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