Method SetInputEvent introduced.
[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 "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