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