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