Updated code for tracking V2 (from Y. Belikov)
[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 #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 "AliV0vertex.h"
28 #endif
29
30 struct GoodVertex {
31   Int_t nlab,plab;
32   Int_t code;
33   Float_t px,py,pz;
34   Float_t x,y,z;
35 };
36 Int_t good_vertices(GoodVertex *gt, Int_t max);
37
38 Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
39    cerr<<"Doing comparison...\n";
40
41    const Double_t V0window=0.05, V0width=0.015; 
42    Double_t V0mass=0.5;
43    switch (code) {
44    case kK0Short:    V0mass=0.497672; break;
45    case kLambda0:    V0mass=1.115683; break;
46    case kLambda0Bar: V0mass=1.115683; break;
47    default: cerr<<"Invalid PDG code !\n"; return 1;
48    }
49
50    TStopwatch timer;
51
52    /*** Load reconstructed vertices ***/
53    TFile *vf=TFile::Open("AliV0vertices.root");
54    if (!vf->IsOpen()) {cerr<<"Can't open AliV0vertices.root !\n"; return 2;}
55    TObjArray varray(1000);
56    TTree *vTree=(TTree*)vf->Get("TreeV");
57    TBranch *branch=vTree->GetBranch("vertices");
58    Int_t nentr=(Int_t)vTree->GetEntries();
59    for (Int_t i=0; i<nentr; i++) {
60        AliV0vertex *iovertex=new AliV0vertex; branch->SetAddress(&iovertex);
61        vTree->GetEvent(i);
62        varray.AddLast(iovertex);
63    }
64
65    /*** Check if the file with the "good" vertices exists ***/
66    GoodVertex gv[1000];
67    Int_t ngood=0;
68    ifstream in("good_vertices");
69    if (in) {
70       cerr<<"Reading good vertices...\n";
71       while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
72                  gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
73                  gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
74          ngood++;
75          cerr<<ngood<<'\r';
76          if (ngood==1000) {
77             cerr<<"Too many good vertices !\n";
78             break;
79          }
80       }
81       if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
82    } else {
83      /*** generate a file with the "good" vertices ***/
84       cerr<<"Marking good vertices (this will take a while)...\n";
85       ngood=good_vertices(gv,1000);
86       ofstream out("good_vertices");
87       if (out) {
88          for (Int_t ngd=0; ngd<ngood; ngd++)            
89             out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
90                  gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
91                  gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
92       } else cerr<<"Can not open file (good_vertices) !\n";
93       out.close();
94    }
95
96    vf->Close();
97
98
99    TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution 
100    hp->SetXTitle("(mrad)"); hp->SetFillColor(2);
101    TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30);
102    hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013); 
103    TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.); 
104    hpt->SetXTitle("(%)"); hpt->SetFillColor(2); 
105
106    TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res. 
107    hx->SetXTitle("(mm)"); hx->SetFillColor(6);
108    TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.);   //y res
109    hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013);
110    TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.);   //z res. 
111    hz->SetXTitle("(mm)"); hz->SetFillColor(6);
112
113
114    Double_t pmin=0.2, pmax=4.2; Int_t nchan=20;
115    TH1F *hgood=new TH1F("hgood","Good Vertices",nchan,pmin,pmax);    
116    TH1F *hfound=new TH1F("hfound","Found Vertices",nchan,pmin,pmax);
117    TH1F *hfake=new TH1F("hfake","Fake Vertices",nchan,pmin,pmax);
118    TH1F *hg=new TH1F("hg","Efficiency for Good Vertices",nchan,pmin,pmax);
119    hg->SetLineColor(4); hg->SetLineWidth(2);
120    TH1F *hf=new TH1F("hf","Probability of Fake Vertices",nchan,pmin,pmax);
121    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
122
123    Double_t mmin=V0mass-V0window, mmax=V0mass+V0window;
124    TH1F *v0s =new TH1F("v0s","V0s Effective Mass",40, mmin, mmax);
125    v0s->SetXTitle("(GeV)"); v0s->SetFillColor(6);
126
127
128    Double_t pxg=0.,pyg=0.,ptg=0.;
129    Int_t nlab=-1, plab=-1;
130    Int_t i;
131    for (i=0; i<nentr; i++) {
132        AliV0vertex *vertex=(AliV0vertex*)varray.UncheckedAt(i);
133        nlab=TMath::Abs(vertex->GetNlabel());
134        plab=TMath::Abs(vertex->GetPlabel());
135
136        vertex->ChangeMassHypothesis(code);
137
138        Double_t mass=vertex->GetEffMass();
139        v0s->Fill(mass);
140
141        Int_t j;
142        for (j=0; j<ngood; j++) {
143           if (gv[j].code != vertex->GetPdgCode()) continue;
144           if (gv[j].nlab == nlab)
145           if (gv[j].plab == plab) break;
146        }
147
148        if (TMath::Abs(mass-V0mass)>V0width) continue;
149
150        Double_t px,py,pz; vertex->GetPxPyPz(px,py,pz);
151        Double_t pt=TMath::Sqrt(px*px+py*py);
152
153        if (j==ngood) {
154           hfake->Fill(pt);
155           cerr<<"Fake vertex: ("<<nlab<<","<<plab<<")\n";
156           continue;
157        }
158
159        pxg=gv[j].px; pyg=gv[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
160        Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px);
161        Double_t lamg=TMath::ATan2(gv[j].pz,ptg), lam=TMath::ATan2(pz,pt);
162        hp->Fill((phi - phig)*1000.);
163        hl->Fill((lam - lamg)*1000.);
164        hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.);
165
166        Double_t x,y,z; vertex->GetXYZ(x,y,z);
167        hx->Fill((x-gv[j].x)*10);
168        hy->Fill((y-gv[j].y)*10);
169        hz->Fill((z-gv[j].z)*10);
170
171        hfound->Fill(ptg);
172        gv[j].nlab=-1;
173
174    }
175    for (i=0; i<ngood; i++) {
176       if (gv[i].code != code) continue;
177       pxg=gv[i].px; pyg=gv[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
178       hgood->Fill(ptg);
179       nlab=gv[i].nlab; plab=gv[i].plab;
180       if (nlab < 0) continue;
181       cerr<<"Vertex ("<<nlab<<','<<plab<<") has not been found !\n";
182    }
183
184    varray.Delete();
185
186    Stat_t ng=hgood->GetEntries();
187    Stat_t nf=hfound->GetEntries();
188
189    cerr<<"Number of found vertices: "<<nentr<<" ("<<nf<<" in the peak)\n";
190    cerr<<"Number of good vertices: "<<ng<<endl;
191
192    if (ng!=0) 
193       cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
194
195    gStyle->SetOptStat(111110);
196    gStyle->SetOptFit(1);
197
198    TCanvas *c1=new TCanvas("c1","",0,0,580,610);
199    c1->Divide(2,2);
200
201    c1->cd(1);
202    gPad->SetFillColor(42); gPad->SetFrameFillColor(10); 
203    //hp->Fit("gaus");
204    hp->Draw();
205    hl->Draw("same"); c1->cd();
206
207    c1->cd(2);
208    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
209    //hpt->Fit("gaus"); c1->cd();
210    hpt->Draw(); c1->cd();
211
212    c1->cd(3);
213    gPad->SetFillColor(42); gPad->SetFrameFillColor(10); 
214    //hx->Fit("gaus"); 
215    hx->Draw(); 
216    hy->Draw("same"); c1->cd();
217
218    c1->cd(4);
219    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
220    //hz->Fit("gaus");
221    hz->Draw();
222
223
224    TCanvas *c2=new TCanvas("c2","",600,0,580,610);
225    c2->Divide(1,2);
226
227    c2->cd(1);
228    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
229    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
230    hg->Divide(hfound,hgood,1,1.,"b");
231    hf->Divide(hfake,hgood,1,1.,"b");
232    hg->SetMaximum(1.4);
233    hg->SetYTitle("V0 reconstruction efficiency");
234    hg->SetXTitle("Pt (GeV/c)");
235    hg->Draw();
236
237    TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
238    line1->Draw("same");
239    TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
240    line2->Draw("same");
241
242    hf->SetFillColor(1);
243    hf->SetFillStyle(3013);
244    hf->SetLineColor(2);
245    hf->SetLineWidth(2);
246    hf->Draw("histsame");
247    TText *text = new TText(0.461176,0.248448,"Fake vertices");
248    text->SetTextSize(0.05);
249    text->Draw();
250    text = new TText(0.453919,1.11408,"Good vertices");
251    text->SetTextSize(0.05);
252    text->Draw();
253
254
255    c2->cd(2);
256    gPad->SetFillColor(42); gPad->SetFrameFillColor(10);
257    v0s->SetXTitle("(GeV/c)"); 
258    v0s->Fit("gaus","","",V0mass-V0width,V0mass+V0width);
259    Double_t max=v0s->GetMaximum();
260    TLine *line3 = new TLine(V0mass-V0width,0.,V0mass-V0width,max);
261    line3->Draw("same");
262    TLine *line4 = new TLine(V0mass+V0width,0.,V0mass+V0width,max);
263    line4->Draw("same");
264
265    timer.Stop(); timer.Print();
266
267    return 0;
268 }
269
270 Int_t good_vertices(GoodVertex *gv, Int_t max) {
271    Int_t nv=0;
272    /*** Get information about the cuts ***/
273    Double_t r2min=0.9*0.9;
274    Double_t r2max=2.9*2.9;
275
276    /*** Get labels of the "good" tracks ***/
277    Double_t dd; Int_t id, label[15000], ngt=0;
278    ifstream in("good_tracks_its");
279    if (!in) {
280      cerr<<"Can't open the file good_tracks_its \n";
281      return nv;
282    }
283    while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) {
284      ngt++;
285      if (ngt>=15000) {
286        cerr<<"Too many good ITS tracks !\n";
287        return nv;
288      }
289    }   
290    if (!in.eof()) {
291       cerr<<"Read error (good_tracks_its) !\n";
292       return nv;
293    }
294
295    /*** Get an access to the kinematics ***/
296    if (gAlice) {delete gAlice; gAlice=0;}
297
298    TFile *file=TFile::Open("galice.root");
299    if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
300    if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
301      cerr<<"gAlice has not been found on galice.root !\n";
302      exit(5);
303    }
304
305    Int_t np=gAlice->GetEvent(0);
306    while (np--) {
307       cerr<<np<<'\r';
308       TParticle *p0=gAlice->Particle(np);
309
310       /*** only these V0s are "good" ***/
311       Int_t code=p0->GetPdgCode();
312       if (code!=kK0Short)
313       if (code!=kLambda0)
314       if (code!=kLambda0Bar) continue; 
315
316       /*** daughter tracks should be "good" ***/
317       Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
318       if (nlab==plab) continue;
319       if (nlab<0) continue;
320       if (plab<0) continue;
321       Int_t i;
322       for (i=0; i<ngt; i++) if (label[i]==nlab) break;
323       if (i==ngt) continue;
324       for (i=0; i<ngt; i++) if (label[i]==plab) break;
325       if (i==ngt) continue;
326
327       /*** fiducial volume ***/
328       TParticle *p=gAlice->Particle(nlab);
329       Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y;
330       if (r2<r2min) continue;
331       if (r2>r2max) continue;
332        
333       if (gAlice->Particle(plab)->GetPDG()->Charge() < 0.) {
334          i=plab; plab=nlab; nlab=i;
335       }
336       
337       gv[nv].code=code;
338       gv[nv].plab=plab; gv[nv].nlab=nlab;
339       gv[nv].px=p0->Px(); gv[nv].py=p0->Py(); gv[nv].pz=p0->Pz();
340       gv[nv].x=x; gv[nv].y=y; gv[nv].z=z;
341       nv++;
342    }
343  
344    delete gAlice; gAlice=0;
345
346    file->Close();  
347
348    return nv;
349 }