BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCorrecSvc Class Reference

#include <DedxCorrecSvc.h>

Inheritance diagram for DedxCorrecSvc:

Public Member Functions

 DedxCorrecSvc (const std::string &name, ISvcLocator *svcloc)
 ~DedxCorrecSvc ()
virtual StatusCode 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 20 of file DedxCorrecSvc.h.

Constructor & Destructor Documentation

◆ DedxCorrecSvc()

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

Definition at line 38 of file DedxCorrecSvc.cxx.

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

Referenced by DedxCorrecSvc().

◆ ~DedxCorrecSvc()

DedxCorrecSvc::~DedxCorrecSvc ( )

Definition at line 50 of file DedxCorrecSvc.cxx.

50{}

Member Function Documentation

◆ CellCorrec()

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

Definition at line 868 of file DedxCorrecSvc.cxx.

869 {
870
871 if ( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0 && m_layerg_flag == 0 &&
872 m_zdep_flag == 0 && m_ggs_flag == 0 )
873 return adc;
874 adc = m_valid[ser] * m_wire_g[ser] * adc;
875 // int lyr = Wire2Lyr(ser);
876 int lyr = ser;
877 double ex = adc;
878 double correct1_ex, correct2_ex, correct3_ex, correct4_ex, correct5_ex;
879
880 if ( m_ggs_flag )
881 {
882 correct1_ex = SaturCorrec( lyr, costheta, adc );
883 ex = correct1_ex;
884 }
885 else { correct1_ex = adc; }
886
887 if ( m_ddg_flag )
888 {
889 correct2_ex = DriftDistCorrec( lyr, dd, correct1_ex );
890 ex = correct2_ex;
891 }
892 else { correct2_ex = correct1_ex; }
893 if ( m_enta_flag )
894 {
895 correct3_ex = DriftDistCorrec( lyr, sinenta, correct2_ex );
896 ex = correct3_ex;
897 }
898 else { correct3_ex = correct2_ex; }
899
900 if ( m_zdep_flag )
901 {
902 correct4_ex = ZdepCorrec( lyr, z, correct3_ex );
903 ex = correct4_ex;
904 }
905 else { correct4_ex = correct3_ex; }
906
907 if ( m_layerg_flag )
908 {
909 correct5_ex = LayerGainCorrec( lyr, correct4_ex );
910 ex = correct5_ex;
911 }
912 else { correct5_ex = correct4_ex; }
913 return ex;
914}
double LayerGainCorrec(int layid, double ex) const
double SaturCorrec(int layid, double costheta, double ex) const
double DriftDistCorrec(int layid, double ddrift, double ex) const
double ZdepCorrec(int layid, double zhit, double ex) const

◆ CosthetaCorrec()

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

Definition at line 594 of file DedxCorrecSvc.cxx.

594 {
595 // cout<<"DedxCorrecSvc::CosthetaCorrec"<<endl;
596 double dedx_cos;
597 // cout<<"costheta vaule is = "<<costheta<< " dedx = " << dedx << endl;
598 if ( fabs( costheta ) > 1.0 ) return dedx;
599
600 const int nbin = cos_k.size() + 1;
601 const double step = 2.00 / nbin;
602 int n = (int)( ( costheta + 1.00 - 0.5 * step ) / step );
603 if ( costheta > ( 1.00 - 0.5 * step ) ) n = nbin - 2;
604
605 if ( costheta > -0.5 * step && costheta <= 0 )
606 dedx_cos = cos_k[n - 1] * costheta + cos_b[n - 1];
607 else if ( costheta > 0 && costheta < 0.5 * step )
608 dedx_cos = cos_k[n + 1] * costheta + cos_b[n + 1];
609 else dedx_cos = cos_k[n] * costheta + cos_b[n];
610
611 // cout << "cos_k[n]="<<cos_k[n] << " cos_b[n]=" << cos_b[n] <<
612 // " dedx_cos=" << dedx_cos << " final dedx=" << dedx/dedx_cos << endl;
613
614 // conside physical edge around 0.93
615 if ( nbin == 80 )
616 { // temporally for only nbin = 80
617 if ( costheta > 0.80 && costheta < 0.95 )
618 {
619 n = 77;
620 if ( costheta < 0.9282 ) n--;
621 if ( costheta < 0.9115 ) n--;
622 if ( costheta < 0.8877 ) n--;
623 if ( costheta < 0.8634 ) n--;
624 if ( costheta < 0.8460 ) n--;
625 if ( costheta < 0.8089 ) n--;
626 dedx_cos = cos_k[n] * costheta + cos_b[n];
627 }
628 else if ( costheta < -0.80 && costheta > -0.95 )
629 {
630 n = 2;
631 if ( costheta > -0.9115 ) n++;
632 if ( costheta > -0.8877 ) n++;
633 if ( costheta > -0.8634 ) n++;
634 if ( costheta > -0.8460 ) n++;
635 if ( costheta > -0.8089 ) n++;
636 dedx_cos = cos_k[n] * costheta + cos_b[n];
637 }
638 }
639
640 if ( dedx_cos > 0 )
641 {
642 dedx_cos = dedx / dedx_cos;
643 return dedx_cos;
644 }
645 else return dedx;
646}
const Int_t n

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ D2I()

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

