]> 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
Fixes in order to eliminate compilation warnings
[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)
11b3ae13 6//
23efe5f1 7Float_t AliITSPid::qtrm(Int_t track)
8{
9 TVector q(*( this->GetVec(track) ));
38f7c306 10 Int_t ml=(Int_t)q(0);
11 if(ml<1)return 0.;
12 if(ml>6)ml=6;
13 float vf[6];
14 Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
15 if(nl==0)return 0.;
23efe5f1 16 switch(nl){
17 case 1:q(6)=q(1); break;
18 case 2:q(6)=(q(1)+q(2))/2.; break;
19 default:
38f7c306 20 for(int fi=0;fi<2;fi++){
21 Int_t swap;
22 do{ swap=0; float qmin=vf[fi];
23 for(int j=fi+1;j<nl;j++)
24 if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
25 }while(swap==1);
26 }
23efe5f1 27 q(6)= (vf[0]+vf[1])/2.;
28 break;
29 }
30 for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
31 return (q(6));
32}
38f7c306 33
79b921e1 34Float_t AliITSPid::qtrm(Float_t qarr[6],Int_t narr)
35{
38f7c306 36 Float_t q[6],qm,qmin;
79b921e1 37 Int_t nl,ml;
38f7c306 38 if(narr>0&&narr<7){ml=narr;}else{return 0;};
39 nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
79b921e1 40 if(nl==0)return 0.;
79b921e1 41 switch(nl){
42 case 1:qm=q[0]; break;
43 case 2:qm=(q[0]+q[1])/2.; break;
44 default:
38f7c306 45 Int_t swap;
46 for(int fi=0;fi<2;fi++){
47 do{ swap=0; qmin=q[fi];
48 for(int j=fi+1;j<nl;j++)
49 if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
50 }while(swap==1);
51 }
52 qm= (q[0]+q[1])/2.;
79b921e1 53 break;
54 }
55 return qm;
56}
11b3ae13 57
58Int_t AliITSPid::wpik(Float_t pm,Float_t q)
23efe5f1 59{
11b3ae13 60 Double_t par[6];
61 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
62 fggpi->SetParameters(par);
23efe5f1 63
11b3ae13 64 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
65 fggka->SetParameters(par);
38f7c306 66
11b3ae13 67 Float_t ppi=fggpi->Eval(q);
68 Float_t pka=fggka->Eval(q);
69 Float_t p=ppi+pka;
70 /*
71 if(!fSilent){
72 fggka->Print();
73 fggpi->Print();
74 if(p>0)cout<<" ppi,pka="<<ppi/p<<" "<<pka/p<<endl;
75 }
76 */
38f7c306 77
11b3ae13 78 if(p>0){
79 ppi=ppi/p;
80 pka=pka/p;
81 fWp=0.; fWpi=ppi; fWk=pka;
82 if( pka>ppi){return fPcode=321;}else{return fPcode=211;}
83 }else{return 0;}
23efe5f1 84}
85//-----------------------------------------------------------
11b3ae13 86Int_t AliITSPid::wpikp(Float_t pm,Float_t q)
23efe5f1 87{
11b3ae13 88 Double_t par[6];
89 for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
90 fggpi->SetParameters(par);
23efe5f1 91
11b3ae13 92 for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
93 fggka->SetParameters(par);
23efe5f1 94
11b3ae13 95 for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
96 fggpr->SetParameters(par);
38f7c306 97
11b3ae13 98 Float_t p,ppi,pka,ppr;
99 if( q>(fggpr->GetParameter(1)+fggpr->GetParameter(2)) )
100 { p=1.0; ppr=1.0; ppi=pka=0.0;
101 }else{
102 ppi=fggpi->Eval(q);
103 pka=fggka->Eval(q);
104 ppr=fggpr->Eval(q);
105 p=ppi+pka+ppr;
106 }
107 if(p>0){
108 ppi=ppi/p;
109 pka=pka/p;
110 ppr=ppr/p;
111 fWp=ppr; fWpi=ppi; fWk=pka;
112 //if(!fSilent)cout<<" ppi,pka,ppr="<<ppi<<" "<<pka<<" "<<ppr<<endl;
38f7c306 113
11b3ae13 114 if( ppi>pka&&ppi>ppr )
115 {return fPcode=211;}
116 else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
117 }
38f7c306 118
11b3ae13 119 }else{return 0;}
23efe5f1 120}
121//-----------------------------------------------------------
122Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
123{
088e0b8d 124 Info("GetPcode","method not implemented - Inputs TClonesArray *%x , Float_t %f",rps,pm);
23efe5f1 125 return 0;
126}
127//-----------------------------------------------------------
79b921e1 128Int_t AliITSPid::GetPcode(AliTPCtrack *track)
129{
130 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
131 Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
132 if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
133 if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
134 Float_t lam=TMath::ATan(par[3]);
135 Float_t pt_1=TMath::Abs(par[4]);
136 Float_t mom=1./(pt_1*TMath::Cos(lam));
137 Float_t dedx=track->GetdEdx();
138 Int_t pcode=GetPcode(dedx/40.,mom);
def162f4 139// cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 140 return pcode?pcode:211;
141 }
142//------------------------------------------------------------
79b921e1 143Int_t AliITSPid::GetPcode(AliITSIOTrack *track)
144{
145 Double_t px,py,pz;
146 px=track->GetPx();
147 py=track->GetPy();
148 pz=track->GetPz();
149 Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
150//???????????????????
def162f4 151 // Float_t dedx=1.0;
152 Float_t dedx=track->GetdEdx();
79b921e1 153//???????????????????
154 Int_t pcode=GetPcode(dedx,mom);
def162f4 155// cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 156return pcode?pcode:211;
157}
79b921e1 158//-----------------------------------------------------------
159Int_t AliITSPid::GetPcode(AliITStrackV2 *track)
160{
161 if(track==0)return 0;
38f7c306 162 // track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
79b921e1 163 track->PropagateTo(3.,0.0028,65.19);
79b921e1 164 track->PropagateToVertex();
165 Double_t xk,par[5]; track->GetExternalParameters(xk,par);
166 Float_t lam=TMath::ATan(par[3]);
167 Float_t pt_1=TMath::Abs(par[4]);
168 Float_t mom=0.;
169 if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
170 Float_t dedx=track->GetdEdx();
def162f4 171// cout<<"lam,pt_1,mom,dedx="<<lam<<","<<pt_1<<","<<mom<<","<<dedx<<endl;
79b921e1 172 Int_t pcode=GetPcode(dedx,mom);
def162f4 173// cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
79b921e1 174return pcode?pcode:211;
175}
176//-----------------------------------------------------------
23efe5f1 177Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
178{
38f7c306 179 fWpi=fWk=fWp=0.; fPcode=0;
23efe5f1 180
11b3ae13 181 if ( pm<=0.400 )
182 { if( q<fCutKa->Eval(pm) )
183 {return pion();}
184 else{ if( q<fCutPr->Eval(pm) )
185 {return kaon();}
186 else{return proton();}
187 }
188 }
189 if ( pm<=0.750 )
190 if ( q>fCutPr->Eval(pm) )
191 {return proton();} else {return wpik(pm,q);};
192 if( pm<=1.10 ){ return wpikp(pm,q); }
23efe5f1 193 return fPcode;
194}
195//-----------------------------------------------------------
196void AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
197 Float_t klo,Float_t khi,Float_t plo,Float_t phi)
198{
199 cut[n][0]=pm;
200 cut[n][1]=pilo;
201 cut[n][2]=pihi;
202 cut[n][3]=klo;
203 cut[n][4]=khi;
204 cut[n][5]=plo;
205 cut[n][6]=phi;
206 return ;
207}
38f7c306 208//------------------------------------------------------------
23efe5f1 209void AliITSPid::SetVec(Int_t ntrack,TVector info)
210{
211TClonesArray& arr=*trs;
212 new( arr[ntrack] ) TVector(info);
213}
214//-----------------------------------------------------------
215TVector* AliITSPid::GetVec(Int_t ntrack)
216{
217TClonesArray& arr=*trs;
218 return (TVector*)arr[ntrack];
219}
220//-----------------------------------------------------------
221void AliITSPid::SetEdep(Int_t track,Float_t Edep)
222{
223 TVector xx(0,11);
224 if( ((TVector*)trs->At(track))->IsValid() )
225 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
226 Int_t j=(Int_t)xx(0); if(j>4)return;
227 xx(++j)=Edep;xx(0)=j;
228 TClonesArray &arr=*trs;
229 new(arr[track])TVector(xx);
230}
231//-----------------------------------------------------------
232void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
233{
234 TVector xx(0,11);
235 if( ((TVector*)trs->At(track))->IsValid() )
236 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
237 xx(10)=Pmom;
238 TClonesArray &arr=*trs;
239 new(arr[track])TVector(xx);
240}
241//-----------------------------------------------------------
242void AliITSPid::SetPcod(Int_t track,Int_t partcode)
243{
244 TVector xx(0,11);
245 if( ((TVector*)trs->At(track))->IsValid() )
246 {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
247 if(xx(11)==0)
248 {xx(11)=partcode; mxtrs++;
249 TClonesArray &arr=*trs;
250 new(arr[track])TVector(xx);
251 }
252}
253//-----------------------------------------------------------
254void AliITSPid::Print(Int_t track)
255{cout<<mxtrs<<" tracks in AliITSPid obj."<<endl;
256 if( ((TVector*)trs->At(track))->IsValid() )
257 {TVector xx( *((TVector*)trs->At(track)) );
258 xx.Print();
259 }
260 else
261 {cout<<"No data for track "<<track<<endl;return;}
262}
263//-----------------------------------------------------------
264void AliITSPid::Tab(void)
265{
266if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
267cout<<"------------------------------------------------------------------------"<<endl;
268cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
269 " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
270cout<<"------------------------------------------------------------------------"<<endl;
271for(Int_t i=0;i<trs->GetEntries();i++)
272{
273 TVector xx( *((TVector*)trs->At(i)) );
274 if( xx.IsValid() && xx(0)>0 )
275 {
276 TVector xx( *((TVector*)trs->At(i)) );
277 if(xx(0)>=2)
278 {
279// 1)Calculate Qtrm
280 xx(6)=(this->qtrm(i));
281
282 }else{
283 xx(6)=xx(1);
284 }
285// 2)Calculate Wpi,Wk,Wp
286 this->GetPcode(xx(6),xx(10)/1000.);
287 xx(7)=GetWpi();
288 xx(8)=GetWk();
289 xx(9)=GetWp();
290// 3)Print table
291 if(xx(0)>0){
def162f4 292// cout<<xx(0)<<" ";
23efe5f1 293 for(Int_t j=1;j<11;j++){
38f7c306 294 if(i<7){ cout.width(7);cout.precision(4);cout<<xx(j);}
295 if(i>7){ cout.width(7);cout.precision(5);cout<<xx(j);}
23efe5f1 296 }
297 cout<<endl;
298 }
299// 4)Update data in TVector
300 TClonesArray &arr=*trs;
301 new(arr[i])TVector(xx);
302 }
303 else
304 {/*cout<<"No data for track "<<i<<endl;*/}
305}// End loop for tracks
306}
307void AliITSPid::Reset(void)
308{
309 for(Int_t i=0;i<trs->GetEntries();i++){
310 TVector xx(0,11);
311 TClonesArray &arr=*trs;
312 new(arr[i])TVector(xx);
313 }
314}
315//-----------------------------------------------------------
316AliITSPid::AliITSPid(Int_t ntrack)
317{
38f7c306 318 fSigmin=0.01;
23efe5f1 319 trs = new TClonesArray("TVector",ntrack);
320 TClonesArray &arr=*trs;
321 for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
322 mxtrs=0;
11b3ae13 323 //
324 fCutKa=new TF1("fcutka","pol4",0.05,0.4);
325 Double_t ka[5]={25.616, -161.59, 408.97, -462.17, 192.86};
326 fCutKa->SetParameters(ka);
327 //
328 fCutPr=new TF1("fcutpr","[0]/x/x+[1]",0.05,1.1);
329 Double_t pr[2]={0.70675,0.4455};
330 fCutPr->SetParameters(pr);
331 //
332 //---------- signal fit ----------
333{//Pions
334fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
335 Double_t parpi_0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
336 1.1405329025e+05, -3.8289694496e+04 };
337 fGGpi[0]->SetParameters(parpi_0);
338fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
339 Double_t parpi_1[10]={ 1.0791668283e-02, 9.7347716496e-01 };
340 fGGpi[1]->SetParameters(parpi_1);
341fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
342 Double_t parpi_2[10]={ 5.8191602279e-04, 9.7285601334e-02 };
343 fGGpi[2]->SetParameters(parpi_2);
344fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
345 Double_t parpi_3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
346 6.2900977606e+03, -2.2935862292e+03 };
347 fGGpi[3]->SetParameters(parpi_3);
348fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
349 Double_t parpi_4[10]={ 9.0419011783e-03, 1.1628922525e+00 };
350 fGGpi[4]->SetParameters(parpi_4);
351fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
352 Double_t parpi_5[10]={ 1.8324872519e-03, 2.1503968838e-01 };
353 fGGpi[5]->SetParameters(parpi_5);
354}//End Pions
355{//Kaons
356fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
357 Double_t parka_0[20]={
358 -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
359 -3.7123338009e+03,1.6003647641e+03 };
360 fGGka[0]->SetParameters(parka_0);
361fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
362 Double_t parka_1[20]={
363 2.5181172905e-01,8.7566001814e-01 };
364 fGGka[1]->SetParameters(parka_1);
365fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
366 Double_t parka_2[20]={
367 8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
368 -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
369 5.7545491696e+01 };
370 fGGka[2]->SetParameters(parka_2);
371}//End Kaons
372{//Protons
373fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
374 Double_t parpr_0[10]={
375 6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
376 -3.5269200901e+03,1.2859128345e+03 };
377 fGGpr[0]->SetParameters(parpr_0);
378fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
379 Double_t parpr_1[10]={
380 9.4970837607e-01,7.3573504201e-01 };
381 fGGpr[1]->SetParameters(parpr_1);
382fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
383 Double_t parpr_2[10]={
384 1.2498403757e-01,2.7845072306e-02 };
385 fGGpr[2]->SetParameters(parpr_2);
386}//End Protons
387 //----------- end fit -----------
388
389 fggpr=new TF1("ggpr","gaus",0.4,1.2);
390 fggpi=new TF1("ggpi","gaus+gaus(3)",0.4,1.2);
391 fggka=new TF1("ggka","gaus",0.4,1.2);
392
393 //-------------------------------------------------
23efe5f1 394const int inf=10;
395// Ncut Pmom pilo pihi klo khi plo phi
396// cut[j] [0] [1] [2] [3] [4] [5] [6]
397//----------------------------------------------------------------
398 SetCut( 1, 0.12 , 0. , 0. , inf , inf , inf , inf );
399 SetCut( 2, 0.20 , 0. , 6.0 , 6.0 , inf , inf , inf );
400 SetCut( 3, 0.30 , 0. , 3.5 , 3.5 , 9.0 , 9.0 , inf );
401 SetCut( 4, 0.41 , 0. , 1.9 , 1.9 , 4.0 , 4.0 , inf );
402//----------------------------------------------------------------
38f7c306 403 SetCut( 5, 0.47 , 0.935, 0.139, 1.738 , 0.498 , 3.5 , inf ); //410-470
404 SetCut( 6, 0.53 , 0.914, 0.136, 1.493 , 0.436 , 3.0 , inf ); //470-530
23efe5f1 405//----------------------------------------------------------------
38f7c306 406 SetCut( 7, 0.59 , 0.895, 0.131, 1.384 , 0.290 , 2.7 , inf ); //530-590
407 SetCut( 8, 0.65 , 0.887, 0.121, 1.167 , 0.287 , 2.5 , inf ); //590-650
408 SetCut( 9, 0.73 , 0.879, 0.120, 1.153 , 0.257 , 2.0 , inf ); //650-730
23efe5f1 409//----------------------------------------------------------------
38f7c306 410 SetCut( 10, 0.83 , 0.880, 0.126, 1.164 , 0.204 , 2.308 , 0.297 ); //730-830
411 SetCut( 11, 0.93 , 0.918, 0.145, 1.164 , 0.204 , 2.00 , 0.168 ); //830-930
412 SetCut( 12, 1.03 , 0.899, 0.128, 1.164 , 0.204 ,1.80 , 0.168);
413 //------------------------ pi,K ---------------------
414 aprob[0][0]=1212; aprob[1][0]=33.; // aprob[0][i] - const for pions,cut[i+5]
415 aprob[0][1]=1022; aprob[1][1]=46.2 ; // aprob[1][i] - kaons
416 //---------------------------------------------------
417 aprob[0][2]= 889.7; aprob[1][2]=66.58; aprob[2][2]=14.53;
418 aprob[0][3]= 686.; aprob[1][3]=88.8; aprob[2][3]=19.27;
419 aprob[0][4]= 697.; aprob[1][4]=125.6; aprob[2][4]=28.67;
420 //------------------------ pi,K,p -------------------
421 aprob[0][5]= 633.7; aprob[1][5]=100.1; aprob[2][5]=37.99; // aprob[2][i] - protons
422 aprob[0][6]= 469.5; aprob[1][6]=20.74; aprob[2][6]=25.43;
423 aprob[0][7]= 355.; aprob[1][7]=
424 355.*(20.74/469.5); aprob[2][7]=34.08;
23efe5f1 425}
11b3ae13 426//End AliITSPid.cxx