[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
  //
  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;
    
    frmska = new TF1("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	//
	75	fSigmin=param.fSigmin;
	76	return (*this);
	77	}
b8def77a	78	//-----------------------------------------------------------
ad350af6	79	Float_t AliTPCPid::Qtrm(Int_t track) const
b8def77a	80	{
ad350af6	81	//
	82	// dummy comment (Boris!!!)
	83	//
b8def77a	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));
ad350af6	107	}// --- End AliTPCPid::Qtrm ---
b8def77a	108
ad350af6	109	Float_t AliTPCPid::Qtrm(Float_t qarr[6],Int_t narr)
b8def77a	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;
ad350af6	132	}// --- End AliTPCPid::Qtrm ---
b8def77a	133
ad350af6	134	Int_t AliTPCPid::Wpik(Int_t nc,Float_t q)
b8def77a	135	{
ad350af6	136	//
	137	// pi-k
	138	//
b8def77a	139	Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
b8def77a	140	Float_t appi,apk;
ad350af6	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];
b8def77a	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;
9dec0e81	153	dpi =TMath::Abs(dqpi);
9dec0e81	154	dk =TMath::Abs(dqk);
b8def77a	155	Double_t dn=appiTMath::Gaus(q,qmpi,sigpi)+apkTMath::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;
ad350af6	159	}else{fWpi=1;return Pion();}
b8def77a	160	Float_t rpik=ppi/(pk+0.0000001);
b8def77a	161	if( !fSilent )
ad350af6	162	cout<<"q,dqpi,dqk, Wpik: ppi,pk,rpik="
b8def77a	163	<<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
	164
	165	fWp=0.; fWpi=ppi; fWk=pk;
ad350af6	166	if( pk>ppi){return Kaon();}else{return Pion();}
b8def77a	167	}
b8def77a	168	//-----------------------------------------------------------
ad350af6	169	Int_t AliTPCPid::Wpikp(Int_t nc,Float_t q)
b8def77a	170	{
ad350af6	171	//
	172	// pi-k-p
	173	//
b8def77a	174	Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
	175	Float_t appi,apk,app;
	176
ad350af6	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];
b8def77a	183
b8def77a	184	//appi = apk = app = 1.;
ad350af6	185	appi = faprob[0][nc-5];
	186	apk = faprob[1][nc-5];
	187	app = faprob[2][nc-5];
b8def77a	188
	189	Double_t dn= appi*TMath::Gaus(q,qmpi,sigpi)
	190	+apkTMath::Gaus(q,qmk,sigk)+appTMath::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	}