Definition at line 2166 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ DipAngCorrec()

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

Definition at line 856 of file DedxCorrecSvc.cxx.

856 {
857 std::cerr << "Error: DipAngCorrec is not implemented yet!" << std::endl;
858 exit( 1 );
859}

Referenced by StandardCorrec().

◆ DocaSinCorrec()

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

Definition at line 809 of file DedxCorrecSvc.cxx.

810 {
811 if ( m_debug ) cout << "DedxCorrecSvc::DocaSinCorrec" << endl;
812 // cout<<"doca = "<<doca<<" eangle = "<<eangle<<endl;
813
814 if ( eangle > 0.25 ) eangle = eangle - 0.5;
815 else if ( eangle < -0.25 ) eangle = eangle + 0.5;
816 int iDoca;
817 int iEAng = 0;
818 doca = doca / doca_norm[layer];
819 iDoca = (Int_t)floor( ( doca - Out_DocaMin ) / Out_Stepdoca );
820 if ( doca < Out_DocaMin ) iDoca = 0;
821 if ( doca > Out_DocaMax ) iDoca = NumSlicesDoca - 1;
822 if ( iDoca >= (Int_t)floor( ( Out_DocaMax - Out_DocaMin ) / Out_Stepdoca ) )
823 iDoca = (Int_t)floor( ( Out_DocaMax - Out_DocaMin ) / Out_Stepdoca ) - 1;
824 if ( iDoca < 0 ) iDoca = 0;
825 if ( m_debug )
826 cout << "iDoca : " << iDoca << " doca : " << doca
827 << " Out_DocaMin : " << Out_DocaMin << " Out_Stepdoca : " << Out_Stepdoca
828 << endl;
829
830 // iEAng = (Int_t)floor((eangle-Phi_Min)/IEangle_step);
831 for ( int i = 0; i < 14; i++ )
832 {
833 // iEAng =0;
834 if ( eangle < Eangle_value_cut[i] || eangle > Eangle_value_cut[i + 1] ) continue;
835 else if ( eangle >= Eangle_value_cut[i] && eangle < Eangle_value_cut[i + 1] )
836 {
837 for ( int k = 0; k < i; k++ ) { iEAng += Eangle_cut_bin[k]; }
838 double eangle_bin_cut_step =
839 ( Eangle_value_cut[i + 1] - Eangle_value_cut[i] ) / Eangle_cut_bin[i];
840 int temp_bin = int( ( eangle - Eangle_value_cut[i] ) / eangle_bin_cut_step );
841 iEAng += temp_bin;
842 }
843 }
844 // cout<<iDoca <<" "<<iEAng<<endl;
845 if ( m_docaeangle[iDoca][iEAng] > 0 )
846 {
847 if ( m_debug && ( iDoca == m_debug_i ) && ( iEAng == m_debug_j ) )
848 cout << "doca: " << doca << " eang" << eangle << " dedx before: " << dedx << endl;
849 dedx = dedx / m_docaeangle[iDoca][iEAng];
850 if ( m_debug && ( iDoca == m_debug_i ) && ( iEAng == m_debug_j ) )
851 cout << "gain: " << m_docaeangle[iDoca][iEAng] << " dedx after: " << dedx << endl;
852 }
853 return dedx;
854}

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ DriftDistCorrec()

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

