//////////////////////////////////////////////////////
/////////// 0 = pion
/////////// 1 = kaon
/////////// 2 = proton
/////////// 3 = electron
/////////////////////////////////////////////////////
#include <math.h>
#include "stdlib.h"
#include "T49Prob.h"
#include "T49EventRoot.h"
#include "T49ParticleRoot.h"
#include "T49Run.h"
#include "PhysConst.h"
//*-*-* ROOT includes
#include "TROOT.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TNtuple.h"
#include "TBranch.h"
#include "TObjArray.h"
#include "TObject.h"
#include "TArrayF.h"
#include "TMatrix.h"
ClassImp(T49Prob)
/*
$Log: T49Prob.C,v $
Revision 1.5 2001/11/12 12:12:03 cblume
Update for ROOT v3.01
* Revision 1.4 1999/11/23 13:54:10 cblume
* Remove files
*
*/
T49Prob::T49Prob()
{
fDedxMatrix = new TMatrix(100,20);
fDedxFitPosMatrix = new TMatrix(12,5);
fDedxFitNegMatrix = new TMatrix(12,5);
}
void T49Prob::ReadDedxPara(Char_t fname[256])
{
FILE *fp;
Int_t i;
Int_t ii;
Int_t row;
Int_t column;
Double_t x;
fp = fopen(fname,"r");
if(fp==NULL)
{
fprintf(stderr,"Error: open cutfile %s\n",fname);
}
printf("read matrix: %s ",fname);
fscanf(fp,"%d",&row);
printf("row: %d ,",row);
fscanf(fp,"%d",&column);
printf("column: %d \n",column);
if (fDedxMatrix) delete fDedxMatrix;
fDedxMatrix = new TMatrix(row,column);
for(i=0;i<row;i++)
{
for(ii=0;ii<column;ii++)
{
fscanf(fp,"%le",&x);
fDedxMatrix->operator()(i,ii) = x;
}
}
}
TArrayD T49Prob::GetAmplitude(Int_t charge,Float_t p)
{
Int_t i;
Int_t ii;
Int_t a;
Double_t A[4];
Double_t amplitude;
TArrayD ampl(4);
for(i=0;i<4;i++)
{
A[i] = 0;
if(charge == 1)
{
for(ii=0;ii<fDedxFitPosMatrix->GetNcols();ii++)
{
a= i*3 ;
A[i] = A[i] + fDedxFitPosMatrix->operator()(a,ii)*pow(p,ii);
}
}
if(charge == -1)
{
for(ii=0;ii<fDedxFitNegMatrix->GetNcols();ii++)
{
a= i*3 ;
A[i] = A[i] + fDedxFitNegMatrix->operator()(a,ii)*pow(p,ii);
}
}
amplitude = A[i];
ampl.AddAt(amplitude,i);
}
return ampl;
}
TArrayD T49Prob::GetMean(Int_t charge,Float_t p)
{
Int_t i;
Int_t ii;
Int_t a;
Double_t M[4];
Double_t meanvalue;
TArrayD mean(4);
for(i=0;i<4;i++)
{
M[i] = 0;
if(charge == 1)
{
for(ii=0;ii<fDedxFitPosMatrix->GetNcols();ii++)
{
a= i*3 ;
M[i] = M[i] + fDedxFitPosMatrix->operator()(a+1,ii)*pow(p,ii);
}
}
if(charge == -1)
{
for(ii=0;ii<fDedxFitNegMatrix->GetNcols();ii++)
{
a= i*3 ;
M[i] = M[i] + fDedxFitNegMatrix->operator()(a+1,ii)*pow(p,ii);
}
}
meanvalue = M[i];
mean.AddAt(meanvalue,i);
}
return mean;
}
TArrayD T49Prob::GetSigma(Int_t charge,Float_t p)
{
Int_t i;
Int_t ii;
Int_t a;
Double_t S[4];
Double_t sigma;
TArrayD sig(4);
for(i=0;i<4;i++)
{
S[i] = 0;
if(charge == 1)
{
for(ii=0;ii<fDedxFitPosMatrix->GetNcols();ii++)
{
a= i*3 ;
S[i] = S[i] + fDedxFitPosMatrix->operator()(a+2,ii)*pow(p,ii);
}
}
if(charge == -1)
{
for(ii=0;ii<fDedxFitNegMatrix->GetNcols();ii++)
{
a= i*3 ;
S[i] = S[i] + fDedxFitNegMatrix->operator()(a+2,ii)*pow(p,ii);
}
}
sigma = S[i];
sig.AddAt(sigma,i);
}
return sig;
}
TArrayD T49Prob::GetDedxAmplitude(Int_t charge,Float_t p,Float_t dedx)
{
Int_t i;
TArrayD A(4);
TArrayD M(4);
TArrayD S(4);
Double_t amplitude;
TArrayD ampl(4);
A = GetAmplitude(charge,p);
M = GetMean(charge,p);
S = GetSigma(charge,p);
for(i=0;i<4;i++)
{
amplitude = A.At(i)*exp(-1*((dedx-M.At(i))*(dedx-M.At(i)))/(2*S.At(i)*S.At(i)));
ampl.AddAt(amplitude,i);
}
return ampl;
}
TArrayD T49Prob::GetDedxProbability(Int_t charge,Float_t p,Float_t dedx)
{
Int_t i;
TArrayD ampl;
Double_t asum;
Double_t prob;
TArrayD DedxProb(5);
ampl = GetDedxAmplitude(charge,p,dedx);
asum = ampl.At(0)+ ampl.At(1)+ ampl.At(2)+ ampl.At(3);
for(i=0;i<4;i++)
{
prob = ampl.At(i)/asum;
DedxProb.AddAt(prob,i);
}
DedxProb.AddAt(asum,4);
return DedxProb;
}
Double_t T49Prob::GetProtonAmplitude(Int_t charge,Float_t p)
{
TArrayD ampl;
ampl = GetAmplitude(charge,p);
return ampl.At(2);
}
Double_t T49Prob::GetKaonAmplitude(Int_t charge,Float_t p)
{
TArrayD ampl;
ampl = GetAmplitude(charge,p);
return ampl.At(1);
}
Double_t T49Prob::GetPionAmplitude(Int_t charge,Float_t p)
{
TArrayD ampl;
ampl = GetAmplitude(charge,p);
return ampl.At(0);
}
Double_t T49Prob::GetElectronAmplitude(Int_t charge,Float_t p)
{
TArrayD ampl;
ampl = GetAmplitude(charge,p);
return ampl.At(3);
}
Double_t T49Prob::GetProtonAmplitude(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD ampl;
ampl = GetDedxAmplitude(charge,p,dedx);
return ampl.At(2);
}
Double_t T49Prob::GetKaonAmplitude(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD ampl;
ampl = GetDedxAmplitude(charge,p,dedx);
return ampl.At(1);
}
Double_t T49Prob::GetPionAmplitude(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD ampl;
ampl = GetDedxAmplitude(charge,p,dedx);
return ampl.At(0);
}
Double_t T49Prob::GetElectronAmplitude(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD ampl;
ampl = GetDedxAmplitude(charge,p,dedx);
return ampl.At(3);
}
Double_t T49Prob::GetProtonProbability(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD prob;
prob = GetDedxProbability(charge,p,dedx);
return prob.At(2);
}
Double_t T49Prob::GetKaonProbability(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD prob;
prob = GetDedxProbability(charge,p,dedx);
return prob.At(1);
}
Double_t T49Prob::GetPionProbability(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD prob;
prob = GetDedxProbability(charge,p,dedx);
return prob.At(0);
}
Double_t T49Prob::GetElectronProbability(Int_t charge,Float_t p,Float_t dedx)
{
TArrayD prob;
prob = GetDedxProbability(charge,p,dedx);
return prob.At(3);
}
void T49Prob::ReadDedxFit(Char_t fname[256])
{
FILE *fp;
Int_t i;
Int_t ii;
Int_t row;
Int_t column;
Double_t x;
fp = fopen(fname,"r");
if(fp==NULL)
{
fprintf(stderr,"Error: open cutfile %s\n",fname);
}
printf("read matrix: %s ",fname);
fscanf(fp,"%d",&row);
printf("row: %d ,",row);
fscanf(fp,"%d",&column);
printf("column: %d \n",column);
//make two lists pos and neg
row = row/2;
if (fDedxFitPosMatrix) delete fDedxFitPosMatrix;
if (fDedxFitNegMatrix) delete fDedxFitNegMatrix;
fDedxFitPosMatrix = new TMatrix(row,column);
fDedxFitNegMatrix = new TMatrix(row,column);
for(i=0;i<(row);i++)
{
for(ii=0;ii<column;ii++)
{
fscanf(fp,"%le",&x);
fDedxFitPosMatrix->operator()(i,ii) = x;
}
}
for(i=0;i<row;i++)
{
for(ii=0;ii<column;ii++)
{
fscanf(fp,"%le",&x);
fDedxFitNegMatrix->operator()(i,ii) = x;
}
}
}
Double_t T49Prob::GetParaPP(Int_t id, Int_t charge, Float_t p)
{
// id = 0 electron
// id = 1 pion
// id = 2 kaon
// id = 3 proton
Float_t Par1 = 3.72;
Float_t Par2 = 1.37;
Float_t Par3 = 4.0;
Float_t Par4 = 1.51;
Float_t Par5 = -2.0;
Double_t mass[4];
Double_t bg;
Int_t part;
Float_t xx, d, beta, gamma, ff;
Float_t BMin, a, b, m;
mass[0] = 0.00000026112;
mass[1] = 0.019488;
mass[2] = 0.24374;
mass[3] = 0.879844;
part = id;
bg = p / sqrt (mass[part]);
BMin = Par1/sqrt(1.0+Par1*Par1);
a = -Par5/2.0*(1.0-BMin*BMin)/TMath::Power(BMin,(Par5+4.0));
b = Par1*Par1*(1.-BMin*BMin*(1.+2./Par5))+log(1.-BMin*BMin);
m = (Par3-Par2)/((Par4/a-b+1.)/2./log(10.)-Par2);
xx = log10(bg);
d = 0.0;
if(xx > Par2)
{
d = xx-Par2+(Par2-Par3)/m;
if(xx < Par3)
{
d = d+TMath::Power((Par3-xx),m)/(m*TMath::Power((Par3-Par2),(m-1.)));
}
}
d = d*2.*log(10.);
beta = bg/sqrt(1.+bg*bg);
gamma = bg/beta;
ff = a*TMath::Power(beta,Par5)*(b+2.*log(gamma)-beta*beta-d);
// now we calculate the truncated mean from the dE/dx
ff = ff*(1.+0.1800*log(ff)-0.0155*log(ff)*log(ff));
return ff;
}
Double_t T49Prob::GetParaPbPb(Int_t id, Int_t charge, Float_t p)
{
// id = 0 electron
// id = 1 pion
// id = 2 kaon
// id = 3 proton
Double_t mass[4];
Double_t x1, x2,p0,bg,ff,d1,d2,d3,dfq;
Int_t part;
Float_t c_par1 = 1.597546;
Float_t d_par1 = 9.8;
Float_t e_par1 = 2.38;
Float_t f_par1 = 0.2;
mass[0] = 0.00000026112;
mass[1] = 0.019488;
mass[2] = 0.24374;
mass[3] = 0.879844;
part = id;
bg = p / sqrt (mass[part]);
x1 = pow(10,(e_par1 - 1.0/3.0 * sqrt(2.0*log(10.0)/(3.0 * f_par1))));
x2 = pow(10,(e_par1 + 2.0/3.0 * sqrt(2.0*log(10.0)/(3.0 * f_par1))));
p0 = c_par1/(d_par1 + 2.0*log(10.0)*e_par1 - 1.0);
if (bg<x1)
dfq = 0;
else
{
dfq = -2.0 * log(bg) + 2.0 * log(10.0) * e_par1;
if(bg<x2)
{
d1 = 2.0/3.0*sqrt(2.0*log(10.0)/(3.0 * f_par1));
d2 = log(bg)/log(10.0);
d3 = pow((e_par1 + d1 - d2),3);
dfq -= f_par1*d3;
}
}
ff = p0*( (1+ bg*bg)/(bg*bg) * ( d_par1 + log(1+(bg*bg)) + dfq)-1.0 + 0.01 * log10(bg) );
return ff;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.