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