Definition at line 714 of file DedxCorrecSvc.cxx.

714 {
715 // cout<<"DedxCorrecSvc::DriftDistCorrec"<<endl;
716 double dedx_ddg=1;
717 // cout<<"m_par_flag = "<<m_par_flag<<endl;
718 // cout<<"dd vaule is = "<<dd<<endl;
719 dd = fabs( dd );
720 if ( m_par_flag == 1 )
721 {
722 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer] * dd + m_ddg[2][layer] * dd * dd +
723 m_ddg[3][layer] * pow( dd, 3 );
724 }
725 else if ( m_par_flag == 0 )
726 {
727 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer] * T1( dd ) + m_ddg[2][layer] * T2( dd ) +
728 m_ddg[3][layer] * T3( dd );
729 }
730 // cout<<"dedx_ddg is = "<<dedx_ddg<<endl;
731 dedx_ddg = dedx / dedx_ddg;
732 if ( dedx_ddg > 0.0 ) return dedx_ddg;
733 else return dedx;
734}

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

◆ EntaCorrec()

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

Definition at line 736 of file DedxCorrecSvc.cxx.

736 {
737 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl;
738 // double dedx_enta;
739 // if(m_par_flag == 1) {
740 // dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta +
741 // m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3);
742 // } else if(m_par_flag == 0) {
743 // dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) +
744 // m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta);
745 // }
746 // dedx_enta = dedx/dedx_enta;
747 // if (dedx_enta>0.0) return dedx_enta;
748 // else return dedx;
749
750 if ( eangle > -0.25 && eangle < 0.25 )
751 {
752 double stepsize = 0.5 / m_venangle.size();
753 int nsize = m_venangle.size();
754 double slope = 0;
755 double offset = 1;
756 double y1 = 0, y2 = 0, x1 = 0, x2 = 0;
757
758 if ( eangle >= -0.25 + 0.5 * stepsize && eangle <= 0.25 - 0.5 * stepsize )
759 {
760 int bin = (int)( ( eangle - ( -0.25 + 0.5 * stepsize ) ) / stepsize );
761 y1 = m_venangle[bin];
762 x1 = -0.25 + 0.5 * stepsize + bin * stepsize;
763 y2 = m_venangle[bin + 1];
764 x2 = -0.25 + 1.5 * stepsize + bin * stepsize;
765 }
766 else if ( eangle <= -0.25 + 0.5 * stepsize )
767 {
768 y1 = m_venangle[0];
769 x1 = -0.25 + 0.5 * stepsize;
770 y2 = m_venangle[1];
771 x2 = -0.25 + 1.5 * stepsize;
772 }
773 else if ( eangle >= 0.25 - 0.5 * stepsize )
774 {
775 y1 = m_venangle[nsize - 2];
776 x1 = 0.25 - 1.5 * stepsize;
777 y2 = m_venangle[nsize - 1];
778 x2 = 0.25 - 0.5 * stepsize;
779 }
780 double kcof = ( y2 - y1 ) / ( x2 - x1 );
781 double bcof = y1 - kcof * x1;
782 double ratio = kcof * eangle + bcof;
783 dedx = dedx / ratio;
784 }
785
786 return dedx;
787}
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ f_larcos()

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

Definition at line 111 of file DedxCorrecSvc.cxx.

111 {
112 if ( nbin != 80 ) return x; // temporally for only nbin = 80
113 double m_absx( fabs( x ) );
114 double m_charge( x / m_absx );
115 if ( m_absx > 0.925 && m_absx <= 0.950 )
116 return 0.9283 * m_charge; // these numbers are from the mean values
117 if ( m_absx > 0.900 && m_absx <= 0.925 )
118 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 return x;
124 std::cerr << "DeDxCorrecSvc::f_larcos: should not reach here!" << std::endl;
125 exit( 1 );
126}
Double_t x[10]

