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