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   // get cuts set by user
 269   Float_t ptMin =  fReaderHeader->GetPtCut();
 270   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 271   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 272   UInt_t  filterMask =  ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
 273
 274   // ----- number of tracks -----
 275   Int_t nTracksStd    = 0;
 276   Short_t mcReaderFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
 277   TClonesArray *mcArray = 0x0;
 278   // check if we have to read from MC branch
 279   if (((AliJetAODReaderHeader*)fReaderHeader)->GetReadStdBranch()) {
 280     if(mcReaderFlag) {
 281       mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 282       if(!mcArray){
 283         Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
 284       }
 285       nTracksStd = mcArray->GetEntriesFast();
 286     }
 287     else {
 288       nTracksStd = fAOD->GetNTracks();
 289       //      printf("no. of standard tracks: %i\n", nTracksStd);
 290     }
 291   }
 292   Int_t nTracksNonStd = 0;
 293   TClonesArray *nonStdTracks = 0x0;
 294   if (((AliJetAODReaderHeader*)fReaderHeader)->GetReadNonStdBranch()) {
 295     nonStdTracks =
 296       (TClonesArray*) fAOD->FindListObject(((AliJetAODReaderHeader*)fReaderHeader)->GetNonStdBranch());
 297     if (nonStdTracks)
 298       nTracksNonStd = nonStdTracks->GetEntries();
 299     //    printf("no. of non-standard tracks: %i\n", nTracksNonStd);
 300   }
 301   Int_t nTracks = nTracksStd + nTracksNonStd;
 302
 303   // temporary storage of signal and pt cut flag
 304   Int_t* sflag  = new Int_t[nTracks];
 305   Int_t* cflag  = new Int_t[nTracks];
 306
 307   // loop over tracks
 308   Int_t aodTrack = 0;
 309   Float_t pt, eta;
 310   TVector3 p3;
 311
 312   // ----- looping over standard branch -----
 313   if (mcArray) {
 314     for (Int_t iTrack = 0; iTrack < nTracksStd; iTrack++) {
 315       AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(iTrack);
 316       if(!track->IsPhysicalPrimary())continue;
 317
 318       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 319       eta = p3.Eta();
 320       if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 321       if(!AcceptAODMCParticle(track,mcReaderFlag))continue;
 322       pt = p3.Pt();
 323
 324       if (aodTrack == 0){
 325         fRef->Delete(); // make sure to delete before placement new...
 326         new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 327       }
 328       new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 329       sflag[aodTrack] = 1;
 330       cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
 331       fRef->Add(track);
 332       aodTrack++;
 333     }
 334   }
 335   else {
 336     for (Int_t iTrack = 0; iTrack < nTracksStd; iTrack++) {
 337       AliAODTrack *track = fAOD->GetTrack(iTrack);
 338       UInt_t status = track->GetStatus();
 339
 340       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 341       p3.SetXYZ(mom[0],mom[1],mom[2]);
 342       pt = p3.Pt();
 343       eta = p3.Eta();
 344       if (status == 0) continue;
 345       if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
 346       if ( (eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 347
 348       if (aodTrack == 0){
 349         fRef->Delete(); // make sure to delete before placement new...
 350         new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 351       }
 352       new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 353       sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
 354       cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
 355       aodTrack++;
 356       fRef->Add(track);
 357     }
 358   }
 359
 360   // ----- reading of non-standard branch -----
 361   for (Int_t iTrack = 0; iTrack < nTracksNonStd; iTrack++) {
 362     AliVParticle *track = dynamic_cast<AliVParticle*> ((*nonStdTracks)[iTrack]);
 363     if (!track)
 364       continue;
 365
 366     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 367     p3.SetXYZ(mom[0],mom[1],mom[2]);
 368     pt = p3.Pt();
 369     eta = p3.Eta();
 370
 371     // which cuts to apply if not AOD track (e.g. MC) ???
 372     AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
 373     if (trackAOD) {
 374       if (trackAOD->GetStatus() == 0)
 375         continue;
 376       if ((filterMask > 0) && !(trackAOD->TestFilterBit(filterMask)))
 377         continue;
 378     }
 379     if ((eta > etaMax) || (eta < etaMin)) continue;      // checking eta cut
 380
 381     if (aodTrack == 0){
 382       fRef->Delete(); // make sure to delete before placement new...
 383       new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
 384     }
 385     new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
 386     sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
 387     cflag[aodTrack] = ( pt > ptMin ) ? 1 : 0;
 388     aodTrack++;
 389     fRef->Add(track);
 390     //    printf("added non-standard track\n");
 391   }
 392
 393   // set the signal flags
 394   fSignalFlag.Set(aodTrack,sflag);
 395   fCutFlag.Set(aodTrack,cflag);
 396
 397   delete [] sflag;
 398   delete [] cflag;
 399
 400   return kTRUE;
 401 }
 402
 403 //__________________________________________________________
 404 void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
 405 {
 406   //
 407   // Set flag to apply MIP correction fApplyMIPCorrection
 408   // - exclusive with fApplyFractionHadronicCorrection
 409   //
 410
 411   fApplyMIPCorrection = val;
 412   if(fApplyMIPCorrection == kTRUE)
 413     {
 414       SetApplyFractionHadronicCorrection(kFALSE);
 415       printf("Enabling MIP Correction \n");
 416     }
 417   else
 418     {
 419       printf("Disabling MIP Correction \n");
 420     }
 421 }
 422
 423 //__________________________________________________________
 424 void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
 425 {
 426   //
 427   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
 428   // - exclusive with fApplyMIPCorrection
 429   //
 430
 431   fApplyFractionHadronicCorrection = val;
 432   if(fApplyFractionHadronicCorrection == kTRUE)
 433     {
 434       SetApplyMIPCorrection(kFALSE);
 435       printf("Enabling Fraction Hadronic Correction \n");
 436     }
 437   else
 438     {
 439       printf("Disabling Fraction Hadronic Correction \n");
 440     }
 441 }
 442
 443 //__________________________________________________________
 444 void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
 445 {
 446   //
 447   // Set value to fFractionHadronicCorrection (default is 0.3)
 448   // apply EMC hadronic correction fApplyFractionHadronicCorrection
 449   // - exclusive with fApplyMIPCorrection
 450   //
 451
 452   fFractionHadronicCorrection = val;
 453   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
 454     {
 455       SetApplyFractionHadronicCorrection(kTRUE);
 456       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
 457     }
 458   else
 459     {
 460       SetApplyFractionHadronicCorrection(kFALSE);
 461     }
 462 }
 463
 464 //____________________________________________________________________________
 465 void AliJetAODReader::CreateTasks(TChain* tree)
 466 {
 467   //
 468   // For reader task initialization
 469   //
 470
 471   fDebug = fReaderHeader->GetDebug();
 472   fDZ = fReaderHeader->GetDZ();
 473   fTree = tree;
 474
 475   // Init EMCAL geometry and create UnitArray object
 476   SetEMCALGeometry();
 477   //  cout << "In create task" << endl;
 478   InitParameters();
 479   InitUnitArray();
 480
 481   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
 482   fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
 483   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
 484   fFillUAFromTracks->SetGeom(fGeom);
 485   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
 486   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 487
 488   if(fDZ)
 489     {
 490       fFillUAFromTracks->SetGrid0(fGrid0);
 491       fFillUAFromTracks->SetGrid1(fGrid1);
 492       fFillUAFromTracks->SetGrid2(fGrid2);
 493       fFillUAFromTracks->SetGrid3(fGrid3);
 494       fFillUAFromTracks->SetGrid4(fGrid4);
 495     }
 496   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
 497   fFillUAFromTracks->SetHadCorrector(fHadCorr);
 498   fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
 499   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
 500   fFillUAFromEMCalDigits->SetGeom(fGeom);
 501   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
 502   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
 503   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
 504   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
 505   fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
 506   // Add the task to global task
 507   fFillUnitArray->Add(fFillUAFromTracks);
 508   fFillUnitArray->Add(fFillUAFromEMCalDigits);
 509   fFillUAFromTracks->SetActive(kFALSE);
 510   fFillUAFromEMCalDigits->SetActive(kFALSE);
 511
 512   cout << "Tasks instantiated at that stage ! " << endl;
 513   cout << "You can loop over events now ! " << endl;
 514
 515 }
 516
 517 //____________________________________________________________________________
 518 Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 519 {
 520   //
 521   // Main function
 522   // Fill the reader part
 523   //
 524
 525   fProcId = procid;
 526   fRefArray = refArray;
 527 //(not used ?)  Int_t nEntRef = fRefArray->GetEntries();
 528 //(not used ?)  Int_t nEntUnit = fUnitArray->GetEntries();
 529
 530   // clear momentum array
 531   ClearArray();
 532   fRef->Clear();
 533
 534   fDebug = fReaderHeader->GetDebug();
 535   fOpt = fReaderHeader->GetDetector();
 536
 537   if(!fAOD) {
 538     return kFALSE;
 539   }
 540
 541   // TPC only or Digits+TPC or Clusters+TPC
 542   if(fOpt%2==!0 && fOpt!=0){
 543     fFillUAFromTracks->SetAOD(fAOD);
 544     fFillUAFromTracks->SetActive(kTRUE);
 545     fFillUAFromTracks->SetUnitArray(fUnitArray);
 546     fFillUAFromTracks->SetRefArray(fRefArray);
 547     fFillUAFromTracks->SetReferences(fRef);
 548     fFillUAFromTracks->SetSignalFlag(fSignalFlag);
 549     fFillUAFromTracks->SetCutFlag(fCutFlag);
 550     fFillUAFromTracks->SetProcId(fProcId);
 551     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
 552     fFillUAFromTracks->Exec("tpc");
 553     if(fOpt==1){
 554       fNumCandidate = fFillUAFromTracks->GetMult();
 555       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 556     }
 557     fSignalFlag = fFillUAFromTracks->GetSignalFlag();
 558     fCutFlag = fFillUAFromTracks->GetCutFlag();
 559   }
 560
 561   // Digits only or Digits+TPC
 562   if(fOpt>=2 && fOpt<=3){
 563     fFillUAFromEMCalDigits->SetAOD(fAOD);
 564     fFillUAFromEMCalDigits->SetActive(kTRUE);
 565     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 566     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 567     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
 568     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 569     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 570     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
 571     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 572     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 573   }
 574
 575   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
 576
 577   return kTRUE;
 578 }
 579
 580 //____________________________________________________________________________
 581 Bool_t AliJetAODReader::SetEMCALGeometry()
 582 {
 583   //
 584   // Set the EMCal Geometry
 585   //
 586
 587   if (!fTree->GetFile())
 588     return kFALSE;
 589
 590   TString geomFile(fTree->GetFile()->GetName());
 591   geomFile.ReplaceAll("AliESDs", "geometry");
 592
 593   // temporary workaround for PROOF bug #18505
 594   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
 595   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
 596
 597   // Define EMCAL geometry to be able to read ESDs
 598   fGeom = AliJetDummyGeo::GetInstance();
 599   if (fGeom == 0)
 600     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
 601
 602   // To be setted to run some AliEMCALGeometry functions
 603   TGeoManager::Import(geomFile);
 604   fGeom->GetTransformationForSM();
 605   printf("\n EMCal Geometry set ! \n");
 606
 607   return kTRUE;
 608
 609 }
 610
 611 //____________________________________________________________________________
 612 void AliJetAODReader::InitParameters()
 613 {
 614   // Initialise parameters
 615   fOpt = fReaderHeader->GetDetector();
 616   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 617   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 618 }
 619
 620 //____________________________________________________________________________
 621 void AliJetAODReader::InitUnitArray()
 622 {
 623   //Initialises unit arrays
 624   Int_t nElements = fTpcGrid->GetNEntries();
 625   Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
 626   if(fArrayInitialised) fUnitArray->Delete();
 627
 628   if(fTpcGrid->GetGridType()==0)
 629     { // Fill the following quantities :
 630       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
 631       // detector flag, in/out jet, pt cut, mass, cluster ID)
 632       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 633         {
 634           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
 635           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 636           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 637           deltaEta = fTpcGrid->GetDeta();
 638           deltaPhi = fTpcGrid->GetDphi();
 639           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 640         }
 641     }
 642
 643   if(fTpcGrid->GetGridType()==1)
 644     {
 645       Int_t nGaps = 0;
 646       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 647
 648       if(fDZ)
 649         {
 650           // Define a grid of cell for the gaps between SM
 651           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 652           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 653           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 654           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 655           fGrid0->SetGridType(0);
 656           fGrid0->SetMatrixIndexes();
 657           fGrid0->SetIndexIJ();
 658           n0 = fGrid0->GetNEntries();
 659           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 660           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 661           fGrid1->SetGridType(0);
 662           fGrid1->SetMatrixIndexes();
 663           fGrid1->SetIndexIJ();
 664           n1 = fGrid1->GetNEntries();
 665           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 666           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 667           fGrid2->SetGridType(0);
 668           fGrid2->SetMatrixIndexes();
 669           fGrid2->SetIndexIJ();
 670           n2 = fGrid2->GetNEntries();
 671           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 672           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 673           fGrid3->SetGridType(0);
 674           fGrid3->SetMatrixIndexes();
 675           fGrid3->SetIndexIJ();
 676           n3 = fGrid3->GetNEntries();
 677           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 678           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 679           fGrid4->SetGridType(0);
 680           fGrid4->SetMatrixIndexes();
 681           fGrid4->SetIndexIJ();
 682           n4 = fGrid4->GetNEntries();
 683
 684           nGaps = n0+n1+n2+n3+n4;
 685
 686         }
 687
 688       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 689         {
 690           if(nBin<fNumUnits)
 691             {
 692               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 693               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 694               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 695               deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 696               deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 697               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 698             }
 699           else {
 700             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 701               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 702               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 703               deltaEta = fTpcGrid->GetDeta();
 704               deltaPhi = fTpcGrid->GetDphi();
 705               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 706             }
 707             else {
 708               if(fDZ) {
 709                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 710                   if(nBin<fNumUnits+nElements+n0)
 711                     {
 712                       phi = eta = 0.;
 713                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
 714                       deltaEta = fGrid0->GetDeta();
 715                       deltaPhi = fGrid0->GetDphi();
 716                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 717                     }
 718                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 719                     {
 720                       phi = eta = 0.;
 721                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
 722                       deltaEta = fGrid1->GetDeta();
 723                       deltaPhi = fGrid1->GetDphi();
 724                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 725                     }
 726                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 727                     {
 728                       phi = eta = 0.;
 729                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
 730                       deltaEta = fGrid2->GetDeta();
 731                       deltaPhi = fGrid2->GetDphi();
 732                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 733                     }
 734                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 735                     {
 736                       phi = eta = 0.;
 737                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
 738                       deltaEta = fGrid3->GetDeta();
 739                       deltaPhi = fGrid3->GetDphi();
 740                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 741                     }
 742                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 743                     {
 744                       phi = eta = 0.;
 745                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
 746                       deltaEta = fGrid4->GetDeta();
 747                       deltaPhi = fGrid4->GetDphi();
 748                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 749                     }
 750                 }
 751               } // end if(fDZ)
 752             } // end else 2
 753           } // end else 1
 754         } // end loop on nBin
 755     } // end grid type == 1
 756   fArrayInitialised = 1;
 757 }
 758