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