ITS/AliITSPid.cxx

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