TPC/AliTPCPid.cxx

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