new macros
[u/mrichter/AliRoot.git] / ITS / EvaluateSPDEffWithTracklets.C
CommitLineData
2065c9d5 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
35void 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
163if(trackleterSPDEff->GetMC()) trackleterSPDEff->SavePredictionMC();
164if(!trackleterSPDEff->WriteHistosToFile()) printf("cannot write histos to file \n");
165trackleterSPDEff->GetPlaneEff()->WriteIntoCDB();
166return;
167}