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