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