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             if(xx(0)>=2)
 426                 {
 427 //       1)Calculate Qtrm
 428                     xx(6)=(this->Qtrm(i));
 429
 430                  }else{
 431                      xx(6)=xx(1);
 432                  }
 433 //       2)Calculate Wpi,Wk,Wp
 434             this->GetPcode(xx(6),xx(10)/1000.);
 435             xx(7)=GetWpi();
 436             xx(8)=GetWk();
 437             xx(9)=GetWp();
 438 //       3)Print table
 439             if(xx(0)>0){
 440                     cout<<xx(0)<<" ";
 441                     for(Int_t j=1;j<11;j++){
 442                         cout.width(7);cout.precision(5);cout<<xx(j);
 443                     }
 444                     cout<<endl;
 445                 }
 446 //        4)Update data in TVector
 447             TClonesArray &arr=*trs;
 448             new(arr[i])TVector(xx);
 449         }
 450     else
 451       {/*cout<<"No data for track "<<i<<endl;*/}
 452 }// End loop for tracks
 453 }
 454 void AliTPCPid::Reset(void)
 455 {
 456   //
 457   // reset
 458   //
 459   for(Int_t i=0;i<trs->GetEntries();i++){
 460     TVector xx(0,11);
 461     TClonesArray &arr=*trs;
 462     new(arr[i])TVector(xx);
 463   }
 464 }
 465 //-----------------------------------------------------------
 466 AliTPCPid::AliTPCPid(Int_t ntrack):TObject(),
 467                   fCutKa(0),
 468                   fCutPr(0),
 469                   fCutKaTune(0.),
 470                   fCutPrTune(0.),
 471                   fSigmin(0.),
 472                   fSilent(0),
 473                   fmxtrs(0),
 474                   trs(0),
 475                   fqtot(0.),
 476                   fWpi(0.),
 477                   fWk(0.),
 478                   fWp(0.),
 479                   fRpik(0.),
 480                   fRppi(0.),
 481                   fRpka(0.),
 482                   fRp(0.),
 483                   fPcode(0)
 484 {
 485     trs = new TClonesArray("TVector",ntrack);
 486     TClonesArray &arr=*trs;
 487     for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
 488
 489
 490     //fCutKa = new TF1("fkaons","[0]/x/x+[1]",0.1,1.2);
 491     //fCutPr = new TF1("fprotons","[0]/x/x +[1]",0.2,1.2);
 492     TF1 *frmska=0;
 493
 494     frmska = new TF1("x_frmska","1.46-7.82*x+16.78*x^2-15.53*x^3+5.24*x^4 ",
 495                 0.1,1.2);
 496     fCutKa = new TF1("fkaons",
 497            "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",
 498            0.1,1.2);
 499     fCutPr = new TF1("fprotons",
 500                      "0.83*(15.32-56.094*x+89.962*x^2-66.1856*x^3+18.4052*x^4)",
 501                      0.2,1.2);
 502     fCutKaTune=1.1; // 0.92;
 503     fCutPrTune=1.0; //0.80;
 504
 505 const int kinf=10;
 506 //         Ncut Pmom   pilo  pihi    klo    khi     plo    phi
 507 //       fcut[j] [0]    [1]    [2]    [3]    [4]     [5]    [6]
 508 //----------------------------------------------------------------
 509     SetCut(  1, 0.120,  0.  ,  0.  , kinf  , kinf   , kinf  , kinf  );
 510     SetCut(  2, 0.200,  0.  ,  6.0 , 6.0  , kinf   , kinf  , kinf  ); //120-200
 511     SetCut(  3, 0.400,  0.  ,  3.5 , 3.5  , 9.0   , 9.0  , kinf  ); //200-400
 512     SetCut(  4, 0.450,  0.  ,  1.9 , 1.9  , 4.0   , 4.0  , kinf  ); //400-450
 513 //----------------------------------------------------------------
 514     SetCut(  5, 0.500, 0.976, 0.108, 1.484 , 0.159  , 3.5  , kinf  );  //450-500
 515     SetCut(  6, 0.550, 0.979, 0.108, 1.376 , 0.145  , 3.0  , kinf  );  //500-550
 516 //----------------------------------------------------------------
 517     SetCut(  7, 0.600, 0.984, 0.111, 1.295 , 0.146 , 2.7  , kinf  );   //550-600
 518     SetCut(  8, 0.650, 0.989, 0.113, 1.239 , 0.141 , 2.5  , kinf  );   //600-650
 519     SetCut(  9, 0.730, 0.995, 0.109, 1.172 , 0.132 , 2.0  , kinf  );   //650-730
 520 //----------------------------------------------------------------
 521     SetCut( 10, 0.830, 1.008, 0.116, 1.117 , 0.134 , 1.703, 0.209 ); //730-830
 522     SetCut( 11, 0.930, 1.019, 0.115, 1.072 , 0.121 , 1.535, 0.215 ); //830-930
 523     SetCut( 12, 1.230, 1.035, 0.117, 1.053 , 0.140  ,1.426, 0.270); //930-1030
 524     //------------------------ pi,K ---------------------
 525     faprob[0][0]=33219.;   faprob[1][0]=1971.;   // faprob[0][i] -    pions
 526     faprob[0][1]=28828.;   faprob[1][1]=1973.; // faprob[1][i] -    kaons
 527     //---------------------------------------------------
 528     faprob[0][2]=24532.;   faprob[1][2]=1932.; faprob[2][2]=1948.;
 529     faprob[0][3]=20797.;   faprob[1][3]=1823.; faprob[2][3]=1970.;
 530     faprob[0][4]=27017.;   faprob[1][4]=2681.; faprob[2][4]=2905.;
 531     //------------------------ pi,K,p -------------------
 532     faprob[0][5]= 24563.;    faprob[1][5]=2816.;  faprob[2][5]=3219.;
 533     faprob[0][6]= 16877.;    faprob[1][6]=2231.;  faprob[2][6]=2746.;
 534     faprob[0][7]= 11557.;    faprob[1][7]=1681;   faprob[2][7]=2190.;
 535
 536     fSilent=kTRUE;
 537 }
 538 //-----------------------------------------------------------
 539
 540
 541