CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCorrecSvc Class Reference

#include <DedxCorrecSvc.h>

+ Inheritance diagram for DedxCorrecSvc:

Public Member Functions

 DedxCorrecSvc (const std::string &name, ISvcLocator *svcloc)
 
 ~DedxCorrecSvc ()
 
virtual StatusCode queryInterface (const InterfaceID &riid, void **ppvUnknown)
 
virtual StatusCode initialize ()
 
virtual StatusCode finalize ()
 
void handle (const Incident &)
 
double RungCorrec (int runNO, double ex) const
 
double WireGainCorrec (int wireid, double ex) const
 
double DriftDistCorrec (int layid, double ddrift, double ex) const
 
double SaturCorrec (int layid, double costheta, double ex) const
 
double CosthetaCorrec (double costheta, double ex) const
 
double T0Correc (double t0, double ex) const
 
double HadronCorrec (double costheta, double dedx) const
 
double ZdepCorrec (int layid, double zhit, double ex) const
 
double EntaCorrec (int layid, double enta, double ex) const
 
double LayerGainCorrec (int layid, double ex) const
 
double DocaSinCorrec (int layid, double doca, double enta, double ex) const
 
double DipAngCorrec (int layid, double costheta, double ex) const
 
double GlobalCorrec (double dedx) const
 
double CellCorrec (int ser, double adc, double dd, double enta, double z, double costheta) const
 
double LayerCorrec (int layer, double z, double costheta, double ex) const
 
double TrkCorrec (double costheta, double dedx) const
 
double StandardCorrec (int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta) const
 
