1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------------------------
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 //-------------------------------------------------------------------------
24 #include <Riostream.h>
26 #include <TLorentzVector.h>
28 #include <TGeoManager.h>
30 #include "AliJetESDReader.h"
31 #include "AliJetESDReaderHeader.h"
33 #include "AliESDtrack.h"
34 //#include "AliEMCALGeometry.h"
35 #include "AliJetDummyGeo.h"
36 #include "AliJetFillUnitArrayTracks.h"
37 #include "AliJetUnitArray.h"
39 ClassImp(AliJetESDReader)
41 AliJetESDReader::AliJetESDReader():
62 //____________________________________________________________________________
64 AliJetESDReader::~AliJetESDReader()
73 //____________________________________________________________________________
75 void AliJetESDReader::OpenInputFiles()
78 fChain = new TChain("esdTree");
80 // get directory and pattern name from the header
81 const char* dirName=fReaderHeader->GetDirectory();
82 const char* pattern=fReaderHeader->GetPattern();
84 // // Add files matching patters to the chain
86 void *dir = gSystem->OpenDirectory(dirName);
87 const char *name = 0x0;
88 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
90 while ((name = gSystem->GetDirEntry(dir))){
91 if (a>=nesd) continue;
93 if (strstr(name,pattern)){
95 sprintf(path,"%s/%s/AliESDs.root",dirName,name);
96 fChain->AddFile(path);
101 gSystem->FreeDirectory(dir);
105 fChain->SetBranchAddress("ESD", &fESD);
107 int nMax = fChain->GetEntries();
109 printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
111 // set number of events in header
112 if (fReaderHeader->GetLastEvent() == -1)
113 fReaderHeader->SetLastEvent(nMax);
115 Int_t nUsr = fReaderHeader->GetLastEvent();
116 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
120 void AliJetESDReader::ConnectTree(TTree* tree, TObject* data) {
122 fChain = (TChain*) tree;
123 fESD = (AliESD*) data;
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);
131 Int_t nUsr = fReaderHeader->GetLastEvent();
132 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
136 //____________________________________________________________________________
138 Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
140 // Fill momentum array
147 // clear momentum array
149 fDebug = fReaderHeader->GetDebug();
150 // get event from chain
151 fChain->GetTree()->GetEntry(event);
157 // get number of tracks in event (for the loop)
158 nt = fESD->GetNumberOfTracks();
159 // temporary storage of signal and pt cut flag
160 Int_t* sflag = new Int_t[nt];
161 Int_t* cflag = new Int_t[nt];
163 // get cuts set by user
164 Float_t ptMin = fReaderHeader->GetPtCut();
165 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
166 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
169 for (Int_t it = 0; it < nt; it++) {
170 AliESDtrack *track = fESD->GetTrack(it);
171 UInt_t status = track->GetStatus();
174 track->GetPxPyPz(mom);
175 p3.SetXYZ(mom[0],mom[1],mom[2]);
177 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
178 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
179 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
180 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
181 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
182 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
184 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
186 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
188 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
190 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
193 // set the signal flags
194 fSignalFlag.Set(goodTrack,sflag);
195 fCutFlag.Set(goodTrack,cflag);
199 if (fTpcGrid || fEmcalGrid) {
202 AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks();
203 fillUAFromTracks->SetReaderHeader(fReaderHeader);
204 fillUAFromTracks->SetMomentumArray(fMomentumArray);
205 fillUAFromTracks->SetTPCGrid(fTpcGrid);
206 fillUAFromTracks->SetEMCalGrid(fEmcalGrid);
207 fillUAFromTracks->SetHadCorrection(fHCorrection);
208 fillUAFromTracks->SetHadCorrector(fHadCorr);
209 fNeta = fillUAFromTracks->GetNeta();
210 fNphi = fillUAFromTracks->GetNphi();
211 fillUAFromTracks->SetActive(kFALSE);
212 // TPC only or Digits+TPC or Clusters+TPC
213 if(fOpt%2==!0 && fOpt!=0) {
214 fillUAFromTracks->SetActive(kTRUE);
215 fillUAFromTracks->SetUnitArray(fUnitArray);
216 fillUAFromTracks->ExecuteTask("tpc");
219 delete fillUAFromTracks;
226 void AliJetESDReader::SetEMCALGeometry()
228 // Define EMCAL geometry to be able to read ESDs
229 fGeom = AliJetDummyGeo::GetInstance();
231 fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
233 // To be setted to run some AliEMCALGeometry functions
234 TGeoManager::Import("geometry.root");
235 fGeom->GetTransformationForSM();
236 printf("\n EMCal Geometry set ! \n");
240 void AliJetESDReader::InitParameters()
242 // Initialise parameters
243 fHCorrection = 0; // For hadron correction
244 fHadCorr = 0; // For hadron correction
245 fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
246 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
249 void AliJetESDReader::InitUnitArray()
251 //Initialises unit arrays
252 Int_t nElements = fTpcGrid->GetNEntries();
253 if(fArrayInitialised) delete [] fUnitArray;
254 if(fTpcGrid->GetGridType()==0)
255 fUnitArray = new AliJetUnitArray[nElements];
256 if(fTpcGrid->GetGridType()==1)
257 fUnitArray = new AliJetUnitArray[fNumUnits+nElements];
258 fArrayInitialised = 1;