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