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