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
 104          delete fSimBuffer;
 105          delete fRecBuffer;
 106
 107          fSimBuffer = 0x0;
 108          fRecBuffer = 0x0;
 109          fCurrentDir++;
 110          continue;
 111        }
 112
 113       Info("ReadNext","Getting event %d",fCurrentEvent);
 114       fTree->GetEvent(fCurrentEvent);
 115       Info("ReadNext","Getting event %d Done",fCurrentEvent);
 116
 117       Int_t retval = 0;
 118       if (fReadRec && fReadSim)
 119        {
 120          retval = ReadRecAndSim();
 121        }
 122       else
 123        {
 124          if (fReadRec) retval = ReadRec();
 125          if (fReadSim) retval = ReadSim();
 126        }
 127
 128       fCurrentEvent++;
 129       if (retval != 0)
 130         {
 131           //something wrong has happend during reading this event, take next
 132           continue;
 133         }
 134
 135       fNEventsRead++;
 136       return retval;//success -> read one event
 137
 138     }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 139
 140   return 1; //no more directories to read
 141
 142
 143 }
 144 /********************************************************************/
 145
 146 Int_t AliReaderAOD::ReadRecAndSim()
 147 {
 148 //Reads raconstructed and simulated data
 149
 150  Info("ReadRecAndSim","Found %d reconstructed tracks and %d simulated particles",
 151        fRecBuffer->GetNumberOfParticles(),fSimBuffer->GetNumberOfParticles());
 152
 153  if (fCuts->GetEntriesFast() == 0x0)
 154   {//if there is no cuts we return pointer to the buffer
 155     if (fEventRec != fRecBuffer)
 156      {
 157        delete fEventRec;
 158        delete fEventSim;
 159      }
 160     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
 161     fEventSim = fSimBuffer;
 162   }
 163  else
 164   {//if there are cuts specified
 165     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
 166      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 167        fEventRec = new AliAOD();
 168        fEventSim = new AliAOD();
 169
 170        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
 171        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
 172      }
 173     else
 174      {//or simply reset it in case it already exists
 175        fEventRec->Reset();
 176        fEventSim->Reset();
 177      }
 178
 179     Int_t npart = fRecBuffer->GetNumberOfParticles();
 180
 181     if (npart != fSimBuffer->GetNumberOfParticles())
 182      {
 183        Error("ReadRecAndSim","There is different number of simulated and reconstructed particles!",
 184                               fSimBuffer->GetNumberOfParticles(),npart);
 185        return 1;
 186      }
 187     for (Int_t i = 0; i < npart; i++)
 188      {
 189        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
 190        AliVAODParticle* psim = fSimBuffer->GetParticle(i);
 191
 192        if (prec == 0x0)
 193         {
 194           Error("ReadRecAndSim","Reconstructed Particle is NULL !!!");
 195           continue;
 196         }
 197        if (psim == 0x0)
 198         {
 199           Error("ReadRecAndSim","Simulated Particle is NULL !!!");
 200           continue;
 201         }
 202
 203        if (Rejected(prec)) continue;//we make cuts only on reconstructed data
 204
 205        fEventRec->AddParticle(prec);
 206        fEventSim->AddParticle( fSimBuffer->GetParticle(i));
 207      }
 208   }
 209
 210  Info("ReadRecAndSim","Read %d reconstructed tracks and %d simulated particles",
 211        fEventRec->GetNumberOfParticles(),fEventSim->GetNumberOfParticles());
 212
 213  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
 214
 215  return 0;
 216 }
 217 /********************************************************************/
 218
 219 Int_t AliReaderAOD::ReadRec()
 220 {
 221 //Reads reconstructed data only
 222
 223  Info("ReadRec","Found %d reconstructed tracks",fRecBuffer->GetNumberOfParticles());
 224
 225  if (fCuts->GetEntriesFast() == 0x0)
 226   {//if there is no cuts we return pointer to the buffer
 227     if (fEventRec != fRecBuffer)
 228      {
 229        delete fEventRec;
 230      }
 231     fEventRec = fRecBuffer;//fEventRec is the pointer that the user gets when he asks about an event
 232   }
 233  else
 234   {//if there are cuts specified
 235     if ( (fEventRec == 0x0) || (fEventRec == fRecBuffer) )
 236      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 237        fEventRec = new AliAOD();
 238
 239        fEventRec->SetParticleClass( fRecBuffer->GetParticleClass() );
 240      }
 241     else
 242      {//or simply reset it in case it already exists
 243        fEventRec->Reset();
 244      }
 245
 246     Int_t npart = fRecBuffer->GetNumberOfParticles();
 247     for (Int_t i = 0; i < npart; i++)
 248      {
 249        AliVAODParticle* prec = fRecBuffer->GetParticle(i);
 250        if (Rejected(prec)) continue;//we make cuts only on simulated data
 251
 252        fEventRec->AddParticle(prec);
 253      }
 254   }
 255
 256  Info("ReadRec","Read %d reconstructed tracks",fEventRec->GetNumberOfParticles());
 257  fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
 258
 259  return 0;
 260 }
 261 /********************************************************************/
 262
 263 Int_t AliReaderAOD::ReadSim()
 264 {
 265 //Reads simulated data only
 266
 267  Info("ReadSim","Found %d simulated particles",fSimBuffer->GetNumberOfParticles());
 268
 269  if (fCuts->GetEntriesFast() == 0x0)
 270   {//if there is no cuts we return pointer to the buffer
 271     if (fEventSim != fSimBuffer)
 272      {
 273        delete fEventSim;
 274      }
 275     fEventSim = fSimBuffer;
 276   }
 277  else
 278   {//if there are cuts specified
 279     if ( (fEventSim == 0x0) || (fEventSim == fSimBuffer) )
 280      {//we need to create a new event, if it is not existing  or it is the same as branch buffer
 281        fEventSim = new AliAOD();
 282
 283        fEventSim->SetParticleClass( fSimBuffer->GetParticleClass() );
 284      }
 285     else
 286      {//or simply reset it in case it already exists
 287        fEventSim->Reset();
 288      }
 289
 290     Int_t npart = fSimBuffer->GetNumberOfParticles();
 291     for (Int_t i = 0; i < npart; i++)
 292      {
 293        AliVAODParticle* prec = fSimBuffer->GetParticle(i);
 294        if (Rejected(prec)) continue;//we make cuts only on simulated data
 295        fEventSim->AddParticle(prec);
 296      }
 297   }
 298
 299  Info("ReadSim","Read %d simulated particles",fEventSim->GetNumberOfParticles());
 300  fTrackCounter->Fill(fEventSim->GetNumberOfParticles());
 301
 302
 303  return 0;
 304 }
 305 /********************************************************************/
 306
 307 Int_t AliReaderAOD::OpenFile(Int_t n)
 308 {
 309 //opens fFile with  tree
 310
 311 // Info("ReadNext","Opening File %d",n);
 312  const TString dirname = GetDirName(n);
 313  if (dirname == "")
 314   {
 315     if (AliVAODParticle::GetDebug() > 2 )
 316       {
 317         Info("OpenFile","Got empty string as a directory name.");
 318       }
 319    return 1;
 320   }
 321
 322  TString filename = dirname +"/"+ fFileName;
 323  fFile = TFile::Open(filename.Data());
 324  if ( fFile == 0x0)
 325   {
 326     Error("OpenFile","Can't open fFile %s",filename.Data());
 327     return 2;
 328   }
 329  if (!fFile->IsOpen())
 330   {
 331     Error("OpenFile","Can't open fFile  %s",filename.Data());
 332     delete fFile;
 333     fFile = 0x0;
 334     return 3;
 335   }
 336
 337  Info("ReadNext","File %s Is Opened, Getting the TREE",filename.Data());
 338
 339  fTree = dynamic_cast<TTree*>(fFile->Get(fgkTreeName));
 340  if (fTree == 0x0)
 341   {
 342     if (AliVAODParticle::GetDebug() > 2 )
 343       {
 344         Info("ReadNext","Can not find TTree object named %s",fgkTreeName.Data());
 345       }
 346     delete fFile;
 347     fFile = 0x0;
 348     return 4;
 349   }
 350
 351 //  Info("ReadNext","Got TREE, Setting branch addresses");
 352
 353   if (fReadRec)
 354    {
 355      TBranch* branch = fTree->GetBranch(fgkReconstructedDataBranchName);
 356      if (branch == 0x0)
 357       {
 358         Error("OpenFile","Can not find branch %s in file %s",
 359               fgkReconstructedDataBranchName.Data(),filename.Data());
 360
 361         delete fTree;
 362         fTree = 0x0;
 363         delete fFile;
 364         fFile = 0x0;
 365         return 5;
 366       }
 367      fTree->SetBranchAddress(fgkReconstructedDataBranchName,&fRecBuffer);
 368    }
 369
 370
 371   if (fReadSim)
 372    {
 373      TBranch* branch = fTree->GetBranch(fgkSimulatedDataBranchName);
 374      if (branch == 0x0)
 375       {
 376         Error("OpenFile","Can not find branch %s in file %s",
 377               fgkSimulatedDataBranchName.Data(),filename.Data());
 378
 379         delete fTree;
 380         fTree = 0x0;
 381         delete fFile;
 382         fFile = 0x0;
 383         return 6;
 384       }
 385      fTree->SetBranchAddress(fgkSimulatedDataBranchName,&fSimBuffer);
 386    }
 387 //  Info("ReadNext","Got TREE, Addresses are set.");
 388 //  Info("ReadNext","Quitting the method.");
 389
 390   return 0;
 391
 392 }
 393 /********************************************************************/
 394
 395 Int_t AliReaderAOD::WriteAOD(AliReader* reader, const char* outfilename, const char* pclassname,  Bool_t /*multcheck*/)
 396 {
 397 //reads tracks from runs and writes them to file
 398   ::Info("AliReaderAOD::Write","________________________________________________________");
 399   ::Info("AliReaderAOD::Write","________________________________________________________");
 400   ::Info("AliReaderAOD::Write","________________________________________________________");
 401
 402   if (reader == 0x0)
 403    {
 404      ::Error("AliReaderAOD::Write","Input Reader is NULL");
 405      return -1;
 406    }
 407   TFile *outfile = TFile::Open(outfilename,"recreate");
 408   if (outfile == 0x0)
 409    {
 410      ::Error("AliReaderAOD::Write","Can not open output file %s",outfilename);
 411      return -1;
 412    }
 413
 414   TTree *tree = new TTree(fgkTreeName,"Tree with tracks");
 415
 416   TBranch *recbranch = 0x0, *simbranch = 0x0;
 417
 418   AliAOD* eventrec = new AliAOD();//must be created before Branch is called. Otherwise clones array is not splitted
 419   AliAOD* eventsim = new AliAOD();//AOD together with fParticles clones array knowing exact type of particles
 420
 421   eventrec->SetParticleClassName(pclassname);
 422   eventsim->SetParticleClassName(pclassname);
 423
 424   AliAOD* recbuffer = eventrec;
 425   AliAOD* simbuffer = eventsim;
 426
 427   if (reader->ReadsRec()) recbranch = tree->Branch(fgkReconstructedDataBranchName,"AliAOD",&recbuffer,32000,99);
 428   if (reader->ReadsSim()) simbranch = tree->Branch(fgkSimulatedDataBranchName,"AliAOD",&simbuffer,32000,99);
 429
 430   reader->Rewind();
 431   while (reader->Next() == kFALSE)
 432    {
 433
 434      if (reader->ReadsRec())
 435       {//here we can get AOD that has different particle type
 436         AliAOD* event = reader->GetEventRec();
 437         if ( eventrec->GetParticleClass() != event->GetParticleClass() )
 438          {//if class type is not what what we whant we copy particles
 439            eventrec->CopyData(event);
 440            recbuffer = eventrec;
 441          }
 442         else
 443          {//else just pointer to event from input reader is passed
 444            recbuffer = event;
 445          }
 446       }
 447
 448      if (reader->ReadsSim())
 449       {
 450         AliAOD* event = reader->GetEventSim();
 451         if ( eventsim->GetParticleClass() != event->GetParticleClass() )
 452          {//if class type is not what what we whant we copy particles
 453            eventsim->CopyData(event);
 454            simbuffer = eventrec;
 455          }
 456         else
 457          {//else just pointer to event from input reader is passed
 458            simbuffer = event;
 459          }
 460       }
 461      tree->Fill();
 462    }
 463
 464   ::Info("AliReaderAOD::Write","Written %d events",tree->GetEntries());
 465   outfile->cd();
 466   tree->Write();
 467
 468   delete eventsim;
 469   delete eventrec;
 470
 471   delete tree;
 472   delete outfile;
 473   return 0;
 474 }
 475