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   fDebug = fReaderHeader->GetDebug();
 173
 174   if (!fAOD) {
 175       return kFALSE;
 176   }
 177
 178   // get number of tracks in event (for the loop)
 179   Int_t nt = fAOD->GetNTracks();
 180   printf("AOD tracks: %5d \t", nt);
 181
 182   // temporary storage of signal and pt cut flag
 183   Int_t* sflag  = new Int_t[nt];
 184   Int_t* cflag  = new Int_t[nt];
 185
 186   // get cuts set by user
 187   Float_t ptMin =  fReaderHeader->GetPtCut();
 188   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 189   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 190   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
 191
 192   //loop over tracks
 193   Int_t aodTrack = 0;
 194   Float_t pt, eta;
 195   TVector3 p3;
 196
 197   for (Int_t it = 0; it < nt; it++) {
 198     AliAODTrack *track = fAOD->GetTrack(it);
 199     UInt_t status = track->GetStatus();
 200
 201     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 202     p3.SetXYZ(mom[0],mom[1],mom[2]);
 203     pt = p3.Pt();
 204     eta = p3.Eta();
 205     if (status == 0) continue;
 206     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
 207     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 208
 209     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 210     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
 211     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
 212     aodTrack++;
 213     fRef->Add(track);
 214   }
 215   printf("Used AOD tracks: %5d \n", aodTrack);
 216   // set the signal flags
 217   fSignalFlag.Set(aodTrack,sflag);
 218   fCutFlag.Set(aodTrack,cflag);
 219
 220   delete [] sflag;
 221   delete [] cflag;
 222
 223   return kTRUE;
 224 }
 225
 226 //__________________________________________________________
 227 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
 228 {
 229   //
 230   // Set flag to apply MIP correction fApplyMIPCorrection
 231   // - exclusive with fApplyFractionHadronicCorrection
 232   //
 233
 234   fApplyMIPCorrection = val;
 235   if(fApplyMIPCorrection == kTRUE)
 236     {
 237       SetApplyFractionHadronicCorrection(kFALSE);
 238       printf("Enabling MIP Correction \n");
 239     }
 240   else
 241     {
 242       printf("Disabling MIP Correction \n");
 243     }
 244 }
 245
 246 //__________________________________________________________
 247 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
 248 {
 249   //
 250   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
 251   // - exclusive with fApplyMIPCorrection
 252   //
 253
 254   fApplyFractionHadronicCorrection = val;
 255   if(fApplyFractionHadronicCorrection == kTRUE)
 256     {
 257       SetApplyMIPCorrection(kFALSE);
 258       printf("Enabling Fraction Hadronic Correction \n");
 259     }
 260   else
 261     {
 262       printf("Disabling Fraction Hadronic Correction \n");
 263     }
 264 }
 265
 266 //__________________________________________________________
 267 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
 268 {
 269   //
 270   // Set value to fFractionHadronicCorrection (default is 0.3)
 271   // apply EMC hadronic correction fApplyFractionHadronicCorrection
 272   // - exclusive with fApplyMIPCorrection
 273   //
 274
 275   fFractionHadronicCorrection = val;
 276   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
 277     {
 278       SetApplyFractionHadronicCorrection(kTRUE);
 279       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
 280     }
 281   else
 282     {
 283       SetApplyFractionHadronicCorrection(kFALSE);
 284     }
 285 }
 286
 287 //____________________________________________________________________________
 288 void AliJetAODReader::CreateTasks(TChain* tree)
 289 {
 290   fDebug = fReaderHeader->GetDebug();
 291   fDZ = fReaderHeader->GetDZ();
 292   fTree = tree;
 293
 294   // Init EMCAL geometry and create UnitArray object
 295   SetEMCALGeometry();
 296   //  cout << "In create task" << endl;
 297   InitParameters();
 298   InitUnitArray();
 299
 300   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
 301   fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
 302   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
 303   fFillUAFromTracks->SetGeom(fGeom);
 304   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
 305   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 306
 307   if(fDZ)
 308     {
 309       fFillUAFromTracks->SetGrid0(fGrid0);
 310       fFillUAFromTracks->SetGrid1(fGrid1);
 311       fFillUAFromTracks->SetGrid2(fGrid2);
 312       fFillUAFromTracks->SetGrid3(fGrid3);
 313       fFillUAFromTracks->SetGrid4(fGrid4);
 314     }
 315   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
 316   fFillUAFromTracks->SetHadCorrector(fHadCorr);
 317   fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
 318   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
 319   fFillUAFromEMCalDigits->SetGeom(fGeom);
 320   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
 321   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
 322   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
 323   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
 324   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
 325
 326   fFillUnitArray->Add(fFillUAFromTracks);
 327   fFillUnitArray->Add(fFillUAFromEMCalDigits);
 328   fFillUAFromTracks->SetActive(kFALSE);
 329   fFillUAFromEMCalDigits->SetActive(kFALSE);
 330
 331   cout << "Tasks instantiated at that stage ! " << endl;
 332   cout << "You can loop over events now ! " << endl;
 333
 334 }
 335
 336 //____________________________________________________________________________
 337 Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 338 {
 339   //
 340   // Main function
 341   // Fill the reader part
 342   //
 343
 344   fProcId = procid;
 345   fRefArray = refArray;
 346 //(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
 347 //(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
 348
 349   // clear momentum array
 350   ClearArray();
 351
 352   fDebug = fReaderHeader->GetDebug();
 353   fOpt = fReaderHeader->GetDetector();
 354
 355   if(!fAOD) {
 356     return kFALSE;
 357   }
 358
 359   // TPC only or Digits+TPC or Clusters+TPC
 360   if(fOpt%2==!0 && fOpt!=0){
 361     fFillUAFromTracks->SetAOD(fAOD);
 362     fFillUAFromTracks->SetActive(kTRUE);
 363     fFillUAFromTracks->SetUnitArray(fUnitArray);
 364     fFillUAFromTracks->SetRefArray(fRefArray);
 365     fFillUAFromTracks->SetProcId(fProcId);
 366     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
 367     fFillUAFromTracks->Exec("tpc");
 368     if(fOpt==1){
 369       fNumCandidate = fFillUAFromTracks->GetMult();
 370       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 371     }
 372   }
 373
 374   // Digits only or Digits+TPC
 375   if(fOpt>=2 && fOpt<=3){
 376     fFillUAFromEMCalDigits->SetAOD(fAOD);
 377     fFillUAFromEMCalDigits->SetActive(kTRUE);
 378     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 379     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 380     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
 381     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 382     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 383     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
 384     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 385     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 386   }
 387
 388   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
 389
 390   return kTRUE;
 391 }
 392
 393 //____________________________________________________________________________
 394 Bool_t AliJetAODReader::SetEMCALGeometry()
 395 {
 396   //
 397   // Set the EMCal Geometry
 398   //
 399
 400   if (!fTree->GetFile())
 401     return kFALSE;
 402
 403   TString geomFile(fTree->GetFile()->GetName());
 404   geomFile.ReplaceAll("AliESDs", "geometry");
 405
 406   // temporary workaround for PROOF bug #18505
 407   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
 408   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
 409
 410   // Define EMCAL geometry to be able to read ESDs
 411   fGeom = AliJetDummyGeo::GetInstance();
 412   if (fGeom == 0)
 413     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
 414
 415   // To be setted to run some AliEMCALGeometry functions
 416   TGeoManager::Import(geomFile);
 417   fGeom->GetTransformationForSM();
 418   printf("\n EMCal Geometry set ! \n");
 419
 420   return kTRUE;
 421
 422 }
 423
 424 //____________________________________________________________________________
 425 void AliJetAODReader::InitParameters()
 426 {
 427   // Initialise parameters
 428   fOpt = fReaderHeader->GetDetector();
 429   //  fHCorrection    = 0;          // For hadron correction
 430   fHadCorr        = 0;          // For hadron correction
 431   if(fEFlag==kFALSE){
 432     if(fOpt==0 || fOpt==1)
 433       fECorrection    = 0;        // For electron correction
 434     else fECorrection = 1;        // For electron correction
 435   }
 436   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 437   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 438 }
 439
 440 //____________________________________________________________________________
 441 void AliJetAODReader::InitUnitArray()
 442 {
 443   //Initialises unit arrays
 444   Int_t nElements = fTpcGrid->GetNEntries();
 445   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
 446   if(fArrayInitialised) fUnitArray->Delete();
 447
 448   if(fTpcGrid->GetGridType()==0)
 449     { // Fill the following quantities :
 450       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
 451       // detector flag, in/out jet, pt cut, mass, cluster ID)
 452       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 453         {
 454           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
 455           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 456           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 457           deltaEta = fTpcGrid->GetDeta();
 458           deltaPhi = fTpcGrid->GetDphi();
 459           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 460         }
 461     }
 462
 463   if(fTpcGrid->GetGridType()==1)
 464     {
 465       Int_t nGaps = 0;
 466       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 467
 468       if(fDZ)
 469         {
 470           // Define a grid of cell for the gaps between SM
 471           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 472           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 473           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 474           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 475           fGrid0->SetGridType(0);
 476           fGrid0->SetMatrixIndexes();
 477           fGrid0->SetIndexIJ();
 478           n0 = fGrid0->GetNEntries();
 479           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 480           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 481           fGrid1->SetGridType(0);
 482           fGrid1->SetMatrixIndexes();
 483           fGrid1->SetIndexIJ();
 484           n1 = fGrid1->GetNEntries();
 485           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 486           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 487           fGrid2->SetGridType(0);
 488           fGrid2->SetMatrixIndexes();
 489           fGrid2->SetIndexIJ();
 490           n2 = fGrid2->GetNEntries();
 491           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 492           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 493           fGrid3->SetGridType(0);
 494           fGrid3->SetMatrixIndexes();
 495           fGrid3->SetIndexIJ();
 496           n3 = fGrid3->GetNEntries();
 497           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 498           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 499           fGrid4->SetGridType(0);
 500           fGrid4->SetMatrixIndexes();
 501           fGrid4->SetIndexIJ();
 502           n4 = fGrid4->GetNEntries();
 503
 504           nGaps = n0+n1+n2+n3+n4;
 505
 506         }
 507
 508       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 509         {
 510           if(nBin<fNumUnits)
 511             {
 512               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 513               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 514               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 515               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 516               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 517               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 518             }
 519           else {
 520             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 521               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 522               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 523               deltaEta = fTpcGrid->GetDeta();
 524               deltaPhi = fTpcGrid->GetDphi();
 525               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 526             }
 527             else {
 528               if(fDZ) {
 529                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 530                   if(nBin<fNumUnits+nElements+n0)
 531                     {
 532                       phi = eta = 0.;
 533                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
 534                       deltaEta = fGrid0->GetDeta();
 535                       deltaPhi = fGrid0->GetDphi();
 536                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 537                     }
 538                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 539                     {
 540                       phi = eta = 0.;
 541                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
 542                       deltaEta = fGrid1->GetDeta();
 543                       deltaPhi = fGrid1->GetDphi();
 544                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 545                     }
 546                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 547                     {
 548                       phi = eta = 0.;
 549                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
 550                       deltaEta = fGrid2->GetDeta();
 551                       deltaPhi = fGrid2->GetDphi();
 552                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 553                     }
 554                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 555                     {
 556                       phi = eta = 0.;
 557                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
 558                       deltaEta = fGrid3->GetDeta();
 559                       deltaPhi = fGrid3->GetDphi();
 560                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 561                     }
 562                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 563                     {
 564                       phi = eta = 0.;
 565                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
 566                       deltaEta = fGrid4->GetDeta();
 567                       deltaPhi = fGrid4->GetDphi();
 568                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 569                     }
 570                 }
 571               } // end if(fDZ)
 572             } // end else 2
 573           } // end else 1
 574         } // end loop on nBin
 575     } // end grid type == 1
 576   fArrayInitialised = 1;
 577 }