[u/mrichter/AliRoot.git] / TPC / 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.                  *
 **************************************************************************/

/* $Id$ */

#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>

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