]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TPC/AliTPCPid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
adding contructors and assignment operator for better usability; removing deprecated...
[u/mrichter/AliRoot.git] / TPC / AliTPCPid.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18#include "AliTPCPid.h"
19#include "TMath.h"
20//#include <TVector.h>
21#include <TF1.h>
22#include <TClonesArray.h>
23//#include "AliITSIOTrack.h"
24#include "AliKalmanTrack.h"
25#include "Riostream.h"
26#include <cassert>
27
28ClassImp(AliTPCPid)
29// Correction 13.01.2003 Z.S.,Dubna
30// 22.01.2003
31//------------------------------------------------------------
32// pid class by B. Batyunya
33// stupid corrections by M.K.
34//
35//_______________________________________________________
36//________________________________________________________
37AliTPCPid::AliTPCPid( const AliTPCPid& r):TObject(r),
38 fCutKa(0),
39 fCutPr(0),
40 fCutKaTune(0.),
41 fCutPrTune(0.),
42 fSigmin(0.),
43 fSilent(0),
44 fmxtrs(0),
45 trs(0),
46 fqtot(0.),
47 fWpi(0.),
48 fWk(0.),
49 fWp(0.),
50 fRpik(0.),
51 fRppi(0.),
52 fRpka(0.),
53 fRp(0.),
54 fPcode(0)
55{
56 // dummy copy constructor
57}
58Float_t AliTPCPid::Qcorr(Float_t xc)
59{
60 //
61 // charge correction
62 //
63 assert(0);
64 Float_t fcorr;
65 fcorr=( 0.766 +0.9692*xc -1.267*xc*xc )*( 1.-TMath::Exp(-xc*64.75) );
66 if(fcorr<=0.1)fcorr=0.1;
67return fqtot/fcorr;
68}
69//__________________________________________________________
70AliTPCPid & AliTPCPid::operator =(const AliTPCPid & param)
71{
72 //
73 // assignment operator - dummy
74 //
75 if(this!=&param){
76 fSigmin=param.fSigmin;
77 }
78 return (*this);
79}
80//-----------------------------------------------------------
81Float_t AliTPCPid::Qtrm(Int_t track) const
82{
83 //
84 // dummy comment (Boris!!!)
85 //
86 TVector q(*( this->GetVec(track) ));
87 Int_t ml=(Int_t)q(0);
88 if(ml<1)return 0.;
89 if(ml>6)ml=6;
90 float vf[6];
91 Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
92 if(nl==0)return 0.;
93 switch(nl){
94 case 1:q(6)=q(1); break;
95 case 2:q(6)=(q(1)+q(2))/2.; break;
96 default:
97 for(int fi=0;fi<2;fi++){
98 Int_t swap;
99 do{ swap=0; float qmin=vf[fi];
100 for(int j=fi+1;j<nl;j++)
101 if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
102 }while(swap==1);
103 }
104 q(6)= (vf[0]+vf[1])/2.;
105 break;
106 }
107 for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
108 return (q(6));
109}// --- End AliTPCPid::Qtrm ---
110
111Float_t AliTPCPid::Qtrm(Float_t qarr[6],Int_t narr)
112{
113 //..................
114 Float_t q[6],qm,qmin;
115 Int_t nl,ml;
116 if(narr>0&&narr<7){ml=narr;}else{return 0;};
117 nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
118 if(nl==0)return 0.;
119 switch(nl){
120 case 1:qm=q[0]; break;
121 case 2:qm=(q[0]+q[1])/2.; break;
122 default:
123 Int_t swap;
124 for(int fi=0;fi<2;fi++){
125 do{ swap=0; qmin=q[fi];
126 for(int j=fi+1;j<nl;j++)
127 if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
128 }while(swap==1);
129 }
130 qm= (q[0]+q[1])/2.;
131 break;
132 }
133 return qm;
134}// --- End AliTPCPid::Qtrm ---
135
136Int_t AliTPCPid::Wpik(Int_t nc,Float_t q)
137{
138 //
139 // pi-k
140 //
141 Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
142 Float_t appi,apk;
143 qmpi =fcut[nc][1];
144 sigpi=fcut[nc][2];
145 qmk =fcut[nc][3];
146 sigk =fcut[nc][4];
147 appi = faprob[0][nc-5];
148 apk = faprob[1][nc-5];
149 if( !fSilent ){
150 cout<<"qmpi,sigpi,qmk,sigk="<<qmpi<<" "<<sigpi<<" "<<qmk<<" "<<sigk<<endl;
151 cout<<"appi,apk="<<appi<<","<<apk<<endl;
152 }
153 Float_t dqpi=(q-qmpi)/sigpi;
154 Float_t dqk =(q-qmk )/sigk;
155 dpi =TMath::Abs(dqpi);
156 dk =TMath::Abs(dqk);
157 Double_t dn=appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk);
158 if(dn>0.){
159 ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
160 pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
161 }else{fWpi=1;return Pion();}
162 Float_t rpik=ppi/(pk+0.0000001);
163 if( !fSilent )
164 cout<<"q,dqpi,dqk, Wpik: ppi,pk,rpik="
165 <<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
166
167 fWp=0.; fWpi=ppi; fWk=pk;
168 if( pk>ppi){return Kaon();}else{return Pion();}
169}
170//-----------------------------------------------------------
171Int_t AliTPCPid::Wpikp(Int_t nc,Float_t q)
172{
173 //
174 // pi-k-p
175 //
176 Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
177 Float_t appi,apk,app;
178
179 qmpi =fcut[nc][1];
180 sigpi=fcut[nc][2];
181 qmk =fcut[nc][3];
182 sigk =fcut[nc][4];
183 qmp =fcut[nc][5];
184 sigp =fcut[nc][6];
185
186 //appi = apk = app = 1.;
187 appi = faprob[0][nc-5];
188 apk = faprob[1][nc-5];
189 app = faprob[2][nc-5];
190
191 Double_t dn= appi*TMath::Gaus(q,qmpi,sigpi)
192 +apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp);
193 if(dn>0.){
194 ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
195 pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
196 pp = app*TMath::Gaus(q,qmp, sigp )/dn;
197 }
198 else{fWpi=1;return Pion();}
199 fWp=pp; fWpi=ppi; fWk=pk;
200 if( !fSilent ){
201cout<<" Wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<"; "<<qmk<<" "<<sigk<<"; "
202 <<qmp<<" "<<sigp<<"; "<<endl;
203cout<<" faprob: "<<appi<<" "<<apk<<" "<<app<<endl;
204cout<<" ppi,pk,pp="<<ppi<<" "<<pk<<" "<<pp<<endl;
205 }
206 if( ppi>pk&&ppi>pp ) { return Pion(); }
207 if(pk>pp){return Kaon();}else{return Proton();}
208}
209//-----------------------------------------------------------
210Int_t AliTPCPid::GetPcode(TClonesArray* /*rps*/,Float_t /*pm*/) const
211{
212 //
213 // dummy ???
214 //
215 return 0;
216}
217//-----------------------------------------------------------
218Int_t AliTPCPid::GetPcode(AliTPCtrack *track)
219{
220 //
221 // get particle code
222 //
223 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
224 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
225 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
226 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
227 Float_t lam=TMath::ATan(par[3]);
228 Float_t pt1=TMath::Abs(par[4]);
229 Float_t mom=1./(pt1*TMath::Cos(lam));
230 Float_t dedx=track->GetdEdx();
231 Int_t pcode=GetPcode(dedx/40.,mom);
232 cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
233 return pcode?pcode:211;
234 }
235/*
236//-----------------------------------------------------------
237Int_t AliTPCPid::GetPcode(AliKalmanTrack *track)
238{
239 //
240 // get particle code
241 //
242 if(track==0)return 0;
243 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
244 track->PropagateTo(3.,0.0028,65.19);
245 track->PropagateToVertex();
246 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
247 Float_t lam=TMath::ATan(par[3]);
248 Float_t pt1=TMath::Abs(par[4]);
249 Float_t mom=0.;
250 if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
251 Float_t dedx=track->GetPIDsignal();
252 cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
253 Int_t pcode=GetPcode(dedx,mom);
254 cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
255return pcode?pcode:211;
256}
257*/
258//-----------------------------------------------------------
259Int_t AliTPCPid::GetPcode(Float_t q,Float_t pm)
260{
261 //
262 // get particle code
263 //
264 fWpi=fWk=fWp=0.; fPcode=0;
265//1)---------------------- 0-120 MeV/c --------------
266 if ( pm<=fcut[1][0] )
267 { fWpi=1.; return Pion(); }
268//2)----------------------120-200 Mev/c ( 29.7.2002 ) -------------
269 if ( pm<=fcut[2][0] )
270 { if( q<fCutKa->Eval(pm) ){fWpi=1.; return Pion(); }
271 else {fWk =1.; return Kaon();} }
272//3)----------------------200-400 Mev/c -------------
273 if ( pm<=fcut[3][0] ){
274 if( q<fCutKa->Eval(pm) )
275 { fWpi=1.;return Pion(); }
276 else
277 { if ( q<=fCutPr->Eval(pm) )
278 {fWk=1.;return Kaon();}
279 else {fWp=1.;return Proton();}
280 }}
281//4)----------------------400-450 Mev/c -------------
282 if ( pm<=fcut[4][0] ){
283 if( q<fCutKaTune*fCutKa->Eval(pm) )
284 { fWpi=1.;return Pion(); }
285 else
286 { if( q<fCutPr->Eval(pm) )
287 {fWk=1.;return Kaon();} else {fWp=1.;return Proton(); }
288 }}
289//5)----------------------450-500 Mev/c -------------
290 if ( pm<=fcut[5][0] ){
291 if ( q>fCutPr->Eval(pm) )
292 {fWp=1.;return Proton();} else {return Wpik(5,q);}}
293//6)----------------------500-550 Mev/c -------------
294 if ( pm<=fcut[6][0] ){
295 if ( q>fCutPr->Eval(pm) )
296 {fWp=1.;return Proton();} else {return Wpik(6,q);}}
297//7)----------------------550-600 Mev/c -------------
298 if ( pm<=fcut[7][0] ){
299 if ( q>fCutPr->Eval(pm) )
300 {fWp=1.;return Proton();} else {return Wpik(7,q);}}
301//8)----------------------600-650 Mev/c -------------
302 if ( pm<=fcut[8][0] ){
303 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
304 else {return Wpik(8,q);};
305//9)----------------------650-730 Mev/c -------------
306 if ( pm<=fcut[9][0] ){
307 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}
308 else {return Wpik(9,q);}}
309//10)----------------------730-830 Mev/c -------------
310 if( pm<=fcut[10][0] ){
311 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
312 else {return Wpik(10,q);};
313//11)----------------------830-930 Mev/c -------------
314 if( pm<=fcut[11][0] ){ return Wpikp(11,q); }
315//12)----------------------930-1030 Mev/c -------------
316 if( pm<=fcut[12][0] )
317 { return Wpikp(12,q); };
318
319 return fPcode;
320}
321//-----------------------------------------------------------
322void AliTPCPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
323 Float_t klo,Float_t khi,Float_t plo,Float_t phi)
324{
325 //
326 // set cuts
327 //
328 fcut[n][0]=pm;
329 fcut[n][1]=pilo;
330 fcut[n][2]=pihi;
331 fcut[n][3]=klo;
332 fcut[n][4]=khi;
333 fcut[n][5]=plo;
334 fcut[n][6]=phi;
335 return ;
336}
337//------------------------------------------------------------
338void AliTPCPid::SetVec(Int_t ntrack,TVector info) const
339{
340 //
341 // new track vector
342 //
343TClonesArray& arr=*trs;
344 new( arr[ntrack] ) TVector(info);
345}
346//-----------------------------------------------------------
347TVector* AliTPCPid::GetVec(Int_t ntrack) const
348{
349 //
350 // get track vector
351 //
352TClonesArray& arr=*trs;
353 return (TVector*)arr[ntrack];
354}
355//-----------------------------------------------------------
356void AliTPCPid::SetEdep(Int_t track,Float_t Edep)
357{
358 //
359 // energy deposit
360 //
361 TVector xx(0,11);
362 if( ((TVector*)trs->At(track))->IsValid() )
363 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
364 Int_t j=(Int_t)xx(0); if(j>4)return;
365 xx(++j)=Edep;xx(0)=j;
366 TClonesArray &arr=*trs;
367 new(arr[track])TVector(xx);
368}
369//-----------------------------------------------------------
370void AliTPCPid::SetPmom(Int_t track,Float_t Pmom)
371{
372 //
373 // set part. momentum
374 //
375 TVector xx(0,11);
376 if( ((TVector*)trs->At(track))->IsValid() )
377 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
378 xx(10)=Pmom;
379 TClonesArray &arr=*trs;
380 new(arr[track])TVector(xx);
381}
382//-----------------------------------------------------------
383void AliTPCPid::SetPcod(Int_t track,Int_t partcode)
384{
385 //
386 // set part. code
387 //
388 TVector xx(0,11);
389 if( ((TVector*)trs->At(track))->IsValid() )
390 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
391 if(xx(11)==0)
392 {xx(11)=partcode; fmxtrs++;
393 TClonesArray &arr=*trs;
394 new(arr[track])TVector(xx);
395 }
396}
397//-----------------------------------------------------------
398void AliTPCPid::PrintPID(Int_t track)
399{
400 //
401 // control print
402 //
403cout<<fmxtrs<<" tracks in AliITSPid obj."<<endl;
404 if( ((TVector*)trs->At(track))->IsValid() )
405 {TVector xx( *((TVector*)trs->At(track)) );
406 xx.Print();
407 }
408 else
409 {cout<<"No data for track "<<track<<endl;return;}
410}
411//-----------------------------------------------------------
412void AliTPCPid::Tab(void)
413{
414 //
415 // fill table
416 //
417if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
418cout<<"------------------------------------------------------------------------"<<endl;
419cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
420 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
421cout<<"------------------------------------------------------------------------"<<endl;
422for(Int_t i=0;i<trs->GetEntries();i++)
423{
424 TVector xx( *((TVector*)trs->At(i)) );
425 if( xx.IsValid() && xx(0)>0 )
426 {
427 if(xx(0)>=2)
428 {
429// 1)Calculate Qtrm
430 xx(6)=(this->Qtrm(i));
431
432 }else{
433 xx(6)=xx(1);
434 }
435// 2)Calculate Wpi,Wk,Wp
436 this->GetPcode(xx(6),xx(10)/1000.);
437 xx(7)=GetWpi();
438 xx(8)=GetWk();
439 xx(9)=GetWp();
440// 3)Print table
441 if(xx(0)>0){
442 cout<<xx(0)<<" ";
443 for(Int_t j=1;j<11;j++){
444 cout.width(7);cout.precision(5);cout<<xx(j);
445 }
446 cout<<endl;
447 }
448// 4)Update data in TVector
449 TClonesArray &arr=*trs;
450 new(arr[i])TVector(xx);
451 }
452 else
453 {/*cout<<"No data for track "<<i<<endl;*/}
454}// End loop for tracks
455}
456void AliTPCPid::Reset(void)
457{
458 //
459 // reset
460 //
461 for(Int_t i=0;i<trs->GetEntries();i++){
462 TVector xx(0,11);
463 TClonesArray &arr=*trs;
464 new(arr[i])TVector(xx);
465 }
466}
467//-----------------------------------------------------------
468AliTPCPid::AliTPCPid(Int_t ntrack):TObject(),
469 fCutKa(0),
470 fCutPr(0),
471 fCutKaTune(0.),
472 fCutPrTune(0.),
473 fSigmin(0.),
474 fSilent(0),
475 fmxtrs(0),
476 trs(0),
477 fqtot(0.),
478 fWpi(0.),
479 fWk(0.),
480 fWp(0.),
481 fRpik(0.),
482 fRppi(0.),
483 fRpka(0.),
484 fRp(0.),
485 fPcode(0)
486{
487 trs = new TClonesArray("TVector",ntrack);
488 TClonesArray &arr=*trs;
489 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
490
491
492 //fCutKa = new TF1("fkaons","[0]/x/x+[1]",0.1,1.2);
493 //fCutPr = new TF1("fprotons","[0]/x/x +[1]",0.2,1.2);
494 //TF1 *frmska=0;
495
496 TF1 frmska("x_frmska","1.46-7.82*x+16.78*x^2-15.53*x^3+5.24*x^4 ",
497 0.1,1.2);
498 fCutKa = new TF1("fkaons",
499 "1.25+0.044/x/x+1.25+0.044*x-13.87*x^2+22.37*x^3-10.05*x^4-2.5*x_frmska",
500 0.1,1.2);
501 fCutPr = new TF1("fprotons",
502 "0.83*(15.32-56.094*x+89.962*x^2-66.1856*x^3+18.4052*x^4)",
503 0.2,1.2);
504 fCutKaTune=1.1; // 0.92;
505 fCutPrTune=1.0; //0.80;
506
507const int kinf=10;
508// Ncut Pmom pilo pihi klo khi plo phi
509// fcut[j] [0] [1] [2] [3] [4] [5] [6]
510//----------------------------------------------------------------
511 SetCut( 1, 0.120, 0. , 0. , kinf , kinf , kinf , kinf );
512 SetCut( 2, 0.200, 0. , 6.0 , 6.0 , kinf , kinf , kinf ); //120-200
513 SetCut( 3, 0.400, 0. , 3.5 , 3.5 , 9.0 , 9.0 , kinf ); //200-400
514 SetCut( 4, 0.450, 0. , 1.9 , 1.9 , 4.0 , 4.0 , kinf ); //400-450
515//----------------------------------------------------------------
516 SetCut( 5, 0.500, 0.976, 0.108, 1.484 , 0.159 , 3.5 , kinf ); //450-500
517 SetCut( 6, 0.550, 0.979, 0.108, 1.376 , 0.145 , 3.0 , kinf ); //500-550
518//----------------------------------------------------------------
519 SetCut( 7, 0.600, 0.984, 0.111, 1.295 , 0.146 , 2.7 , kinf ); //550-600
520 SetCut( 8, 0.650, 0.989, 0.113, 1.239 , 0.141 , 2.5 , kinf ); //600-650
521 SetCut( 9, 0.730, 0.995, 0.109, 1.172 , 0.132 , 2.0 , kinf ); //650-730
522//----------------------------------------------------------------
523 SetCut( 10, 0.830, 1.008, 0.116, 1.117 , 0.134 , 1.703, 0.209 ); //730-830
524 SetCut( 11, 0.930, 1.019, 0.115, 1.072 , 0.121 , 1.535, 0.215 ); //830-930
525 SetCut( 12, 1.230, 1.035, 0.117, 1.053 , 0.140 ,1.426, 0.270); //930-1030
526 //------------------------ pi,K ---------------------
527 faprob[0][0]=33219.; faprob[1][0]=1971.; // faprob[0][i] - pions
528 faprob[0][1]=28828.; faprob[1][1]=1973.; // faprob[1][i] - kaons
529 //---------------------------------------------------
530 faprob[0][2]=24532.; faprob[1][2]=1932.; faprob[2][2]=1948.;
531 faprob[0][3]=20797.; faprob[1][3]=1823.; faprob[2][3]=1970.;
532 faprob[0][4]=27017.; faprob[1][4]=2681.; faprob[2][4]=2905.;
533 //------------------------ pi,K,p -------------------
534 faprob[0][5]= 24563.; faprob[1][5]=2816.; faprob[2][5]=3219.;
535 faprob[0][6]= 16877.; faprob[1][6]=2231.; faprob[2][6]=2746.;
536 faprob[0][7]= 11557.; faprob[1][7]=1681; faprob[2][7]=2190.;
537
538 fSilent=kTRUE;
539}
540//-----------------------------------------------------------
541
542
543