ITS/AliITSPid.cxx

   1 #include "AliITSPid.h"
   2 #include "TMath.h"
   3 #include <Riostream.h>
   4 #include <TClonesArray.h>
   5 #include <TVector.h>
   6 #include "AliTPCtrack.h"
   7 #include "AliITStrackV2.h"
   8 #include <TF1.h>
   9
  10 ClassImp(AliITSPid)
  11 //
  12 Float_t AliITSPid::Qtrm(Int_t track)
  13 {
  14     TVector q(*( this->GetVec(track)  ));
  15     Int_t ml=(Int_t)q(0);
  16     if(ml<1)return 0.;
  17     if(ml>6)ml=6;
  18     float vf[6];
  19     Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
  20     if(nl==0)return 0.;
  21     switch(nl){
  22       case  1:q(6)=q(1); break;
  23       case  2:q(6)=(q(1)+q(2))/2.; break;
  24       default:
  25         for(int fi=0;fi<2;fi++){
  26           Int_t swap;
  27           do{ swap=0; float qmin=vf[fi];
  28           for(int j=fi+1;j<nl;j++)
  29             if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
  30           }while(swap==1);
  31         }
  32         q(6)= (vf[0]+vf[1])/2.;
  33         break;
  34     }
  35     for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
  36     return (q(6));
  37 }
  38
  39 Float_t AliITSPid::Qtrm(Float_t qarr[6],Int_t narr)
  40 {
  41   Float_t q[6],qm,qmin;
  42   Int_t nl,ml;
  43   if(narr>0&&narr<7){ml=narr;}else{return 0;};
  44   nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
  45   if(nl==0)return 0.;
  46     switch(nl){
  47       case  1:qm=q[0]; break;
  48       case  2:qm=(q[0]+q[1])/2.; break;
  49       default:
  50         Int_t swap;
  51         for(int fi=0;fi<2;fi++){
  52           do{ swap=0; qmin=q[fi];
  53           for(int j=fi+1;j<nl;j++)
  54             if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
  55           }while(swap==1);
  56         }
  57         qm= (q[0]+q[1])/2.;
  58         break;
  59     }
  60     return qm;
  61 }
  62
  63 Int_t   AliITSPid::Wpik(Float_t pm,Float_t q)
  64 {
  65   Double_t par[6];
  66   for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
  67   fggpi->SetParameters(par);
  68
  69   for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
  70   fggka->SetParameters(par);
  71
  72   Float_t ppi=fggpi->Eval(q);
  73   Float_t pka=fggka->Eval(q);
  74   Float_t p=ppi+pka;
  75   /*
  76   if(!fSilent){
  77     fggka->Print();
  78     fggpi->Print();
  79     if(p>0)cout<<" ppi,pka="<<ppi/p<<"  "<<pka/p<<endl;
  80   }
  81   */
  82
  83   if(p>0){
  84     ppi=ppi/p;
  85     pka=pka/p;
  86     fWp=0.; fWpi=ppi; fWk=pka;
  87     if( pka>ppi){return fPcode=321;}else{return fPcode=211;}
  88   }else{return 0;}
  89 }
  90 //-----------------------------------------------------------
  91 Int_t   AliITSPid::Wpikp(Float_t pm,Float_t q)
  92 {
  93   Double_t par[6];
  94   for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
  95   fggpi->SetParameters(par);
  96
  97   for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
  98   fggka->SetParameters(par);
  99
 100   for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
 101   fggpr->SetParameters(par);
 102
 103   Float_t p,ppi,pka,ppr;
 104   if( q>(fggpr->GetParameter(1)+fggpr->GetParameter(2)) )
 105       { p=1.0; ppr=1.0; ppi=pka=0.0;
 106     }else{
 107     ppi=fggpi->Eval(q);
 108     pka=fggka->Eval(q);
 109     ppr=fggpr->Eval(q);
 110     p=ppi+pka+ppr;
 111   }
 112   if(p>0){
 113     ppi=ppi/p;
 114     pka=pka/p;
 115     ppr=ppr/p;
 116     fWp=ppr; fWpi=ppi; fWk=pka;
 117     //if(!fSilent)cout<<" ppi,pka,ppr="<<ppi<<"  "<<pka<<" "<<ppr<<endl;
 118
 119    if( ppi>pka&&ppi>ppr )
 120            {return fPcode=211;}
 121    else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
 122    }
 123
 124   }else{return 0;}
 125 }
 126 //-----------------------------------------------------------
 127 Int_t   AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
 128 {
 129   Info("GetPcode","method not implemented - Inputs TClonesArray *%x , Float_t %f",rps,pm);
 130     return 0;
 131 }
 132 //-----------------------------------------------------------
 133 Int_t   AliITSPid::GetPcode(AliTPCtrack *track)
 134 {
 135       Double_t xk,par[5]; track->GetExternalParameters(xk,par);
 136       Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
 137       if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 138       if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 139       Float_t lam=TMath::ATan(par[3]);
 140       Float_t pt1=TMath::Abs(par[4]);
 141       Float_t mom=1./(pt1*TMath::Cos(lam));
 142       Float_t dedx=track->GetdEdx();
 143     Int_t pcode=GetPcode(dedx/40.,mom);
 144 //    cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
 145     return pcode?pcode:211;
 146     }
 147 //------------------------------------------------------------
 148 Int_t   AliITSPid::GetPcode(AliITSIOTrack *track)
 149 {
 150     Double_t px,py,pz;
 151     px=track->GetPx();
 152     py=track->GetPy();
 153     pz=track->GetPz();
 154     Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
 155 //???????????????????
 156     // Float_t dedx=1.0;
 157     Float_t dedx=track->GetdEdx();
 158 //???????????????????
 159     Int_t pcode=GetPcode(dedx,mom);
 160 //    cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
 161 return pcode?pcode:211;
 162 }
 163 //-----------------------------------------------------------
 164 Int_t   AliITSPid::GetPcode(AliITStrackV2 *track)
 165 {
 166   if(track==0)return 0;
 167   //      track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
 168       track->PropagateTo(3.,0.0028,65.19);
 169       track->PropagateToVertex();
 170     Double_t xk,par[5]; track->GetExternalParameters(xk,par);
 171     Float_t lam=TMath::ATan(par[3]);
 172     Float_t pt1=TMath::Abs(par[4]);
 173     Float_t mom=0.;
 174     if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
 175     Float_t dedx=track->GetdEdx();
 176 //    cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
 177     Int_t pcode=GetPcode(dedx,mom);
 178 //    cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
 179 return pcode?pcode:211;
 180 }
 181 //-----------------------------------------------------------
 182 Int_t   AliITSPid::GetPcode(Float_t q,Float_t pm)
 183 {
 184     fWpi=fWk=fWp=0.;     fPcode=0;
 185
 186     if ( pm<=0.400 )
 187         { if( q<fCutKa->Eval(pm) )
 188             {return Pion();}
 189         else{ if( q<fCutPr->Eval(pm) )
 190                 {return Kaon();}
 191             else{return Proton();}
 192             }
 193         }
 194     if ( pm<=0.750 )
 195         if ( q>fCutPr->Eval(pm)  )
 196             {return Proton();} else {return Wpik(pm,q);};
 197     if( pm<=1.10 ){ return Wpikp(pm,q); }
 198     return fPcode;
 199 }
 200 //-----------------------------------------------------------
 201 void    AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
 202                         Float_t klo,Float_t khi,Float_t plo,Float_t phi)
 203 {
 204     fCut[n][0]=pm;
 205     fCut[n][1]=pilo;
 206     fCut[n][2]=pihi;
 207     fCut[n][3]=klo;
 208     fCut[n][4]=khi;
 209     fCut[n][5]=plo;
 210     fCut[n][6]=phi;
 211     return ;
 212 }
 213 //------------------------------------------------------------
 214 void AliITSPid::SetVec(Int_t ntrack,TVector info)
 215 {
 216 TClonesArray& arr=*fTrs;
 217     new( arr[ntrack] ) TVector(info);
 218 }
 219 //-----------------------------------------------------------
 220 TVector* AliITSPid::GetVec(Int_t ntrack)
 221 {
 222 TClonesArray& arr=*fTrs;
 223     return (TVector*)arr[ntrack];
 224 }
 225 //-----------------------------------------------------------
 226 void AliITSPid::SetEdep(Int_t track,Float_t Edep)
 227 {
 228     TVector xx(0,11);
 229     if( ((TVector*)fTrs->At(track))->IsValid() )
 230         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
 231     Int_t j=(Int_t)xx(0); if(j>4)return;
 232     xx(++j)=Edep;xx(0)=j;
 233     TClonesArray &arr=*fTrs;
 234     new(arr[track])TVector(xx);
 235 }
 236 //-----------------------------------------------------------
 237 void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
 238 {
 239     TVector xx(0,11);
 240     if( ((TVector*)fTrs->At(track))->IsValid() )
 241         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
 242     xx(10)=Pmom;
 243     TClonesArray &arr=*fTrs;
 244     new(arr[track])TVector(xx);
 245 }
 246 //-----------------------------------------------------------
 247 void AliITSPid::SetPcod(Int_t track,Int_t partcode)
 248 {
 249     TVector xx(0,11);
 250     if( ((TVector*)fTrs->At(track))->IsValid() )
 251         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
 252     if(xx(11)==0)
 253         {xx(11)=partcode; fMxtrs++;
 254         TClonesArray &arr=*fTrs;
 255         new(arr[track])TVector(xx);
 256         }
 257 }
 258 //-----------------------------------------------------------
 259 void AliITSPid::Print(Int_t track)
 260 {cout<<fMxtrs<<" tracks in AliITSPid obj."<<endl;
 261     if( ((TVector*)fTrs->At(track))->IsValid() )
 262         {TVector xx( *((TVector*)fTrs->At(track)) );
 263          xx.Print();
 264          }
 265     else
 266         {cout<<"No data for track "<<track<<endl;return;}
 267 }
 268 //-----------------------------------------------------------
 269 void AliITSPid::Tab(void)
 270 {
 271 if(fTrs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
 272 cout<<"------------------------------------------------------------------------"<<endl;
 273 cout<<"Nq"<<"   q1  "<<"   q2  "<<"   q3  "<<"   q4  "<<"   q5   "<<
 274       " Qtrm    "    <<"  Wpi  "<<"  Wk   "<<"  Wp  "<<"Pmom  "<<endl;
 275 cout<<"------------------------------------------------------------------------"<<endl;
 276 for(Int_t i=0;i<fTrs->GetEntries();i++)
 277 {
 278   TVector xx( *((TVector*)fTrs->At(i)) );
 279     if( xx.IsValid() && xx(0)>0 )
 280         {
 281             TVector xx( *((TVector*)fTrs->At(i)) );
 282             if(xx(0)>=2)
 283                 {
 284 //       1)Calculate Qtrm
 285                     xx(6)=(this->Qtrm(i));
 286
 287                  }else{
 288                      xx(6)=xx(1);
 289                  }
 290 //       2)Calculate Wpi,Wk,Wp
 291             this->GetPcode(xx(6),xx(10)/1000.);
 292             xx(7)=GetWpi();
 293             xx(8)=GetWk();
 294             xx(9)=GetWp();
 295 //       3)Print table
 296             if(xx(0)>0){
 297 //                  cout<<xx(0)<<" ";
 298                     for(Int_t j=1;j<11;j++){
 299                       if(i<7){ cout.width(7);cout.precision(4);cout<<xx(j);}
 300                       if(i>7){ cout.width(7);cout.precision(5);cout<<xx(j);}
 301                     }
 302                     cout<<endl;
 303                 }
 304 //        4)Update data in TVector
 305             TClonesArray &arr=*fTrs;
 306             new(arr[i])TVector(xx);
 307         }
 308     else
 309       {/*cout<<"No data for track "<<i<<endl;*/}
 310 }// End loop for tracks
 311 }
 312 void AliITSPid::Reset(void)
 313 {
 314   for(Int_t i=0;i<fTrs->GetEntries();i++){
 315     TVector xx(0,11);
 316     TClonesArray &arr=*fTrs;
 317     new(arr[i])TVector(xx);
 318   }
 319 }
 320 //-----------------------------------------------------------
 321 AliITSPid::AliITSPid(Int_t ntrack)
 322 {
 323   fSigmin=0.01;
 324     fTrs = new TClonesArray("TVector",ntrack);
 325     TClonesArray &arr=*fTrs;
 326     for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
 327     fMxtrs=0;
 328     //
 329     fCutKa=new TF1("fcutka","pol4",0.05,0.4);
 330     Double_t ka[5]={25.616, -161.59, 408.97, -462.17, 192.86};
 331     fCutKa->SetParameters(ka);
 332     //
 333     fCutPr=new TF1("fcutpr","[0]/x/x+[1]",0.05,1.1);
 334     Double_t pr[2]={0.70675,0.4455};
 335     fCutPr->SetParameters(pr);
 336     //
 337     //---------- signal fit ----------
 338 {//Pions
 339 fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
 340   Double_t parpi0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
 341    1.1405329025e+05, -3.8289694496e+04  };
 342   fGGpi[0]->SetParameters(parpi0);
 343 fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
 344   Double_t parpi1[10]={ 1.0791668283e-02, 9.7347716496e-01  };
 345   fGGpi[1]->SetParameters(parpi1);
 346 fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
 347   Double_t parpi2[10]={ 5.8191602279e-04, 9.7285601334e-02  };
 348   fGGpi[2]->SetParameters(parpi2);
 349 fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
 350   Double_t parpi3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
 351    6.2900977606e+03, -2.2935862292e+03  };
 352   fGGpi[3]->SetParameters(parpi3);
 353 fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
 354   Double_t parpi4[10]={ 9.0419011783e-03, 1.1628922525e+00  };
 355   fGGpi[4]->SetParameters(parpi4);
 356 fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
 357   Double_t parpi5[10]={ 1.8324872519e-03, 2.1503968838e-01  };
 358   fGGpi[5]->SetParameters(parpi5);
 359 }//End Pions
 360 {//Kaons
 361 fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
 362   Double_t parka0[20]={
 363   -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
 364   -3.7123338009e+03,1.6003647641e+03  };
 365   fGGka[0]->SetParameters(parka0);
 366 fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
 367   Double_t parka1[20]={
 368   2.5181172905e-01,8.7566001814e-01  };
 369   fGGka[1]->SetParameters(parka1);
 370 fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
 371   Double_t parka2[20]={
 372   8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
 373   -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
 374   5.7545491696e+01  };
 375   fGGka[2]->SetParameters(parka2);
 376 }//End Kaons
 377 {//Protons
 378 fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
 379   Double_t parpr0[10]={
 380   6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
 381   -3.5269200901e+03,1.2859128345e+03  };
 382   fGGpr[0]->SetParameters(parpr0);
 383 fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
 384   Double_t parpr1[10]={
 385   9.4970837607e-01,7.3573504201e-01  };
 386   fGGpr[1]->SetParameters(parpr1);
 387 fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
 388   Double_t parpr2[10]={
 389   1.2498403757e-01,2.7845072306e-02  };
 390   fGGpr[2]->SetParameters(parpr2);
 391 }//End Protons
 392     //----------- end fit -----------
 393
 394     fggpr=new TF1("ggpr","gaus",0.4,1.2);
 395     fggpi=new TF1("ggpi","gaus+gaus(3)",0.4,1.2);
 396     fggka=new TF1("ggka","gaus",0.4,1.2);
 397
 398     //-------------------------------------------------
 399 const int kInf=10;
 400 //         Ncut Pmom   pilo  pihi    klo    khi     plo    phi
 401 //       cut[j] [0]    [1]    [2]    [3]    [4]     [5]    [6]
 402 //----------------------------------------------------------------
 403     SetCut(  1, 0.12 ,  0.  ,  0.  , kInf  , kInf   , kInf  , kInf  );
 404     SetCut(  2, 0.20 ,  0.  ,  6.0 , 6.0  , kInf   , kInf  , kInf  );
 405     SetCut(  3, 0.30 ,  0.  ,  3.5 , 3.5  , 9.0   , 9.0  , kInf  );
 406     SetCut(  4, 0.41 ,  0.  ,  1.9 , 1.9  , 4.0   , 4.0  , kInf  );
 407 //----------------------------------------------------------------
 408     SetCut(  5, 0.47 , 0.935, 0.139, 1.738 , 0.498  , 3.5  , kInf  );  //410-470
 409     SetCut(  6, 0.53 , 0.914, 0.136, 1.493 , 0.436  , 3.0  , kInf  );  //470-530
 410 //----------------------------------------------------------------
 411     SetCut(  7, 0.59 , 0.895, 0.131, 1.384 , 0.290 , 2.7  , kInf  );    //530-590
 412     SetCut(  8, 0.65 , 0.887, 0.121, 1.167 , 0.287 , 2.5  , kInf  );     //590-650
 413     SetCut(  9, 0.73 , 0.879, 0.120, 1.153 , 0.257 , 2.0  , kInf  );     //650-730
 414 //----------------------------------------------------------------
 415     SetCut( 10, 0.83 , 0.880, 0.126, 1.164 , 0.204 , 2.308 , 0.297 );       //730-830
 416     SetCut( 11, 0.93 , 0.918, 0.145, 1.164 , 0.204 , 2.00 , 0.168 );        //830-930
 417     SetCut( 12, 1.03 , 0.899, 0.128, 1.164 , 0.204  ,1.80 , 0.168);
 418     //------------------------ pi,K ---------------------
 419     fAprob[0][0]=1212;     fAprob[1][0]=33.;   // fAprob[0][i] - const for pions,cut[i+5]
 420     fAprob[0][1]=1022;     fAprob[1][1]=46.2 ; // fAprob[1][i] -           kaons
 421     //---------------------------------------------------
 422     fAprob[0][2]= 889.7;   fAprob[1][2]=66.58; fAprob[2][2]=14.53;
 423     fAprob[0][3]= 686.;    fAprob[1][3]=88.8;  fAprob[2][3]=19.27;
 424     fAprob[0][4]= 697.;    fAprob[1][4]=125.6; fAprob[2][4]=28.67;
 425     //------------------------ pi,K,p -------------------
 426     fAprob[0][5]= 633.7;   fAprob[1][5]=100.1;   fAprob[2][5]=37.99;   // fAprob[2][i] -  protons
 427     fAprob[0][6]= 469.5;   fAprob[1][6]=20.74;   fAprob[2][6]=25.43;
 428     fAprob[0][7]= 355.;    fAprob[1][7]=
 429                           355.*(20.74/469.5);  fAprob[2][7]=34.08;
 430 }
 431 //End AliITSPid.cxx