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