StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
rdTree.C
1 // read TTree
2 
3 rdTree(int neve=100, TString Tname0="mc_eveID1", int flag=0, float Emax=40.){
4  gROOT->LoadMacro("bfc.C");
5  bfc(0,"fzin sim_T gen_T","mc_eveID6.fzd");
6 
7  gSystem->Load("EEmc.so");
8  // gStyle->SetPalette(1,0);
9  gStyle->SetOptStat(111111);
10 
11  TString Tname;
12 
13  // Tname="mc_one_pid3eta30.pt1.75";
14  // Tname="/star/data22/MC/balewski/eemc_sim2002/new/"+Tname0;
15  // Tname="/star/data22/MC/balewski/eemc_sim2002b/"+Tname0;
16  Tname=Tname0;
17 
18  printf("read upto %d events from file=%s.root\n",neve,Tname.Data());
19  TFile *f = new TFile(Tname+".root");
20  assert(f->IsOpen());
21  TTree *t4 = (TTree*)f->Get("t4");
22 
23  // create a pointer to an event object. This will be used
24  // to read the branch values.
25  EEevent *event = new EEevent();
26 
27  // if (gROOT->IsBatch()) return; new TBrowser(); t4->StartViewer(); return;
28 
29  // allocate memory for needed branches
30  TClonesArray *secA=new TClonesArray("EEsectorDst",1000);
31 
32  // set the branch address
33  TBranch *BRsec = t4->GetBranch("Sec");
34 
35  BRsec->SetAddress(&secA);
36 
37  Int_t nevent = (Int_t)t4->GetEntries();
38  printf("Total events in TTree=%d\n",nevent);
39 
40  const int nh=6;
41  TH2F *h[nh];
42  h[0]=new TH2F("pr1","EE Pre-1 #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5);
43  h[1]=new TH2F("pr2","EE Pre-2 #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5);
44  h[2]=new TH2F("tw","EE towers #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5);
45  h[3]=new TH2F("Post","EE Post #eta-Bin vs #phi-sector",60,1.,13.,12,0.5,12.5);
46  h[4]=new TH2F;
47  h[5]=new TH2F;
48 
49  TH2F *hde[nh];
50  hde[0]=new TH2F("deD1","EE energy loss : Pre-1 vs. Tower (GeV)",50,.0,20.,50,0.,1.);
51  hde[1]=new TH2F("deD2","EE energy loss : Post vs. Tower (GeV)",50,.0,20.,50,0.,.5);
52  hde[2]=new TH2F("deD3","EE energy loss : Pre-1 vs. Post (GeV)",50,.0,1.,50,0.,.5);
53  hde[3]=new TH2F;
54  hde[4]=new TH2F;
55  hde[5]=new TH2F;
56 
57 
58  TH1F *hde1[nh];
59  hde1[0]=new TH1F("dePr1","EE energy deposit : Pre-1 (GeV)",500,.0,Emax/150.);
60  hde1[1]=new TH1F("dePr2","EE energy deposit : Pre-2 (GeV)",500,.0,Emax/150.);
61  hde1[2]=new TH1F("deTw","EE energy deposit : Tower (GeV)",500,.0,Emax/7.);
62  hde1[3]=new TH1F("dePo","EE energy deposit : Post (GeV)",200,.0,Emax/150.);
63  hde1[4]=new TH1F("deU","EE energy deposit : smdU (GeV)",300,.0,Emax/70.);
64  hde1[5]=new TH1F("deV","EE energy deposit : smdV (GeV)",300,.0,Emax/70.);
65  TH1F *hdeSum=new TH1F("deT+S","EE energy deposit : Tower+smdU+V (GeV)",600,.0,Emax/7.);
66 
67  int i;
68  for(i=0;i<nh;i++) {
69  h[i]->GetXaxis()->SetTitle("#phi sectors [1-12]");
70  h[i]->GetYaxis()->SetTitle("#eta bins (bin1 #rightarrow #eta=2.0)");
71  h[i]->GetZaxis()->SetTitle("Energy deposit (GeV)");
72  }
73 
74  // ........... LOOP OVER EVENTS ............
75 
76  for (Int_t ie=0;ie<nevent;ie++) {
77  if(ie>=neve) break;
78  int i;
79  //for(i=0;i<nh;i++) h[i]->Reset(); // clear eve-by-eve histo
80 
81  //read this branch only
82  BRsec->GetEntry(ie);
83 
84  // if(ie%20==0) printf("\n\iEve=%d eveID=%d, eveType=%d, nSec=%d with data :\n",ie,eveID,eveType,secA->GetEntries());
85  if(ie%20==0) printf("\n\iEve=%d nSec=%d with data \n",ie,secA->GetEntries());
86 
87  const int nz=6;
88  float enerS[nz]; //4 Tower~like detectors +2 Smd
89  {int i; for(i=0;i<nz;i++) enerS[i]=0;}
90 
91  int is;
92  for(is=0;is<secA->GetEntries();is++) {
93  EEsectorDst *sec=(EEsectorDst*)secA->At(is);
94  if(ie<1) sec->print();
95 
96  TClonesArray *hitA;
97  int ih;
98 
99  TClonesArray *hitAA[]={sec->getPre1Hits(),sec->getPre2Hits(),sec->getTwHits(),sec->getPostHits(),sec->getSmdUHits(),sec->getSmdVHits()};
100  int iz;
101  for(iz=0;iz<4;iz++) { // over tower/pre/post
102  hitA=hitAA[iz];
103  if(ie<1) printf(" sectorID=%d iz=%d nHit=%d :\n",sec->getID(),iz,hitA->GetEntries());
104  for(ih=0;ih<hitA->GetEntries();ih++) {
105  char sub;
106  int eta;
107  float ener;
108  EEtwHitDst *hit=(EEtwHitDst*)hitA->At(ih);
109  hit->get(sub,eta,ener);
110  if(ie<1) printf(" ih=%d sec=%d sub=%c etaBin=%d ener=%f\n",ih, sec->getID(), sub, eta,ener);
111  float x=sec->getID()+0.1+0.2*(sub-'A');
112  h[iz]->Fill(x,eta,ener);
113  enerS[iz]+=ener;
114  // hit->print();
115  }
116  }// end of loop over pre1/2/Tw/Post
117 
118 
119  for(iz=4;iz<6;iz++) { // over tower/pre/post
120  hitA=hitAA[iz];
121  if(ie<1) printf(" sectorID=%d iz=%d nHit=%d :\n",sec->getID(),iz,hitA->GetEntries());
122  for(ih=0;ih<hitA->GetEntries();ih++) {
123  EEsmdHitDst *hit2=(EEsmdHitDst*)hitA->At(ih);
124  int strip;
125  float ener;
126  hit2->get(strip,ener);
127  if(ie<1) printf(" ih=%d strip=%d etaBin=%d ener=%f\n",ih, sec->getID(), strip,ener);
128  // h[iz]->Fill(x,eta,ener);
129  enerS[iz]+=ener;
130  // hit->print();
131  }
132  }// end of loop over pre1/2/Tw/Post
133 
134 
135  }// end of loop over sector
136 
137  hde[0]->Fill(enerS[2],enerS[0]);
138  hde[1]->Fill(enerS[2],enerS[3]);
139  hde[2]->Fill(enerS[3],enerS[0]);
140 
141  for(i=0;i<nz;i++) hde1[i]->Fill(enerS[i]);
142 
143  hdeSum->Fill(enerS[2]+enerS[4]+enerS[5]);
144 
145  if(ie<1) for(i=0;i<nz;i++) printf("enerS[%d]=%f\n",i,enerS[i]);
146 
147  }// end of loop over events
148 
149 
150  printf("Total events in B TTree=%d\n",nevent);
151 
152  TFile fh(Tname0+".hist.root","recreate");
153 
154  hde1[2]->Fit("gaus");// tottal tower energy
155  hde1[2]->GetFunction("gaus")->SetLineColor(kRed);
156 
157  hdeSum->Fit("gaus");// tottal tower energy
158  hdeSum->GetFunction("gaus")->SetLineColor(kRed);
159 
160  for(i=0;i<nz;i++) {
161  h[i]->Write();
162  hde[i]->Write();
163  hde1[i]->Write();
164  }
165  hdeSum->Write();
166  fh.ls();
167 
168  TCanvas *can = new TCanvas("cx","main plot",600,800);
169  can->Range(0,0,1,1);
170 
171  TPad *pad0 = new TPad("pad0", "apd0",0.0,0.90,1.,1.);
172  pad0->Draw();
173  pad0->cd();
174  pad0->Range(0,0,1,1);
175 
176  TLatex *tlx=new TLatex(0.12,0.5,strstr(Tname0.Data(),"/mc_"));
177  tlx->SetTextSize(0.4);
178  tlx->Draw();
179 
180  can->cd();
181  TPad *c1 = new TPad("pad1", "apd1",0.0,0.0,1.,.90);
182  c1->Draw();
183 
184  c1->Divide(2,4);
185  for(i=0;i<nz;i++) {
186  c1->cd(i+1);
187  hde1[i]->Draw();
188  // if(i<=3)
189  gPad->SetLogy();
190  }
191  c1->cd(7);
192  hdeSum->Draw();
193  gPad->SetLogy();
194 
195 
196  if(flag) {
197  TString psFile=Tname0+".ps";
198  can->Print(psFile.Data());
199  printf("## cp %s target ",psFile.Data());
200  }
201  return;
202  c1=new TCanvas();
203  c1->Divide(2,2);
204  for(i=0;i<3;i++) {
205  c1->cd(i+1);
206  hde[i]->Draw("box");
207  }
208  c1->cd(4);
209  h[2]->Draw("box");
210  if(flag)c1->Print("fig2.ps");
211 
212 
213  c1=new TCanvas();
214  c1->Divide(2,2);
215  printf("Content of all used event DISPLAY\n");
216  for(i=0;i<nh;i++) {
217  c1->cd(i+1);
218  h[i]->Draw("LEGO2");
219  h[i]->Print();
220  }
221  if(flag) c1->Print("fig3.ps");
222 
223 }
224 
Generated by   doxygen 1.8.5