double StandardHitCorrec (int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const
 
double StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const
 
double PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const
 
double D2I (const double &cosTheta, const double &D) const
 
double I2D (const double &cosTheta, const double &I) const
 
void set_flag (int calib_F)
 
double f_larcos (double x, int nbin)
 
virtual double RungCorrec (int runNO, double ex) const =0
 
virtual double WireGainCorrec (int wireid, double ex) const =0
 
virtual double DriftDistCorrec (int layid, double ddrift, double ex) const =0
 
virtual double SaturCorrec (int layid, double costheta, double ex) const =0
 
virtual double EntaCorrec (int layid, double enta, double ex) const =0
 
virtual double ZdepCorrec (int layer, double z, double dedx) const =0
 
virtual double LayerGainCorrec (int layid, double dedx) const =0
 
virtual double GlobalCorrec (double dedx) const =0
 
virtual double CellCorrec (int ser, double adc, double dd, double enta, double z, double theta) const =0
 
virtual double LayerCorrec (int layer, double z, double costheta, double ex) const =0
 
virtual double TrkCorrec (double costheta, double dedx) const =0
 
virtual double StandardCorrec (int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta) const =0
 
virtual double StandardHitCorrec (int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0
 
virtual double StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const =0
 
virtual double PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const =0
 
virtual void set_flag (int calib_F)=0
 

Additional Inherited Members

- Static Public Member Functions inherited from IDedxCorrecSvc
static const InterfaceID & interfaceID ()
 

Detailed Description

Definition at line 19 of file DedxCorrecSvc.h.

Constructor & Destructor Documentation

◆ DedxCorrecSvc()

DedxCorrecSvc::DedxCorrecSvc ( const std::string &  name,
ISvcLocator *  svcloc 
)

Definition at line 39 of file DedxCorrecSvc.cxx.

39 :
40 geosvc(0),Service (name, svcloc) {
41 declareProperty("Run",m_run=1);
42 declareProperty("init",m_init=1);
43 declareProperty("par_flag",m_par_flag=0);
44 declareProperty("Debug",m_debug=false);
45 declareProperty("DebugI",m_debug_i=1);
46 declareProperty("DebugJ",m_debug_j=1);
47 m_initConstFlg = false;
48 // cout<<"DedxCorrecSvc::DedxCorrecSvc"<<endl;
49 }

◆ ~DedxCorrecSvc()

DedxCorrecSvc::~DedxCorrecSvc ( )

Definition at line 51 of file DedxCorrecSvc.cxx.

51 {
52}

Member Function Documentation

◆ CellCorrec()

double DedxCorrecSvc::CellCorrec ( int  ser,
double  adc,
double  dd,
double  enta,
double  z,
double  costheta 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 842 of file DedxCorrecSvc.cxx.

843 {
844
845 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
846 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
847 m_ggs_flag == 0) return adc;
848 adc = m_valid[ser]*m_wire_g[ser]*adc;
849 // int lyr = Wire2Lyr(ser);
850 int lyr = ser;
851 double ex = adc;
852 double correct1_ex , correct2_ex, correct3_ex, correct4_ex, correct5_ex;
853
854 if( m_ggs_flag) {
855 correct1_ex = SaturCorrec( lyr, costheta, adc );
856 ex = correct1_ex;
857 } else {
858 correct1_ex = adc;
859 }
860
861 if( m_ddg_flag) {
862 correct2_ex = DriftDistCorrec( lyr, dd, correct1_ex );
863 ex = correct2_ex;
864 } else {
865 correct2_ex =correct1_ex;
866 }
867 if( m_enta_flag) {
868 correct3_ex = DriftDistCorrec( lyr, sinenta, correct2_ex );
869 ex = correct3_ex;
870 } else {
871 correct3_ex =correct2_ex;
872 }
873
874 if( m_zdep_flag) {
875 correct4_ex = ZdepCorrec( lyr, z, correct3_ex );
876 ex = correct4_ex;
877 } else {
878 correct4_ex = correct3_ex;
879 }
880
881 if( m_layerg_flag) {
882 correct5_ex = LayerGainCorrec( lyr, correct4_ex );
883 ex = correct5_ex;
884 } else {
885 correct5_ex = correct4_ex;
886 }
887 return ex;
888
889}
double LayerGainCorrec(int layid, double ex) const
double SaturCorrec(int layid, double costheta, double ex) const
double DriftDistCorrec(int layid, double ddrift, double ex) const
double ZdepCorrec(int layid, double zhit, double ex) const

◆ CosthetaCorrec()

double DedxCorrecSvc::CosthetaCorrec ( double  costheta,
double  ex 
) const

Definition at line 565 of file DedxCorrecSvc.cxx.

565 {
566 //cout<<"DedxCorrecSvc::CosthetaCorrec"<<endl;
567 double dedx_cos;
568 //cout<<"costheta vaule is = "<<costheta<< " dedx = " << dedx << endl;
569 if(fabs(costheta)>1.0) return dedx;
570
571 const int nbin=cos_k.size()+1;
572 const double step=2.00/nbin;
573 int n= (int)((costheta+1.00-0.5*step)/step);
574 if(costheta>(1.00-0.5*step))
575 n=nbin-2;
576
577 if(costheta>-0.5*step && costheta<=0)
578 dedx_cos=cos_k[n-1]*costheta+cos_b[n-1];
579 else if(costheta>0 && costheta<0.5*step)
580 dedx_cos=cos_k[n+1]*costheta+cos_b[n+1];
581 else
582 dedx_cos=cos_k[n]*costheta+cos_b[n];
583
584 //cout << "cos_k[n]="<<cos_k[n] << " cos_b[n]=" << cos_b[n] <<
585 // " dedx_cos=" << dedx_cos << " final dedx=" << dedx/dedx_cos << endl;
586
587 //conside physical edge around 0.93
588 if(nbin==80){ // temporally for only nbin = 80
589 if(costheta>0.80 && costheta<0.95){
590 n = 77;
591 if(costheta<0.9282) n--;
592 if(costheta<0.9115) n--;
593 if(costheta<0.8877) n--;
594 if(costheta<0.8634) n--;
595 if(costheta<0.8460) n--;
596 if(costheta<0.8089) n--;
597 dedx_cos=cos_k[n]*costheta+cos_b[n];
598 }
599 else if(costheta<-0.80 && costheta>-0.95){
600 n = 2;
601 if(costheta>-0.9115) n++;
602 if(costheta>-0.8877) n++;
603 if(costheta>-0.8634) n++;
604 if(costheta>-0.8460) n++;
605 if(costheta>-0.8089) n++;
606 dedx_cos=cos_k[n]*costheta+cos_b[n];
607 }
608 }
609
610 if(dedx_cos>0){
611 dedx_cos = dedx/dedx_cos;
612 return dedx_cos;
613 }
614 else return dedx;
615}
const Int_t n

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ D2I()

double DedxCorrecSvc::D2I ( const double &  cosTheta,
const double &  D 
) const

Definition at line 2153 of file DedxCorrecSvc.cxx.

2154{
2155 //cout<<"DedxCorrecSvc::D2I"<<endl;
2156 double absCosTheta = fabs(cosTheta);
2157 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2158 double chargeDensity = D/projection;
2159 double numerator = 1 + m_alpha*chargeDensity;
2160 double denominator = 1 + m_gamma*chargeDensity;
2161 double I = D;
2162
2163 //if(denominator > 0.1)
2164
2165 I = D * m_ratio* numerator/denominator;
2166// cout << "m_alpha " << m_alpha << endl;
2167// cout << "m_gamma " << m_gamma << endl;
2168// cout << "m_delta " << m_delta << endl;
2169// cout << "m_power " << m_power << endl;
2170// cout << "m_ratio " << m_ratio << endl;
2171 return I;
2172}
const DifComplex I

Referenced by HadronCorrec().

◆ DipAngCorrec()

double DedxCorrecSvc::DipAngCorrec ( int  layid,
double  costheta,
double  ex 
) const

Definition at line 828 of file DedxCorrecSvc.cxx.

828 {
829}

Referenced by StandardCorrec().

◆ DocaSinCorrec()

double DedxCorrecSvc::DocaSinCorrec ( int  layid,
double  doca,
double  enta,
double  ex 
) const

Definition at line 785 of file DedxCorrecSvc.cxx.

785 {
786 if(m_debug) cout<<"DedxCorrecSvc::DocaSinCorrec"<<endl;
787 // cout<<"doca = "<<doca<<" eangle = "<<eangle<<endl;
788
789 if(eangle>0.25) eangle = eangle -0.5;
790 else if(eangle < -0.25) eangle = eangle +0.5;
791 int iDoca;
792 int iEAng = 0;
793 doca = doca/doca_norm[layer];
794 iDoca =(Int_t)floor((doca-Out_DocaMin)/Out_Stepdoca);
795 if(doca<Out_DocaMin) iDoca=0;
796 if(doca>Out_DocaMax) iDoca=NumSlicesDoca-1;
797 if(iDoca >=(Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca) )
798 iDoca = (Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca)-1;
799 if(iDoca<0) iDoca = 0;
800 if(m_debug) cout<<"iDoca : "<<iDoca<<" doca : "<<doca<<" Out_DocaMin : "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<endl;
801
802 //iEAng = (Int_t)floor((eangle-Phi_Min)/IEangle_step);
803 for(int i =0; i<14; i++){
804 //iEAng =0;
805 if(eangle <Eangle_value_cut[i] || eangle > Eangle_value_cut[i+1]) continue;
806 else if(eangle>= Eangle_value_cut[i] && eangle < Eangle_value_cut[i+1]) {
807 for(int k =0; k<i; k++){
808 iEAng+= Eangle_cut_bin[k];
809 }
810 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
811 int temp_bin = int((eangle-Eangle_value_cut[i])/eangle_bin_cut_step);
812 iEAng+= temp_bin;
813 }
814 }
815 //cout<<iDoca <<" "<<iEAng<<endl;
816 if(m_docaeangle[iDoca][iEAng]>0) {
817 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
818 cout << "doca: " << doca << " eang" << eangle << " dedx before: " << dedx << endl;
819 dedx = dedx/m_docaeangle[iDoca][iEAng];
820 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
821 cout << "gain: " << m_docaeangle[iDoca][iEAng] << " dedx after: " << dedx << endl;
822 }
823 return dedx;
824}
#define NumSlicesDoca
#define Out_DocaMin
#define Out_DocaMax

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ DriftDistCorrec()

double DedxCorrecSvc::DriftDistCorrec ( int  layid,
double  ddrift,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 694 of file DedxCorrecSvc.cxx.

694 {
695 // cout<<"DedxCorrecSvc::DriftDistCorrec"<<endl;
696 double dedx_ddg;
697 //cout<<"m_par_flag = "<<m_par_flag<<endl;
698 //cout<<"dd vaule is = "<<dd<<endl;
699 dd = fabs(dd);
700 if(m_par_flag == 1) {
701 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*dd +
702 m_ddg[2][layer]*dd*dd + m_ddg[3][layer]*pow(dd,3);
703 } else if(m_par_flag == 0) {
704 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*T1(dd) +
705 m_ddg[2][layer]*T2(dd) + m_ddg[3][layer]*T3(dd);
706 }
707 //cout<<"dedx_ddg is = "<<dedx_ddg<<endl;
708 dedx_ddg = dedx/dedx_ddg;
709 if (dedx_ddg>0.0) return dedx_ddg;
710 else return dedx;
711}

Referenced by CellCorrec(), StandardCorrec(), and StandardHitCorrec().

◆ EntaCorrec()

double DedxCorrecSvc::EntaCorrec ( int  layid,
double  enta,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 715 of file DedxCorrecSvc.cxx.

715 {
716 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl;
717// double dedx_enta;
718// if(m_par_flag == 1) {
719// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta +
720// m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3);
721// } else if(m_par_flag == 0) {
722// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) +
723// m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta);
724// }
725// dedx_enta = dedx/dedx_enta;
726// if (dedx_enta>0.0) return dedx_enta;
727// else return dedx;
728
729 if(eangle>-0.25 && eangle<0.25){
730 double stepsize= 0.5/m_venangle.size();
731 int nsize= m_venangle.size();
732 double slope=0;
733 double offset=1;
734 double y1=0,y2=0,x1=0,x2=0;
735
736 if(eangle>=-0.25+0.5*stepsize && eangle<=0.25-0.5*stepsize){
737 int bin = (int)( (eangle-(-0.25+0.5*stepsize))/stepsize);
738 y1=m_venangle[bin];
739 x1=-0.25+0.5*stepsize+bin*stepsize;
740 y2=m_venangle[bin+1];
741 x2=-0.25+1.5*stepsize+bin*stepsize;
742 }
743 else if(eangle<=-0.25+0.5*stepsize){
744 y1=m_venangle[0];
745 x1=-0.25+0.5*stepsize;
746 y2=m_venangle[1];
747 x2=-0.25+1.5*stepsize;
748 }
749 else if( eangle>=0.25-0.5*stepsize){
750 y1=m_venangle[nsize-2];
751 x1=0.25-1.5*stepsize;
752 y2=m_venangle[nsize-1];
753 x2=0.25-0.5*stepsize;
754 }
755 double kcof= (y2-y1)/(x2-x1);
756 double bcof= y1-kcof*x1;
757 double ratio = kcof*eangle+bcof;
758 dedx=dedx/ratio;
759 }
760
761 return dedx;
762}
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ f_larcos()

double DedxCorrecSvc::f_larcos ( double  x,
int  nbin 
)

Definition at line 121 of file DedxCorrecSvc.cxx.

121 {
122 if(nbin!=80) return x; // temporally for only nbin = 80
123 double m_absx(fabs(x));
124 double m_charge(x/m_absx);
125 if(m_absx>0.925 && m_absx<=0.950) return 0.9283*m_charge; // these numbers are from the mean values
126 if(m_absx>0.900 && m_absx<=0.925) return 0.9115*m_charge; // of each bin in cos(theta) distribution
127 if(m_absx>0.875 && m_absx<=0.900) return 0.8877*m_charge;
128 if(m_absx>0.850 && m_absx<=0.875) return 0.8634*m_charge;
129 if(m_absx>0.825 && m_absx<=0.850) return 0.8460*m_charge;
130 if(m_absx>0.800 && m_absx<=0.825) return 0.8089*m_charge;
131}
Double_t x[10]

◆ finalize()

StatusCode DedxCorrecSvc::finalize ( )
virtual

Definition at line 94 of file DedxCorrecSvc.cxx.

94 {
95 // cout<<"DedxCorrecSvc::finalize()"<<endl;
96 MsgStream log(messageService(), name());
97 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq;
98 return StatusCode::SUCCESS;
99}

Referenced by main().

◆ GlobalCorrec()

double DedxCorrecSvc::GlobalCorrec ( double  dedx) const
virtual

Implements IDedxCorrecSvc.

Definition at line 832 of file DedxCorrecSvc.cxx.

832 {
833 if( m_mdcg_flag == 0 ) return dedx;
834 double calib_ex = dedx;
835 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
836 return calib_ex;
837}

Referenced by StandardCorrec().

◆ HadronCorrec()

double DedxCorrecSvc::HadronCorrec ( double  costheta,
double  dedx 
) const

Definition at line 651 of file DedxCorrecSvc.cxx.

651 {
652 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl;
653 //constant 1.00 in the following function is the mean dedx of normalized electrons.
654 dedx=dedx/550.00;
655 return D2I(costheta, I2D(costheta,1.00)/1.00*dedx)*550;
656}
double D2I(const double &cosTheta, const double &D) const
double I2D(const double &cosTheta, const double &I) const

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ handle()

void DedxCorrecSvc::handle ( const Incident &  inc)

Definition at line 101 of file DedxCorrecSvc.cxx.

101 {
102 // cout<<"DedxCorrecSvc::handle"<<endl;
103 MsgStream log( messageService(), name() );
104 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
105
106 if ( inc.type() == "NewRun" ){
107 log << MSG::DEBUG << "New Run" << endreq;
108
109 if( ! m_initConstFlg )
110 {
111 if( m_init == 0 ) { init_param(); }
112 else if( m_init ==1 ) { init_param_svc(); }
113 /* if( ! init_param() ){
114 log << MSG::ERROR
115 << "can not initilize Mdc Calib Constants" << endreq;
116 }*/
117 }
118 }
119}

◆ I2D()

double DedxCorrecSvc::I2D ( const double &  cosTheta,
const double &  I 
) const

Definition at line 2175 of file DedxCorrecSvc.cxx.

2176{
2177 // cout<<" DedxCorrecSvc::I2D"<<endl;
2178 double absCosTheta = fabs(cosTheta);
2179 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2180
2181 // Do the quadratic to compute d of the electron
2182 double a = m_alpha / projection;
2183 double b = 1 - m_gamma / projection*(I/m_ratio);
2184 double c = -I/m_ratio;
2185
2186 if(b==0 && a==0){
2187 cout<<"both a and b coefficiants for hadron correction are 0"<<endl;
2188 return I;
2189 }
2190
2191 double D = a != 0? (-b + sqrt(b*b - 4.0*a*c))/(2.0*a) : -c/b;
2192 if(D<0){
2193 cout<<"D is less 0! will try another solution"<<endl;
2194 D=a != 0? (-b - sqrt(b*b + 4.0*a*c))/(2.0*a) : -c/b;
2195 if(D<0){
2196 cout<<"D is still less 0! just return uncorrectecd value"<<endl;
2197 return I;
2198 }
2199 }
2200
2201 return D;
2202}

Referenced by HadronCorrec().

◆ initialize()

StatusCode DedxCorrecSvc::initialize ( )
virtual

Definition at line 64 of file DedxCorrecSvc.cxx.

64 {
65 // cout<<"DedxCorrecSvc::initialize"<<endl;
66 MsgStream log(messageService(), name());
67 log << MSG::INFO << name() << "DedxCorrecSvc::initialize()" << endreq;
68
69 StatusCode sc = Service::initialize();
70 if( sc.isFailure() ) return sc;
71
72 IIncidentSvc* incsvc;
73 sc = service("IncidentSvc", incsvc);
74 int priority = 100;
75 if( sc.isSuccess() ){
76 //incsvc -> addListener(this, "BeginEvent", priority);
77 incsvc -> addListener(this, "NewRun", priority);
78 }
79
80 StatusCode scgeo = service("MdcGeomSvc", geosvc);
81 if(scgeo.isFailure() ) {
82 log << MSG::ERROR << "GeoSvc failing!!!!!!!" << endreq;
83 return scgeo;
84 }
85
86 StatusCode scmgn = service ("MagneticFieldSvc",m_pIMF);
87 if(scmgn!=StatusCode::SUCCESS) {
88 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
89 }
90
91 return StatusCode::SUCCESS;
92}

Referenced by main().

◆ LayerCorrec()

double DedxCorrecSvc::LayerCorrec ( int  layer,
double  z,
double  costheta,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 892 of file DedxCorrecSvc.cxx.

892 {
893 //cout<<"DedxCorrecSvc::LayerCorrec"<<endl;
894
895 if( m_zdep_flag == 0 && m_ggs_flag == 0 ) return ex;
896
897 double calib_ex = ZdepCorrec( layer, z, ex );
898 if( m_ggs_flag != 0 ) calib_ex = SaturCorrec( layer, costheta, calib_ex );
899 return calib_ex;
900
901}

◆ LayerGainCorrec()

double DedxCorrecSvc::LayerGainCorrec ( int  layid,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 658 of file DedxCorrecSvc.cxx.

658 {
659 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl;
660 if( m_layer_g[layid] > 0 ){
661 double ch_dedx = dedx/m_layer_g[layid];
662 return ch_dedx;
663 }
664 else if( m_layer_g[layid] == 0 ){
665 return dedx;
666 }
667 else return 0;
668}

Referenced by CellCorrec(), and StandardCorrec().

◆ PathL()

double DedxCorrecSvc::PathL ( int  ntpFlag,
const Dedx_Helix hel,
int  layer,
int  cellid,
double  z 
) const
virtual

‍***********----------------—***************---------------—****************‍//

‍***********----------------—***************---------------—****************‍//

assumed the first half circle

‍***********----------------—***************---------------—****************‍//

‍***********----------------—***************---------------—****************‍//

in the Local Helix case, phi must be small

Implements IDedxCorrecSvc.

Definition at line 1489 of file DedxCorrecSvc.cxx.

1489 {
1490
1491 double length = geosvc->Layer(layer)->Length();
1492 int genlay = geosvc->Layer(layer)->Gen();
1493 double East_lay_X = geosvc->GeneralLayer(genlay)->SxEast();
1494 double East_lay_Y = geosvc->GeneralLayer(genlay)->SyEast();
1495 double East_lay_Z = geosvc->GeneralLayer(genlay)->SzEast();
1496 double West_lay_X = geosvc->GeneralLayer(genlay)->SxWest();
1497 double West_lay_Y = geosvc->GeneralLayer(genlay)->SyWest();
1498 double West_lay_Z = geosvc->GeneralLayer(genlay)->SzWest();
1499
1500 HepPoint3D East_origin(East_lay_X, East_lay_Y, East_lay_Z);
1501 HepPoint3D West_origin(West_lay_X, West_lay_Y, West_lay_Z);
1502 Hep3Vector wire = (CLHEP::Hep3Vector)East_origin - (CLHEP::Hep3Vector)West_origin;
1503 HepPoint3D piovt_z =(z*10+length/2 )/length * wire + West_origin;
1504 piovt_z = piovt_z*0.1; // conversion of the units(mm=>cm)
1505
1506
1507 //-------------------------------temp -------------------------------//
1508 HepPoint3D piv(hel.pivot());
1509 //-------------------------------temp -------------------------------//
1510
1511 double dr0 = hel.a()[0];
1512 double phi0 = hel.a()[1];
1513 double kappa = hel.a()[2];
1514 double dz0 = hel.a()[3];
1515 double tanl = hel.a()[4];
1516
1517 // Choose the local field !!
1518 double Bz = 1000*(m_pIMF->getReferField());
1519 double ALPHA_loc = 1000/(2.99792458*Bz);
1520
1521 // Choose the local field !!
1522 //double Bz = 1.0; // will be obtained from magnetic field service
1523 //double ALPHA_loc = 1000/(2.99792458*Bz);
1524 //double ALPHA_loc = 333.564095;
1525 int charge = ( kappa >= 0 )? 1 : -1;
1526 double rho = ALPHA_loc/kappa;
1527 double pt = fabs( 1.0/kappa );
1528 double lambda = atan( tanl );
1529 double theta = M_PI_2- lambda;
1530 double sintheta = sin(M_PI_2-atan(tanl));
1531 // theta = 2.0*M_PI*theta/360.;
1532 double phi = fmod(phi0 + M_PI*4, M_PI*2);
1533 double csf0 = cos(phi);
1534 double snf0 = (1. - csf0) * (1. + csf0);
1535 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1536 if(phi > M_PI) snf0 = - snf0;
1537 //if(ntpFlag>0)
1538 //cout<<"rho = "<<rho<<" hel.dr() + rho = "<<hel.dr() + rho<<endl;
1539
1540 //-------------------------------temp -------------------------------//
1541 double x_c = piv.x() + ( hel.dr() + rho )*csf0;
1542 double y_c = piv.y() + ( hel.dr() + rho )*snf0;
1543 double z_c = piv.z() + hel.dz();
1544 HepPoint3D ccenter(x_c, y_c, 0.0);
1545 double m_c_perp(ccenter.perp());
1546 Hep3Vector m_c_unit((HepPoint3D)ccenter.unit());
1547 //-------------------------------temp -------------------------------//
1548
1549 ////change to boost coordinate system
1550 double x_c_boost = x_c - piovt_z.x();
1551 double y_c_boost = y_c - piovt_z.y();
1552 double z_c_boost = z_c - piovt_z.z();
1553 HepPoint3D ccenter_boost(x_c_boost, y_c_boost, 0.0);
1554 double m_c_perp_boost(ccenter_boost.perp());
1555 //if(ntpFlag>0)
1556 //cout<<" ccenter = "<<ccenter<<" m_c_perp ="<<m_c_perp<<endl;
1557 Hep3Vector m_c_unit_boost((HepPoint3D)ccenter_boost.unit());
1558
1559 //phi region estimation
1560 double phi_io[2];
1561 HepPoint3D IO = (-1)*charge*m_c_unit;
1562 double dphi0 = fmod( IO.phi()+4*M_PI,2*M_PI ) - phi;
1563 double IO_phi = fmod( IO.phi()+4*M_PI,2*M_PI );
1564 //double dphi0_bak = IO_phi - phi;
1565 //if(dphi0 != 0)
1566 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1567 if( dphi0 > M_PI ) dphi0 -= 2*M_PI;
1568 else if( dphi0 < -M_PI ) dphi0 += 2*M_PI;
1569 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1570 //cout<<"charge is = "<<charge<<endl;
1571 phi_io[0] = -(1+charge)*M_PI_2 - charge*dphi0;
1572 //phi_io[0] = -(1+charge)*M_PI_2 + charge*dphi0;
1573 phi_io[1] = phi_io[0]+1.5*M_PI;
1574 //cout<<"phi_io[0] is : "<<phi_io[0]<<" phi_io[1]:"<<phi_io[1]<<endl;
1575 double m_crio[2];
1576 double m_zb, m_zf, Calpha;
1577
1578 // retrieve Mdc geometry information
1579 double rcsiz1= geosvc->Layer(layer)->RCSiz1();
1580 double rcsiz2= geosvc->Layer(layer)->RCSiz2();
1581 double rlay= geosvc->Layer(layer)->Radius();
1582 int ncell= geosvc->Layer(layer)->NCell();
1583 double phioffset= geosvc->Layer(layer)->Offset();
1584 float shift= geosvc->Layer(layer)->Shift();
1585 double slant= geosvc->Layer(layer)->Slant();
1586 //double length = geosvc->Layer(layer)->Length();
1587 int type = geosvc->Layer(layer)->Sup()->Type();
1588
1589 ///***********-------------------***************------------------****************//
1590 //check the value if same //
1591 ///***********-------------------***************------------------****************//
1592 int w0id = geosvc->Layer(layer)->Wirst();
1593 int wid = w0id + cellid;
1594 HepPoint3D backkward = geosvc->Wire(wid)->BWirePos();
1595 HepPoint3D forward = geosvc->Wire(wid)->FWirePos();
1596 double x_lay_backward = geosvc->Wire(layer, cellid)->Backward().x();
1597 double y_lay_backward = geosvc->Wire(layer, cellid)->Backward().y();
1598 double x_lay_forward = geosvc->Wire(layer, cellid)->Forward().x();
1599 double y_lay_forward = geosvc->Wire(layer, cellid)->Forward().y();
1600 double r_lay_backward = sqrt(x_lay_backward*x_lay_backward+y_lay_backward*y_lay_backward);
1601 double r_lay_forward = sqrt(x_lay_forward*x_lay_forward+y_lay_forward*y_lay_forward);
1602 double r_lay_use = ((z*10+length/2)/length)*(r_lay_backward-r_lay_forward) + r_lay_forward;
1603 /*for(int i=0; i<43; i++){
1604 double r_up= geosvc->Layer(i)->RCSiz1();
1605 double r_down= geosvc->Layer(i)->RCSiz2();
1606 cout<<i<<" "<<r_up<<" "<<r_down<<endl;
1607 }*/
1608 // shift = shift*type;
1609 // cout<< "type "<< type <<endl;
1610 r_lay_use = 0.1*r_lay_use;
1611 rcsiz1 = 0.1*rcsiz1;
1612 rcsiz2 = 0.1*rcsiz2;
1613 rlay = 0.1*rlay;
1614 length = 0.1*length;
1615 m_zb = 0.5*length;
1616 m_zf = -0.5*length;
1617 m_crio[0] = rlay - rcsiz1;
1618 m_crio[1] = rlay + rcsiz2;
1619 //conversion of the units(mm=>cm)
1620 int sign = -1; ///assumed the first half circle
1621 int epflag[2];
1622 Hep3Vector iocand[2];
1623 Hep3Vector cell_IO[2];
1624 double dphi, downin;
1625 Hep3Vector zvector;
1626 //if (ntpFlag>0) cout<<"z = "<<z<<" , m_zb = "<<m_zb<<" , m_zf = "<<m_zf<<endl;
1627 if( type ) {
1628 downin = (z*z-m_zb*m_zb)*pow(tan(slant),2);
1629 m_crio[0] = sqrt(m_crio[0]*m_crio[0]+downin);
1630 m_crio[1] = sqrt(m_crio[1]*m_crio[1]+downin);
1631 }
1632
1633 //int stced[2];
1634
1635 for( int i = 0; i < 2; i++ ) {
1636 double cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[i]*m_crio[i] - rho*rho;
1637 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[i] );
1638 if(fabs(cos_alpha)>1&&i==0) return(-1.0);
1639 if(fabs(cos_alpha)>1&&i==1) {
1640 cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[0]*m_crio[0] - rho*rho;
1641 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[0] );
1642 Calpha = 2.0*M_PI-acos( cos_alpha );
1643 } else {
1644 Calpha = acos( cos_alpha );
1645 }
1646 epflag[i] = 0;
1647 iocand[i] = m_c_unit_boost;
1648 iocand[i].rotateZ( charge*sign*Calpha );
1649 iocand[i]*= m_crio[i];
1650 //if (ntpFlag>0) cout<<"iocand corridate is : "<<iocand[i]<<endl;
1651
1652 ///***********-------------------***************------------------****************//
1653 //compare with standard coordinate system results //
1654 ///***********-------------------***************------------------****************//
1655 //-------------------------------temp-------------------------//
1656 iocand[i] = iocand[i]+ piovt_z; //change to standard coordinate system
1657 //iocand[i].y() = iocand[i].y() + piovt_z.y();
1658 //if (ntpFlag>0) cout<<i<<setw(10)<<iocand[i].x()<<setw(10)<<iocand[i].y()<<endl;
1659 //------------------------------temp -------------------------//
1660
1661 double xx = iocand[i].x() - x_c;
1662 double yy = iocand[i].y() - y_c;
1663 //dphi = atan2(yy, xx) - phi0 - M_PI_2*(1+charge);
1664 dphi = atan2(yy, xx) - phi0 - M_PI_2*(1-charge);
1665 dphi = fmod( dphi + 8.0*M_PI, 2*M_PI );
1666
1667 if( dphi < phi_io[0] ) {
1668 dphi += 2*M_PI;
1669 }
1670 else if( phi_io[1] < dphi ) {
1671 dphi -= 2*M_PI;
1672 }
1673 ///in the Local Helix case, phi must be small
1674
1675 //Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl);
1676 Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl-piovt_z.z());
1677 //if (ntpFlag>0) cout<<"z_c-rho*dphi*tanl : "<<z_c-rho*dphi*tanl<<endl;
1678 cell_IO[i] = iocand[i];
1679 cell_IO[i] += zvector;
1680 //---------------------------------temp ---------------------------------//
1681 //cell_IO[i] = hel.x(dphi);//compare with above results
1682 //---------------------------------temp ---------------------------------//
1683
1684 double xcio, ycio, phip;
1685 ///// z region check active volume is between zb+2. and zf-2. [cm]
1686 /*
1687 float delrin = 2.0;
1688 if( m_zf-delrin > zvector.z() ){
1689 phip = z_c - m_zb + delrin;
1690 phip = phip/( rho*tanl );
1691 phip = phip + phi0;
1692 xcio = x_c - rho*cos( phip );
1693 ycio = y_c - rho*sin( phip );
1694 cell_IO[i].setX( xcio );
1695 cell_IO[i].setY( ycio );
1696 cell_IO[i].setZ( m_zb+delrin );
1697 epflag[i] = 1;
1698 }
1699 else if( m_zb+delrin < zvector.z() ){
1700 phip = z_c - m_zf-delrin;
1701 phip = phip/( rho*tanl );
1702 phip = phip + phi0;
1703 xcio = x_c - rho*cos( phip );
1704 ycio = y_c - rho*sin( phip );
1705 cell_IO[i].setX( xcio );
1706 cell_IO[i].setY( ycio );
1707 cell_IO[i].setZ( m_zf-delrin );
1708 epflag[i] = 1;
1709 }
1710 else{
1711 */
1712 // cell_IO[i] = iocand;
1713 // cell_IO[i] += zvector;
1714 // }
1715 //dphi = dphi -M_PI;
1716 cell_IO[i] = hel.x(dphi);
1717 //if (ntpFlag>0) { cout<<"cell_IO corridate : "<<cell_IO[i]<<endl;}
1718 // if(i==0) cout<<"zhit value is = "<<z<<endl;
1719 }
1720
1721 // path length estimation
1722 // phi calculation
1723 Hep3Vector cl = cell_IO[1] - cell_IO[0];
1724
1725 //float ch_tphi, ch_tdphi;
1726 double ch_theta;
1727 double ch_dphi;
1728 double ch_ltrk = 0;
1729 double ch_ltrk_rp = 0;
1730 ch_dphi = cl.perp()*0.5/(ALPHA_loc*pt);
1731 ch_dphi = 2.0 * asin( ch_dphi );
1732 ch_ltrk_rp = ch_dphi*ALPHA_loc*pt;
1733 double rpi_path = sqrt(cl.x()*cl.x()+cl.y()*cl.y());
1734 ch_ltrk = sqrt( ch_ltrk_rp*ch_ltrk_rp + cl.z()*cl.z() );
1735 double path = ch_ltrk_rp/ sintheta;
1736 ch_theta = cl.theta();
1737 /* if (ntpFlag>0)
1738 {
1739 // if((ch_ltrk_rp-rpi_path)>0.001 || (ch_ltrk-path)>0.001)
1740 cout<<"ch_ltrk_rp : "<<ch_ltrk_rp<<" rpi_path: "<<rpi_path<<endl;
1741 cout<<"ch_ltrk : "<<ch_ltrk<<" path:"<<path<<endl;
1742 }
1743 */
1744 /*
1745 if( ch_theta < theta*0.85 || 1.15*theta < ch_theta)
1746 ch_ltrk *= -1; /// miss calculation
1747
1748 if( epflag[0] == 1 || epflag [1] == 1 )
1749 ch_ltrk *= -1; /// invalid region for dE/dx or outside of Mdc
1750 */
1751 // judge how many cells are traversed and deal with different situations
1752 double phibin, phi_in, phi_out, inlow, inup, outlow, outup, gap, phi1, phi2, phi_mid, phi_midin, phi_midout;
1753 int cid_in, cid_out;
1754 double inlow_z, inup_z, outlow_z, outup_z, gap_z, phi1_z, phi2_z, phi_mid_z, phi_midin_z, phi_midout_z;
1755 //cache sampl in each cell for this layer
1756
1757 std::vector<double> sampl;
1758 sampl.resize(ncell);
1759 for(int k=0; k<ncell; k++) {
1760 sampl[k] = -1.0;
1761 }
1762
1763 cid_in = cid_out = -1;
1764 phi_in = cell_IO[0].phi();
1765 phi_out = cell_IO[1].phi();
1766 //phi_in = iocand[0].phi();
1767 //phi_out = iocand[1].phi();
1768 phi_in = fmod( phi_in+4*M_PI,2*M_PI );
1769 phi_out = fmod( phi_out+4*M_PI,2*M_PI );
1770 phibin = 2.0*M_PI/ncell;
1771 // gap = fabs(phi_out-phi_in);
1772
1773 //determine the in/out cell id
1774 Hep3Vector cell_mid=0.5*(cell_IO[0]+cell_IO[1]);
1775 //Hep3Vector cell_mid=0.5*(iocand[0]+iocand[0]);
1776 //cout<<cell_mid<<endl;
1777 //double stereophi = shift*phibin*(0.5-cell_mid.z()/length);
1778 //phioffset = phioffset+stereophi;
1779 double stphi[2], phioff[2];
1780 stphi[0] = shift*phibin*(0.5-cell_IO[0].z()/length);
1781 stphi[1] = shift*phibin*(0.5-cell_IO[1].z()/length);
1782 //stphi[0] = shift*phibin*(0.5-z/length);
1783 //stphi[1] = shift*phibin*(0.5-z/length);
1784 phioff[0] = phioffset+stphi[0];
1785 phioff[1] = phioffset+stphi[1];
1786
1787 for(int i=0; i<ncell; i++) {
1788
1789 double x_lay_backward_cell = geosvc->Wire(layer, i)->Backward().x()*0.1;
1790 double y_lay_backward_cell = geosvc->Wire(layer, i)->Backward().y()*0.1;
1791 double x_lay_forward_cell = geosvc->Wire(layer, i)->Forward().x()*0.1;
1792 double y_lay_forward_cell = geosvc->Wire(layer, i)->Forward().y()*0.1;
1793 //if(ntpFlag>0) cout<<layer<<setw(10)<<i<<setw(10)<<x_lay_forward<<setw(10)<<y_lay_forward<<setw(10)<<x_lay_backward<<setw(10)<<y_lay_backward<<setw(10)<<r_lay_forward<<setw(10)<<endl;
1794 //Hep3Vector lay_backward, lay_forward;
1795 Hep3Vector lay_backward(x_lay_backward_cell, y_lay_backward_cell, 0);
1796 Hep3Vector lay_forward(x_lay_forward_cell, y_lay_forward_cell, 0);
1797 Hep3Vector Cell_z[2];
1798 Cell_z[0] = ((cell_IO[0].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1799 Cell_z[1] = ((cell_IO[1].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1800 double z_phi[2];
1801 z_phi[0] = Cell_z[0].phi();
1802 z_phi[1] = Cell_z[1].phi();
1803 z_phi[0] = fmod( z_phi[0]+4*M_PI,2*M_PI );
1804 z_phi[1] = fmod( z_phi[1]+4*M_PI,2*M_PI );
1805 //double backward_phi = lay_backward.phi();
1806 //double forward_phi = lay_forward.phi();
1807 //backward_phi = fmod( backward_phi+4*M_PI,2*M_PI );
1808 //forward_phi = fmod( forward_phi+4*M_PI,2*M_PI );
1809 //if(ntpFlag>0) cout<<"backward_phi cal : "<<atan2(y_lay_backward,x_lay_backward)<<" forward_phi : "<<atan2(y_lay_forward,x_lay_forward)<<endl;
1810 //if(ntpFlag>0) cout<<layer<<" backward_phi : "<<backward_phi<<" forward_phi : "<<forward_phi<<endl;
1811
1812 //double z_phi[2];
1813 //z_phi[0] = ((cell_IO[0].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1814 //z_phi[1] = ((cell_IO[1].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1815 //if(ntpFlag>0) cout<<"phi z : "<<z_phi[0]<<" "<<z_phi[1]<<endl;
1816 inlow_z = z_phi[0] - phibin*0.5;
1817 inup_z = z_phi[0] + phibin*0.5;
1818 outlow_z = z_phi[1] - phibin*0.5;
1819 outup_z = z_phi[1] + phibin*0.5;
1820 inlow_z = fmod( inlow_z+4*M_PI,2*M_PI );
1821 inup_z = fmod( inup_z+4*M_PI,2*M_PI );
1822 outlow_z = fmod( outlow_z+4*M_PI,2*M_PI );
1823 outup_z = fmod( outup_z+4*M_PI,2*M_PI );
1824
1825
1826 // for stereo layer
1827 inlow = phioff[0]+phibin*i-phibin*0.5;
1828 inup = phioff[0]+phibin*i+phibin*0.5;
1829 outlow = phioff[1]+phibin*i-phibin*0.5;
1830 outup = phioff[1]+phibin*i+phibin*0.5;
1831 inlow = fmod( inlow+4*M_PI,2*M_PI );
1832 inup = fmod( inup+4*M_PI,2*M_PI );
1833 outlow = fmod( outlow+4*M_PI,2*M_PI );
1834 outup = fmod( outup+4*M_PI,2*M_PI );
1835
1836#ifdef YDEBUG
1837 if(ntpFlag > 0) cout << "shift " << shift
1838 <<" phi_in " << phi_in << " phi_out " << phi_out
1839 << " inup "<< inup << " inlow " << inlow
1840 << " outup "<< outup << " outlow " << outlow << endl;
1841#endif
1842
1843 /*if(phi_in>=inlow&&phi_in<inup) cid_in = i;
1844 if(phi_out>=outlow&&phi_out<outup) cid_out = i;
1845 if ( inlow>inup) {
1846 if((phi_in>=inlow&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup)) cid_in = i;
1847 }
1848 if ( outlow>outup) {
1849 if((phi_out>=outlow&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup)) cid_out = i;
1850 }*/
1851
1852 if(phi_in>=inlow_z&&phi_in<inup_z) cid_in = i;
1853 if(phi_out>=outlow_z&&phi_out<outup_z) cid_out = i;
1854 if ( inlow_z>inup_z) {
1855 if((phi_in>=inlow_z&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup_z)) cid_in = i;
1856 }
1857 if ( outlow_z>outup_z) {
1858 if((phi_out>=outlow_z&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup_z)) cid_out = i;
1859 }
1860 }
1861
1862 phi_midin = phi_midout = phi1 = phi2 = -999.0;
1863 gap = -999.0;
1864 //only one cell traversed
1865 //if(ntpFlag > 0) cout<<"cid_in = "<<cid_in<<" cid_out = "<<cid_out<<endl;
1866 if(cid_in == -1 || cid_out == -1) return -1;
1867
1868 if( cid_in == cid_out) {
1869 sampl[cid_in]= ch_ltrk;
1870 //return ch_ltrk;
1871 return sampl[cellid];
1872 } else if(cid_in < cid_out) {
1873 //in cell id less than out cell id
1874 //deal with the special case crossing x axis
1875 if( cid_out-cid_in>ncell/2 ) {
1876
1877 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1878 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
1879 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
1880 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
1881 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
1882 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1883 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1884 Hep3Vector Cell_z[2];
1885 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1886 double phi_cin_z = Cell_z[0].phi();
1887 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
1888 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
1889 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
1890 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
1891 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
1892 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1893 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1894 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1895 double phi_cout_z = Cell_z[1].phi();
1896 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
1897
1898 phi_midin_z = phi_cin_z-phibin*0.5;
1899 phi_midout_z = phi_cout_z+phibin*0.5;
1900 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
1901 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
1902 phi1_z = phi_midout_z-phi_out;
1903 phi2_z = phi_in-phi_midin_z;
1904 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
1905 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
1906 gap_z = phi1_z+phi2_z+(ncell-1-cid_out+cid_in)*phibin;
1907 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
1908 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
1909 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
1910 for( int jj = cid_out+1; jj<ncell; jj++) {
1911 sampl[jj]=phibin/gap_z*ch_ltrk;
1912 }
1913 for( int jj = 0; jj<cid_in; jj++) {
1914 sampl[jj]=phibin/gap_z*ch_ltrk;
1915 }
1916
1917 /*
1918 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1919 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
1920 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
1921 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1922 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1923 phi1 = phi_midout-phi_out;
1924 phi2 = phi_in-phi_midin;
1925 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1926 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1927 gap = phi1+phi2+(ncell-1-cid_out+cid_in)*phibin;
1928 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1929 sampl[cid_in]=phi2/gap*ch_ltrk;
1930 sampl[cid_out]=phi1/gap*ch_ltrk;
1931 for( int jj = cid_out+1; jj<ncell; jj++) {
1932 sampl[jj]=phibin/gap*ch_ltrk;
1933 }
1934 for( int jj = 0; jj<cid_in; jj++) {
1935 sampl[jj]=phibin/gap*ch_ltrk;
1936 }*/
1937 /*
1938 cout<<"LLLLLLLLLLLLL"<<endl;
1939 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
1940 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
1941 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
1942 << phi_out << " phi_midout " << phi_midout <<endl;
1943 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
1944 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
1945 */
1946 } else {
1947 //normal case
1948 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
1949 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
1950 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
1951 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
1952 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1953 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1954 Hep3Vector Cell_z[2];
1955 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1956 double phi_cin_z = Cell_z[0].phi();
1957 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
1958 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
1959 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
1960 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
1961 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
1962 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1963 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1964 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1965 double phi_cout_z = Cell_z[1].phi();
1966 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
1967
1968 phi_midin_z = phi_cin_z+phibin*0.5;
1969 phi_midout_z = phi_cout_z-phibin*0.5;
1970 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
1971 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
1972 phi1_z = phi_midin_z-phi_in;
1973 phi2_z = phi_out-phi_midout_z;
1974 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
1975 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
1976 gap_z = phi1_z+phi2_z+(cid_out-cid_in-1)*phibin;
1977 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
1978 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
1979 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
1980 for( int jj = cid_in+1; jj<cid_out; jj++) {
1981 sampl[jj]=phibin/gap_z*ch_ltrk;
1982 }
1983 //normal case
1984 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
1985 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
1986 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1987 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1988 phi1 = phi_midin-phi_in;
1989 phi2 = phi_out-phi_midout;
1990 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1991 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1992 gap = phi1+phi2+(cid_out-cid_in-1)*phibin;
1993 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1994 sampl[cid_in]=phi1/gap*ch_ltrk;
1995 sampl[cid_out]=phi2/gap*ch_ltrk;
1996 for( int jj = cid_in+1; jj<cid_out; jj++) {
1997 sampl[jj]=phibin/gap*ch_ltrk;
1998 }*/
1999 }
2000
2001 } else if(cid_in > cid_out) {
2002 //in cell id greater than out cell id
2003 //deal with the special case crossing x axis
2004 if( cid_in-cid_out>ncell/2 ) {
2005 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
2006 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
2007 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
2008 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
2009 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2010 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2011 Hep3Vector Cell_z[2];
2012 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2013 double phi_cin_z = Cell_z[0].phi();
2014 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
2015 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
2016 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
2017 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
2018 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
2019 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2020 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2021 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2022 double phi_cout_z = Cell_z[1].phi();
2023 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
2024
2025 phi_midin_z = phi_cin_z+phibin*0.5;
2026 phi_midout_z = phi_cout_z-phibin*0.5;
2027 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
2028 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
2029 phi1_z = phi_midin_z-phi_in;
2030 phi2_z = phi_out-phi_midout_z;
2031 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
2032 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
2033 gap_z = phi1_z+phi2_z+(ncell-1-cid_in+cid_out)*phibin;
2034 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
2035 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
2036 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
2037 for( int jj = cid_in+1; jj<ncell; jj++) {
2038 sampl[jj]=phibin/gap_z*ch_ltrk;
2039 }
2040 for( int jj = 0; jj<cid_out; jj++) {
2041 sampl[jj]=phibin/gap_z*ch_ltrk;
2042 }
2043
2044 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
2045 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
2046 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
2047 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2048 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2049 phi1 = phi_midin-phi_in;
2050 phi2 = phi_out-phi_midout;
2051 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2052 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2053 gap = phi1+phi2+(ncell-1-cid_in+cid_out)*phibin;
2054 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2055 sampl[cid_out]=phi2/gap*ch_ltrk;
2056 sampl[cid_in]=phi1/gap*ch_ltrk;
2057 for( int jj = cid_in+1; jj<ncell; jj++) {
2058 sampl[jj]=phibin/gap*ch_ltrk;
2059 }
2060 for( int jj = 0; jj<cid_out; jj++) {
2061 sampl[jj]=phibin/gap*ch_ltrk;
2062 }*/
2063 /*
2064 cout<<"LLLLLLLLLLLLL"<<endl;
2065 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2066 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2067 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2068 << phi_out << " phi_midout " << phi_midout <<endl;
2069 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2070 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
2071 */
2072 } else {
2073 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1;
2074 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1;
2075 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1;
2076 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1;
2077 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2078 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2079 Hep3Vector Cell_z[2];
2080 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2081 double phi_cin_z = Cell_z[0].phi();
2082 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI );
2083 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1;
2084 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1;
2085 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1;
2086 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1;
2087 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2088 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2089 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2090 double phi_cout_z = Cell_z[1].phi();
2091 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI );
2092
2093 phi_midin_z = phi_cin_z-phibin*0.5;
2094 phi_midout_z = phi_cout_z+phibin*0.5;
2095 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI );
2096 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI );
2097 phi1_z = phi_midout_z-phi_out;
2098 phi2_z = phi_in-phi_midin_z;
2099 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI);
2100 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI);
2101 gap_z = phi1_z+phi2_z+(cid_in-cid_out-1)*phibin;
2102 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI);
2103 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
2104 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
2105 for( int jj = cid_out+1; jj<cid_in; jj++) {
2106 sampl[jj]=phibin/gap_z*ch_ltrk;
2107 }
2108
2109 //normal case
2110 /*phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
2111 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
2112 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2113 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2114 phi1 = phi_midout-phi_out;
2115 phi2 = phi_in-phi_midin;
2116 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2117 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2118 gap = phi1+phi2+(cid_in-cid_out-1)*phibin;
2119 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2120 sampl[cid_in]=phi2/gap*ch_ltrk;
2121 sampl[cid_out]=phi1/gap*ch_ltrk;
2122 for( int jj = cid_out+1; jj<cid_in; jj++) {
2123 sampl[jj]=phibin/gap*ch_ltrk;
2124 }*/
2125 }
2126 }
2127
2128#ifdef YDEBUG
2129 if(sampl[cellid]<0.0) {
2130 if(cid_in!=cid_out) cout<<"?????????"<<endl;
2131 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2132 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2133
2134 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2135 << phi_out << " phi_midout " << phi_midout <<endl;
2136 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2137 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
2138
2139
2140 for(int l=0; l<ncell; l++) {
2141 if(sampl[l]>=0.0)
2142 cout<<"picked cellid " << l << " sampling length "<< sampl[l]<< endl;
2143 }
2144 }
2145#endif
2146 return sampl[cellid];
2147}
double tan(const BesAngle a)
Definition: BesAngle.h:216
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
double length
Double_t phi2
Double_t phi1
#define M_PI
Definition: TConstant.h:4
const HepVector & a(void) const
returns helix parameters.
Definition: Dedx_Helix.h:238
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Dedx_Helix.cxx:209
double dz(void) const
Definition: Dedx_Helix.h:220
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Dedx_Helix.h:184
double dr(void) const
returns an element of parameters.
Definition: Dedx_Helix.h:202
virtual double getReferField()=0
virtual const MdcGeoLayer *const Layer(unsigned id)=0
virtual const MdcGeoWire *const Wire(unsigned id)=0
virtual const MdcGeoGeneral *const GeneralLayer(unsigned id)=0
double SzWest(void) const
double SxEast(void) const
double SyWest(void) const
double SxWest(void) const
double SzEast(void) const
double SyEast(void) const
double Radius(void) const
Definition: MdcGeoLayer.h:160
double Slant(void) const
Definition: MdcGeoLayer.h:158
int Gen(void) const
Definition: MdcGeoLayer.h:175
double RCSiz2(void) const
Definition: MdcGeoLayer.h:163
double Length(void) const
Definition: MdcGeoLayer.h:161
MdcGeoSuper * Sup(void) const
Definition: MdcGeoLayer.h:174
double Shift(void) const
Definition: MdcGeoLayer.h:167
double RCSiz1(void) const
Definition: MdcGeoLayer.h:162
int NCell(void) const
Definition: MdcGeoLayer.h:165
int Wirst(void) const
Definition: MdcGeoLayer.h:157
double Offset(void) const
Definition: MdcGeoLayer.h:166
int Type(void) const
Definition: MdcGeoSuper.h:35
HepPoint3D FWirePos(void) const
Definition: MdcGeoWire.h:131
HepPoint3D BWirePos(void) const
Definition: MdcGeoWire.h:130
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128

