source: ELYSE/HEAD/scripts/ROOT/analysis.C @ 286

Last change on this file since 286 was 286, checked in by campagne, 17 years ago

ELYSE sauvegarde provisoire (JEC)
File size: 9.4 KB
Line 
1//
2//  This script shows how to read, with ROOT, a ELYSE.root
3// file having the "TTree within TTree" of the ELYSE_batch
4// progrem (see Analysis.cxx).
5//
6//  Usage :
7//    OS> <setup ROOT>
8//  ( OS> root batch_dummy.C )
9//    OS> root analyis.C
10//
11//  G.Barrand
12//
13
14//////////////////////////////////////////////////////////////////////////////
15bool plot(
16 TH1D& aHisto1D
17,TH2D& aHisto2D
18) 
19//////////////////////////////////////////////////////////////////////////////
20//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
21{
22  TCanvas* plotter = new TCanvas("canvas","",10,10,800,600);
23  plotter->Divide(1,2); 
24
25  plotter->cd(1);
26  aHisto1D.Draw();
27
28  plotter->cd(2);
29  aHisto2D.Draw();
30
31  plotter->Update();
32
33  return true;
34}
35//////////////////////////////////////////////////////////////////////////////
36void analysis() 
37//////////////////////////////////////////////////////////////////////////////
38//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
39{
40
41  ////////////////////////////////////////////////////////
42  // Create histograms ///////////////////////////////////
43  ////////////////////////////////////////////////////////
44  TH1D* hits_times = new TH1D("hits_times","Hits times",100,0,3000);
45
46  TH2D* digits_time_pe = 
47    new TH2D("digits_pe_time","Digits PE time",100,0,3000,100,0,10);
48
49  int hit_count = 0;
50  bool dump = true;
51  //bool dump = false;
52
53  ////////////////////////////////////////////////////////
54  /// Read data an file histos ///////////////////////////
55  ////////////////////////////////////////////////////////
56  TFile* f = new TFile("ELYSE.root");
57  TTree* tEvent = (TTree*)f->Get("Event");
58
59 {TObjArray* brs = tEvent->GetListOfBranches();
60  for(int i=0;i<brs->GetEntries();i++) {
61    TBranch* b = (TBranch*)brs->At(i);
62    std::cout << " Event branch : " << b->GetName()
63              << std::endl;   
64  }}
65
66  //declare Non subTuple part
67  Int_t eventId,inputEvtId, interMode, vtxVol;
68  tEvent->SetBranchAddress("eventId",&eventId);
69  tEvent->SetBranchAddress("inputEvtId",&inputEvtId);
70  tEvent->SetBranchAddress("interMode",&interMode);
71  tEvent->SetBranchAddress("vtxVol",&vtxVol);
72
73  Int_t nPart,leptonIndex, protonIndex;
74  tEvent->SetBranchAddress("nPart",&nPart);
75  tEvent->SetBranchAddress("leptonIndex",&leptonIndex);
76  tEvent->SetBranchAddress("protonIndex",&protonIndex);
77
78  Int_t nHits; 
79  tEvent->SetBranchAddress("nHits",&nHits);
80
81  Int_t nDigits;
82  Double_t sumPE;
83  tEvent->SetBranchAddress("nDigits",&nDigits);
84  tEvent->SetBranchAddress("sumPE",&sumPE);
85
86  Int_t nEvent = (Int_t)tEvent->GetEntries(); 
87  //nEvent = 1;
88  //nEvent = 10000;
89  std::cout << " nEvents = " << nEvent << std::endl;
90
91  for (Int_t i=0; i<nEvent; ++i) {
92
93    TTree* Event_vtxPos = new TTree();
94    tEvent->SetBranchAddress("vtxPos",&Event_vtxPos);
95
96    TTree* Event_track = new TTree();
97    tEvent->SetBranchAddress("track",&Event_track);
98
99    TTree* Event_hit = new TTree();
100    tEvent->SetBranchAddress("hit",&Event_hit);
101   
102    TTree* Event_digit = new TTree();
103    tEvent->SetBranchAddress("digit",&Event_digit);
104
105    tEvent->GetEntry(i);
106
107    //{TObjArray* brs = Event_track->GetListOfBranches();
108    //for(int i=0;i<brs->GetEntries();i++) {
109    //  TBranch* b = (TBranch*)brs->At(i);
110    //  std::cout << " Event_track branch : " << b->GetName() << std::endl;   
111    //  }}
112
113    // Bind sub tuple variables
114    Double_t vx,vy,vz;
115    Event_vtxPos->SetBranchAddress("x",&vx);
116    Event_vtxPos->SetBranchAddress("y",&vy);
117    Event_vtxPos->SetBranchAddress("z",&vz);
118
119    Int_t pId,parent;
120    Float_t timeStart;
121    Double_t mass, pTot, ETot;
122    Int_t startVol, stopVol;
123    Event_track->SetBranchAddress("pId",&pId);
124    Event_track->SetBranchAddress("parent",&parent);
125    Event_track->SetBranchAddress("timeStart",&timeStart);
126    Event_track->SetBranchAddress("mass",&mass);
127    Event_track->SetBranchAddress("pTot",&pTot);
128    Event_track->SetBranchAddress("ETot",&ETot);
129    Event_track->SetBranchAddress("startVol",&startVol);
130    Event_track->SetBranchAddress("stopVol",&stopVol);
131
132    Int_t tubeId_hit; //JEC 16/1/06
133    Int_t totalPE;   
134    Event_hit->SetBranchAddress("tubeId",&tubeId_hit); //JEC 16/1/06   
135    Event_hit->SetBranchAddress("totalPE",&totalPE);
136
137    Int_t tubeId;
138    Double_t digit_pe, digit_time;
139    Event_digit->SetBranchAddress("tubeId",&tubeId);
140    Event_digit->SetBranchAddress("pe",&digit_pe);
141    Event_digit->SetBranchAddress("time",&digit_time);
142
143    if(dump)
144    std::cout << ">>>>>>>>>>>>> Event{" << i << "}: "
145              << " evt Id " << eventId
146    //        << " evt Input Id " << inputEvtId
147    //        << "\n interaction mode " << interMode
148    //        << " start in volume " << vtxVol << "\n"
149              <<" #tracks: " << nPart
150              <<" #hits: " << nHits
151              <<" #digits: " << nDigits
152              << std::endl;
153
154    Int_t nVtx = (Int_t)Event_vtxPos->GetEntries();
155    Int_t nTracks = (Int_t)Event_track->GetEntries();
156    Int_t nTubeHits = (Int_t)Event_hit->GetEntries();
157    Int_t nTubeDigits = (Int_t)Event_digit->GetEntries();
158
159    if(dump)
160    std::cout << "Verif :" 
161              << " nVtx = " << nVtx
162              << " nTracks = " << nTracks
163              << " nTube Hits = " << nTubeHits
164              << " nTube Digits = " << nTubeDigits
165              << std::endl;
166
167    for (Int_t j=0; j<nTracks; ++j) {
168
169      //Sub tuples of Event_tracks :
170      TTree* Track_direction = new TTree();
171      Event_track->SetBranchAddress("direction",&Track_direction);
172      TTree* Track_momentum = new TTree();
173      Event_track->SetBranchAddress("momentum",&Track_momentum);
174      TTree* Track_startPos = new TTree();
175      Event_track->SetBranchAddress("startPos",&Track_startPos);
176      TTree* Track_stopPos = new TTree();
177      Event_track->SetBranchAddress("stopPos",&Track_stopPos);
178
179      Event_track->GetEntry(j);
180
181      Double_t dx,dy,dz;
182      Track_direction->SetBranchAddress("dx",&dx);
183      Track_direction->SetBranchAddress("dy",&dy);
184      Track_direction->SetBranchAddress("dz",&dz);
185
186      Double_t px,py,pz;
187      Track_momentum->SetBranchAddress("px",&px);
188      Track_momentum->SetBranchAddress("py",&py);
189      Track_momentum->SetBranchAddress("pz",&pz);
190
191      Double_t start_x,start_y,start_z;
192      Track_startPos->SetBranchAddress("x",&start_x);
193      Track_startPos->SetBranchAddress("y",&start_y);
194      Track_startPos->SetBranchAddress("z",&start_z);
195
196      Double_t stop_x,stop_y,stop_z;
197      Track_stopPos->SetBranchAddress("x",&stop_x);
198      Track_stopPos->SetBranchAddress("y",&stop_y);
199      Track_stopPos->SetBranchAddress("z",&stop_z);
200
201      // One entry only :
202      //std::cout << "debug : "
203      //          << " " << Track_direction->GetEntries()
204      //          << " " << Track_momentum->GetEntries()
205      //          << " " << Track_startPos->GetEntries()
206      //          << " " << Track_stopPos->GetEntries()
207      //          << std::endl;
208      if(Track_direction->GetEntries()==1) Track_direction->GetEntry(0);
209      if(Track_momentum->GetEntries()==1)   Track_momentum->GetEntry(0);
210      if(Track_startPos->GetEntries()==1) Track_startPos->GetEntry(0);
211      if(Track_stopPos->GetEntries()==1) Track_stopPos->GetEntry(0);
212
213      if(dump)
214      std::cout << "----> Tk{"<<j<<"}: " 
215                << " pId " << pId
216                << " parent " << parent
217                << " creation time " << timeStart
218                << " Volumes " << startVol << " " << stopVol << "\n"
219                << " Start Pos (" << start_x << "," << start_y << "," << start_z << ")\n"
220                << " Stop Pos (" << stop_x << "," << stop_y << "," << stop_z << ")\n"
221                << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n"
222                << " m " << mass
223                << " ETot " << ETot
224                << " pTot " << pTot
225                << " px,py,pz " << px << " " << py << " " << pz << "\n"
226                << std::endl;
227
228      delete Track_direction;
229      delete Track_momentum;
230      delete Track_startPos;
231      delete Track_stopPos;
232
233    }//loop on Tracks
234
235    //--------
236    // The Hits
237    //--------
238
239
240    for (Int_t k=0; k<nTubeHits; ++k) {
241     
242      TTree* Hit_pe = new TTree();
243      Event_hit->SetBranchAddress("pe",&Hit_pe);
244 
245      Event_hit->GetEntry(k);
246
247      Float_t hit_time;
248      //Float_t trk_length; //NV 13/6/06
249      Hit_pe->SetBranchAddress("time",&hit_time);
250      //Hit_pe->SetBranchAddress("length",&trk_length);
251
252      //JEC 16/1/06 add the tubeId_hit info
253      if(dump)
254      std::cout << "----> Hit{"<<k<<"}: tube[" << tubeId_hit << "] total #PE " << totalPE << std::endl;
255
256      for (Int_t ki=0; ki<Hit_pe->GetEntries(); ++ki) {
257        Hit_pe->GetEntry(ki);
258
259        if(dump)
260        std::cout << "<" << hit_time << ">" 
261          /* << "<" << trk_length << ">" */ ;
262
263        hits_times->Fill(hit_time);
264
265        hit_count++;
266      }
267      if(dump)
268      std::cout << std::endl;
269
270      delete Hit_pe;
271    }//Loop on Hits
272
273    //--------
274    // The Digits
275    //--------
276    for (Int_t l=0; l<nTubeDigits; ++l) {
277      Event_digit->GetEntry(l);
278     
279      if(dump)
280      std::cout << "----> Digit{"<<l<<"}: " 
281                << "tube[" << tubeId << "] = " 
282                << " pe: " << digit_pe
283                << " time: " << digit_time
284                << std::endl;
285     
286      digits_time_pe->Fill(digit_time,digit_pe);
287
288    }//Loop on Digits
289
290    delete Event_vtxPos;
291    delete Event_track;
292    delete Event_hit;
293    delete Event_digit;
294
295  }//loop on event
296
297  if(dump)
298  std::cout << " nEvents = " << nEvent << " hits = " << hit_count << std::endl;
299
300  ////////////////////////////////////////////////////////
301  // Plot histograms /////////////////////////////////////
302  ////////////////////////////////////////////////////////
303  plot(*hits_times,*digits_time_pe);
304
305}//read_event
Note: See TracBrowser for help on using the repository browser.