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