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