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