MC info can now be written also into root file for running with Task on Proof (G...
[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
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("TrackletsMCpred.root");
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 }