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():
63 //____________________________________________________________________________
65 AliJetESDReader::~AliJetESDReader()
74 //____________________________________________________________________________
76 void AliJetESDReader::OpenInputFiles()
79 fChain = new TChain("esdTree");
81 // get directory and pattern name from the header
82 const char* dirName=fReaderHeader->GetDirectory();
83 const char* pattern=fReaderHeader->GetPattern();
85 // // Add files matching patters to the chain
87 void *dir = gSystem->OpenDirectory(dirName);
88 const char *name = 0x0;
89 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
91 while ((name = gSystem->GetDirEntry(dir))){
92 if (a>=nesd) continue;
94 if (strstr(name,pattern)){
96 sprintf(path,"%s/%s/AliESDs.root",dirName,name);
97 fChain->AddFile(path);
102 gSystem->FreeDirectory(dir);
106 fChain->SetBranchAddress("ESD", &fESD);
108 int nMax = fChain->GetEntries();
110 printf("\nTotal number of events in chain= %d",nMax);
112 // set number of events in header
113 if (fReaderHeader->GetLastEvent() == -1)
114 fReaderHeader->SetLastEvent(nMax);
116 Int_t nUsr = fReaderHeader->GetLastEvent();
117 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
121 void AliJetESDReader::ConnectTree(TTree* tree) {
123 fChain = (TChain*) tree;
125 fChain->SetBranchAddress("ESD", &fESD);
126 Int_t nMax = fChain->GetEntries();
127 printf("\nTotal number of events in chain= %5d", nMax);
128 // set number of events in header
129 if (fReaderHeader->GetLastEvent() == -1)
130 fReaderHeader->SetLastEvent(nMax);
132 Int_t nUsr = fReaderHeader->GetLastEvent();
133 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
137 //____________________________________________________________________________
139 Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
141 // Fill momentum array
148 // clear momentum array
150 fDebug = fReaderHeader->GetDebug();
152 // get event from chain
153 fChain->GetTree()->GetEntry(event);
159 // get number of tracks in event (for the loop)
160 nt = fESD->GetNumberOfTracks();
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];
165 // get cuts set by user
166 Float_t ptMin = fReaderHeader->GetPtCut();
167 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
168 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
171 for (Int_t it = 0; it < nt; it++) {
172 AliESDtrack *track = fESD->GetTrack(it);
173 UInt_t status = track->GetStatus();
176 track->GetPxPyPz(mom);
177 p3.SetXYZ(mom[0],mom[1],mom[2]);
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
186 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
188 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
190 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
192 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
195 // set the signal flags
196 fSignalFlag.Set(goodTrack,sflag);
197 fCutFlag.Set(goodTrack,cflag);
201 AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks();
202 fillUAFromTracks->SetReaderHeader(fReaderHeader);
203 fillUAFromTracks->SetMomentumArray(fMomentumArray);
204 fillUAFromTracks->SetTPCGrid(fTpcGrid);
205 fillUAFromTracks->SetEMCalGrid(fEmcalGrid);
206 fillUAFromTracks->SetHadCorrection(fHCorrection);
207 fillUAFromTracks->SetHadCorrector(fHadCorr);
208 fNeta = fillUAFromTracks->GetNeta();
209 fNphi = fillUAFromTracks->GetNphi();
210 fillUAFromTracks->SetActive(kFALSE);
211 // TPC only or Digits+TPC or Clusters+TPC
212 if(fOpt%2==!0 && fOpt!=0){
213 fillUAFromTracks->SetActive(kTRUE);
214 fillUAFromTracks->SetUnitArray(fUnitArray);
215 fillUAFromTracks->ExecuteTask("tpc");
217 delete fillUAFromTracks;
222 void AliJetESDReader::SetEMCALGeometry()
224 // Define EMCAL geometry to be able to read ESDs
225 fGeom = AliJetDummyGeo::GetInstance();
227 fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
229 // To be setted to run some AliEMCALGeometry functions
230 TGeoManager::Import("geometry.root");
231 fGeom->GetTransformationForSM();
233 printf("\n EMCal Geometry set ! \n");
237 void AliJetESDReader::InitParameters()
239 // Initialise parameters
240 fHCorrection = 0; // For hadron correction
241 fHadCorr = 0; // For hadron correction
242 fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
243 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
246 void AliJetESDReader::InitUnitArray()
248 //Initialises unit arrays
249 Int_t nElements = fTpcGrid->GetNEntries();
250 if(fArrayInitialised) delete [] fUnitArray;
251 if(fTpcGrid->GetGridType()==0)
252 fUnitArray = new AliJetUnitArray[nElements];
253 if(fTpcGrid->GetGridType()==1)
254 fUnitArray = new AliJetUnitArray[fNumUnits+nElements];
255 fArrayInitialised = 1;