ad350af6	196	else{fWpi=1;return Pion();}
b8def77a	197	fWp=pp; fWpi=ppi; fWk=pk;
b8def77a	198	if( !fSilent ){
ad350af6	199	cout<<" Wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<"; "<<qmk<<" "<<sigk<<"; "
b8def77a	200	<<qmp<<" "<<sigp<<"; "<<endl;
ad350af6	201	cout<<" faprob: "<<appi<<" "<<apk<<" "<<app<<endl;
b8def77a	202	cout<<" ppi,pk,pp="<<ppi<<" "<<pk<<" "<<pp<<endl;
b8def77a	203	}
ad350af6	204	if( ppi>pk&&ppi>pp ) { return Pion(); }
ad350af6	205	if(pk>pp){return Kaon();}else{return Proton();}
b8def77a	206	}
b8def77a	207	//-----------------------------------------------------------
ad350af6	208	Int_t AliTPCPid::GetPcode(TClonesArray* /rps/,Float_t /pm/) const
b8def77a	209	{
ad350af6	210	//
	211	// dummy ???
	212	//
b8def77a	213	return 0;
	214	}
	215	//-----------------------------------------------------------
	216	Int_t AliTPCPid::GetPcode(AliTPCtrack *track)
	217	{
ad350af6	218	//
	219	// get particle code
	220	//
b8def77a	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]);
ad350af6	226	Float_t pt1=TMath::Abs(par[4]);
ad350af6	227	Float_t mom=1./(pt1*TMath::Cos(lam));
b8def77a	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	}
c84a5e9e	233	/*
b8def77a	234	//-----------------------------------------------------------
c300c64c	235	Int_t AliTPCPid::GetPcode(AliKalmanTrack *track)
b8def77a	236	{
ad350af6	237	//
	238	// get particle code
	239	//
b8def77a	240	if(track==0)return 0;
	241	// track->Propagate(track->GetAlpha(),3.,0.1/65.191.848,0.11.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]);
ad350af6	246	Float_t pt1=TMath::Abs(par[4]);
b8def77a	247	Float_t mom=0.;
ad350af6	248	if( (pt1TMath::Cos(lam))!=0. ){ mom=1./(pt1TMath::Cos(lam)); }else{mom=0.;};
bbc6cd2c	249	Float_t dedx=track->GetPIDsignal();
ad350af6	250	cout<<"lam,pt1,mom,dedx="<<lam<<","<<pt1<<","<<mom<<","<<dedx<<endl;
b8def77a	251	Int_t pcode=GetPcode(dedx,mom);
	252	cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
	253	return pcode?pcode:211;
	254	}
c84a5e9e	255	*/
b8def77a	256	//-----------------------------------------------------------
	257	Int_t AliTPCPid::GetPcode(Float_t q,Float_t pm)
	258	{
ad350af6	259	//
	260	// get particle code
	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	{
ad350af6	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;
b8def77a	333	return ;
	334	}
	335	//------------------------------------------------------------
ad350af6	336	void AliTPCPid::SetVec(Int_t ntrack,TVector info) const
b8def77a	337	{
ad350af6	338	//
	339	// new track vector
	340	//
b8def77a	341	TClonesArray& arr=*trs;
	342	new( arr[ntrack] ) TVector(info);
	343	}
	344	//-----------------------------------------------------------
ad350af6	345	TVector* AliTPCPid::GetVec(Int_t ntrack) const
b8def77a	346	{
ad350af6	347	//
	348	// get track vector
	349	//
b8def77a	350	TClonesArray& arr=*trs;
	351	return (TVector*)arr[ntrack];
	352	}
	353	//-----------------------------------------------------------
	354	void AliTPCPid::SetEdep(Int_t track,Float_t Edep)
	355	{
ad350af6	356	//
	357	// energy deposit
	358	//
b8def77a	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	{
ad350af6	370	//
	371	// set part. momentum
	372	//
b8def77a	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	{
ad350af6	383	//
	384	// set part. code
	385	//
b8def77a	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)
ad350af6	390	{xx(11)=partcode; fmxtrs++;
b8def77a	391	TClonesArray &arr=*trs;
	392	new(arr[track])TVector(xx);
	393	}
	394	}
	395	//-----------------------------------------------------------
7360a026	396	void AliTPCPid::PrintPID(Int_t track)
ad350af6	397	{
	398	//
	399	// control print
	400	//
	401	cout<<fmxtrs<<" tracks in AliITSPid obj."<<endl;
b8def77a	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	{
ad350af6	412	//
	413	// fill table
	414	//
b8def77a	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	{
b8def77a	425	if(xx(0)>=2)
	426	{
	427	// 1)Calculate Qtrm
ad350af6	428	xx(6)=(this->Qtrm(i));
b8def77a	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	{
ad350af6	456	//
	457	// reset
	458	//
b8def77a	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	//-----------------------------------------------------------
179c6296	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)
b8def77a	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);
179c6296	488
b8def77a	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);
ad350af6	492	TF1 *frmska=0;
176aff27	493
ad350af6	494	frmska = new TF1("x_frmska","1.46-7.82x+16.78x^2-15.53x^3+5.24x^4 ",
b8def77a	495	0.1,1.2);
	496	fCutKa = new TF1("fkaons",
	497	"1.25+0.044/x/x+1.25+0.044x-13.87x^2+22.37x^3-10.05x^4-2.5*x_frmska",
	498	0.1,1.2);
	499	fCutPr = new TF1("fprotons",
	500	"0.83(15.32-56.094x+89.962x^2-66.1856x^3+18.4052*x^4)",
	501	0.2,1.2);
	502	fCutKaTune=1.1; // 0.92;
	503	fCutPrTune=1.0; //0.80;
	504
ad350af6	505	const int kinf=10;
b8def77a	506	// Ncut Pmom pilo pihi klo khi plo phi
ad350af6	507	// fcut[j] [0] [1] [2] [3] [4] [5] [6]
b8def77a	508	//----------------------------------------------------------------
ad350af6	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
b8def77a	513	//----------------------------------------------------------------
ad350af6	514	SetCut( 5, 0.500, 0.976, 0.108, 1.484 , 0.159 , 3.5 , kinf ); //450-500
ad350af6	515	SetCut( 6, 0.550, 0.979, 0.108, 1.376 , 0.145 , 3.0 , kinf ); //500-550
b8def77a	516	//----------------------------------------------------------------
ad350af6	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
b8def77a	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 ---------------------
ad350af6	525	faprob[0][0]=33219.; faprob[1][0]=1971.; // faprob[0][i] - pions
ad350af6	526	faprob[0][1]=28828.; faprob[1][1]=1973.; // faprob[1][i] - kaons
b8def77a	527	//---------------------------------------------------
ad350af6	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.;
b8def77a	531	//------------------------ pi,K,p -------------------
ad350af6	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.;
b8def77a	535
	536	fSilent=kTRUE;
	537	}
	538	//-----------------------------------------------------------
	539
	540
	541