Changed AliRunLoader::GetRunLoader() into AliRunLoader::Instance()
[u/mrichter/AliRoot.git] / 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");
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 }