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