◆ finalize()

StatusCode DedxCorrecSvc::finalize ( )
virtual

Definition at line 83 of file DedxCorrecSvc.cxx.

83 {
84 // cout<<"DedxCorrecSvc::finalize()"<<endl;
85 MsgStream log( msgSvc(), name() );
86 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endmsg;
87 return StatusCode::SUCCESS;
88}
IMessageSvc * msgSvc()

Referenced by main().

◆ GlobalCorrec()

double DedxCorrecSvc::GlobalCorrec ( double dedx) const

Definition at line 861 of file DedxCorrecSvc.cxx.

861 {
862 if ( m_mdcg_flag == 0 ) return dedx;
863 double calib_ex = dedx;
864 if ( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
865 return calib_ex;
866}

Referenced by StandardCorrec().

◆ HadronCorrec()

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

Definition at line 672 of file DedxCorrecSvc.cxx.

672 {
673 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl;
674 // constant 1.00 in the following function is the mean dedx of normalized electrons.
675 dedx = dedx / 550.00;
676 return D2I( costheta, I2D( costheta, 1.00 ) / 1.00 * dedx ) * 550;
677}
double D2I(const double &cosTheta, const double &D) const
double I2D(const double &cosTheta, const double &I) const

Referenced by StandardCorrec(), and StandardTrackCorrec().

◆ handle()

void DedxCorrecSvc::handle ( const Incident & inc)

Definition at line 90 of file DedxCorrecSvc.cxx.

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

◆ I2D()

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

Definition at line 2187 of file DedxCorrecSvc.cxx.

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

Referenced by HadronCorrec().

◆ initialize()

StatusCode DedxCorrecSvc::initialize ( )
virtual

Definition at line 52 of file DedxCorrecSvc.cxx.

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

Referenced by main().

◆ LayerCorrec()

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

Definition at line 916 of file DedxCorrecSvc.cxx.

916 {
917 // cout<<"DedxCorrecSvc::LayerCorrec"<<endl;
918
919 if ( m_zdep_flag == 0 && m_ggs_flag == 0 ) return ex;
920
921 double calib_ex = ZdepCorrec( layer, z, ex );
922 if ( m_ggs_flag != 0 ) calib_ex = SaturCorrec( layer, costheta, calib_ex );
923 return calib_ex;
924}

◆ LayerGainCorrec()

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

Definition at line 678 of file DedxCorrecSvc.cxx.

678 {
679 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl;
680 if ( m_layer_g[layid] > 0 )
681 {
682 double ch_dedx = dedx / m_layer_g[layid];
683 return ch_dedx;
684 }
685 else if ( m_layer_g[layid] == 0 ) { return dedx; }
686 else return 0;
687}

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

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

◆ RungCorrec()

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

Definition at line 689 of file DedxCorrecSvc.cxx.

689 {
690 // cout<<"DedxCorrecSvc::RungCorrec"<<endl;
691 double dedx_rung;
692 int run_ture = 0;
693 // cout << " runNO : "<<runNO << " evtNO " << evtNO << endl;
694
695 for ( int j = 0; j < 100000; j++ )
696 {
697 // for(int j=0;j<10;j++) {
698 if ( ( runNO == m_rung[2][j] ) && ( evtNO >= m_rung[4][j] ) && ( evtNO <= m_rung[5][j] ) )
699 {
700 dedx_rung = dedx / m_rung[0][j];
701 // cout <<"j " << j << "runno " << m_rung[2][j] << " evtNO " << evtNO << " evtfrom "
702 // << m_rung[4][j] << " evtto " << m_rung[5][j] << " rung " << m_rung[0][j] <<endl;
703 run_ture = 1;
704 return dedx_rung;
705 }
706 }
707 if ( run_ture == 0 )
708 { cout << "Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to " << runNO << endl; }
709
710 std::cerr << __FILE__ << ":" << __LINE__ << " Should not reach here!" << std::endl;
711 exit( 1 );
712}

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

◆ SaturCorrec()

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

Definition at line 573 of file DedxCorrecSvc.cxx.

573 {
574 // cout<<"DedxCorrecSvc::SaturCorrec"<<endl;
575 double dedx_ggs=1;
576 // cout<<"costheta vaule is = "<<costheta<<endl;
577 costheta = fabs( costheta );
578 if ( m_par_flag == 1 )
579 {
580 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer] * costheta +
581 m_ggs[2][layer] * pow( costheta, 2 ) + m_ggs[3][layer] * pow( costheta, 3 );
582 }
583 else if ( m_par_flag == 0 )
584 {
585 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer] * T1( costheta ) +
586 m_ggs[2][layer] * T2( costheta ) + m_ggs[3][layer] * T3( costheta );
587 }
588 // cout<<"dedx_ggs is = "<<dedx_ggs<<endl;
589 dedx_ggs = dedx / dedx_ggs;
590 if ( dedx_ggs > 0.0 ) return dedx_ggs;
591 else return dedx;
592}

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

