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