StRoot: StEEmcUtil/macros/plTw.C Source File

00001 plTw( int pid=14, char *name="deT+S") {   //type of particles, histo name
00002 
00003   const int ng=4;
00004   TGraphErrors *hg[ng];
00005   int nH=InitGraph(hg,ng);
00006   printf("%d graphs initialized\n",nH);
00007   int i;
00008   for(i=0;i<ng;i++) printf("i=%d, add=%0x\n",i,hg[i]);
00009   
00010   gStyle->SetOptStat(111111);
00011   gStyle->SetOptFit(111111);
00012   
00013   //pt bins
00014   float pt[]={0.5,1,2,3,4,5,10,20,30,40,50}; int npt=11;
00015   //float pt[]={0.5,1,10,50}; int npt=4;
00016   
00017   //eta bins 
00018   float eta[]={1.2,1.75};
00019   int neta=2;
00020 
00021   int i=0,k=0;
00022   
00023   for(i=0;i<npt;i++) 
00024         for(k=0;k<neta;k++)
00025       {
00026         float Etot=pt[i]*cosh(eta[k]); // approximate, mass ignored
00027         char fname[1000];
00028         sprintf(fname,"/star/data22/MC/balewski/eemc/sim2002/f/mc_pid%d_pt%.1f_eta%4.2f",pid,pt[i],eta[k]);
00029         //  printf( "file%s Etot=%f  eta=%f\n",fname,Etot,eta[k]);// continue;
00030         // E=pT/sin(theta)
00031         doJob(fname,hg,Etot,k,pid,name);
00032         //      if(i>1) break;
00033         //      return;
00034       }
00035   
00036   //hg[0]->Print();
00037   // hg[2]->Print();
00038   //    return;
00039   char tit[100];
00040   sprintf(tit,"%s-pid%d",name,pid);
00041   
00042   int kkk=0;
00043   if (pid<4) kkk=1;
00044   TCanvas *can = new TCanvas(tit,tit,600,400+400*kkk);
00045   if(kkk)  { can->Divide(1,2); can->cd(1);}
00046    hgDrawE(hg,pid,name);
00047    
00048   if(kkk)  { can->cd(2);
00049   hgDrawSE(hg+2,pid,name);
00050   }
00051   can->Print();
00052   
00053 }
00054 
00055 //--------------------------------------------
00056 //--------------------------------------------
00057 //--------------------------------------------
00058 
00059 void doJob(char *fname0,TGraphErrors **hg, float Etot, int ieta,int pid,char *hname) {
00060   printf("job--> %s Etot=%f  ieta=%d hist=%s\n",fname0,Etot,ieta,hname);
00061   assert(ieta>=0 && ieta<=1);
00062 
00063   TString fname=fname0;
00064   TFile *dir=new TFile(fname+".hist.root");
00065   // dir->ls();
00066   
00067   if(!dir->IsOpen()){ printf("Open failed, take next\n"); return;}
00068   
00069   TH1F* h1=0;
00070   double x,y,ey;
00071   int n=-1;
00072    
00073   h1=(TH1F*) dir->Get(hname);  
00074   assert(h1);
00075   
00076   TF1 *ff=h1->GetFunction("gaus");  assert(ff);
00077 
00078   // total energy  ........................
00079   x=Etot;
00080   
00081   if(pid<4) { // use relative energy for gam & ele
00082     y=ff->GetParameter(1)/Etot;
00083     ey=ff->GetParError(1)/Etot;
00084   }
00085   
00086   if(pid>4) { // use relative MEAN & RMS for others
00087     y=h1->GetMean();
00088     ey=h1->GetRMS();
00089   }
00090   
00091   
00092   TGraphErrors *gr=hg[0+ieta];
00093   n=gr->GetN();
00094   gr->SetPoint(n,x,y);
00095   gr->SetPointError(n,0.,ey);
00096   
00097     if(pid>4) return;
00098 
00099   // energy resolution ........................
00100   x=1/sqrt(Etot);
00101   y=ff->GetParameter(2)/Etot;
00102   ey=ff->GetParError(2)/Etot;
00103   printf("Etot=%f x=%f y=%f ey=%g\n",Etot,x,y);
00104 
00105   TGraphErrors *gr=hg[2+ieta];
00106   n=gr->GetN();
00107   gr->SetPoint(n,x,y);
00108   gr->SetPointError(n,0.,ey);
00109   
00110 }
00111 
00112 
00113 //------------------------------------------------
00114 //------------------------------------------------
00115 //------------------------------------------------
00116 
00117 int InitGraph( TGraphErrors **hg,int ng){
00118   char *name[]={"reco E or E/E, eta=1.20"
00119                 ,"reco E or E/E, eta=1.75"
00120                 ,"reco sigE , eta=1.20"
00121                 ,"reco sigE , eta=1.75"
00122   }
00123   
00124   int j;
00125   
00126   {int j; for(j=0;j<ng;j++) hg[j]=0;  }
00127 
00128   int kk=0;
00129   for(j=0;j<ng;j++) {
00130     int i=j; // tmp
00131     TGraphErrors *gr=new TGraphErrors();
00132     hg[kk]=gr;
00133     gr->SetMarkerStyle(28);
00134     gr->SetMarkerColor(kGreen);
00135     gr->SetName(name[i]);
00136     int icol=kBlue, isym=28;
00137     if(j%2==0 )isym=25;
00138     if(j%2==1 )isym=24;
00139 
00140     if(j%2==0) icol=kRed;
00141     if(j%2==1) icol=kGreen;
00142  
00143     gr->SetMarkerStyle(isym);
00144     gr->SetMarkerColor(icol);
00145     gr->SetLineColor(icol);
00146     kk++;
00147   }
00148   printf("Initialized %d Graphs\n",kk);
00149   return kk;
00150 }
00151 
00152 
00153 //------------------------------------------------
00154 //------------------------------------------------
00155 //------------------------------------------------
00156 
00157 
00158 hgDrawE(TGraphErrors **hg, int pid, char *name) {   //type of particles, histo name
00159   char *pname[]={"gamma","xx","electron","xx","xx","muon","xx","xx","pion-","xx","xx","xx","xx","proton"};
00160   
00161   TGraphErrors *gr=0;
00162   gr=hg[0];
00163   gr->Draw("AP");
00164   TString titHead;
00165 
00166   if(pid>4 && strstr(name,"deT+S")){
00167     titHead="mean #pm RMS E_{TOWER+SMD} (GeV)";
00168     gr->SetMinimum( -0.2);
00169  }
00170   if(pid>4 && strstr(name,"deTw")){
00171     titHead="mean #pm RMS E_{TOWER only} (GeV)";
00172     gr->SetMinimum( -0.2);
00173   }
00174 
00175   if(pid<4 && strstr(name,"deT+S")) {
00176     gr->SetMaximum( 0.054);  gr->SetMinimum( 0.048);
00177     titHead="Gauss E_{REC}/E_{INC}  #pm #sigma/E_{INC}  for  TOWER+SMD ";
00178   }
00179   if(pid<4 && strstr(name,"deTw")) {
00180     gr->SetMaximum( 0.042);  gr->SetMinimum( 0.037);
00181     titHead="Gauss E_{REC}/E_{INC}  #pm #sigma/E_{INC}  for  TOWER only ";
00182   }
00183 
00184   gPad->SetLogx();
00185   gr->GetXaxis()->SetLimits(0.7,170.);
00186   gr->GetXaxis()->SetTitle("Incident Energy (GeV)");
00187   
00188   char tit[1000];
00189   sprintf(tit," %s;    PID=%d\n",titHead.Data(),pid);
00190   gr->SetTitle(tit);
00191   hg[1]->Draw("P");
00192 
00193   lg=new TLegend(0.3,0.65,0.55,0.8);
00194   TString head="Incident    ";
00195   head+=pname[pid-1];
00196   lg->SetHeader(head);
00197   lg->AddEntry(hg[0],"  #eta=1.20","LP");
00198   lg->AddEntry(hg[1],"  #eta=1.75","LP");
00199   lg->Draw();
00200 
00201 
00202   
00203 }
00204 
00205 
00206 
00207 
00208 hgDrawSE(TGraphErrors **hg, int pid, char *name) {   //type of particles
00209   char *pname[]={"gamma","xx","electron","xx","xx","muon","xx","xx","pion-","xx","xx","xx","xx","proton"};
00210   
00211   TGraphErrors *gr=0;
00212   gr=hg[0];
00213   gr->Draw("AP");
00214   gr->Fit("pol1");
00215 
00216   TString titHead;
00217   if(strstr(name,"deT+S")) titHead="  #sigma/ E_{INC} for TOWER+SMD";
00218   if(strstr(name,"deTw")) titHead="  #sigma/ E_{INC} for TOWER only";
00219 
00220   gr->SetMinimum( 0.0);
00221 
00222   if(pid>4){
00223     gr->SetMaximum( 0.004); 
00224   }
00225 
00226   if(pid<4){
00227     gr->SetMaximum( 0.01); 
00228   }
00229 
00230   gr->GetXaxis()->SetTitle("1 /  #sqrt{E_{INC}/GeV}");
00231   
00232   char tit[1000];
00233   sprintf(tit,"Gauss %s;    PID=%d\n",titHead.Data(),pid);
00234   gr->SetTitle(tit);
00235 
00236 
00237   hg[1]->Draw("P");
00238 
00239   lg=new TLegend(0.6,0.15,0.85,0.3);
00240   TString head="Incident    ";
00241   head+=pname[pid-1];
00242   lg->SetHeader(head);
00243   lg->AddEntry(hg[0],"  #eta=1.20 +FIT","LP");
00244   lg->AddEntry(hg[1],"  #eta=1.75","LP");
00245   lg->Draw();
00246 
00247 
00248   
00249 }
00250 
00251 
00252