A more elaborate approach for associating the reconstructed and MC pileup vertices
[u/mrichter/AliRoot.git] / ITS / UPGRADE / macros / QA / AliITSUComparisonPileup.C
CommitLineData
394f72ea 1/****************************************************************************
2 * This macro calculates the efficiency of pileup reconstruction. *
3 * Works only for events generated with the AliGenPileup generator. *
4 * *
5 * Before running, load the ITSU libraries: *
6 * gSystem->Load("libITSUpgradeBase");gSystem->Load("libITSUpgradeRec"); *
7 * *
8 * Origin: I.Belikov, IPHC, Iouri.Belikov@iphc.cnrs.fr *
9 ****************************************************************************/
10
11#if !defined(__CINT__) || defined(__MAKECINT__)
12 #include <Riostream.h>
13 #include <TMath.h>
14 #include <TTree.h>
15 #include <TParticle.h>
c4b8f157 16 #include <TParticlePDG.h>
394f72ea 17 #include <TCanvas.h>
18 #include <TFile.h>
c4b8f157 19 #include <TLine.h>
394f72ea 20 #include <TROOT.h>
21
22 #include "AliStack.h"
23 #include "AliHeader.h"
24 #include "AliGenCocktailEventHeader.h"
25 #include "AliRunLoader.h"
26 #include "AliRun.h"
27 #include "AliESDEvent.h"
28 #include "AliESDtrack.h"
29#endif
30
31Int_t GoodPileupVertices(const Char_t *dir=".");
32
33extern AliRun *gAlice;
34extern TROOT *gROOT;
35
36Int_t AliITSUComparisonPileup(const Char_t *dir=".") {
37 ::Info("AliITSUComparisonPileup.C","Doing comparison...");
38
39 // **** Book histogramms
e4607b46 40 Int_t nb=35;
41 Float_t min=0, max=70.;
c4b8f157 42 TH2F *h2spd=(TH2F*)gROOT->FindObject("h2spd");
43 if (!h2spd)
44 h2spd=new TH2F("h2spd",";Good vertices;Reconstructed vertices",
45 nb,min,max, nb,min,max);
46 h2spd->SetLineColor(2);
47 TH2F *h2trk=(TH2F*)gROOT->FindObject("h2trk");
48 if (!h2trk)
49 h2trk=new TH2F("h2trk",";Good vertices;Reconstructed vertices",
50 nb,min,max, nb,min,max);
51 h2trk->SetLineColor(4);
52
53 nb=100;
54 min=-0.03; max=0.03;
55 TH1F *hzspd=(TH1F*)gROOT->FindObject("hzspd");
56 if (!hzspd)
57 hzspd=new TH1F("hzspd","Resolution in Z;#DeltaZ (cm);",nb,min,max);
58 hzspd->SetLineColor(2);
59 TH1F *hztrk=(TH1F*)gROOT->FindObject("hztrk");
60 if (!hztrk)
61 hztrk=new TH1F("hztrk","Resolution in Z;#DeltaZ (cm);",nb,min,max);
62 hztrk->SetLineColor(4);
63
64
65 nb=30;
66 min=-10.; max=10.;
67 TH1F *hgood=(TH1F*)gROOT->FindObject("hgood");
68 if (!hgood)
69 hgood=new TH1F("hgood",";Z (cm);",nb,min,max);
70 TH1F *hfoundspd=(TH1F*)gROOT->FindObject("hfoundspd");
71 if (!hfoundspd)
72 hfoundspd=new TH1F("hfoundspd",";Z (cm);",nb,min,max);
73 TH1F *heffspd=(TH1F*)gROOT->FindObject("heffspd");
74 if (!heffspd)
75 heffspd=new TH1F("heffspd","Efficiency;Z (cm);Efficiency",nb,min,max);
76 heffspd->SetLineColor(2);
77 TH1F *hfoundtrk=(TH1F*)gROOT->FindObject("hfoundtrk");
78 if (!hfoundtrk)
79 hfoundtrk=new TH1F("hfoundtrk",";Z (cm);",nb,min,max);
80 TH1F *hefftrk=(TH1F*)gROOT->FindObject("hefftrk");
81 if (!hefftrk)
82 hefftrk=new TH1F("hefftrk",";Z (cm);Efficiency",nb,min,max);
83 hefftrk->SetLineColor(4);
394f72ea 84
85
86 // **** Generate a rerefence file with reconstructable vertices
87 Char_t fname[100];
88 sprintf(fname,"%s/GoodPileupVertices.root",dir);
89 TFile *refFile=TFile::Open(fname,"old");
90 if (!refFile || !refFile->IsOpen()) {
91 ::Info("AliITSUComparisonPileup.C",
92 "Marking good pileup vertices (will take a while)...");
93 if (GoodPileupVertices(dir)) {
94 ::Error("AliITSUComparisonPileup.C","Can't generate the reference file !");
95 return 1;
96 }
97 }
98 refFile=TFile::Open(fname,"old");
99 if (!refFile || !refFile->IsOpen()) {
100 ::Error("AliITSUComparisonPileup.C","Can't open the reference file !");
101 return 1;
102 }
103 TTree *refTree=(TTree*)refFile->Get("refTree");
104 if (!refTree) {
105 ::Error("AliITSUComparisonPileup.C","Can't get the reference tree !");
106 return 2;
107 }
108 TBranch *branch=refTree->GetBranch("Vertices");
109 if (!branch) {
110 ::Error("AliITSUComparisonPileup.C","Can't get the vertex branch !");
111 return 3;
112 }
113 TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
114 branch->SetAddress(&refs);
115
116
117 // **** Open the ESD
118 sprintf(fname,"%s/AliESDs.root",dir);
119 TFile *ef=TFile::Open(fname);
120 if ((!ef)||(!ef->IsOpen())) {
121 ::Error("AliITSUComparisonPileup.C","Can't open AliESDs.root !");
122 return 4;
123 }
124 AliESDEvent* event = new AliESDEvent();
125 TTree* esdTree = (TTree*) ef->Get("esdTree");
126 if (!esdTree) {
127 ::Error("AliITSUComparisonPileup.C", "no ESD tree found");
128 return 6;
129 }
130 event->ReadFromTree(esdTree);
131
132
c4b8f157 133 //******* Loop over reconstructed events *********
394f72ea 134 Int_t e=0;
135 while (esdTree->GetEvent(e)) {
136 cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
c4b8f157 137 TClonesArray *verticesSPD=event->GetPileupVerticesSPD();
138 Int_t nfoundSPD=verticesSPD->GetEntries();
139 TClonesArray *verticesTRK=event->GetPileupVerticesTracks();
140 Int_t nfoundTRK=verticesTRK->GetEntries();
394f72ea 141
142 if (refTree->GetEvent(e++)==0) {
143 cerr<<"No reconstructable vertices for this event !\n";
144 continue;
145 }
146 Int_t ngood=refs->GetEntriesFast();
c4b8f157 147 cout<<"Found SPD: "<<nfoundSPD<<" good: "<<ngood<<endl;
148
149 h2spd->Fill(ngood,nfoundSPD);
150 h2trk->Fill(ngood,nfoundTRK);
151
152 for (Int_t g=0; g<ngood; g++) {
e4607b46 153 Bool_t Associate(const AliESDVertex *g, const AliESDVertex *f,
154 const AliESDEvent *esd);
c4b8f157 155 const AliESDVertex *vtxg=(AliESDVertex*)refs->UncheckedAt(g);
156 Double_t zg=vtxg->GetZv();
157 hgood->Fill(zg);
158
159 const AliESDVertex *vtxf=0;
160 Double_t zf=0.;
161 Int_t f=0;
162 for (; f<nfoundSPD; f++) {
163 vtxf=(AliESDVertex*)verticesSPD->UncheckedAt(f);
164 if (!vtxf->GetStatus()) continue;
e4607b46 165 if (!Associate(vtxg,vtxf,event)) continue;
c4b8f157 166 zf=vtxf->GetZv();
c4b8f157 167 break;
168 }
169 if (f>=nfoundSPD) {
170 vtxf=event->GetPrimaryVertexSPD();
171 if (!vtxf->GetStatus()) goto trk;
e4607b46 172 if (!Associate(vtxg,vtxf,event)) goto trk;
c4b8f157 173 zf=vtxf->GetZv();
c4b8f157 174 }
175 hfoundspd->Fill(zg);
176 hzspd->Fill(zf-zg);
177
178 trk:
179 for (f=0; f<nfoundTRK; f++) {
180 vtxf=(AliESDVertex*)verticesTRK->UncheckedAt(f);
181 if (!vtxf->GetStatus()) continue;
e4607b46 182 if (!Associate(vtxg,vtxf,event)) continue;
c4b8f157 183 zf=vtxf->GetZv();
c4b8f157 184 break;
185 }
186 if (f>=nfoundTRK) {
187 vtxf=event->GetPrimaryVertexTracks();
188 if (!vtxf->GetStatus()) continue;
e4607b46 189 if (!Associate(vtxg,vtxf,event)) continue;
c4b8f157 190 zf=vtxf->GetZv();
c4b8f157 191 }
192 hfoundtrk->Fill(zg);
193 hztrk->Fill(zf-zg);
394f72ea 194
c4b8f157 195 }
394f72ea 196
c4b8f157 197 } //***** End of the loop over reconstructed events
394f72ea 198
199 delete event;
200 delete esdTree;
201 ef->Close();
202
203 refFile->Close();
204
c4b8f157 205 TCanvas *c1=new TCanvas("c1","",0,0,700,1000);
206 c1->Divide(1,3);
207 c1->cd(1);
208 gPad->SetGridx(); gPad->SetGridy();
209 h2spd->Draw("box");
210 h2trk->Draw("boxsame");
e4607b46 211 TLine *l=new TLine(0,0,60,60);
c4b8f157 212 l->Draw("same");
213
214 c1->cd(2);
215 gPad->SetGridx(); gPad->SetGridy();
216 heffspd->Divide(hfoundspd,hgood,1,1,"b");
e4607b46 217 heffspd->SetMinimum(0.); heffspd->SetMaximum(1.2);
c4b8f157 218 heffspd->Draw("hist");
219 hefftrk->Divide(hfoundtrk,hgood,1,1,"b");
220 hefftrk->Draw("histsame");
221
222 c1->cd(3);
223 gPad->SetGridx(); gPad->SetGridy();
e4607b46 224 hztrk->Draw();
225 hzspd->Draw("same");
c4b8f157 226
227 TFile fc("AliITSUComparisonPileup.root","RECREATE");
394f72ea 228 c1->Write();
229 fc.Close();
230
231 return 0;
232}
233
e4607b46 234Bool_t
235Associate(const AliESDVertex *g,const AliESDVertex *f,const AliESDEvent *esd) {
236 UShort_t *idxg=g->GetIndices(); Int_t ng=g->GetNIndices();
237 UShort_t *idxf=f->GetIndices(); Int_t nf=f->GetNIndices();
238
239 if (nf==0) {
240 // SPD vertex
241 Double_t zg=g->GetZv();
242 Double_t zf=f->GetZv();
243 if (TMath::Abs(zf-zg)>2e-2) return kFALSE;
244 return kTRUE;
245 }
246 // TRK vertex
247 Int_t nass=0;
248 for (Int_t i=0; i<ng; i++) {
249 UShort_t labg=idxg[i];
250 for (Int_t j=0; j<nf; j++) {
251 const AliESDtrack *t=esd->GetTrack(idxf[j]);
252 UShort_t labf=TMath::Abs(t->GetLabel());
253 if (labg != labf) continue;
254 nass++;
255 break;
256 }
257 }
258
259 if (Float_t(nass)/ng > 0.5) return kTRUE;
260 return kFALSE;
261}
394f72ea 262
263Int_t GoodPileupVertices(const Char_t *dir) {
e4607b46 264 Int_t FindContributors(Float_t tz, AliStack *stack, UShort_t *idx);
394f72ea 265 if (gAlice) {
266 delete AliRunLoader::Instance();
267 delete gAlice;//if everything was OK here it is already NULL
268 gAlice = 0x0;
269 }
270
271 Char_t fname[100];
272 sprintf(fname,"%s/galice.root",dir);
273
274 AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
275 if (!rl) {
276 ::Error("GoodPileupVertices","Can't start session !");
277 return 1;
278 }
279
280 rl->LoadgAlice();
281 rl->LoadHeader();
282 rl->LoadKinematics();
283
284
285 Int_t nev=rl->GetNumberOfEvents();
286 ::Info("GoodPileupVertices","Number of events : %d\n",nev);
287
288 sprintf(fname,"%s/GoodPileupVertices.root",dir);
289 TFile *refFile=TFile::Open(fname,"recreate");
290 TClonesArray dummy("AliESDVertex",100), *refs=&dummy;
291 TTree refTree("refTree","Tree with the reconstructable vertices");
292 refTree.Branch("Vertices",&refs);
293
294 //******** Loop over generated events
295 for (Int_t e=0; e<nev; e++) {
296 rl->GetEvent(e); refFile->cd();
c4b8f157 297 AliStack* stack = rl->Stack();
394f72ea 298
299 AliHeader *ah=rl->GetHeader();
300 AliGenCocktailEventHeader *cock=
301 (AliGenCocktailEventHeader*)ah->GenEventHeader();
302 TList *headers=cock->GetHeaders();
303 const Int_t nvtx=headers->GetEntries();
c4b8f157 304 const Int_t np=ah->GetNtrack();
305 cout<<"Event "<<e<<" Number of vertices, particles: "
306 <<nvtx<<' '<<np<<endl;
394f72ea 307
308 Int_t nv=0;
309 for (Int_t v=0; v<nvtx; v++) {
310 AliGenEventHeader *h=(AliGenEventHeader *)headers->At(v);
311 TArrayF vtx(3); h->PrimaryVertex(vtx);
c4b8f157 312 Float_t t=h->InteractionTime();
e4607b46 313 UShort_t *idx=new UShort_t[np];
314 Int_t ntrk=FindContributors(t,stack,idx);
c4b8f157 315 if (ntrk<3) continue;
394f72ea 316 AliESDVertex *vertex=new ((*refs)[nv]) AliESDVertex();
317 vertex->SetXv(vtx[0]);
318 vertex->SetYv(vtx[1]);
319 vertex->SetZv(vtx[2]);
c4b8f157 320 vertex->SetNContributors(ntrk);
e4607b46 321 vertex->SetIndices(ntrk,idx);
322 delete idx;
394f72ea 323 nv++;
324 }
325 refTree.Fill();
326 refs->Clear();
327 } //*** end of the loop over generated events
328
329 refTree.Write();
330 refFile->Close();
331
332 delete rl;
333 return 0;
334}
335
e4607b46 336Int_t FindContributors(Float_t tz, AliStack *stack, UShort_t *idx) {
c4b8f157 337 Int_t ntrk=0;
e4607b46 338 Int_t np=stack->GetNtrack();
c4b8f157 339 for (Int_t i=0; i<np; i++) {
340 if (!stack->IsPhysicalPrimary(i)) continue;
341 TParticle *part=stack->Particle(i);
342 if (!part) continue;
343 TParticlePDG *partPDG = part->GetPDG();
344 if (!partPDG) continue;
345 if (TMath::Abs(partPDG->Charge())<1e-10) continue;
e4607b46 346 if (part->Pt()<0.5) continue;
c4b8f157 347 Float_t dt=0.5*(tz-part->T())/(tz+part->T());
348 if (TMath::Abs(dt)>1e-5) continue;
e4607b46 349 idx[ntrk]=i;
c4b8f157 350 ntrk++;
351 }
352 return ntrk;
353}