BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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)
 
 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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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
 
virtual double RungCorrec (int runNO, int evtNO, 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, int evtNO, 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, int evtNO, 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, int evtNO, 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 ()
 
static const InterfaceID & interfaceID ()
 

Detailed Description

Definition at line 19 of file InstallArea/include/DedxCorrecSvc/DedxCorrecSvc/DedxCorrecSvc.h.

Constructor & Destructor Documentation

◆ DedxCorrecSvc() [1/2]

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

◆ ~DedxCorrecSvc() [1/2]

DedxCorrecSvc::~DedxCorrecSvc ( )

Definition at line 51 of file DedxCorrecSvc.cxx.

51 {
52}

◆ DedxCorrecSvc() [2/2]

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

◆ ~DedxCorrecSvc() [2/2]

DedxCorrecSvc::~DedxCorrecSvc ( )

Member Function Documentation

◆ CellCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 846 of file DedxCorrecSvc.cxx.

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

◆ CellCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ CosthetaCorrec() [1/2]

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

Definition at line 567 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ CosthetaCorrec() [2/2]

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

◆ D2I() [1/2]

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

Definition at line 2159 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ D2I() [2/2]

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

◆ DipAngCorrec() [1/2]

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

Definition at line 832 of file DedxCorrecSvc.cxx.

832 {
833}

Referenced by StandardCorrec().

◆ DipAngCorrec() [2/2]

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

◆ DocaSinCorrec() [1/2]

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

Definition at line 789 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ DocaSinCorrec() [2/2]

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

◆ DriftDistCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 698 of file DedxCorrecSvc.cxx.

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

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

◆ DriftDistCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ EntaCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 719 of file DedxCorrecSvc.cxx.

719 {
720 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl;
721// double dedx_enta;
722// if(m_par_flag == 1) {
723// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta +
724// m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3);
725// } else if(m_par_flag == 0) {
726// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) +
727// m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta);
728// }
729// dedx_enta = dedx/dedx_enta;
730// if (dedx_enta>0.0) return dedx_enta;
731// else return dedx;
732
733 if(eangle>-0.25 && eangle<0.25){
734 double stepsize= 0.5/m_venangle.size();
735 int nsize= m_venangle.size();
736 double slope=0;
737 double offset=1;
738 double y1=0,y2=0,x1=0,x2=0;
739
740 if(eangle>=-0.25+0.5*stepsize && eangle<=0.25-0.5*stepsize){
741 int bin = (int)( (eangle-(-0.25+0.5*stepsize))/stepsize);
742 y1=m_venangle[bin];
743 x1=-0.25+0.5*stepsize+bin*stepsize;
744 y2=m_venangle[bin+1];
745 x2=-0.25+1.5*stepsize+bin*stepsize;
746 }
747 else if(eangle<=-0.25+0.5*stepsize){
748 y1=m_venangle[0];
749 x1=-0.25+0.5*stepsize;
750 y2=m_venangle[1];
751 x2=-0.25+1.5*stepsize;
752 }
753 else if( eangle>=0.25-0.5*stepsize){
754 y1=m_venangle[nsize-2];
755 x1=0.25-1.5*stepsize;
756 y2=m_venangle[nsize-1];
757 x2=0.25-0.5*stepsize;
758 }
759 double kcof= (y2-y1)/(x2-x1);
760 double bcof= y1-kcof*x1;
761 double ratio = kcof*eangle+bcof;
762 dedx=dedx/ratio;
763 }
764
765 return dedx;
766}
bin
*******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
slope
void slope()
Definition: slope.cxx:12

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ EntaCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ f_larcos() [1/2]

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}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57

◆ f_larcos() [2/2]

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

◆ finalize() [1/2]

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().

◆ finalize() [2/2]

virtual StatusCode DedxCorrecSvc::finalize ( )
virtual

◆ GlobalCorrec() [1/2]

double DedxCorrecSvc::GlobalCorrec ( double  dedx) const
virtual

Implements IDedxCorrecSvc.

Definition at line 836 of file DedxCorrecSvc.cxx.

836 {
837 if( m_mdcg_flag == 0 ) return dedx;
838 double calib_ex = dedx;
839 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
840 return calib_ex;
841}

Referenced by StandardCorrec().

◆ GlobalCorrec() [2/2]

double DedxCorrecSvc::GlobalCorrec ( double  dedx) const
virtual

Implements IDedxCorrecSvc.

◆ HadronCorrec() [1/2]

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

Definition at line 653 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ HadronCorrec() [2/2]

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

◆ handle() [1/2]

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}

◆ handle() [2/2]

void DedxCorrecSvc::handle ( const Incident &  )

◆ I2D() [1/2]

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

Definition at line 2181 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ I2D() [2/2]

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

◆ initialize() [1/2]

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().

◆ initialize() [2/2]

virtual StatusCode DedxCorrecSvc::initialize ( )
virtual

◆ LayerCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 896 of file DedxCorrecSvc.cxx.

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

◆ LayerCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ LayerGainCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 660 of file DedxCorrecSvc.cxx.

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

Referenced by CellCorrec(), and StandardCorrec().

◆ LayerGainCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ PathL() [1/2]

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 1495 of file DedxCorrecSvc.cxx.

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

◆ PathL() [2/2]

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

Implements IDedxCorrecSvc.

◆ queryInterface() [1/2]

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}

◆ queryInterface() [2/2]

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

◆ RungCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 674 of file DedxCorrecSvc.cxx.

674 {
675 //cout<<"DedxCorrecSvc::RungCorrec"<<endl;
676 double dedx_rung;
677 int run_ture =0;
678 // cout << " runNO : "<<runNO << " evtNO " << evtNO << endl;
679
680 for(int j=0;j<100000;j++) {
681 // for(int j=0;j<10;j++) {
682 if((runNO == m_rung[2][j]) && (evtNO >= m_rung[4][j]) && (evtNO <= m_rung[5][j])) {
683 dedx_rung = dedx/m_rung[0][j];
684 // cout <<"j " << j << "runno " << m_rung[2][j] << " evtNO " << evtNO << " evtfrom " << m_rung[4][j] << " evtto " << m_rung[5][j] << " rung " << m_rung[0][j] <<endl;
685 run_ture =1;
686 return dedx_rung;
687 }
688 }
689 if(run_ture ==0)
690 {
691 cout<<"Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to "<<runNO<<endl;
692 exit(1);
693 }
694}

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

◆ RungCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ SaturCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 547 of file DedxCorrecSvc.cxx.

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

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

◆ SaturCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ set_flag() [1/2]

void DedxCorrecSvc::set_flag ( int  calib_F)
virtual

Implements IDedxCorrecSvc.

Definition at line 1137 of file DedxCorrecSvc.cxx.

1137 {
1138 // cout<<"DedxCorrecSvc::set_flag"<<endl;
1139 cout<<"calib_F is = "<<calib_F<<endl;
1140 if(calib_F<0){ m_enta_flag = 0; calib_F = abs(calib_F); }
1141 else m_enta_flag = 1;
1142 m_rung_flag = ( calib_F & 1 )? 1 : 0;
1143 m_wireg_flag = ( calib_F & 2 )? 1 : 0;
1144 m_dcasin_flag = ( calib_F & 4 )? 1 : 0;
1145 m_ggs_flag = ( calib_F & 8 )? 1 : 0;
1146 m_ddg_flag = 0;
1147 //m_ddg_flag = ( calib_F & 8 )? 1 : 0;
1148 m_t0_flag = ( calib_F & 16 )? 1 : 0;
1149 m_sat_flag = ( calib_F & 32 )? 1 : 0;
1150 m_layerg_flag = ( calib_F & 64 )? 1 : 0;
1151 //m_dcasin_flag = ( calib_F & 128 )? 1 : 0;
1152 m_dip_flag = ( calib_F & 256 )? 1 : 0;
1153 m_mdcg_flag = ( calib_F & 512 )? 1 : 0;
1154}

◆ set_flag() [2/2]

void DedxCorrecSvc::set_flag ( int  calib_F)
virtual

Implements IDedxCorrecSvc.

◆ StandardCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 923 of file DedxCorrecSvc.cxx.

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

◆ StandardCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ StandardHitCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 997 of file DedxCorrecSvc.cxx.

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

◆ StandardHitCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ StandardTrackCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 1062 of file DedxCorrecSvc.cxx.

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

◆ StandardTrackCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ T0Correc() [1/2]

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

Definition at line 621 of file DedxCorrecSvc.cxx.

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

Referenced by StandardTrackCorrec().

◆ T0Correc() [2/2]

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

◆ TrkCorrec() [1/2]

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

dEdx calib. for each track

Implements IDedxCorrecSvc.

Definition at line 908 of file DedxCorrecSvc.cxx.

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

◆ TrkCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ WireGainCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 535 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ WireGainCorrec() [2/2]

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

Implements IDedxCorrecSvc.

◆ ZdepCorrec() [1/2]

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

Implements IDedxCorrecSvc.

Definition at line 770 of file DedxCorrecSvc.cxx.

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

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

◆ ZdepCorrec() [2/2]

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

Implements IDedxCorrecSvc.

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