[u/mrichter/AliRoot.git] / TPC / Attic / AliTPCPid.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/// \class AliTPCPid

#include "AliTPCPid.h"
#include "TMath.h"
//#include <TVector.h>

#include <TF1.h>
#include <TClonesArray.h>
//#include "AliITSIOTrack.h"
#include "AliKalmanTrack.h"
#include "Riostream.h"
#include <cassert>

using std::cout;
using std::endl;
/// \cond CLASSIMP
ClassImp(AliTPCPid)
/// \endcond
// Correction 13.01.2003 Z.S.,Dubna
//            22.01.2003
//------------------------------------------------------------
// pid class by B. Batyunya
// stupid corrections by M.K.
//
//_______________________________________________________
//________________________________________________________
AliTPCPid::AliTPCPid( const AliTPCPid& r):TObject(r),
		  fCutKa(0),
		  fCutPr(0),
		  fCutKaTune(0.),
		  fCutPrTune(0.),
		  fSigmin(0.),
		  fSilent(0),
		  fmxtrs(0),
		  trs(0),
		  fqtot(0.),
		  fWpi(0.),
		  fWk(0.),
		  fWp(0.),
		  fRpik(0.),
		  fRppi(0.),
		  fRpka(0.),
		  fRp(0.),
		  fPcode(0)
{
  // dummy copy constructor
}
Float_t AliTPCPid::Qcorr(Float_t xc)
{
  /// charge correction

    assert(0);
  Float_t fcorr;
  fcorr=( 0.766 +0.9692*xc -1.267*xc*xc )*( 1.-TMath::Exp(-xc*64.75) );
  if(fcorr<=0.1)fcorr=0.1;
return fqtot/fcorr;
}
//__________________________________________________________
AliTPCPid & AliTPCPid::operator =(const AliTPCPid & param)
{
  /// assignment operator - dummy

  if(this!=&param){
     fSigmin=param.fSigmin;
  }
  return (*this);
}
//-----------------------------------------------------------
Float_t AliTPCPid::Qtrm(Int_t track) const
{
  /// dummy comment (Boris!!!)

    TVector q(*( this->GetVec(track)  ));
    Int_t ml=(Int_t)q(0);
    if(ml<1)return 0.;
    if(ml>6)ml=6;
    float vf[6];
    Int_t nl=0; for(Int_t i=0;i<ml;i++){if(q(i)>fSigmin){vf[nl]=q(i+1);nl++;}}
    if(nl==0)return 0.;
    switch(nl){
      case  1:q(6)=q(1); break;
      case  2:q(6)=(q(1)+q(2))/2.; break;
      default:
	for(int fi=0;fi<2;fi++){
	  Int_t swap;
	  do{ swap=0; float qmin=vf[fi];
	  for(int j=fi+1;j<nl;j++)
	    if(qmin>vf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;};
	  }while(swap==1);
	}
        q(6)= (vf[0]+vf[1])/2.;
        break;
    }
    for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
    return (q(6));
}// --- End AliTPCPid::Qtrm ---

Float_t AliTPCPid::Qtrm(Float_t qarr[6],Int_t narr)
{
  /// ..................

  Float_t q[6],qm,qmin;
  Int_t nl,ml;
  if(narr>0&&narr<7){ml=narr;}else{return 0;};
  nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>fSigmin){q[nl]=qarr[i];nl++;}}
  if(nl==0)return 0.;
    switch(nl){
      case  1:qm=q[0]; break;
      case  2:qm=(q[0]+q[1])/2.; break;
      default:
	Int_t swap;
	for(int fi=0;fi<2;fi++){
	  do{ swap=0; qmin=q[fi];
	  for(int j=fi+1;j<nl;j++)
	    if(qmin>q[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;};
	  }while(swap==1);
	}
        qm= (q[0]+q[1])/2.;
        break;
    }
    return qm;
}// --- End AliTPCPid::Qtrm  ---

