ANALYSIS/AliReaderAOD.cxx

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