733b26ac2073bc554158312d736a0704be24d1b0
[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 // --- Standard library ---
24 #include <Riostream.h>
25
26 // --- ROOT system ---
27 #include <TSystem.h>
28 #include <TStopwatch.h>
29 #include <TLorentzVector.h>
30 #include <TVector3.h>
31 #include <TTask.h>
32 #include <TGeoManager.h>
33 #include <assert.h>
34 #include <TRefArray.h>
35 #include <TMath.h>
36 #include <TChain.h>
37
38
39 // --- AliRoot header files ---
40 #include "AliJetESDReader.h"
41 #include "AliJetESDReaderHeader.h"
42 #include "AliESDEvent.h"
43 #include "AliESD.h"
44 #include "AliESDtrack.h"
45 #include "AliJetDummyGeo.h"
46 #include "AliJetFillUnitArrayTracks.h"
47 #include "AliJetFillUnitArrayEMCalDigits.h"
48 #include "AliJetUnitArray.h"
49
50 ClassImp(AliJetESDReader)
51
52 AliJetESDReader::AliJetESDReader():
53   AliJetReader(),  
54   fGeom(0),
55   fChain(0x0),
56   fESD(0x0),
57   fHadCorr(0x0),
58   fTpcGrid(0x0),
59   fEmcalGrid(0x0),
60   fGrid0(0),
61   fGrid1(0),
62   fGrid2(0),
63   fGrid3(0),
64   fGrid4(0),
65   fPtCut(0),
66   fHCorrection(0),
67   fNumUnits(0),
68   fDebug(0),
69   fMass(0),
70   fSign(0),
71   fNIn(0),
72   fOpt(0),
73   fDZ(0),
74   fNeta(0),
75   fNphi(0),
76   fArrayInitialised(0)
77 {
78   // Constructor 
79 }
80
81 //____________________________________________________________________________
82 AliJetESDReader::~AliJetESDReader()
83 {
84   // Destructor
85     delete fChain;
86     delete fESD;
87     delete fTpcGrid;
88     delete fEmcalGrid;
89     if(fDZ)
90       {
91         delete fGrid0;
92         delete fGrid1;
93         delete fGrid2;
94         delete fGrid3;
95         delete fGrid4;
96       }
97 }
98
99 //____________________________________________________________________________
100 void AliJetESDReader::OpenInputFiles()
101 {
102   // chain for the ESDs
103   fChain   = new TChain("esdTree");
104
105   // get directory and pattern name from the header
106    const char* dirName=fReaderHeader->GetDirectory();
107    const char* pattern=fReaderHeader->GetPattern();
108     
109    // Add files matching patters to the chain
110    void *dir  = gSystem->OpenDirectory(dirName);
111    const char *name = 0x0;
112    int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
113    int a = 0;
114    while ((name = gSystem->GetDirEntry(dir))){
115        if (a>=nesd) continue;
116        
117        if (strstr(name,pattern)){
118          printf("Adding %s\n",name);
119          char path[256];
120          //        sprintf(path,"%s/%s/AliESDs.root",dirName,name);
121          sprintf(path,"%s/%s/AliESDs.root",dirName,name);
122          fChain->AddFile(path);
123          a++;
124        }
125    }
126   
127   gSystem->FreeDirectory(dir);
128   
129   fESD = new AliESDEvent();
130   fESD->ReadFromTree(fChain);
131
132   int nMax = fChain->GetEntries(); 
133
134   printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
135   
136   // set number of events in header
137   if (fReaderHeader->GetLastEvent() == -1)
138     fReaderHeader->SetLastEvent(nMax);
139   else {
140     Int_t nUsr = fReaderHeader->GetLastEvent();
141     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
142   }
143 }
144
145 //____________________________________________________________________________
146 void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
147     // Connect the tree
148      fESD   = (AliESDEvent*) esd;
149 }
150
151 //____________________________________________________________________________
152 Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
153 {
154   // Fill momentum array
155
156   Int_t goodTrack = 0;
157   Int_t nt = 0;
158   Float_t pt, eta;
159   TVector3 p3;
160   
161   // clear momentum array
162   ClearArray();
163   fDebug = fReaderHeader->GetDebug();
164   fOpt = fReaderHeader->GetDetector();
165
166   if (!fESD) {
167       return kFALSE;
168   }
169
170   // get number of tracks in event (for the loop)
171   nt = fESD->GetNumberOfTracks();
172   printf("Fill Momentum Array %5d  \n", nt);
173   
174   // temporary storage of signal and pt cut flag
175   Int_t* sflag  = new Int_t[nt];
176   Int_t* cflag  = new Int_t[nt];
177   
178   // get cuts set by user
179   Float_t ptMin =  fReaderHeader->GetPtCut();
180   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
181   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
182   
183   //loop over tracks
184   for (Int_t it = 0; it < nt; it++) {
185       AliESDtrack *track = fESD->GetTrack(it);
186       UInt_t status = track->GetStatus();
187       
188       Double_t mom[3];
189       track->GetPxPyPz(mom);
190       p3.SetXYZ(mom[0],mom[1],mom[2]);
191       pt = p3.Pt();
192       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
193       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
194       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly() 
195           && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
196       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly() 
197           && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
198       eta = p3.Eta();
199       if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
200       
201       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
202       sflag[goodTrack]=0;
203       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
204       cflag[goodTrack]=0;
205       if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
206       goodTrack++;
207   }
208   // set the signal flags
209   fSignalFlag.Set(goodTrack,sflag);
210   fCutFlag.Set(goodTrack,cflag);
211   return kTRUE;
212
213 }
214
215 //____________________________________________________________________________
216 void AliJetESDReader::CreateTasks()
217 {
218   fDebug = fReaderHeader->GetDebug();
219   fDZ = fReaderHeader->GetDZ();
220
221   // Init EMCAL geometry and create UnitArray object
222   SetEMCALGeometry();
223   InitParameters();
224   InitUnitArray();
225
226   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
227   fFillUAFromTracks = new AliJetFillUnitArrayTracks(); 
228   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
229   fFillUAFromTracks->SetGeom(fGeom);
230   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
231   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
232
233   if(fDZ)
234     {
235       fFillUAFromTracks->SetGrid0(fGrid0);
236       fFillUAFromTracks->SetGrid1(fGrid1);
237       fFillUAFromTracks->SetGrid2(fGrid2);
238       fFillUAFromTracks->SetGrid3(fGrid3);
239       fFillUAFromTracks->SetGrid4(fGrid4);
240     }
241   fFillUAFromTracks->SetHadCorrection(fHCorrection);
242   fFillUAFromTracks->SetHadCorrector(fHadCorr);
243   fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
244   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
245   fFillUAFromEMCalDigits->SetGeom(fGeom);
246   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
247   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
248   //      fFillUnitArray->Add(fFillUAFromTracks);
249   fFillUnitArray->Add(fFillUAFromEMCalDigits);
250   fFillUAFromTracks->SetActive(kFALSE);
251   fFillUAFromEMCalDigits->SetActive(kFALSE);
252
253   cout << "Tasks instantiated at that stage ! " << endl;
254   cout << "You can loop over events now ! " << endl;
255    
256 }
257
258 //____________________________________________________________________________
259 //void AliJetESDReader::ExecTasks(Int_t event)
260 Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/)
261 {
262   // clear momentum array
263   Int_t nEntRef = fRefArray->GetEntries();
264
265   for(Int_t i=0; i<nEntRef; i++)
266     {      
267       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitTrackID(0);
268       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.);
269       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller);
270       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc);
271       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet);
272     }
273
274   ClearArray();
275
276   fDebug = fReaderHeader->GetDebug();
277   fOpt = fReaderHeader->GetDetector();
278   //  InitParameters();
279
280   if(!fESD) {
281     return kFALSE;
282   }
283   
284   /*
285   // get event from chain
286   // For TSelectors
287   //  fChain->GetTree()->GetEntry(event);
288   // For interactive process
289   //  fChain->GetEntry(event);
290   fChain->GetEvent(event);
291   */
292
293   // TPC only or Digits+TPC or Clusters+TPC
294   if(fOpt%2==!0 && fOpt!=0){ 
295     fFillUAFromTracks->SetESD(fESD);
296     fFillUAFromTracks->SetActive(kTRUE);
297     fFillUAFromTracks->SetUnitArray(fUnitArray);
298     fFillUAFromTracks->SetRefArray(fRefArray);
299     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed  !!!
300                                                  // Incompatibility with Andrei's analysis framework
301     fFillUAFromTracks->Exec("tpc");
302     if(fOpt==1){
303       fNumCandidate = fFillUAFromTracks->GetMult();
304       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
305     }
306   }
307
308   // Digits only or Digits+TPC  
309   if(fOpt>=2 && fOpt<=3){
310     fFillUAFromEMCalDigits->SetESD(fESD);
311     fFillUAFromEMCalDigits->SetActive(kTRUE);
312     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
313     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
314     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
315     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
316     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!!
317     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
318     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
319   }
320
321   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
322
323   return kTRUE;
324 }
325
326 //____________________________________________________________________________
327 void AliJetESDReader::SetEMCALGeometry()
328 {
329   // Define EMCAL geometry to be able to read ESDs
330     fGeom = AliJetDummyGeo::GetInstance();
331     if (fGeom == 0)
332         fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
333
334     // To be setted to run some AliEMCALGeometry functions
335     TGeoManager::Import("geometry.root");
336     //    fGeom->GetTransformationForSM();  
337     printf("\n EMCal Geometry set ! \n");
338
339 }
340
341
342 //____________________________________________________________________________  
343 void AliJetESDReader::InitParameters()
344 {
345     // Initialise parameters
346     fHCorrection    = 0;                 // For hadron correction
347     fHadCorr        = 0;                 // For hadron correction
348     fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
349     if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
350 }
351
352 //____________________________________________________________________________
353 void AliJetESDReader::InitUnitArray()
354 {
355   //Initialises unit arrays
356   Int_t nElements = fTpcGrid->GetNEntries();
357   Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
358   if(fArrayInitialised) fUnitArray->Delete();
359
360   if(fTpcGrid->GetGridType()==0)
361     { // Fill the following quantities :
362       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi, 
363       // detector flag, in/out jet, pt cut, mass, cluster ID)
364       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
365         {
366           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
367           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
368           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
369           Deta = fTpcGrid->GetDeta();
370           Dphi = fTpcGrid->GetDphi();
371           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
372         }
373     }
374
375   if(fTpcGrid->GetGridType()==1)
376     {
377       Int_t nGaps = 0;
378       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
379
380       if(fDZ)
381         {
382           // Define a grid of cell for the gaps between SM
383           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
384           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
385           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
386           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
387           fGrid0->SetGridType(0);
388           fGrid0->SetMatrixIndexes();
389           fGrid0->SetIndexIJ();
390           n0 = fGrid0->GetNEntries();
391           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
392           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
393           fGrid1->SetGridType(0);
394           fGrid1->SetMatrixIndexes();
395           fGrid1->SetIndexIJ();
396           n1 = fGrid1->GetNEntries();
397           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
398           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
399           fGrid2->SetGridType(0);
400           fGrid2->SetMatrixIndexes();
401           fGrid2->SetIndexIJ();
402           n2 = fGrid2->GetNEntries();
403           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
404           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
405           fGrid3->SetGridType(0);  
406           fGrid3->SetMatrixIndexes();
407           fGrid3->SetIndexIJ();
408           n3 = fGrid3->GetNEntries();
409           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
410           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
411           fGrid4->SetGridType(0);
412           fGrid4->SetMatrixIndexes();
413           fGrid4->SetIndexIJ();
414           n4 = fGrid4->GetNEntries();
415
416           if(fDebug>1) 
417             {
418               cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
419               cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl;
420               cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl;
421               cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl;
422               cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl;
423             }
424           
425           nGaps = n0+n1+n2+n3+n4;
426
427         }
428
429       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
430         {
431           if(nBin<fNumUnits)
432             {
433               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
434               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
435               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
436               Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
437               Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
438               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
439             } 
440           else {
441             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
442               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
443               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
444               Deta = fTpcGrid->GetDeta();
445               Dphi = fTpcGrid->GetDphi();
446               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
447             }
448             else {
449               if(fDZ) {
450                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
451                   if(nBin<fNumUnits+nElements+n0)
452                     {
453                       Float_t phi = eta = 0.;
454                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
455                       Deta = fGrid0->GetDeta(); 
456                       Dphi = fGrid0->GetDphi(); 
457                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
458                     }
459                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
460                     {
461                       Float_t phi = eta = 0.;
462                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
463                       Deta = fGrid1->GetDeta(); 
464                       Dphi = fGrid1->GetDphi(); 
465                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
466                     }
467                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
468                     {
469                       Float_t phi = eta = 0.;
470                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
471                       Deta = fGrid2->GetDeta(); 
472                       Dphi = fGrid2->GetDphi(); 
473                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
474                     }
475                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
476                     {
477                       Float_t phi = eta = 0.;
478                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
479                       Deta = fGrid3->GetDeta(); 
480                       Dphi = fGrid3->GetDphi(); 
481                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
482                     }
483                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
484                     {
485                       Float_t phi = eta = 0.;
486                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
487                       Deta = fGrid4->GetDeta(); 
488                       Dphi = fGrid4->GetDphi(); 
489                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
490                     }
491                 }
492               } // end if(fDZ)
493             } // end else 2
494           } // end else 1
495         } // end loop on nBin
496     } // end grid type == 1
497   fArrayInitialised = 1;
498 }
499
500
501