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