Effective C++ warnings
[u/mrichter/AliRoot.git] / ITS / AliITSPid.cxx
1 //------------------------------------------------------------//
2 // Class for identification of pions,kaons and protons in ITS //
3 // Prior particles population (probabilities) are taken from  //
4 // Hijing event generator.                                   //
5 //------------------------------------------------------------//
6 // #include <stdlib.h>
7 #include "AliITSPid.h"
8 //#include "TMath.h"
9 #include <Riostream.h>
10 #include <TClonesArray.h>
11 #include <TVector.h>
12 #include "AliKalmanTrack.h"
13 #include "AliITSIOTrack.h"
14 #include "AliITStrackV2.h"
15 #include <TF1.h>
16
17 ClassImp(AliITSPid)
18
19
20 AliITSPid::AliITSPid(const AliITSPid &source) : TObject(source){
21     // Copy constructor. This is a function which is not allowed to be
22     // done: protected It exits with an error.
23     // Inputs:
24     //      AliITSPid &source  An AliITSPid class.
25     // Outputs:
26     //      none.
27     // Return:
28     //      none.
29
30   Error("AliITSPid","You are not allowed to make a copy of the AliITSPid\n");
31 }
32 //______________________________________________________________________
33 AliITSPid& AliITSPid::operator=(const AliITSPid& /* source */){
34     // Assignment operator. This is a function which is not allowed to be
35     // done to the ITS. It exits with an error.
36     // Inputs:
37     //      AliITSPid &source  An AliITSPid class.
38     // Outputs:
39     //      none.
40     // Return:
41     //      none.
42
43   Error("=operator","Assignment operator not allowed for this class\n");
44   return *this;
45 }
46
47 //
48   Float_t AliITSPid::Qtrm(Int_t track) {
49 //
50 // This calculates truncated mean signal for given track.
51 //
52     TVector q(*( this->GetVec(track)  ));
53     Int_t ml=(Int_t)q(0);
54     if(ml<1)return 0.;
55     if(ml>6)ml=6;
56     float vf[6];
57     Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
58     if(nl==0)return 0.;
59     switch(nl){
60       case  1:q(6)=q(1); break;
61       case  2:q(6)=(q(1)+q(2))/2.; break;
62       default:
63         for(int fi=0;fi<2;fi++){
64           Int_t swap;
65           do{ swap=0; float qmin=vf[fi];
66           for(int j=fi+1;j<nl;j++)
67             if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
68           }while(swap==1);
69         }
70         q(6)= (vf[0]+vf[1])/2.;
71         break;
72     }
73     for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
74     return (q(6));
75 }
76
77 Float_t AliITSPid::Qtrm(Float_t qarr[6],Int_t narr) const{
78 //
79 //This calculates truncated mean signal for given signal set.
80 //
81   Float_t q[6],qm,qmin;
82   Int_t nl,ml;
83   if(narr>0&&narr<7){ml=narr;}else{return 0;};
84   nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
85   if(nl==0)return 0.;
86     switch(nl){
87       case  1:qm=q[0]; break;
88       case  2:qm=(q[0]+q[1])/2.; break;
89       default:
90         Int_t swap;
91         for(int fi=0;fi<2;fi++){
92           do{ swap=0; qmin=q[fi];
93           for(int j=fi+1;j<nl;j++)
94             if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
95           }while(swap==1);
96         }
97         qm= (q[0]+q[1])/2.;
98         break;
99     }
100     return qm;
101 }
102
103 Int_t   AliITSPid::Wpik(Float_t pm,Float_t q){
104   //Calcutates probabilityes of pions and kaons
105   //Returns particle code for dominant probability.
106   Double_t par[6];
107   for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
108   fggpi->SetParameters(par);
109
110   for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
111   fggka->SetParameters(par);
112
113   Float_t ppi=fggpi->Eval(q);
114   Float_t pka=fggka->Eval(q);
115   Float_t p=ppi+pka;
116   /*
117   if(!fSilent){
118     fggka->Print();
119     fggpi->Print();
120     if(p>0)cout<<" ppi,pka="<<ppi/p<<"  "<<pka/p<<endl;
121   }
122   */
123
124   if(p>0){
125     ppi=ppi/p; 
126     pka=pka/p;
127     fWp=0.; fWpi=ppi; fWk=pka;
128     if( pka>ppi){return fPcode=321;}else{return fPcode=211;}
129   }else{return 0;}
130 }
131 //-----------------------------------------------------------
132 Int_t   AliITSPid::Wpikp(Float_t pm,Float_t q){
133   //
134   //Calcutates probabilityes of pions,kaons and protons.
135   //Returns particle code for dominant probability.
136   Double_t par[6];
137   for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
138   fggpi->SetParameters(par);
139
140   for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
141   fggka->SetParameters(par);
142
143   for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
144   fggpr->SetParameters(par);
145
146   Float_t p,ppi,pka,ppr;
147   if( q>(fggpr->GetParameter(1)+fggpr->GetParameter(2)) )
148       { p=1.0; ppr=1.0; ppi=pka=0.0;
149     }else{ 
150     ppi=fggpi->Eval(q);
151     pka=fggka->Eval(q);
152     ppr=fggpr->Eval(q);
153     p=ppi+pka+ppr;
154   }
155   if(p>0){
156     ppi=ppi/p; 
157     pka=pka/p;
158     ppr=ppr/p;
159     fWp=ppr; fWpi=ppi; fWk=pka;
160     //if(!fSilent)cout<<" ppi,pka,ppr="<<ppi<<"  "<<pka<<" "<<ppr<<endl;
161
162    if( ppi>pka&&ppi>ppr )
163            {return fPcode=211;}
164    else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
165    }
166
167   }else{return 0;}
168 }
169 //-----------------------------------------------------------
170 Int_t   AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
171 {
172   //Returns particle code
173   Info("GetPcode","method not implemented - Inputs TClonesArray *%x , Float_t %f",rps,pm); 
174     return 0;    
175 }
176 //-----------------------------------------------------------
177 Int_t   AliITSPid::GetPcode(AliKalmanTrack *track)
178 {
179   //Returns particle code for given track.
180       Double_t xk,par[5]; track->GetExternalParameters(xk,par);
181       Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
182       if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
183       if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
184       Float_t lam=TMath::ATan(par[3]); 
185       Float_t pt1=TMath::Abs(par[4]);
186       Float_t mom=1./(pt1*TMath::Cos(lam));
187       Float_t dedx=track->GetPIDsignal();
188     Int_t pcode=GetPcode(dedx/40.,mom);
189 //    cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
190     return pcode?pcode:211;
191     }
192 //------------------------------------------------------------
193 Int_t   AliITSPid::GetPcode(AliITSIOTrack *track)
194 {
195   //Returns particle code for given track(V1).
196     Double_t px,py,pz;
197     px=track->GetPx();
198     py=track->GetPy();
199     pz=track->GetPz();
200     Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
201 //???????????????????
202     // Float_t dedx=1.0;
203     Float_t dedx=track->GetdEdx();
204 //???????????????????    
205     Int_t pcode=GetPcode(dedx,mom);
206 //    cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
207 return pcode?pcode:211;
208 }
209 //-----------------------------------------------------------
210 Int_t   AliITSPid::GetPcode(AliITStrackV2 *track)
211 {
212 //Returns particle code for given track(V2).
213   if(track==0)return 0;
214   //      track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
215       track->PropagateTo(3.,0.0028,65.19);
216       track->PropagateToVertex();
217     Double_t xk,par[5]; track->GetExternalParameters(xk,par);
218     Float_t lam=TMath::ATan(par[3]);
219     Float_t pt1=TMath::Abs(par[4]);
220     Float_t mom=0.;
221     if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
222     Float_t dedx=track->GetdEdx();
223 //    cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
224     Int_t pcode=GetPcode(dedx,mom);
225 //    cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
226 return pcode?pcode:211;
227 }
228 //-----------------------------------------------------------
229 Int_t   AliITSPid::GetPcode(Float_t q,Float_t pm)
230 {
231   //Returns particle code for given signal and momentum.
232     fWpi=fWk=fWp=0.;     fPcode=0;
233
234     if ( pm<=0.400 )
235         { if( q<fCutKa->Eval(pm) )
236             {return Pion();}
237         else{ if( q<fCutPr->Eval(pm) )
238                 {return Kaon();}
239             else{return Proton();}
240             } 
241         }
242     if ( pm<=0.750 )
243         if ( q>fCutPr->Eval(pm)  )
244             {return Proton();} else {return Wpik(pm,q);};
245     if( pm<=1.10 ){ return Wpikp(pm,q); }
246     return fPcode;    
247 }
248 //-----------------------------------------------------------
249 void    AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
250                         Float_t klo,Float_t khi,Float_t plo,Float_t phi)
251 {
252   // The cut-table initializer method.
253     fCut[n][0]=pm;
254     fCut[n][1]=pilo;
255     fCut[n][2]=pihi;
256     fCut[n][3]=klo;
257     fCut[n][4]=khi;
258     fCut[n][5]=plo;
259     fCut[n][6]=phi;
260     return ;    
261 }
262 //------------------------------------------------------------
263 void AliITSPid::SetVec(Int_t ntrack,const TVector& info) const
264 {
265   //Store track info in tracls table
266 TClonesArray& arr=*fTrs;
267     new( arr[ntrack] ) TVector(info);
268 }
269 //-----------------------------------------------------------
270 TVector* AliITSPid::GetVec(Int_t ntrack) const
271 {
272   //Get given track from track table 
273 TClonesArray& arr=*fTrs;
274     return (TVector*)arr[ntrack];
275 }
276 //-----------------------------------------------------------
277 void AliITSPid::SetEdep(Int_t track,Float_t Edep)
278 {
279   //Set dEdx for given track
280     TVector xx(0,11);
281     if( ((TVector*)fTrs->At(track))->IsValid() )
282         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
283     Int_t j=(Int_t)xx(0); if(j>4)return;
284     xx(++j)=Edep;xx(0)=j;
285     TClonesArray &arr=*fTrs;
286     new(arr[track])TVector(xx);
287 }
288 //-----------------------------------------------------------
289 void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
290 {
291   //Set momentum for given track
292     TVector xx(0,11);
293     if( ((TVector*)fTrs->At(track))->IsValid() )
294         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
295     xx(10)=Pmom;
296     TClonesArray &arr=*fTrs;
297     new(arr[track])TVector(xx);
298 }
299 //-----------------------------------------------------------
300 void AliITSPid::SetPcod(Int_t track,Int_t partcode)
301 {
302   //Set particle code for given track
303     TVector xx(0,11);
304     if( ((TVector*)fTrs->At(track))->IsValid() )
305         {TVector yy( *((TVector*)fTrs->At(track)) );xx=yy; }
306     if(xx(11)==0)
307         {xx(11)=partcode; fMxtrs++;
308         TClonesArray &arr=*fTrs;
309         new(arr[track])TVector(xx);
310         }
311 }
312 //-----------------------------------------------------------
313 void AliITSPid::Print(Int_t track)
314 {
315   //Prints information for given track
316 cout<<fMxtrs<<" tracks in AliITSPid obj."<<endl;
317     if( ((TVector*)fTrs->At(track))->IsValid() )
318         {TVector xx( *((TVector*)fTrs->At(track)) );
319          xx.Print();
320          }
321     else 
322         {cout<<"No data for track "<<track<<endl;return;}
323 }
324 //-----------------------------------------------------------
325 void AliITSPid::Tab(void)
326 {
327   //Make PID for tracks stored in tracks table
328 if(fTrs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
329 cout<<"------------------------------------------------------------------------"<<endl;
330 cout<<"Nq"<<"   q1  "<<"   q2  "<<"   q3  "<<"   q4  "<<"   q5   "<<
331       " Qtrm    "    <<"  Wpi  "<<"  Wk   "<<"  Wp  "<<"Pmom  "<<endl;
332 cout<<"------------------------------------------------------------------------"<<endl;
333 for(Int_t i=0;i<fTrs->GetEntries();i++)
334 {
335   TVector xx( *((TVector*)fTrs->At(i)) );     
336     if( xx.IsValid() && xx(0)>0 )
337         {
338             TVector xx( *((TVector*)fTrs->At(i)) );
339             if(xx(0)>=2)
340                 {
341 //       1)Calculate Qtrm       
342                     xx(6)=(this->Qtrm(i));
343
344                  }else{
345                      xx(6)=xx(1);
346                  }
347 //       2)Calculate Wpi,Wk,Wp
348             this->GetPcode(xx(6),xx(10)/1000.);
349             xx(7)=GetWpi();
350             xx(8)=GetWk();
351             xx(9)=GetWp();
352 //       3)Print table
353             if(xx(0)>0){
354 //                  cout<<xx(0)<<" ";
355                     for(Int_t j=1;j<11;j++){
356                       if(i<7){ cout.width(7);cout.precision(4);cout<<xx(j);}
357                       if(i>7){ cout.width(7);cout.precision(5);cout<<xx(j);}
358                     }
359                     cout<<endl;
360                 }
361 //        4)Update data in TVector
362             TClonesArray &arr=*fTrs;
363             new(arr[i])TVector(xx);      
364         }
365     else 
366       {/*cout<<"No data for track "<<i<<endl;*/}
367 }// End loop for tracks
368 }
369 void AliITSPid::Reset(void)
370 {
371   //Reset tracks table
372   for(Int_t i=0;i<fTrs->GetEntries();i++){
373     TVector xx(0,11);
374     TClonesArray &arr=*fTrs;
375     new(arr[i])TVector(xx);
376   }
377 }
378 //-----------------------------------------------------------
379 AliITSPid::AliITSPid(Int_t ntrack)
380 {
381   //Constructor for AliITSPid class
382   fSigmin=0.01;
383     fTrs = new TClonesArray("TVector",ntrack);
384     TClonesArray &arr=*fTrs;
385     for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
386     fMxtrs=0;
387     //   
388     fCutKa=new TF1("fcutka","pol4",0.05,0.4);
389     Double_t ka[5]={25.616, -161.59, 408.97, -462.17, 192.86};
390     fCutKa->SetParameters(ka);
391     //
392     fCutPr=new TF1("fcutpr","[0]/x/x+[1]",0.05,1.1);
393     Double_t pr[2]={0.70675,0.4455};
394     fCutPr->SetParameters(pr);
395     //
396     //---------- signal fit ----------
397 {//Pions
398 fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
399   Double_t parpi0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
400    1.1405329025e+05, -3.8289694496e+04  };
401   fGGpi[0]->SetParameters(parpi0);
402 fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
403   Double_t parpi1[10]={ 1.0791668283e-02, 9.7347716496e-01  };
404   fGGpi[1]->SetParameters(parpi1);
405 fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
406   Double_t parpi2[10]={ 5.8191602279e-04, 9.7285601334e-02  };
407   fGGpi[2]->SetParameters(parpi2);
408 fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
409   Double_t parpi3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
410    6.2900977606e+03, -2.2935862292e+03  };
411   fGGpi[3]->SetParameters(parpi3);
412 fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
413   Double_t parpi4[10]={ 9.0419011783e-03, 1.1628922525e+00  };
414   fGGpi[4]->SetParameters(parpi4);
415 fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
416   Double_t parpi5[10]={ 1.8324872519e-03, 2.1503968838e-01  };
417   fGGpi[5]->SetParameters(parpi5);
418 }//End Pions
419 {//Kaons
420 fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
421   Double_t parka0[20]={
422   -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
423   -3.7123338009e+03,1.6003647641e+03  };
424   fGGka[0]->SetParameters(parka0);
425 fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
426   Double_t parka1[20]={
427   2.5181172905e-01,8.7566001814e-01  };
428   fGGka[1]->SetParameters(parka1);
429 fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
430   Double_t parka2[20]={
431   8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
432   -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
433   5.7545491696e+01  };
434   fGGka[2]->SetParameters(parka2);
435 }//End Kaons
436 {//Protons
437 fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
438   Double_t parpr0[10]={
439   6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
440   -3.5269200901e+03,1.2859128345e+03  };
441   fGGpr[0]->SetParameters(parpr0);
442 fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
443   Double_t parpr1[10]={
444   9.4970837607e-01,7.3573504201e-01  };
445   fGGpr[1]->SetParameters(parpr1);
446 fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
447   Double_t parpr2[10]={
448   1.2498403757e-01,2.7845072306e-02  };
449   fGGpr[2]->SetParameters(parpr2);
450 }//End Protons
451     //----------- end fit -----------
452
453     fggpr=new TF1("ggpr","gaus",0.4,1.2);
454     fggpi=new TF1("ggpi","gaus+gaus(3)",0.4,1.2);
455     fggka=new TF1("ggka","gaus",0.4,1.2);
456
457     //-------------------------------------------------
458 const int kInf=10;
459 //         Ncut Pmom   pilo  pihi    klo    khi     plo    phi
460 //       cut[j] [0]    [1]    [2]    [3]    [4]     [5]    [6]
461 //----------------------------------------------------------------
462     SetCut(  1, 0.12 ,  0.  ,  0.  , kInf  , kInf   , kInf  , kInf  );
463     SetCut(  2, 0.20 ,  0.  ,  6.0 , 6.0  , kInf   , kInf  , kInf  );
464     SetCut(  3, 0.30 ,  0.  ,  3.5 , 3.5  , 9.0   , 9.0  , kInf  );
465     SetCut(  4, 0.41 ,  0.  ,  1.9 , 1.9  , 4.0   , 4.0  , kInf  );
466 //----------------------------------------------------------------
467     SetCut(  5, 0.47 , 0.935, 0.139, 1.738 , 0.498  , 3.5  , kInf  );  //410-470
468     SetCut(  6, 0.53 , 0.914, 0.136, 1.493 , 0.436  , 3.0  , kInf  );  //470-530
469 //----------------------------------------------------------------    
470     SetCut(  7, 0.59 , 0.895, 0.131, 1.384 , 0.290 , 2.7  , kInf  );    //530-590
471     SetCut(  8, 0.65 , 0.887, 0.121, 1.167 , 0.287 , 2.5  , kInf  );     //590-650
472     SetCut(  9, 0.73 , 0.879, 0.120, 1.153 , 0.257 , 2.0  , kInf  );     //650-730
473 //----------------------------------------------------------------    
474     SetCut( 10, 0.83 , 0.880, 0.126, 1.164 , 0.204 , 2.308 , 0.297 );       //730-830
475     SetCut( 11, 0.93 , 0.918, 0.145, 1.164 , 0.204 , 2.00 , 0.168 );        //830-930
476     SetCut( 12, 1.03 , 0.899, 0.128, 1.164 , 0.204  ,1.80 , 0.168);
477     //------------------------ pi,K ---------------------
478     fAprob[0][0]=1212;     fAprob[1][0]=33.;   // fAprob[0][i] - const for pions,cut[i+5] 
479     fAprob[0][1]=1022;     fAprob[1][1]=46.2 ; // fAprob[1][i] -           kaons
480     //---------------------------------------------------
481     fAprob[0][2]= 889.7;   fAprob[1][2]=66.58; fAprob[2][2]=14.53;
482     fAprob[0][3]= 686.;    fAprob[1][3]=88.8;  fAprob[2][3]=19.27;   
483     fAprob[0][4]= 697.;    fAprob[1][4]=125.6; fAprob[2][4]=28.67;
484     //------------------------ pi,K,p -------------------
485     fAprob[0][5]= 633.7;   fAprob[1][5]=100.1;   fAprob[2][5]=37.99;   // fAprob[2][i] -  protons
486     fAprob[0][6]= 469.5;   fAprob[1][6]=20.74;   fAprob[2][6]=25.43;
487     fAprob[0][7]= 355.;    fAprob[1][7]=
488                           355.*(20.74/469.5);  fAprob[2][7]=34.08;
489 }
490 //End AliITSPid.cxx