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