JETAN/AliJetESDReader.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 //-------------------------------------------------------------------------
  17 // Jet ESD Reader
  18 // ESD reader for jet analysis
  19 // Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
  20 //          Magali Estienne <magali.estienne@IReS.in2p3.fr>
  21 //-------------------------------------------------------------------------
  22
  23
  24 #include <Riostream.h>
  25 #include <TSystem.h>
  26 #include <TLorentzVector.h>
  27 #include <TVector3.h>
  28 #include <TGeoManager.h>
  29
  30 #include "AliJetESDReader.h"
  31 #include "AliJetESDReaderHeader.h"
  32 #include "AliESDEvent.h"
  33 #include "AliESDtrack.h"
  34 //#include "AliEMCALGeometry.h"
  35 #include "AliJetDummyGeo.h"
  36 #include "AliJetFillUnitArrayTracks.h"
  37 #include "AliJetUnitArray.h"
  38
  39 ClassImp(AliJetESDReader)
  40
  41 AliJetESDReader::AliJetESDReader():
  42     AliJetReader(),
  43     fGeom(0),
  44     fChain(0x0),
  45     fESD(0x0),
  46     fHadCorr(0x0),
  47     fTpcGrid(0x0),
  48     fEmcalGrid(0x0),
  49     fPtCut(0),
  50     fHCorrection(0),
  51     fNumUnits(0),
  52     fDebug(0),
  53     fNIn(0),
  54     fOpt(0),
  55     fNeta(0),
  56     fNphi(0),
  57     fArrayInitialised(0)
  58 {
  59   // Constructor
  60 }
  61
  62 //____________________________________________________________________________
  63
  64 AliJetESDReader::~AliJetESDReader()
  65 {
  66   // Destructor
  67     delete fChain;
  68     delete fESD;
  69     delete fTpcGrid;
  70     delete fEmcalGrid;
  71 }
  72
  73 //____________________________________________________________________________
  74
  75 void AliJetESDReader::OpenInputFiles()
  76 {
  77   // chain for the ESDs
  78   fChain   = new TChain("esdTree");
  79
  80   // get directory and pattern name from the header
  81    const char* dirName=fReaderHeader->GetDirectory();
  82    const char* pattern=fReaderHeader->GetPattern();
  83
  84 //   // Add files matching patters to the chain
  85
  86    void *dir  = gSystem->OpenDirectory(dirName);
  87    const char *name = 0x0;
  88    int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
  89    int a = 0;
  90    while ((name = gSystem->GetDirEntry(dir))){
  91        if (a>=nesd) continue;
  92
  93        if (strstr(name,pattern)){
  94            char path[256];
  95            sprintf(path,"%s/%s/AliESDs.root",dirName,name);
  96            fChain->AddFile(path);
  97            a++;
  98        }
  99    }
 100
 101   gSystem->FreeDirectory(dir);
 102
 103
 104   fESD = 0;
 105   fChain->SetBranchAddress("ESD",    &fESD);
 106
 107   int nMax = fChain->GetEntries();
 108
 109   printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
 110
 111   // set number of events in header
 112   if (fReaderHeader->GetLastEvent() == -1)
 113     fReaderHeader->SetLastEvent(nMax);
 114   else {
 115     Int_t nUsr = fReaderHeader->GetLastEvent();
 116     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 117   }
 118 }
 119
 120 void AliJetESDReader::ConnectTree(TTree* tree, TObject* data) {
 121     // Connect the tree
 122      fChain = (TChain*)      tree;
 123      fESD   = (AliESDEvent*) data;
 124
 125      Int_t nMax = fChain->GetEntries();
 126      printf("\n AliJetESDReader: Total number of events in chain= %5d \n", nMax);
 127      // set number of events in header
 128      if (fReaderHeader->GetLastEvent() == -1)
 129          fReaderHeader->SetLastEvent(nMax);
 130      else {
 131          Int_t nUsr = fReaderHeader->GetLastEvent();
 132          fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 133      }
 134 }
 135
 136 //____________________________________________________________________________
 137
 138 Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
 139 {
 140   // Fill momentum array
 141
 142   Int_t goodTrack = 0;
 143   Int_t nt = 0;
 144   Float_t pt, eta;
 145   TVector3 p3;
 146
 147   // clear momentum array
 148   ClearArray();
 149   fDebug = fReaderHeader->GetDebug();
 150   // get event from chain
 151   // fChain->GetTree()->GetEntry(event);
 152
 153   if (!fESD) {
 154       return kFALSE;
 155   }
 156
 157   // get number of tracks in event (for the loop)
 158   nt = fESD->GetNumberOfTracks();
 159   printf("Fill Momentum Array %5d ", nt);
 160
 161   // temporary storage of signal and pt cut flag
 162   Int_t* sflag  = new Int_t[nt];
 163   Int_t* cflag  = new Int_t[nt];
 164
 165   // get cuts set by user
 166   Float_t ptMin =  fReaderHeader->GetPtCut();
 167   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 168   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 169
 170   //loop over tracks
 171   for (Int_t it = 0; it < nt; it++) {
 172       AliESDtrack *track = fESD->GetTrack(it);
 173       UInt_t status = track->GetStatus();
 174
 175       Double_t mom[3];
 176       track->GetPxPyPz(mom);
 177       p3.SetXYZ(mom[0],mom[1],mom[2]);
 178       pt = p3.Pt();
 179       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
 180       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
 181       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 182           && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
 183       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 184           && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
 185       eta = p3.Eta();
 186       if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
 187
 188       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
 189       sflag[goodTrack]=0;
 190       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 191       cflag[goodTrack]=0;
 192       if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
 193       goodTrack++;
 194   }
 195   // set the signal flags
 196   fSignalFlag.Set(goodTrack,sflag);
 197   fCutFlag.Set(goodTrack,cflag);
 198
 199 //
 200 //
 201   if (fTpcGrid || fEmcalGrid) {
 202       SetEMCALGeometry();
 203       InitParameters();
 204       AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks();
 205       fillUAFromTracks->SetReaderHeader(fReaderHeader);
 206       fillUAFromTracks->SetMomentumArray(fMomentumArray);
 207       fillUAFromTracks->SetTPCGrid(fTpcGrid);
 208       fillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 209       fillUAFromTracks->SetHadCorrection(fHCorrection);
 210       fillUAFromTracks->SetHadCorrector(fHadCorr);
 211       fNeta = fillUAFromTracks->GetNeta();
 212       fNphi = fillUAFromTracks->GetNphi();
 213       fillUAFromTracks->SetActive(kFALSE);
 214       // TPC only or Digits+TPC or Clusters+TPC
 215       if(fOpt%2==!0 && fOpt!=0) {
 216           fillUAFromTracks->SetActive(kTRUE);
 217           fillUAFromTracks->SetUnitArray(fUnitArray);
 218           fillUAFromTracks->ExecuteTask("tpc");
 219       }
 220
 221       delete fillUAFromTracks;
 222   }
 223
 224   return kTRUE;
 225 }
 226
 227
 228 void AliJetESDReader::SetEMCALGeometry()
 229 {
 230   // Define EMCAL geometry to be able to read ESDs
 231     fGeom = AliJetDummyGeo::GetInstance();
 232     if (fGeom == 0)
 233         fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
 234
 235     // To be setted to run some AliEMCALGeometry functions
 236     TGeoManager::Import("geometry.root");
 237     fGeom->GetTransformationForSM();
 238     printf("\n EMCal Geometry set ! \n");
 239
 240 }
 241
 242 void AliJetESDReader::InitParameters()
 243 {
 244     // Initialise parameters
 245     fHCorrection    = 0;                 // For hadron correction
 246     fHadCorr        = 0;                 // For hadron correction
 247     fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 248     if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 249 }
 250
 251 void AliJetESDReader::InitUnitArray()
 252 {
 253   //Initialises unit arrays
 254     Int_t nElements = fTpcGrid->GetNEntries();
 255     if(fArrayInitialised) delete [] fUnitArray;
 256     if(fTpcGrid->GetGridType()==0)
 257         fUnitArray = new AliJetUnitArray[nElements];
 258     if(fTpcGrid->GetGridType()==1)
 259         fUnitArray = new AliJetUnitArray[fNumUnits+nElements];
 260     fArrayInitialised = 1;
 261 }
 262
 263