ANALYSIS/AliReaderAOD.cxx

   1 #include "AliReaderAOD.h"
   2 //______________________________________________________________________________
   3 ////////////////////////////////////////////////////////////////////////////////
   4 //                                                                            //
   5 // class AliReaderAOD                                                         //
   6 //                                                                            //
   7 // Reader and Writer for AOD format.                                          //
   8 // AODs are stored in a tree named by the variable fgkTreeName.               //
   9 // There is stored 1 or 2 branches. Each of them stores AOD objects           //
  10 // First branch is named by the variable fgkReconstructedDataBranchName       //
  11 // ("reconstructed.") and keeps reconstructed data.                           //
  12 // Second branch is called by the variable fgkSimulatedDataBranchName         //
  13 // ("simulated.") and stores Monte carlo truth. If both branches are present  //
  14 // AODs are parallel, i.e. nth particle in one branch corresponds to the nth  //
  15 // particle in the other one.                                                 //
  16 //                                                                            //
  17 // Since we accept different formats of particles that are stored in AODs     //
  18 // reader must take care of that fact: clean buffer if the next file contains //
  19 // different particle type.                                                   //
  20 //                                                                            //
  21 // Piotr.Skowronski@cern.ch                                                   //
  22 //                                                                            //
  23 ////////////////////////////////////////////////////////////////////////////////
  24
  25 ClassImp(AliReaderAOD)
  26
  27 #include <TError.h>
  28 #include <TFile.h>
  29 #include <TTree.h>
  30 #include <TH1.h>
  31 #include "AliAOD.h"
  32
  33
  34 const TString AliReaderAOD::fgkTreeName("TAOD");
  35 const TString AliReaderAOD::fgkReconstructedDataBranchName("reconstructed.");
  36 const TString AliReaderAOD::fgkSimulatedDataBranchName("simulated.");
  37
  38 AliReaderAOD::AliReaderAOD(const Char_t* aodfilename):
  39  fFileName(aodfilename),
  40  fReadSim(kFALSE),
  41  fReadRec(kTRUE),
  42  fTree(0x0),
  43  fFile(0x0),
  44  fSimBuffer(0x0),
  45  fRecBuffer(0x0)
  46 {
  47   //ctor
  48 }
  49 /********************************************************************/
  50
  51 AliReaderAOD::~AliReaderAOD()
  52 {
  53 //dtor
  54   if (fEventSim == fSimBuffer )
  55    {
  56     fEventSim = 0x0;
  57     fEventRec = 0x0;
  58    }
  59   delete fSimBuffer;
  60   delete fRecBuffer;
  61
  62   delete fTree;
  63   delete fFile;
  64 }
  65 /********************************************************************/
  66
  67 void AliReaderAOD::Rewind()
  68 {
  69 //Rewinds reading
  70   delete fTree;
  71   fTree = 0x0;
  72   delete fFile;
  73   fFile = 0x0;
  74   fCurrentDir = 0;
  75   fNEventsRead= 0;
  76 }
  77 /********************************************************************/
  78 Int_t AliReaderAOD::ReadNext()
  79 {
  80 //Reads next event
  81
  82   Info("ReadNext","Entered");
  83   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
  84     {
  85       if (fFile == 0x0)
  86        {
  87          Int_t openfailed = OpenFile(fCurrentDir);//rl is opened here
  88          if (openfailed)
  89            {
  90              //Error("ReadNext","Error Occured while opening directory number %d",fCurrentDir);
  91              fCurrentDir++;
  92              continue;
  93            }
  94          fCurrentEvent = 0;
  95        }
  96       //Tree must exist because OpenFile would reuturn error in the other case
  97       if ( fCurrentEvent >= fTree->GetEntries() )
  98        {
  99          delete fTree;
 100          fTree = 0x0;
 101          delete fFile;
 102          fFile = 0x0;
 103          fSimBuffer = 0x0;
 104          fRecBuffer = 0x0;
 105          fCurrentDir++;
 106          continue;
 107        }
 108
 109       Info("ReadNext","Getting event %d",fCurrentEvent);
 110       fTree->GetEvent(fCurrentEvent);
 111       Info("ReadNext","Getting event %d Done",fCurrentEvent);
 112
 113       Int_t retval = 0;
 114       if (fReadRec && fReadSim)
 115        {
 116          retval = ReadRecAndSim();
 117        }
 118       else
 119        {
 120          if (fReadRec) retval = ReadRec();
 121          if (fReadSim) retval = ReadSim();
 122        }
 123
 124       fCurrentEvent++;
 125       if (retval == 0) fNEventsRead++;
 126
 127       return retval;//success -> read one event
 128
 129     }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 130
 131   return 1; //no more directories to read
 132
 133
 134 }
 135 /********************************************************************/
 136
 137 Int_t AliReaderAOD::ReadRecAndSim()
 138 {
 139 //Reads raconstructed and simulated data
 140
 141  Info("ReadRecAndSim","Found %d reconstructed tracks and %d simulated particles",
 142        fRecBuffer->GetNumberOfParticles(),fSimBuffer->GetNumberOfParticles());
 143
 144  if (fCuts->GetEntriesFast() == 0x0)
 145   {//if there is no cuts we return pointer to the buffer
 146     if (fEventRec != fRecBuffer)
 147      {
 148        delete fEventRec;
 149        delete fEventSim;
 150      }
 151     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
 152     fEventSim = fSimBuffer;
 153   }
 154  else
 155   {//if there are cuts specified
 156     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
 157      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 158        fEventRec = new AliAOD();
 159        fEventSim = new AliAOD();
 160
 161        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
 162        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
 163      }
 164     else
 165      {//or simply reset it in case it already exists
 166        fEventRec->Reset();
 167        fEventSim->Reset();
 168      }
 169
 170     Int_t npart = fRecBuffer->GetNumberOfParticles();
 171     for (Int_t i = 0; i < npart; i++)
 172      {
 173        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
 174        if (Rejected(prec)) continue;//we make cuts only on simulated data
 175
 176        fEventRec->AddParticle(prec);
 177        fEventSim->AddParticle( fSimBuffer->GetParticle(i));
 178      }
 179   }
 180
 181  Info("ReadRecAndSim","Read %d reconstructed tracks and %d simulated particles",
 182        fEventRec->GetNumberOfParticles(),fEventSim->GetNumberOfParticles());
 183
 184  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
 185
 186  return 0;
 187 }
 188 /********************************************************************/
 189
 190 Int_t AliReaderAOD::ReadRec()
 191 {
 192 //Reads reconstructed data only
 193
 194  Info("ReadRec","Found %d reconstructed tracks",fRecBuffer->GetNumberOfParticles());
 195
 196  if (fCuts->GetEntriesFast() == 0x0)
 197   {//if there is no cuts we return pointer to the buffer
 198     if (fEventRec != fRecBuffer)
 199      {
 200        delete fEventRec;
 201      }
 202     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
 203   }
 204  else
 205   {//if there are cuts specified
 206     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
 207      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 208        fEventRec = new AliAOD();
 209
 210        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
 211      }
 212     else
 213      {//or simply reset it in case it already exists
 214        fEventRec->Reset();
 215      }
 216
 217     Int_t npart = fRecBuffer->GetNumberOfParticles();
 218     for (Int_t i = 0; i < npart; i++)
 219      {
 220        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
 221        if (Rejected(prec)) continue;//we make cuts only on simulated data
 222
 223        fEventRec->AddParticle(prec);
 224      }
 225   }
 226
 227  Info("ReadRec","Read %d reconstructed tracks",fEventRec->GetNumberOfParticles());
 228  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
 229
 230  return 0;
 231 }
 232 /********************************************************************/
 233
 234 Int_t AliReaderAOD::ReadSim()
 235 {
 236 //Reads simulated data only
 237
 238  Info("ReadSim","Found %d simulated particles",fSimBuffer->GetNumberOfParticles());
 239
 240  if (fCuts->GetEntriesFast() == 0x0)
 241   {//if there is no cuts we return pointer to the buffer
 242     if (fEventSim != fSimBuffer)
 243      {
 244        delete fEventSim;
 245      }
 246     fEventSim = fSimBuffer;
 247   }
 248  else
 249   {//if there are cuts specified
 250     if ( (fEventSim == 0x0) || (fEventSim == fSimBuffer) )
 251      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 252        fEventSim = new AliAOD();
 253
 254        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
 255      }
 256     else
 257      {//or simply reset it in case it already exists
 258        fEventSim->Reset();
 259      }
 260
 261     Int_t npart = fSimBuffer->GetNumberOfParticles();
 262     for (Int_t i = 0; i < npart; i++)
 263      {
 264        AliVAODParticle* prec = fSimBuffer->GetParticle(i);
 265        if (Rejected(prec)) continue;//we make cuts only on simulated data
 266        fEventSim->AddParticle(prec);
 267      }
 268   }
 269
 270  Info("ReadSim","Read %d simulated particles",fEventSim->GetNumberOfParticles());
 271  fTrackCounter->Fill(fEventSim->GetNumberOfParticles());
 272
 273
 274  return 0;
 275 }
 276 /********************************************************************/
 277
 278 Int_t AliReaderAOD::OpenFile(Int_t n)
 279 {
 280 //opens fFile with  tree
 281
 282 // Info("ReadNext","Opening File %d",n);
 283  const TString dirname = GetDirName(n);
 284  if (dirname == "")
 285   {
 286     if (AliVAODParticle::GetDebug() > 2 )
 287       {
 288         Info("OpenFile","Got empty string as a directory name.");
 289       }
 290    return 1;
 291   }
 292
 293  TString filename = dirname +"/"+ fFileName;
 294  fFile = TFile::Open(filename.Data());
 295  if ( fFile == 0x0)
 296   {
 297     Error("OpenFile","Can't open fFile %s",filename.Data());
 298     return 2;
 299   }
 300  if (!fFile->IsOpen())
 301   {
 302     Error("OpenFile","Can't open fFile  %s",filename.Data());
 303     delete fFile;
 304     fFile = 0x0;
 305     return 3;
 306   }
 307
 308  Info("ReadNext","File Is Opened, Getting the TREE");
 309
 310  fTree = dynamic_cast<TTree*>(fFile->Get(fgkTreeName));
 311  if (fTree == 0x0)
 312   {
 313     if (AliVAODParticle::GetDebug() > 2 )
 314       {
 315         Info("ReadNext","Can not find TTree object named %s",fgkTreeName.Data());
 316       }
 317     delete fFile;
 318     fFile = 0x0;
 319     return 4;
 320   }
 321
 322 //  Info("ReadNext","Got TREE, Setting branch addresses");
 323
 324   if (fReadRec)
 325    {
 326      TBranch* branch = fTree->GetBranch(fgkReconstructedDataBranchName);
 327      if (branch == 0x0)
 328       {
 329         Error("OpenFile","Can not find branch %s in file %s",
 330               fgkReconstructedDataBranchName.Data(),filename.Data());
 331
 332         delete fTree;
 333         fTree = 0x0;
 334         delete fFile;
 335         fFile = 0x0;
 336         return 5;
 337       }
 338      fTree->SetBranchAddress(fgkReconstructedDataBranchName,&fRecBuffer);
 339    }
 340
 341
 342   if (fReadSim)
 343    {
 344      TBranch* branch = fTree->GetBranch(fgkSimulatedDataBranchName);
 345      if (branch == 0x0)
 346       {
 347         Error("OpenFile","Can not find branch %s in file %s",
 348               fgkSimulatedDataBranchName.Data(),filename.Data());
 349
 350         delete fTree;
 351         fTree = 0x0;
 352         delete fFile;
 353         fFile = 0x0;
 354         return 6;
 355       }
 356      fTree->SetBranchAddress(fgkSimulatedDataBranchName,&fSimBuffer);
 357    }
 358 //  Info("ReadNext","Got TREE, Addresses are set.");
 359 //  Info("ReadNext","Quitting the method.");
 360
 361   return 0;
 362
 363 }
 364 /********************************************************************/
 365
 366 Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, const char* pclassname,  Bool_t /*multcheck*/)
 367 {
 368 //reads tracks from runs and writes them to file
 369   ::Info("AliReaderAOD::Write","________________________________________________________");
 370   ::Info("AliReaderAOD::Write","________________________________________________________");
 371   ::Info("AliReaderAOD::Write","________________________________________________________");
 372
 373   if (reader == 0x0)
 374    {
 375      ::Error("AliReaderAOD::Write","Input Reader is NULL");
 376      return -1;
 377    }
 378   TFile *outfile = TFile::Open(outfilename,"recreate");
 379   if (outfile == 0x0)
 380    {
 381      ::Error("AliReaderAOD::Write","Can not open output file %s",outfilename);
 382      return -1;
 383    }
 384
 385   TTree *tree = new TTree(fgkTreeName,"Tree with tracks");
 386
 387   TBranch *recbranch = 0x0, *simbranch = 0x0;
 388
 389   AliAOD* eventrec = new AliAOD();//must be created before Branch is called. Otherwise clones array is not splitted
 390   AliAOD* eventsim = new AliAOD();//AOD together with fParticles clones array knowing exact type of particles
 391
 392   eventrec->SetParticleClassName(pclassname);
 393   eventsim->SetParticleClassName(pclassname);
 394
 395   AliAOD* recbuffer = eventrec;
 396   AliAOD* simbuffer = eventsim;
 397
 398   if (reader->ReadsRec()) recbranch = tree->Branch(fgkReconstructedDataBranchName,"AliAOD",&recbuffer,32000,99);
 399   if (reader->ReadsSim()) simbranch = tree->Branch(fgkSimulatedDataBranchName,"AliAOD",&simbuffer,32000,99);
 400
 401   reader->Rewind();
 402   while (reader->Next() == kFALSE)
 403    {
 404
 405      if (reader->ReadsRec())
 406       {//here we can get AOD that has different particle type
 407         AliAOD* event = reader->GetEventRec();
 408         if ( eventrec->GetParticleClass() != event->GetParticleClass() )
 409          {//if class type is not what what we whant we copy particles
 410            eventrec->CopyData(event);
 411            recbuffer = eventrec;
 412          }
 413         else
 414          {//else just pointer to event from input reader is passed
 415            recbuffer = event;
 416          }
 417         recbuffer->GetParticle(0)->Print();
 418       }
 419
 420      if (reader->ReadsSim())
 421       {
 422         AliAOD* event = reader->GetEventSim();
 423         if ( eventsim->GetParticleClass() != event->GetParticleClass() )
 424          {//if class type is not what what we whant we copy particles
 425            eventsim->CopyData(event);
 426            simbuffer = eventrec;
 427          }
 428         else
 429          {//else just pointer to event from input reader is passed
 430            simbuffer = event;
 431          }
 432         simbuffer->GetParticle(0)->Print();
 433       }
 434      tree->Fill();
 435    }
 436
 437   ::Info("AliReaderAOD::Write","Written %d events",tree->GetEntries());
 438   outfile->cd();
 439   tree->Write();
 440
 441   delete eventsim;
 442   delete eventrec;
 443
 444   delete tree;
 445   delete outfile;
 446   return 0;
 447 }
 448