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