Read events are enumerated from 1
[u/mrichter/AliRoot.git] / ANALYSIS / AliReader.cxx
1 #include "AliReader.h"
2 //_________________________________________________________________________
3 ///////////////////////////////////////////////////////////////////////////
4 //
5 // class AliReader
6 //
7 // Reader Base class 
8 // Reads particles and tracks and
9 // puts them to the AliAOD objects and eventually, if needed, buffers AliAODs in AliAODRun(s)
10 //
11 // User loops over events calling method Next. In case of success this method returns 0.
12 // In case of error or if there is no more events to read, non-0 value is returned
13 //
14 // Reading can be rewound to the beginning using method Rewind.
15 //
16 // Tracks are read to the fEventRec (contains reconstructed tracks) 
17 // and fEventSim (corresponding MC simulated data) data members,
18 // that are of the type AliAOD. 
19 //
20 // If a given reader has ability of reading both, reconstructed and simulated data, 
21 // these are structured in AODs so a "n'th" simulated particle 
22 // (the one stored in the fEventSim at slot n) 
23 // corresponds to the n'th reconstructed track (the one stored in the fEventRec at slot n).
24 //
25 // The same reconstructed track can be present more than ones in the AOD,
26 // but with a different PID. In this case
27 // pointer to the corresponding MC simulated particles is also present more than ones.
28 // This situation happens if you want to read all particles 
29 // with PID probability of being , e.g.,  pion higher than 60%
30 // and being kaon higher than 40%. Than, if a given track has probability Ppid(pi)=52% and Ppid(K)=48% 
31 // than it is read twise.
32 //
33 // Provides functionality for both buffering and non-buffering reading
34 // This can be switched on/off via method SetEventBuffering(bool)
35 // The main method that inheriting classes need to implement is ReadNext()
36 // that read next event in queue.
37 //
38 // The others are:
39 // Bool_t  ReadsRec() const; specifies if reader is able to read simulated particles
40 // Bool_t  ReadsSim() const; specifies if reader is able to read reconstructed tracks
41 // void    Rewind(); rewind reading to the beginning
42 //
43 // This class provides full functionality for reading from many sources
44 // User can provide TObjArray of TObjStrings (SetDirs method or via parameter 
45 // in the constructor) which desribes paths of directories to search data in.
46 // If none specified current directory is searched.
47 // 
48 // Piotr.Skowronski@cern.ch
49 //
50 ///////////////////////////////////////////////////////////////////////////
51
52 #include <TString.h>
53 #include <TObjString.h>
54 #include <TObjArray.h>
55 #include <TClass.h>
56 #include <TRandom.h>
57 #include <TH1.h>
58
59 #include "AliAODParticleCut.h"
60 #include "AliAOD.h"
61 #include "AliAODRun.h"
62  
63 ClassImp(AliReader)
64 //pure virtual
65     
66 /*************************************************************************************/
67
68 AliReader::AliReader():
69  fCuts(new TObjArray()),
70  fDirs(0x0),
71  fCurrentEvent(0),
72  fCurrentDir(0),
73  fNEventsRead(0),
74  fEventRec(0x0),
75  fEventSim(0x0),
76  fRunSim(0x0),
77  fRunRec(0x0),
78  fIsRead(kFALSE),
79  fBufferEvents(kFALSE),
80  fBlend(kFALSE),
81  fFirst(0),
82  fLast(0),
83  fTrackCounter(0x0)
84 {
85 //constructor
86 }
87 /*************************************************************************************/
88
89 AliReader::AliReader(TObjArray* dirs):
90  fCuts(new TObjArray()),
91  fDirs(dirs),
92  fCurrentEvent(0),
93  fCurrentDir(0),
94  fNEventsRead(0),
95  fEventRec(0x0),
96  fEventSim(0x0),
97  fRunSim(0x0),
98  fRunRec(0x0),
99  fIsRead(kFALSE),
100  fBufferEvents(kFALSE),
101  fBlend(kFALSE),
102  fFirst(0),
103  fLast(0),
104  fTrackCounter(0x0)
105 {
106 //ctor with array of directories to read as parameter
107 }
108 /*************************************************************************************/
109 AliReader::AliReader(const AliReader& in):
110  TNamed(in),
111  fCuts((in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0),
112  fDirs((in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0),
113  fCurrentEvent(0),
114  fCurrentDir(0),
115  fNEventsRead(0),
116  fEventRec(0x0),
117  fEventSim(0x0),
118  fRunSim(0x0),
119  fRunRec(0x0),
120  fIsRead(kFALSE),
121  fBufferEvents(in.fBufferEvents),
122  fBlend(in.fBlend),
123  fFirst(in.fFirst),
124  fLast(in.fLast),
125  fTrackCounter(0x0)
126 {
127  //cpy constructor
128 }
129
130 AliReader::~AliReader()
131 {
132 //destructor
133  if(fCuts)
134   {
135    fCuts->SetOwner();
136    delete fCuts;
137   }
138  delete fEventSim;
139  delete fEventRec;
140  delete fTrackCounter;
141 }
142 /*************************************************************************************/
143
144 AliReader& AliReader::operator=(const AliReader& in)
145 {
146   //Assigment operator
147  if (this == &in) return *this;  
148  TNamed::operator=( (const TNamed&)in );
149   
150  fCuts = (in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0;
151  fDirs = (in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0;
152  fCurrentEvent = 0;
153  fCurrentDir = 0;
154  fNEventsRead = 0;
155  fEventRec = 0x0;
156  fEventSim = 0x0;
157  fRunSim = 0x0;
158  fRunRec = 0x0;
159  fIsRead = kFALSE;
160  fBufferEvents = in.fBufferEvents;
161  fBlend = in.fBlend;
162  fFirst = in.fFirst;
163  fLast = in.fLast;
164  fTrackCounter = 0x0;
165  return *this;  
166 }
167 /*************************************************************************************/
168
169 Int_t AliReader::Next()
170 {
171 //moves to next event
172
173   //if asked to read up to event nb. fLast, and it is overcome, report no more events
174   if ((fNEventsRead >= fLast) && (fLast > 0) ) return kTRUE;
175   
176   if (fTrackCounter == 0x0)//create Track Counter
177    {
178      fTrackCounter = new TH1I("trackcounter","Track Counter",20000,0,20000);
179      fTrackCounter->SetDirectory(0x0);
180    }
181   
182   do //if asked to read from event fFirst, rewind to it
183    {
184     if ( ReadNext() == kTRUE) //if no more evets, return it
185       return kTRUE;
186    }while (fNEventsRead < fFirst);
187    
188   //here we have event
189   
190   if (fBlend) Blend();//Mix particles order 
191   
192   if (fBufferEvents)//store events if buffering is on
193    {
194      if ( ReadsRec() && fEventRec) 
195        fRunRec->SetEvent(fNEventsRead-1-fFirst,fEventRec);
196      if ( ReadsSim() && fEventSim)
197        fRunSim->SetEvent(fNEventsRead-1-fFirst,fEventSim);
198    }
199   return kFALSE;
200 }
201 /*************************************************************************************/
202
203 void AliReader::AddParticleCut(AliAODParticleCut* cut)
204 {
205  //sets the new cut. MAKES A COPY OF THE CUT !!!!
206  
207   if (!cut) //if cut is NULL return with error
208    {
209     Error("AddParticleType","NULL pointers are not accepted any more.\nIf You want to accept all particles of this type, set an empty cut ");
210     return;
211    }
212   AliAODParticleCut *c = (AliAODParticleCut*)cut->Clone();
213   fCuts->Add(c);
214 }
215 /********************************************************************/
216
217 AliAOD* AliReader::GetEventSim(Int_t n)
218  {
219  //returns Nth event with simulated particles
220   if (ReadsSim() == kFALSE)
221    {
222      Error("GetParticleEvent","This reader is not able to provide simulated particles.");
223      return 0;
224    } 
225    
226   if (!fIsRead) 
227    { 
228     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
229     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
230     
231     if (Read(fRunSim,fRunRec))
232      {
233        Error("GetParticleEvent","Error in reading");
234        return 0x0;
235      }
236     else fIsRead = kTRUE;
237    }
238   return fRunSim->GetEvent(n);
239  }
240 /********************************************************************/
241
242 AliAOD* AliReader::GetEventRec(Int_t n)
243  {
244  //returns Nth event with reconstructed tracks
245   if (ReadsRec() == kFALSE)
246    {
247      Error("GetTrackEvent","This reader is not able to provide recosntructed tracks.");
248      return 0;
249    } 
250   if (!fIsRead) 
251    {
252     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
253     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
254     
255     if(Read(fRunSim,fRunRec))
256      {
257        Error("GetTrackEvent","Error in reading");
258        return 0x0;
259      }
260     else fIsRead = kTRUE;
261    }
262   return fRunRec->GetEvent(n);
263  }
264 /********************************************************************/
265
266 Int_t AliReader::GetNumberOfSimEvents()
267  {
268  //returns number of events of particles
269   if (ReadsSim() == kFALSE)
270    {
271      Error("GetNumberOfPartEvents","This reader is not able to provide simulated particles.");
272      return 0;
273    } 
274    
275   if (!fIsRead) 
276    {
277     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
278     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
279     
280     if (Read(fRunSim,fRunRec))
281      {
282        Error("GetNumberOfPartEvents","Error in reading");
283        return 0;
284      }
285     else fIsRead = kTRUE;
286    }
287    return fRunSim->GetNumberOfEvents();
288  }
289 /********************************************************************/
290  
291 Int_t AliReader::GetNumberOfRecEvents()
292  {
293  //returns number of events of tracks
294   if (ReadsRec() == kFALSE)
295    {
296      Error("GetNumberOfTrackEvents","This reader is not able to provide recosntructed tracks.");
297      return 0;
298    } 
299   if (!fIsRead)
300    {
301      if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
302      if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
303      
304      if(Read(fRunSim,fRunRec))
305       {
306         Error("GetNumberOfTrackEvents","Error in reading");
307         return 0;
308       }
309      else fIsRead = kTRUE;
310    }
311   return fRunRec->GetNumberOfEvents();
312  }
313 /********************************************************************/
314
315 Int_t AliReader::Read(AliAODRun* particles, AliAODRun *tracks)
316 {
317  //reads data and puts put to the particles and tracks objects
318  //reurns 0 if everything is OK
319  //
320   Info("Read","");
321   
322   if ( ReadsSim() && (particles == 0x0) ) //check if an object is instatiated
323    {
324      Error("Read"," particles object must be instatiated before passing it to the reader");
325      return 1;
326    }
327   if ( ReadsRec() && (tracks == 0x0) )  //check if an object is instatiated
328    {
329      Error("Read"," tracks object must be instatiated before passing it to the reader");
330      return 1;
331    }
332    
333   if (ReadsSim()) particles->Reset();//clear runs == delete all old events
334   if (ReadsRec()) tracks->Reset();
335   
336   Rewind();
337   
338   Int_t i = 0;
339   while(Next() == kFALSE)
340    {
341      if (ReadsRec()) tracks->SetEvent(i,fEventRec);
342      if (ReadsSim()) particles->SetEvent(i,fEventSim);
343      i++;
344    }
345   return 0;
346 }      
347 /*************************************************************************************/
348
349 Bool_t AliReader::Rejected(AliVAODParticle* p)
350 {
351  //Method examines whether particle meets all cut and particle type criteria
352   
353    if(p==0x0)//of corse we not pass NULL pointers
354     {
355      Warning("Rejected()","No Pasaran! We never accept NULL pointers");
356      return kTRUE;
357     }
358    //if no particle is specified, we pass all particles
359    //excluding NULL pointers, of course
360   if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
361   for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
362    {
363      AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i));
364      if(!cut.Rejected(p)) return kFALSE;  //accepted
365    }
366    
367   return kTRUE;//not accepted
368 }
369 /*************************************************************************************/
370
371 Bool_t  AliReader::Rejected(Int_t pid)
372 {
373 //this method checks if any of existing cuts accepts this pid particles
374 //or any cuts accepts all particles
375
376  if(pid == 0)
377   return kTRUE;
378
379  if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
380   
381  for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)   
382    {
383      AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i));
384      //if some of cuts accepts all particles or some accepts particles of this type, accept
385      if ( (cut.GetPID() == 0) || (cut.GetPID() == pid) ) return kFALSE; 
386    }
387  return kTRUE;
388 }
389 /*************************************************************************************/
390
391 TString AliReader::GetDirName(Int_t entry)
392 {
393 //returns directory name of next one to read
394   TString  retval;//return value
395   if (fDirs ==  0x0)
396    { 
397      if (entry == 0)
398       {
399        retval = ".";
400        return retval;
401       }
402      else
403       {
404        return retval;
405       }  
406    }
407   
408   
409   if ( (entry >= fDirs->GetEntries()) || (entry < 0))//if out of bounds return empty string
410    {                                            //note that entry==0 is accepted even if array is empty (size=0)
411     if ( (fDirs->GetEntries() == 0) && (entry == 0) )
412       { 
413         retval = ".";
414         return retval;
415       }
416      if (AliVAODParticle::GetDebug() > 0)
417       {
418         Warning("GetDirName","Index %d out of bounds",entry);
419       }
420      return retval;
421    }
422
423
424   TClass *objclass = fDirs->At(entry)->IsA();
425   TClass *stringclass = TObjString::Class();
426
427   TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
428
429   if(dir == 0x0)
430    {
431      Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
432      return retval;
433    }
434   if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
435   retval = dir->String();
436   return retval;
437 }
438 /*************************************************************************************/
439
440 void AliReader::Blend()
441 {
442   //randomly change positions of the particles after reading
443   //is used to check if some distributions (of many particle properties) 
444   //depend on the order of particles
445   //(tracking gives particles Pt sorted)
446   Int_t npart = 0;
447
448   if (fEventSim ) 
449    {
450      npart = fEventSim->GetNumberOfParticles();
451     } 
452   else 
453     if (fEventRec ) 
454      {
455         npart = fEventRec->GetNumberOfParticles();
456      }
457     else
458      {
459        return;
460      }
461   for (Int_t i = 2;  i < npart; i++)
462    {
463      Int_t with = gRandom->Integer(i);
464 //     Info("Blend","%d %d",i, with);
465      if (fEventSim) fEventSim->SwapParticles(i,with);
466      if (fEventRec) fEventRec->SwapParticles(i,with);
467    }
468 }
469 /*************************************************************************************/
470
471 void AliReader::WriteTrackCounter() const
472 {
473  //writes the counter histogram
474  
475  if (fTrackCounter) fTrackCounter->Write(0,TObject::kOverwrite);
476  else 
477   {
478     Warning("WriteTrackCounter","Counter is NULL");
479   }
480 }