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