JETAN/AliJetESDReader.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 ESD Reader
  18 // ESD reader for jet analysis
  19 // Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
  20 //          Magali Estienne <magali.estienne@IReS.in2p3.fr>
  21 //-------------------------------------------------------------------------
  22
  23 // --- Standard library ---
  24 #include <Riostream.h>
  25
  26 // --- ROOT system ---
  27 #include <TSystem.h>
  28 #include <TStopwatch.h>
  29 #include <TLorentzVector.h>
  30 #include <TVector3.h>
  31 #include <TTask.h>
  32 #include <TGeoManager.h>
  33 #include <assert.h>
  34 #include <TRefArray.h>
  35 #include <TMath.h>
  36 #include <TChain.h>
  37
  38
  39 // --- AliRoot header files ---
  40 #include "AliJetESDReader.h"
  41 #include "AliJetESDReaderHeader.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESD.h"
  44 #include "AliESDtrack.h"
  45 #include "AliJetDummyGeo.h"
  46 #include "AliJetFillUnitArrayTracks.h"
  47 #include "AliJetFillUnitArrayEMCalDigits.h"
  48 #include "AliJetUnitArray.h"
  49
  50 ClassImp(AliJetESDReader)
  51
  52 AliJetESDReader::AliJetESDReader():
  53   AliJetReader(),
  54   fGeom(0),
  55   fChain(0x0),
  56   fESD(0x0),
  57   fHadCorr(0x0),
  58   fTpcGrid(0x0),
  59   fEmcalGrid(0x0),
  60   fGrid0(0),
  61   fGrid1(0),
  62   fGrid2(0),
  63   fGrid3(0),
  64   fGrid4(0),
  65   fPtCut(0),
  66   fHCorrection(0),
  67   fNumUnits(0),
  68   fDebug(0),
  69   fMass(0),
  70   fSign(0),
  71   fNIn(0),
  72   fOpt(0),
  73   fDZ(0),
  74   fNeta(0),
  75   fNphi(0),
  76   fArrayInitialised(0)
  77 {
  78   // Constructor
  79 }
  80
  81 //____________________________________________________________________________
  82 AliJetESDReader::~AliJetESDReader()
  83 {
  84   // Destructor
  85     delete fChain;
  86     delete fESD;
  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 void AliJetESDReader::OpenInputFiles()
 101 {
 102   // chain for the ESDs
 103   fChain   = new TChain("esdTree");
 104
 105   // get directory and pattern name from the header
 106    const char* dirName=fReaderHeader->GetDirectory();
 107    const char* pattern=fReaderHeader->GetPattern();
 108
 109    // Add files matching patters to the chain
 110    void *dir  = gSystem->OpenDirectory(dirName);
 111    const char *name = 0x0;
 112    int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
 113    int a = 0;
 114    while ((name = gSystem->GetDirEntry(dir))){
 115        if (a>=nesd) continue;
 116
 117        if (strstr(name,pattern)){
 118          printf("Adding %s\n",name);
 119          char path[256];
 120          //        sprintf(path,"%s/%s/AliESDs.root",dirName,name);
 121          sprintf(path,"%s/%s/AliESDs.root",dirName,name);
 122          fChain->AddFile(path);
 123          a++;
 124        }
 125    }
 126
 127   gSystem->FreeDirectory(dir);
 128
 129   fESD = new AliESDEvent();
 130   fESD->ReadFromTree(fChain);
 131
 132   int nMax = fChain->GetEntries();
 133
 134   printf("\n AliJetESDReader: 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 void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
 147     // Connect the tree
 148      fESD   = (AliESDEvent*) esd;
 149 }
 150
 151 //____________________________________________________________________________
 152 Bool_t AliJetESDReader::FillMomentumArray(Int_t /*event*/)
 153 {
 154   // Fill momentum array
 155
 156   Int_t goodTrack = 0;
 157   Int_t nt = 0;
 158   Float_t pt, eta;
 159   TVector3 p3;
 160
 161   // clear momentum array
 162   ClearArray();
 163   fDebug = fReaderHeader->GetDebug();
 164   fOpt = fReaderHeader->GetDetector();
 165
 166   if (!fESD) {
 167       return kFALSE;
 168   }
 169
 170   // get number of tracks in event (for the loop)
 171   nt = fESD->GetNumberOfTracks();
 172   printf("Fill Momentum Array %5d  \n", nt);
 173
 174   // temporary storage of signal and pt cut flag
 175   Int_t* sflag  = new Int_t[nt];
 176   Int_t* cflag  = new Int_t[nt];
 177
 178   // get cuts set by user
 179   Float_t ptMin =  fReaderHeader->GetPtCut();
 180   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 181   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 182
 183   //loop over tracks
 184   for (Int_t it = 0; it < nt; it++) {
 185       AliESDtrack *track = fESD->GetTrack(it);
 186       UInt_t status = track->GetStatus();
 187
 188       Double_t mom[3];
 189       track->GetPxPyPz(mom);
 190       p3.SetXYZ(mom[0],mom[1],mom[2]);
 191       pt = p3.Pt();
 192       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
 193       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
 194       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 195           && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
 196       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 197           && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
 198       eta = p3.Eta();
 199       if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
 200
 201       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
 202       sflag[goodTrack]=0;
 203       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 204       cflag[goodTrack]=0;
 205       if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
 206       goodTrack++;
 207   }
 208   // set the signal flags
 209   fSignalFlag.Set(goodTrack,sflag);
 210   fCutFlag.Set(goodTrack,cflag);
 211   return kTRUE;
 212
 213 }
 214
 215 //____________________________________________________________________________
 216 void AliJetESDReader::CreateTasks()
 217 {
 218   fDebug = fReaderHeader->GetDebug();
 219   fDZ = fReaderHeader->GetDZ();
 220
 221   // Init EMCAL geometry and create UnitArray object
 222   SetEMCALGeometry();
 223   InitParameters();
 224   InitUnitArray();
 225
 226   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
 227   fFillUAFromTracks = new AliJetFillUnitArrayTracks();
 228   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
 229   fFillUAFromTracks->SetGeom(fGeom);
 230   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
 231   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 232
 233   if(fDZ)
 234     {
 235       fFillUAFromTracks->SetGrid0(fGrid0);
 236       fFillUAFromTracks->SetGrid1(fGrid1);
 237       fFillUAFromTracks->SetGrid2(fGrid2);
 238       fFillUAFromTracks->SetGrid3(fGrid3);
 239       fFillUAFromTracks->SetGrid4(fGrid4);
 240     }
 241   fFillUAFromTracks->SetHadCorrection(fHCorrection);
 242   fFillUAFromTracks->SetHadCorrector(fHadCorr);
 243   fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
 244   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
 245   fFillUAFromEMCalDigits->SetGeom(fGeom);
 246   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
 247   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
 248   //      fFillUnitArray->Add(fFillUAFromTracks);
 249   fFillUnitArray->Add(fFillUAFromEMCalDigits);
 250   fFillUAFromTracks->SetActive(kFALSE);
 251   fFillUAFromEMCalDigits->SetActive(kFALSE);
 252
 253   cout << "Tasks instantiated at that stage ! " << endl;
 254   cout << "You can loop over events now ! " << endl;
 255
 256 }
 257
 258 //____________________________________________________________________________
 259 //void AliJetESDReader::ExecTasks(Int_t event)
 260 Bool_t AliJetESDReader::ExecTasks(Int_t /*event*/)
 261 {
 262   // clear momentum array
 263   Int_t nEntRef = fRefArray->GetEntries();
 264
 265   for(Int_t i=0; i<nEntRef; i++)
 266     {
 267       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitTrackID(0);
 268       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitEnergy(0.);
 269       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitCutFlag(kPtSmaller);
 270       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitDetectorFlag(kTpc);
 271       ((AliJetUnitArray*)fRefArray->At(i))->SetUnitFlag(kOutJet);
 272     }
 273
 274   ClearArray();
 275
 276   fDebug = fReaderHeader->GetDebug();
 277   fOpt = fReaderHeader->GetDetector();
 278   //  InitParameters();
 279
 280   if(!fESD) {
 281     return kFALSE;
 282   }
 283
 284   /*
 285   // get event from chain
 286   // For TSelectors
 287   //  fChain->GetTree()->GetEntry(event);
 288   // For interactive process
 289   //  fChain->GetEntry(event);
 290   fChain->GetEvent(event);
 291   */
 292
 293   // TPC only or Digits+TPC or Clusters+TPC
 294   if(fOpt%2==!0 && fOpt!=0){
 295     fFillUAFromTracks->SetESD(fESD);
 296     fFillUAFromTracks->SetActive(kTRUE);
 297     fFillUAFromTracks->SetUnitArray(fUnitArray);
 298     fFillUAFromTracks->SetRefArray(fRefArray);
 299     //    fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed  !!!
 300                                                  // Incompatibility with Andrei's analysis framework
 301     fFillUAFromTracks->Exec("tpc");
 302     if(fOpt==1){
 303       fNumCandidate = fFillUAFromTracks->GetMult();
 304       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 305     }
 306   }
 307
 308   // Digits only or Digits+TPC
 309   if(fOpt>=2 && fOpt<=3){
 310     fFillUAFromEMCalDigits->SetESD(fESD);
 311     fFillUAFromEMCalDigits->SetActive(kTRUE);
 312     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 313     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 314     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 315     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 316     fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily changed !!!
 317     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 318     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 319   }
 320
 321   //  fFillUnitArray->ExecuteTask(); // => Temporarily commented
 322
 323   return kTRUE;
 324 }
 325
 326 //____________________________________________________________________________
 327 void AliJetESDReader::SetEMCALGeometry()
 328 {
 329   // Define EMCAL geometry to be able to read ESDs
 330     fGeom = AliJetDummyGeo::GetInstance();
 331     if (fGeom == 0)
 332         fGeom = AliJetDummyGeo::GetInstance("SHISH_77_TRD1_2X2_FINAL_110DEG","EMCAL");
 333
 334     // To be setted to run some AliEMCALGeometry functions
 335     TGeoManager::Import("geometry.root");
 336     //    fGeom->GetTransformationForSM();
 337     printf("\n EMCal Geometry set ! \n");
 338
 339 }
 340
 341
 342 //____________________________________________________________________________
 343 void AliJetESDReader::InitParameters()
 344 {
 345     // Initialise parameters
 346     fHCorrection    = 0;                 // For hadron correction
 347     fHadCorr        = 0;                 // For hadron correction
 348     fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 349     if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 350 }
 351
 352 //____________________________________________________________________________
 353 void AliJetESDReader::InitUnitArray()
 354 {
 355   //Initialises unit arrays
 356   Int_t nElements = fTpcGrid->GetNEntries();
 357   Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
 358   if(fArrayInitialised) fUnitArray->Delete();
 359
 360   if(fTpcGrid->GetGridType()==0)
 361     { // Fill the following quantities :
 362       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
 363       // detector flag, in/out jet, pt cut, mass, cluster ID)
 364       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 365         {
 366           //      fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
 367           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 368           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 369           Deta = fTpcGrid->GetDeta();
 370           Dphi = fTpcGrid->GetDphi();
 371           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 372         }
 373     }
 374
 375   if(fTpcGrid->GetGridType()==1)
 376     {
 377       Int_t nGaps = 0;
 378       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 379
 380       if(fDZ)
 381         {
 382           // Define a grid of cell for the gaps between SM
 383           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 384           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 385           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 386           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 387           fGrid0->SetGridType(0);
 388           fGrid0->SetMatrixIndexes();
 389           fGrid0->SetIndexIJ();
 390           n0 = fGrid0->GetNEntries();
 391           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 392           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 393           fGrid1->SetGridType(0);
 394           fGrid1->SetMatrixIndexes();
 395           fGrid1->SetIndexIJ();
 396           n1 = fGrid1->GetNEntries();
 397           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 398           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 399           fGrid2->SetGridType(0);
 400           fGrid2->SetMatrixIndexes();
 401           fGrid2->SetIndexIJ();
 402           n2 = fGrid2->GetNEntries();
 403           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 404           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 405           fGrid3->SetGridType(0);
 406           fGrid3->SetMatrixIndexes();
 407           fGrid3->SetIndexIJ();
 408           n3 = fGrid3->GetNEntries();
 409           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 410           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 411           fGrid4->SetGridType(0);
 412           fGrid4->SetMatrixIndexes();
 413           fGrid4->SetIndexIJ();
 414           n4 = fGrid4->GetNEntries();
 415
 416           if(fDebug>1)
 417             {
 418               cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
 419               cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl;
 420               cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl;
 421               cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl;
 422               cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl;
 423             }
 424
 425           nGaps = n0+n1+n2+n3+n4;
 426
 427         }
 428
 429       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 430         {
 431           if(nBin<fNumUnits)
 432             {
 433               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 434               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 435               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 436               Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 437               Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 438               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 439             }
 440           else {
 441             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 442               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 443               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 444               Deta = fTpcGrid->GetDeta();
 445               Dphi = fTpcGrid->GetDphi();
 446               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 447             }
 448             else {
 449               if(fDZ) {
 450                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 451                   if(nBin<fNumUnits+nElements+n0)
 452                     {
 453                       Float_t phi = eta = 0.;
 454                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
 455                       Deta = fGrid0->GetDeta();
 456                       Dphi = fGrid0->GetDphi();
 457                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 458                     }
 459                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 460                     {
 461                       Float_t phi = eta = 0.;
 462                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
 463                       Deta = fGrid1->GetDeta();
 464                       Dphi = fGrid1->GetDphi();
 465                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 466                     }
 467                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 468                     {
 469                       Float_t phi = eta = 0.;
 470                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
 471                       Deta = fGrid2->GetDeta();
 472                       Dphi = fGrid2->GetDphi();
 473                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 474                     }
 475                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 476                     {
 477                       Float_t phi = eta = 0.;
 478                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
 479                       Deta = fGrid3->GetDeta();
 480                       Dphi = fGrid3->GetDphi();
 481                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 482                     }
 483                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 484                     {
 485                       Float_t phi = eta = 0.;
 486                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
 487                       Deta = fGrid4->GetDeta();
 488                       Dphi = fGrid4->GetDphi();
 489                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,0.,0.,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,0.,-1);
 490                     }
 491                 }
 492               } // end if(fDZ)
 493             } // end else 2
 494           } // end else 1
 495         } // end loop on nBin
 496     } // end grid type == 1
 497   fArrayInitialised = 1;
 498 }
 499
 500
 501