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 <TGliteXmlEventlist.h>
  60
  61 #include "AliAODParticleCut.h"
  62 #include "AliAOD.h"
  63 #include "AliAODRun.h"
  64
  65 ClassImp(AliReader)
  66 //pure virtual
  67
  68 /*************************************************************************************/
  69
  70 AliReader::AliReader():
  71  fEventList(0x0),
  72  fCuts(new TObjArray()),
  73  fDirs(0x0),
  74  fCurrentEvent(0),
  75  fCurrentDir(0),
  76  fNEventsRead(0),
  77  fEventRec(0x0),
  78  fEventSim(0x0),
  79  fRunSim(0x0),
  80  fRunRec(0x0),
  81  fIsRead(kFALSE),
  82  fBufferEvents(kFALSE),
  83  fBlend(kFALSE),
  84  fFirst(0),
  85  fLast(0),
  86  fTrackCounter(0x0)
  87 {
  88 //constructor
  89 }
  90 /*************************************************************************************/
  91
  92 AliReader::AliReader(TObjArray* dirs):
  93  fEventList(0x0),
  94  fCuts(new TObjArray()),
  95  fDirs(dirs),
  96  fCurrentEvent(0),
  97  fCurrentDir(0),
  98  fNEventsRead(0),
  99  fEventRec(0x0),
 100  fEventSim(0x0),
 101  fRunSim(0x0),
 102  fRunRec(0x0),
 103  fIsRead(kFALSE),
 104  fBufferEvents(kFALSE),
 105  fBlend(kFALSE),
 106  fFirst(0),
 107  fLast(0),
 108  fTrackCounter(0x0)
 109 {
 110 //ctor with array of directories to read as parameter
 111 }
 112 /*************************************************************************************/
 113 AliReader::AliReader(const AliReader& in):
 114  TNamed(in),
 115  fEventList((in.fEventList)?(TGliteXmlEventlist*)in.fEventList->Clone():0x0),
 116  fCuts((in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0),
 117  fDirs((in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0),
 118  fCurrentEvent(0),
 119  fCurrentDir(0),
 120  fNEventsRead(0),
 121  fEventRec(0x0),
 122  fEventSim(0x0),
 123  fRunSim(0x0),
 124  fRunRec(0x0),
 125  fIsRead(kFALSE),
 126  fBufferEvents(in.fBufferEvents),
 127  fBlend(in.fBlend),
 128  fFirst(in.fFirst),
 129  fLast(in.fLast),
 130  fTrackCounter(0x0)
 131 {
 132  //cpy constructor
 133 }
 134
 135 AliReader::~AliReader()
 136 {
 137 //destructor
 138  if(fCuts)
 139   {
 140    fCuts->SetOwner();
 141    delete fCuts;
 142   }
 143  delete fEventSim;
 144  delete fEventRec;
 145  delete fTrackCounter;
 146  delete fEventList;
 147 }
 148 /*************************************************************************************/
 149
 150 AliReader& AliReader::operator=(const AliReader& in)
 151 {
 152   //Assigment operator
 153  if (this == &in) return *this;
 154  TNamed::operator=( (const TNamed&)in );
 155
 156  fCuts = (in.fCuts)?(TObjArray*)in.fCuts->Clone():0x0;
 157  fDirs = (in.fDirs)?(TObjArray*)in.fDirs->Clone():0x0;
 158  fCurrentEvent = 0;
 159  fCurrentDir = 0;
 160  fNEventsRead = 0;
 161  fEventRec = 0x0;
 162  fEventSim = 0x0;
 163  fRunSim = 0x0;
 164  fRunRec = 0x0;
 165  fIsRead = kFALSE;
 166  fBufferEvents = in.fBufferEvents;
 167  fBlend = in.fBlend;
 168  fFirst = in.fFirst;
 169  fLast = in.fLast;
 170  fTrackCounter = 0x0;
 171  return *this;
 172 }
 173 /*************************************************************************************/
 174
 175 Int_t AliReader::Next()
 176 {
 177 //moves to next event
 178
 179   //if asked to read up to event nb. fLast, and it is overcome, report no more events
 180   if ((fNEventsRead >= fLast) && (fLast > 0) ) return kTRUE;
 181
 182   if (fTrackCounter == 0x0)//create Track Counter
 183    {
 184      fTrackCounter = new TH1I("trackcounter","Track Counter",20000,0,20000);
 185      fTrackCounter->SetDirectory(0x0);
 186    }
 187
 188   do //if asked to read from event fFirst, rewind to it
 189    {
 190     if ( ReadNext() == kTRUE) //if no more evets, return it
 191       return kTRUE;
 192    }while (fNEventsRead < fFirst);
 193
 194   //here we have event
 195
 196   if (fBlend) Blend();//Mix particles order
 197
 198   if (fBufferEvents)//store events if buffering is on
 199    {
 200      if ( ReadsRec() && fEventRec)
 201        fRunRec->SetEvent(fNEventsRead-1-fFirst,fEventRec);
 202      if ( ReadsSim() && fEventSim)
 203        fRunSim->SetEvent(fNEventsRead-1-fFirst,fEventSim);
 204    }
 205   return kFALSE;
 206 }
 207 /*************************************************************************************/
 208
 209 void AliReader::AddParticleCut(AliAODParticleCut* cut)
 210 {
 211  //sets the new cut. MAKES A COPY OF THE CUT !!!!
 212
 213   if (!cut) //if cut is NULL return with error
 214    {
 215     Error("AddParticleType","NULL pointers are not accepted any more.\nIf You want to accept all particles of this type, set an empty cut ");
 216     return;
 217    }
 218   AliAODParticleCut *c = (AliAODParticleCut*)cut->Clone();
 219   fCuts->Add(c);
 220 }
 221 /********************************************************************/
 222
 223 AliAOD* AliReader::GetEventSim(Int_t n)
 224  {
 225  //returns Nth event with simulated particles
 226   if (ReadsSim() == kFALSE)
 227    {
 228      Error("GetParticleEvent","This reader is not able to provide simulated particles.");
 229      return 0;
 230    }
 231
 232   if (!fIsRead)
 233    {
 234     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
 235     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
 236
 237     if (Read(fRunSim,fRunRec))
 238      {
 239        Error("GetParticleEvent","Error in reading");
 240        return 0x0;
 241      }
 242     else fIsRead = kTRUE;
 243    }
 244   return fRunSim->GetEvent(n);
 245  }
 246 /********************************************************************/
 247
 248 AliAOD* AliReader::GetEventRec(Int_t n)
 249  {
 250  //returns Nth event with reconstructed tracks
 251   if (ReadsRec() == kFALSE)
 252    {
 253      Error("GetTrackEvent","This reader is not able to provide recosntructed tracks.");
 254      return 0;
 255    }
 256   if (!fIsRead)
 257    {
 258     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
 259     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
 260
 261     if(Read(fRunSim,fRunRec))
 262      {
 263        Error("GetTrackEvent","Error in reading");
 264        return 0x0;
 265      }
 266     else fIsRead = kTRUE;
 267    }
 268   return fRunRec->GetEvent(n);
 269  }
 270 /********************************************************************/
 271
 272 Int_t AliReader::GetNumberOfSimEvents()
 273  {
 274  //returns number of events of particles
 275   if (ReadsSim() == kFALSE)
 276    {
 277      Error("GetNumberOfPartEvents","This reader is not able to provide simulated particles.");
 278      return 0;
 279    }
 280
 281   if (!fIsRead)
 282    {
 283     if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
 284     if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
 285
 286     if (Read(fRunSim,fRunRec))
 287      {
 288        Error("GetNumberOfPartEvents","Error in reading");
 289        return 0;
 290      }
 291     else fIsRead = kTRUE;
 292    }
 293    return fRunSim->GetNumberOfEvents();
 294  }
 295 /********************************************************************/
 296
 297 Int_t AliReader::GetNumberOfRecEvents()
 298  {
 299  //returns number of events of tracks
 300   if (ReadsRec() == kFALSE)
 301    {
 302      Error("GetNumberOfTrackEvents","This reader is not able to provide recosntructed tracks.");
 303      return 0;
 304    }
 305   if (!fIsRead)
 306    {
 307      if (ReadsSim() && (fRunSim == 0x0)) fRunSim = new AliAODRun();
 308      if (ReadsRec() && (fRunRec == 0x0)) fRunRec = new AliAODRun();
 309
 310      if(Read(fRunSim,fRunRec))
 311       {
 312         Error("GetNumberOfTrackEvents","Error in reading");
 313         return 0;
 314       }
 315      else fIsRead = kTRUE;
 316    }
 317   return fRunRec->GetNumberOfEvents();
 318  }
 319 /********************************************************************/
 320
 321 Int_t AliReader::Read(AliAODRun* particles, AliAODRun *tracks)
 322 {
 323  //reads data and puts put to the particles and tracks objects
 324  //reurns 0 if everything is OK
 325  //
 326   Info("Read","");
 327
 328   if ( ReadsSim() && (particles == 0x0) ) //check if an object is instatiated
 329    {
 330      Error("Read"," particles object must be instatiated before passing it to the reader");
 331      return 1;
 332    }
 333   if ( ReadsRec() && (tracks == 0x0) )  //check if an object is instatiated
 334    {
 335      Error("Read"," tracks object must be instatiated before passing it to the reader");
 336      return 1;
 337    }
 338
 339   if (ReadsSim()) particles->Reset();//clear runs == delete all old events
 340   if (ReadsRec()) tracks->Reset();
 341
 342   Rewind();
 343
 344   Int_t i = 0;
 345   while(Next() == kFALSE)
 346    {
 347      if (ReadsRec()) tracks->SetEvent(i,fEventRec);
 348      if (ReadsSim()) particles->SetEvent(i,fEventSim);
 349      i++;
 350    }
 351   return 0;
 352 }
 353 /*************************************************************************************/
 354
 355 Bool_t AliReader::Rejected(AliVAODParticle* p)
 356 {
 357  //Method examines whether particle meets all cut and particle type criteria
 358
 359    if(p==0x0)//of corse we not pass NULL pointers
 360     {
 361      Warning("Rejected()","No Pasaran! We never accept NULL pointers");
 362      return kTRUE;
 363     }
 364    //if no particle is specified, we pass all particles
 365    //excluding NULL pointers, of course
 366   if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
 367   for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)
 368    {
 369      AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i));
 370      if(!cut.Rejected(p)) return kFALSE;  //accepted
 371    }
 372
 373   return kTRUE;//not accepted
 374 }
 375 /*************************************************************************************/
 376
 377 Bool_t  AliReader::Rejected(Int_t pid)
 378 {
 379 //this method checks if any of existing cuts accepts this pid particles
 380 //or any cuts accepts all particles
 381
 382  if(pid == 0)
 383   return kTRUE;
 384
 385  if ( fCuts->GetEntriesFast() == 0 ) return kFALSE; //if no cut specified accept all particles
 386
 387  for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)
 388    {
 389      AliAODParticleCut &cut = *((AliAODParticleCut*)fCuts->At(i));
 390      //if some of cuts accepts all particles or some accepts particles of this type, accept
 391      if ( (cut.GetPID() == 0) || (cut.GetPID() == pid) ) return kFALSE;
 392    }
 393  return kTRUE;
 394 }
 395 /*************************************************************************************/
 396
 397 TString AliReader::GetDirName(Int_t entry)
 398 {
 399 //returns directory name of next one to read
 400   TString  retval;//return value
 401   if (fDirs ==  0x0)
 402    {
 403      if (entry == 0)
 404       {
 405        retval = ".";
 406        return retval;
 407       }
 408      else
 409       {
 410        return retval;
 411       }
 412    }
 413
 414
 415   if ( (entry >= fDirs->GetEntries()) || (entry < 0))//if out of bounds return empty string
 416    {                                            //note that entry==0 is accepted even if array is empty (size=0)
 417     if ( (fDirs->GetEntries() == 0) && (entry == 0) )
 418       {
 419         retval = ".";
 420         return retval;
 421       }
 422      if (AliVAODParticle::GetDebug() > 0)
 423       {
 424         Warning("GetDirName","Index %d out of bounds",entry);
 425       }
 426      return retval;
 427    }
 428
 429
 430   TClass *objclass = fDirs->At(entry)->IsA();
 431   TClass *stringclass = TObjString::Class();
 432
 433   TObjString *dir = (TObjString*)objclass->DynamicCast(stringclass,fDirs->At(entry));
 434
 435   if(dir == 0x0)
 436    {
 437      Error("GetDirName","Object in TObjArray is not a TObjString or its descendant");
 438      return retval;
 439    }
 440   if (gDebug > 0) Info("GetDirName","Returned ok %s",dir->String().Data());
 441   retval = dir->String();
 442   return retval;
 443 }
 444 /*************************************************************************************/
 445
 446 void AliReader::Blend()
 447 {
 448   //randomly change positions of the particles after reading
 449   //is used to check if some distributions (of many particle properties)
 450   //depend on the order of particles
 451   //(tracking gives particles Pt sorted)
 452   Int_t npart = 0;
 453
 454   if (fEventSim )
 455    {
 456      npart = fEventSim->GetNumberOfParticles();
 457     }
 458   else
 459     if (fEventRec )
 460      {
 461         npart = fEventRec->GetNumberOfParticles();
 462      }
 463     else
 464      {
 465        return;
 466      }
 467   for (Int_t i = 2;  i < npart; i++)
 468    {
 469      Int_t with = gRandom->Integer(i);
 470 //     Info("Blend","%d %d",i, with);
 471      if (fEventSim) fEventSim->SwapParticles(i,with);
 472      if (fEventRec) fEventRec->SwapParticles(i,with);
 473    }
 474 }
 475 /*************************************************************************************/
 476
 477 void AliReader::WriteTrackCounter() const
 478 {
 479  //writes the counter histogram
 480
 481  if (fTrackCounter) fTrackCounter->Write(0,TObject::kOverwrite);
 482  else
 483   {
 484     Warning("WriteTrackCounter","Counter is NULL");
 485   }
 486 }