]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/oldmacros/dedxanal.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated selection in ReadFromTracks()
[u/mrichter/AliRoot.git] / ITS / oldmacros / dedxanal.C
CommitLineData
08cb00f2 1//=========================================
2int HMax=80,PMax=6.;
3int NStat=15000;
4Float_t pmin,pmax;
5
6char tit[]="--------------------------";
7TH1F pions=TH1F("Qpi",tit,100,0.5,PMax);
8TH1F *qpi =&pions;
9
10TH2F qplotP, qplotKa, qplotPi, qplotE;
11TH1F *q1plot;
12TH2F *qplot;
13void dedxhis(Int_t pcod=0,Float_t pmin=0,Float_t pmax=1300)
14{
15 if(!q1plot)q1plot = new TH1F("Qhis","Particle charge",100,0.5,PMax);
16 if(!qplot){ qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
17 // TH2F qplotP, qplotKa, qplotPi, qplotE;
18 qplot->Copy(qplotP);
19 qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlack); qplotP.SetMarkerSize(.2);
20 qplot->Copy(qplotKa);
21 qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.2);
22 qplot->Copy(qplotPi);
23 qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); qplotPi.SetMarkerSize(.2);
24 qplot->Copy(qplotE);
25 qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.2);
26 }else{
27 qplotP.Reset();
28 qplotPi.Reset(); qplotKa.Reset(); qplotE.Reset();
29 }
30// -------------------------------------------------------------
31 Float_t Qtr,Pmom;
32 for(Int_t i=0;i<NStat;i++){
33 TVector tabv(*( pid->GetVec(i) ));
34 Qtr=tabv(6);Pmom=tabv(10);
35 if(pcod==0)
36 { if(tabv(0)>=1 ){
37 //qplot->Fill( Pmom,Qtr );
38 if(tabv(11)==2212)qplotP.Fill( Pmom,Qtr );
39 if(tabv(11)== 321)qplotKa.Fill( Pmom,Qtr );
40 if(tabv(11)== 211)qplotPi.Fill( Pmom,Qtr );
41 if(tabv(11)== 11)qplotE.Fill( Pmom,Qtr );
42 }
43 }
44 else { if(tabv(11)==pcod && tabv(0)>=1 )qplot->Fill( Pmom,Qtr );}
45
46 //if( tabv(0)>=2 && tabv(11)==pcod )qplot->Fill( Pmom,Qtr );
47 //if(tabv(0)>=1)qplot->Fill( Pmom,Qtr );
48 if(Pmom>pmin&&Pmom<pmax)q1plot->Fill(Qtr);
49 }
50// ------------------------------------------------------------
51// c1=new TCanvas("c1"," ",10,10,700,500);
52// pad1 =new TPad("p1"," " ,0 , 0,0,1,21);
53// pad1 =new TPad("p1"," " ,0 , 0,.3,1,21);
54// pad2 =new TPad("p2"," " ,.35, 0, 1,1,0);
55// pad2->Draw(); pad2->cd();
56//.....................
57 //if(pcod==0){ qplotP.Draw(); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same"); }
58
59 if(pcod==0){
60
61 //qplotP->GetXaxis()->SetTitleSize(0.05);
62 //qplotP->GetYaxis()->SetTitleSize(0.05);
63 //gPad->SetFillColor(kWhite); // b.b.
64 gStyle->SetOptStat(0);
65 qplotP->SetXTitle("(Mev/c)");
66 qplotP->SetYTitle("(mips)");
67 qplotP.Draw(); qplotKa.Draw("same"); qplotPi.Draw("same");
68 qplotE.Draw("same");
69 TText *text = new TText(800.,11.,"Protons");
70 text->SetTextSize(0.05);
71 text->SetTextColor(kBlack);
72 text->Draw();
73 TText *text = new TText(800.,10.,"Kaons");
74 text->SetTextSize(0.05);
75 text->SetTextColor(kRed);
76 text->Draw();
77 TText *text = new TText(800.,9.,"Pions");
78 text->SetTextSize(0.05);
79 text->SetTextColor(kBlue);
80 text->Draw();
81 TText *text = new TText(800.,8.,"Electrons");
82 text->SetTextSize(0.05);
83 text->SetTextColor(kGreen);
84 text->Draw();
85
86 }
87 else{
88 qplot->Draw();}
89 c1->Range(0,0,1300,10);
90 gStyle->SetLineColor(kRed);
91 gStyle->SetLineWidth(2);
92 TLine *lj[3],*lk[3];
93 for(Int_t j=0;j<3;j++){
94 Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
95 x1=pid->cut[j+1][0]; x2=pid->cut[j+2][0];
96 y1=y2=pid->cut[j+2][2];
97 lj[j]=new TLine(1000*x1,y1,1000*x2,y2);
98 //lj[j]->Draw();
99 if(j==0){yy1=10.;}else{yy1=lj[j-1]->GetY1();}
100 yy2=lj[j]->GetY1();
101 xx1=xx2=x1;
102 lk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2);
103 //lk[j]->Draw();
104 }
105 //Draw pions-kaons cuts.
106//..................................
107 TLine *mj[7],*mk[7];
108 for(Int_t j=0;j<7;j++){
109 Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
110 x1=pid->cut[j+2][0]; x2=pid->cut[j+3][0];
111 y1=y2=pid->cut[j+3][5];
112 mj[j]=new TLine(1000*x1,y1,1000*x2,y2);
113 //mj[j]->Draw();
114 if(j==0){yy1=10.;}else{yy1=mj[j-1]->GetY1();}
115 yy2=mj[j]->GetY1();
116 xx1=xx2=x1;
117 mk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2);
118 //mk[j]->Draw();
119 }
120 //Draw kaons-protons cuts.
121 gStyle->SetLineWidth(1.);
122}
123
124void qhispi(Float_t pmin=0,Float_t pmax=1300)
125{
126 char tit[]="--------------------------";
127 sprintf(tit,"%s%.1f %.1f","Pmom ",pmin,pmax);
128// TH1F *qpi = new TH1F("Qpi",tit,100,0.5,PMax);
129// TH1F *nhit = new TH1F("Nhit",tit,100,0,6);
130 TH2F *nhit = new TH2F("Nhit",tit,100,0,6,100,0,1300);
131// -------------------------------------------------------------
132 Float_t Qtr,Pmom;
133 for(Int_t i=0;i<NStat;i++){
134 TVector tabv(*( pid->GetVec(i) ));
135 Qtr=tabv(6);Pmom=tabv(10);
136 if(Pmom>pmin&&Pmom<pmax&&abs(tabv(11))==211)qpi->Fill(Qtr);
137// if(tabv(11)==2212)nhit->Fill(tabv(0));
138 if(tabv(11)==211)nhit->Fill(tabv(0),Pmom);
139 }
140// ------------------------------------------------------------
141 int pistat=qpi->GetEntries();
142 //qpi->SetMaximum(pistat+1);
143 qpi->SetMaximum(HMax);
144 qpi->SetFillColor(42);
145 qpi->Fit("gaus");
146 TF1 *fun=qpi->GetFunction("gaus");
147 fun->SetLineWidth(.1);
148 qpi->SetLineWidth(.1);
149 qpi->Draw();
150//---------------------------
151 if(0){
152 gaus1=new TF1("g1","gaus",.5,2);
153 qpi->SetMaximum(-1111);
154 gaus1->SetParameter(0, qpi->GetMaximum() );
155 qpi->SetMaximum(HMax);
156 gaus1->SetParameter(1,1);
157 gaus1->SetParameter(2,.12);
158 gaus1->SetLineWidth(.1);
159 gaus1->Draw("same");
160 }
161}
162
163void qhiska(Float_t pmin=0,Float_t pmax=1300)
164{
165 TH1F *qka = new TH1F("Qka","Particle charge",100,0.5,PMax);
166// -------------------------------------------------------------
167 Float_t Qtr,Pmom;
168 for(Int_t i=0;i<NStat;i++){
169 TVector tabv(*( pid->GetVec(i) ));
170 Qtr=tabv(6);Pmom=tabv(10);
171// qplot->Fill( Pmom,Qtr );
172 if(Pmom>pmin&&Pmom<pmax&&tabv(11)==321)qka->Fill(Qtr);
173 }
174// ------------------------------------------------------------
175 qka->SetFillColor(0);
176 qka->SetLineColor(kRed);
177//...................
178 qka->Fit("gaus","0");
179 TF1 *fun=qka->GetFunction("gaus");
180 fun->SetLineWidth(.1);
181 qka->SetLineWidth(.1);
182 qka->Draw("same");
183//...................
184// qka->Draw("same");
185// gaus_k(pmin,pmax,qka);
186
187}
188
189void qhisp(Float_t pmin=0,Float_t pmax=1300)
190{
191 TH1F *qpr = new TH1F("Qpr","Particle charge",100,0.5,PMax);
192// -------------------------------------------------------------
193 qpr->Clear();
194 Float_t Qtr,Pmom;
195 for(Int_t i=0;i<NStat;i++){
196 TVector tabv(*( pid->GetVec(i) ));
197 Qtr=tabv(6);Pmom=tabv(10);
198// qplot->Fill( Pmom,Qtr );
199 if(Pmom>pmin&&Pmom<pmax&&tabv(11)==2212)qpr->Fill(Qtr);
200 }
201// ------------------------------------------------------------
202 qpr->SetFillColor(16);
203//...................
204 qpr->Fit("gaus","0");
205 TF1 *fun=qpr->GetFunction("gaus");
206 fun->SetLineWidth(.1);
207 qpr->SetLineWidth(.1);
208 qpr->Draw("same");
209//...................
210// qpr->Draw("same");
211// gaus_p(pmin,pmax,qpr);
212}
213//---------------------------------------
214
215void gaus_k(Float_t pmin=450,Float_t pmax=470,TH1F *qka){
216 Float_t xmean,xsig;
217 if(pmin==410&&pmax==470){xmean=pid->cut[5][3]; xsig=pid->cut[5][4];}
218 else
219 if(pmin==470&&pmax==530){xmean=pid->cut[6][3]; xsig=pid->cut[6][4];}
220 else
221 if(pmin==730&&pmax==830){xmean=pid->cut[10][3]; xsig=pid->cut[10][4];}
222 else
223 if(pmin==830&&pmax==930){xmean=pid->cut[11][3]; xsig=pid->cut[11][4];}
224 else
225 if(pmin==930&&pmax==1030){xmean=pid->cut[12][3]; xsig=pid->cut[12][4];}
226 else{return;}
227 gaus1=new TF1("g1","gaus",xmean-3*xsig,xmean+3*xsig);
228 qka->SetMaximum(-1111);
229 gaus1->SetParameter(0, qka->GetMaximum() );
230 qka->SetMaximum(HMax);
231 gaus1->SetParameter(1,xmean);
232 gaus1->SetParameter(2,xsig);
233 gaus1->SetLineWidth(.1);
234 gaus1->Draw("same");
235}
236//---------------------------------------
237void gaus_p(Float_t pmin=730,Float_t pmax=830,TH1F *qp){
238 Float_t xmean,xsig;
239 if(pmin==730&&pmax==830){xmean=pid->cut[10][5]; xsig=pid->cut[10][6];}
240 else
241 if(pmin==830&&pmax==930){xmean=pid->cut[11][5]; xsig=pid->cut[11][6];}
242 else
243 if(pmin==930&&pmax==1030){xmean=pid->cut[12][5]; xsig=pid->cut[12][6];}
244 else{return;}
245 gaus1=new TF1("g1","gaus",xmean-3*xsig,xmean+3*xsig);
246 qp->SetMaximum(-1111);
247 gaus1->SetParameter(0, qp->GetMaximum() );
248 qp->SetMaximum(HMax);
249 gaus1->SetParameter(1,xmean);
250 gaus1->SetParameter(2,xsig);
251 gaus1->SetLineWidth(.1);
252 gaus1->Draw("same");
253}
254
255//----b.b.---------------------------------
256void fitpi(Float_t pmin=0,Float_t pmax=1300,char *tit)
257{
258 cout<<"fitpi: NStat, PMax, pmin, pmax ="<<NStat<<","<<PMax<<","<<pmin<<","<<pmax<<endl;
259 TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
260// -------------------------------------------------------------
261 Float_t Qtr,Pmom;
262 for(Int_t i=0;i<NStat;i++){
263 TVector tabv(*( pid->GetVec(i) ));
264 Qtr=tabv(6);Pmom=tabv(10);
265 if(Pmom>pmin&&Pmom<pmax&&tabv(11)==211)q1fit->Fill(Qtr);
266 }
267 q1fit->SetFillColor(0);
268 q1fit->SetLineColor(kRed);
269 q1fit->Draw();
270 q1fit->Fit("gaus");
271}
272
273//---------------------------------------
274void fitka(Float_t pmin=0,Float_t pmax=1300,char *tit)
275{
276 //TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
277 TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,10.);
278// -------------------------------------------------------------
279 Float_t Qtr,Pmom;
280 for(Int_t i=0;i<NStat;i++){
281 TVector tabv(*( pid->GetVec(i) ));
282 Qtr=tabv(6);Pmom=tabv(10);
283 if(Pmom>pmin&&Pmom<pmax&&tabv(11)==321)q1fit->Fill(Qtr);
284 }
285 q1fit->SetFillColor(0);
286 q1fit->SetLineColor(kRed);
287 q1fit->Draw();
288 q1fit->Fit("gaus");
289}
290//---------------------------------------
291void fitp(Float_t pmin=0,Float_t pmax=1300,char *tit)
292{
293 TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
294// -------------------------------------------------------------
295 Float_t Qtr,Pmom;
296 for(Int_t i=0;i<NStat;i++){
297 TVector tabv(*( pid->GetVec(i) ));
298 Qtr=tabv(6);Pmom=tabv(10);
299 if(Pmom>pmin&&Pmom<pmax&&tabv(11)==2212)q1fit->Fill(Qtr);
300 }
301 q1fit->SetFillColor(0);
302 q1fit->Draw();
303 q1fit->Fit("gaus");
304}
305//---------------------------------------
306void effall(){
307
308 eff(211,211,"",3);
309 eff(321,321,"same",2);
310 eff(2212,2212,"same",1);
311}
312//---------------------------------------
313void eff(int pc=211,int pc2=211,char *opt="",int color=2)
314{
315 const Int_t HSize=100;
316 TH1F *efpi = new TH1F("Effpi","Eff of PID",HSize,0.025,1400);
317 TH1F *pmom = new TH1F("Pmom" ,"Eff of PID",HSize,0.025,1400);
318 TH1F *heff = new TH1F("Eff%" ,"Eff of PID",HSize,0.025,1400);
319// -------------------------------------------------------------
320 Float_t Qtr,Pmom;
321 Int_t Pcode;
322 for(Int_t i=0;i<NStat;i++){
323 TVector tabv(*( pid->GetVec(i) ));
324 Qtr=tabv(6);Pmom=tabv(10);
325 Pcode= ( pid->GetPcode(Qtr,Pmom/1000.) );
326 if(tabv(0)>=4){
327 if(Pcode==pc&&tabv(11)==pc2) efpi->Fill(Pmom);
328 if(tabv(11)==pc2) pmom->Fill(Pmom);
329 }
330 }
331// ------------------------------------------------------------
332 for(int i=1;i<=HSize;i++){
333 if( pmom->fArray[i] > 0 )
334 heff->fArray[i] = color * ( efpi->fArray[i] )/( pmom->fArray[i] );
335 }
336 heff->SetLineColor(color);
337 if(color==1)heff->SetLineWidth(2)else heff->SetLineWidth(.5);
338 heff->Draw(opt);
339return;
340 pmom->SetFillColor(16); pmom->Draw();
341 efpi->SetFillColor(0); //42
342 efpi->Draw("same");
343}
344//---------------------------------------
345Int_t NRange=0;
346 Float_t xpmin[]={410,470,730,830,930};
347 Float_t xpmax[]={470,530,830,930,1030};
348void all(){
349 PMax=3; HMax=500; NStat=99000;
350 Float_t pmin,pmax;
351 pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
352 qhispi(pmin,pmax); qhisp(pmin,pmax); qhiska(pmin,pmax);
353}
354//---------------------------------------
355void fitkall(){
356 PMax=3;
357 Float_t pmin,pmax;
358 //pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
359 pmin=200;pmax=300; // b.b.
360 char str[]=" ";
361 sprintf(str,"Kaons %d - %d MeV/c",pmin,pmax);
362 gStyle->SetOptFit();
363 fitka(pmin,pmax,str);
364}
365
366//--b.b.-------------------------------------
367void fitpiall(){
368 PMax=3;
369 NRange=3; // b.b.
370 Float_t pmin,pmax;
371 cout<<"fitpiall: NRange ="<<NRange<<endl;
372 //pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
373 pmin=100;pmax=200; // b.b.
374 char str[]=" ";
375 sprintf(str,"Pions %d - %d MeV/c",pmin,pmax);
376 gStyle->SetOptFit();
377 cout<<"fitpiall: pmin, pmax ="<<pmin<<","<<pmax<<endl;
378 fitpi(pmin,pmax,str);
379}
380
381//---------------------------------------
382void fitpall(){
383 PMax=3;
384 Float_t pmin,pmax;
385 if(NRange==0)NRange=2;
386 pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=2;
387 char str[]=" ";
388 sprintf(str,"Protons %d - %d MeV/c",pmin,pmax);
389 gStyle->SetOptFit();
390 fitp(pmin,pmax,str);
391}
392//-------------
393void newcuts(){
394//
395
396 pid->SetCut(3,.3, 0 , 2.5, 2.5 , 9, 9, 10 ); //200-300
397
398 pid->SetCut(5,.47, 1 , 0.12, 1.98 , 1.17 , 2.5, 10 );//410-470
399 pid->SetCut(6,.53, 1 , 0.12, 1.73 , 0.15 , 2.5, 10 );//470-530
400 pid->SetCut(7,.59, 0 , 0, 1.18 , 1.125 , 2.5, 10 );//530-590
401
402 pid->SetCut(8,.65, 0 , 0, 1.18, 1.125 , 2.3, 10 );//590-650
403}
404//---------------------------------------
405void qhisall(Float_t pmin=0.25,Float_t pmax=.700)
406//void qhisall()
407{
408 qhispi(pmin,pmax);
409 qhisp(pmin,pmax);
410 qhiska(pmin,pmax);
411}
412//--------------------------------------
413void pcode(){
414 TH1F *his = new TH1F("pcode","tit",100,10,2300);
415Float_t Pcod;
416for(Int_t i=0;i<NStat;i++){
417 TVector tabv(*( pid->GetVec(i) ));
418 Pcod=tabv(11); if(Pcod>0) his->Fill(Pcod);
419}
420his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
421}
422
423void signal(){
424 TH1F *hisR = new TH1F("Rec signal","tit",100,0,13);
425 TH1F *hisH = new TH1F("Hit signal","tit",100,0,13);
426 Float_t xx;
427 for(Int_t i=0;i<NStat;i++){
428 TVector tabv(*( pid->GetVec(i) ));
429 xx=tabv(1); if(tabv(0)>0) hisR->Fill(xx);
430 }
431 hisR->SetFillColor(0);hisR->SetLineColor(kRed);hisR->Draw();
432 hisH->SetFillColor(0);hisH->SetLineColor(kBlue);hisH->Draw("same");
433}
434//-------------------------------------
435
436void pmom(){
437 TH1F *his = new TH1F("pmom","tit",100,0,1000);
438 Float_t xx;
439 for(Int_t i=0;i<NStat;i++){
440 TVector tabv(*( pid->GetVec(i) ));
441 xx=tabv(10); if(tabv(0)>0) his->Fill(xx);
442 }
443 his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
444}
445
446// Print next track
447int jtrack=0;
448void track(){
449 if(jtrack==NStat)jtrack=0;
450 for(Int_t i=jtrack;i<NStat;i++){
451 TVector tabv(*( pid->GetVec(i) ));
452 if(tabv(0)>0)
453 { jtrack=i;break; }
454 }
455cout<<"Track No "<<jtrack<<endl;
456pid->Print(jtrack++);
457}
458
459// Fill histogram for track number
460//--------------------------------
461void tracks(){
462 TH1F *his = new TH1F("tracks","tit",100,0,15000);
463 Float_t xx;
464 for(Int_t i=0;i<NStat;i++){
465 TVector tabv(*( pid->GetVec(i) ));
466 if(tabv(0)>0)
467 { his->Fill(i); }
468 }
469 his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
470}
471// Fill pid table with reconstructed tracks
472
473
474#ifndef __CINT__
475 #include <iostream.h>
476 #include <fstream.h>
477
478 #include "AliRun.h"
479 #include "AliITS.h"
480 #include "AliITSgeom.h"
481 #include "AliITStrackerV2.h"
482 #include "AliITStrackV2.h"
483 #include "AliITSclusterV2.h"
484
485 #include "TFile.h"
486 #include "TTree.h"
487 #include "TH1.h"
488 #include "TObjArray.h"
489 #include "TStyle.h"
490 #include "TCanvas.h"
491 #include "TLine.h"
492 #include "TText.h"
493#endif
494
495struct GoodTrackITS {
496 Int_t lab;
497 Int_t code;
498 Float_t px,py,pz;
499 Float_t x,y,z;
500};
501
502
503
504void filltab_tracks(){
505
506 cerr<<"Filling of track table...\n";
507
508 Int_t event=0;
509Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event);
510
511 const Int_t MAX=15000;
512 Int_t nentr=0; TObjArray tarray(2000);
513 {/* Load tracks */
514 TFile *tf=TFile::Open("AliITStracksV2.root");
515 if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
516 char tname[100]; sprintf(tname,"TreeT_ITS_%d",event);
517 TTree *tracktree=(TTree*)tf->Get(tname);
518 if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
519 TBranch *tbranch=tracktree->GetBranch("tracks");
520 nentr=(Int_t)tracktree->GetEntries();
521 for (Int_t i=0; i<nentr; i++) {
522 AliITStrackV2 *iotrack=new AliITStrackV2;
523 tbranch->SetAddress(&iotrack);
524 tracktree->GetEvent(i);
525 tarray.AddLast(iotrack);
526 }
527 delete tracktree; //Thanks to Mariana Bondila
528 tf->Close();
529 }
530
531 /* Generate a list of "good" tracks */
532 GoodTrackITS gt[MAX];
533 Int_t ngood=0;
534 Float_t ConvRadAng=180./TMath::Pi();
535 ifstream in("good_tracks_its");
536 if (in) {
537 cerr<<"Reading good tracks...\n";
538 while (in>>gt[ngood].lab>>gt[ngood].code>>
539 gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
540 gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
541 ngood++;
542 cerr<<ngood<<'\r';
543 if (ngood==MAX) {
544 cerr<<"Too many good tracks !\n";
545 break;
546 }
547 }
548 if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
549 } else {
550 cerr<<"Marking good tracks (this will take a while)...\n";
551 ngood=good_tracks_its(gt,MAX,event);
552 ofstream out("good_tracks_its");
553 if (out) {
554 for (Int_t ngd=0; ngd<ngood; ngd++)
555 out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
556 <<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
557 <<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
558 } else cerr<<"Can not open file (good_tracks_its) !\n";
559 out.close();
560 }
561 cerr<<"Number of good tracks : "<<ngood<<endl;
562 TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
563 TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
564 TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
565 hpt->SetFillColor(2);
566 TH1F *hmpt=new TH1F("hmpt","Transverse impact parameter",30,-300,300);
567 hmpt->SetFillColor(6);
568 TH1F *hz=new TH1F("hz","Longitudinal impact parameter",30,-300,300);
569 //hmpt->SetFillColor(6);
570
571 AliITStrackV2 *trk=(AliITStrackV2*)tarray.UncheckedAt(0);
572 Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
573 Double_t pmax=6.0+pmin;
574
575 TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
576 TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
577 TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
578 TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
579 hg->SetLineColor(4); hg->SetLineWidth(2);
580 TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax);
581 hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
582
583 TH1F *hptw=new TH1F("hptw","Weghted pt",30,pmax,pmin);
584
585 while (ngood--) {
586 Int_t lab=gt[ngood].lab, tlab=-1;
587 Double_t pxg=gt[ngood].px, pyg=gt[ngood].py, pzg=gt[ngood].pz;
588 Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
589 Double_t pg=1000.*TMath::Sqrt(pxg*pxg+pyg*pyg+pzg*pzg); // b.b.
590
591
592 if (ptg<pmin) continue;
593
594 hgood->Fill(ptg);
595
596 AliITStrackV2 *track=0;
597 Int_t j;
598 for (j=0; j<nentr; j++) {
599 track=(AliITStrackV2*)tarray.UncheckedAt(j);
600 tlab=track->GetLabel();
601 if (lab==TMath::Abs(tlab)) break;
602 }
603 if (j==nentr) {
604 cerr<<"Track "<<lab<<" was not found !\n";
605 continue;
606 }
607 //track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
608 track->PropagateTo(3.,0.0028,65.19);
609
610 track->PropagateToVertex();
611
612 if (lab==tlab) hfound->Fill(ptg);
613 else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
614
615 Double_t xv,par[5]; track->GetExternalParameters(xv,par);
616 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
617 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
618 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
619 Float_t lam=TMath::ATan(par[3]);
620 Float_t pt_1=TMath::Abs(par[4]);
621 Double_t phig=TMath::ATan2(pyg,pxg);
622 //hp->Fill((phi - phig)*1000.);
623
624 Double_t lamg=TMath::ATan2(pzg,ptg);
625 //hl->Fill((lam - lamg)*1000.);
626
627 Double_t d=10000*track->GetD();
628 //hmpt->Fill(d);
629
630 //hptw->Fill(ptg,TMath::Abs(d));
631
632 Double_t z=10000*track->GetZ();
633 //hz->Fill(z);
634
635
636
637// ------------------------------------------------ b.b. -----
638 Int_t nev=0;
639 Float_t mom=1000.*(1./(pt_1*TMath::Cos(lam)));
640 Float_t dedx=track->GetdEdx();
641 //hep->Fill(mom,dedx,1.);
642 Int_t pcode=gt[ngood].code;
643 pid->SetEdep(10000*nev+ngood,dedx);
644 pid->SetPmom(10000*nev+ngood,mom);
645 pid->SetPcod(10000*nev+ngood,abs(pcode));
646 //cout<<"!!!!! pcode, pmod, dedx ="<<pcode<<","<<mom<<","<<dedx<<endl;
647
648 }
649
650 pid->Tab();
651
652 Stat_t ng=hgood->GetEntries();
653 //cerr<<"Good tracks "<<ng<<endl;
654 Stat_t nf=hfound->GetEntries();
655 Stat_t nfak=hfake->GetEntries();
656 if (ng!=0)
657 cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
658 //delete gAlice; //b.b.
659 cout<<"Done with nfound(good): "<<nf<<" + nfake: "<<nfak<<" = "<<nf+nfak<<" tracks."<<endl;
660// -------------------- b.b. -------------------
661
662 //return 0;
663}
664
665Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event) {
666 if (gAlice) {delete gAlice; gAlice=0;}
667
668 TFile *file=TFile::Open("galice.root");
669 if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
670 if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
671 cerr<<"gAlice have not been found on galice.root !\n";
672 exit(5);
673 }
674
675 Int_t np=gAlice->GetEvent(event);
676
677 Int_t *good=new Int_t[np];
678 Int_t k;
679 for (k=0; k<np; k++) good[k]=0;
680
681 AliITS *ITS=(AliITS*)gAlice->GetDetector("ITS");
682 if (!ITS) {
683 cerr<<"can't get ITS !\n"; exit(8);
684 }
685 AliITSgeom *geom=ITS->GetITSgeom();
686 if (!geom) {
687 cerr<<"can't get ITS geometry !\n"; exit(9);
688 }
689
690 TFile *cf=TFile::Open("AliITSclustersV2.root");
691 if (!cf->IsOpen()){
692 cerr<<"Can't open AliITSclustersV2.root !\n"; exit(6);
693 }
694 char cname[100]; sprintf(cname,"TreeC_ITS_%d",event);
695 TTree *cTree=(TTree*)cf->Get(cname);
696 if (!cTree) {
697 cerr<<"Can't get cTree !\n"; exit(7);
698 }
699 TBranch *branch=cTree->GetBranch("Clusters");
700 if (!branch) {
701 cerr<<"Can't get clusters branch !\n"; exit(8);
702 }
703 TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
704 branch->SetAddress(&clusters);
705
706 Int_t entr=(Int_t)cTree->GetEntries();
707 for (k=0; k<entr; k++) {
708 cTree->GetEvent(k);
709 Int_t ncl=clusters->GetEntriesFast(); if (ncl==0) continue;
710 Int_t lay,lad,det; geom->GetModuleId(k,lay,lad,det);
711 if (lay<1 || lay>6) {
712 cerr<<"wrong layer !\n"; exit(10);
713 }
714 while (ncl--) {
715 AliITSclusterV2 *pnt=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
716 Int_t l0=pnt->GetLabel(0);
717 Int_t l1=pnt->GetLabel(1);
718 Int_t l2=pnt->GetLabel(2);
719 Int_t mask=1<<(lay-1);
720 if (l0>=0) good[l0]|=mask;
721 if (l1>=0) good[l1]|=mask;
722 if (l2>=0) good[l2]|=mask;
723 }
724 }
725 clusters->Delete(); delete clusters;
726 delete cTree; //Thanks to Mariana Bondila
727 cf->Close();
728
729 ifstream in("good_tracks_tpc");
730 if (!in) {
731 cerr<<"can't get good_tracks_tpc !\n"; exit(11);
732 }
733 Int_t nt=0;
734 Double_t px,py,pz,x,y,z;
735 Int_t code,lab;
736 while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
737 if (good[lab] != 0x3F) continue;
738 TParticle *p = (TParticle*)gAlice->Particle(lab);
739 gt[nt].lab=lab;
740 gt[nt].code=p->GetPdgCode();
741 //**** px py pz - in global coordinate system
742 gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
743 gt[nt].x=gt[nt].y=gt[nt].z=0.;
744 nt++;
745 if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
746 }
747
748 delete[] good;
749
750 delete gAlice; gAlice=0;
751 file->Close();
752
753 return nt;
754}
755
756
757//----------------------------------------
758void filltab(){
759 cout<<"Fill tab..";
760 int nev=0;
761 int track,pcode;
762 float signal,pmom;
763 for(int t=0;t<10000;t++){
764 track=t; signal=5.; pmom=1200.; pcode=321;
765 pid->SetEdep(10000*nev+track,signal);
766 pid->SetPmom(10000*nev+track,pmom);
767 pid->SetPcod(10000*nev+track,abs(pcode));
768 }
769 pid->Tab();
770 cout<<"Done."<<endl;
771}
772//
773void loadpid(){
774 if(pid==0){
775 TFile *f=new TFile("pidhit.root");
776 AliITSPid *pid=(AliITSPid*)f->Get("AliITSPid");
777 if(pid)
778 {cout<<"Load PID object from PIDHIT.ROOT file"<<endl;}
779 else{cout<<"ERROR load PID object from PIDHIT.ROOT file"<<endl;}
780 }
781 pid->Print(0);
782}
783
784void quit(){
785 gROOT->Reset();
786 gROOT->ProcessLine(".q");
787}
788//---------------------------------------
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804