◆ set_flag()

void DedxCorrecSvc::set_flag ( int calib_F)

Definition at line 1120 of file DedxCorrecSvc.cxx.

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

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

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

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

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

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

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

◆ T0Correc()

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

Definition at line 648 of file DedxCorrecSvc.cxx.

648 {
649 // cout<<"DedxCorrecSvc::T0Correc"<<endl;
650 double dedx_t0;
651 const int nbin = t0_k.size() + 1;
652 if ( nbin != 35 ) cout << "there should be 35 bins for t0 correction, check it!" << endl;
653
654 int n = 0;
655 if ( t0 < 575 ) n = 0;
656 else if ( t0 < 610 ) n = 1;
657 else if ( t0 >= 610 && t0 < 1190 ) { n = (int)( ( t0 - 610.0 ) / 20.0 ) + 2; }
658 else if ( t0 >= 1190 && t0 < 1225 ) n = 31;
659 else if ( t0 >= 1225 && t0 < 1275 ) n = 32;
660 else if ( t0 >= 1275 ) n = 33;
661
662 dedx_t0 = t0_k[n] * t0 + t0_b[n];
663
664 if ( dedx_t0 > 0 )
665 {
666 dedx_t0 = dedx / dedx_t0;
667 return dedx_t0;
668 }
669 else return dedx;
670}

Referenced by StandardTrackCorrec().

◆ TrkCorrec()

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

dEdx calib. for each track

Definition at line 926 of file DedxCorrecSvc.cxx.

926 {
927 // cout<<"DedxCorrecSvc::TrkCorrec"<<endl;
928 if ( m_mdcg_flag == 0 ) return dedx;
929 /// dEdx calib. for each track
930 double calib_ex = dedx;
931
932 double run_const = m_dedx_gain;
933 if ( run_const > 0 && m_mdcg_flag != 0 ) calib_ex /= run_const;
934
935 return calib_ex;
936}

◆ WireGainCorrec()

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

Definition at line 563 of file DedxCorrecSvc.cxx.

563 {
564 if ( m_wire_g[wireid] > 0 )
565 {
566 double ch_dedx = ( ex / m_wire_g[wireid] ) * m_valid[wireid];
567 return ch_dedx;
568 }
569 else if ( m_wire_g[wireid] == 0 ) { return ex; }
570 else return 0;
571}

Referenced by StandardCorrec(), and StandardHitCorrec().

◆ ZdepCorrec()

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

Definition at line 789 of file DedxCorrecSvc.cxx.

789 {
790 // cout<<"DedxCorrecSvc::ZdepCorrec"<<endl;
791 double dedx_zdep=1;
792 if ( m_par_flag == 1 )
793 {
794 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer] * z + m_zdep[2][layer] * z * z +
795 m_zdep[3][layer] * pow( z, 3 );
796 }
797 else if ( m_par_flag == 0 )
798 {
799 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer] * T1( z ) + m_zdep[2][layer] * T2( z ) +
800 m_zdep[3][layer] * T3( z );
801 }
802 dedx_zdep = dedx / dedx_zdep;
803 if ( dedx_zdep > 0.0 ) return dedx_zdep;
804 else return dedx;
805
806 // return dedx;
807}

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


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