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 }