ITS/EvaluateSPDEffWithTracklets.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2
   3 #include <Riostream.h>
   4 #include <TFile.h>
   5 #include <TTree.h>
   6 #include <TBranch.h>
   7 #include <TROOT.h>
   8 #include <TObjArray.h>
   9 #include <TGeoManager.h>
  10 #include <TProfile.h>
  11
  12 #include "AliESD.h"
  13 #include "AliESDVertex.h"
  14 #include "AliESDEvent.h"
  15
  16 #include "AliRun.h"
  17 #include "AliRunLoader.h"
  18
  19 #include "AliCDBManager.h"
  20 #include "AliGeomManager.h"
  21 #include "AliITS.h"
  22 #include "AliITSgeom.h"
  23 #include "AliITSLoader.h"
  24 #include "AliITSRecPoint.h"
  25
  26 #include "AliITSTrackleterSPDEff.h"
  27 #include "AliITSPlaneEffSPD.h"
  28
  29 #include <AliHeader.h>
  30 #include <AliGenEventHeader.h>
  31 #include <AliStack.h>
  32
  33 #endif
  34
  35 void EvaluateSPDEffWithTracklets(Char_t* dir=".", Bool_t mc=kTRUE, Bool_t bckg=kFALSE) {
  36
  37   Char_t str[256];
  38
  39   AliITSTrackleterSPDEff* trackleterSPDEff = new AliITSTrackleterSPDEff();
  40 // outer layer (estrapolation)
  41   trackleterSPDEff->SetPhiWindow(0.07);
  42   trackleterSPDEff->SetZetaWindow(0.4);
  43 // inner layer (interpolation)
  44   trackleterSPDEff->SetPhiWindowL1(0.10);
  45   trackleterSPDEff->SetZetaWindowL1(0.6);
  46 //
  47   trackleterSPDEff->SetUpdateOncePerEventPlaneEff();
  48 // Study the residual background: reflect outer RecPoints
  49   if(bckg) trackleterSPDEff->SetReflectClusterAroundZAxisForLayer(1,kTRUE);
  50 //
  51 // this special setting for MC
  52   if(mc) trackleterSPDEff->SetMC();
  53   if(trackleterSPDEff->GetMC()) trackleterSPDEff->SetUseOnlyStableParticle();
  54 //
  55 // this for having histograms (both from base class and the new ones)
  56   trackleterSPDEff->SetHistOn();
  57 //
  58 //
  59   const Int_t minCont=3;
  60   const Bool_t VtxMC=kFALSE;
  61 //
  62   // Defining pointers
  63   AliRunLoader* runLoader;
  64
  65     if (gAlice) {
  66       delete gAlice->GetRunLoader();
  67       delete gAlice;
  68       gAlice=0;
  69     }
  70
  71     sprintf(str,"%s/galice.root",dir);
  72     runLoader = AliRunLoader::Open(str);
  73     runLoader->LoadgAlice();
  74     gAlice = runLoader->GetAliRun();
  75
  76     runLoader->LoadKinematics("read");
  77     runLoader->LoadTrackRefs("read");
  78     Int_t retval = runLoader->LoadHeader();
  79     if (retval){
  80       cerr<<"LoadHeader returned error"<<endl;
  81       return;
  82     }
  83
  84     // open the new ESD file
  85     sprintf(str,"%s/AliESDs.root",dir);
  86
  87     TFile inFile(str, "READ");
  88     TTree *esdTree = (TTree*)inFile.Get("esdTree");
  89     AliESDEvent *esd = new AliESDEvent();
  90     esd->ReadFromTree(esdTree);
  91
  92     // Set OfflineConditionsDataBase if needed
  93     AliCDBManager* man = AliCDBManager::Instance();
  94     if (!man->IsDefaultStorageSet()) {
  95       printf("Setting a local default storage and run number 0\n");
  96       man->SetDefaultStorage("local://$ALICE_ROOT");
  97       man->SetRun(0);
  98     }
  99     else {
 100       printf("Using deafult storage \n");
 101     }
 102     // retrives geometry
 103     if (!gGeoManager) {
 104       sprintf(str,"%s/geometry.root",dir);
 105       AliGeomManager::LoadGeometry(str);
 106     }
 107     AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
 108     if (!ITSloader){
 109       cerr<<"ITS loader not found"<<endl;
 110       return;
 111     }
 112     ITSloader->LoadRecPoints("read");
 113
 114     // getting number of events
 115     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
 116
 117     // loop over events
 118     for (Int_t iev=0; iev<nEvents; iev++) {
 119
 120       runLoader->GetEvent(iev);
 121       // read events
 122       esdTree->GetEvent(iev);
 123
 124       // get the ESD vertex
 125       const AliESDVertex* vtxESD = esd->GetVertex();
 126       Double_t esdvtx[3];
 127       vtxESD->GetXYZ(esdvtx);
 128       Int_t ncont=vtxESD->GetNContributors();
 129       if(ncont <= minCont) continue;
 130       Float_t ESDvtx[3];
 131       ESDvtx[0]=esdvtx[0];
 132       ESDvtx[1]=esdvtx[1];
 133       ESDvtx[2]=esdvtx[2];
 134
 135       // get the MC vertex
 136       TArrayF vertex(3);
 137       runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 138
 139      // Read the generated particles
 140      AliStack *pStack=0x0; TTree *tRefTree=0x0;
 141      if (trackleterSPDEff->GetMC()) {
 142        pStack= runLoader->Stack();
 143        tRefTree= runLoader->TreeTR();
 144      }
 145
 146      TTree *itsClusterTree = ITSloader->TreeR();
 147
 148      if(!VtxMC) {
 149       if (ESDvtx[2]!=0.) {
 150         if(trackleterSPDEff->GetMC()) trackleterSPDEff->Reconstruct(itsClusterTree, ESDvtx, ESDvtx, pStack,tRefTree);
 151         else trackleterSPDEff->Reconstruct(itsClusterTree, ESDvtx, ESDvtx); }
 152      }
 153      else {
 154        Float_t vtx[3]={0.,0.,vertex[2]};
 155        if(trackleterSPDEff->GetMC()) trackleterSPDEff->Reconstruct(itsClusterTree, vtx, vtx, pStack,tRefTree);
 156      }
 157
 158    } // end loop over events
 159
 160    runLoader->UnloadAll();
 161    delete runLoader;
 162
 163 if(trackleterSPDEff->GetMC()) trackleterSPDEff->SavePredictionMC();
 164 if(!trackleterSPDEff->WriteHistosToFile()) printf("cannot write histos to file \n");
 165 //trackleterSPDEff->GetPlaneEff()->WriteIntoCDB();
 166 const char* name="AliITSPlaneEffSPDtracklet.root";
 167 TFile* pefile = TFile::Open(name, "RECREATE");
 168 Int_t nb=trackleterSPDEff->GetPlaneEff()->Write();
 169 if(nb>0) printf("Writing PlaneEfficiency to file %s\n",name);
 170 pefile->Close();
 171 return;
 172 }