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       //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
 124 void 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
 163 void 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
 189 void 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
 215 void 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 //---------------------------------------
 237 void 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.---------------------------------
 256 void 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 //---------------------------------------
 274 void 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 //---------------------------------------
 291 void 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 //---------------------------------------
 306 void effall(){
 307
 308       eff(211,211,"",3);
 309       eff(321,321,"same",2);
 310      eff(2212,2212,"same",1);
 311 }
 312 //---------------------------------------
 313 void 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);
 339 return;
 340     pmom->SetFillColor(16);    pmom->Draw();
 341     efpi->SetFillColor(0); //42
 342     efpi->Draw("same");
 343 }
 344 //---------------------------------------
 345 Int_t NRange=0;
 346     Float_t xpmin[]={410,470,730,830,930};
 347     Float_t xpmax[]={470,530,830,930,1030};
 348 void 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 //---------------------------------------
 355 void 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.-------------------------------------
 367 void 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 //---------------------------------------
 382 void 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 //-------------
 393 void 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 //---------------------------------------
 405 void 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 //--------------------------------------
 413 void pcode(){
 414     TH1F *his =  new TH1F("pcode","tit",100,10,2300);
 415 Float_t Pcod;
 416 for(Int_t i=0;i<NStat;i++){
 417     TVector tabv(*( pid->GetVec(i) ));
 418     Pcod=tabv(11); if(Pcod>0)  his->Fill(Pcod);
 419 }
 420 his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
 421 }
 422
 423 void 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
 436 void 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
 447 int jtrack=0;
 448 void 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     }
 455 cout<<"Track No "<<jtrack<<endl;
 456 pid->Print(jtrack++);
 457 }
 458
 459 // Fill histogram for track number
 460 //--------------------------------
 461 void 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
 495 struct 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
 504 void filltab_tracks(){
 505
 506   cerr<<"Filling of track table...\n";
 507
 508   Int_t event=0;
 509 Int_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
 665 Int_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 //----------------------------------------
 758 void 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 //
 773 void 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
 784 void 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