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     fAOD(0x0),
  50     fRef(new TRefArray),
  51     fDebug(0),
  52     fOpt(0),
  53     fGeom(0),
  54     fHadCorr(0x0),
  55     fTpcGrid(0x0),
  56     fEmcalGrid(0x0),
  57     fGrid0(0),
  58     fGrid1(0),
  59     fGrid2(0),
  60     fGrid3(0),
  61     fGrid4(0),
  62     fPtCut(0),
  63     fApplyElectronCorrection(kFALSE),
  64     fApplyMIPCorrection(kTRUE),
  65     fApplyFractionHadronicCorrection(kFALSE),
  66     fFractionHadronicCorrection(0.3),
  67     fNumUnits(0),
  68     fMass(0),
  69     fSign(0),
  70     fNIn(0),
  71     fDZ(0),
  72     fNeta(0),
  73     fNphi(0),
  74     fRefArray(0x0),
  75     fProcId(kFALSE)
  76 {
  77   // Constructor
  78 }
  79
  80 //____________________________________________________________________________
  81
  82 AliJetAODReader::~AliJetAODReader()
  83 {
  84   // Destructor
  85     delete fAOD;
  86     delete fRef;
  87     delete fTpcGrid;
  88     delete fEmcalGrid;
  89     if(fDZ)
  90       {
  91         delete fGrid0;
  92         delete fGrid1;
  93         delete fGrid2;
  94         delete fGrid3;
  95         delete fGrid4;
  96       }
  97
  98 }
  99
 100 //____________________________________________________________________________
 101
 102 void AliJetAODReader::OpenInputFiles()
 103 {
 104     // Open the necessary input files
 105     // chain for the AODs
 106   fChain   = new TChain("aodTree");
 107
 108   // get directory and pattern name from the header
 109    const char* dirName=fReaderHeader->GetDirectory();
 110    const char* pattern=fReaderHeader->GetPattern();
 111
 112 //   // Add files matching patters to the chain
 113
 114    void *dir  = gSystem->OpenDirectory(dirName);
 115    const char *name = 0x0;
 116    int naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
 117    int a = 0;
 118    while ((name = gSystem->GetDirEntry(dir))){
 119        if (a>=naod) continue;
 120
 121        if (strstr(name,pattern)){
 122            char path[256];
 123            sprintf(path,"%s/%s/aod.root",dirName,name);
 124            fChain->AddFile(path);
 125            a++;
 126        }
 127    }
 128
 129   gSystem->FreeDirectory(dir);
 130
 131   fAOD = 0;
 132   fChain->SetBranchAddress("AOD",&fAOD);
 133
 134   int nMax = fChain->GetEntries();
 135
 136   printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
 137
 138   // set number of events in header
 139   if (fReaderHeader->GetLastEvent() == -1)
 140     fReaderHeader->SetLastEvent(nMax);
 141   else {
 142     Int_t nUsr = fReaderHeader->GetLastEvent();
 143     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 144   }
 145 }
 146
 147 //____________________________________________________________________________
 148
 149 void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
 150     // Connect the tree
 151     // For AOD reader it's needed only to set the number of events
 152      fChain = (TChain*)      tree;
 153
 154      Int_t nMax = fChain->GetEntries();
 155      printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
 156      // set number of events in header
 157      if (fReaderHeader->GetLastEvent() == -1)
 158          fReaderHeader->SetLastEvent(nMax);
 159      else {
 160          Int_t nUsr = fReaderHeader->GetLastEvent();
 161          fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 162      }
 163 }
 164
 165 //____________________________________________________________________________
 166
 167 Bool_t AliJetAODReader::FillMomentumArray()
 168 {
 169   // Clear momentum array
 170   ClearArray();
 171   fRef->Clear();
 172
 173   fDebug = fReaderHeader->GetDebug();
 174
 175   if (!fAOD) {
 176       return kFALSE;
 177   }
 178
 179   // get number of tracks in event (for the loop)
 180   Int_t nt = fAOD->GetNTracks();
 181   printf("AOD tracks: %5d \t", 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   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
 192
 193   //loop over tracks
 194   Int_t aodTrack = 0;
 195   Float_t pt, eta;
 196   TVector3 p3;
 197
 198   for (Int_t it = 0; it < nt; it++) {
 199     AliAODTrack *track = fAOD->GetTrack(it);
 200     UInt_t status = track->GetStatus();
 201
 202     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 203     p3.SetXYZ(mom[0],mom[1],mom[2]);
 204     pt = p3.Pt();
 205     eta = p3.Eta();
 206     if (status == 0) continue;
 207     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
 208     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 209
 210     if (aodTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 211
 212     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 213     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
 214     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
 215     aodTrack++;
 216     fRef->Add(track);
 217   }
 218   printf("Used AOD tracks: %5d \n", aodTrack);
 219   // set the signal flags
 220   fSignalFlag.Set(aodTrack,sflag);
 221   fCutFlag.Set(aodTrack,cflag);
 222
 223   delete [] sflag;
 224   delete [] cflag;
 225
 226   return kTRUE;
 227 }
 228
 229 //__________________________________________________________
 230 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
 231 {
 232   //
 233   // Set flag to apply MIP correction fApplyMIPCorrection
 234   // - exclusive with fApplyFractionHadronicCorrection
 235   //
 236
 237   fApplyMIPCorrection = val;
 238   if(fApplyMIPCorrection == kTRUE)
 239     {
 240       SetApplyFractionHadronicCorrection(kFALSE);
 241       printf("Enabling MIP Correction \n");
 242     }
 243   else
 244     {
 245       printf("Disabling MIP Correction \n");
 246     }
 247 }
 248
 249 //__________________________________________________________
 250 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
 251 {
 252   //
 253   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
 254   // - exclusive with fApplyMIPCorrection
 255   //
 256
 257   fApplyFractionHadronicCorrection = val;
 258   if(fApplyFractionHadronicCorrection == kTRUE)
 259     {
 260       SetApplyMIPCorrection(kFALSE);
 261       printf("Enabling Fraction Hadronic Correction \n");
 262     }
 263   else
 264     {
 265       printf("Disabling Fraction Hadronic Correction \n");
 266     }
 267 }
 268
 269 //__________________________________________________________
 270 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
 271 {
 272   //
 273   // Set value to fFractionHadronicCorrection (default is 0.3)
 274   // apply EMC hadronic correction fApplyFractionHadronicCorrection
 275   // - exclusive with fApplyMIPCorrection
 276   //
 277
 278   fFractionHadronicCorrection = val;
 279   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
 280     {
 281       SetApplyFractionHadronicCorrection(kTRUE);
 282       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
 283     }
 284   else
 285     {
 286       SetApplyFractionHadronicCorrection(kFALSE);
 287     }
 288 }
 289
 290 //____________________________________________________________________________
 291 void AliJetAODReader::CreateTasks(TChain* tree)
 292 {
 293   //
 294   // For reader task initialization
 295   //
 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 }