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