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 "AliAODMCParticle.h"
  40 #include "AliJetDummyGeo.h"
  41 #include "AliJetAODFillUnitArrayTracks.h"
  42 #include "AliJetAODFillUnitArrayEMCalDigits.h"
  43 #include "AliJetHadronCorrectionv1.h"
  44 #include "AliJetUnitArray.h"
  45
  46 ClassImp(AliJetAODReader)
  47
  48 AliJetAODReader::AliJetAODReader():
  49     AliJetReader(),
  50     fAOD(0x0),
  51     fRef(new TRefArray),
  52     fDebug(0),
  53     fOpt(0),
  54     fGeom(0),
  55     fHadCorr(0x0),
  56     fTpcGrid(0x0),
  57     fEmcalGrid(0x0),
  58     fGrid0(0),
  59     fGrid1(0),
  60     fGrid2(0),
  61     fGrid3(0),
  62     fGrid4(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            fChain->AddFile(Form("%s/%s/aod.root",dirName,name));
 123            a++;
 124        }
 125    }
 126
 127   gSystem->FreeDirectory(dir);
 128
 129   fAOD = 0;
 130   fChain->SetBranchAddress("AOD",&fAOD);
 131
 132   int nMax = fChain->GetEntries();
 133
 134   printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
 135
 136   // set number of events in header
 137   if (fReaderHeader->GetLastEvent() == -1)
 138     fReaderHeader->SetLastEvent(nMax);
 139   else {
 140     Int_t nUsr = fReaderHeader->GetLastEvent();
 141     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 142   }
 143 }
 144
 145 //____________________________________________________________________________
 146
 147 void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
 148     // Connect the tree
 149     // For AOD reader it's needed only to set the number of events
 150      fChain = (TChain*)      tree;
 151
 152      Int_t nMax = fChain->GetEntries();
 153      printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
 154      // set number of events in header
 155      if (fReaderHeader->GetLastEvent() == -1)
 156          fReaderHeader->SetLastEvent(nMax);
 157      else {
 158          Int_t nUsr = fReaderHeader->GetLastEvent();
 159          fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 160      }
 161 }
 162
 163
 164 Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){
 165
 166   //
 167   // filter for charge and for charged and neutral, no detector
 168   // response yet
 169   // physical priamries already selected
 170
 171   Int_t   pdg     = TMath::Abs(mcP->GetPdgCode());
 172
 173   // exclude neutrinos anyway
 174   if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
 175
 176   if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE;
 177   if(flag==AliJetAODReaderHeader::kChargedMC){
 178     if(mcP->Charge()==0)return kFALSE;
 179     return kTRUE;
 180   }
 181
 182   return kTRUE;
 183
 184 }
 185
 186
 187 Bool_t AliJetAODReader::FillMomentumArrayMC(){
 188
 189   //
 190   // This routine fetches the MC particles from the AOD
 191   // Depending on the flag all particles except neurinos are use
 192   // or only charged particles
 193   //
 194
 195   TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 196   if(!mcArray){
 197     Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
 198     return kFALSE;
 199   }
 200
 201   Int_t nMC = mcArray->GetEntriesFast();
 202
 203   // get number of tracks in event (for the loop)
 204   if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC);
 205
 206   // temporary storage of signal and pt cut flag
 207   Int_t* sflag  = new Int_t[nMC];
 208   Int_t* cflag  = new Int_t[nMC];
 209
 210   // get cuts set by user
 211   Float_t ptMin =  fReaderHeader->GetPtCut();
 212   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 213   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 214   Int_t mcTrack = 0;
 215   Float_t pt, eta;
 216   TVector3 p3;
 217
 218   Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
 219
 220
 221   for (Int_t it = 0; it < nMC; it++) {
 222     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
 223     if(!track->IsPhysicalPrimary())continue;
 224
 225     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 226     eta = p3.Eta();
 227     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 228     if(!AcceptAODMCParticle(track,readerFlag))continue;
 229     pt = p3.Pt();
 230
 231
 232
 233
 234     if (mcTrack == 0){
 235       fRef->Delete(); // make sure to delete before placement new...
 236       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 237     }
 238     new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag());
 239     sflag[mcTrack] = 1;
 240     cflag[mcTrack] = ( pt > ptMin ) ? 1: 0;
 241     fRef->Add(track);
 242     mcTrack++;
 243   }
 244   if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack);
 245   // set the signal flags
 246   fSignalFlag.Set(mcTrack,sflag);
 247   fCutFlag.Set(mcTrack,cflag);
 248
 249   delete [] sflag;
 250   delete [] cflag;
 251
 252   return kTRUE;
 253 }
 254
 255 //____________________________________________________________________________
 256
 257 Bool_t AliJetAODReader::FillMomentumArray()
 258 {
 259   // Clear momentum array
 260   ClearArray();
 261   fRef->Clear();
 262   fDebug = fReaderHeader->GetDebug();
 263
 264   if (!fAOD) {
 265       return kFALSE;
 266   }
 267
 268   if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
 269     return FillMomentumArrayMC();
 270   }
 271
 272   // get number of tracks in event (for the loop)
 273   Int_t nt = fAOD->GetNTracks();
 274   if(fDebug>0)printf("AOD tracks: %5d \t", nt);
 275
 276   // temporary storage of signal and pt cut flag
 277   Int_t* sflag  = new Int_t[nt];
 278   Int_t* cflag  = new Int_t[nt];
 279
 280   // get cuts set by user
 281   Float_t ptMin =  fReaderHeader->GetPtCut();
 282   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 283   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 284   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
 285
 286   //loop over tracks
 287   Int_t aodTrack = 0;
 288   Float_t pt, eta;
 289   TVector3 p3;
 290
 291   for (Int_t it = 0; it < nt; it++) {
 292     AliAODTrack *track = fAOD->GetTrack(it);
 293     UInt_t status = track->GetStatus();
 294
 295     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 296     p3.SetXYZ(mom[0],mom[1],mom[2]);
 297     pt = p3.Pt();
 298     eta = p3.Eta();
 299     if (status == 0) continue;
 300     if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
 301     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 302
 303     if (aodTrack == 0){
 304       fRef->Delete(); // make sure to delete before placement new...
 305       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 306     }
 307     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 308     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
 309     cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
 310     aodTrack++;
 311     fRef->Add(track);
 312   }
 313   if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack);
 314   // set the signal flags
 315   fSignalFlag.Set(aodTrack,sflag);
 316   fCutFlag.Set(aodTrack,cflag);
 317
 318   delete [] sflag;
 319   delete [] cflag;
 320
 321   return kTRUE;
 322 }
 323
 324 //__________________________________________________________
 325 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
 326 {
 327   //
 328   // Set flag to apply MIP correction fApplyMIPCorrection
 329   // - exclusive with fApplyFractionHadronicCorrection
 330   //
 331
 332   fApplyMIPCorrection = val;
 333   if(fApplyMIPCorrection == kTRUE)
 334     {
 335       SetApplyFractionHadronicCorrection(kFALSE);
 336       printf("Enabling MIP Correction \n");
 337     }
 338   else
 339     {
 340       printf("Disabling MIP Correction \n");
 341     }
 342 }
 343
 344 //__________________________________________________________
 345 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
 346 {
 347   //
 348   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
 349   // - exclusive with fApplyMIPCorrection
 350   //
 351
 352   fApplyFractionHadronicCorrection = val;
 353   if(fApplyFractionHadronicCorrection == kTRUE)
 354     {
 355       SetApplyMIPCorrection(kFALSE);
 356       printf("Enabling Fraction Hadronic Correction \n");
 357     }
 358   else
 359     {
 360       printf("Disabling Fraction Hadronic Correction \n");
 361     }
 362 }
 363
 364 //__________________________________________________________
 365 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
 366 {
 367   //
 368   // Set value to fFractionHadronicCorrection (default is 0.3)
 369   // apply EMC hadronic correction fApplyFractionHadronicCorrection
 370   // - exclusive with fApplyMIPCorrection
 371   //
 372
 373   fFractionHadronicCorrection = val;
 374   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
 375     {
 376       SetApplyFractionHadronicCorrection(kTRUE);
 377       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
 378     }
 379   else
 380     {
 381       SetApplyFractionHadronicCorrection(kFALSE);
 382     }
 383 }
 384
 385 //____________________________________________________________________________
 386 void AliJetAODReader::CreateTasks(TChain* tree)
 387 {
 388   //
 389   // For reader task initialization
 390   //
 391
 392   fDebug = fReaderHeader->GetDebug();
 393   fDZ = fReaderHeader->GetDZ();
 394   fTree = tree;
 395
 396   // Init EMCAL geometry and create UnitArray object
 397   SetEMCALGeometry();
 398   //  cout << "In create task" << endl;
 399   InitParameters();
 400   InitUnitArray();
 401
 402   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
 403   fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
 404   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
 405   fFillUAFromTracks->SetGeom(fGeom);
 406   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
 407   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 408
 409   if(fDZ)
 410     {
 411       fFillUAFromTracks->SetGrid0(fGrid0);
 412       fFillUAFromTracks->SetGrid1(fGrid1);
 413       fFillUAFromTracks->SetGrid2(fGrid2);
 414       fFillUAFromTracks->SetGrid3(fGrid3);
 415       fFillUAFromTracks->SetGrid4(fGrid4);
 416     }
 417   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
 418   fFillUAFromTracks->SetHadCorrector(fHadCorr);
 419   fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
 420   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
 421   fFillUAFromEMCalDigits->SetGeom(fGeom);
 422   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
 423   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
 424   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
 425   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
 426   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
 427   // Add the task to global task
 428   fFillUnitArray->Add(fFillUAFromTracks);
 429   fFillUnitArray->Add(fFillUAFromEMCalDigits);
 430   fFillUAFromTracks->SetActive(kFALSE);
 431   fFillUAFromEMCalDigits->SetActive(kFALSE);
 432
 433   cout << "Tasks instantiated at that stage ! " << endl;
 434   cout << "You can loop over events now ! " << endl;
 435
 436 }
 437
 438 //____________________________________________________________________________
 439 Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 440 {
 441   //
 442   // Main function
 443   // Fill the reader part
 444   //
 445
 446   fProcId = procid;
 447   fRefArray = refArray;
 448 //(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
 449 //(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
 450
 451   // clear momentum array
 452   ClearArray();
 453   fRef->Clear();
 454
 455   fDebug = fReaderHeader->GetDebug();
 456   fOpt = fReaderHeader->GetDetector();
 457
 458   if(!fAOD) {
 459     return kFALSE;
 460   }
 461
 462   // TPC only or Digits+TPC or Clusters+TPC
 463   if(fOpt%2==!0 && fOpt!=0){
 464     fFillUAFromTracks->SetAOD(fAOD);
 465     fFillUAFromTracks->SetActive(kTRUE);
 466     fFillUAFromTracks->SetUnitArray(fUnitArray);
 467     fFillUAFromTracks->SetRefArray(fRefArray);
 468     fFillUAFromTracks->SetReferences(fRef);
 469     fFillUAFromTracks->SetSignalFlag(fSignalFlag);
 470     fFillUAFromTracks->SetCutFlag(fCutFlag);
 471     fFillUAFromTracks->SetProcId(fProcId);
 472     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
 473     fFillUAFromTracks->Exec("tpc");
 474     if(fOpt==1){
 475       fNumCandidate = fFillUAFromTracks->GetMult();
 476       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 477     }
 478     fSignalFlag = fFillUAFromTracks->GetSignalFlag();
 479     fCutFlag = fFillUAFromTracks->GetCutFlag();
 480   }
 481
 482   // Digits only or Digits+TPC
 483   if(fOpt>=2 && fOpt<=3){
 484     fFillUAFromEMCalDigits->SetAOD(fAOD);
 485     fFillUAFromEMCalDigits->SetActive(kTRUE);
 486     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 487     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 488     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
 489     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 490     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 491     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
 492     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 493     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 494   }
 495
 496   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
 497
 498   return kTRUE;
 499 }
 500
 501 //____________________________________________________________________________
 502 Bool_t AliJetAODReader::SetEMCALGeometry()
 503 {
 504   //
 505   // Set the EMCal Geometry
 506   //
 507
 508   if (!fTree->GetFile())
 509     return kFALSE;
 510
 511   TString geomFile(fTree->GetFile()->GetName());
 512   geomFile.ReplaceAll("AliESDs", "geometry");
 513
 514   // temporary workaround for PROOF bug #18505
 515   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
 516   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
 517
 518   // Define EMCAL geometry to be able to read ESDs
 519   fGeom = AliJetDummyGeo::GetInstance();
 520   if (fGeom == 0)
 521     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
 522
 523   // To be setted to run some AliEMCALGeometry functions
 524   TGeoManager::Import(geomFile);
 525   fGeom->GetTransformationForSM();
 526   printf("\n EMCal Geometry set ! \n");
 527
 528   return kTRUE;
 529
 530 }
 531
 532 //____________________________________________________________________________
 533 void AliJetAODReader::InitParameters()
 534 {
 535   // Initialise parameters
 536   fOpt = fReaderHeader->GetDetector();
 537   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 538   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 539 }
 540
 541 //____________________________________________________________________________
 542 void AliJetAODReader::InitUnitArray()
 543 {
 544   //Initialises unit arrays
 545   Int_t nElements = fTpcGrid->GetNEntries();
 546   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
 547   if(fArrayInitialised) fUnitArray->Delete();
 548
 549   if(fTpcGrid->GetGridType()==0)
 550     { // Fill the following quantities :
 551       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
 552       // detector flag, in/out jet, pt cut, mass, cluster ID)
 553       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 554         {
 555           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
 556           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 557           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 558           deltaEta = fTpcGrid->GetDeta();
 559           deltaPhi = fTpcGrid->GetDphi();
 560           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 561         }
 562     }
 563
 564   if(fTpcGrid->GetGridType()==1)
 565     {
 566       Int_t nGaps = 0;
 567       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 568
 569       if(fDZ)
 570         {
 571           // Define a grid of cell for the gaps between SM
 572           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 573           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 574           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 575           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 576           fGrid0->SetGridType(0);
 577           fGrid0->SetMatrixIndexes();
 578           fGrid0->SetIndexIJ();
 579           n0 = fGrid0->GetNEntries();
 580           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 581           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 582           fGrid1->SetGridType(0);
 583           fGrid1->SetMatrixIndexes();
 584           fGrid1->SetIndexIJ();
 585           n1 = fGrid1->GetNEntries();
 586           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 587           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 588           fGrid2->SetGridType(0);
 589           fGrid2->SetMatrixIndexes();
 590           fGrid2->SetIndexIJ();
 591           n2 = fGrid2->GetNEntries();
 592           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 593           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 594           fGrid3->SetGridType(0);
 595           fGrid3->SetMatrixIndexes();
 596           fGrid3->SetIndexIJ();
 597           n3 = fGrid3->GetNEntries();
 598           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 599           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 600           fGrid4->SetGridType(0);
 601           fGrid4->SetMatrixIndexes();
 602           fGrid4->SetIndexIJ();
 603           n4 = fGrid4->GetNEntries();
 604
 605           nGaps = n0+n1+n2+n3+n4;
 606
 607         }
 608
 609       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 610         {
 611           if(nBin<fNumUnits)
 612             {
 613               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 614               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 615               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 616               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 617               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 618               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 619             }
 620           else {
 621             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 622               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 623               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 624               deltaEta = fTpcGrid->GetDeta();
 625               deltaPhi = fTpcGrid->GetDphi();
 626               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 627             }
 628             else {
 629               if(fDZ) {
 630                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 631                   if(nBin<fNumUnits+nElements+n0)
 632                     {
 633                       phi = eta = 0.;
 634                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
 635                       deltaEta = fGrid0->GetDeta();
 636                       deltaPhi = fGrid0->GetDphi();
 637                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 638                     }
 639                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 640                     {
 641                       phi = eta = 0.;
 642                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
 643                       deltaEta = fGrid1->GetDeta();
 644                       deltaPhi = fGrid1->GetDphi();
 645                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 646                     }
 647                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 648                     {
 649                       phi = eta = 0.;
 650                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
 651                       deltaEta = fGrid2->GetDeta();
 652                       deltaPhi = fGrid2->GetDphi();
 653                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 654                     }
 655                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 656                     {
 657                       phi = eta = 0.;
 658                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
 659                       deltaEta = fGrid3->GetDeta();
 660                       deltaPhi = fGrid3->GetDphi();
 661                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 662                     }
 663                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 664                     {
 665                       phi = eta = 0.;
 666                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
 667                       deltaEta = fGrid4->GetDeta();
 668                       deltaPhi = fGrid4->GetDphi();
 669                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 670                     }
 671                 }
 672               } // end if(fDZ)
 673             } // end else 2
 674           } // end else 1
 675         } // end loop on nBin
 676     } // end grid type == 1
 677   fArrayInitialised = 1;
 678 }
 679