1c7ade22de981b9f2bebec8be627b7f153f1c1b1
[u/mrichter/AliRoot.git] / JETAN / AliJetAODReader.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 AOD Reader 
18 // AOD reader for jet analysis
19 // This is the reader which must be used if the jet analysis task
20 // is executed after the ESD filter task, in order to read its output
21 //
22 // Author: Davide Perrino <davide.perrino@cern.ch>
23 //------------------------------------------------------------------------- 
24
25
26 #include <Riostream.h>
27 #include <TSystem.h>
28 #include <TLorentzVector.h>
29 #include <TVector3.h>
30 #include <TChain.h>
31 #include <TFile.h>
32 #include <TTask.h>
33 #include <TGeoManager.h>
34
35 #include "AliJetAODReader.h"
36 #include "AliJetAODReaderHeader.h"
37 #include "AliAODEvent.h"
38 #include "AliAODTrack.h"
39 #include "AliAODMCParticle.h"
40 #include "AliJetDummyGeo.h"
41 #include "AliJetAODFillUnitArrayTracks.h"
42 #include "AliJetAODFillUnitArrayEMCalDigits.h"
43 #include "AliJetHadronCorrectionv1.h"
44 #include "AliJetUnitArray.h"
45
46 ClassImp(AliJetAODReader)
47
48 AliJetAODReader::AliJetAODReader():
49     AliJetReader(),
50     fAOD(0x0),
51     fRef(new TRefArray),
52     fDebug(0),
53     fOpt(0),
54     fGeom(0),
55     fHadCorr(0x0),
56     fTpcGrid(0x0),
57     fEmcalGrid(0x0),
58     fGrid0(0),
59     fGrid1(0),
60     fGrid2(0),
61     fGrid3(0),
62     fGrid4(0),
63     fPtCut(0),
64     fApplyElectronCorrection(kFALSE),
65     fApplyMIPCorrection(kTRUE),
66     fApplyFractionHadronicCorrection(kFALSE),
67     fFractionHadronicCorrection(0.3),
68     fNumUnits(0),
69     fMass(0),
70     fSign(0),
71     fNIn(0),
72     fDZ(0),
73     fNeta(0),
74     fNphi(0),
75     fRefArray(0x0),
76     fProcId(kFALSE)
77 {
78   // Constructor    
79 }
80
81 //____________________________________________________________________________
82
83 AliJetAODReader::~AliJetAODReader()
84 {
85   // Destructor
86     delete fAOD;
87     delete fRef;
88     delete fTpcGrid;
89     delete fEmcalGrid;
90     if(fDZ)
91       {
92         delete fGrid0;
93         delete fGrid1;
94         delete fGrid2;
95         delete fGrid3;
96         delete fGrid4;
97       }
98
99 }
100
101 //____________________________________________________________________________
102
103 void AliJetAODReader::OpenInputFiles()
104 {
105     // Open the necessary input files
106     // chain for the AODs
107   fChain   = new TChain("aodTree");
108
109   // get directory and pattern name from the header
110    const char* dirName=fReaderHeader->GetDirectory();
111    const char* pattern=fReaderHeader->GetPattern();
112
113 //   // Add files matching patters to the chain
114
115    void *dir  = gSystem->OpenDirectory(dirName);
116    const char *name = 0x0;
117    int naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
118    int a = 0;
119    while ((name = gSystem->GetDirEntry(dir))){
120        if (a>=naod) continue;
121        
122        if (strstr(name,pattern)){
123            char path[256];
124            sprintf(path,"%s/%s/aod.root",dirName,name);
125            fChain->AddFile(path);
126            a++;
127        }
128    }
129   
130   gSystem->FreeDirectory(dir);
131   
132   fAOD = 0;
133   fChain->SetBranchAddress("AOD",&fAOD);
134   
135   int nMax = fChain->GetEntries(); 
136
137   printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
138   
139   // set number of events in header
140   if (fReaderHeader->GetLastEvent() == -1)
141     fReaderHeader->SetLastEvent(nMax);
142   else {
143     Int_t nUsr = fReaderHeader->GetLastEvent();
144     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
145   }
146 }
147
148 //____________________________________________________________________________
149
150 void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
151     // Connect the tree
152     // For AOD reader it's needed only to set the number of events
153      fChain = (TChain*)      tree;
154      
155      Int_t nMax = fChain->GetEntries(); 
156      printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
157      // set number of events in header
158      if (fReaderHeader->GetLastEvent() == -1)
159          fReaderHeader->SetLastEvent(nMax);
160      else {
161          Int_t nUsr = fReaderHeader->GetLastEvent();
162          fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
163      }
164 }
165
166
167 Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){
168
169   //
170   // filter for charge and for charged and neutral, no detector
171   // response yet
172   // physical priamries already selected
173
174   Int_t   pdg     = TMath::Abs(mcP->GetPdgCode());
175
176   // exclude neutrinos anyway   
177   if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
178
179   if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE;
180   if(flag==AliJetAODReaderHeader::kChargedMC){
181     if(mcP->Charge()==0)return kFALSE;
182     return kTRUE;
183   }
184
185   return kTRUE;
186
187 }
188  
189
190 Bool_t AliJetAODReader::FillMomentumArrayMC(){
191
192   // 
193   // This routine fetches the MC particles from the AOD
194   // Depending on the flag all particles except neurinos are use
195   // or only charged particles
196   //
197
198   TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
199   if(!mcArray){
200     Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
201     return kFALSE;
202   }
203   
204   Int_t nMC = mcArray->GetEntriesFast();
205   
206   // get number of tracks in event (for the loop)
207   if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC);
208   
209   // temporary storage of signal and pt cut flag
210   Int_t* sflag  = new Int_t[nMC];
211   Int_t* cflag  = new Int_t[nMC];
212   
213   // get cuts set by user
214   Float_t ptMin =  fReaderHeader->GetPtCut();
215   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
216   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
217   Int_t mcTrack = 0;
218   Float_t pt, eta;
219   TVector3 p3;
220
221   Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
222
223
224   for (Int_t it = 0; it < nMC; it++) {
225     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
226     if(!track->IsPhysicalPrimary())continue;
227
228     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
229     eta = p3.Eta();
230     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
231     if(!AcceptAODMCParticle(track,readerFlag))continue;
232     pt = p3.Pt();
233
234
235     
236
237     if (mcTrack == 0){
238       fRef->Delete(); // make sure to delete before placement new...
239       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
240     }
241     new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag());
242     sflag[mcTrack] = 1;
243     cflag[mcTrack] = ( pt > ptMin ) ? 1: 0;
244     fRef->Add(track);
245     mcTrack++;
246   }
247   if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack);
248   // set the signal flags
249   fSignalFlag.Set(mcTrack,sflag);
250   fCutFlag.Set(mcTrack,cflag);
251
252   delete [] sflag;
253   delete [] cflag;
254
255   return kTRUE;
256 }
257
258 //____________________________________________________________________________
259
260 Bool_t AliJetAODReader::FillMomentumArray()
261 {
262   // Clear momentum array
263   ClearArray();
264   fRef->Clear();
265   fDebug = fReaderHeader->GetDebug();
266
267   if (!fAOD) {
268       return kFALSE;
269   }
270
271   if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
272     return FillMomentumArrayMC();
273   }
274    
275   // get number of tracks in event (for the loop)
276   Int_t nt = fAOD->GetNTracks();
277   if(fDebug>0)printf("AOD tracks: %5d \t", nt);
278   
279   // temporary storage of signal and pt cut flag
280   Int_t* sflag  = new Int_t[nt];
281   Int_t* cflag  = new Int_t[nt];
282   
283   // get cuts set by user
284   Float_t ptMin =  fReaderHeader->GetPtCut();
285   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
286   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
287   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
288
289   //loop over tracks
290   Int_t aodTrack = 0;
291   Float_t pt, eta;
292   TVector3 p3;
293
294   for (Int_t it = 0; it < nt; it++) {
295     AliAODTrack *track = fAOD->GetTrack(it);
296     UInt_t status = track->GetStatus();
297     
298     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
299     p3.SetXYZ(mom[0],mom[1],mom[2]);
300     pt = p3.Pt();
301     eta = p3.Eta();
302     if (status == 0) continue;
303     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
304     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
305
306     if (aodTrack == 0){
307       fRef->Delete(); // make sure to delete before placement new...
308       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
309     }
310     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
311     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
312     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
313     aodTrack++;
314     fRef->Add(track);
315   }
316   if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack);
317   // set the signal flags
318   fSignalFlag.Set(aodTrack,sflag);
319   fCutFlag.Set(aodTrack,cflag);
320
321   delete [] sflag;
322   delete [] cflag;
323   
324   return kTRUE;
325 }
326
327 //__________________________________________________________
328 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
329 {
330   //
331   // Set flag to apply MIP correction fApplyMIPCorrection
332   // - exclusive with fApplyFractionHadronicCorrection
333   //
334
335   fApplyMIPCorrection = val;
336   if(fApplyMIPCorrection == kTRUE)
337     {
338       SetApplyFractionHadronicCorrection(kFALSE);
339       printf("Enabling MIP Correction \n");
340     }
341   else
342     {
343       printf("Disabling MIP Correction \n");
344     }
345 }
346
347 //__________________________________________________________
348 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
349 {
350   //
351   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
352   // - exclusive with fApplyMIPCorrection
353   //
354
355   fApplyFractionHadronicCorrection = val;
356   if(fApplyFractionHadronicCorrection == kTRUE)
357     {
358       SetApplyMIPCorrection(kFALSE);
359       printf("Enabling Fraction Hadronic Correction \n");
360     }
361   else
362     {
363       printf("Disabling Fraction Hadronic Correction \n");
364     }
365 }
366
367 //__________________________________________________________
368 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
369 {
370   //
371   // Set value to fFractionHadronicCorrection (default is 0.3)
372   // apply EMC hadronic correction fApplyFractionHadronicCorrection
373   // - exclusive with fApplyMIPCorrection
374   //
375
376   fFractionHadronicCorrection = val;
377   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
378     {
379       SetApplyFractionHadronicCorrection(kTRUE);
380       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
381     }
382   else
383     {
384       SetApplyFractionHadronicCorrection(kFALSE);
385     }
386 }
387
388 //____________________________________________________________________________
389 void AliJetAODReader::CreateTasks(TChain* tree)
390 {
391   //
392   // For reader task initialization
393   //
394
395   fDebug = fReaderHeader->GetDebug();
396   fDZ = fReaderHeader->GetDZ();
397   fTree = tree;
398
399   // Init EMCAL geometry and create UnitArray object
400   SetEMCALGeometry();
401   //  cout << "In create task" << endl;
402   InitParameters();
403   InitUnitArray();
404
405   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
406   fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
407   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
408   fFillUAFromTracks->SetGeom(fGeom);
409   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
410   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
411
412   if(fDZ)
413     {
414       fFillUAFromTracks->SetGrid0(fGrid0);
415       fFillUAFromTracks->SetGrid1(fGrid1);
416       fFillUAFromTracks->SetGrid2(fGrid2);
417       fFillUAFromTracks->SetGrid3(fGrid3);
418       fFillUAFromTracks->SetGrid4(fGrid4);
419     }
420   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
421   fFillUAFromTracks->SetHadCorrector(fHadCorr);
422   fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
423   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
424   fFillUAFromEMCalDigits->SetGeom(fGeom);
425   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
426   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
427   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
428   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
429   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
430   // Add the task to global task
431   fFillUnitArray->Add(fFillUAFromTracks);
432   fFillUnitArray->Add(fFillUAFromEMCalDigits);
433   fFillUAFromTracks->SetActive(kFALSE);
434   fFillUAFromEMCalDigits->SetActive(kFALSE);
435
436   cout << "Tasks instantiated at that stage ! " << endl;
437   cout << "You can loop over events now ! " << endl;
438
439 }
440
441 //____________________________________________________________________________
442 Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
443 {
444   //
445   // Main function
446   // Fill the reader part
447   //
448
449   fProcId = procid;
450   fRefArray = refArray;
451 //(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
452 //(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
453
454   // clear momentum array
455   ClearArray();
456   fRef->Clear();
457
458   fDebug = fReaderHeader->GetDebug();
459   fOpt = fReaderHeader->GetDetector();
460
461   if(!fAOD) {
462     return kFALSE;
463   }
464
465   // TPC only or Digits+TPC or Clusters+TPC
466   if(fOpt%2==!0 && fOpt!=0){
467     fFillUAFromTracks->SetAOD(fAOD);
468     fFillUAFromTracks->SetActive(kTRUE);
469     fFillUAFromTracks->SetUnitArray(fUnitArray);
470     fFillUAFromTracks->SetRefArray(fRefArray);
471     fFillUAFromTracks->SetReferences(fRef);
472     fFillUAFromTracks->SetSignalFlag(fSignalFlag);
473     fFillUAFromTracks->SetCutFlag(fCutFlag);
474     fFillUAFromTracks->SetProcId(fProcId);
475     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
476     fFillUAFromTracks->Exec("tpc");
477     if(fOpt==1){
478       fNumCandidate = fFillUAFromTracks->GetMult();
479       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
480     }
481     fSignalFlag = fFillUAFromTracks->GetSignalFlag();
482     fCutFlag = fFillUAFromTracks->GetCutFlag();
483   }
484
485   // Digits only or Digits+TPC
486   if(fOpt>=2 && fOpt<=3){
487     fFillUAFromEMCalDigits->SetAOD(fAOD);
488     fFillUAFromEMCalDigits->SetActive(kTRUE);
489     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
490     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
491     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
492     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
493     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
494     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
495     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
496     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
497   }
498
499   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
500
501   return kTRUE;
502 }
503
504 //____________________________________________________________________________
505 Bool_t AliJetAODReader::SetEMCALGeometry()
506 {
507   // 
508   // Set the EMCal Geometry
509   //
510   
511   if (!fTree->GetFile()) 
512     return kFALSE;
513
514   TString geomFile(fTree->GetFile()->GetName());
515   geomFile.ReplaceAll("AliESDs", "geometry");
516   
517   // temporary workaround for PROOF bug #18505
518   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
519   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
520
521   // Define EMCAL geometry to be able to read ESDs
522   fGeom = AliJetDummyGeo::GetInstance();
523   if (fGeom == 0)
524     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
525   
526   // To be setted to run some AliEMCALGeometry functions
527   TGeoManager::Import(geomFile);
528   fGeom->GetTransformationForSM();  
529   printf("\n EMCal Geometry set ! \n");
530   
531   return kTRUE;
532   
533 }
534
535 //____________________________________________________________________________
536 void AliJetAODReader::InitParameters()
537 {
538   // Initialise parameters
539   fOpt = fReaderHeader->GetDetector();
540   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
541   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
542 }
543
544 //____________________________________________________________________________
545 void AliJetAODReader::InitUnitArray()
546 {
547   //Initialises unit arrays
548   Int_t nElements = fTpcGrid->GetNEntries();
549   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
550   if(fArrayInitialised) fUnitArray->Delete();
551
552   if(fTpcGrid->GetGridType()==0)
553     { // Fill the following quantities :
554       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
555       // detector flag, in/out jet, pt cut, mass, cluster ID)
556       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
557         {
558           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
559           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
560           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
561           deltaEta = fTpcGrid->GetDeta();
562           deltaPhi = fTpcGrid->GetDphi();
563           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
564         }
565     }
566
567   if(fTpcGrid->GetGridType()==1)
568     {
569       Int_t nGaps = 0;
570       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
571
572       if(fDZ)
573         {
574           // Define a grid of cell for the gaps between SM
575           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
576           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
577           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
578           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
579           fGrid0->SetGridType(0);
580           fGrid0->SetMatrixIndexes();
581           fGrid0->SetIndexIJ();
582           n0 = fGrid0->GetNEntries();
583           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
584           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
585           fGrid1->SetGridType(0);
586           fGrid1->SetMatrixIndexes();
587           fGrid1->SetIndexIJ();
588           n1 = fGrid1->GetNEntries();
589           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
590           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
591           fGrid2->SetGridType(0);
592           fGrid2->SetMatrixIndexes();
593           fGrid2->SetIndexIJ();
594           n2 = fGrid2->GetNEntries();
595           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
596           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
597           fGrid3->SetGridType(0);
598           fGrid3->SetMatrixIndexes();
599           fGrid3->SetIndexIJ();
600           n3 = fGrid3->GetNEntries();
601           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
602           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
603           fGrid4->SetGridType(0);
604           fGrid4->SetMatrixIndexes();
605           fGrid4->SetIndexIJ();
606           n4 = fGrid4->GetNEntries();
607
608           nGaps = n0+n1+n2+n3+n4;
609
610         }
611
612       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
613         {
614           if(nBin<fNumUnits)
615             {
616               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
617               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
618               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
619               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
620               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
621               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
622             } 
623           else {
624             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
625               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
626               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
627               deltaEta = fTpcGrid->GetDeta();
628               deltaPhi = fTpcGrid->GetDphi();
629               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
630             }
631             else {
632               if(fDZ) {
633                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
634                   if(nBin<fNumUnits+nElements+n0)
635                     {
636                       phi = eta = 0.;
637                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
638                       deltaEta = fGrid0->GetDeta(); 
639                       deltaPhi = fGrid0->GetDphi(); 
640                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
641                     }
642                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
643                     {
644                       phi = eta = 0.;
645                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
646                       deltaEta = fGrid1->GetDeta(); 
647                       deltaPhi = fGrid1->GetDphi(); 
648                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
649                     }
650                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
651                     {
652                       phi = eta = 0.;
653                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
654                       deltaEta = fGrid2->GetDeta(); 
655                       deltaPhi = fGrid2->GetDphi(); 
656                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
657                     }
658                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
659                     {
660                       phi = eta = 0.;
661                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
662                       deltaEta = fGrid3->GetDeta(); 
663                       deltaPhi = fGrid3->GetDphi(); 
664                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
665                     }
666                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
667                     {
668                       phi = eta = 0.;
669                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
670                       deltaEta = fGrid4->GetDeta(); 
671                       deltaPhi = fGrid4->GetDphi(); 
672                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
673                     }
674                 }
675               } // end if(fDZ)
676             } // end else 2
677           } // end else 1
678         } // end loop on nBin
679     } // end grid type == 1
680   fArrayInitialised = 1;
681 }
682