Coding violations corrected.
[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) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
238     new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag());
239     sflag[mcTrack] = 1;
240     cflag[mcTrack] = ( pt > ptMin ) ? 1: 0;
241     fRef->Add(track);
242     mcTrack++;
243   }
244   if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack);
245   // set the signal flags
246   fSignalFlag.Set(mcTrack,sflag);
247   fCutFlag.Set(mcTrack,cflag);
248
249   delete [] sflag;
250   delete [] cflag;
251
252   return kTRUE;
253 }
254
255 //____________________________________________________________________________
256
257 Bool_t AliJetAODReader::FillMomentumArray()
258 {
259   // Clear momentum array
260   ClearArray();
261   fRef->Clear();
262   fDebug = fReaderHeader->GetDebug();
263
264   if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
265     return FillMomentumArrayMC();
266   }
267   
268
269   
270   if (!fAOD) {
271       return kFALSE;
272   }
273   
274   // get number of tracks in event (for the loop)
275   Int_t nt = fAOD->GetNTracks();
276   if(fDebug>0)printf("AOD tracks: %5d \t", nt);
277   
278   // temporary storage of signal and pt cut flag
279   Int_t* sflag  = new Int_t[nt];
280   Int_t* cflag  = new Int_t[nt];
281   
282   // get cuts set by user
283   Float_t ptMin =  fReaderHeader->GetPtCut();
284   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
285   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
286   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
287
288   //loop over tracks
289   Int_t aodTrack = 0;
290   Float_t pt, eta;
291   TVector3 p3;
292
293   for (Int_t it = 0; it < nt; it++) {
294     AliAODTrack *track = fAOD->GetTrack(it);
295     UInt_t status = track->GetStatus();
296     
297     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
298     p3.SetXYZ(mom[0],mom[1],mom[2]);
299     pt = p3.Pt();
300     eta = p3.Eta();
301     if (status == 0) continue;
302     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
303     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
304
305     if (aodTrack == 0){
306       fRef->Delete(); // make sure to delete before placement new...
307       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
308     }
309     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
310     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
311     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
312     aodTrack++;
313     fRef->Add(track);
314   }
315   if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack);
316   // set the signal flags
317   fSignalFlag.Set(aodTrack,sflag);
318   fCutFlag.Set(aodTrack,cflag);
319
320   delete [] sflag;
321   delete [] cflag;
322   
323   return kTRUE;
324 }
325
326 //__________________________________________________________
327 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
328 {
329   //
330   // Set flag to apply MIP correction fApplyMIPCorrection
331   // - exclusive with fApplyFractionHadronicCorrection
332   //
333
334   fApplyMIPCorrection = val;
335   if(fApplyMIPCorrection == kTRUE)
336     {
337       SetApplyFractionHadronicCorrection(kFALSE);
338       printf("Enabling MIP Correction \n");
339     }
340   else
341     {
342       printf("Disabling MIP Correction \n");
343     }
344 }
345
346 //__________________________________________________________
347 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
348 {
349   //
350   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
351   // - exclusive with fApplyMIPCorrection
352   //
353
354   fApplyFractionHadronicCorrection = val;
355   if(fApplyFractionHadronicCorrection == kTRUE)
356     {
357       SetApplyMIPCorrection(kFALSE);
358       printf("Enabling Fraction Hadronic Correction \n");
359     }
360   else
361     {
362       printf("Disabling Fraction Hadronic Correction \n");
363     }
364 }
365
366 //__________________________________________________________
367 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
368 {
369   //
370   // Set value to fFractionHadronicCorrection (default is 0.3)
371   // apply EMC hadronic correction fApplyFractionHadronicCorrection
372   // - exclusive with fApplyMIPCorrection
373   //
374
375   fFractionHadronicCorrection = val;
376   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
377     {
378       SetApplyFractionHadronicCorrection(kTRUE);
379       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
380     }
381   else
382     {
383       SetApplyFractionHadronicCorrection(kFALSE);
384     }
385 }
386
387 //____________________________________________________________________________
388 void AliJetAODReader::CreateTasks(TChain* tree)
389 {
390   //
391   // For reader task initialization
392   //
393
394   fDebug = fReaderHeader->GetDebug();
395   fDZ = fReaderHeader->GetDZ();
396   fTree = tree;
397
398   // Init EMCAL geometry and create UnitArray object
399   SetEMCALGeometry();
400   //  cout << "In create task" << endl;
401   InitParameters();
402   InitUnitArray();
403
404   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
405   fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
406   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
407   fFillUAFromTracks->SetGeom(fGeom);
408   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
409   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
410
411   if(fDZ)
412     {
413       fFillUAFromTracks->SetGrid0(fGrid0);
414       fFillUAFromTracks->SetGrid1(fGrid1);
415       fFillUAFromTracks->SetGrid2(fGrid2);
416       fFillUAFromTracks->SetGrid3(fGrid3);
417       fFillUAFromTracks->SetGrid4(fGrid4);
418     }
419   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
420   fFillUAFromTracks->SetHadCorrector(fHadCorr);
421   fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
422   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
423   fFillUAFromEMCalDigits->SetGeom(fGeom);
424   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
425   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
426   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
427   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
428   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
429   // Add the task to global task
430   fFillUnitArray->Add(fFillUAFromTracks);
431   fFillUnitArray->Add(fFillUAFromEMCalDigits);
432   fFillUAFromTracks->SetActive(kFALSE);
433   fFillUAFromEMCalDigits->SetActive(kFALSE);
434
435   cout << "Tasks instantiated at that stage ! " << endl;
436   cout << "You can loop over events now ! " << endl;
437
438 }
439
440 //____________________________________________________________________________
441 Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
442 {
443   //
444   // Main function
445   // Fill the reader part
446   //
447
448   fProcId = procid;
449   fRefArray = refArray;
450 //(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
451 //(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
452
453   // clear momentum array
454   ClearArray();
455   fRef->Clear();
456
457   fDebug = fReaderHeader->GetDebug();
458   fOpt = fReaderHeader->GetDetector();
459
460   if(!fAOD) {
461     return kFALSE;
462   }
463
464   // TPC only or Digits+TPC or Clusters+TPC
465   if(fOpt%2==!0 && fOpt!=0){
466     fFillUAFromTracks->SetAOD(fAOD);
467     fFillUAFromTracks->SetActive(kTRUE);
468     fFillUAFromTracks->SetUnitArray(fUnitArray);
469     fFillUAFromTracks->SetRefArray(fRefArray);
470     fFillUAFromTracks->SetReferences(fRef);
471     fFillUAFromTracks->SetSignalFlag(fSignalFlag);
472     fFillUAFromTracks->SetCutFlag(fCutFlag);
473     fFillUAFromTracks->SetProcId(fProcId);
474     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
475     fFillUAFromTracks->Exec("tpc");
476     if(fOpt==1){
477       fNumCandidate = fFillUAFromTracks->GetMult();
478       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
479     }
480     fSignalFlag = fFillUAFromTracks->GetSignalFlag();
481     fCutFlag = fFillUAFromTracks->GetCutFlag();
482   }
483
484   // Digits only or Digits+TPC
485   if(fOpt>=2 && fOpt<=3){
486     fFillUAFromEMCalDigits->SetAOD(fAOD);
487     fFillUAFromEMCalDigits->SetActive(kTRUE);
488     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
489     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
490     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
491     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
492     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
493     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
494     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
495     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
496   }
497
498   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
499
500   return kTRUE;
501 }
502
503 //____________________________________________________________________________
504 Bool_t AliJetAODReader::SetEMCALGeometry()
505 {
506   // 
507   // Set the EMCal Geometry
508   //
509   
510   if (!fTree->GetFile()) 
511     return kFALSE;
512
513   TString geomFile(fTree->GetFile()->GetName());
514   geomFile.ReplaceAll("AliESDs", "geometry");
515   
516   // temporary workaround for PROOF bug #18505
517   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
518   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
519
520   // Define EMCAL geometry to be able to read ESDs
521   fGeom = AliJetDummyGeo::GetInstance();
522   if (fGeom == 0)
523     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
524   
525   // To be setted to run some AliEMCALGeometry functions
526   TGeoManager::Import(geomFile);
527   fGeom->GetTransformationForSM();  
528   printf("\n EMCal Geometry set ! \n");
529   
530   return kTRUE;
531   
532 }
533
534 //____________________________________________________________________________
535 void AliJetAODReader::InitParameters()
536 {
537   // Initialise parameters
538   fOpt = fReaderHeader->GetDetector();
539   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
540   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
541 }
542
543 //____________________________________________________________________________
544 void AliJetAODReader::InitUnitArray()
545 {
546   //Initialises unit arrays
547   Int_t nElements = fTpcGrid->GetNEntries();
548   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
549   if(fArrayInitialised) fUnitArray->Delete();
550
551   if(fTpcGrid->GetGridType()==0)
552     { // Fill the following quantities :
553       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
554       // detector flag, in/out jet, pt cut, mass, cluster ID)
555       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
556         {
557           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
558           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
559           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
560           deltaEta = fTpcGrid->GetDeta();
561           deltaPhi = fTpcGrid->GetDphi();
562           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
563         }
564     }
565
566   if(fTpcGrid->GetGridType()==1)
567     {
568       Int_t nGaps = 0;
569       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
570
571       if(fDZ)
572         {
573           // Define a grid of cell for the gaps between SM
574           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
575           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
576           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
577           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
578           fGrid0->SetGridType(0);
579           fGrid0->SetMatrixIndexes();
580           fGrid0->SetIndexIJ();
581           n0 = fGrid0->GetNEntries();
582           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
583           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
584           fGrid1->SetGridType(0);
585           fGrid1->SetMatrixIndexes();
586           fGrid1->SetIndexIJ();
587           n1 = fGrid1->GetNEntries();
588           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
589           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
590           fGrid2->SetGridType(0);
591           fGrid2->SetMatrixIndexes();
592           fGrid2->SetIndexIJ();
593           n2 = fGrid2->GetNEntries();
594           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
595           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
596           fGrid3->SetGridType(0);
597           fGrid3->SetMatrixIndexes();
598           fGrid3->SetIndexIJ();
599           n3 = fGrid3->GetNEntries();
600           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
601           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
602           fGrid4->SetGridType(0);
603           fGrid4->SetMatrixIndexes();
604           fGrid4->SetIndexIJ();
605           n4 = fGrid4->GetNEntries();
606
607           nGaps = n0+n1+n2+n3+n4;
608
609         }
610
611       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++) 
612         {
613           if(nBin<fNumUnits)
614             {
615               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry 
616               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
617               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
618               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
619               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
620               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
621             } 
622           else {
623             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
624               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
625               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
626               deltaEta = fTpcGrid->GetDeta();
627               deltaPhi = fTpcGrid->GetDphi();
628               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
629             }
630             else {
631               if(fDZ) {
632                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
633                   if(nBin<fNumUnits+nElements+n0)
634                     {
635                       phi = eta = 0.;
636                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
637                       deltaEta = fGrid0->GetDeta(); 
638                       deltaPhi = fGrid0->GetDphi(); 
639                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
640                     }
641                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
642                     {
643                       phi = eta = 0.;
644                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
645                       deltaEta = fGrid1->GetDeta(); 
646                       deltaPhi = fGrid1->GetDphi(); 
647                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
648                     }
649                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
650                     {
651                       phi = eta = 0.;
652                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
653                       deltaEta = fGrid2->GetDeta(); 
654                       deltaPhi = fGrid2->GetDphi(); 
655                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
656                     }
657                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
658                     {
659                       phi = eta = 0.;
660                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
661                       deltaEta = fGrid3->GetDeta(); 
662                       deltaPhi = fGrid3->GetDphi(); 
663                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
664                     }
665                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
666                     {
667                       phi = eta = 0.;
668                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
669                       deltaEta = fGrid4->GetDeta(); 
670                       deltaPhi = fGrid4->GetDphi(); 
671                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
672                     }
673                 }
674               } // end if(fDZ)
675             } // end else 2
676           } // end else 1
677         } // end loop on nBin
678     } // end grid type == 1
679   fArrayInitialised = 1;
680 }
681