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