BOSS 7.1.2
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 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)
 

Detailed Description

Definition at line 18 of file DedxCorrecSvc.h.

Constructor & Destructor Documentation

◆ DedxCorrecSvc()

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

Definition at line 40 of file DedxCorrecSvc.cxx.

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

◆ ~DedxCorrecSvc()

DedxCorrecSvc::~DedxCorrecSvc ( )

Definition at line 52 of file DedxCorrecSvc.cxx.

52 {
53}

Member Function Documentation

◆ CellCorrec()

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

Definition at line 838 of file DedxCorrecSvc.cxx.

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

◆ CosthetaCorrec()

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

Definition at line 559 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ D2I()

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

Definition at line 2151 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ DipAngCorrec()

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

Definition at line 824 of file DedxCorrecSvc.cxx.

824 {
825}

Referenced by StandardCorrec().

◆ DocaSinCorrec()

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

Definition at line 781 of file DedxCorrecSvc.cxx.

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

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ DriftDistCorrec()

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

Definition at line 690 of file DedxCorrecSvc.cxx.

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

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

◆ EntaCorrec()

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

Definition at line 711 of file DedxCorrecSvc.cxx.

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

◆ f_larcos()

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

Definition at line 113 of file DedxCorrecSvc.cxx.

113 {
114 if(nbin!=80) return x; // temporally for only nbin = 80
115 double m_absx(fabs(x));
116 double m_charge(x/m_absx);
117 if(m_absx>0.925 && m_absx<=0.950) return 0.9283*m_charge; // these numbers are from the mean values
118 if(m_absx>0.900 && m_absx<=0.925) return 0.9115*m_charge; // of each bin in cos(theta) distribution
119 if(m_absx>0.875 && m_absx<=0.900) return 0.8877*m_charge;
120 if(m_absx>0.850 && m_absx<=0.875) return 0.8634*m_charge;
121 if(m_absx>0.825 && m_absx<=0.850) return 0.8460*m_charge;
122 if(m_absx>0.800 && m_absx<=0.825) return 0.8089*m_charge;
123}
x
Double_t x[10]
Definition DataBase/tau_mode.c:57

◆ finalize()

StatusCode DedxCorrecSvc::finalize ( )
virtual

Definition at line 86 of file DedxCorrecSvc.cxx.

86 {
87 // cout<<"DedxCorrecSvc::finalize()"<<endl;
88 MsgStream log(messageService(), name());
89 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq;
90 return StatusCode::SUCCESS;
91}

Referenced by main().

◆ GlobalCorrec()

double DedxCorrecSvc::GlobalCorrec ( double dedx) const

Definition at line 828 of file DedxCorrecSvc.cxx.

828 {
829 if( m_mdcg_flag == 0 ) return dedx;
830 double calib_ex = dedx;
831 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
832 return calib_ex;
833}

Referenced by StandardCorrec().

◆ HadronCorrec()

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

Definition at line 645 of file DedxCorrecSvc.cxx.

645 {
646 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl;
647 //constant 1.00 in the following function is the mean dedx of normalized electrons.
648 dedx=dedx/550.00;
649 return D2I(costheta, I2D(costheta,1.00)/1.00*dedx)*550;
650}
DedxCorrecSvc::D2I
double D2I(const double &cosTheta, const double &D) const
Definition DedxCorrecSvc.cxx:2151
DedxCorrecSvc::I2D
double I2D(const double &cosTheta, const double &I) const
Definition DedxCorrecSvc.cxx:2173

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ handle()

void DedxCorrecSvc::handle ( const Incident & inc)

Definition at line 93 of file DedxCorrecSvc.cxx.

93 {
94 // cout<<"DedxCorrecSvc::handle"<<endl;
95 MsgStream log( messageService(), name() );
96 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
97
98 if ( inc.type() == "NewRun" ){
99 log << MSG::DEBUG << "New Run" << endreq;
100
101 if( ! m_initConstFlg )
102 {
103 if( m_init == 0 ) { init_param(); }
104 else if( m_init ==1 ) { init_param_svc(); }
105 /* if( ! init_param() ){
106 log << MSG::ERROR
107 << "can not initilize Mdc Calib Constants" << endreq;
108 }*/
109 }
110 }
111}

◆ I2D()

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

Definition at line 2173 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ initialize()

StatusCode DedxCorrecSvc::initialize ( )
virtual

Definition at line 56 of file DedxCorrecSvc.cxx.

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

Referenced by main().

◆ LayerCorrec()

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

Definition at line 888 of file DedxCorrecSvc.cxx.

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

◆ LayerGainCorrec()

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

Definition at line 652 of file DedxCorrecSvc.cxx.

652 {
653 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl;
654 if( m_layer_g[layid] > 0 ){
655 double ch_dedx = dedx/m_layer_g[layid];
656 return ch_dedx;
657 }
658 else if( m_layer_g[layid] == 0 ){
659 return dedx;
660 }
661 else return 0;
662}

Referenced by CellCorrec(), and StandardCorrec().

◆ PathL()

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

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

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

assumed the first half circle

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

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

in the Local Helix case, phi must be small

Definition at line 1487 of file DedxCorrecSvc.cxx.

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

◆ RungCorrec()

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

Definition at line 666 of file DedxCorrecSvc.cxx.

666 {
667 //cout<<"DedxCorrecSvc::RungCorrec"<<endl;
668 double dedx_rung;
669 int run_ture =0;
670 // cout << " runNO : "<<runNO << " evtNO " << evtNO << endl;
671
672 for(int j=0;j<100000;j++) {
673 // for(int j=0;j<10;j++) {
674 if((runNO == m_rung[2][j]) && (evtNO >= m_rung[4][j]) && (evtNO <= m_rung[5][j])) {
675 dedx_rung = dedx/m_rung[0][j];
676 // 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;
677 run_ture =1;
678 return dedx_rung;
679 }
680 }
681 if(run_ture ==0)
682 {
683 cout<<"Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to "<<runNO<<endl;
684 exit(1);
685 }
686}

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

◆ SaturCorrec()

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

Definition at line 539 of file DedxCorrecSvc.cxx.

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

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

◆ set_flag()

void DedxCorrecSvc::set_flag ( int calib_F)

Definition at line 1129 of file DedxCorrecSvc.cxx.

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

◆ StandardCorrec()

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

Definition at line 915 of file DedxCorrecSvc.cxx.

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

◆ StandardHitCorrec()

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

Definition at line 989 of file DedxCorrecSvc.cxx.

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

◆ StandardTrackCorrec()

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

Definition at line 1054 of file DedxCorrecSvc.cxx.

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

◆ T0Correc()

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

Definition at line 613 of file DedxCorrecSvc.cxx.

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

Referenced by StandardTrackCorrec().

◆ TrkCorrec()

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

dEdx calib. for each track

Definition at line 900 of file DedxCorrecSvc.cxx.

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

◆ WireGainCorrec()

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

Definition at line 527 of file DedxCorrecSvc.cxx.

527 {
528 if( m_wire_g[wireid] > 0 ){
529 double ch_dedx = (ex/m_wire_g[wireid])*m_valid[wireid];
530 return ch_dedx;
531 }
532 else if( m_wire_g[wireid] == 0 ){
533 return ex;
534 }
535 else return 0;
536}

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ ZdepCorrec()

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

Definition at line 762 of file DedxCorrecSvc.cxx.

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

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

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