Int_t	AliTPCPid::Wpik(Int_t nc,Float_t q)
{
  ///  pi-k

    Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
    Float_t appi,apk;
    qmpi =fcut[nc][1];
    sigpi=fcut[nc][2];
    qmk  =fcut[nc][3];
    sigk =fcut[nc][4];
    appi = faprob[0][nc-5];
    apk  = faprob[1][nc-5];
    if( !fSilent ){
    cout<<"qmpi,sigpi,qmk,sigk="<<qmpi<<"  "<<sigpi<<"  "<<qmk<<"  "<<sigk<<endl;
    cout<<"appi,apk="<<appi<<","<<apk<<endl;
    }
    Float_t dqpi=(q-qmpi)/sigpi;
    Float_t dqk =(q-qmk )/sigk;
    dpi =TMath::Abs(dqpi);
    dk  =TMath::Abs(dqk);
    Double_t dn=appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk);
    if(dn>0.){
      ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
      pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
    }else{fWpi=1;return Pion();}
    Float_t rpik=ppi/(pk+0.0000001); 
    if( !fSilent )
	cout<<"q,dqpi,dqk, Wpik: ppi,pk,rpik="
	<<q<<"  "<<dqpi<<"  "<<dqk<<"  "<<ppi<<"  "<<pk<<"  "<<rpik<<endl;

    fWp=0.; fWpi=ppi; fWk=pk;
    if( pk>ppi){return Kaon();}else{return Pion();}
}
//-----------------------------------------------------------
Int_t	AliTPCPid::Wpikp(Int_t nc,Float_t q)
{
  /// pi-k-p

   Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
   Float_t appi,apk,app;

    qmpi =fcut[nc][1];
    sigpi=fcut[nc][2];
    qmk  =fcut[nc][3];
    sigk =fcut[nc][4];
    qmp  =fcut[nc][5];
    sigp =fcut[nc][6];

    //appi = apk = app = 1.;
    appi = faprob[0][nc-5];
    apk  = faprob[1][nc-5];
    app  = faprob[2][nc-5];

    Double_t dn= appi*TMath::Gaus(q,qmpi,sigpi)
      +apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp);
    if(dn>0.){
    ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn;
    pk = apk*TMath::Gaus(q,qmk, sigk )/dn;
    pp = app*TMath::Gaus(q,qmp, sigp )/dn;
     }
    else{fWpi=1;return Pion();}
    fWp=pp; fWpi=ppi; fWk=pk;
    if( !fSilent ){
cout<<" Wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<";   "<<qmk<<" "<<sigk<<";   "
    <<qmp<<" "<<sigp<<"; "<<endl;
cout<<" faprob: "<<appi<<"  "<<apk<<"  "<<app<<endl;
cout<<" ppi,pk,pp="<<ppi<<"  "<<pk<<"  "<<pp<<endl;
     }
    if( ppi>pk&&ppi>pp )  { return Pion(); }
    if(pk>pp){return Kaon();}else{return Proton();}
}
//-----------------------------------------------------------
Int_t	AliTPCPid::GetPcode(TClonesArray* /*rps*/,Float_t /*pm*/) const
{
  /// dummy ???

    return 0;    
}
//-----------------------------------------------------------
Int_t   AliTPCPid::GetPcode(AliTPCtrack *track)
{
  /// get particle code

      Double_t xk,par[5]; track->GetExternalParameters(xk,par);
      Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
      if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
      if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
      Float_t lam=TMath::ATan(par[3]); 
      Float_t pt1=TMath::Abs(par[4]);
      Float_t mom=1./(pt1*TMath::Cos(lam));
      Float_t dedx=track->GetdEdx();
    Int_t pcode=GetPcode(dedx/40.,mom);
    cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
    return pcode?pcode:211;
    }
