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