JETAN/AliJetESDmcReader.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 // Jet ESD Reader
  17 // ESD reader for jet analysis (it reads the esd and the MC trees)
  18 // Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
  19
  20 #include <Riostream.h>
  21 #include <TSystem.h>
  22 #include <TLorentzVector.h>
  23 #include <TVector3.h>
  24 #include <TPDGCode.h>
  25 #include <TParticle.h>
  26 #include <TParticlePDG.h>
  27 #include <TMath.h>
  28 #include "AliJetESDmcReader.h"
  29 #include "AliJetESDReaderHeader.h"
  30 #include "AliESD.h"
  31 #include "AliESDtrack.h"
  32 #include "AliHeader.h"
  33
  34 ClassImp(AliJetESDmcReader)
  35
  36 AliJetESDmcReader::AliJetESDmcReader():
  37     AliJetESDReader(),
  38     fChainMC(0x0),
  39     fAliHeader(0)
  40 {
  41   // Constructor
  42 }
  43
  44 //____________________________________________________________________________
  45
  46 AliJetESDmcReader::~AliJetESDmcReader()
  47 {
  48   // Destructor
  49     delete fChainMC;
  50 }
  51
  52 //____________________________________________________________________________
  53
  54 void AliJetESDmcReader::OpenInputFiles()
  55
  56 {
  57   // chain for the ESDs
  58   fChain   = new TChain("esdTree");
  59   fChainMC = new TChain("mcStackTree");
  60
  61   // get directory and pattern name from the header
  62   const char* dirName=fReaderHeader->GetDirectory();
  63   const char* pattern=fReaderHeader->GetPattern();
  64
  65   // Add files matching patters to the chain
  66   void *dir  = gSystem->OpenDirectory(dirName);
  67   const char *name = 0x0;
  68   int nesd = ((AliJetESDReaderHeader*)fReaderHeader)->GetNesd();
  69   int a = 0;
  70   while ((name = gSystem->GetDirEntry(dir))){
  71     if (a>=nesd) continue;
  72     if (strstr(name,pattern)){
  73       printf("Adding %s\n",name);
  74       char path[256];
  75       sprintf(path,"%s/%s",dirName,name);
  76       fChain->AddFile(path,-1);
  77       fChainMC->AddFile(path,-1);
  78       a++;
  79     }
  80   }
  81   printf("%d ESDs added\n",a);
  82
  83   gSystem->FreeDirectory(dir);
  84   fChain->SetBranchAddress("ESD",    &fESD);
  85   fChainMC->SetBranchAddress("Header", &fAliHeader);
  86   fChainMC->SetBranchAddress("Stack",  &fArrayMC);
  87
  88   int nMax = fChain->GetEntries();
  89   printf("\nTotal number of events in chain= %d\n",nMax);
  90
  91   // set number of events in header
  92   if (fReaderHeader->GetLastEvent() == -1)
  93     fReaderHeader->SetLastEvent(nMax);
  94   else {
  95     Int_t nUsr = fReaderHeader->GetLastEvent();
  96     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
  97   }
  98 }
  99
 100 //____________________________________________________________________________
 101
 102 Bool_t AliJetESDmcReader::FillMomentumArray(Int_t event)
 103 {
 104   // Fill momentum array
 105
 106   Int_t goodTrack = 0;
 107   Int_t nt = 0;
 108   Int_t pdgCode = 0;
 109   Float_t pt, eta;
 110   Float_t energyMC, pxMC, pyMC, pzMC, ptMC; // Monte Carlo
 111   TVector3 p3;
 112
 113   // clear momentum array
 114   ClearArray();
 115
 116   // get event from chain
 117   fChain->GetEntry(event);
 118   fChainMC->GetEntry(event);
 119
 120   // get number of tracks in event (for the loop)
 121   nt = fESD->GetNumberOfTracks();
 122
 123  // temporary storage of signal and pt cut flag
 124   Int_t* sflag  = new Int_t[nt];
 125   Int_t* cflag  = new Int_t[nt];
 126
 127   // get cuts set by user
 128   Float_t ptMin = fReaderHeader->GetPtCut();
 129   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 130   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 131
 132   //loop over tracks
 133   for (Int_t it = 0; it < nt; it++) {
 134       AliESDtrack *track = fESD->GetTrack(it);
 135       UInt_t status = track->GetStatus();
 136
 137       Double_t mom[3];
 138       track->GetPxPyPz(mom);
 139       p3.SetXYZ(mom[0],mom[1],mom[2]);
 140       pt = p3.Pt();
 141       if (((status & AliESDtrack::kITSrefit) == 0) ||
 142           ((status & AliESDtrack::kTPCrefit) == 0)) continue;    // quality check
 143       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 144           && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
 145       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 146           && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
 147       eta = p3.Eta();
 148       if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
 149
 150       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
 151
 152       sflag[goodTrack]=0;
 153       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 154       cflag[goodTrack]=0;
 155       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 156
 157       goodTrack++;
 158
 159       // Monte Carlo information
 160       Int_t label = TMath::Abs(track->GetLabel());
 161       TClonesArray &arrayMC = *fArrayMC;
 162       TParticle *part = (TParticle*)arrayMC[label]; //particle
 163       ptMC = part->Pt();
 164       pdgCode = part->GetPdgCode();
 165       energyMC = part->Energy();
 166       pxMC = part->Px(); pyMC = part->Py(); pzMC = part->Pz();
 167   }
 168   // set the signal flags
 169   fSignalFlag.Set(goodTrack,sflag);
 170   fCutFlag.Set(goodTrack,cflag);
 171   return kTRUE;
 172
 173 }
 174
 175