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