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