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