Correction in output format.
[u/mrichter/AliRoot.git] / JETAN / AliJetESDReader.cxx
CommitLineData
99e5fe42 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 **************************************************************************/
1994c1b0 15
16//-------------------------------------------------------------------------
99e5fe42 17// Jet ESD Reader
18// ESD reader for jet analysis
b45b0c92 19// Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
20// Magali Estienne <magali.estienne@IReS.in2p3.fr>
1994c1b0 21//-------------------------------------------------------------------------
22
99e5fe42 23
24#include <Riostream.h>
25#include <TSystem.h>
26#include <TLorentzVector.h>
83a444b1 27#include <TVector3.h>
b45b0c92 28#include <TGeoManager.h>
29
99e5fe42 30#include "AliJetESDReader.h"
31#include "AliJetESDReaderHeader.h"
32#include "AliESD.h"
33#include "AliESDtrack.h"
5b81d1ce 34//#include "AliEMCALGeometry.h"
35#include "AliJetDummyGeo.h"
b45b0c92 36#include "AliJetFillUnitArrayTracks.h"
37#include "AliJetUnitArray.h"
99e5fe42 38
39ClassImp(AliJetESDReader)
40
1b7d5d7e 41AliJetESDReader::AliJetESDReader():
b45b0c92 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)
99e5fe42 58{
59 // Constructor
99e5fe42 60}
61
62//____________________________________________________________________________
63
64AliJetESDReader::~AliJetESDReader()
65{
66 // Destructor
b45b0c92 67 delete fChain;
68 delete fESD;
69 delete fTpcGrid;
70 delete fEmcalGrid;
99e5fe42 71}
72
73//____________________________________________________________________________
74
75void AliJetESDReader::OpenInputFiles()
99e5fe42 76{
77 // chain for the ESDs
78 fChain = new TChain("esdTree");
99e5fe42 79
80 // get directory and pattern name from the header
7d0f353c 81 const char* dirName=fReaderHeader->GetDirectory();
82 const char* pattern=fReaderHeader->GetPattern();
99e5fe42 83
7d0f353c 84// // Add files matching patters to the chain
85
86 void *dir = gSystem->OpenDirectory(dirName);
87 const char *name = 0x0;
b45b0c92 88 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
7d0f353c 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 }
1994c1b0 100
101 gSystem->FreeDirectory(dir);
7d0f353c 102
103
104 fESD = 0;
1994c1b0 105 fChain->SetBranchAddress("ESD", &fESD);
106
99e5fe42 107 int nMax = fChain->GetEntries();
7d0f353c 108
5a50fd9a 109 printf("\nTotal number of events in chain= %d \n",nMax);
1994c1b0 110
99e5fe42 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
7d0f353c 120void AliJetESDReader::ConnectTree(TTree* tree) {
b45b0c92 121 // Connect the tree
7d0f353c 122 fChain = (TChain*) tree;
123
124 fChain->SetBranchAddress("ESD", &fESD);
125 Int_t nMax = fChain->GetEntries();
5a50fd9a 126 printf("\nTotal number of events in chain= %5d \n", nMax);
7d0f353c 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
99e5fe42 136//____________________________________________________________________________
137
7d0f353c 138Bool_t AliJetESDReader::FillMomentumArray(Int_t event)
99e5fe42 139{
1994c1b0 140 // Fill momentum array
141
99e5fe42 142 Int_t goodTrack = 0;
143 Int_t nt = 0;
83a444b1 144 Float_t pt, eta;
145 TVector3 p3;
1994c1b0 146
99e5fe42 147 // clear momentum array
69b18641 148 ClearArray();
149 fDebug = fReaderHeader->GetDebug();
99e5fe42 150 // get event from chain
7d0f353c 151 fChain->GetTree()->GetEntry(event);
5a50fd9a 152
7d0f353c 153 if (!fESD) {
154 return kFALSE;
155 }
1994c1b0 156
99e5fe42 157 // get number of tracks in event (for the loop)
158 nt = fESD->GetNumberOfTracks();
1994c1b0 159 // temporary storage of signal and pt cut flag
83a444b1 160 Int_t* sflag = new Int_t[nt];
161 Int_t* cflag = new Int_t[nt];
1994c1b0 162
99e5fe42 163 // get cuts set by user
7d0f353c 164 Float_t ptMin = fReaderHeader->GetPtCut();
83a444b1 165 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
166 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
1994c1b0 167
99e5fe42 168 //loop over tracks
169 for (Int_t it = 0; it < nt; it++) {
7d0f353c 170 AliESDtrack *track = fESD->GetTrack(it);
171 UInt_t status = track->GetStatus();
172
173 Double_t mom[3];
174 track->GetPxPyPz(mom);
175 p3.SetXYZ(mom[0],mom[1],mom[2]);
176 pt = p3.Pt();
b45b0c92 177 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
178 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
7d0f353c 179 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
b45b0c92 180 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
7d0f353c 181 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
b45b0c92 182 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
7d0f353c 183 eta = p3.Eta();
b45b0c92 184 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
7d0f353c 185
186 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
187 sflag[goodTrack]=0;
188 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
189 cflag[goodTrack]=0;
69b18641 190 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
7d0f353c 191 goodTrack++;
99e5fe42 192 }
8011d399 193 // set the signal flags
83a444b1 194 fSignalFlag.Set(goodTrack,sflag);
195 fCutFlag.Set(goodTrack,cflag);
b45b0c92 196
197//
198//
69b18641 199 if (fTpcGrid || fEmcalGrid) {
200 SetEMCALGeometry();
201 InitParameters();
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");
217 }
218
219 delete fillUAFromTracks;
b45b0c92 220 }
69b18641 221
7d0f353c 222 return kTRUE;
99e5fe42 223}
224
b45b0c92 225
226void AliJetESDReader::SetEMCALGeometry()
227{
228 // Define EMCAL geometry to be able to read ESDs
5b81d1ce 229 fGeom = AliJetDummyGeo::GetInstance();
b45b0c92 230 if (fGeom == 0)
5b81d1ce 231 fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
b45b0c92 232
233 // To be setted to run some AliEMCALGeometry functions
234 TGeoManager::Import("geometry.root");
235 fGeom->GetTransformationForSM();
b45b0c92 236 printf("\n EMCal Geometry set ! \n");
237
238}
7d0f353c 239
b45b0c92 240void AliJetESDReader::InitParameters()
241{
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");
247}
248
249void AliJetESDReader::InitUnitArray()
250{
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;
259}
7d0f353c 260
99e5fe42 261