/*
//-----------------------------------------------------------
Int_t   AliTPCPid::GetPcode(AliKalmanTrack *track)
{
  //
  // get particle code
  //
  if(track==0)return 0;
  //      track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
      track->PropagateTo(3.,0.0028,65.19);
      track->PropagateToVertex();
    Double_t xk,par[5]; track->GetExternalParameters(xk,par);
    Float_t lam=TMath::ATan(par[3]);
    Float_t pt1=TMath::Abs(par[4]);
    Float_t mom=0.;
    if( (pt1*TMath::Cos(lam))!=0. ){ mom=1./(pt1*TMath::Cos(lam)); }else{mom=0.;};
    Float_t dedx=track->GetPIDsignal();
    cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
    Int_t pcode=GetPcode(dedx,mom);
    cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
return pcode?pcode:211;
}
*/
//-----------------------------------------------------------
Int_t	AliTPCPid::GetPcode(Float_t q,Float_t pm)
{
  /// get particle code

    fWpi=fWk=fWp=0.;     fPcode=0;
//1)---------------------- 0-120 MeV/c --------------
    if ( pm<=fcut[1][0] )
	{ fWpi=1.; return Pion(); }
//2)----------------------120-200 Mev/c ( 29.7.2002 ) ------------- 
    if ( pm<=fcut[2][0] )
	{ if( q<fCutKa->Eval(pm) ){fWpi=1.; return Pion(); } 
	                     else {fWk =1.; return  Kaon();} }
//3)----------------------200-400 Mev/c -------------
    if ( pm<=fcut[3][0] ){
	if( q<fCutKa->Eval(pm) )
		{ fWpi=1.;return Pion(); }
	     else
		{ if ( q<=fCutPr->Eval(pm) ) 
			     {fWk=1.;return Kaon();} 
			else {fWp=1.;return Proton();}
		}}
//4)----------------------400-450 Mev/c -------------
    if ( pm<=fcut[4][0] ){
	if( q<fCutKaTune*fCutKa->Eval(pm) )
		{ fWpi=1.;return Pion(); }
	    else
		{ if( q<fCutPr->Eval(pm) ) 
                  {fWk=1.;return Kaon();} else {fWp=1.;return Proton(); }
                }}
//5)----------------------450-500 Mev/c -------------
    if ( pm<=fcut[5][0] ){
	if ( q>fCutPr->Eval(pm) )
	  {fWp=1.;return Proton();} else {return Wpik(5,q);}}
//6)----------------------500-550 Mev/c -------------
    if ( pm<=fcut[6][0] ){
	if ( q>fCutPr->Eval(pm) )
	  {fWp=1.;return Proton();} else {return Wpik(6,q);}}
//7)----------------------550-600 Mev/c -------------
    if ( pm<=fcut[7][0] ){
	if ( q>fCutPr->Eval(pm) )
	  {fWp=1.;return Proton();} else {return Wpik(7,q);}}
//8)----------------------600-650 Mev/c -------------
    if ( pm<=fcut[8][0] ){
      if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
                                     else {return Wpik(8,q);};
//9)----------------------650-730 Mev/c -------------
    if ( pm<=fcut[9][0] ){
      if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}
      else {return Wpik(9,q);}}
//10)----------------------730-830 Mev/c -------------
    if( pm<=fcut[10][0] ){
      if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();}}
                                     else {return Wpik(10,q);};
//11)----------------------830-930 Mev/c -------------
    if( pm<=fcut[11][0] ){ return Wpikp(11,q); }
