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