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 "AliJetHadronCorrection.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     fPtCut(0),
  64     fApplyElectronCorrection(kFALSE),
  65     fApplyMIPCorrection(kTRUE),
  66     fApplyFractionHadronicCorrection(kFALSE),
  67     fFractionHadronicCorrection(0.3),
  68     fNumUnits(0),
  69     fMass(0),
  70     fSign(0),
  71     fNIn(0),
  72     fDZ(0),
  73     fNeta(0),
  74     fNphi(0),
  75     fRefArray(0x0),
  76     fProcId(kFALSE)
  77 {
  78   // Constructor
  79 }
  80
  81 //____________________________________________________________________________
  82
  83 AliJetAODReader::~AliJetAODReader()
  84 {
  85   // Destructor
  86     delete fAOD;
  87     delete fRef;
  88     delete fTpcGrid;
  89     delete fEmcalGrid;
  90     if(fDZ)
  91       {
  92         delete fGrid0;
  93         delete fGrid1;
  94         delete fGrid2;
  95         delete fGrid3;
  96         delete fGrid4;
  97       }
  98
  99 }
 100
 101 //____________________________________________________________________________
 102
 103 void AliJetAODReader::OpenInputFiles()
 104 {
 105     // Open the necessary input files
 106     // chain for the AODs
 107   fChain   = new TChain("aodTree");
 108
 109   // get directory and pattern name from the header
 110    const char* dirName=fReaderHeader->GetDirectory();
 111    const char* pattern=fReaderHeader->GetPattern();
 112
 113 //   // Add files matching patters to the chain
 114
 115    void *dir  = gSystem->OpenDirectory(dirName);
 116    const char *name = 0x0;
 117    int naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
 118    int a = 0;
 119    while ((name = gSystem->GetDirEntry(dir))){
 120        if (a>=naod) continue;
 121
 122        if (strstr(name,pattern)){
 123            char path[256];
 124            sprintf(path,"%s/%s/aod.root",dirName,name);
 125            fChain->AddFile(path);
 126            a++;
 127        }
 128    }
 129
 130   gSystem->FreeDirectory(dir);
 131
 132   fAOD = 0;
 133   fChain->SetBranchAddress("AOD",&fAOD);
 134
 135   int nMax = fChain->GetEntries();
 136
 137   printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
 138
 139   // set number of events in header
 140   if (fReaderHeader->GetLastEvent() == -1)
 141     fReaderHeader->SetLastEvent(nMax);
 142   else {
 143     Int_t nUsr = fReaderHeader->GetLastEvent();
 144     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 145   }
 146 }
 147
 148 //____________________________________________________________________________
 149
 150 void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
 151     // Connect the tree
 152     // For AOD reader it's needed only to set the number of events
 153      fChain = (TChain*)      tree;
 154
 155      Int_t nMax = fChain->GetEntries();
 156      printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
 157      // set number of events in header
 158      if (fReaderHeader->GetLastEvent() == -1)
 159          fReaderHeader->SetLastEvent(nMax);
 160      else {
 161          Int_t nUsr = fReaderHeader->GetLastEvent();
 162          fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 163      }
 164 }
 165
 166
 167 Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){
 168
 169   //
 170   // filter for charge and for charged and neutral, no detector
 171   // response yet
 172   // physical priamries already selected
 173
 174   Int_t   pdg     = TMath::Abs(mcP->GetPdgCode());
 175
 176   // exclude neutrinos anyway
 177   if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
 178
 179   if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE;
 180   if(flag==AliJetAODReaderHeader::kChargedMC){
 181     if(mcP->Charge()==0)return kFALSE;
 182     return kTRUE;
 183   }
 184
 185   return kTRUE;
 186
 187 }
 188
 189
 190 Bool_t AliJetAODReader::FillMomentumArrayMC(){
 191
 192   //
 193   // This routine fetches the MC particles from the AOD
 194   // Depending on the flag all particles except neurinos are use
 195   // or only charged particles
 196   //
 197
 198   TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 199   if(!mcArray){
 200     Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
 201     return kFALSE;
 202   }
 203
 204   Int_t nMC = mcArray->GetEntriesFast();
 205
 206   // get number of tracks in event (for the loop)
 207   if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC);
 208
 209   // temporary storage of signal and pt cut flag
 210   Int_t* sflag  = new Int_t[nMC];
 211   Int_t* cflag  = new Int_t[nMC];
 212
 213   // get cuts set by user
 214   Float_t ptMin =  fReaderHeader->GetPtCut();
 215   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 216   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 217   Int_t mcTrack = 0;
 218   Float_t pt, eta;
 219   TVector3 p3;
 220
 221   Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
 222
 223
 224   for (Int_t it = 0; it < nMC; it++) {
 225     AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
 226     if(!track->IsPhysicalPrimary())continue;
 227
 228     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 229     eta = p3.Eta();
 230     if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 231     if(!AcceptAODMCParticle(track,readerFlag))continue;
 232     pt = p3.Pt();
 233
 234
 235
 236
 237     if (mcTrack == 0) new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 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(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
 265     return FillMomentumArrayMC();
 266   }
 267
 268
 269
 270   if (!fAOD) {
 271       return kFALSE;
 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) 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
 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
 454   fDebug = fReaderHeader->GetDebug();
 455   fOpt = fReaderHeader->GetDetector();
 456
 457   if(!fAOD) {
 458     return kFALSE;
 459   }
 460
 461   // TPC only or Digits+TPC or Clusters+TPC
 462   if(fOpt%2==!0 && fOpt!=0){
 463     fFillUAFromTracks->SetAOD(fAOD);
 464     fFillUAFromTracks->SetActive(kTRUE);
 465     fFillUAFromTracks->SetUnitArray(fUnitArray);
 466     fFillUAFromTracks->SetRefArray(fRefArray);
 467     fFillUAFromTracks->SetProcId(fProcId);
 468     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
 469     fFillUAFromTracks->Exec("tpc");
 470     if(fOpt==1){
 471       fNumCandidate = fFillUAFromTracks->GetMult();
 472       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 473     }
 474   }
 475
 476   // Digits only or Digits+TPC
 477   if(fOpt>=2 && fOpt<=3){
 478     fFillUAFromEMCalDigits->SetAOD(fAOD);
 479     fFillUAFromEMCalDigits->SetActive(kTRUE);
 480     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 481     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 482     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
 483     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 484     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 485     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
 486     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 487     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 488   }
 489
 490   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
 491
 492   return kTRUE;
 493 }
 494
 495 //____________________________________________________________________________
 496 Bool_t AliJetAODReader::SetEMCALGeometry()
 497 {
 498   //
 499   // Set the EMCal Geometry
 500   //
 501
 502   if (!fTree->GetFile())
 503     return kFALSE;
 504
 505   TString geomFile(fTree->GetFile()->GetName());
 506   geomFile.ReplaceAll("AliESDs", "geometry");
 507
 508   // temporary workaround for PROOF bug #18505
 509   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
 510   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
 511
 512   // Define EMCAL geometry to be able to read ESDs
 513   fGeom = AliJetDummyGeo::GetInstance();
 514   if (fGeom == 0)
 515     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
 516
 517   // To be setted to run some AliEMCALGeometry functions
 518   TGeoManager::Import(geomFile);
 519   fGeom->GetTransformationForSM();
 520   printf("\n EMCal Geometry set ! \n");
 521
 522   return kTRUE;
 523
 524 }
 525
 526 //____________________________________________________________________________
 527 void AliJetAODReader::InitParameters()
 528 {
 529   // Initialise parameters
 530   fOpt = fReaderHeader->GetDetector();
 531   //  fHCorrection    = 0;          // For hadron correction
 532   fHadCorr        = 0;          // For hadron correction
 533   if(fEFlag==kFALSE){
 534     if(fOpt==0 || fOpt==1)
 535       fECorrection    = 0;        // For electron correction
 536     else fECorrection = 1;        // For electron correction
 537   }
 538   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 539   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 540 }
 541
 542 //____________________________________________________________________________
 543 void AliJetAODReader::InitUnitArray()
 544 {
 545   //Initialises unit arrays
 546   Int_t nElements = fTpcGrid->GetNEntries();
 547   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
 548   if(fArrayInitialised) fUnitArray->Delete();
 549
 550   if(fTpcGrid->GetGridType()==0)
 551     { // Fill the following quantities :
 552       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
 553       // detector flag, in/out jet, pt cut, mass, cluster ID)
 554       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 555         {
 556           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
 557           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 558           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 559           deltaEta = fTpcGrid->GetDeta();
 560           deltaPhi = fTpcGrid->GetDphi();
 561           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 562         }
 563     }
 564
 565   if(fTpcGrid->GetGridType()==1)
 566     {
 567       Int_t nGaps = 0;
 568       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 569
 570       if(fDZ)
 571         {
 572           // Define a grid of cell for the gaps between SM
 573           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 574           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 575           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 576           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 577           fGrid0->SetGridType(0);
 578           fGrid0->SetMatrixIndexes();
 579           fGrid0->SetIndexIJ();
 580           n0 = fGrid0->GetNEntries();
 581           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 582           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 583           fGrid1->SetGridType(0);
 584           fGrid1->SetMatrixIndexes();
 585           fGrid1->SetIndexIJ();
 586           n1 = fGrid1->GetNEntries();
 587           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 588           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 589           fGrid2->SetGridType(0);
 590           fGrid2->SetMatrixIndexes();
 591           fGrid2->SetIndexIJ();
 592           n2 = fGrid2->GetNEntries();
 593           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 594           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 595           fGrid3->SetGridType(0);
 596           fGrid3->SetMatrixIndexes();
 597           fGrid3->SetIndexIJ();
 598           n3 = fGrid3->GetNEntries();
 599           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 600           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 601           fGrid4->SetGridType(0);
 602           fGrid4->SetMatrixIndexes();
 603           fGrid4->SetIndexIJ();
 604           n4 = fGrid4->GetNEntries();
 605
 606           nGaps = n0+n1+n2+n3+n4;
 607
 608         }
 609
 610       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 611         {
 612           if(nBin<fNumUnits)
 613             {
 614               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 615               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 616               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 617               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 618               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 619               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 620             }
 621           else {
 622             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 623               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 624               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 625               deltaEta = fTpcGrid->GetDeta();
 626               deltaPhi = fTpcGrid->GetDphi();
 627               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 628             }
 629             else {
 630               if(fDZ) {
 631                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 632                   if(nBin<fNumUnits+nElements+n0)
 633                     {
 634                       phi = eta = 0.;
 635                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
 636                       deltaEta = fGrid0->GetDeta();
 637                       deltaPhi = fGrid0->GetDphi();
 638                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 639                     }
 640                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 641                     {
 642                       phi = eta = 0.;
 643                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
 644                       deltaEta = fGrid1->GetDeta();
 645                       deltaPhi = fGrid1->GetDphi();
 646                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 647                     }
 648                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 649                     {
 650                       phi = eta = 0.;
 651                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
 652                       deltaEta = fGrid2->GetDeta();
 653                       deltaPhi = fGrid2->GetDphi();
 654                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 655                     }
 656                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 657                     {
 658                       phi = eta = 0.;
 659                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
 660                       deltaEta = fGrid3->GetDeta();
 661                       deltaPhi = fGrid3->GetDphi();
 662                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 663                     }
 664                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 665                     {
 666                       phi = eta = 0.;
 667                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
 668                       deltaEta = fGrid4->GetDeta();
 669                       deltaPhi = fGrid4->GetDphi();
 670                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 671                     }
 672                 }
 673               } // end if(fDZ)
 674             } // end else 2
 675           } // end else 1
 676         } // end loop on nBin
 677     } // end grid type == 1
 678   fArrayInitialised = 1;
 679 }
 680