ca28c5f5 |
1 | /**************************************************************************** |
2 | * Very important, delicate and rather obscure macro. * |
3 | * * |
4 | * Creates list of "findable" cascades, * |
5 | * calculates efficiency, resolutions etc. * |
6 | * * |
7 | * Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr * |
8 | ****************************************************************************/ |
9 | |
10 | #ifndef __CINT__ |
11 | #include <iostream.h> |
12 | #include <fstream.h> |
13 | |
14 | #include "TH1.h" |
15 | #include "TFile.h" |
16 | #include "TTree.h" |
17 | #include "TObjArray.h" |
18 | #include "TStyle.h" |
19 | #include "TCanvas.h" |
20 | #include "TLine.h" |
21 | #include "TText.h" |
22 | #include "TParticle.h" |
23 | #include "TStopwatch.h" |
24 | |
25 | #include "AliRun.h" |
26 | #include "AliPDG.h" |
27 | #include "AliCascadeVertex.h" |
28 | #endif |
29 | |
30 | struct GoodCascade { |
31 | Int_t nlab,plab; // V0's daughter labels |
32 | Int_t blab; // Bachelor label |
33 | Int_t code; |
34 | Float_t px,py,pz; |
35 | Float_t x,y,z; |
36 | }; |
37 | Int_t good_cascades(GoodCascade *gt, Int_t max); |
38 | |
39 | Int_t AliCascadeComparison(Int_t code=3312) { |
40 | //code= 3312; //kXiMinus |
41 | //code=-3312; //kXiPlusBar |
42 | //code= 3334; //kOmegaMinus |
43 | //code=-3334; //kOmegaPlusBar |
44 | |
45 | cerr<<"Doing comparison...\n"; |
46 | |
47 | const Double_t cascadeWindow=0.05, cascadeWidth=0.015; |
04b2a5f1 |
48 | Double_t ptncut=0.12, ptpcut=0.33, kine0cut=0.003; |
49 | Double_t ptbcut=0.11, kinecut=0.002; |
50 | Double_t cascadeMass=1.32131; |
ca28c5f5 |
51 | switch (code) { |
52 | case kXiMinus: |
04b2a5f1 |
53 | break; |
54 | case kXiPlusBar: |
55 | ptncut=0.33; ptpcut=0.12; |
56 | break; |
ca28c5f5 |
57 | case kOmegaMinus: |
04b2a5f1 |
58 | cascadeMass=1.67245; |
59 | kine0cut=0.001; |
60 | ptbcut=0.22; kinecut=0.006; |
61 | break; |
62 | case kOmegaPlusBar: |
63 | cascadeMass=1.67245; |
64 | kine0cut=0.001; |
65 | ptncut=0.33; ptpcut=0.12; |
66 | ptbcut=0.22; kinecut=0.006; |
67 | break; |
ca28c5f5 |
68 | default: cerr<<"Invalid PDG code !\n"; return 1; |
69 | } |
70 | |
71 | TStopwatch timer; |
72 | |
73 | /*** Load reconstructed cascades ***/ |
74 | TFile *cf=TFile::Open("AliCascadeVertices.root"); |
75 | if (!cf->IsOpen()){cerr<<"Can't open AliCascadeVertices.root !\n";return 2;} |
76 | TObjArray carray(1000); |
77 | TTree *cTree=(TTree*)cf->Get("TreeCasc"); |
78 | TBranch *branch=cTree->GetBranch("cascades"); |
79 | Int_t nentr=(Int_t)cTree->GetEntries(); |
80 | for (Int_t i=0; i<nentr; i++) { |
81 | AliCascadeVertex *iovertex=new AliCascadeVertex; |
82 | branch->SetAddress(&iovertex); |
83 | cTree->GetEvent(i); |
84 | carray.AddLast(iovertex); |
85 | } |
86 | delete cTree; |
87 | cf->Close(); |
88 | |
89 | /*** Check if the file with the "good" cascades exists ***/ |
90 | GoodCascade gc[100]; |
91 | Int_t ngood=0; |
92 | ifstream in("good_cascades"); |
93 | if (in) { |
94 | cerr<<"Reading good cascades...\n"; |
95 | while (in>>gc[ngood].nlab>>gc[ngood].plab>> |
96 | gc[ngood].blab>>gc[ngood].code>> |
97 | gc[ngood].px>>gc[ngood].py>>gc[ngood].pz>> |
98 | gc[ngood].x >>gc[ngood].y >>gc[ngood].z) { |
99 | ngood++; |
100 | cerr<<ngood<<'\r'; |
101 | if (ngood==100) { |
102 | cerr<<"Too many good cascades !\n"; |
103 | break; |
104 | } |
105 | } |
106 | if (!in.eof()) cerr<<"Read error (good_cascades) !\n"; |
107 | } else { |
108 | /*** generate a file with the "good" cascades ***/ |
109 | cerr<<"Marking good cascades (this will take a while)...\n"; |
110 | ngood=good_cascades(gc,100); |
111 | ofstream out("good_cascades"); |
112 | if (out) { |
113 | for (Int_t ngd=0; ngd<ngood; ngd++) |
114 | out<<gc[ngd].nlab<<' '<<gc[ngd].plab<<' '<< |
115 | gc[ngd].blab<<' '<<gc[ngd].code<<' '<< |
116 | gc[ngd].px<<' '<<gc[ngd].py<<' '<<gc[ngd].pz<<' '<< |
117 | gc[ngd].x <<' '<<gc[ngd].y <<' '<<gc[ngd].z <<endl; |
118 | } else cerr<<"Can not open file (good_cascades) !\n"; |
119 | out.close(); |
120 | } |
121 | |
122 | /*** create some histograms ***/ |
123 | TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution |
124 | hp->SetXTitle("(mrad)"); hp->SetFillColor(2); |
125 | TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30); |
126 | hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013); |
127 | TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.); |
128 | hpt->SetXTitle("(%)"); hpt->SetFillColor(2); |
129 | |
130 | TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res. |
131 | hx->SetXTitle("(mm)"); hx->SetFillColor(6); |
132 | TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.); //y res |
133 | hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013); |
134 | TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.); //z res. |
135 | hz->SetXTitle("(mm)"); hz->SetFillColor(6); |
136 | |
137 | Double_t pmin=0.2, pmax=4.2; Int_t nchan=20; |
138 | TH1F *hgood=new TH1F("hgood","Good Cascades",nchan,pmin,pmax); |
139 | TH1F *hfound=new TH1F("hfound","Found Cascades",nchan,pmin,pmax); |
140 | TH1F *hfake=new TH1F("hfake","Fake Cascades",nchan,pmin,pmax); |
141 | TH1F *hg=new TH1F("hg","Efficiency for Good Cascades",nchan,pmin,pmax); |
142 | hg->SetLineColor(4); hg->SetLineWidth(2); |
143 | TH1F *hf=new TH1F("hf","Probability of Fake Cascades",nchan,pmin,pmax); |
144 | hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); |
145 | |
146 | Double_t mmin=cascadeMass-cascadeWindow, mmax=cascadeMass+cascadeWindow; |
147 | TH1F *cs =new TH1F("cs","Cascade Effective Mass",40, mmin, mmax); |
04b2a5f1 |
148 | cs->SetXTitle("(GeV)"); |
149 | cs->SetLineColor(4); cs->SetLineWidth(4); |
150 | TH1F *csf =new TH1F("csf","Fake Cascade Effective Mass",40, mmin, mmax); |
151 | csf->SetXTitle("(GeV)"); csf->SetFillColor(6); |
ca28c5f5 |
152 | |
153 | Double_t pxg=0.,pyg=0.,ptg=0.; |
154 | Int_t nlab=-1, plab=-1, blab=-1; |
155 | Int_t i; |
156 | for (i=0; i<nentr; i++) { |
157 | AliCascadeVertex *cascade=(AliCascadeVertex*)carray.UncheckedAt(i); |
04b2a5f1 |
158 | nlab=TMath::Abs(cascade->GetNlabel()); |
159 | plab=TMath::Abs(cascade->GetPlabel()); |
160 | blab=TMath::Abs(cascade->GetBlabel()); |
161 | |
162 | /** Kinematical cuts **/ |
163 | Double_t pxn,pyn,pzn; cascade->GetNPxPyPz(pxn,pyn,pzn); |
164 | Double_t ptn=TMath::Sqrt(pxn*pxn + pyn*pyn); |
165 | if (ptn < ptncut) continue; |
166 | Double_t pxp,pyp,pzp; cascade->GetPPxPyPz(pxp,pyp,pzp); |
167 | Double_t ptp=TMath::Sqrt(pxp*pxp + pyp*pyp); |
168 | if (ptp < ptpcut) continue; |
169 | Double_t pxb,pyb,pzb; cascade->GetBPxPyPz(pxb,pyb,pzb); |
170 | Double_t ptb=TMath::Sqrt(pxb*pxb + pyb*pyb); |
171 | if (ptb < ptbcut) continue; |
172 | Double_t kine0; |
173 | Double_t kine=cascade->ChangeMassHypothesis(kine0,code); |
174 | if (TMath::Abs(kine0)>kine0cut) continue; |
175 | //if (TMath::Abs(kine)>kinecut) continue; |
ca28c5f5 |
176 | |
177 | Double_t mass=cascade->GetEffMass(); |
178 | cs->Fill(mass); |
04b2a5f1 |
179 | csf->Fill(mass); |
ca28c5f5 |
180 | |
181 | if (TMath::Abs(mass-cascadeMass)>cascadeWidth) continue; |
182 | |
183 | Int_t j; |
184 | for (j=0; j<ngood; j++) { |
185 | if (gc[j].code != cascade->GetPdgCode()) continue; |
186 | if (gc[j].nlab == nlab) |
187 | if (gc[j].plab == plab) |
188 | if (gc[j].blab == blab) break; |
189 | } |
190 | |
191 | Double_t px,py,pz; cascade->GetPxPyPz(px,py,pz); |
192 | Double_t pt=TMath::Sqrt(px*px+py*py); |
193 | |
194 | if (j==ngood) { |
195 | hfake->Fill(pt); |
196 | cerr<<"Fake cascade: ("<<nlab<<","<<plab<<","<<blab<<")\n"; |
197 | continue; |
198 | } |
04b2a5f1 |
199 | csf->Fill(mass,-1); |
ca28c5f5 |
200 | |
201 | pxg=gc[j].px; pyg=gc[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg); |
202 | Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px); |
203 | Double_t lamg=TMath::ATan2(gc[j].pz,ptg), lam=TMath::ATan2(pz,pt); |
204 | hp->Fill((phi - phig)*1000.); |
205 | hl->Fill((lam - lamg)*1000.); |
206 | hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.); |
207 | |
208 | Double_t x,y,z; cascade->GetXYZ(x,y,z); |
209 | hx->Fill((x-gc[j].x)*10); |
210 | hy->Fill((y-gc[j].y)*10); |
211 | hz->Fill((z-gc[j].z)*10); |
212 | |
213 | hfound->Fill(ptg); |
214 | gc[j].nlab=-1; |
215 | |
216 | } |
217 | for (i=0; i<ngood; i++) { |
218 | if (gc[i].code != code) continue; |
219 | pxg=gc[i].px; pyg=gc[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg); |
220 | hgood->Fill(ptg); |
221 | nlab=gc[i].nlab; plab=gc[i].plab; blab=gc[i].blab; |
222 | if (nlab < 0) continue; |
223 | cerr<<"Cascade ("<<nlab<<','<<plab<<","<<blab<<") has not been found !\n"; |
224 | } |
225 | |
226 | carray.Delete(); |
227 | |
228 | Stat_t ng=hgood->GetEntries(); |
229 | Stat_t nf=hfound->GetEntries(); |
230 | |
231 | cerr<<"Number of found cascades: "<<nentr<<" ("<<nf<<" in the peak)\n"; |
232 | cerr<<"Number of good cascades: "<<ng<<endl; |
233 | |
234 | if (ng!=0) |
235 | cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n"; |
236 | |
237 | gStyle->SetOptStat(111110); |
238 | gStyle->SetOptFit(1); |
239 | |
240 | TCanvas *c1=new TCanvas("c1","",0,0,580,610); |
241 | c1->Divide(2,2); |
242 | |
243 | c1->cd(1); |
244 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
245 | //hp->Fit("gaus"); |
246 | hp->Draw(); |
247 | hl->Draw("same"); c1->cd(); |
248 | |
249 | c1->cd(2); |
250 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
251 | //hpt->Fit("gaus"); c1->cd(); |
252 | hpt->Draw(); c1->cd(); |
253 | |
254 | c1->cd(3); |
255 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
256 | //hx->Fit("gaus"); |
257 | hx->Draw(); |
258 | hy->Draw("same"); c1->cd(); |
259 | |
260 | c1->cd(4); |
261 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
262 | //hz->Fit("gaus"); |
263 | hz->Draw(); |
264 | |
265 | |
266 | TCanvas *c2=new TCanvas("c2","",600,0,580,610); |
267 | c2->Divide(1,2); |
268 | |
269 | c2->cd(1); |
270 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
271 | hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2(); |
272 | hg->Divide(hfound,hgood,1,1.,"b"); |
273 | hf->Divide(hfake,hgood,1,1.,"b"); |
274 | hg->SetMaximum(1.4); |
275 | hg->SetYTitle("Cascade reconstruction efficiency"); |
276 | hg->SetXTitle("Pt (GeV/c)"); |
277 | hg->Draw(); |
278 | |
279 | TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4); |
280 | line1->Draw("same"); |
281 | TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4); |
282 | line2->Draw("same"); |
283 | |
284 | hf->SetFillColor(1); |
285 | hf->SetFillStyle(3013); |
286 | hf->SetLineColor(2); |
287 | hf->SetLineWidth(2); |
288 | hf->Draw("histsame"); |
289 | TText *text = new TText(0.461176,0.248448,"Fake cascades"); |
290 | text->SetTextSize(0.05); |
291 | text->Draw(); |
292 | text = new TText(0.453919,1.11408,"Good cascades"); |
293 | text->SetTextSize(0.05); |
294 | text->Draw(); |
295 | |
296 | |
297 | c2->cd(2); |
298 | gPad->SetFillColor(42); gPad->SetFrameFillColor(10); |
ca28c5f5 |
299 | //cs->Fit("gaus","","",cascadeMass-cascadeWidth,cascadeMass+cascadeWidth); |
300 | cs->Draw(); |
04b2a5f1 |
301 | csf->Draw("same"); |
ca28c5f5 |
302 | Double_t max=cs->GetMaximum(); |
303 | TLine *line3 = |
304 | new TLine(cascadeMass-cascadeWidth,0.,cascadeMass-cascadeWidth,max); |
305 | line3->Draw("same"); |
306 | TLine *line4 = |
307 | new TLine(cascadeMass+cascadeWidth,0.,cascadeMass+cascadeWidth,max); |
308 | line4->Draw("same"); |
309 | |
310 | timer.Stop(); timer.Print(); |
311 | |
312 | return 0; |
313 | } |
314 | |
315 | |
316 | |
317 | Int_t good_cascades(GoodCascade *gc, Int_t max) { |
318 | Int_t nc=0; |
319 | /*** Get information about the cuts ***/ |
320 | Double_t r2min=0.9*0.9; |
321 | Double_t r2max=2.9*2.9; |
322 | |
323 | /*** Get labels of the "good" tracks ***/ |
324 | Double_t dd; Int_t id, label[15000], ngt=0; |
325 | ifstream in("good_tracks_its"); |
326 | if (!in) { |
327 | cerr<<"Can't open the file good_tracks_its \n"; |
328 | return nc; |
329 | } |
330 | while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) { |
331 | ngt++; |
332 | if (ngt>=15000) { |
333 | cerr<<"Too many good ITS tracks !\n"; |
334 | return nc; |
335 | } |
336 | } |
337 | if (!in.eof()) { |
338 | cerr<<"Read error (good_tracks_its) !\n"; |
339 | return nc; |
340 | } |
341 | |
342 | /*** Get an access to the kinematics ***/ |
343 | if (gAlice) {delete gAlice; gAlice=0;} |
344 | |
345 | TFile *file=TFile::Open("galice.root"); |
346 | if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);} |
347 | if (!(gAlice=(AliRun*)file->Get("gAlice"))) { |
348 | cerr<<"gAlice has not been found on galice.root !\n"; |
349 | exit(5); |
350 | } |
351 | |
352 | Int_t np=gAlice->GetEvent(0); |
353 | while (np--) { |
354 | cerr<<np<<'\r'; |
355 | TParticle *cp=gAlice->Particle(np); |
356 | |
357 | /*** only these cascades are "good" ***/ |
358 | Int_t code=cp->GetPdgCode(); |
359 | if (code!=kXiMinus) if (code!=kXiPlusBar) |
360 | if (code!=kOmegaMinus) if (code!=kOmegaPlusBar) continue; |
361 | |
362 | /*** daughter tracks must be "good" ***/ |
363 | Int_t v0lab=cp->GetFirstDaughter(), blab=cp->GetLastDaughter(); |
364 | if (v0lab==blab) continue; |
365 | if (v0lab<0) continue; |
366 | if (blab<0) continue; |
367 | |
368 | TParticle *p0=gAlice->Particle(v0lab); |
369 | TParticle *bp=gAlice->Particle(blab); |
370 | if ((p0->GetPdgCode()!=kLambda0) && (p0->GetPdgCode()!=kLambda0Bar)) { |
371 | TParticle *p=p0; p0=bp; bp=p; |
372 | Int_t i=v0lab; v0lab=blab; blab=i; |
373 | if ((p0->GetPdgCode()!=kLambda0) && (p0->GetPdgCode()!=kLambda0Bar)) |
374 | continue; |
375 | } |
376 | |
377 | /** is the bachelor "good" ? **/ |
378 | Int_t i; |
379 | for (i=0; i<ngt; i++) if (label[i]==blab) break; |
380 | if (i==ngt) continue; |
381 | |
382 | /** is the V0 "good" ? **/ |
383 | Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter(); |
384 | if (nlab==plab) continue; |
385 | if (nlab<0) continue; |
386 | if (plab<0) continue; |
387 | |
388 | for (i=0; i<ngt; i++) if (label[i]==nlab) break; |
389 | if (i==ngt) continue; |
390 | for (i=0; i<ngt; i++) if (label[i]==plab) break; |
391 | if (i==ngt) continue; |
392 | |
393 | /*** fiducial volume ***/ |
394 | Double_t x=bp->Vx(), y=bp->Vy(), r2=x*x+y*y; //bachelor |
395 | if (r2<r2min) continue; |
396 | if (r2>r2max) continue; |
397 | TParticle *pp=gAlice->Particle(plab); |
398 | x=pp->Vx(); y=pp->Vy(); r2=x*x+y*y; //V0 |
399 | if (r2<r2min) continue; |
400 | if (r2>r2max) continue; |
401 | |
402 | if (gAlice->Particle(plab)->GetPDG()->Charge() < 0.) { |
403 | i=plab; plab=nlab; nlab=i; |
404 | } |
405 | |
406 | gc[nc].code=code; |
407 | gc[nc].plab=plab; gc[nc].nlab=nlab; gc[nc].blab=blab; |
408 | gc[nc].px=cp->Px(); gc[nc].py=cp->Py(); gc[nc].pz=cp->Pz(); |
409 | gc[nc].x=bp->Vx(); gc[nc].y=bp->Vy(); gc[nc].z=bp->Vz(); |
410 | nc++; |
411 | |
412 | } |
413 | |
414 | |
415 | return nc; |
416 | } |