//12)----------------------930-1030 Mev/c -------------
    if( pm<=fcut[12][0] )
     { return Wpikp(12,q); };

    return fPcode;    
}
//-----------------------------------------------------------
void	AliTPCPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
			Float_t klo,Float_t khi,Float_t plo,Float_t phi)
{
  /// set cuts

    fcut[n][0]=pm;
    fcut[n][1]=pilo;
    fcut[n][2]=pihi;
    fcut[n][3]=klo;
    fcut[n][4]=khi;
    fcut[n][5]=plo;
    fcut[n][6]=phi;
    return ;    
}
//------------------------------------------------------------
void AliTPCPid::SetVec(Int_t ntrack,TVector info) const
{
  /// new track vector

TClonesArray& arr=*trs;
    new( arr[ntrack] ) TVector(info);
}
//-----------------------------------------------------------
TVector* AliTPCPid::GetVec(Int_t ntrack) const
{
  /// get track vector

TClonesArray& arr=*trs;
    return (TVector*)arr[ntrack];
}
//-----------------------------------------------------------
void AliTPCPid::SetEdep(Int_t track,Float_t Edep)
{
  /// energy deposit

    TVector xx(0,11);
    if( ((TVector*)trs->At(track))->IsValid() )
	{TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
    Int_t j=(Int_t)xx(0); if(j>4)return;
    xx(++j)=Edep;xx(0)=j;
    TClonesArray &arr=*trs;
    new(arr[track])TVector(xx);
}
//-----------------------------------------------------------
void AliTPCPid::SetPmom(Int_t track,Float_t Pmom)
{
  /// set part. momentum

    TVector xx(0,11);
    if( ((TVector*)trs->At(track))->IsValid() )
	{TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
    xx(10)=Pmom;
    TClonesArray &arr=*trs;
    new(arr[track])TVector(xx);
}
//-----------------------------------------------------------
void AliTPCPid::SetPcod(Int_t track,Int_t partcode)
{
  /// set part. code

    TVector xx(0,11);
    if( ((TVector*)trs->At(track))->IsValid() )
	{TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
    if(xx(11)==0)
	{xx(11)=partcode; fmxtrs++;
	TClonesArray &arr=*trs;
	new(arr[track])TVector(xx);
	}
}
//-----------------------------------------------------------
void AliTPCPid::PrintPID(Int_t track)
{
  /// control print

cout<<fmxtrs<<" tracks in AliITSPid obj."<<endl;
    if( ((TVector*)trs->At(track))->IsValid() )
	{TVector xx( *((TVector*)trs->At(track)) );
	 xx.Print();
	 }
    else 
	{cout<<"No data for track "<<track<<endl;return;}
}
//-----------------------------------------------------------
void AliTPCPid::Tab(void)
{
  /// fill table

if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
cout<<"------------------------------------------------------------------------"<<endl;
cout<<"Nq"<<"   q1  "<<"   q2  "<<"   q3  "<<"   q4  "<<"   q5   "<<
      " Qtrm    "    <<"  Wpi  "<<"  Wk   "<<"  Wp  "<<"Pmom  "<<endl;
cout<<"------------------------------------------------------------------------"<<endl;
for(Int_t i=0;i<trs->GetEntries();i++)
{
  TVector xx( *((TVector*)trs->At(i)) );     
    if( xx.IsValid() && xx(0)>0 )
	{
	    if(xx(0)>=2)
		{
//       1)Calculate Qtrm	
		    xx(6)=(this->Qtrm(i));

	    	 }else{
		     xx(6)=xx(1);
		 }
//	 2)Calculate Wpi,Wk,Wp
  	    this->GetPcode(xx(6),xx(10)/1000.);
	    xx(7)=GetWpi();
	    xx(8)=GetWk();
	    xx(9)=GetWp();
//       3)Print table
	    if(xx(0)>0){
		    cout<<xx(0)<<" ";
		    for(Int_t j=1;j<11;j++){
		        cout.width(7);cout.precision(5);cout<<xx(j);
		    }
		    cout<<endl;
		}
//	  4)Update data in TVector
	    TClonesArray &arr=*trs;
	    new(arr[i])TVector(xx);	 
	}
    else 
      {/*cout<<"No data for track "<<i<<endl;*/}
}// End loop for tracks
}
void AliTPCPid::Reset(void)
{
  /// reset

  for(Int_t i=0;i<trs->GetEntries();i++){
    TVector xx(0,11);
    TClonesArray &arr=*trs;
    new(arr[i])TVector(xx);
  }
}
//-----------------------------------------------------------
AliTPCPid::AliTPCPid(Int_t ntrack):TObject(),
		  fCutKa(0),
		  fCutPr(0),
		  fCutKaTune(0.),
		  fCutPrTune(0.),
		  fSigmin(0.),
		  fSilent(0),
		  fmxtrs(0),
		  trs(0),
		  fqtot(0.),
		  fWpi(0.),
		  fWk(0.),
		  fWp(0.),
		  fRpik(0.),
		  fRppi(0.),
		  fRpka(0.),
		  fRp(0.),
		  fPcode(0)
{
    trs = new TClonesArray("TVector",ntrack);
    TClonesArray &arr=*trs;
    for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);


    //fCutKa = new TF1("fkaons","[0]/x/x+[1]",0.1,1.2);
    //fCutPr = new TF1("fprotons","[0]/x/x +[1]",0.2,1.2);
    //TF1 *frmska=0;
    
    TF1 frmska("x_frmska","1.46-7.82*x+16.78*x^2-15.53*x^3+5.24*x^4 ",
	        0.1,1.2);
    fCutKa = new TF1("fkaons",
	   "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",
	   0.1,1.2);
    fCutPr = new TF1("fprotons",
                     "0.83*(15.32-56.094*x+89.962*x^2-66.1856*x^3+18.4052*x^4)",
		     0.2,1.2); 
    fCutKaTune=1.1; // 0.92; 
    fCutPrTune=1.0; //0.80;
    
const int kinf=10;
//         Ncut Pmom   pilo  pihi    klo    khi     plo    phi
//       fcut[j] [0]    [1]    [2]    [3]    [4]     [5]    [6]
//----------------------------------------------------------------
    SetCut(  1, 0.120,  0.  ,  0.  , kinf  , kinf   , kinf  , kinf  );
    SetCut(  2, 0.200,  0.  ,  6.0 , 6.0  , kinf   , kinf  , kinf  ); //120-200
    SetCut(  3, 0.400,  0.  ,  3.5 , 3.5  , 9.0   , 9.0  , kinf  ); //200-400
    SetCut(  4, 0.450,  0.  ,  1.9 , 1.9  , 4.0   , 4.0  , kinf  ); //400-450
//----------------------------------------------------------------
    SetCut(  5, 0.500, 0.976, 0.108, 1.484 , 0.159  , 3.5  , kinf  );  //450-500
    SetCut(  6, 0.550, 0.979, 0.108, 1.376 , 0.145  , 3.0  , kinf  );  //500-550
//----------------------------------------------------------------    
    SetCut(  7, 0.600, 0.984, 0.111, 1.295 , 0.146 , 2.7  , kinf  );   //550-600
    SetCut(  8, 0.650, 0.989, 0.113, 1.239 , 0.141 , 2.5  , kinf  );   //600-650
    SetCut(  9, 0.730, 0.995, 0.109, 1.172 , 0.132 , 2.0  , kinf  );   //650-730
//----------------------------------------------------------------    
    SetCut( 10, 0.830, 1.008, 0.116, 1.117 , 0.134 , 1.703, 0.209 ); //730-830
    SetCut( 11, 0.930, 1.019, 0.115, 1.072 , 0.121 , 1.535, 0.215 ); //830-930
    SetCut( 12, 1.230, 1.035, 0.117, 1.053 , 0.140  ,1.426, 0.270); //930-1030
    //------------------------ pi,K ---------------------
    faprob[0][0]=33219.;   faprob[1][0]=1971.;   // faprob[0][i] -    pions 
    faprob[0][1]=28828.;   faprob[1][1]=1973.; // faprob[1][i] -    kaons
    //---------------------------------------------------
    faprob[0][2]=24532.;   faprob[1][2]=1932.; faprob[2][2]=1948.;
    faprob[0][3]=20797.;   faprob[1][3]=1823.; faprob[2][3]=1970.;   
    faprob[0][4]=27017.;   faprob[1][4]=2681.; faprob[2][4]=2905.;
    //------------------------ pi,K,p -------------------
    faprob[0][5]= 24563.;    faprob[1][5]=2816.;  faprob[2][5]=3219.;  
    faprob[0][6]= 16877.;    faprob[1][6]=2231.;  faprob[2][6]=2746.;
    faprob[0][7]= 11557.;    faprob[1][7]=1681;   faprob[2][7]=2190.;

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