StRoot: StEEmcUtil/macros/rdTree.C Source File

00001 // read TTree
00002 
00003 rdTree(int neve=100,  TString Tname0="mc_eveID1", int flag=0, float Emax=40.){  
00004  gROOT->LoadMacro("bfc.C"); 
00005   bfc(0,"fzin sim_T gen_T","mc_eveID6.fzd");
00006 
00007   gSystem->Load("EEmc.so"); 
00008   //  gStyle->SetPalette(1,0);
00009   gStyle->SetOptStat(111111);
00010 
00011  TString Tname;
00012 
00013  // Tname="mc_one_pid3eta30.pt1.75";
00014  // Tname="/star/data22/MC/balewski/eemc_sim2002/new/"+Tname0;
00015  // Tname="/star/data22/MC/balewski/eemc_sim2002b/"+Tname0;
00016  Tname=Tname0;
00017 
00018  printf("read upto %d events from file=%s.root\n",neve,Tname.Data());
00019  TFile *f = new TFile(Tname+".root");
00020  assert(f->IsOpen());
00021  TTree *t4 = (TTree*)f->Get("t4");
00022   
00023   // create a pointer to an event object. This will be used
00024   // to read the branch values.
00025   EEevent *event = new EEevent();
00026   
00027   //  if (gROOT->IsBatch()) return;  new TBrowser();  t4->StartViewer();  return;
00028 
00029   // allocate memory for needed branches 
00030   TClonesArray *secA=new TClonesArray("EEsectorDst",1000);
00031 
00032   // set the branch address
00033   TBranch *BRsec = t4->GetBranch("Sec");
00034 
00035   BRsec->SetAddress(&secA);
00036   
00037   Int_t nevent = (Int_t)t4->GetEntries();
00038   printf("Total events in TTree=%d\n",nevent);
00039 
00040   const int nh=6;
00041   TH2F *h[nh];
00042   h[0]=new TH2F("pr1","EE Pre-1  #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5); 
00043   h[1]=new TH2F("pr2","EE Pre-2  #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5); 
00044   h[2]=new TH2F("tw","EE towers  #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5); 
00045   h[3]=new TH2F("Post","EE Post  #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5); 
00046   h[4]=new TH2F;
00047   h[5]=new TH2F;
00048 
00049   TH2F *hde[nh];
00050   hde[0]=new TH2F("deD1","EE energy loss : Pre-1  vs. Tower  (GeV)",50,.0,20.,50,0.,1.); 
00051   hde[1]=new TH2F("deD2","EE energy loss : Post  vs. Tower  (GeV)",50,.0,20.,50,0.,.5); 
00052   hde[2]=new TH2F("deD3","EE energy loss : Pre-1  vs. Post  (GeV)",50,.0,1.,50,0.,.5); 
00053   hde[3]=new TH2F;
00054   hde[4]=new TH2F;
00055   hde[5]=new TH2F;
00056 
00057 
00058   TH1F *hde1[nh];
00059   hde1[0]=new TH1F("dePr1","EE energy deposit : Pre-1 (GeV)",500,.0,Emax/150.);
00060   hde1[1]=new TH1F("dePr2","EE energy deposit : Pre-2 (GeV)",500,.0,Emax/150.);
00061   hde1[2]=new TH1F("deTw","EE energy deposit : Tower (GeV)",500,.0,Emax/7.);
00062   hde1[3]=new TH1F("dePo","EE energy deposit : Post (GeV)",200,.0,Emax/150.);
00063   hde1[4]=new TH1F("deU","EE energy deposit : smdU (GeV)",300,.0,Emax/70.);
00064   hde1[5]=new TH1F("deV","EE energy deposit : smdV (GeV)",300,.0,Emax/70.);
00065   TH1F *hdeSum=new TH1F("deT+S","EE energy deposit : Tower+smdU+V (GeV)",600,.0,Emax/7.);
00066 
00067   int i;
00068   for(i=0;i<nh;i++) {
00069     h[i]->GetXaxis()->SetTitle("#phi sectors [1-12]");
00070     h[i]->GetYaxis()->SetTitle("#eta bins (bin1 #rightarrow #eta=2.0)");
00071     h[i]->GetZaxis()->SetTitle("Energy deposit (GeV)");
00072   }
00073 
00074   // ...........  LOOP OVER EVENTS ............
00075 
00076   for (Int_t ie=0;ie<nevent;ie++) {    
00077     if(ie>=neve) break;
00078     int i;
00079     //for(i=0;i<nh;i++) h[i]->Reset(); // clear eve-by-eve histo  
00080     
00081     //read this branch only
00082     BRsec->GetEntry(ie);  
00083     
00084     // if(ie%20==0) printf("\n\iEve=%d eveID=%d, eveType=%d, nSec=%d with data :\n",ie,eveID,eveType,secA->GetEntries());
00085    if(ie%20==0) printf("\n\iEve=%d  nSec=%d with data \n",ie,secA->GetEntries());
00086     
00087     const int nz=6;
00088     float enerS[nz]; //4 Tower~like detectors +2 Smd
00089     {int i; for(i=0;i<nz;i++) enerS[i]=0;}
00090 
00091     int is;
00092     for(is=0;is<secA->GetEntries();is++) {
00093       EEsectorDst *sec=(EEsectorDst*)secA->At(is);
00094       if(ie<1) sec->print();
00095       
00096       TClonesArray *hitA;
00097       int ih;
00098   
00099       TClonesArray *hitAA[]={sec->getPre1Hits(),sec->getPre2Hits(),sec->getTwHits(),sec->getPostHits(),sec->getSmdUHits(),sec->getSmdVHits()};
00100       int iz;
00101       for(iz=0;iz<4;iz++) { // over tower/pre/post
00102         hitA=hitAA[iz];
00103         if(ie<1)    printf("  sectorID=%d  iz=%d nHit=%d :\n",sec->getID(),iz,hitA->GetEntries());      
00104         for(ih=0;ih<hitA->GetEntries();ih++) {
00105           char sub;
00106           int eta;
00107           float ener;
00108           EEtwHitDst *hit=(EEtwHitDst*)hitA->At(ih);
00109           hit->get(sub,eta,ener);
00110          if(ie<1) printf("    ih=%d sec=%d sub=%c etaBin=%d ener=%f\n",ih, sec->getID(), sub, eta,ener);
00111           float x=sec->getID()+0.1+0.2*(sub-'A');
00112           h[iz]->Fill(x,eta,ener);
00113           enerS[iz]+=ener;
00114           // hit->print();      
00115         }
00116       }// end of loop over pre1/2/Tw/Post       
00117       
00118      
00119       for(iz=4;iz<6;iz++) { // over tower/pre/post
00120         hitA=hitAA[iz];
00121          if(ie<1) printf("  sectorID=%d  iz=%d nHit=%d :\n",sec->getID(),iz,hitA->GetEntries());
00122         for(ih=0;ih<hitA->GetEntries();ih++) {
00123           EEsmdHitDst *hit2=(EEsmdHitDst*)hitA->At(ih);
00124           int strip;
00125           float ener;
00126           hit2->get(strip,ener);
00127           if(ie<1) printf("    ih=%d strip=%d etaBin=%d ener=%f\n",ih, sec->getID(), strip,ener);
00128           //      h[iz]->Fill(x,eta,ener);
00129           enerS[iz]+=ener;
00130           // hit->print();      
00131         }
00132       }// end of loop over pre1/2/Tw/Post       
00133       
00134       
00135     }// end of loop over sector
00136     
00137     hde[0]->Fill(enerS[2],enerS[0]);
00138     hde[1]->Fill(enerS[2],enerS[3]);
00139     hde[2]->Fill(enerS[3],enerS[0]);
00140 
00141     for(i=0;i<nz;i++) hde1[i]->Fill(enerS[i]);
00142     
00143     hdeSum->Fill(enerS[2]+enerS[4]+enerS[5]);
00144 
00145     if(ie<1) for(i=0;i<nz;i++) printf("enerS[%d]=%f\n",i,enerS[i]);
00146 
00147     }// end of loop over events
00148     
00149 
00150   printf("Total events in B TTree=%d\n",nevent);
00151 
00152   TFile fh(Tname0+".hist.root","recreate");
00153   
00154   hde1[2]->Fit("gaus");// tottal tower energy
00155   hde1[2]->GetFunction("gaus")->SetLineColor(kRed);
00156 
00157   hdeSum->Fit("gaus");// tottal tower energy
00158   hdeSum->GetFunction("gaus")->SetLineColor(kRed);
00159 
00160   for(i=0;i<nz;i++) {
00161     h[i]->Write();
00162     hde[i]->Write();
00163     hde1[i]->Write();
00164   }    
00165   hdeSum->Write();
00166   fh.ls();
00167 
00168   TCanvas *can = new TCanvas("cx","main plot",600,800);
00169   can->Range(0,0,1,1);
00170     
00171   TPad *pad0 = new TPad("pad0", "apd0",0.0,0.90,1.,1.);
00172   pad0->Draw();
00173   pad0->cd();
00174   pad0->Range(0,0,1,1);
00175     
00176   TLatex *tlx=new TLatex(0.12,0.5,strstr(Tname0.Data(),"/mc_"));
00177   tlx->SetTextSize(0.4);
00178   tlx->Draw();
00179     
00180   can->cd();
00181   TPad *c1 = new TPad("pad1", "apd1",0.0,0.0,1.,.90);
00182   c1->Draw();
00183   
00184   c1->Divide(2,4);
00185   for(i=0;i<nz;i++) {
00186     c1->cd(i+1);
00187     hde1[i]->Draw();
00188     //    if(i<=3)
00189     gPad->SetLogy();
00190   }
00191   c1->cd(7);
00192   hdeSum->Draw();
00193   gPad->SetLogy();
00194   
00195 
00196   if(flag) {
00197     TString psFile=Tname0+".ps";
00198     can->Print(psFile.Data());
00199     printf("## cp %s target ",psFile.Data());
00200   }
00201   return;
00202   c1=new TCanvas();
00203   c1->Divide(2,2);
00204   for(i=0;i<3;i++) {
00205     c1->cd(i+1);
00206     hde[i]->Draw("box");
00207   }
00208     c1->cd(4);
00209     h[2]->Draw("box");
00210     if(flag)c1->Print("fig2.ps");
00211     
00212     
00213     c1=new TCanvas();
00214     c1->Divide(2,2);   
00215     printf("Content of all used event DISPLAY\n");
00216     for(i=0;i<nh;i++) {
00217       c1->cd(i+1);
00218       h[i]->Draw("LEGO2");
00219       h[i]->Print();
00220     }
00221     if(flag) c1->Print("fig3.ps");
00222   
00223 }
00224