◆ queryInterface()

StatusCode DedxCorrecSvc::queryInterface ( const InterfaceID &  riid,
void **  ppvUnknown 
)
virtual

Definition at line 54 of file DedxCorrecSvc.cxx.

54 {
55 //cout<<"DedxCorrecSvc::queryInterface"<<endl;
56 if( IID_IDedxCorrecSvc.versionMatch(riid) ){
57 *ppvInterface = static_cast<IDedxCorrecSvc*> (this);
58 } else{
59 return Service::queryInterface(riid, ppvInterface);
60 }
61 return StatusCode::SUCCESS;
62}

◆ RungCorrec()

double DedxCorrecSvc::RungCorrec ( int  runNO,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 672 of file DedxCorrecSvc.cxx.

672 {
673 // cout<<"DedxCorrecSvc::RungCorrec"<<endl;
674 double dedx_rung;
675 int run_ture =0;
676 //cout<<"N_run : "<<N_run<<" runNO : "<<runNO<<endl;
677
678 for(int j=0;j<10000;j++) {
679 if(runNO == m_rung[2][j]) {
680 dedx_rung = dedx/m_rung[0][j];
681 run_ture =1;
682 return dedx_rung;
683 }
684 }
685 if(run_ture ==0)
686 {
687 cout<<"Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to "<<runNO<<endl;
688 exit(1);
689 }
690}

Referenced by StandardCorrec(), StandardHitCorrec(), and StandardTrackCorrec().

◆ SaturCorrec()

double DedxCorrecSvc::SaturCorrec ( int  layid,
double  costheta,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 545 of file DedxCorrecSvc.cxx.

545 {
546 //cout<<"DedxCorrecSvc::SaturCorrec"<<endl;
547 double dedx_ggs;
548 //cout<<"costheta vaule is = "<<costheta<<endl;
549 costheta = fabs(costheta);
550 if(m_par_flag == 1) {
551 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*costheta +
552 m_ggs[2][layer]*pow(costheta,2) + m_ggs[3][layer]*pow(costheta,3);
553 } else if(m_par_flag == 0) {
554 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*T1(costheta) +
555 m_ggs[2][layer]*T2(costheta) + m_ggs[3][layer]*T3(costheta);
556 }
557 //cout<<"dedx_ggs is = "<<dedx_ggs<<endl;
558 dedx_ggs = dedx/dedx_ggs;
559 if(dedx_ggs>0.0) return dedx_ggs;
560 else return dedx;
561}

Referenced by CellCorrec(), LayerCorrec(), StandardCorrec(), and StandardHitCorrec().

◆ set_flag()

void DedxCorrecSvc::set_flag ( int  calib_F)
virtual

Implements IDedxCorrecSvc.

Definition at line 1131 of file DedxCorrecSvc.cxx.

1131 {
1132 // cout<<"DedxCorrecSvc::set_flag"<<endl;
1133 cout<<"calib_F is = "<<calib_F<<endl;
1134 if(calib_F<0){ m_enta_flag = 0; calib_F = abs(calib_F); }
1135 else m_enta_flag = 1;
1136 m_rung_flag = ( calib_F & 1 )? 1 : 0;
1137 m_wireg_flag = ( calib_F & 2 )? 1 : 0;
1138 m_dcasin_flag = ( calib_F & 4 )? 1 : 0;
1139 m_ggs_flag = ( calib_F & 8 )? 1 : 0;
1140 m_ddg_flag = 0;
1141 //m_ddg_flag = ( calib_F & 8 )? 1 : 0;
1142 m_t0_flag = ( calib_F & 16 )? 1 : 0;
1143 m_sat_flag = ( calib_F & 32 )? 1 : 0;
1144 m_layerg_flag = ( calib_F & 64 )? 1 : 0;
1145 //m_dcasin_flag = ( calib_F & 128 )? 1 : 0;
1146 m_dip_flag = ( calib_F & 256 )? 1 : 0;
1147 m_mdcg_flag = ( calib_F & 512 )? 1 : 0;
1148}
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212

◆ StandardCorrec()

double DedxCorrecSvc::StandardCorrec ( int  runFlag,
int  ntpFlag,
int  runNO,
double  pathl,
int  wid,
int  layid,
double  adc,
double  dd,
double  eangle,
double  z,
double  costheta 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 919 of file DedxCorrecSvc.cxx.

919 {
920 // cout<<"DedxCorrecSvc::StandardCorrec"<<endl;
921 //int layid = MdcID::layer(mdcid);
922 //int localwid = MdcID::wire(mdcid);
923 //int w0id = geosvc->Layer(layid)->Wirst();
924 //int wid = w0id + localwid;
925 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
926 ////pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
927 double ex = adc;
928 if( runNO<0 ) return ex;
929 ex = ex*1.5/pathl;
930 //double ex = adc/pathl;
931 //if( runNO<0 ) return ex;
932 if( ntpFlag ==0) return ex;
933 //double ex = adc/pathl;
934 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_dcasin_flag==0 && m_ddg_flag == 0
935 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
936 m_ggs_flag == 0) return ex;
937
938 if(m_rung_flag) {
939 ex = RungCorrec( runNO, ex );
940 }
941
942 if( m_wireg_flag) {
943 ex = WireGainCorrec(wid, ex);
944 }
945
946 if( m_dcasin_flag) {
947 if(runFlag<3) {
948 ex = DriftDistCorrec( layid, dd, ex );
949 }
950 else{ ex = DocaSinCorrec(layid, dd, eangle, ex);}
951 }
952
953 if(m_enta_flag && runFlag>=3){
954 ex = EntaCorrec(layid, eangle, ex);
955 }
956
957 // ddg is not use anymore, it's replaced by DocaSin
958 if(m_ddg_flag) {
959 ex = DriftDistCorrec( layid, dd, ex );
960 }
961
962 if(m_ggs_flag) {
963 if(runFlag <3 ){
964 ex = SaturCorrec( layid, costheta, ex );
965 }
966 else{ ex = CosthetaCorrec( costheta, ex );}
967 // Staur is OLD, Costheta is NEW.
968 }
969
970 if( m_sat_flag) {
971 ex = HadronCorrec( costheta, ex );
972 }
973
974 if( m_zdep_flag) {
975 ex = ZdepCorrec( layid, z, ex );
976 }
977
978 if( m_layerg_flag) {
979 ex = LayerGainCorrec( layid, ex );
980 }
981
982 if( m_dip_flag){
983 ex = DipAngCorrec(layid, costheta, ex);
984 }
985
986 if( m_mdcg_flag) {
987 ex = GlobalCorrec( ex );
988 }
989 return ex;
990}
double CosthetaCorrec(double costheta, double ex) const
double RungCorrec(int runNO, double ex) const
double EntaCorrec(int layid, double enta, double ex) const
double WireGainCorrec(int wireid, double ex) const
double HadronCorrec(double costheta, double dedx) const
double DocaSinCorrec(int layid, double doca, double enta, double ex) const
double DipAngCorrec(int layid, double costheta, double ex) const
double GlobalCorrec(double dedx) const

◆ StandardHitCorrec()

double DedxCorrecSvc::StandardHitCorrec ( int  calib_rec_Flag,
int  runFlag,
int  ntpFlag,
int  runNO,
double  pathl,
int  wid,
int  layid,
double  adc,
double  dd,
double  eangle,
double  costheta 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 993 of file DedxCorrecSvc.cxx.

993 {
994 if(m_debug) cout<<"DedxCorrecSvc::StandardHitCorrec"<<endl;
995 //int layid = MdcID::layer(mdcid);
996 //int localwid = MdcID::wire(mdcid);
997 //int w0id = geosvc->Layer(layid)->Wirst();
998 //int wid = w0id + localwid;
999 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
1000 //cout<<"DedxCorrecSvc::StandardHitCorrec pathl= "<<pathl<<endl;
1001 //pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
1002 double ex = adc;
1003 if( runNO<0 ) return ex;
1004 ex = ex*1.5/pathl;
1005 //if( runNO<0 ) return ex;
1006 //double ex = adc/pathl;
1007 if( ntpFlag ==0) return ex;
1008 //if(ntpFlag>0) cout<<"dE/dx value after path correction: "<<ex<<endl;
1009 //double ex = adc/pathl;
1010 //cout<<m_rung_flag << " "<<m_wireg_flag<<" "<<m_ddg_flag<<" "<<m_ggs_flag<<endl;
1011 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
1012 && m_layerg_flag == 0 && m_zdep_flag == 0 && m_dcasin_flag==0
1013 && m_ggs_flag == 0 && m_enta_flag==0) return ex;
1014
1015 if(m_rung_flag) {
1016 ex = RungCorrec( runNO, ex );
1017 }
1018 //if(ntpFlag>0) cout<<"dE/dx value after run by run correction: "<<ex<<endl;
1019
1020 if( m_wireg_flag) {
1021 ex = WireGainCorrec(wid, ex);
1022 }
1023 //if(ntpFlag>0) cout<<"dE/dx value after wire gain correction: "<<ex<<endl;
1024
1025 if( m_dcasin_flag) {
1026 if(runFlag<3){
1027 ex = DriftDistCorrec( layid, dd, ex );
1028 }
1029 else{
1030 //cout<<layid<<" "<<dd<<" "<<eangle<<" "<<ex<<endl;
1031 ex = DocaSinCorrec(layid, dd, eangle, ex);
1032 }
1033 }
1034
1035 // 1D entrance angle correction
1036 if(m_enta_flag && runFlag>=3){
1037 ex = EntaCorrec(layid, eangle, ex);
1038 }
1039
1040 // ddg is not used anymore
1041 if( m_ddg_flag) {
1042 ex = DriftDistCorrec( layid, dd, ex );
1043 }
1044 //if(ntpFlag>0) cout<<"dE/dx value after ddg correction: "<<ex<<endl;
1045
1046 if(m_ggs_flag) {
1047 if(runFlag <3 ){
1048 ex = SaturCorrec( layid, costheta, ex );
1049 }
1050 else{ ex = ex;} // do not do the cos(theta) correction at hit level
1051 }
1052 //if(ntpFlag>0) cout<<"dE/dx value after costheta correction: "<<ex<<endl;
1053 return ex;
1054}

◆ StandardTrackCorrec()

double DedxCorrecSvc::StandardTrackCorrec ( int  calib_rec_Flag,
int  typFlag,
int  ntpFlag,
int  runNO,
double  ex,
double  costheta,
double  t0 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 1058 of file DedxCorrecSvc.cxx.

1058 {
1059 if(m_debug) cout<<"DedxCorrecSvc::StandardTrackCorrec"<<endl;
1060 if( runNO<0 ) return ex;
1061 if( ntpFlag ==0) return ex;
1062 if( runFlag <3) return ex;
1063 if( calib_rec_Flag ==1){
1064 ex = CosthetaCorrec( costheta, ex );
1065 if(m_t0_flag) ex= T0Correc(t0, ex);
1066 return ex;
1067 }
1068
1069 //if(ntpFlag>0) cout<<"trcak value before costheta correction: "<<ex<<endl;
1070 if( m_ggs_flag) {
1071 ex = CosthetaCorrec( costheta, ex );
1072 }
1073 //if(ntpFlag>0) cout<<"trcak value after costheta correction: "<<ex<<endl;
1074 if(m_t0_flag){
1075 if(runFlag==3) {ex= T0Correc(t0, ex);}
1076 else if(runFlag>3) {ex=ex;}
1077 }
1078 //if(ntpFlag>0) cout<<"trcak value after all correction: "<<ex<<endl;
1079 if( m_sat_flag) {
1080 ex = HadronCorrec( costheta, ex );
1081 }
1082
1083 if(m_rung_flag && calib_rec_Flag ==2){
1084 double rungain_temp =RungCorrec( runNO, ex )/ex;
1085 ex = ex/rungain_temp;
1086 }
1087
1088 //if(ntpFlag>0) cout<<"trcak value no run gain correction: "<<ex<<endl;
1089 return ex;
1090
1091}
double T0Correc(double t0, double ex) const

◆ T0Correc()

double DedxCorrecSvc::T0Correc ( double  t0,
double  ex 
) const

Definition at line 619 of file DedxCorrecSvc.cxx.

619 {
620 // cout<<"DedxCorrecSvc::T0Correc"<<endl;
621 double dedx_t0;
622 const int nbin=t0_k.size()+1;
623 if(nbin!=35)
624 cout<<"there should be 35 bins for t0 correction, check it!"<<endl;
625
626 int n=0;
627 if(t0<575)
628 n=0;
629 else if(t0<610)
630 n=1;
631 else if(t0>=610 && t0<1190){
632 n=(int)( (t0-610.0)/20.0 ) + 2;
633 }
634 else if(t0>=1190 && t0<1225)
635 n=31;
636 else if(t0>=1225 && t0<1275)
637 n=32;
638 else if(t0>=1275)
639 n=33;
640
641 dedx_t0=t0_k[n]*t0+t0_b[n];
642
643 if(dedx_t0>0){
644 dedx_t0 = dedx/dedx_t0;
645 return dedx_t0;
646 }
647 else return dedx;
648}

Referenced by StandardTrackCorrec().

◆ TrkCorrec()

double DedxCorrecSvc::TrkCorrec ( double  costheta,
double  dedx 
) const
virtual

dEdx calib. for each track

Implements IDedxCorrecSvc.

Definition at line 904 of file DedxCorrecSvc.cxx.

904 {
905 // cout<<"DedxCorrecSvc::TrkCorrec"<<endl;
906 if( m_mdcg_flag == 0 ) return dedx;
907 ///dEdx calib. for each track
908 double calib_ex = dedx;
909
910 double run_const = m_dedx_gain;
911 if( run_const > 0 && m_mdcg_flag != 0 ) calib_ex /= run_const;
912
913 return calib_ex;
914
915}

◆ WireGainCorrec()

double DedxCorrecSvc::WireGainCorrec ( int  wireid,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 533 of file DedxCorrecSvc.cxx.

533 {
534 if( m_wire_g[wireid] > 0 ){
535 double ch_dedx = (ex/m_wire_g[wireid])*m_valid[wireid];
536 return ch_dedx;
537 }
538 else if( m_wire_g[wireid] == 0 ){
539 return ex;
540 }
541 else return 0;
542}

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ ZdepCorrec()

double DedxCorrecSvc::ZdepCorrec ( int  layid,
double  zhit,
double  ex 
) const
virtual

Implements IDedxCorrecSvc.

Definition at line 766 of file DedxCorrecSvc.cxx.

766 {
767 // cout<<"DedxCorrecSvc::ZdepCorrec"<<endl;
768 double dedx_zdep;
769 if(m_par_flag == 1) {
770 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*z +
771 m_zdep[2][layer]*z*z + m_zdep[3][layer]*pow(z,3);
772 } else if(m_par_flag == 0) {
773 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*T1(z) +
774 m_zdep[2][layer]*T2(z) + m_zdep[3][layer]*T3(z);
775 }
776 dedx_zdep = dedx/dedx_zdep;
777 if (dedx_zdep>0.0) return dedx_zdep;
778 else return dedx;
779
780 //return dedx;
781}

Referenced by CellCorrec(), LayerCorrec(), and StandardCorrec().


The documentation for this class was generated from the following files: