bugfix
[u/mrichter/AliRoot.git] / ITS / AliITSPid.cxx
CommitLineData
23efe5f1 1#include "AliITSPid.h"
2#include "TMath.h"
4ae5bbc4 3#include <Riostream.h>
e75c69cd 4#include <TClonesArray.h>
5#include <TVector.h>
6#include "AliTPCtrack.h"
7#include "AliITStrackV2.h"
8#include <TF1.h>
9
23efe5f1 10ClassImp(AliITSPid)
11b3ae13 11//
e75c69cd 12Float_t AliITSPid::Qtrm(Int_t track)
23efe5f1 13{
14 TVector q(*( this->GetVec(track) ));
38f7c306 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.;
23efe5f1 21 switch(nl){
22 case 1:q(6)=q(1); break;
23 case 2:q(6)=(q(1)+q(2))/2.; break;
24 default:
38f7c306 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 }
23efe5f1 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}
38f7c306 38
e75c69cd 39Float_t AliITSPid::Qtrm(Float_t qarr[6],Int_t narr)
79b921e1 40{
38f7c306 41 Float_t q[6],qm,qmin;
79b921e1 42 Int_t nl,ml;
38f7c306 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++;}}
79b921e1 45 if(nl==0)return 0.;
79b921e1 46 switch(nl){
47 case 1:qm=q[0]; break;
48 case 2:qm=(q[0]+q[1])/2.; break;
49 default:
38f7c306 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.;
79b921e1 58 break;
59 }
60 return qm;
61}
11b3ae13 62
e75c69cd 63Int_t AliITSPid::Wpik(Float_t pm,Float_t q)
23efe5f1 64{
11b3ae13 65 Double_t par[6];
66 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
67 fggpi->SetParameters(par);
23efe5f1 68
11b3ae13 69 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
70 fggka->SetParameters(par);
38f7c306 71
11b3ae13 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 */
38f7c306 82
11b3ae13 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;}
23efe5f1 89}
90//-----------------------------------------------------------
e75c69cd 91Int_t AliITSPid::Wpikp(Float_t pm,Float_t q)
23efe5f1 92{
11b3ae13 93 Double_t par[6];
94 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
95 fggpi->SetParameters(par);
23efe5f1 96
11b3ae13 97 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
98 fggka->SetParameters(par);
23efe5f1 99
11b3ae13 100 for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
101 fggpr->SetParameters(par);
38f7c306 102
11b3ae13 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;
38f7c306 118
11b3ae13 119 if( ppi>pka&&ppi>ppr )
120 {return fPcode=211;}
121 else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
122 }
38f7c306 123
11b3ae13 124 }else{return 0;}
23efe5f1 125}
126//-----------------------------------------------------------
127Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
128{
088e0b8d 129 Info("GetPcode","method not implemented - Inputs TClonesArray *%x , Float_t %f",rps,pm);
23efe5f1 130 return 0;
131}
132//-----------------------------------------------------------
79b921e1 133Int_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]);
e75c69cd 140 Float_t pt1=TMath::Abs(par[4]);
141 Float_t mom=1./(pt1*TMath::Cos(lam));
79b921e1 142 Float_t dedx=track->GetdEdx();
143 Int_t pcode=GetPcode(dedx/40.,mom);
def162f4 144// cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 145 return pcode?pcode:211;
146 }
147//------------------------------------------------------------
79b921e1 148Int_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//???????????????????
def162f4 156 // Float_t dedx=1.0;
157 Float_t dedx=track->GetdEdx();
79b921e1 158//???????????????????
159 Int_t pcode=GetPcode(dedx,mom);
def162f4 160// cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 161return pcode?pcode:211;
162}
79b921e1 163//-----------------------------------------------------------
164Int_t AliITSPid::GetPcode(AliITStrackV2 *track)
165{
166 if(track==0)return 0;
38f7c306 167 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
79b921e1 168 track->PropagateTo(3.,0.0028,65.19);
79b921e1 169 track->PropagateToVertex();
170 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
171 Float_t lam=TMath::ATan(par[3]);
e75c69cd 172 Float_t pt1=TMath::Abs(par[4]);
79b921e1 173 Float_t mom=0.;
e75c69cd 174 if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
79b921e1 175 Float_t dedx=track->GetdEdx();
e75c69cd 176// cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
79b921e1 177 Int_t pcode=GetPcode(dedx,mom);
def162f4 178// cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 179return pcode?pcode:211;
180}
181//-----------------------------------------------------------
23efe5f1 182Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
183{
38f7c306 184 fWpi=fWk=fWp=0.; fPcode=0;
23efe5f1 185
11b3ae13 186 if ( pm<=0.400 )
187 { if( q<fCutKa->Eval(pm) )
e75c69cd 188 {return Pion();}
11b3ae13 189 else{ if( q<fCutPr->Eval(pm) )
e75c69cd 190 {return Kaon();}
191 else{return Proton();}
11b3ae13 192 }
193 }
194 if ( pm<=0.750 )
195 if ( q>fCutPr->Eval(pm) )
e75c69cd 196 {return Proton();} else {return Wpik(pm,q);};
197 if( pm<=1.10 ){ return Wpikp(pm,q); }
23efe5f1 198 return fPcode;
199}
200//-----------------------------------------------------------
201void 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{
e75c69cd 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;
23efe5f1 211 return ;
212}
38f7c306 213//------------------------------------------------------------
23efe5f1 214void AliITSPid::SetVec(Int_t ntrack,TVector info)
215{
e75c69cd 216TClonesArray& arr=*fTrs;
23efe5f1 217 new( arr[ntrack] ) TVector(info);
218}
219//-----------------------------------------------------------
220TVector* AliITSPid::GetVec(Int_t ntrack)
221{
e75c69cd 222TClonesArray& arr=*fTrs;
23efe5f1 223 return (TVector*)arr[ntrack];
224}
225//-----------------------------------------------------------
226void AliITSPid::SetEdep(Int_t track,Float_t Edep)
227{
228 TVector xx(0,11);
e75c69cd 229 if( ((TVector*)fTrs->At(track))->IsValid() )
230 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
23efe5f1 231 Int_t j=(Int_t)xx(0); if(j>4)return;
232 xx(++j)=Edep;xx(0)=j;
e75c69cd 233 TClonesArray &arr=*fTrs;
23efe5f1 234 new(arr[track])TVector(xx);
235}
236//-----------------------------------------------------------
237void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
238{
239 TVector xx(0,11);
e75c69cd 240 if( ((TVector*)fTrs->At(track))->IsValid() )
241 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
23efe5f1 242 xx(10)=Pmom;
e75c69cd 243 TClonesArray &arr=*fTrs;
23efe5f1 244 new(arr[track])TVector(xx);
245}
246//-----------------------------------------------------------
247void AliITSPid::SetPcod(Int_t track,Int_t partcode)
248{
249 TVector xx(0,11);
e75c69cd 250 if( ((TVector*)fTrs->At(track))->IsValid() )
251 {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
23efe5f1 252 if(xx(11)==0)
e75c69cd 253 {xx(11)=partcode; fMxtrs++;
254 TClonesArray &arr=*fTrs;
23efe5f1 255 new(arr[track])TVector(xx);
256 }
257}
258//-----------------------------------------------------------
259void AliITSPid::Print(Int_t track)
e75c69cd 260{cout<<fMxtrs<<" tracks in AliITSPid obj."<<endl;
261 if( ((TVector*)fTrs->At(track))->IsValid() )
262 {TVector xx( *((TVector*)fTrs->At(track)) );
23efe5f1 263 xx.Print();
264 }
265 else
266 {cout<<"No data for track "<<track<<endl;return;}
267}
268//-----------------------------------------------------------
269void AliITSPid::Tab(void)
270{
e75c69cd 271if(fTrs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
23efe5f1 272cout<<"------------------------------------------------------------------------"<<endl;
273cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
274 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
275cout<<"------------------------------------------------------------------------"<<endl;
e75c69cd 276for(Int_t i=0;i<fTrs->GetEntries();i++)
23efe5f1 277{
e75c69cd 278 TVector xx( *((TVector*)fTrs->At(i)) );
23efe5f1 279 if( xx.IsValid() && xx(0)>0 )
280 {
e75c69cd 281 TVector xx( *((TVector*)fTrs->At(i)) );
23efe5f1 282 if(xx(0)>=2)
283 {
284// 1)Calculate Qtrm
e75c69cd 285 xx(6)=(this->Qtrm(i));
23efe5f1 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){
def162f4 297// cout<<xx(0)<<" ";
23efe5f1 298 for(Int_t j=1;j<11;j++){
38f7c306 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);}
23efe5f1 301 }
302 cout<<endl;
303 }
304// 4)Update data in TVector
e75c69cd 305 TClonesArray &arr=*fTrs;
23efe5f1 306 new(arr[i])TVector(xx);
307 }
308 else
309 {/*cout<<"No data for track "<<i<<endl;*/}
310}// End loop for tracks
311}
312void AliITSPid::Reset(void)
313{
e75c69cd 314 for(Int_t i=0;i<fTrs->GetEntries();i++){
23efe5f1 315 TVector xx(0,11);
e75c69cd 316 TClonesArray &arr=*fTrs;
23efe5f1 317 new(arr[i])TVector(xx);
318 }
319}
320//-----------------------------------------------------------
321AliITSPid::AliITSPid(Int_t ntrack)
322{
38f7c306 323 fSigmin=0.01;
e75c69cd 324 fTrs = new TClonesArray("TVector",ntrack);
325 TClonesArray &arr=*fTrs;
23efe5f1 326 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
e75c69cd 327 fMxtrs=0;
11b3ae13 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
339fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
e75c69cd 340 Double_t parpi0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
11b3ae13 341 1.1405329025e+05, -3.8289694496e+04 };
e75c69cd 342 fGGpi[0]->SetParameters(parpi0);
11b3ae13 343fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
e75c69cd 344 Double_t parpi1[10]={ 1.0791668283e-02, 9.7347716496e-01 };
345 fGGpi[1]->SetParameters(parpi1);
11b3ae13 346fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
e75c69cd 347 Double_t parpi2[10]={ 5.8191602279e-04, 9.7285601334e-02 };
348 fGGpi[2]->SetParameters(parpi2);
11b3ae13 349fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
e75c69cd 350 Double_t parpi3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
11b3ae13 351 6.2900977606e+03, -2.2935862292e+03 };
e75c69cd 352 fGGpi[3]->SetParameters(parpi3);
11b3ae13 353fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
e75c69cd 354 Double_t parpi4[10]={ 9.0419011783e-03, 1.1628922525e+00 };
355 fGGpi[4]->SetParameters(parpi4);
11b3ae13 356fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
e75c69cd 357 Double_t parpi5[10]={ 1.8324872519e-03, 2.1503968838e-01 };
358 fGGpi[5]->SetParameters(parpi5);
11b3ae13 359}//End Pions
360{//Kaons
361fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
e75c69cd 362 Double_t parka0[20]={
11b3ae13 363 -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
364 -3.7123338009e+03,1.6003647641e+03 };
e75c69cd 365 fGGka[0]->SetParameters(parka0);
11b3ae13 366fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
e75c69cd 367 Double_t parka1[20]={
11b3ae13 368 2.5181172905e-01,8.7566001814e-01 };
e75c69cd 369 fGGka[1]->SetParameters(parka1);
11b3ae13 370fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
e75c69cd 371 Double_t parka2[20]={
11b3ae13 372 8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
373 -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
374 5.7545491696e+01 };
e75c69cd 375 fGGka[2]->SetParameters(parka2);
11b3ae13 376}//End Kaons
377{//Protons
378fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
e75c69cd 379 Double_t parpr0[10]={
11b3ae13 380 6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
381 -3.5269200901e+03,1.2859128345e+03 };
e75c69cd 382 fGGpr[0]->SetParameters(parpr0);
11b3ae13 383fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
e75c69cd 384 Double_t parpr1[10]={
11b3ae13 385 9.4970837607e-01,7.3573504201e-01 };
e75c69cd 386 fGGpr[1]->SetParameters(parpr1);
11b3ae13 387fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
e75c69cd 388 Double_t parpr2[10]={
11b3ae13 389 1.2498403757e-01,2.7845072306e-02 };
e75c69cd 390 fGGpr[2]->SetParameters(parpr2);
11b3ae13 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 //-------------------------------------------------
e75c69cd 399const int kInf=10;
23efe5f1 400// Ncut Pmom pilo pihi klo khi plo phi
401// cut[j] [0] [1] [2] [3] [4] [5] [6]
402//----------------------------------------------------------------
e75c69cd 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 );
23efe5f1 407//----------------------------------------------------------------
e75c69cd 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
23efe5f1 410//----------------------------------------------------------------
e75c69cd 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
23efe5f1 414//----------------------------------------------------------------
38f7c306 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 ---------------------
e75c69cd 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
38f7c306 421 //---------------------------------------------------
e75c69cd 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;
38f7c306 425 //------------------------ pi,K,p -------------------
e75c69cd 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;
23efe5f1 430}
11b3ae13 431//End AliITSPid.cxx