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