//_HADES_CLASS_DESCRIPTION
////////////////////////////////////////////////////////////////////////////
//*-- AUTHOR : J. Markert
////////////////////////////////////////////////////////////////////////////
// HMdcDeDx
//
// This transformation class calculates the "pseudo dEdx" from t2-t1 (time above threshold)
// of all fired drift cells included in one HMdcSeg. The transformation is performed doing a
// normalization of the measured t2-t1 with impact angle and minimum distance to wire
// (MDCCAL2 parametrization) and the algorithm to normalize the single measurements and
// average over all cells included in the segment.
//-------------------------------------------------------------------------
// Calculation of dEdx:
// Float_t calcDeDx(HMdcSeg*,Float_t*,Float_t*,UChar_t*,Float_t*,UChar_t*,
// Int_t vers=2,Int_t mod=2, Bool_t useTruncMean=kTRUE, Float_t truncMeanWindow=-99.)
// calculates dedx of a MDC (or 2) segment by doing normalization for
// impact angle and distance from wire for the fired cells of
// the segment. The measurements are sorted in decreasing order,
// the mean and sigma is calculated. Afterwards the truncated mean
// is calculated and returned as dedx. Mean,sigma,number of wires before
// truncating,truncated mean,truncated mean sigma and number of wires
// after truncating can be obtained by the returned pointers.
// Different methods can be selected:
// vers==0 : Input filling from inner HMdcSeg
// vers==1 : Input filling from outer HMdcSeg
// vers==2 : Input filling from inner+outer HMdcSeg (default)
// mod==0 : calc dedx from 1st module in segment
// mod==1 : calc dedx from 2nd module in segment
// mod==2 : calc dedx from whole segment (default)
// useTruncMean: kTRUE (default) apply truncated Mean
// truncMeanWindow (unit: SIGMA RMS of mean TOT): -99 (default) use standard window
// returns -99 if nothing was calculated
//-------------------------------------------------------------------------
// Settings of the truncated Mean method:
// The truncated mean is calculated according to the set window (in sigma units)
// (setWindow(Float_t window)) arround the mean. For the calculation of the truncated
// mean only drift cells with a dEdx inside the window arround the mean (calculated from
// all drift cells of the segment) will be taken into account.
// With setMinimumWires(Int_t minwires) the algorithm can be forced to keep
// a minimum number of wires. The method will stop removing drift cells from the active
// sample if the minimum number is reached.
////////////////////////////////////////////////////////////////////////////
#include "hmdcdedx2.h"
#include "hpario.h"
#include "hdetpario.h"
#include "hmessagemgr.h"
#include "hparamlist.h"
#include "hmessagemgr.h"
#include "hmdcdetector.h"
#include "hspectrometer.h"
#include "hmatrixcategory.h"
#include "hlinearcategory.h"
#include "hrecevent.h"
#include "hmdcdef.h"
#include "hmdctrackddef.h"
#include "hmdccal1.h"
#include "hmdcseg.h"
#include "hmdchit.h"
#include "hmdcclusinf.h"
#include "hmdcclusfit.h"
#include "hmdcwirefit.h"
#include "hmdcsizescells.h"
#include "hpidphysicsconstants.h"
#include "hgeantkine.h"
#include "TStyle.h"
#include "TMath.h"
#include "TRandom.h"
#include "TH1.h"
#include <stdlib.h>
ClassImp(HMdcDeDx2)
Float_t HMdcDeDx2::MaxMdcMinDist[4] = {4.,4.,8.,9.};
HMdcDeDx2::HMdcDeDx2(const char* name,const char* title,
const char* context)
: HParCond(name,title,context)
{
//
strcpy(detName,"Mdc");
clear();
loccal.set(4,0,0,0,0);
minimumWires=5;
isInitialized=kFALSE;
setWindow(1.);
measurements.Set(MAX_WIRES);
measurements.Reset(-99);
module =2;
method =2;
useCalibration=kTRUE;
debug=kFALSE;
for(Int_t i=0;i<4;i++){
MinDistStep[i]=MaxMdcMinDist[i]/(Float_t) N_DIST;
}
AngleStep=5.;
memset(&parMax[0][0][0][0],0,sizeof(Double_t)*6*4*N_ANGLE*N_DIST);
}
HMdcDeDx2::~HMdcDeDx2()
{
// destructor
}
void HMdcDeDx2::clear()
{
hefr =0;
status=kFALSE;
resetInputVersions();
changed=kFALSE;
}
void HMdcDeDx2::printParam(TString opt)
{
// prints the parameters of HMdcDeDx2 to the screen.
cout<<"HMdcDeDx2:"<<endl;
if(opt.CompareTo("par")==0||opt.CompareTo("all")==0)
{
cout<<"par:"<<endl;
for(Int_t s=0;s<6;s++){
for(Int_t m=0;m<4;m++){
for(Int_t a=0;a<N_ANGLE;a++){
for(Int_t d=0;d<N_DIST;d++){
for(Int_t i=0;i<N_PARAM;i++){
printf("s %i m %i angle %2i dist %2i %1.15e %1.15e %1.15e %1.15en",
s,m,a,d,
par[s][m][a][d][0],par[s][m][a][d][1],par[s][m][a][d][2],par[s][m][a][d][3]);
}
}
}
}
}
}
if(opt.CompareTo("parMax")==0||opt.CompareTo("all")==0)
{
cout<<"parMax:"<<endl;
for(Int_t s=0;s<6;s++){
for(Int_t m=0;m<4;m++){
for(Int_t a=0;a<N_ANGLE;a++){
for(Int_t d=0;d<N_DIST;d++){
printf("s %i m %i angle %2i dist %2i %7.3fn",
s,m,a,d,
parMax[s][m][a][d]);
}
}
}
}
}
if(opt.CompareTo("shiftpar")==0||opt.CompareTo("all")==0)
{
cout<<"shift par:"<<endl;
for(Int_t s=0;s<6;s++){
for(Int_t m=0;m<4;m++){
for(Int_t a=0;a<N_ANGLE;a++){
for(Int_t d=0;d<N_DIST;d++){
for(Int_t i=0;i<N_PARAM;i++){
printf("s %i m %i angle %2i dist %2i %1.15e %1.15e %1.15e %1.15en",
s,m,a,d,
shiftpar[s][m][a][d][0],shiftpar[s][m][a][d][1],shiftpar[s][m][a][d][2],shiftpar[s][m][a][d][3]);
}
}
}
}
}
}
cout<<"window:"<<endl;
printf("window %7.3f n",window);
cout<<"hefr:"<<endl;
printf("hefr %7.3f n",hefr);
}
Bool_t HMdcDeDx2::init(HParIo* inp,Int_t* set)
{
// intitializes the container from an input
HDetParIo* input=inp->getDetParIo("HCondParIo");
if (input) return (input->init(this,set));
return kFALSE;
}
Bool_t HMdcDeDx2::createMaxPar(Bool_t print)
{
// create the maximum allowed value of ToT to keep
// the numerical calulations inside the range of double.
// The procedure loops over all s,m,abin,dbin to set
// the boundary automatically. This function needs
// the parameter array for the fit functions to be initialized
// before. The boundary value are stored inside the container.
Double_t dEdX;
Double_t ToT;
for(Int_t s=0;s<6;s++){
for(Int_t m=0;m<4;m++){
for(Int_t a=0;a<N_ANGLE;a++){
for(Int_t d=0;d<N_DIST;d++){
Double_t* p= &par[0][m][a][d][0];
ToT=1.;
dEdX=TMath::Power(10.,(TMath::Power((ToT-p[0])/p[1],(1./p[2]))))-p[3];
while(TMath::Finite(dEdX))
{
ToT+=1.;
dEdX=TMath::Power(10.,(TMath::Power((ToT-p[0])/p[1],(1./p[2]))))-p[3];
if(!TMath::Finite(dEdX)) {
parMax[s][m][a][d]=ToT-1;
break;
}
}
if(print)
{
cout<<"s "<<s
<<" m "<<m
<<" a "<<a
<<" d "<<d
<<" ToTmax " <<parMax[s][m][a][d]
<<endl;
}
}
}
}
}
return kTRUE;
}
Bool_t HMdcDeDx2::initContainer()
{
if(!isInitialized)
{
catcal =(HCategory*)((HRecEvent*)(gHades->getCurrentEvent()))->getCategory(catMdcCal1);
if(!catcal) { Error("initContainer()","HMdcCal1 Category not found in current Event!"); exit(1);}
sizescells=HMdcSizesCells::getExObject();
if(!sizescells)
{
Error("init()","NO HMDCSIZESCELLS CONTAINER IN INPUT!");
return kFALSE;
}
//if(!sizescells->hasChanged()) sizescells->initContainer();
isInitialized=kTRUE;
}
return kTRUE;
}
Int_t HMdcDeDx2::write(HParIo* output)
{
// writes the container to an output
HDetParIo* out=output->getDetParIo("HCondParIo");
if (out) return out->write(this);
return -1;
}
void HMdcDeDx2::putParams(HParamList* l)
{
// Puts all params of HMdcDeDx2 to the parameter list of
// HParamList (which ist used by the io);
if (!l) return;
l->addBinary("par" ,&par [0][0][0][0][0],6*4*N_ANGLE*N_DIST*N_PARAM );
l->addBinary("parMax" ,&parMax [0][0][0][0] ,6*4*N_ANGLE*N_DIST );
l->addBinary("shiftpar" ,&shiftpar[0][0][0][0][0],6*4*N_ANGLE*N_DIST*N_SHIFT_PARAM);
l->add ("window" ,window);
l->add ("hefr" ,hefr);
}
Bool_t HMdcDeDx2::getParams(HParamList* l)
{
if (!l) return kFALSE;
if (!( l->fillBinary("par" ,&par [0][0][0][0][0],6*4*N_ANGLE*N_DIST*N_PARAM ))) return kFALSE;
if (!( l->fillBinary("parMax" ,&parMax [0][0][0][0] ,6*4*N_ANGLE*N_DIST ))) return kFALSE;
if (!( l->fillBinary("shiftpar" ,&shiftpar[0][0][0][0][0],6*4*N_ANGLE*N_DIST*N_SHIFT_PARAM))) return kFALSE;
if (!( l->fill ("window" ,&window ))) return kFALSE;
if (!( l->fill ("hefr" ,&hefr ))) return kFALSE;
return kTRUE;
}
void HMdcDeDx2::sort(Int_t nHits)
{
// Puts the measurement values into decreasing order.
Int_t nhit=nHits;
if(nHits>MAX_WIRES) nhit=MAX_WIRES;
register Int_t a,b,c;
Int_t exchange;
Double_t val;
for(a=0;a<nhit-1;++a)
{
exchange = 0;
c = a;
val = measurements[a];
for(b=a+1;b<nhit;++b)
{
if(measurements[b]>val) // : < increasing :> decreasing
{
c = b;
val = measurements[b];
exchange = 1;
}
}
if(exchange)
{
measurements[c] = measurements[a];
measurements[a] = val;
}
}
}
Float_t HMdcDeDx2::calcDeDx(HMdcSeg* seg[2],
Float_t* meanold,Float_t* sigmaold,UChar_t* nwire,
Float_t* sigmanew,UChar_t* nwiretrunc,
Int_t vers,Int_t mod,
Bool_t useTruncMean,
Float_t truncMeanWindow)
{
// calculates dedx of a MDC segment by doing normalization for
// impact angle and distance from wire for the fired cells of
// the segment. The measurements are sorted in decreasing order,
// the mean and sigma is calculated. Afterwards the truncated mean
// according to the set window (in sigma units) arround the mean
// is calculated and returned as dedx. Mean,sigma,number of wires before
// truncating,truncated mean,truncated mean sigma and number of wires
// after truncating can be obtained by the returned pointers.
// Different methods can be selected:
// vers==0 : Input filling from inner HMdcSeg
// vers==1 : Input filling from outer HMdcSeg
// vers==2 : Input filling from inner+outer HMdcSeg (default)
// mod==0 : calc dedx from 1st module in segment
// mod==1 : calc dedx from 2nd module in segment
// mod==2 : calc dedx from whole segment (default)
// useTruncMean: kTRUE (default) apply truncated Mean
// truncMeanWindow (unit: SIGMA RMS of mean TOT): -99 (default) use standard window
// returns -99 if nothing was calculated
//to be shure initialized values are returned
*meanold =-1.;
*sigmaold =-1.;
*nwire = 0;
*sigmanew =-1.;
*nwiretrunc= 0;
Float_t dedx =-99.;
Float_t mean =-99.;
Float_t sigma =-99.;
UChar_t nWire = 0;
UChar_t nWireTrunc = 0;
if(vers>=0&&vers<=2) method=vers;
if(seg[0]==0&& (method==0||method==2))
{
Error("HMdcDeDx2::calcDeDx()","ZERO POINTER FOR inner HMdcSeg RECEIVED!");
return mean;
}
if(seg[1]==0&& (method==1))
{
return mean;
}
if(mod>=0&&mod<3)
{
module=mod;
}
else
{
Warning("calcDeDx()","Unknown module type! Will use both modules of Segment!");
module=2;
}
//---------------------------------------------------
// getting input
measurements.Reset(-99);
nWire=fillingInput(seg);
if(debug) {
cout<<"calcDeDx() from fillingInput(): "
<<"nw "<<((Int_t)nWire)
<<" methode "<<method
<<" module "<<module
<<endl;
}
if(nWire==0) return mean;
*nwire=nWire;
if(nWire>0) mean = TMath::Mean(nWire,measurements.GetArray(), NULL);
*meanold=mean;
//---------------------------------------------------
// calculating sigma
if(nWire>=1)
{
sigma = TMath::RMS(nWire,measurements.GetArray());
*sigmaold = sigma;
//--------------------truncating measurements outside +- x*sigma
if (useTruncMean) {
nWireTrunc = select(mean,sigma,nWire,truncMeanWindow);
} else {
sort(nWire);
nWireTrunc=nWire;
}
*nwiretrunc = nWireTrunc;
//--------------------calculating truncated mean
dedx = TMath::Mean(nWireTrunc,measurements.GetArray(), NULL);
//--------------------calculating new sigma
sigma = TMath::RMS(nWireTrunc,measurements.GetArray());
*sigmanew = sigma;
}else{
Warning("calcDeDx()","nWire <=1 : skipped %i and %i with t2<=-998 !",
ctskipmod1,ctskipmod1);
}
//---------------------------------------------------
// now calibrate mean value of ToT to dEdX
Int_t s=0;
if(seg[0]) s=seg[0]->getSec();
dedx=toTTodEdX(s,ref_MOD,ref_ANGLE,ref_DIST,dedx);
if(dedx>1000) {
if(debug){cout<<"calcDeDx() : dEdX>1000 -->"<<dedx<<endl;}
dedx=1000.;
}
//---------------------------------------------------
return dedx;
}
UChar_t HMdcDeDx2::fillingInput(HMdcSeg* seg[2])
{
// Fills array of measurements from HMdcSeg and returns the number of fired wires in Segment
// The single measurements are normalized by impact angle and distance from wire. Only the first
// wires < MAX_WIRES are accepted.
ctskipmod0=0;
ctskipmod1=0;
Float_t t1,t2;
Double_t alpha,mindist;
Double_t x1,y1,z1,x2,y2,z2;
UChar_t nWire =0;
//----------------------------------------getting input--------------------------------------------------
Int_t low=0;
Int_t up =2;
if (method==0) {low=0;up=1;} // inner seg
else if(method==1) {low=1;up=2;} // outer seg
else if(method==2) {low=0;up=2;} // both seg
for(Int_t io=low;io<up;io++)
{
if(seg[io]==0) continue;
for(Int_t l=0;l<12;l++)
{
if(module==0&&l>5) continue; // fill for 1st module only
if(module==1&&l<5) continue; // fill for 2st module only
Int_t nCell=seg[io]->getNCells(l);
for(Int_t i=0;i<nCell;i++)
{
loccal[0]=seg[io]->getSec();
Int_t ioseg= seg[io]->getIOSeg();
if(ioseg==0){
(l<6)? loccal[1]=0 : loccal[1]=1;
}else if(ioseg==1){
(l<6)? loccal[1]=2 : loccal[1]=3;
}
(l<6)? loccal[2]=l : loccal[2]=l-6;
loccal[3]=seg[io]->getCell(l,i);
calcSegPoints(seg[io],x1,y1,z1,x2,y2,z2);
(*sizescells)[loccal[0]][loccal[1]][loccal[2]].getImpact(x1,y1,z1,x2,y2,z2,loccal[3],alpha,mindist);
HMdcCal1* cal1 =(HMdcCal1*)catcal->getObject(loccal);
if(cal1!=0)
{
t1 =cal1->getTime1();
t2 =cal1->getTime2();
if(t2!=-998&&t2!=-999&&
t1!=-998&&t1!=-999&&
TMath::Finite(t1) &&
TMath::Finite(t2)
)
{ // some times t1 or t2 can be missing (-999,-998)
// calculating mean value
if(nWire<MAX_WIRES-1)
{
if(useCalibration)measurements[nWire]=normalize(loccal[0],loccal[1],alpha,mindist,t2-t1);
else measurements[nWire]=t2-t1;
nWire++;
}else{
Warning("fillingInput()","Number of wires in Segment > MAX_WIRES! Skipped!");
}
} else {
if(loccal[1]==0||loccal[1]==2)ctskipmod0++;
if(loccal[1]==1||loccal[1]==3)ctskipmod1++;}
}else{
Warning("calcDeDx()","ZERO pointer recieved for cal1 object!");
}
}
}
}
return nWire;
}
UChar_t HMdcDeDx2::select(Float_t mean,Float_t sigma,UChar_t nWire,Float_t wind)
{
// loops over the measurement array with nWire values and puts values to
// -99 if the measurements are outside the wanted window arround the mean value mean
// (set with setWindow() (sigma units)).
// The procedure keeps at least a number of minimum wires set with setMinimumWires().
UChar_t nWTrunc=nWire;
UChar_t count=0;
UChar_t minW=(UChar_t)getMinimumWires();
if(nWTrunc<minW)
{
if(debug){cout<<"select : skipp because nWT<minW "<<endl;}
return nWTrunc;
}
//---------------------------------------------------
// decide about the window for truncated Mean method
Float_t tempWindow=getWindow();
if(wind>0) tempWindow = wind;
//---------------------------------------------------
//--------------------sorting measurements in decreasing order
sort(nWire);
if(debug){
cout<<"---------------------------------------------------"<<endl;
cout<<"measurements before truncation :"<<endl;
for(Int_t i=0;i<nWire;i++){
cout<<i<<" "<<measurements[i]<<endl;
}
}
Bool_t fail_high=kFALSE;
Bool_t fail_low =kFALSE;
while(nWTrunc>minW &&
count<nWire )
{
//------------------------------------------------------------
// starting from highest val
if(!fail_high && fabs(measurements[count]-mean)>(tempWindow*sigma))
{ // if measurement outside the window
measurements[count]=-99;
nWTrunc--;
} else fail_high=kTRUE;
//------------------------------------------------------------
//------------------------------------------------------------
// taking lowest val next
if(nWTrunc>minW)
{
if(!fail_low && fabs(measurements[nWire-1-count]-mean)>(tempWindow*sigma))
{ // if measurement outside the window
measurements[nWire-1-count]=-99;
nWTrunc--;
}
else fail_low =kTRUE;
} else fail_low =kTRUE;
//------------------------------------------------------------
count++;
if(fail_low && fail_high) break;
}
sort(nWire);
if(debug){
cout<<"---------------------------------------------------"<<endl;
cout<<"measurements after truncation :"<<endl;
for(Int_t i=0;i<nWTrunc;i++){
cout<<i<<" "<<measurements[i]<<endl;
}
cout<<"method "<<method
<<" window "<<tempWindow
<<" mean "<<mean
<<" meanN "<<TMath::Mean(nWTrunc,measurements.GetArray())
<<" sig "<<sigma
<<" sigN "<<TMath::RMS(nWTrunc,measurements.GetArray())
<<" w "<<((Int_t)nWire)
<<" trunc w "<<((Int_t)nWTrunc)
<<endl;
}
return nWTrunc;
}
void HMdcDeDx2::setFuncPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t* p,Int_t size)
{
// set the values for the dedx functions by s, m, angle bin, dist bin. Needs pointer
// to parameter array[N_PARAM]
if(size==N_PARAM&&
abin>=0&&abin<N_ANGLE&&
dbin>=0&&dbin<N_DIST &&
s>=0&&s<6&&
m>=0&&m<4){
for(Int_t i=0;i<size;i++) par[s][m][abin][dbin][i]=p[i];
} else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncPar(Double_t* p)
{
// set all values for the dedx functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST][N_PARAM]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
for(Int_t i=0;i<N_PARAM;i++) {
Int_t maxS=4*N_ANGLE*N_DIST*N_PARAM;
Int_t maxM=N_ANGLE*N_DIST*N_PARAM;
Int_t maxA=N_DIST*N_PARAM;
par[s][m][a][d][i]=p[s*maxS+m*maxM+a*maxA+d*N_PARAM+i];
}
}
}
}
}
}
void HMdcDeDx2::getFuncPar(Double_t* p)
{
// set all values for the dedx functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST][N_PARAM]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
for(Int_t i=0;i<N_PARAM;i++) {
Int_t maxS=4*N_ANGLE*N_DIST*N_PARAM;
Int_t maxM=N_ANGLE*N_DIST*N_PARAM;
Int_t maxA=N_DIST*N_PARAM;
p[s*maxS+m*maxM+a*maxA+d*N_PARAM+i]=par[s][m][a][d][i];
}
}
}
}
}
}
void HMdcDeDx2::setFuncMaxPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t p)
{
// set the values for the dedx functions by s, m, angle bin, dist bin.
if(abin>=0&&abin<N_ANGLE&&
dbin>=0&&dbin<N_DIST &&
s>=0&&s<6&&
m>=0&&m<4){
parMax[s][m][abin][dbin]=p;
} else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncMaxPar(Double_t* p)
{
// set all values for the dedx functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
Int_t maxS=4*N_ANGLE*N_DIST;
Int_t maxM=N_ANGLE*N_DIST;
Int_t maxA=N_DIST;
parMax[s][m][a][d]=p[s*maxS+m*maxM+a*maxA+d];
}
}
}
}
}
void HMdcDeDx2::getFuncMaxPar(Double_t* p)
{
// set all values for the dedx functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
Int_t maxS=4*N_ANGLE*N_DIST;
Int_t maxM=N_ANGLE*N_DIST;
Int_t maxA=N_DIST;
p[s*maxS+m*maxM+a*maxA+d]=parMax[s][m][a][d];
}
}
}
}
}
void HMdcDeDx2::setFuncWidthPar(Int_t s,Int_t m,Int_t abin,Int_t dbin,Double_t* p,Int_t size)
{
// set the values fpr the dedx width functions by s, m, angle bin, dist bin. Needs pointer
// to parameter array[N_SHIFT_PARAM]
if(size==N_SHIFT_PARAM&&
abin>=0&&abin<N_ANGLE&&
dbin>=0&&dbin<N_DIST &&
s>=0&&s<6&&
m>=0&&m<4){
for(Int_t i=0;i<size;i++) shiftpar[s][m][abin][dbin][i]=p[i];
} else Error("SetFuncPar(...)","array indices out of range!");
}
void HMdcDeDx2::setFuncWidthPar(Double_t* p)
{
// set all values for the dedx width functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST][N_SHIFT_PARAM]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
for(Int_t i=0;i<N_SHIFT_PARAM;i++) {
Int_t maxS=4*N_ANGLE*N_DIST*N_SHIFT_PARAM;
Int_t maxM=N_ANGLE*N_DIST*N_SHIFT_PARAM;
Int_t maxA=N_DIST*N_SHIFT_PARAM;
shiftpar[s][m][a][d][i]=p[s*maxS+m*maxM+a*maxA+d*N_SHIFT_PARAM+i];
}
}
}
}
}
}
void HMdcDeDx2::getFuncWidthPar(Double_t* p)
{
// set all values for the dedx width functions. Needs pointer
// to parameter array[6][4][N_ANGLE][N_DIST][N_SHIFT_PARAM]
for(Int_t s=0;s<6;s++) {
for(Int_t m=0;m<4;m++) {
for(Int_t a=0;a<N_ANGLE;a++) {
for(Int_t d=0;d<N_DIST;d++) {
for(Int_t i=0;i<N_SHIFT_PARAM;i++) {
Int_t maxS=4*N_ANGLE*N_DIST*N_SHIFT_PARAM;
Int_t maxM=N_ANGLE*N_DIST*N_SHIFT_PARAM;
Int_t maxA=N_DIST*N_SHIFT_PARAM;
p[s*maxS+m*maxM+a*maxA+d*N_SHIFT_PARAM+i]=shiftpar[s][m][a][d][i];
}
}
}
}
}
}
void HMdcDeDx2::calcSegPoints(HMdcSeg * seg,Double_t& x1, Double_t& y1, Double_t& z1, Double_t& x2, Double_t& y2, Double_t& z2)
{
// calculates 2 coordinate points from segment
Double_t teta=seg->getTheta();
Double_t phi =seg->getPhi();
Double_t z =seg->getZ();
Double_t r =seg->getR();
Double_t pi2 =acos(-1.)/2.;
Double_t X=r*cos(phi+pi2);
Double_t Y=r*sin(phi+pi2);
Double_t Z=z;
x1=X;
y1=Y;
z1=Z;
x2=X+cos(phi)*sin(teta);
y2=Y+sin(phi)*sin(teta);
z2=Z+cos(teta);
}
void HMdcDeDx2::findBin(Int_t m,Double_t* angle,Double_t* mindist,Int_t* abin,Int_t* dbin)
{
Double_t a=*angle;
Double_t d=*mindist;
if(d>MaxMdcMinDist[m]) d = MaxMdcMinDist[m];
if(a>90) a=90.;
if(a<0) a= 0.;
(*abin)=(Int_t)((90.-a)/AngleStep);
(*dbin)=(Int_t)(d/MinDistStep[m]);
if((*abin)==18) (*abin)=17;
if((*dbin)==N_DIST) (*dbin)=N_DIST-1;
}
Double_t HMdcDeDx2::toTSigma(Int_t s,Int_t m,Double_t angle,Double_t mindist,Int_t shift)
{
// returns asymmetric width of ToT for single drift cell
// shift == 0 (default) : returns 0
// ==-1 : low bound width of ToT distribution (sigma)
// == 1 : upper bound width of ToT distribution (sigma)
Int_t abin=0;
Int_t dbin=0;
findBin(m,&angle,&mindist,&abin,&dbin);
// shift is half FWHM = 2.355/2=1.1775 sigma
// we have to recalulate sigma here
if (shift<0) return (shiftpar[s][m][abin][dbin][0]/1.1775);
else if (shift>0) return (shiftpar[s][m][abin][dbin][1]/1.1775);
else return 0.;
}
Double_t HMdcDeDx2::toTTodEdX(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t ToT)
{
// calculation of ToT -> dEdX for single drift cell
Int_t abin=0;
Int_t dbin=0;
findBin(m,&angle,&mindist,&abin,&dbin);
Double_t* p= &par[s][m][abin][dbin][0];
//-----------------------------------------------------
if(ToT>parMax[s][m][abin][dbin]) ToT=parMax[s][m][abin][dbin];
Double_t dEdX=TMath::Power(10.,(TMath::Power((ToT-p[0])/p[1],(1./p[2]))))-p[3];
if(!TMath::Finite(dEdX)){
Warning("totTodEdX()","dEdX out of range: s %i m %i abin %i dbin %i ToT %1.15e dEdX %1.15e !",s,m,abin,dbin,ToT,dEdX);
dEdX=-1;
}
return dEdX;
}
Double_t HMdcDeDx2::dEdXToToT(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t dEdX)
{
// calculation of dEdX -> ToT for single drift cell
Int_t abin=0;
Int_t dbin=0;
findBin(m,&angle,&mindist,&abin,&dbin);
Double_t* p= &par[s][m][abin][dbin][0];
Double_t ToT=p[0]+p[1]*TMath::Power((TMath::Log10(dEdX+p[3])),p[2]);
if(!TMath::Finite(ToT)){
Warning("dEdXToToT()","ToT out of range: s %i m %i abin %i dbin %i ToT %1.15e dEdX %1.15e !",s,m,abin,dbin,ToT,dEdX);
ToT=-1;
}
return ToT;
}
Double_t HMdcDeDx2::normalize(Int_t s,Int_t m,Double_t angle,Double_t mindist,Double_t ToT)
{
// calibrate ToT :
// dEdX=TotTodEdX(t2-t1) ----> dEdXToToT(dEdX) for reference module,angle,distbin
if(ToT<0) return ToT;
Double_t dEdX=toTTodEdX(s,m,angle,mindist,ToT);
if(dEdX>=0) ToT=dEdXToToT(s,ref_MOD,ref_ANGLE,ref_DIST,dEdX);
return ToT;
}
//----------------------dEdx-----------------------------------------------
Double_t HMdcDeDx2::energyLoss(Int_t id,Double_t p,Double_t hefr)
{
// Calculates the dEdx (MeV 1/g cm2) of an particle with GEANT ID id
// and momentum p (MeV) for He/i-butan gas mixture with He fraction hefr
// (He (hefr) / i-butan (1-hefr)). The fomular is taken from PDG and doesn't
// include the density correction term. The values for the mean excitation
// energy are taken from Sauli.
if(p==0) return -1;
if(hefr<0.||hefr>1.) return -1;
Double_t mass =HPidPhysicsConstants::mass(id);
if(mass==0) return -1;
// keep function value inside boundary
if (mass< 150 &&p<5 ) p=5 ;
else if (mass>=150 &&mass<1000&&p<20) p=20;
else if (mass>=1000&&mass<2000&&p<35) p=35;
else if (mass>=2000 &&p<55) p=55;
// Z and A
Double_t Z_gas =2.*hefr+(1.-hefr)*34.;
Double_t A_gas =4.*hefr+(1.-hefr)*58.;
// I_0 and I
Double_t I_0_gas=24.6*hefr+(1.-hefr)*10.8;
Double_t I2 =pow(I_0_gas*Z_gas*(1.e-6),2); // sauli
//Double_t I2 =pow(16.*pow(Z_gas,0.9),2); //C.Lippmann thesis
Double_t K =0.307075; // MeV cm2 PDG, 4*pi*N_{A}*r_{e}2*m_{e}2*c2
Double_t mass2 =pow(mass,2);
Double_t m_ec2 =HPidPhysicsConstants::mass(3);
Double_t z2 =pow((Float_t)HPidPhysicsConstants::charge(id),2);
Double_t p2 =pow(p,2);
Double_t beta2 =1./((mass2/p2)+1.);
Double_t gamma2 =1./(1-beta2);
Double_t gamma =sqrt(gamma2);
Double_t Tmax =(2.*m_ec2*beta2*gamma2)/(1.+ 2.*gamma*(m_ec2/mass)+pow((m_ec2/mass),2));
Double_t term1 =K*z2*(Z_gas/A_gas)*(1./beta2);
Double_t term2 =((2.*m_ec2*beta2*gamma2*Tmax)/I2);
Double_t dedx =term1 * (0.5*log(term2)-beta2);
return dedx;
}
TGraph* HMdcDeDx2::energyLossGraph(Int_t id,Double_t hefr,TString opt,Bool_t exchange,Int_t markerstyle,Int_t markercolor,Float_t markersize)
{
// creates a TGraph for particle with GEANT ID id and
// and He/i-butan gas mixture with He fraction hefr
// (He (hefr) / i-butan (1-hefr) dEdx vs p,beta,1/beta2 or betagamma
// depending on opt (p,beta,1/beta2,betagamma). exchange=kTRUE
// will exchange x-axis and y-axis.
Double_t pmin=50;
Double_t pmax=100000.;
Int_t np =100;
Double_t ptot=0;
Double_t xarray[np];
Double_t yarray[np];
Int_t vers=0;
Double_t mass=HPidPhysicsConstants::mass(id);
if(opt.CompareTo("p" )==0)vers=0;
else if(opt.CompareTo("beta" )==0)vers=1;
else if(opt.CompareTo("1/beta2" )==0)vers=2;
else if(opt.CompareTo("betagamma")==0)vers=3;
else {cout<<"HMdcDedx::calcDeDxGraph():unknow option!"<<endl;}
for(Int_t i=1;i<=np;i++)
{
ptot=pmin*pow((pmax/pmin),(i-1)/(Double_t)(np-1.));
if(vers==0){xarray[i-1]=ptot;}
if(vers==1){xarray[i-1]=sqrt(1./((pow(mass,2)/pow(ptot,2))+1.));}
if(vers==2){xarray[i-1]=((pow(mass,2)/pow(ptot,2))+1.);}
if(vers==3){xarray[i-1]=(ptot/mass);}
yarray[i-1]=HMdcDeDx2::energyLoss(id,ptot,hefr);
}
Char_t name[300];
sprintf(name,"dedx_%s_He_%5.1f_ibutane_%5.1f",HPidPhysicsConstants::pid(id),hefr*100,(1-hefr)*100);
TGraph* g=0;
if(!exchange)g=new TGraph(np,xarray,yarray);
else g=new TGraph(np,yarray,xarray);
g->SetName(name);
g->SetMarkerStyle(markerstyle);
g->SetMarkerSize (markersize);
g->SetMarkerColor(markercolor);
return g;
}
Double_t HMdcDeDx2::scaledTimeAboveThreshold(HGeantKine* kine,
Double_t p,
Float_t t1, Float_t t2, Float_t* t2err,
Int_t s,Int_t m,
Double_t alpha,Double_t mindist)
{
// calculated the t2 of drift cell measurent
// from dedx of a given particle with momentum p.
// The assumption is that all drift cell measurements
// belonging to one track can be treated independent.
// return t2 and smeared error to pointer t2err for an
// single measurement. If the result t2-t1 would
// be < 0 a random value for t2 0-20 ns larger than
// t1(already including error) will be created.
// The smeared error will be returned to the pointer t2err.
//------------------------------------------------
// check input
if(!kine) {
Error("caledTimeAboveThreshold()","retrieved ZERO pointer for kine object!");
return -99;
}
if(p<0) {
Error("caledTimeAboveThreshold()","momentum = %5.3f is wrong!",p);
return -99;
}
if(t1<-997) {
Error("caledTimeAboveThreshold()","t1 = %5.3f is wrong!",t1);
return -99;
}
if(t2<-997) {
Error("caledTimeAboveThreshold()","t2 = %5.3f is wrong!",t1);
return -99;
}
//------------------------------------------------
//------------------------------------------------
// get particle ID from kine object
Int_t pTr=-1,pID=0;
kine->getParticle(pTr,pID);
//------------------------------------------------
//------------------------------------------------
// get mass and charge of particle for dedx
// calculation
Double_t mass =HPidPhysicsConstants::mass (pID);
Int_t charge=HPidPhysicsConstants::charge(pID);
if(charge==0||mass==0)
{
Warning("HMdcDeDx2::scaledEnergyLoss()","id %i %s mass %7.5f charge %i p %7.5f skipped!"
,pID,HPidPhysicsConstants::pid(pID),mass,charge,p);
return t2;
}
//------------------------------------------------
//------------------------------------------------
// calulate analytical dedx of particle
Double_t dedx = energyLoss(pID,p,getHeFraction());
if(!TMath::Finite(dedx)){
Error("caledTimeAboveThreshold()","dedx either NaN or Inf!");
return t2;
}
//------------------------------------------------
//------------------------------------------------
// recalculate dedx to mean value of Time over
// threshold for the given drift chamber
Double_t ToTmean= dEdXToToT(s,m,alpha,mindist,dedx);
if(!TMath::Finite(ToTmean)||dedx<0){
Error("caledTimeAboveThreshold()","ToT either NaN or Inf!");
cout<<"s " <<s
<<" m " <<m
<<" angle " <<alpha
<<" dist " <<mindist
<<" dedx " <<dedx
<<" pid " <<pID
<<" mass " <<mass
<<" charge "<<charge
<<" mom " <<p
<<endl;
return t2;
}
//------------------------------------------------
//------------------------------------------------
// generate some realistic smearing arround
// the ToT mean value for single measurement
// The distribution of ToT is asymmetric.
// the smearing has to be done for left and right half
// different.
Double_t sigL = toTSigma (s,m,alpha,mindist,-1);
Double_t sigR = toTSigma (s,m,alpha,mindist, 1);
Double_t smear;
if(gRandom->Rndm()>0.5){
smear = TMath::Abs(gRandom->Gaus(0., sigR));
} else{
smear =-TMath::Abs(gRandom->Gaus(0., sigL));
}
if(t2err)*t2err=smear;
//------------------------------------------------
//------------------------------------------------
// return t2 mean of an single measurement. If the
// result t2-t1 would be < 0 a random value for
// t2 0-20 ns larger than t1 will be created.
t2 = ToTmean+t1;
if(t2<=t1) t2=t1+gRandom->Rndm()*20.;
if(debug){
cout<<"scaledEnergyLoss() "
<<" s "<<s
<<" m "<<m
<<" a "<<alpha
<<" d "<<mindist
<<" id "<<pID
<<" mass "<<mass
<<" charge "<<charge
<<" mom "<<p
<<" dedx "<<dedx
<<" t1 "<<t1
<<" t2 "<<t2
<<" t2err "<<smear
<<" tot "<<ToTmean
<<" sigL "<<sigL
<<" sigR "<<sigR
<<endl;
}
return t2;
}
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.