]> 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
Script supersceeded by AliForwarddNdetaTask.C and
[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 fSigmin=param.fSigmin;
76 return (*this);
77}
78//-----------------------------------------------------------
79Float_t AliTPCPid::Qtrm(Int_t track) const
80{
81 //
82 // dummy comment (Boris!!!)
83 //
84 TVector q(*( this->GetVec(track) ));
85 Int_t ml=(Int_t)q(0);
86 if(ml<1)return 0.;
87 if(ml>6)ml=6;
88 float vf[6];
89 Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
90 if(nl==0)return 0.;
91 switch(nl){
92 case 1:q(6)=q(1); break;
93 case 2:q(6)=(q(1)+q(2))/2.; break;
94 default:
95 for(int fi=0;fi<2;fi++){
96 Int_t swap;
97 do{ swap=0; float qmin=vf[fi];
98 for(int j=fi+1;j<nl;j++)
99 if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
100 }while(swap==1);
101 }
102 q(6)= (vf[0]+vf[1])/2.;
103 break;
104 }
105 for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
106 return (q(6));
107}// --- End AliTPCPid::Qtrm ---
108
109Float_t AliTPCPid::Qtrm(Float_t qarr[6],Int_t narr)
110{
111 //..................
112 Float_t q[6],qm,qmin;
113 Int_t nl,ml;
114 if(narr>0&&narr<7){ml=narr;}else{return 0;};
115 nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
116 if(nl==0)return 0.;
117 switch(nl){
118 case 1:qm=q[0]; break;
119 case 2:qm=(q[0]+q[1])/2.; break;
120 default:
121 Int_t swap;
122 for(int fi=0;fi<2;fi++){
123 do{ swap=0; qmin=q[fi];
124 for(int j=fi+1;j<nl;j++)
125 if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
126 }while(swap==1);
127 }
128 qm= (q[0]+q[1])/2.;
129 break;
130 }
131 return qm;
132}// --- End AliTPCPid::Qtrm ---
133
134Int_t AliTPCPid::Wpik(Int_t nc,Float_t q)
135{
136 //
137 // pi-k
138 //
139 Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
140 Float_t appi,apk;
141 qmpi =fcut[nc][1];
142 sigpi=fcut[nc][2];
143 qmk =fcut[nc][3];
144 sigk =fcut[nc][4];
145 appi = faprob[0][nc-5];
146 apk = faprob[1][nc-5];
147 if( !fSilent ){
148 cout<<"qmpi,sigpi,qmk,sigk="<<qmpi<<" "<<sigpi<<" "<<qmk<<" "<<sigk<<endl;
149 cout<<"appi,apk="<<appi<<","<<apk<<endl;
150 }
151 Float_t dqpi=(q-qmpi)/sigpi;
152 Float_t dqk =(q-qmk )/sigk;
153 dpi =TMath::Abs(dqpi);
154 dk =TMath::Abs(dqk);
155 Double_t dn=appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk);
156 if(dn>0.){
157 ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
158 pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
159 }else{fWpi=1;return Pion();}
160 Float_t rpik=ppi/(pk+0.0000001);
161 if( !fSilent )
162 cout<<"q,dqpi,dqk, Wpik: ppi,pk,rpik="
163 <<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
164
165 fWp=0.; fWpi=ppi; fWk=pk;
166 if( pk>ppi){return Kaon();}else{return Pion();}
167}
168//-----------------------------------------------------------
169Int_t AliTPCPid::Wpikp(Int_t nc,Float_t q)
170{
171 //
172 // pi-k-p
173 //
174 Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
175 Float_t appi,apk,app;
176
177 qmpi =fcut[nc][1];
178 sigpi=fcut[nc][2];
179 qmk =fcut[nc][3];
180 sigk =fcut[nc][4];
181 qmp =fcut[nc][5];
182 sigp =fcut[nc][6];
183
184 //appi = apk = app = 1.;
185 appi = faprob[0][nc-5];
186 apk = faprob[1][nc-5];
187 app = faprob[2][nc-5];
188
189 Double_t dn= appi*TMath::Gaus(q,qmpi,sigpi)
190 +apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp);
191 if(dn>0.){
192 ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
193 pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
194 pp = app*TMath::Gaus(q,qmp, sigp )/dn;
195 }
196 else{fWpi=1;return Pion();}
197 fWp=pp; fWpi=ppi; fWk=pk;
198 if( !fSilent ){
199cout<<" Wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<"; "<<qmk<<" "<<sigk<<"; "
200 <<qmp<<" "<<sigp<<"; "<<endl;
201cout<<" faprob: "<<appi<<" "<<apk<<" "<<app<<endl;
202cout<<" ppi,pk,pp="<<ppi<<" "<<pk<<" "<<pp<<endl;
203 }
204 if( ppi>pk&&ppi>pp ) { return Pion(); }
205 if(pk>pp){return Kaon();}else{return Proton();}
206}
207//-----------------------------------------------------------
208Int_t AliTPCPid::GetPcode(TClonesArray* /*rps*/,Float_t /*pm*/) const
209{
210 //
211 // dummy ???
212 //
213 return 0;
214}
215//-----------------------------------------------------------
216Int_t AliTPCPid::GetPcode(AliTPCtrack *track)
217{
218 //
219 // get particle code
220 //
221 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
222 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
223 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
224 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
225 Float_t lam=TMath::ATan(par[3]);
226 Float_t pt1=TMath::Abs(par[4]);
227 Float_t mom=1./(pt1*TMath::Cos(lam));
228 Float_t dedx=track->GetdEdx();
229 Int_t pcode=GetPcode(dedx/40.,mom);
230 cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
231 return pcode?pcode:211;
232 }
233/*
234//-----------------------------------------------------------
235Int_t AliTPCPid::GetPcode(AliKalmanTrack *track)
236{
237 //
238 // get particle code
239 //
240 if(track==0)return 0;
241 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
242 track->PropagateTo(3.,0.0028,65.19);
243 track->PropagateToVertex();
244 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
245 Float_t lam=TMath::ATan(par[3]);
246 Float_t pt1=TMath::Abs(par[4]);
247 Float_t mom=0.;
248 if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
249 Float_t dedx=track->GetPIDsignal();
250 cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
251 Int_t pcode=GetPcode(dedx,mom);
252 cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
253return pcode?pcode:211;
254}
255*/
256//-----------------------------------------------------------
257Int_t AliTPCPid::GetPcode(Float_t q,Float_t pm)
258{
259 //
260 // get particle code
261 //
262 fWpi=fWk=fWp=0.; fPcode=0;
263//1)---------------------- 0-120 MeV/c --------------
264 if ( pm<=fcut[1][0] )
265 { fWpi=1.; return Pion(); }
266//2)----------------------120-200 Mev/c ( 29.7.2002 ) -------------
267 if ( pm<=fcut[2][0] )
268 { if( q<fCutKa->Eval(pm) ){fWpi=1.; return Pion(); }
269 else {fWk =1.; return Kaon();} }
270//3)----------------------200-400 Mev/c -------------
271 if ( pm<=fcut[3][0] ){
272 if( q<fCutKa->Eval(pm) )
273 { fWpi=1.;return Pion(); }
274 else
275 { if ( q<=fCutPr->Eval(pm) )
276 {fWk=1.;return Kaon();}
277 else {fWp=1.;return Proton();}
278 }}
279//4)----------------------400-450 Mev/c -------------
280 if ( pm<=fcut[4][0] ){
281 if( q<fCutKaTune*fCutKa->Eval(pm) )
282 { fWpi=1.;return Pion(); }
283 else
284 { if( q<fCutPr->Eval(pm) )
285 {fWk=1.;return Kaon();} else {fWp=1.;return Proton(); }
286 }}
287//5)----------------------450-500 Mev/c -------------
288 if ( pm<=fcut[5][0] ){
289 if ( q>fCutPr->Eval(pm) )
290 {fWp=1.;return Proton();} else {return Wpik(5,q);}}
291//6)----------------------500-550 Mev/c -------------
292 if ( pm<=fcut[6][0] ){
293 if ( q>fCutPr->Eval(pm) )
294 {fWp=1.;return Proton();} else {return Wpik(6,q);}}
295//7)----------------------550-600 Mev/c -------------
296 if ( pm<=fcut[7][0] ){
297 if ( q>fCutPr->Eval(pm) )
298 {fWp=1.;return Proton();} else {return Wpik(7,q);}}
299//8)----------------------600-650 Mev/c -------------
300 if ( pm<=fcut[8][0] ){
301 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
302 else {return Wpik(8,q);};
303//9)----------------------650-730 Mev/c -------------
304 if ( pm<=fcut[9][0] ){
305 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}
306 else {return Wpik(9,q);}}
307//10)----------------------730-830 Mev/c -------------
308 if( pm<=fcut[10][0] ){
309 if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
310 else {return Wpik(10,q);};
311//11)----------------------830-930 Mev/c -------------
312 if( pm<=fcut[11][0] ){ return Wpikp(11,q); }
313//12)----------------------930-1030 Mev/c -------------
314 if( pm<=fcut[12][0] )
315 { return Wpikp(12,q); };
316
317 return fPcode;
318}
319//-----------------------------------------------------------
320void AliTPCPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
321 Float_t klo,Float_t khi,Float_t plo,Float_t phi)
322{
323 //
324 // set cuts
325 //
326 fcut[n][0]=pm;
327 fcut[n][1]=pilo;
328 fcut[n][2]=pihi;
329 fcut[n][3]=klo;
330 fcut[n][4]=khi;
331 fcut[n][5]=plo;
332 fcut[n][6]=phi;
333 return ;
334}
335//------------------------------------------------------------
336void AliTPCPid::SetVec(Int_t ntrack,TVector info) const
337{
338 //
339 // new track vector
340 //
341TClonesArray& arr=*trs;
342 new( arr[ntrack] ) TVector(info);
343}
344//-----------------------------------------------------------
345TVector* AliTPCPid::GetVec(Int_t ntrack) const
346{
347 //
348 // get track vector
349 //
350TClonesArray& arr=*trs;
351 return (TVector*)arr[ntrack];
352}
353//-----------------------------------------------------------
354void AliTPCPid::SetEdep(Int_t track,Float_t Edep)
355{
356 //
357 // energy deposit
358 //
359 TVector xx(0,11);
360 if( ((TVector*)trs->At(track))->IsValid() )
361 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
362 Int_t j=(Int_t)xx(0); if(j>4)return;
363 xx(++j)=Edep;xx(0)=j;
364 TClonesArray &arr=*trs;
365 new(arr[track])TVector(xx);
366}
367//-----------------------------------------------------------
368void AliTPCPid::SetPmom(Int_t track,Float_t Pmom)
369{
370 //
371 // set part. momentum
372 //
373 TVector xx(0,11);
374 if( ((TVector*)trs->At(track))->IsValid() )
375 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
376 xx(10)=Pmom;
377 TClonesArray &arr=*trs;
378 new(arr[track])TVector(xx);
379}
380//-----------------------------------------------------------
381void AliTPCPid::SetPcod(Int_t track,Int_t partcode)
382{
383 //
384 // set part. code
385 //
386 TVector xx(0,11);
387 if( ((TVector*)trs->At(track))->IsValid() )
388 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
389 if(xx(11)==0)
390 {xx(11)=partcode; fmxtrs++;
391 TClonesArray &arr=*trs;
392 new(arr[track])TVector(xx);
393 }
394}
395//-----------------------------------------------------------
396void AliTPCPid::PrintPID(Int_t track)
397{
398 //
399 // control print
400 //
401cout<<fmxtrs<<" tracks in AliITSPid obj."<<endl;
402 if( ((TVector*)trs->At(track))->IsValid() )
403 {TVector xx( *((TVector*)trs->At(track)) );
404 xx.Print();
405 }
406 else
407 {cout<<"No data for track "<<track<<endl;return;}
408}
409//-----------------------------------------------------------
410void AliTPCPid::Tab(void)
411{
412 //
413 // fill table
414 //
415if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
416cout<<"------------------------------------------------------------------------"<<endl;
417cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
418 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
419cout<<"------------------------------------------------------------------------"<<endl;
420for(Int_t i=0;i<trs->GetEntries();i++)
421{
422 TVector xx( *((TVector*)trs->At(i)) );
423 if( xx.IsValid() && xx(0)>0 )
424 {
425 if(xx(0)>=2)
426 {
427// 1)Calculate Qtrm
428 xx(6)=(this->Qtrm(i));
429
430 }else{
431 xx(6)=xx(1);
432 }
433// 2)Calculate Wpi,Wk,Wp
434 this->GetPcode(xx(6),xx(10)/1000.);
435 xx(7)=GetWpi();
436 xx(8)=GetWk();
437 xx(9)=GetWp();
438// 3)Print table
439 if(xx(0)>0){
440 cout<<xx(0)<<" ";
441 for(Int_t j=1;j<11;j++){
442 cout.width(7);cout.precision(5);cout<<xx(j);
443 }
444 cout<<endl;
445 }
446// 4)Update data in TVector
447 TClonesArray &arr=*trs;
448 new(arr[i])TVector(xx);
449 }
450 else
451 {/*cout<<"No data for track "<<i<<endl;*/}
452}// End loop for tracks
453}
454void AliTPCPid::Reset(void)
455{
456 //
457 // reset
458 //
459 for(Int_t i=0;i<trs->GetEntries();i++){
460 TVector xx(0,11);
461 TClonesArray &arr=*trs;
462 new(arr[i])TVector(xx);
463 }
464}
465//-----------------------------------------------------------
466AliTPCPid::AliTPCPid(Int_t ntrack):TObject(),
467 fCutKa(0),
468 fCutPr(0),
469 fCutKaTune(0.),
470 fCutPrTune(0.),
471 fSigmin(0.),
472 fSilent(0),
473 fmxtrs(0),
474 trs(0),
475 fqtot(0.),
476 fWpi(0.),
477 fWk(0.),
478 fWp(0.),
479 fRpik(0.),
480 fRppi(0.),
481 fRpka(0.),
482 fRp(0.),
483 fPcode(0)
484{
485 trs = new TClonesArray("TVector",ntrack);
486 TClonesArray &arr=*trs;
487 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
488
489
490 //fCutKa = new TF1("fkaons","[0]/x/x+[1]",0.1,1.2);
491 //fCutPr = new TF1("fprotons","[0]/x/x +[1]",0.2,1.2);
492 TF1 *frmska=0;
493
494 frmska = new TF1("x_frmska","1.46-7.82*x+16.78*x^2-15.53*x^3+5.24*x^4 ",
495 0.1,1.2);
496 fCutKa = new TF1("fkaons",
497 "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",
498 0.1,1.2);
499 fCutPr = new TF1("fprotons",
500 "0.83*(15.32-56.094*x+89.962*x^2-66.1856*x^3+18.4052*x^4)",
501 0.2,1.2);
502 fCutKaTune=1.1; // 0.92;
503 fCutPrTune=1.0; //0.80;
504
505const int kinf=10;
506// Ncut Pmom pilo pihi klo khi plo phi
507// fcut[j] [0] [1] [2] [3] [4] [5] [6]
508//----------------------------------------------------------------
509 SetCut( 1, 0.120, 0. , 0. , kinf , kinf , kinf , kinf );
510 SetCut( 2, 0.200, 0. , 6.0 , 6.0 , kinf , kinf , kinf ); //120-200
511 SetCut( 3, 0.400, 0. , 3.5 , 3.5 , 9.0 , 9.0 , kinf ); //200-400
512 SetCut( 4, 0.450, 0. , 1.9 , 1.9 , 4.0 , 4.0 , kinf ); //400-450
513//----------------------------------------------------------------
514 SetCut( 5, 0.500, 0.976, 0.108, 1.484 , 0.159 , 3.5 , kinf ); //450-500
515 SetCut( 6, 0.550, 0.979, 0.108, 1.376 , 0.145 , 3.0 , kinf ); //500-550
516//----------------------------------------------------------------
517 SetCut( 7, 0.600, 0.984, 0.111, 1.295 , 0.146 , 2.7 , kinf ); //550-600
518 SetCut( 8, 0.650, 0.989, 0.113, 1.239 , 0.141 , 2.5 , kinf ); //600-650
519 SetCut( 9, 0.730, 0.995, 0.109, 1.172 , 0.132 , 2.0 , kinf ); //650-730
520//----------------------------------------------------------------
521 SetCut( 10, 0.830, 1.008, 0.116, 1.117 , 0.134 , 1.703, 0.209 ); //730-830
522 SetCut( 11, 0.930, 1.019, 0.115, 1.072 , 0.121 , 1.535, 0.215 ); //830-930
523 SetCut( 12, 1.230, 1.035, 0.117, 1.053 , 0.140 ,1.426, 0.270); //930-1030
524 //------------------------ pi,K ---------------------
525 faprob[0][0]=33219.; faprob[1][0]=1971.; // faprob[0][i] - pions
526 faprob[0][1]=28828.; faprob[1][1]=1973.; // faprob[1][i] - kaons
527 //---------------------------------------------------
528 faprob[0][2]=24532.; faprob[1][2]=1932.; faprob[2][2]=1948.;
529 faprob[0][3]=20797.; faprob[1][3]=1823.; faprob[2][3]=1970.;
530 faprob[0][4]=27017.; faprob[1][4]=2681.; faprob[2][4]=2905.;
531 //------------------------ pi,K,p -------------------
532 faprob[0][5]= 24563.; faprob[1][5]=2816.; faprob[2][5]=3219.;
533 faprob[0][6]= 16877.; faprob[1][6]=2231.; faprob[2][6]=2746.;
534 faprob[0][7]= 11557.; faprob[1][7]=1681; faprob[2][7]=2190.;
535
536 fSilent=kTRUE;
537}
538//-----------------------------------------------------------
539
540
541