]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSPid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New methods and data member added by M. Horner.
[u/mrichter/AliRoot.git] / ITS / AliITSPid.cxx
CommitLineData
23efe5f1 1#include "AliITSPid.h"
2#include "TMath.h"
79b921e1 3#include "AliITSIOTrack.h"
4ae5bbc4 4#include <Riostream.h>
23efe5f1 5ClassImp(AliITSPid)
6
23efe5f1 7Float_t AliITSPid::qcorr(Float_t xc)
8{
9 assert(0);
10 Float_t fcorr;
11 fcorr=( 0.766 +0.9692*xc -1.267*xc*xc )*( 1.-TMath::Exp(-xc*64.75) );
12 if(fcorr<=0.1)fcorr=0.1;
13return qtot/fcorr;
14}
38f7c306 15
23efe5f1 16Float_t AliITSPid::qtrm(Int_t track)
17{
18 TVector q(*( this->GetVec(track) ));
38f7c306 19 Int_t ml=(Int_t)q(0);
20 if(ml<1)return 0.;
21 if(ml>6)ml=6;
22 float vf[6];
23 Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
24 if(nl==0)return 0.;
23efe5f1 25 switch(nl){
26 case 1:q(6)=q(1); break;
27 case 2:q(6)=(q(1)+q(2))/2.; break;
28 default:
38f7c306 29 for(int fi=0;fi<2;fi++){
30 Int_t swap;
31 do{ swap=0; float qmin=vf[fi];
32 for(int j=fi+1;j<nl;j++)
33 if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
34 }while(swap==1);
35 }
23efe5f1 36 q(6)= (vf[0]+vf[1])/2.;
37 break;
38 }
39 for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
40 return (q(6));
41}
38f7c306 42
79b921e1 43Float_t AliITSPid::qtrm(Float_t qarr[6],Int_t narr)
44{
38f7c306 45 Float_t q[6],qm,qmin;
79b921e1 46 Int_t nl,ml;
38f7c306 47 if(narr>0&&narr<7){ml=narr;}else{return 0;};
48 nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
79b921e1 49 if(nl==0)return 0.;
79b921e1 50 switch(nl){
51 case 1:qm=q[0]; break;
52 case 2:qm=(q[0]+q[1])/2.; break;
53 default:
38f7c306 54 Int_t swap;
55 for(int fi=0;fi<2;fi++){
56 do{ swap=0; qmin=q[fi];
57 for(int j=fi+1;j<nl;j++)
58 if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
59 }while(swap==1);
60 }
61 qm= (q[0]+q[1])/2.;
79b921e1 62 break;
63 }
64 return qm;
65}
23efe5f1 66//-----------------------------------------------------------
79b921e1 67//inline
23efe5f1 68Int_t AliITSPid::wpik(Int_t nc,Float_t q)
69{
38f7c306 70 // return pion();
23efe5f1 71
72 Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
38f7c306 73 Float_t appi,apk;
23efe5f1 74 qmpi =cut[nc][1];
75 sigpi=cut[nc][2];
76 qmk =cut[nc][3];
77 sigk =cut[nc][4];
38f7c306 78 appi = aprob[0][nc-5];
79 apk = aprob[1][nc-5];
def162f4 80// cout<<"qmpi,sigpi,qmk,sigk="<<qmpi<<" "<<sigpi<<" "<<qmk<<" "<<sigk<<endl;
81// cout<<"appi,apk="<<appi<<","<<apk<<endl;
23efe5f1 82 Float_t dqpi=(q-qmpi)/sigpi;
83 Float_t dqk =(q-qmk )/sigk;
84 if( dqk<-1. )return pion();
6311f063 85 dpi =TMath::Abs(dqpi);
86 dk =TMath::Abs(dqk);
38f7c306 87 //ppi =1.- TMath::Erf(dpi); // +0.5;
88 //pk =1.- TMath::Erf(dk); // +0.5;
89 ppi=appi*TMath::Gaus(q,qmpi,sigpi)
90 /(appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk));
91 pk = apk*TMath::Gaus(q,qmk, sigk )
92 /(appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk));
93
def162f4 94// Float_t rpik=ppi/(pk+0.0000001);
95// cout<<"q,dqpi,dqk, wpik: ppi,pk,rpik="
96// <<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
38f7c306 97
98 fWp=0.; fWpi=ppi; fWk=pk;
99 if( pk>ppi){return kaon();}else{return pion();}
23efe5f1 100}
101//-----------------------------------------------------------
79b921e1 102//inline
23efe5f1 103Int_t AliITSPid::wpikp(Int_t nc,Float_t q)
104{
38f7c306 105 Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
106 Float_t appi,apk,app;
23efe5f1 107
108 qmpi =cut[nc][1];
109 sigpi=cut[nc][2];
110 qmk =cut[nc][3];
111 sigk =cut[nc][4];
112 qmp =cut[nc][5];
113 sigp =cut[nc][6];
114
38f7c306 115 //appi = apk = app = 1.;
116 appi = aprob[0][nc-5];
117 apk = aprob[1][nc-5];
118 app = aprob[2][nc-5];
23efe5f1 119
38f7c306 120 //ppi =TMath::Erf( TMath::Abs((q-qmpi)/sigpi) )+0.5;
121 //pka =TMath::Erf( TMath::Abs((q-qmk )/sigk) )+0.5;
122 //ppr =TMath::Erf( TMath::Abs((q-qmp )/sigp) )+0.5;
123
124 ppi=appi*TMath::Gaus(q,qmpi,sigpi)
125 /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
126 pk = apk*TMath::Gaus(q,qmk, sigk )
127 /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
128 pp = app*TMath::Gaus(q,qmp, sigp )
129 /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
130
131 fWp=pp; fWpi=ppi; fWk=pk;
132
def162f4 133//cout<<" wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<"; "<<qmk<<" "<<sigk<<"; "
134// <<qmp<<" "<<sigp<<"; "<<endl;
135// cout<<" aprob: "<<appi<<" "<<apk<<" "<<app<<endl;
136//cout<<" ppi,pk,pp="<<ppi<<" "<<pk<<" "<<pp<<endl;
38f7c306 137
138 if( ppi>pk&&ppi>pp ) { return pion(); }
139 if(pk>pp){return kaon();}else{return proton();}
23efe5f1 140}
141//-----------------------------------------------------------
142Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
143{
144 return 0;
145}
146//-----------------------------------------------------------
79b921e1 147Int_t AliITSPid::GetPcode(AliTPCtrack *track)
148{
149 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
150 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
151 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
152 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
153 Float_t lam=TMath::ATan(par[3]);
154 Float_t pt_1=TMath::Abs(par[4]);
155 Float_t mom=1./(pt_1*TMath::Cos(lam));
156 Float_t dedx=track->GetdEdx();
157 Int_t pcode=GetPcode(dedx/40.,mom);
def162f4 158// cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 159 return pcode?pcode:211;
160 }
161//------------------------------------------------------------
79b921e1 162Int_t AliITSPid::GetPcode(AliITSIOTrack *track)
163{
164 Double_t px,py,pz;
165 px=track->GetPx();
166 py=track->GetPy();
167 pz=track->GetPz();
168 Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
169//???????????????????
def162f4 170 // Float_t dedx=1.0;
171 Float_t dedx=track->GetdEdx();
79b921e1 172//???????????????????
173 Int_t pcode=GetPcode(dedx,mom);
def162f4 174// cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 175return pcode?pcode:211;
176}
79b921e1 177//-----------------------------------------------------------
178Int_t AliITSPid::GetPcode(AliITStrackV2 *track)
179{
180 if(track==0)return 0;
38f7c306 181 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
79b921e1 182 track->PropagateTo(3.,0.0028,65.19);
79b921e1 183 track->PropagateToVertex();
184 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
185 Float_t lam=TMath::ATan(par[3]);
186 Float_t pt_1=TMath::Abs(par[4]);
187 Float_t mom=0.;
188 if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
189 Float_t dedx=track->GetdEdx();
def162f4 190// cout<<"lam,pt_1,mom,dedx="<<lam<<","<<pt_1<<","<<mom<<","<<dedx<<endl;
79b921e1 191 Int_t pcode=GetPcode(dedx,mom);
def162f4 192// cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 193return pcode?pcode:211;
194}
195//-----------------------------------------------------------
23efe5f1 196Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
197{
38f7c306 198 fWpi=fWk=fWp=0.; fPcode=0;
199//1)---------------------- 0-120 MeV/c --------------
23efe5f1 200 if ( pm<=cut[1][0] )
201 { return pion(); }
38f7c306 202//2)----------------------120-200 Mev/c -------------
23efe5f1 203 if ( pm<=cut[2][0] )
204 { if( q<cut[2][2] ){ return pion(); } else { return kaon();} }
38f7c306 205//3)----------------------200-300 Mev/c -------------
23efe5f1 206 if ( pm<=cut[3][0] )
207 if( q<cut[3][2])
208 { return pion(); }
209 else
210 { if ( q<=cut[3][5] ) {return kaon();} else {return proton();}}
38f7c306 211//4)----------------------300-410 Mev/c -------------
23efe5f1 212 if ( pm<=cut[4][0] )
213 if( q<cut[4][2] )
214 { return pion(); }
215 else
216 { if( q<=cut[4][4] ) {return kaon();} else {return proton(); }}
38f7c306 217//5)----------------------410-470 Mev/c -------------
23efe5f1 218 if ( pm<=cut[5][0] )
219 if ( q>cut[5][5] ) {return proton();} else {return wpik(5,q);};
38f7c306 220//6)----------------------470-530 Mev/c -------------
23efe5f1 221 if ( pm<=cut[6][0] )
222 if ( q>cut[6][5] ) {return proton();} else {return wpik(6,q);};
38f7c306 223//7)----------------------530-590 Mev/c -------------
23efe5f1 224 if ( pm<=cut[7][0] )
38f7c306 225 if ( q<=cut[7][5] ) {return wpik(7,q);} else {return proton();};
226//8)----------------------590-650 Mev/c -------------
23efe5f1 227 if ( pm<=cut[8][0] )
38f7c306 228 if ( q<=cut[8][5] ) {return wpik(8,q);} else {return proton();};
229//9)----------------------650-730 Mev/c -------------
23efe5f1 230 if ( pm<=cut[9][0] )
38f7c306 231 if ( q<=cut[9][5] ) {return wpik(9,q);} else {return proton();};
232//10)----------------------730-830 Mev/c -------------
23efe5f1 233 if( pm<=cut[10][0] ){ return wpikp(10,q); }
38f7c306 234//11)----------------------830-930 Mev/c -------------
23efe5f1 235 if( pm<=cut[11][0] ){ return wpikp(11,q); }
38f7c306 236//12)----------------------930-1030 Mev/c -------------
23efe5f1 237 if( pm<=cut[12][0] ){ return wpikp(12,q); }
238
239 return fPcode;
240}
241//-----------------------------------------------------------
242void AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
243 Float_t klo,Float_t khi,Float_t plo,Float_t phi)
244{
245 cut[n][0]=pm;
246 cut[n][1]=pilo;
247 cut[n][2]=pihi;
248 cut[n][3]=klo;
249 cut[n][4]=khi;
250 cut[n][5]=plo;
251 cut[n][6]=phi;
252 return ;
253}
38f7c306 254//------------------------------------------------------------
23efe5f1 255void AliITSPid::SetVec(Int_t ntrack,TVector info)
256{
257TClonesArray& arr=*trs;
258 new( arr[ntrack] ) TVector(info);
259}
260//-----------------------------------------------------------
261TVector* AliITSPid::GetVec(Int_t ntrack)
262{
263TClonesArray& arr=*trs;
264 return (TVector*)arr[ntrack];
265}
266//-----------------------------------------------------------
267void AliITSPid::SetEdep(Int_t track,Float_t Edep)
268{
269 TVector xx(0,11);
270 if( ((TVector*)trs->At(track))->IsValid() )
271 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
272 Int_t j=(Int_t)xx(0); if(j>4)return;
273 xx(++j)=Edep;xx(0)=j;
274 TClonesArray &arr=*trs;
275 new(arr[track])TVector(xx);
276}
277//-----------------------------------------------------------
278void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
279{
280 TVector xx(0,11);
281 if( ((TVector*)trs->At(track))->IsValid() )
282 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
283 xx(10)=Pmom;
284 TClonesArray &arr=*trs;
285 new(arr[track])TVector(xx);
286}
287//-----------------------------------------------------------
288void AliITSPid::SetPcod(Int_t track,Int_t partcode)
289{
290 TVector xx(0,11);
291 if( ((TVector*)trs->At(track))->IsValid() )
292 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
293 if(xx(11)==0)
294 {xx(11)=partcode; mxtrs++;
295 TClonesArray &arr=*trs;
296 new(arr[track])TVector(xx);
297 }
298}
299//-----------------------------------------------------------
300void AliITSPid::Print(Int_t track)
301{cout<<mxtrs<<" tracks in AliITSPid obj."<<endl;
302 if( ((TVector*)trs->At(track))->IsValid() )
303 {TVector xx( *((TVector*)trs->At(track)) );
304 xx.Print();
305 }
306 else
307 {cout<<"No data for track "<<track<<endl;return;}
308}
309//-----------------------------------------------------------
310void AliITSPid::Tab(void)
311{
312if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
313cout<<"------------------------------------------------------------------------"<<endl;
314cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
315 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
316cout<<"------------------------------------------------------------------------"<<endl;
317for(Int_t i=0;i<trs->GetEntries();i++)
318{
319 TVector xx( *((TVector*)trs->At(i)) );
320 if( xx.IsValid() && xx(0)>0 )
321 {
322 TVector xx( *((TVector*)trs->At(i)) );
323 if(xx(0)>=2)
324 {
325// 1)Calculate Qtrm
326 xx(6)=(this->qtrm(i));
327
328 }else{
329 xx(6)=xx(1);
330 }
331// 2)Calculate Wpi,Wk,Wp
332 this->GetPcode(xx(6),xx(10)/1000.);
333 xx(7)=GetWpi();
334 xx(8)=GetWk();
335 xx(9)=GetWp();
336// 3)Print table
337 if(xx(0)>0){
def162f4 338// cout<<xx(0)<<" ";
23efe5f1 339 for(Int_t j=1;j<11;j++){
38f7c306 340 if(i<7){ cout.width(7);cout.precision(4);cout<<xx(j);}
341 if(i>7){ cout.width(7);cout.precision(5);cout<<xx(j);}
23efe5f1 342 }
343 cout<<endl;
344 }
345// 4)Update data in TVector
346 TClonesArray &arr=*trs;
347 new(arr[i])TVector(xx);
348 }
349 else
350 {/*cout<<"No data for track "<<i<<endl;*/}
351}// End loop for tracks
352}
353void AliITSPid::Reset(void)
354{
355 for(Int_t i=0;i<trs->GetEntries();i++){
356 TVector xx(0,11);
357 TClonesArray &arr=*trs;
358 new(arr[i])TVector(xx);
359 }
360}
361//-----------------------------------------------------------
362AliITSPid::AliITSPid(Int_t ntrack)
363{
38f7c306 364 fSigmin=0.01;
23efe5f1 365 trs = new TClonesArray("TVector",ntrack);
366 TClonesArray &arr=*trs;
367 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
368 mxtrs=0;
369const int inf=10;
370// Ncut Pmom pilo pihi klo khi plo phi
371// cut[j] [0] [1] [2] [3] [4] [5] [6]
372//----------------------------------------------------------------
373 SetCut( 1, 0.12 , 0. , 0. , inf , inf , inf , inf );
374 SetCut( 2, 0.20 , 0. , 6.0 , 6.0 , inf , inf , inf );
375 SetCut( 3, 0.30 , 0. , 3.5 , 3.5 , 9.0 , 9.0 , inf );
376 SetCut( 4, 0.41 , 0. , 1.9 , 1.9 , 4.0 , 4.0 , inf );
377//----------------------------------------------------------------
38f7c306 378 SetCut( 5, 0.47 , 0.935, 0.139, 1.738 , 0.498 , 3.5 , inf ); //410-470
379 SetCut( 6, 0.53 , 0.914, 0.136, 1.493 , 0.436 , 3.0 , inf ); //470-530
23efe5f1 380//----------------------------------------------------------------
38f7c306 381 SetCut( 7, 0.59 , 0.895, 0.131, 1.384 , 0.290 , 2.7 , inf ); //530-590
382 SetCut( 8, 0.65 , 0.887, 0.121, 1.167 , 0.287 , 2.5 , inf ); //590-650
383 SetCut( 9, 0.73 , 0.879, 0.120, 1.153 , 0.257 , 2.0 , inf ); //650-730
23efe5f1 384//----------------------------------------------------------------
38f7c306 385 SetCut( 10, 0.83 , 0.880, 0.126, 1.164 , 0.204 , 2.308 , 0.297 ); //730-830
386 SetCut( 11, 0.93 , 0.918, 0.145, 1.164 , 0.204 , 2.00 , 0.168 ); //830-930
387 SetCut( 12, 1.03 , 0.899, 0.128, 1.164 , 0.204 ,1.80 , 0.168);
388 //------------------------ pi,K ---------------------
389 aprob[0][0]=1212; aprob[1][0]=33.; // aprob[0][i] - const for pions,cut[i+5]
390 aprob[0][1]=1022; aprob[1][1]=46.2 ; // aprob[1][i] - kaons
391 //---------------------------------------------------
392 aprob[0][2]= 889.7; aprob[1][2]=66.58; aprob[2][2]=14.53;
393 aprob[0][3]= 686.; aprob[1][3]=88.8; aprob[2][3]=19.27;
394 aprob[0][4]= 697.; aprob[1][4]=125.6; aprob[2][4]=28.67;
395 //------------------------ pi,K,p -------------------
396 aprob[0][5]= 633.7; aprob[1][5]=100.1; aprob[2][5]=37.99; // aprob[2][i] - protons
397 aprob[0][6]= 469.5; aprob[1][6]=20.74; aprob[2][6]=25.43;
398 aprob[0][7]= 355.; aprob[1][7]=
399 355.*(20.74/469.5); aprob[2][7]=34.08;
23efe5f1 400}
401//-----------------------------------------------------------
402