364c8b308e65de6bdbada97aab129f76dfa02621
[u/mrichter/AliRoot.git] / 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 "AliESD.h"
33 #include "AliESDtrack.h"
34 #include "AliEMCALGeometry.h"
35 #include "AliJetFillUnitArrayTracks.h"
36 #include "AliJetUnitArray.h"
37
38 ClassImp(AliJetESDReader)
39
40 AliJetESDReader::AliJetESDReader():
41     AliJetReader(),  
42     fGeom(0),
43     fChain(0x0),
44     fESD(0x0),
45     fHadCorr(0x0),
46     fTpcGrid(0x0),
47     fEmcalGrid(0x0),
48     fPtCut(0),
49     fHCorrection(0),
50     fNumUnits(0),
51     fDebug(0),
52     fNIn(0),
53     fOpt(0),
54     fNeta(0),
55     fNphi(0),
56     fArrayInitialised(0) 
57 {
58   // Constructor    
59     SetEMCALGeometry();
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("\nTotal number of events in chain= %d",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) {
121     // Connect the tree
122      fChain = (TChain*) tree;
123      
124      fChain->SetBranchAddress("ESD",    &fESD);
125      Int_t nMax = fChain->GetEntries(); 
126      printf("\nTotal number of events in chain= %5d", 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    InitParameters();
151   // get event from chain
152   fChain->GetTree()->GetEntry(event);
153
154   if (!fESD) {
155       return kFALSE;
156   }
157   
158   // get number of tracks in event (for the loop)
159   nt = fESD->GetNumberOfTracks();
160   // temporary storage of signal and pt cut flag
161   Int_t* sflag  = new Int_t[nt];
162   Int_t* cflag  = new Int_t[nt];
163   
164   // get cuts set by user
165   Float_t ptMin =  fReaderHeader->GetPtCut();
166   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
167   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
168   
169   //loop over tracks
170   for (Int_t it = 0; it < nt; it++) {
171       AliESDtrack *track = fESD->GetTrack(it);
172       UInt_t status = track->GetStatus();
173       
174       Double_t mom[3];
175       track->GetPxPyPz(mom);
176       p3.SetXYZ(mom[0],mom[1],mom[2]);
177       pt = p3.Pt();
178       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
179       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
180       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
181           && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
182       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
183           && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
184       eta = p3.Eta();
185       if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
186       
187       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
188       sflag[goodTrack]=0;
189       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
190       cflag[goodTrack]=0;
191       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
192       goodTrack++;
193   }
194   // set the signal flags
195   fSignalFlag.Set(goodTrack,sflag);
196   fCutFlag.Set(goodTrack,cflag);
197
198 //
199 //
200   AliJetFillUnitArrayTracks *fillUAFromTracks = new AliJetFillUnitArrayTracks(); 
201   fillUAFromTracks->SetReaderHeader(fReaderHeader);
202   fillUAFromTracks->SetMomentumArray(fMomentumArray);
203   fillUAFromTracks->SetTPCGrid(fTpcGrid);
204   fillUAFromTracks->SetEMCalGrid(fEmcalGrid);
205   fillUAFromTracks->SetHadCorrection(fHCorrection);
206   fillUAFromTracks->SetHadCorrector(fHadCorr);
207   fNeta = fillUAFromTracks->GetNeta();
208   fNphi = fillUAFromTracks->GetNphi();
209   fillUAFromTracks->SetActive(kFALSE);
210   // TPC only or Digits+TPC or Clusters+TPC
211   if(fOpt%2==!0 && fOpt!=0){ 
212     fillUAFromTracks->SetActive(kTRUE);
213     fillUAFromTracks->SetUnitArray(fUnitArray);
214     fillUAFromTracks->ExecuteTask("tpc");
215   }
216   delete fillUAFromTracks;
217   return kTRUE;
218 }
219
220
221 void AliJetESDReader::SetEMCALGeometry()
222 {
223   // Define EMCAL geometry to be able to read ESDs
224     fGeom = AliEMCALGeometry::GetInstance();
225     if (fGeom == 0)
226         fGeom = AliEMCALGeometry::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
227     
228     // To be setted to run some AliEMCALGeometry functions
229     TGeoManager::Import("geometry.root");
230     fGeom->GetTransformationForSM();  
231     
232     printf("\n EMCal Geometry set ! \n");
233
234 }
235   
236 void AliJetESDReader::InitParameters()
237 {
238     // Initialise parameters
239     fHCorrection    = 0;                 // For hadron correction
240     fHadCorr        = 0;                 // For hadron correction
241     fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
242     if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
243 }
244
245 void AliJetESDReader::InitUnitArray()
246 {
247   //Initialises unit arrays
248     Int_t nElements = fTpcGrid->GetNEntries();
249     if(fArrayInitialised) delete [] fUnitArray;
250     if(fTpcGrid->GetGridType()==0)
251         fUnitArray = new AliJetUnitArray[nElements];
252     if(fTpcGrid->GetGridType()==1)
253         fUnitArray = new AliJetUnitArray[fNumUnits+nElements];
254     fArrayInitialised = 1;
255 }
256
257