1 //------------------------------------------------------------//
2 // Class for identification of pions,kaons and protons in ITS //
3 // Prior particles population (probabilities) are taken from //
4 // Hijing event generator. //
5 //------------------------------------------------------------//
10 #include <TClonesArray.h>
12 #include "AliKalmanTrack.h"
13 #include "AliITSIOTrack.h"
14 #include "AliITStrackV2.h"
20 AliITSPid::AliITSPid(const AliITSPid &source) : TObject(source){
21 // Copy constructor. This is a function which is not allowed to be
22 // done: protected It exits with an error.
24 // AliITSPid &source An AliITSPid class.
30 Error("AliITSPid","You are not allowed to make a copy of the AliITSPid\n");
32 //______________________________________________________________________
33 AliITSPid& AliITSPid::operator=(const AliITSPid& /* source */){
34 // Assignment operator. This is a function which is not allowed to be
35 // done to the ITS. It exits with an error.
37 // AliITSPid &source An AliITSPid class.
43 Error("=operator","Assignment operator not allowed for this class\n");
48 Float_t AliITSPid::Qtrm(Int_t track) {
50 // This calculates truncated mean signal for given track.
52 TVector q(*( this->GetVec(track) ));
57 Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
60 case 1:q(6)=q(1); break;
61 case 2:q(6)=(q(1)+q(2))/2.; break;
63 for(int fi=0;fi<2;fi++){
65 do{ swap=0; float qmin=vf[fi];
66 for(int j=fi+1;j<nl;j++)
67 if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
70 q(6)= (vf[0]+vf[1])/2.;
73 for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
77 Float_t AliITSPid::Qtrm(Float_t qarr[6],Int_t narr) const{
79 //This calculates truncated mean signal for given signal set.
83 if(narr>0&&narr<7){ml=narr;}else{return 0;};
84 nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
87 case 1:qm=q[0]; break;
88 case 2:qm=(q[0]+q[1])/2.; break;
91 for(int fi=0;fi<2;fi++){
92 do{ swap=0; qmin=q[fi];
93 for(int j=fi+1;j<nl;j++)
94 if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
103 Int_t AliITSPid::Wpik(Float_t pm,Float_t q){
104 //Calcutates probabilityes of pions and kaons
105 //Returns particle code for dominant probability.
107 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
108 fggpi->SetParameters(par);
110 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
111 fggka->SetParameters(par);
113 Float_t ppi=fggpi->Eval(q);
114 Float_t pka=fggka->Eval(q);
120 if(p>0)cout<<" ppi,pka="<<ppi/p<<" "<<pka/p<<endl;
127 fWp=0.; fWpi=ppi; fWk=pka;
128 if( pka>ppi){return fPcode=321;}else{return fPcode=211;}
131 //-----------------------------------------------------------
132 Int_t AliITSPid::Wpikp(Float_t pm,Float_t q){
134 //Calcutates probabilityes of pions,kaons and protons.
135 //Returns particle code for dominant probability.
137 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
138 fggpi->SetParameters(par);
140 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
141 fggka->SetParameters(par);
143 for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
144 fggpr->SetParameters(par);
146 Float_t p,ppi,pka,ppr;
147 if( q>(fggpr->GetParameter(1)+fggpr->GetParameter(2)) )
148 { p=1.0; ppr=1.0; ppi=pka=0.0;
159 fWp=ppr; fWpi=ppi; fWk=pka;
160 //if(!fSilent)cout<<" ppi,pka,ppr="<<ppi<<" "<<pka<<" "<<ppr<<endl;
162 if( ppi>pka&&ppi>ppr )
164 else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
169 //-----------------------------------------------------------
170 Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
172 //Returns particle code
173 Info("GetPcode","method not implemented - Inputs TClonesArray *%x , Float_t %f",rps,pm);
176 //-----------------------------------------------------------
177 Int_t AliITSPid::GetPcode(AliKalmanTrack *track)
179 //Returns particle code for given track.
180 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
181 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
182 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
183 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
184 Float_t lam=TMath::ATan(par[3]);
185 Float_t pt1=TMath::Abs(par[4]);
186 Float_t mom=1./(pt1*TMath::Cos(lam));
187 Float_t dedx=track->GetPIDsignal();
188 Int_t pcode=GetPcode(dedx/40.,mom);
189 // cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
190 return pcode?pcode:211;
192 //------------------------------------------------------------
193 Int_t AliITSPid::GetPcode(AliITSIOTrack *track)
195 //Returns particle code for given track(V1).
200 Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
201 //???????????????????
203 Float_t dedx=track->GetdEdx();
204 //???????????????????
205 Int_t pcode=GetPcode(dedx,mom);
206 // cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
207 return pcode?pcode:211;
209 //-----------------------------------------------------------
210 Int_t AliITSPid::GetPcode(AliITStrackV2 *track)
212 //Returns particle code for given track(V2).
213 if(track==0)return 0;
214 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
215 track->PropagateTo(3.,0.0028,65.19);
216 track->PropagateToVertex();
217 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
218 Float_t lam=TMath::ATan(par[3]);
219 Float_t pt1=TMath::Abs(par[4]);
221 if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
222 Float_t dedx=track->GetdEdx();
223 // cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
224 Int_t pcode=GetPcode(dedx,mom);
225 // cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
226 return pcode?pcode:211;
228 //-----------------------------------------------------------
229 Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
231 //Returns particle code for given signal and momentum.
232 fWpi=fWk=fWp=0.; fPcode=0;
235 { if( q<fCutKa->Eval(pm) )
237 else{ if( q<fCutPr->Eval(pm) )
239 else{return Proton();}
243 if ( q>fCutPr->Eval(pm) )
244 {return Proton();} else {return Wpik(pm,q);};
245 if( pm<=1.10 ){ return Wpikp(pm,q); }
248 //-----------------------------------------------------------
249 void AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
250 Float_t klo,Float_t khi,Float_t plo,Float_t phi)
252 // The cut-table initializer method.
262 //------------------------------------------------------------
263 void AliITSPid::SetVec(Int_t ntrack,const TVector& info) const
265 //Store track info in tracls table
266 TClonesArray& arr=*fTrs;
267 new( arr[ntrack] ) TVector(info);
269 //-----------------------------------------------------------
270 TVector* AliITSPid::GetVec(Int_t ntrack) const
272 //Get given track from track table
273 TClonesArray& arr=*fTrs;
274 return (TVector*)arr[ntrack];
276 //-----------------------------------------------------------
277 void AliITSPid::SetEdep(Int_t track,Float_t Edep)
279 //Set dEdx for given track
281 if( ((TVector*)fTrs->At(track))->IsValid() )
282 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
283 Int_t j=(Int_t)xx(0); if(j>4)return;
284 xx(++j)=Edep;xx(0)=j;
285 TClonesArray &arr=*fTrs;
286 new(arr[track])TVector(xx);
288 //-----------------------------------------------------------
289 void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
291 //Set momentum for given track
293 if( ((TVector*)fTrs->At(track))->IsValid() )
294 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
296 TClonesArray &arr=*fTrs;
297 new(arr[track])TVector(xx);
299 //-----------------------------------------------------------
300 void AliITSPid::SetPcod(Int_t track,Int_t partcode)
302 //Set particle code for given track
304 if( ((TVector*)fTrs->At(track))->IsValid() )
305 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
307 {xx(11)=partcode; fMxtrs++;
308 TClonesArray &arr=*fTrs;
309 new(arr[track])TVector(xx);
312 //-----------------------------------------------------------
313 void AliITSPid::Print(Int_t track)
315 //Prints information for given track
316 cout<<fMxtrs<<" tracks in AliITSPid obj."<<endl;
317 if( ((TVector*)fTrs->At(track))->IsValid() )
318 {TVector xx( *((TVector*)fTrs->At(track)) );
322 {cout<<"No data for track "<<track<<endl;return;}
324 //-----------------------------------------------------------
325 void AliITSPid::Tab(void)
327 //Make PID for tracks stored in tracks table
328 if(fTrs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
329 cout<<"------------------------------------------------------------------------"<<endl;
330 cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
331 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
332 cout<<"------------------------------------------------------------------------"<<endl;
333 for(Int_t i=0;i<fTrs->GetEntries();i++)
335 TVector xx( *((TVector*)fTrs->At(i)) );
336 if( xx.IsValid() && xx(0)>0 )
338 TVector xx( *((TVector*)fTrs->At(i)) );
342 xx(6)=(this->Qtrm(i));
347 // 2)Calculate Wpi,Wk,Wp
348 this->GetPcode(xx(6),xx(10)/1000.);
355 for(Int_t j=1;j<11;j++){
356 if(i<7){ cout.width(7);cout.precision(4);cout<<xx(j);}
357 if(i>7){ cout.width(7);cout.precision(5);cout<<xx(j);}
361 // 4)Update data in TVector
362 TClonesArray &arr=*fTrs;
363 new(arr[i])TVector(xx);
366 {/*cout<<"No data for track "<<i<<endl;*/}
367 }// End loop for tracks
369 void AliITSPid::Reset(void)
372 for(Int_t i=0;i<fTrs->GetEntries();i++){
374 TClonesArray &arr=*fTrs;
375 new(arr[i])TVector(xx);
378 //-----------------------------------------------------------
379 AliITSPid::AliITSPid(Int_t ntrack)
381 //Constructor for AliITSPid class
383 fTrs = new TClonesArray("TVector",ntrack);
384 TClonesArray &arr=*fTrs;
385 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
388 fCutKa=new TF1("fcutka","pol4",0.05,0.4);
389 Double_t ka[5]={25.616, -161.59, 408.97, -462.17, 192.86};
390 fCutKa->SetParameters(ka);
392 fCutPr=new TF1("fcutpr","[0]/x/x+[1]",0.05,1.1);
393 Double_t pr[2]={0.70675,0.4455};
394 fCutPr->SetParameters(pr);
396 //---------- signal fit ----------
398 fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
399 Double_t parpi0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
400 1.1405329025e+05, -3.8289694496e+04 };
401 fGGpi[0]->SetParameters(parpi0);
402 fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
403 Double_t parpi1[10]={ 1.0791668283e-02, 9.7347716496e-01 };
404 fGGpi[1]->SetParameters(parpi1);
405 fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
406 Double_t parpi2[10]={ 5.8191602279e-04, 9.7285601334e-02 };
407 fGGpi[2]->SetParameters(parpi2);
408 fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
409 Double_t parpi3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
410 6.2900977606e+03, -2.2935862292e+03 };
411 fGGpi[3]->SetParameters(parpi3);
412 fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
413 Double_t parpi4[10]={ 9.0419011783e-03, 1.1628922525e+00 };
414 fGGpi[4]->SetParameters(parpi4);
415 fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
416 Double_t parpi5[10]={ 1.8324872519e-03, 2.1503968838e-01 };
417 fGGpi[5]->SetParameters(parpi5);
420 fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
421 Double_t parka0[20]={
422 -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
423 -3.7123338009e+03,1.6003647641e+03 };
424 fGGka[0]->SetParameters(parka0);
425 fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
426 Double_t parka1[20]={
427 2.5181172905e-01,8.7566001814e-01 };
428 fGGka[1]->SetParameters(parka1);
429 fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
430 Double_t parka2[20]={
431 8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
432 -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
434 fGGka[2]->SetParameters(parka2);
437 fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
438 Double_t parpr0[10]={
439 6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
440 -3.5269200901e+03,1.2859128345e+03 };
441 fGGpr[0]->SetParameters(parpr0);
442 fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
443 Double_t parpr1[10]={
444 9.4970837607e-01,7.3573504201e-01 };
445 fGGpr[1]->SetParameters(parpr1);
446 fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
447 Double_t parpr2[10]={
448 1.2498403757e-01,2.7845072306e-02 };
449 fGGpr[2]->SetParameters(parpr2);
451 //----------- end fit -----------
453 fggpr=new TF1("ggpr","gaus",0.4,1.2);
454 fggpi=new TF1("ggpi","gaus+gaus(3)",0.4,1.2);
455 fggka=new TF1("ggka","gaus",0.4,1.2);
457 //-------------------------------------------------
459 // Ncut Pmom pilo pihi klo khi plo phi
460 // cut[j] [0] [1] [2] [3] [4] [5] [6]
461 //----------------------------------------------------------------
462 SetCut( 1, 0.12 , 0. , 0. , kInf , kInf , kInf , kInf );
463 SetCut( 2, 0.20 , 0. , 6.0 , 6.0 , kInf , kInf , kInf );
464 SetCut( 3, 0.30 , 0. , 3.5 , 3.5 , 9.0 , 9.0 , kInf );
465 SetCut( 4, 0.41 , 0. , 1.9 , 1.9 , 4.0 , 4.0 , kInf );
466 //----------------------------------------------------------------
467 SetCut( 5, 0.47 , 0.935, 0.139, 1.738 , 0.498 , 3.5 , kInf ); //410-470
468 SetCut( 6, 0.53 , 0.914, 0.136, 1.493 , 0.436 , 3.0 , kInf ); //470-530
469 //----------------------------------------------------------------
470 SetCut( 7, 0.59 , 0.895, 0.131, 1.384 , 0.290 , 2.7 , kInf ); //530-590
471 SetCut( 8, 0.65 , 0.887, 0.121, 1.167 , 0.287 , 2.5 , kInf ); //590-650
472 SetCut( 9, 0.73 , 0.879, 0.120, 1.153 , 0.257 , 2.0 , kInf ); //650-730
473 //----------------------------------------------------------------
474 SetCut( 10, 0.83 , 0.880, 0.126, 1.164 , 0.204 , 2.308 , 0.297 ); //730-830
475 SetCut( 11, 0.93 , 0.918, 0.145, 1.164 , 0.204 , 2.00 , 0.168 ); //830-930
476 SetCut( 12, 1.03 , 0.899, 0.128, 1.164 , 0.204 ,1.80 , 0.168);
477 //------------------------ pi,K ---------------------
478 fAprob[0][0]=1212; fAprob[1][0]=33.; // fAprob[0][i] - const for pions,cut[i+5]
479 fAprob[0][1]=1022; fAprob[1][1]=46.2 ; // fAprob[1][i] - kaons
480 //---------------------------------------------------
481 fAprob[0][2]= 889.7; fAprob[1][2]=66.58; fAprob[2][2]=14.53;
482 fAprob[0][3]= 686.; fAprob[1][3]=88.8; fAprob[2][3]=19.27;
483 fAprob[0][4]= 697.; fAprob[1][4]=125.6; fAprob[2][4]=28.67;
484 //------------------------ pi,K,p -------------------
485 fAprob[0][5]= 633.7; fAprob[1][5]=100.1; fAprob[2][5]=37.99; // fAprob[2][i] - protons
486 fAprob[0][6]= 469.5; fAprob[1][6]=20.74; fAprob[2][6]=25.43;
487 fAprob[0][7]= 355.; fAprob[1][7]=
488 355.*(20.74/469.5); fAprob[2][7]=34.08;