Loading...
Searching...
No Matches
DedxCorrecSvc Class Reference
#include <DedxCorrecSvc.h>
Inheritance diagram for DedxCorrecSvc:
Public Member Functions | |
| DedxCorrecSvc (const std::string &name, ISvcLocator *svcloc) | |
| ~DedxCorrecSvc () | |
| virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvUnknown) |
| virtual StatusCode | initialize () |
| virtual StatusCode | finalize () |
| void | handle (const Incident &) |
| double | RungCorrec (int runNO, double ex) const |
| double | WireGainCorrec (int wireid, double ex) const |
| double | DriftDistCorrec (int layid, double ddrift, double ex) const |
| double | SaturCorrec (int layid, double costheta, double ex) const |
| double | CosthetaCorrec (double costheta, double ex) const |
| double | T0Correc (double t0, double ex) const |
| double | HadronCorrec (double costheta, double dedx) const |
| double | ZdepCorrec (int layid, double zhit, double ex) const |
| double | EntaCorrec (int layid, double enta, double ex) const |
| double | LayerGainCorrec (int layid, double ex) const |
| double | DocaSinCorrec (int layid, double doca, double enta, double ex) const |
| double | DipAngCorrec (int layid, double costheta, double ex) const |
| double | GlobalCorrec (double dedx) const |
| double | CellCorrec (int ser, double adc, double dd, double enta, double z, double costheta) const |
| double | LayerCorrec (int layer, double z, double costheta, double ex) const |
| double | TrkCorrec (double costheta, double dedx) const |
| double | StandardCorrec (int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta) const |
| double | StandardHitCorrec (int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const |
| double | StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const |
| double | PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const |
| double | D2I (const double &cosTheta, const double &D) const |
| double | I2D (const double &cosTheta, const double &I) const |
| void | set_flag (int calib_F) |
| double | f_larcos (double x, int nbin) |
| virtual double | RungCorrec (int runNO, double ex) const =0 |
| virtual double | WireGainCorrec (int wireid, double ex) const =0 |
| virtual double | DriftDistCorrec (int layid, double ddrift, double ex) const =0 |
| virtual double | SaturCorrec (int layid, double costheta, double ex) const =0 |
| virtual double | EntaCorrec (int layid, double enta, double ex) const =0 |
| virtual double | ZdepCorrec (int layer, double z, double dedx) const =0 |
| virtual double | LayerGainCorrec (int layid, double dedx) const =0 |
| virtual double | GlobalCorrec (double dedx) const =0 |
| virtual double | CellCorrec (int ser, double adc, double dd, double enta, double z, double theta) const =0 |
| virtual double | LayerCorrec (int layer, double z, double costheta, double ex) const =0 |
| virtual double | TrkCorrec (double costheta, double dedx) const =0 |
| virtual double | StandardCorrec (int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double z, double costheta) const =0 |
| virtual double | StandardHitCorrec (int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0 |
| virtual double | StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const =0 |
| virtual double | PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const =0 |
| virtual void | set_flag (int calib_F)=0 |
Additional Inherited Members | |
| Static Public Member Functions inherited from IDedxCorrecSvc | |
| static const InterfaceID & | interfaceID () |
Detailed Description
Definition at line 19 of file DedxCorrecSvc.h.
Constructor & Destructor Documentation
◆ DedxCorrecSvc()
| DedxCorrecSvc::DedxCorrecSvc | ( | const std::string & | name, |
| ISvcLocator * | svcloc | ||
| ) |
Definition at line 39 of file DedxCorrecSvc.cxx.
39 :
40 geosvc(0),Service (name, svcloc) {
41 declareProperty("Run",m_run=1);
42 declareProperty("init",m_init=1);
43 declareProperty("par_flag",m_par_flag=0);
44 declareProperty("Debug",m_debug=false);
45 declareProperty("DebugI",m_debug_i=1);
46 declareProperty("DebugJ",m_debug_j=1);
47 m_initConstFlg = false;
48 // cout<<"DedxCorrecSvc::DedxCorrecSvc"<<endl;
49 }
◆ ~DedxCorrecSvc()
| DedxCorrecSvc::~DedxCorrecSvc | ( | ) |
Definition at line 51 of file DedxCorrecSvc.cxx.
51 {
52}
Member Function Documentation
◆ CellCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 842 of file DedxCorrecSvc.cxx.
843 {
844
845 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
846 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
847 m_ggs_flag == 0) return adc;
848 adc = m_valid[ser]*m_wire_g[ser]*adc;
849 // int lyr = Wire2Lyr(ser);
850 int lyr = ser;
851 double ex = adc;
852 double correct1_ex , correct2_ex, correct3_ex, correct4_ex, correct5_ex;
853
854 if( m_ggs_flag) {
855 correct1_ex = SaturCorrec( lyr, costheta, adc );
856 ex = correct1_ex;
857 } else {
858 correct1_ex = adc;
859 }
860
861 if( m_ddg_flag) {
862 correct2_ex = DriftDistCorrec( lyr, dd, correct1_ex );
863 ex = correct2_ex;
864 } else {
865 correct2_ex =correct1_ex;
866 }
867 if( m_enta_flag) {
868 correct3_ex = DriftDistCorrec( lyr, sinenta, correct2_ex );
869 ex = correct3_ex;
870 } else {
871 correct3_ex =correct2_ex;
872 }
873
874 if( m_zdep_flag) {
875 correct4_ex = ZdepCorrec( lyr, z, correct3_ex );
876 ex = correct4_ex;
877 } else {
878 correct4_ex = correct3_ex;
879 }
880
881 if( m_layerg_flag) {
882 correct5_ex = LayerGainCorrec( lyr, correct4_ex );
883 ex = correct5_ex;
884 } else {
885 correct5_ex = correct4_ex;
886 }
887 return ex;
888
889}
double LayerGainCorrec(int layid, double ex) const
Definition: DedxCorrecSvc.cxx:658
double SaturCorrec(int layid, double costheta, double ex) const
Definition: DedxCorrecSvc.cxx:545
double DriftDistCorrec(int layid, double ddrift, double ex) const
Definition: DedxCorrecSvc.cxx:694
double ZdepCorrec(int layid, double zhit, double ex) const
Definition: DedxCorrecSvc.cxx:766
◆ CosthetaCorrec()
| double DedxCorrecSvc::CosthetaCorrec | ( | double | costheta, |
| double | ex | ||
| ) | const |
Definition at line 565 of file DedxCorrecSvc.cxx.
565 {
566 //cout<<"DedxCorrecSvc::CosthetaCorrec"<<endl;
567 double dedx_cos;
568 //cout<<"costheta vaule is = "<<costheta<< " dedx = " << dedx << endl;
569 if(fabs(costheta)>1.0) return dedx;
570
571 const int nbin=cos_k.size()+1;
572 const double step=2.00/nbin;
574 if(costheta>(1.00-0.5*step))
575 n=nbin-2;
576
577 if(costheta>-0.5*step && costheta<=0)
579 else if(costheta>0 && costheta<0.5*step)
581 else
583
584 //cout << "cos_k[n]="<<cos_k[n] << " cos_b[n]=" << cos_b[n] <<
585 // " dedx_cos=" << dedx_cos << " final dedx=" << dedx/dedx_cos << endl;
586
587 //conside physical edge around 0.93
588 if(nbin==80){ // temporally for only nbin = 80
589 if(costheta>0.80 && costheta<0.95){
590 n = 77;
598 }
599 else if(costheta<-0.80 && costheta>-0.95){
600 n = 2;
607 }
608 }
609
610 if(dedx_cos>0){
611 dedx_cos = dedx/dedx_cos;
612 return dedx_cos;
613 }
614 else return dedx;
615}
Referenced by StandardCorrec(), and StandardTrackCorrec().
◆ D2I()
| double DedxCorrecSvc::D2I | ( | const double & | cosTheta, |
| const double & | D | ||
| ) | const |
Definition at line 2153 of file DedxCorrecSvc.cxx.
2154{
2155 //cout<<"DedxCorrecSvc::D2I"<<endl;
2156 double absCosTheta = fabs(cosTheta);
2157 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2158 double chargeDensity = D/projection;
2159 double numerator = 1 + m_alpha*chargeDensity;
2160 double denominator = 1 + m_gamma*chargeDensity;
2162
2163 //if(denominator > 0.1)
2164
2165 I = D * m_ratio* numerator/denominator;
2166// cout << "m_alpha " << m_alpha << endl;
2167// cout << "m_gamma " << m_gamma << endl;
2168// cout << "m_delta " << m_delta << endl;
2169// cout << "m_power " << m_power << endl;
2170// cout << "m_ratio " << m_ratio << endl;
2172}
const DifComplex I
Referenced by HadronCorrec().
◆ DipAngCorrec()
| double DedxCorrecSvc::DipAngCorrec | ( | int | layid, |
| double | costheta, | ||
| double | ex | ||
| ) | const |
◆ DocaSinCorrec()
| double DedxCorrecSvc::DocaSinCorrec | ( | int | layid, |
| double | doca, | ||
| double | enta, | ||
| double | ex | ||
| ) | const |
Definition at line 785 of file DedxCorrecSvc.cxx.
785 {
786 if(m_debug) cout<<"DedxCorrecSvc::DocaSinCorrec"<<endl;
787 // cout<<"doca = "<<doca<<" eangle = "<<eangle<<endl;
788
789 if(eangle>0.25) eangle = eangle -0.5;
790 else if(eangle < -0.25) eangle = eangle +0.5;
791 int iDoca;
792 int iEAng = 0;
793 doca = doca/doca_norm[layer];
794 iDoca =(Int_t)floor((doca-Out_DocaMin)/Out_Stepdoca);
799 if(iDoca<0) iDoca = 0;
800 if(m_debug) cout<<"iDoca : "<<iDoca<<" doca : "<<doca<<" Out_DocaMin : "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<endl;
801
802 //iEAng = (Int_t)floor((eangle-Phi_Min)/IEangle_step);
803 for(int i =0; i<14; i++){
804 //iEAng =0;
805 if(eangle <Eangle_value_cut[i] || eangle > Eangle_value_cut[i+1]) continue;
806 else if(eangle>= Eangle_value_cut[i] && eangle < Eangle_value_cut[i+1]) {
807 for(int k =0; k<i; k++){
808 iEAng+= Eangle_cut_bin[k];
809 }
810 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
811 int temp_bin = int((eangle-Eangle_value_cut[i])/eangle_bin_cut_step);
812 iEAng+= temp_bin;
813 }
814 }
815 //cout<<iDoca <<" "<<iEAng<<endl;
816 if(m_docaeangle[iDoca][iEAng]>0) {
817 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
818 cout << "doca: " << doca << " eang" << eangle << " dedx before: " << dedx << endl;
819 dedx = dedx/m_docaeangle[iDoca][iEAng];
820 if(m_debug && (iDoca==m_debug_i) && (iEAng==m_debug_j))
821 cout << "gain: " << m_docaeangle[iDoca][iEAng] << " dedx after: " << dedx << endl;
822 }
823 return dedx;
824}
Referenced by StandardCorrec(), and StandardHitCorrec().
◆ DriftDistCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 694 of file DedxCorrecSvc.cxx.
694 {
695 // cout<<"DedxCorrecSvc::DriftDistCorrec"<<endl;
696 double dedx_ddg;
697 //cout<<"m_par_flag = "<<m_par_flag<<endl;
698 //cout<<"dd vaule is = "<<dd<<endl;
699 dd = fabs(dd);
700 if(m_par_flag == 1) {
701 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*dd +
702 m_ddg[2][layer]*dd*dd + m_ddg[3][layer]*pow(dd,3);
703 } else if(m_par_flag == 0) {
704 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*T1(dd) +
705 m_ddg[2][layer]*T2(dd) + m_ddg[3][layer]*T3(dd);
706 }
707 //cout<<"dedx_ddg is = "<<dedx_ddg<<endl;
708 dedx_ddg = dedx/dedx_ddg;
709 if (dedx_ddg>0.0) return dedx_ddg;
710 else return dedx;
711}
Referenced by CellCorrec(), StandardCorrec(), and StandardHitCorrec().
◆ EntaCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 715 of file DedxCorrecSvc.cxx.
715 {
716 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl;
717// double dedx_enta;
718// if(m_par_flag == 1) {
719// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta +
720// m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3);
721// } else if(m_par_flag == 0) {
722// dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) +
723// m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta);
724// }
725// dedx_enta = dedx/dedx_enta;
726// if (dedx_enta>0.0) return dedx_enta;
727// else return dedx;
728
729 if(eangle>-0.25 && eangle<0.25){
730 double stepsize= 0.5/m_venangle.size();
731 int nsize= m_venangle.size();
732 double slope=0;
733 double offset=1;
734 double y1=0,y2=0,x1=0,x2=0;
735
736 if(eangle>=-0.25+0.5*stepsize && eangle<=0.25-0.5*stepsize){
738 y1=m_venangle[bin];
739 x1=-0.25+0.5*stepsize+bin*stepsize;
740 y2=m_venangle[bin+1];
741 x2=-0.25+1.5*stepsize+bin*stepsize;
742 }
743 else if(eangle<=-0.25+0.5*stepsize){
744 y1=m_venangle[0];
745 x1=-0.25+0.5*stepsize;
746 y2=m_venangle[1];
747 x2=-0.25+1.5*stepsize;
748 }
749 else if( eangle>=0.25-0.5*stepsize){
750 y1=m_venangle[nsize-2];
751 x1=0.25-1.5*stepsize;
752 y2=m_venangle[nsize-1];
753 x2=0.25-0.5*stepsize;
754 }
755 double kcof= (y2-y1)/(x2-x1);
756 double bcof= y1-kcof*x1;
757 double ratio = kcof*eangle+bcof;
758 dedx=dedx/ratio;
759 }
760
761 return dedx;
762}
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
Referenced by StandardCorrec(), and StandardHitCorrec().
◆ f_larcos()
| double DedxCorrecSvc::f_larcos | ( | double | x, |
| int | nbin | ||
| ) |
Definition at line 121 of file DedxCorrecSvc.cxx.
121 {
123 double m_absx(fabs(x));
124 double m_charge(x/m_absx);
125 if(m_absx>0.925 && m_absx<=0.950) return 0.9283*m_charge; // these numbers are from the mean values
126 if(m_absx>0.900 && m_absx<=0.925) return 0.9115*m_charge; // of each bin in cos(theta) distribution
127 if(m_absx>0.875 && m_absx<=0.900) return 0.8877*m_charge;
128 if(m_absx>0.850 && m_absx<=0.875) return 0.8634*m_charge;
129 if(m_absx>0.825 && m_absx<=0.850) return 0.8460*m_charge;
130 if(m_absx>0.800 && m_absx<=0.825) return 0.8089*m_charge;
131}
◆ finalize()
|
virtual |
Definition at line 94 of file DedxCorrecSvc.cxx.
94 {
95 // cout<<"DedxCorrecSvc::finalize()"<<endl;
96 MsgStream log(messageService(), name());
97 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq;
98 return StatusCode::SUCCESS;
99}
Referenced by main().
◆ GlobalCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 832 of file DedxCorrecSvc.cxx.
832 {
833 if( m_mdcg_flag == 0 ) return dedx;
834 double calib_ex = dedx;
835 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain;
836 return calib_ex;
837}
Referenced by StandardCorrec().
◆ HadronCorrec()
| double DedxCorrecSvc::HadronCorrec | ( | double | costheta, |
| double | dedx | ||
| ) | const |
Definition at line 651 of file DedxCorrecSvc.cxx.
651 {
652 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl;
653 //constant 1.00 in the following function is the mean dedx of normalized electrons.
654 dedx=dedx/550.00;
656}
double D2I(const double &cosTheta, const double &D) const
Definition: DedxCorrecSvc.cxx:2153
double I2D(const double &cosTheta, const double &I) const
Definition: DedxCorrecSvc.cxx:2175
Referenced by StandardCorrec(), and StandardTrackCorrec().
◆ handle()
| void DedxCorrecSvc::handle | ( | const Incident & | inc | ) |
Definition at line 101 of file DedxCorrecSvc.cxx.
101 {
102 // cout<<"DedxCorrecSvc::handle"<<endl;
103 MsgStream log( messageService(), name() );
104 log << MSG::DEBUG << "handle: " << inc.type() << endreq;
105
106 if ( inc.type() == "NewRun" ){
107 log << MSG::DEBUG << "New Run" << endreq;
108
109 if( ! m_initConstFlg )
110 {
111 if( m_init == 0 ) { init_param(); }
112 else if( m_init ==1 ) { init_param_svc(); }
113 /* if( ! init_param() ){
114 log << MSG::ERROR
115 << "can not initilize Mdc Calib Constants" << endreq;
116 }*/
117 }
118 }
119}
◆ I2D()
| double DedxCorrecSvc::I2D | ( | const double & | cosTheta, |
| const double & | I | ||
| ) | const |
Definition at line 2175 of file DedxCorrecSvc.cxx.
2176{
2177 // cout<<" DedxCorrecSvc::I2D"<<endl;
2178 double absCosTheta = fabs(cosTheta);
2179 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire
2180
2181 // Do the quadratic to compute d of the electron
2182 double a = m_alpha / projection;
2185
2186 if(b==0 && a==0){
2187 cout<<"both a and b coefficiants for hadron correction are 0"<<endl;
2189 }
2190
2191 double D = a != 0? (-b + sqrt(b*b - 4.0*a*c))/(2.0*a) : -c/b;
2192 if(D<0){
2193 cout<<"D is less 0! will try another solution"<<endl;
2194 D=a != 0? (-b - sqrt(b*b + 4.0*a*c))/(2.0*a) : -c/b;
2195 if(D<0){
2196 cout<<"D is still less 0! just return uncorrectecd value"<<endl;
2198 }
2199 }
2200
2201 return D;
2202}
Referenced by HadronCorrec().
◆ initialize()
|
virtual |
Definition at line 64 of file DedxCorrecSvc.cxx.
64 {
65 // cout<<"DedxCorrecSvc::initialize"<<endl;
66 MsgStream log(messageService(), name());
67 log << MSG::INFO << name() << "DedxCorrecSvc::initialize()" << endreq;
68
69 StatusCode sc = Service::initialize();
70 if( sc.isFailure() ) return sc;
71
72 IIncidentSvc* incsvc;
73 sc = service("IncidentSvc", incsvc);
74 int priority = 100;
75 if( sc.isSuccess() ){
76 //incsvc -> addListener(this, "BeginEvent", priority);
77 incsvc -> addListener(this, "NewRun", priority);
78 }
79
80 StatusCode scgeo = service("MdcGeomSvc", geosvc);
81 if(scgeo.isFailure() ) {
82 log << MSG::ERROR << "GeoSvc failing!!!!!!!" << endreq;
83 return scgeo;
84 }
85
86 StatusCode scmgn = service ("MagneticFieldSvc",m_pIMF);
87 if(scmgn!=StatusCode::SUCCESS) {
88 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
89 }
90
91 return StatusCode::SUCCESS;
92}
Referenced by main().
◆ LayerCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 892 of file DedxCorrecSvc.cxx.
892 {
893 //cout<<"DedxCorrecSvc::LayerCorrec"<<endl;
894
895 if( m_zdep_flag == 0 && m_ggs_flag == 0 ) return ex;
896
899 return calib_ex;
900
901}
◆ LayerGainCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 658 of file DedxCorrecSvc.cxx.
658 {
659 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl;
660 if( m_layer_g[layid] > 0 ){
661 double ch_dedx = dedx/m_layer_g[layid];
662 return ch_dedx;
663 }
664 else if( m_layer_g[layid] == 0 ){
665 return dedx;
666 }
667 else return 0;
668}
Referenced by CellCorrec(), and StandardCorrec().
◆ PathL()
|
virtual |
***********----------------—***************---------------—****************//
***********----------------—***************---------------—****************//
assumed the first half circle
***********----------------—***************---------------—****************//
***********----------------—***************---------------—****************//
in the Local Helix case, phi must be small
Implements IDedxCorrecSvc.
Definition at line 1489 of file DedxCorrecSvc.cxx.
1489 {
1490
1499
1500 HepPoint3D East_origin(East_lay_X, East_lay_Y, East_lay_Z);
1501 HepPoint3D West_origin(West_lay_X, West_lay_Y, West_lay_Z);
1502 Hep3Vector wire = (CLHEP::Hep3Vector)East_origin - (CLHEP::Hep3Vector)West_origin;
1503 HepPoint3D piovt_z =(z*10+length/2 )/length * wire + West_origin;
1504 piovt_z = piovt_z*0.1; // conversion of the units(mm=>cm)
1505
1506
1507 //-------------------------------temp -------------------------------//
1509 //-------------------------------temp -------------------------------//
1510
1516
1517 // Choose the local field !!
1519 double ALPHA_loc = 1000/(2.99792458*Bz);
1520
1521 // Choose the local field !!
1522 //double Bz = 1.0; // will be obtained from magnetic field service
1523 //double ALPHA_loc = 1000/(2.99792458*Bz);
1524 //double ALPHA_loc = 333.564095;
1525 int charge = ( kappa >= 0 )? 1 : -1;
1526 double rho = ALPHA_loc/kappa;
1527 double pt = fabs( 1.0/kappa );
1528 double lambda = atan( tanl );
1529 double theta = M_PI_2- lambda;
1531 // theta = 2.0*M_PI*theta/360.;
1534 double snf0 = (1. - csf0) * (1. + csf0);
1535 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
1537 //if(ntpFlag>0)
1538 //cout<<"rho = "<<rho<<" hel.dr() + rho = "<<hel.dr() + rho<<endl;
1539
1540 //-------------------------------temp -------------------------------//
1544 HepPoint3D ccenter(x_c, y_c, 0.0);
1545 double m_c_perp(ccenter.perp());
1546 Hep3Vector m_c_unit((HepPoint3D)ccenter.unit());
1547 //-------------------------------temp -------------------------------//
1548
1549 ////change to boost coordinate system
1550 double x_c_boost = x_c - piovt_z.x();
1551 double y_c_boost = y_c - piovt_z.y();
1552 double z_c_boost = z_c - piovt_z.z();
1553 HepPoint3D ccenter_boost(x_c_boost, y_c_boost, 0.0);
1554 double m_c_perp_boost(ccenter_boost.perp());
1555 //if(ntpFlag>0)
1556 //cout<<" ccenter = "<<ccenter<<" m_c_perp ="<<m_c_perp<<endl;
1557 Hep3Vector m_c_unit_boost((HepPoint3D)ccenter_boost.unit());
1558
1559 //phi region estimation
1560 double phi_io[2];
1561 HepPoint3D IO = (-1)*charge*m_c_unit;
1564 //double dphi0_bak = IO_phi - phi;
1565 //if(dphi0 != 0)
1566 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1569 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl;
1570 //cout<<"charge is = "<<charge<<endl;
1571 phi_io[0] = -(1+charge)*M_PI_2 - charge*dphi0;
1572 //phi_io[0] = -(1+charge)*M_PI_2 + charge*dphi0;
1573 phi_io[1] = phi_io[0]+1.5*M_PI;
1574 //cout<<"phi_io[0] is : "<<phi_io[0]<<" phi_io[1]:"<<phi_io[1]<<endl;
1575 double m_crio[2];
1576 double m_zb, m_zf, Calpha;
1577
1578 // retrieve Mdc geometry information
1586 //double length = geosvc->Layer(layer)->Length();
1588
1589 ///***********-------------------***************------------------****************//
1590 //check the value if same //
1591 ///***********-------------------***************------------------****************//
1593 int wid = w0id + cellid;
1600 double r_lay_backward = sqrt(x_lay_backward*x_lay_backward+y_lay_backward*y_lay_backward);
1601 double r_lay_forward = sqrt(x_lay_forward*x_lay_forward+y_lay_forward*y_lay_forward);
1602 double r_lay_use = ((z*10+length/2)/length)*(r_lay_backward-r_lay_forward) + r_lay_forward;
1603 /*for(int i=0; i<43; i++){
1604 double r_up= geosvc->Layer(i)->RCSiz1();
1605 double r_down= geosvc->Layer(i)->RCSiz2();
1606 cout<<i<<" "<<r_up<<" "<<r_down<<endl;
1607 }*/
1608 // shift = shift*type;
1609 // cout<< "type "<< type <<endl;
1610 r_lay_use = 0.1*r_lay_use;
1611 rcsiz1 = 0.1*rcsiz1;
1612 rcsiz2 = 0.1*rcsiz2;
1613 rlay = 0.1*rlay;
1614 length = 0.1*length;
1615 m_zb = 0.5*length;
1616 m_zf = -0.5*length;
1617 m_crio[0] = rlay - rcsiz1;
1618 m_crio[1] = rlay + rcsiz2;
1619 //conversion of the units(mm=>cm)
1620 int sign = -1; ///assumed the first half circle
1621 int epflag[2];
1622 Hep3Vector iocand[2];
1623 Hep3Vector cell_IO[2];
1624 double dphi, downin;
1625 Hep3Vector zvector;
1626 //if (ntpFlag>0) cout<<"z = "<<z<<" , m_zb = "<<m_zb<<" , m_zf = "<<m_zf<<endl;
1627 if( type ) {
1628 downin = (z*z-m_zb*m_zb)*pow(tan(slant),2);
1629 m_crio[0] = sqrt(m_crio[0]*m_crio[0]+downin);
1630 m_crio[1] = sqrt(m_crio[1]*m_crio[1]+downin);
1631 }
1632
1633 //int stced[2];
1634
1635 for( int i = 0; i < 2; i++ ) {
1636 double cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[i]*m_crio[i] - rho*rho;
1637 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[i] );
1638 if(fabs(cos_alpha)>1&&i==0) return(-1.0);
1639 if(fabs(cos_alpha)>1&&i==1) {
1640 cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[0]*m_crio[0] - rho*rho;
1641 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[0] );
1642 Calpha = 2.0*M_PI-acos( cos_alpha );
1643 } else {
1644 Calpha = acos( cos_alpha );
1645 }
1646 epflag[i] = 0;
1647 iocand[i] = m_c_unit_boost;
1648 iocand[i].rotateZ( charge*sign*Calpha );
1649 iocand[i]*= m_crio[i];
1650 //if (ntpFlag>0) cout<<"iocand corridate is : "<<iocand[i]<<endl;
1651
1652 ///***********-------------------***************------------------****************//
1653 //compare with standard coordinate system results //
1654 ///***********-------------------***************------------------****************//
1655 //-------------------------------temp-------------------------//
1656 iocand[i] = iocand[i]+ piovt_z; //change to standard coordinate system
1657 //iocand[i].y() = iocand[i].y() + piovt_z.y();
1658 //if (ntpFlag>0) cout<<i<<setw(10)<<iocand[i].x()<<setw(10)<<iocand[i].y()<<endl;
1659 //------------------------------temp -------------------------//
1660
1661 double xx = iocand[i].x() - x_c;
1662 double yy = iocand[i].y() - y_c;
1663 //dphi = atan2(yy, xx) - phi0 - M_PI_2*(1+charge);
1664 dphi = atan2(yy, xx) - phi0 - M_PI_2*(1-charge);
1666
1667 if( dphi < phi_io[0] ) {
1668 dphi += 2*M_PI;
1669 }
1670 else if( phi_io[1] < dphi ) {
1671 dphi -= 2*M_PI;
1672 }
1673 ///in the Local Helix case, phi must be small
1674
1675 //Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl);
1676 Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl-piovt_z.z());
1677 //if (ntpFlag>0) cout<<"z_c-rho*dphi*tanl : "<<z_c-rho*dphi*tanl<<endl;
1678 cell_IO[i] = iocand[i];
1679 cell_IO[i] += zvector;
1680 //---------------------------------temp ---------------------------------//
1681 //cell_IO[i] = hel.x(dphi);//compare with above results
1682 //---------------------------------temp ---------------------------------//
1683
1684 double xcio, ycio, phip;
1685 ///// z region check active volume is between zb+2. and zf-2. [cm]
1686 /*
1687 float delrin = 2.0;
1688 if( m_zf-delrin > zvector.z() ){
1689 phip = z_c - m_zb + delrin;
1690 phip = phip/( rho*tanl );
1691 phip = phip + phi0;
1692 xcio = x_c - rho*cos( phip );
1693 ycio = y_c - rho*sin( phip );
1694 cell_IO[i].setX( xcio );
1695 cell_IO[i].setY( ycio );
1696 cell_IO[i].setZ( m_zb+delrin );
1697 epflag[i] = 1;
1698 }
1699 else if( m_zb+delrin < zvector.z() ){
1700 phip = z_c - m_zf-delrin;
1701 phip = phip/( rho*tanl );
1702 phip = phip + phi0;
1703 xcio = x_c - rho*cos( phip );
1704 ycio = y_c - rho*sin( phip );
1705 cell_IO[i].setX( xcio );
1706 cell_IO[i].setY( ycio );
1707 cell_IO[i].setZ( m_zf-delrin );
1708 epflag[i] = 1;
1709 }
1710 else{
1711 */
1712 // cell_IO[i] = iocand;
1713 // cell_IO[i] += zvector;
1714 // }
1715 //dphi = dphi -M_PI;
1716 cell_IO[i] = hel.x(dphi);
1717 //if (ntpFlag>0) { cout<<"cell_IO corridate : "<<cell_IO[i]<<endl;}
1718 // if(i==0) cout<<"zhit value is = "<<z<<endl;
1719 }
1720
1721 // path length estimation
1722 // phi calculation
1723 Hep3Vector cl = cell_IO[1] - cell_IO[0];
1724
1725 //float ch_tphi, ch_tdphi;
1726 double ch_theta;
1727 double ch_dphi;
1728 double ch_ltrk = 0;
1729 double ch_ltrk_rp = 0;
1730 ch_dphi = cl.perp()*0.5/(ALPHA_loc*pt);
1731 ch_dphi = 2.0 * asin( ch_dphi );
1732 ch_ltrk_rp = ch_dphi*ALPHA_loc*pt;
1733 double rpi_path = sqrt(cl.x()*cl.x()+cl.y()*cl.y());
1734 ch_ltrk = sqrt( ch_ltrk_rp*ch_ltrk_rp + cl.z()*cl.z() );
1735 double path = ch_ltrk_rp/ sintheta;
1736 ch_theta = cl.theta();
1737 /* if (ntpFlag>0)
1738 {
1739 // if((ch_ltrk_rp-rpi_path)>0.001 || (ch_ltrk-path)>0.001)
1740 cout<<"ch_ltrk_rp : "<<ch_ltrk_rp<<" rpi_path: "<<rpi_path<<endl;
1741 cout<<"ch_ltrk : "<<ch_ltrk<<" path:"<<path<<endl;
1742 }
1743 */
1744 /*
1745 if( ch_theta < theta*0.85 || 1.15*theta < ch_theta)
1746 ch_ltrk *= -1; /// miss calculation
1747
1748 if( epflag[0] == 1 || epflag [1] == 1 )
1749 ch_ltrk *= -1; /// invalid region for dE/dx or outside of Mdc
1750 */
1751 // judge how many cells are traversed and deal with different situations
1752 double phibin, phi_in, phi_out, inlow, inup, outlow, outup, gap, phi1, phi2, phi_mid, phi_midin, phi_midout;
1753 int cid_in, cid_out;
1754 double inlow_z, inup_z, outlow_z, outup_z, gap_z, phi1_z, phi2_z, phi_mid_z, phi_midin_z, phi_midout_z;
1755 //cache sampl in each cell for this layer
1756
1757 std::vector<double> sampl;
1758 sampl.resize(ncell);
1759 for(int k=0; k<ncell; k++) {
1760 sampl[k] = -1.0;
1761 }
1762
1763 cid_in = cid_out = -1;
1764 phi_in = cell_IO[0].phi();
1765 phi_out = cell_IO[1].phi();
1766 //phi_in = iocand[0].phi();
1767 //phi_out = iocand[1].phi();
1770 phibin = 2.0*M_PI/ncell;
1771 // gap = fabs(phi_out-phi_in);
1772
1773 //determine the in/out cell id
1774 Hep3Vector cell_mid=0.5*(cell_IO[0]+cell_IO[1]);
1775 //Hep3Vector cell_mid=0.5*(iocand[0]+iocand[0]);
1776 //cout<<cell_mid<<endl;
1777 //double stereophi = shift*phibin*(0.5-cell_mid.z()/length);
1778 //phioffset = phioffset+stereophi;
1779 double stphi[2], phioff[2];
1780 stphi[0] = shift*phibin*(0.5-cell_IO[0].z()/length);
1781 stphi[1] = shift*phibin*(0.5-cell_IO[1].z()/length);
1782 //stphi[0] = shift*phibin*(0.5-z/length);
1783 //stphi[1] = shift*phibin*(0.5-z/length);
1784 phioff[0] = phioffset+stphi[0];
1785 phioff[1] = phioffset+stphi[1];
1786
1787 for(int i=0; i<ncell; i++) {
1788
1793 //if(ntpFlag>0) cout<<layer<<setw(10)<<i<<setw(10)<<x_lay_forward<<setw(10)<<y_lay_forward<<setw(10)<<x_lay_backward<<setw(10)<<y_lay_backward<<setw(10)<<r_lay_forward<<setw(10)<<endl;
1794 //Hep3Vector lay_backward, lay_forward;
1795 Hep3Vector lay_backward(x_lay_backward_cell, y_lay_backward_cell, 0);
1796 Hep3Vector lay_forward(x_lay_forward_cell, y_lay_forward_cell, 0);
1797 Hep3Vector Cell_z[2];
1798 Cell_z[0] = ((cell_IO[0].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1799 Cell_z[1] = ((cell_IO[1].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward;
1800 double z_phi[2];
1801 z_phi[0] = Cell_z[0].phi();
1802 z_phi[1] = Cell_z[1].phi();
1805 //double backward_phi = lay_backward.phi();
1806 //double forward_phi = lay_forward.phi();
1807 //backward_phi = fmod( backward_phi+4*M_PI,2*M_PI );
1808 //forward_phi = fmod( forward_phi+4*M_PI,2*M_PI );
1809 //if(ntpFlag>0) cout<<"backward_phi cal : "<<atan2(y_lay_backward,x_lay_backward)<<" forward_phi : "<<atan2(y_lay_forward,x_lay_forward)<<endl;
1810 //if(ntpFlag>0) cout<<layer<<" backward_phi : "<<backward_phi<<" forward_phi : "<<forward_phi<<endl;
1811
1812 //double z_phi[2];
1813 //z_phi[0] = ((cell_IO[0].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1814 //z_phi[1] = ((cell_IO[1].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi;
1815 //if(ntpFlag>0) cout<<"phi z : "<<z_phi[0]<<" "<<z_phi[1]<<endl;
1816 inlow_z = z_phi[0] - phibin*0.5;
1817 inup_z = z_phi[0] + phibin*0.5;
1818 outlow_z = z_phi[1] - phibin*0.5;
1819 outup_z = z_phi[1] + phibin*0.5;
1824
1825
1826 // for stereo layer
1827 inlow = phioff[0]+phibin*i-phibin*0.5;
1828 inup = phioff[0]+phibin*i+phibin*0.5;
1829 outlow = phioff[1]+phibin*i-phibin*0.5;
1830 outup = phioff[1]+phibin*i+phibin*0.5;
1835
1836#ifdef YDEBUG
1837 if(ntpFlag > 0) cout << "shift " << shift
1838 <<" phi_in " << phi_in << " phi_out " << phi_out
1839 << " inup "<< inup << " inlow " << inlow
1840 << " outup "<< outup << " outlow " << outlow << endl;
1841#endif
1842
1843 /*if(phi_in>=inlow&&phi_in<inup) cid_in = i;
1844 if(phi_out>=outlow&&phi_out<outup) cid_out = i;
1845 if ( inlow>inup) {
1846 if((phi_in>=inlow&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup)) cid_in = i;
1847 }
1848 if ( outlow>outup) {
1849 if((phi_out>=outlow&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup)) cid_out = i;
1850 }*/
1851
1852 if(phi_in>=inlow_z&&phi_in<inup_z) cid_in = i;
1853 if(phi_out>=outlow_z&&phi_out<outup_z) cid_out = i;
1854 if ( inlow_z>inup_z) {
1856 }
1857 if ( outlow_z>outup_z) {
1859 }
1860 }
1861
1863 gap = -999.0;
1864 //only one cell traversed
1865 //if(ntpFlag > 0) cout<<"cid_in = "<<cid_in<<" cid_out = "<<cid_out<<endl;
1866 if(cid_in == -1 || cid_out == -1) return -1;
1867
1868 if( cid_in == cid_out) {
1869 sampl[cid_in]= ch_ltrk;
1870 //return ch_ltrk;
1871 return sampl[cellid];
1872 } else if(cid_in < cid_out) {
1873 //in cell id less than out cell id
1874 //deal with the special case crossing x axis
1875 if( cid_out-cid_in>ncell/2 ) {
1876
1877 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1882 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1883 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1884 Hep3Vector Cell_z[2];
1885 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1886 double phi_cin_z = Cell_z[0].phi();
1892 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1893 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1894 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1895 double phi_cout_z = Cell_z[1].phi();
1897
1898 phi_midin_z = phi_cin_z-phibin*0.5;
1899 phi_midout_z = phi_cout_z+phibin*0.5;
1902 phi1_z = phi_midout_z-phi_out;
1903 phi2_z = phi_in-phi_midin_z;
1906 gap_z = phi1_z+phi2_z+(ncell-1-cid_out+cid_in)*phibin;
1908 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
1909 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
1910 for( int jj = cid_out+1; jj<ncell; jj++) {
1911 sampl[jj]=phibin/gap_z*ch_ltrk;
1912 }
1913 for( int jj = 0; jj<cid_in; jj++) {
1914 sampl[jj]=phibin/gap_z*ch_ltrk;
1915 }
1916
1917 /*
1918 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
1919 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
1920 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
1921 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1922 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1923 phi1 = phi_midout-phi_out;
1924 phi2 = phi_in-phi_midin;
1925 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1926 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1927 gap = phi1+phi2+(ncell-1-cid_out+cid_in)*phibin;
1928 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1929 sampl[cid_in]=phi2/gap*ch_ltrk;
1930 sampl[cid_out]=phi1/gap*ch_ltrk;
1931 for( int jj = cid_out+1; jj<ncell; jj++) {
1932 sampl[jj]=phibin/gap*ch_ltrk;
1933 }
1934 for( int jj = 0; jj<cid_in; jj++) {
1935 sampl[jj]=phibin/gap*ch_ltrk;
1936 }*/
1937 /*
1938 cout<<"LLLLLLLLLLLLL"<<endl;
1939 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
1940 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
1941 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
1942 << phi_out << " phi_midout " << phi_midout <<endl;
1943 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
1944 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
1945 */
1946 } else {
1947 //normal case
1952 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
1953 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
1954 Hep3Vector Cell_z[2];
1955 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
1956 double phi_cin_z = Cell_z[0].phi();
1962 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
1963 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
1964 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
1965 double phi_cout_z = Cell_z[1].phi();
1967
1968 phi_midin_z = phi_cin_z+phibin*0.5;
1969 phi_midout_z = phi_cout_z-phibin*0.5;
1972 phi1_z = phi_midin_z-phi_in;
1973 phi2_z = phi_out-phi_midout_z;
1976 gap_z = phi1_z+phi2_z+(cid_out-cid_in-1)*phibin;
1978 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
1979 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
1980 for( int jj = cid_in+1; jj<cid_out; jj++) {
1981 sampl[jj]=phibin/gap_z*ch_ltrk;
1982 }
1983 //normal case
1984 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
1985 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
1986 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
1987 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
1988 phi1 = phi_midin-phi_in;
1989 phi2 = phi_out-phi_midout;
1990 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
1991 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
1992 gap = phi1+phi2+(cid_out-cid_in-1)*phibin;
1993 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
1994 sampl[cid_in]=phi1/gap*ch_ltrk;
1995 sampl[cid_out]=phi2/gap*ch_ltrk;
1996 for( int jj = cid_in+1; jj<cid_out; jj++) {
1997 sampl[jj]=phibin/gap*ch_ltrk;
1998 }*/
1999 }
2000
2001 } else if(cid_in > cid_out) {
2002 //in cell id greater than out cell id
2003 //deal with the special case crossing x axis
2004 if( cid_in-cid_out>ncell/2 ) {
2009 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2010 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2011 Hep3Vector Cell_z[2];
2012 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2013 double phi_cin_z = Cell_z[0].phi();
2019 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2020 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2021 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2022 double phi_cout_z = Cell_z[1].phi();
2024
2025 phi_midin_z = phi_cin_z+phibin*0.5;
2026 phi_midout_z = phi_cout_z-phibin*0.5;
2029 phi1_z = phi_midin_z-phi_in;
2030 phi2_z = phi_out-phi_midout_z;
2033 gap_z = phi1_z+phi2_z+(ncell-1-cid_in+cid_out)*phibin;
2035 sampl[cid_out]=phi2_z/gap_z*ch_ltrk;
2036 sampl[cid_in]=phi1_z/gap_z*ch_ltrk;
2037 for( int jj = cid_in+1; jj<ncell; jj++) {
2038 sampl[jj]=phibin/gap_z*ch_ltrk;
2039 }
2040 for( int jj = 0; jj<cid_out; jj++) {
2041 sampl[jj]=phibin/gap_z*ch_ltrk;
2042 }
2043
2044 // if(gap>=M_PI) gap = 2.0*M_PI-gap;
2045 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5;
2046 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5;
2047 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2048 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2049 phi1 = phi_midin-phi_in;
2050 phi2 = phi_out-phi_midout;
2051 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2052 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2053 gap = phi1+phi2+(ncell-1-cid_in+cid_out)*phibin;
2054 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2055 sampl[cid_out]=phi2/gap*ch_ltrk;
2056 sampl[cid_in]=phi1/gap*ch_ltrk;
2057 for( int jj = cid_in+1; jj<ncell; jj++) {
2058 sampl[jj]=phibin/gap*ch_ltrk;
2059 }
2060 for( int jj = 0; jj<cid_out; jj++) {
2061 sampl[jj]=phibin/gap*ch_ltrk;
2062 }*/
2063 /*
2064 cout<<"LLLLLLLLLLLLL"<<endl;
2065 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2066 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2067 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2068 << phi_out << " phi_midout " << phi_midout <<endl;
2069 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2070 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl;
2071 */
2072 } else {
2077 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0);
2078 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0);
2079 Hep3Vector Cell_z[2];
2080 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin;
2081 double phi_cin_z = Cell_z[0].phi();
2087 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0);
2088 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0);
2089 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout;
2090 double phi_cout_z = Cell_z[1].phi();
2092
2093 phi_midin_z = phi_cin_z-phibin*0.5;
2094 phi_midout_z = phi_cout_z+phibin*0.5;
2097 phi1_z = phi_midout_z-phi_out;
2098 phi2_z = phi_in-phi_midin_z;
2101 gap_z = phi1_z+phi2_z+(cid_in-cid_out-1)*phibin;
2103 sampl[cid_in]=phi2_z/gap_z*ch_ltrk;
2104 sampl[cid_out]=phi1_z/gap_z*ch_ltrk;
2105 for( int jj = cid_out+1; jj<cid_in; jj++) {
2106 sampl[jj]=phibin/gap_z*ch_ltrk;
2107 }
2108
2109 //normal case
2110 /*phi_midin = phioff[0]+phibin*cid_in-phibin*0.5;
2111 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5;
2112 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI );
2113 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI );
2114 phi1 = phi_midout-phi_out;
2115 phi2 = phi_in-phi_midin;
2116 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI);
2117 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI);
2118 gap = phi1+phi2+(cid_in-cid_out-1)*phibin;
2119 gap = fmod(gap+2.0*M_PI,2.0*M_PI);
2120 sampl[cid_in]=phi2/gap*ch_ltrk;
2121 sampl[cid_out]=phi1/gap*ch_ltrk;
2122 for( int jj = cid_out+1; jj<cid_in; jj++) {
2123 sampl[jj]=phibin/gap*ch_ltrk;
2124 }*/
2125 }
2126 }
2127
2128#ifdef YDEBUG
2129 if(sampl[cellid]<0.0) {
2130 if(cid_in!=cid_out) cout<<"?????????"<<endl;
2131 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift
2132 << " cid_in " << cid_in << " cid_out " << cid_out << endl;
2133
2134 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out "
2135 << phi_out << " phi_midout " << phi_midout <<endl;
2136 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 "
2138
2139
2140 for(int l=0; l<ncell; l++) {
2141 if(sampl[l]>=0.0)
2142 cout<<"picked cellid " << l << " sampling length "<< sampl[l]<< endl;
2143 }
2144 }
2145#endif
2146 return sampl[cellid];
2147}
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Dedx_Helix.cxx:209
virtual double getReferField()=0
virtual const MdcGeoLayer *const Layer(unsigned id)=0
virtual const MdcGeoWire *const Wire(unsigned id)=0
virtual const MdcGeoGeneral *const GeneralLayer(unsigned id)=0
◆ queryInterface()
|
virtual |
Definition at line 54 of file DedxCorrecSvc.cxx.
54 {
55 //cout<<"DedxCorrecSvc::queryInterface"<<endl;
56 if( IID_IDedxCorrecSvc.versionMatch(riid) ){
58 } else{
59 return Service::queryInterface(riid, ppvInterface);
60 }
61 return StatusCode::SUCCESS;
62}
Definition: IDedxCorrecSvc.h:13
◆ RungCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 672 of file DedxCorrecSvc.cxx.
672 {
673 // cout<<"DedxCorrecSvc::RungCorrec"<<endl;
674 double dedx_rung;
675 int run_ture =0;
676 //cout<<"N_run : "<<N_run<<" runNO : "<<runNO<<endl;
677
678 for(int j=0;j<10000;j++) {
679 if(runNO == m_rung[2][j]) {
680 dedx_rung = dedx/m_rung[0][j];
681 run_ture =1;
682 return dedx_rung;
683 }
684 }
685 if(run_ture ==0)
686 {
687 cout<<"Warning!!! in DedxCorrecSvc::RungCorrec(): no rungain to "<<runNO<<endl;
688 exit(1);
689 }
690}
Referenced by StandardCorrec(), StandardHitCorrec(), and StandardTrackCorrec().
◆ SaturCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 545 of file DedxCorrecSvc.cxx.
545 {
546 //cout<<"DedxCorrecSvc::SaturCorrec"<<endl;
547 double dedx_ggs;
548 //cout<<"costheta vaule is = "<<costheta<<endl;
549 costheta = fabs(costheta);
550 if(m_par_flag == 1) {
551 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*costheta +
552 m_ggs[2][layer]*pow(costheta,2) + m_ggs[3][layer]*pow(costheta,3);
553 } else if(m_par_flag == 0) {
554 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*T1(costheta) +
555 m_ggs[2][layer]*T2(costheta) + m_ggs[3][layer]*T3(costheta);
556 }
557 //cout<<"dedx_ggs is = "<<dedx_ggs<<endl;
558 dedx_ggs = dedx/dedx_ggs;
559 if(dedx_ggs>0.0) return dedx_ggs;
560 else return dedx;
561}
Referenced by CellCorrec(), LayerCorrec(), StandardCorrec(), and StandardHitCorrec().
◆ set_flag()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 1131 of file DedxCorrecSvc.cxx.
1131 {
1132 // cout<<"DedxCorrecSvc::set_flag"<<endl;
1133 cout<<"calib_F is = "<<calib_F<<endl;
1135 else m_enta_flag = 1;
1136 m_rung_flag = ( calib_F & 1 )? 1 : 0;
1137 m_wireg_flag = ( calib_F & 2 )? 1 : 0;
1138 m_dcasin_flag = ( calib_F & 4 )? 1 : 0;
1139 m_ggs_flag = ( calib_F & 8 )? 1 : 0;
1140 m_ddg_flag = 0;
1141 //m_ddg_flag = ( calib_F & 8 )? 1 : 0;
1142 m_t0_flag = ( calib_F & 16 )? 1 : 0;
1143 m_sat_flag = ( calib_F & 32 )? 1 : 0;
1144 m_layerg_flag = ( calib_F & 64 )? 1 : 0;
1145 //m_dcasin_flag = ( calib_F & 128 )? 1 : 0;
1146 m_dip_flag = ( calib_F & 256 )? 1 : 0;
1147 m_mdcg_flag = ( calib_F & 512 )? 1 : 0;
1148}
◆ StandardCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 919 of file DedxCorrecSvc.cxx.
919 {
920 // cout<<"DedxCorrecSvc::StandardCorrec"<<endl;
921 //int layid = MdcID::layer(mdcid);
922 //int localwid = MdcID::wire(mdcid);
923 //int w0id = geosvc->Layer(layid)->Wirst();
924 //int wid = w0id + localwid;
925 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
926 ////pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
927 double ex = adc;
928 if( runNO<0 ) return ex;
929 ex = ex*1.5/pathl;
930 //double ex = adc/pathl;
931 //if( runNO<0 ) return ex;
932 if( ntpFlag ==0) return ex;
933 //double ex = adc/pathl;
934 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_dcasin_flag==0 && m_ddg_flag == 0
935 && m_layerg_flag == 0 && m_zdep_flag == 0 &&
936 m_ggs_flag == 0) return ex;
937
938 if(m_rung_flag) {
939 ex = RungCorrec( runNO, ex );
940 }
941
942 if( m_wireg_flag) {
943 ex = WireGainCorrec(wid, ex);
944 }
945
946 if( m_dcasin_flag) {
947 if(runFlag<3) {
948 ex = DriftDistCorrec( layid, dd, ex );
949 }
951 }
952
953 if(m_enta_flag && runFlag>=3){
954 ex = EntaCorrec(layid, eangle, ex);
955 }
956
957 // ddg is not use anymore, it's replaced by DocaSin
958 if(m_ddg_flag) {
959 ex = DriftDistCorrec( layid, dd, ex );
960 }
961
962 if(m_ggs_flag) {
963 if(runFlag <3 ){
964 ex = SaturCorrec( layid, costheta, ex );
965 }
967 // Staur is OLD, Costheta is NEW.
968 }
969
970 if( m_sat_flag) {
971 ex = HadronCorrec( costheta, ex );
972 }
973
974 if( m_zdep_flag) {
975 ex = ZdepCorrec( layid, z, ex );
976 }
977
978 if( m_layerg_flag) {
979 ex = LayerGainCorrec( layid, ex );
980 }
981
982 if( m_dip_flag){
983 ex = DipAngCorrec(layid, costheta, ex);
984 }
985
986 if( m_mdcg_flag) {
987 ex = GlobalCorrec( ex );
988 }
989 return ex;
990}
double CosthetaCorrec(double costheta, double ex) const
Definition: DedxCorrecSvc.cxx:565
double RungCorrec(int runNO, double ex) const
Definition: DedxCorrecSvc.cxx:672
double EntaCorrec(int layid, double enta, double ex) const
Definition: DedxCorrecSvc.cxx:715
double WireGainCorrec(int wireid, double ex) const
Definition: DedxCorrecSvc.cxx:533
double HadronCorrec(double costheta, double dedx) const
Definition: DedxCorrecSvc.cxx:651
double DocaSinCorrec(int layid, double doca, double enta, double ex) const
Definition: DedxCorrecSvc.cxx:785
double DipAngCorrec(int layid, double costheta, double ex) const
Definition: DedxCorrecSvc.cxx:828
◆ StandardHitCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 993 of file DedxCorrecSvc.cxx.
993 {
994 if(m_debug) cout<<"DedxCorrecSvc::StandardHitCorrec"<<endl;
995 //int layid = MdcID::layer(mdcid);
996 //int localwid = MdcID::wire(mdcid);
997 //int w0id = geosvc->Layer(layid)->Wirst();
998 //int wid = w0id + localwid;
999 //double pathl = PathL(ntpFlag, hel, layid, localwid, z);
1000 //cout<<"DedxCorrecSvc::StandardHitCorrec pathl= "<<pathl<<endl;
1001 //pathl= PathLCosmic(hel, layid, localwid, z, sigmaz );
1002 double ex = adc;
1003 if( runNO<0 ) return ex;
1004 ex = ex*1.5/pathl;
1005 //if( runNO<0 ) return ex;
1006 //double ex = adc/pathl;
1007 if( ntpFlag ==0) return ex;
1008 //if(ntpFlag>0) cout<<"dE/dx value after path correction: "<<ex<<endl;
1009 //double ex = adc/pathl;
1010 //cout<<m_rung_flag << " "<<m_wireg_flag<<" "<<m_ddg_flag<<" "<<m_ggs_flag<<endl;
1011 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0
1012 && m_layerg_flag == 0 && m_zdep_flag == 0 && m_dcasin_flag==0
1013 && m_ggs_flag == 0 && m_enta_flag==0) return ex;
1014
1015 if(m_rung_flag) {
1016 ex = RungCorrec( runNO, ex );
1017 }
1018 //if(ntpFlag>0) cout<<"dE/dx value after run by run correction: "<<ex<<endl;
1019
1020 if( m_wireg_flag) {
1021 ex = WireGainCorrec(wid, ex);
1022 }
1023 //if(ntpFlag>0) cout<<"dE/dx value after wire gain correction: "<<ex<<endl;
1024
1025 if( m_dcasin_flag) {
1026 if(runFlag<3){
1027 ex = DriftDistCorrec( layid, dd, ex );
1028 }
1029 else{
1030 //cout<<layid<<" "<<dd<<" "<<eangle<<" "<<ex<<endl;
1031 ex = DocaSinCorrec(layid, dd, eangle, ex);
1032 }
1033 }
1034
1035 // 1D entrance angle correction
1036 if(m_enta_flag && runFlag>=3){
1037 ex = EntaCorrec(layid, eangle, ex);
1038 }
1039
1040 // ddg is not used anymore
1041 if( m_ddg_flag) {
1042 ex = DriftDistCorrec( layid, dd, ex );
1043 }
1044 //if(ntpFlag>0) cout<<"dE/dx value after ddg correction: "<<ex<<endl;
1045
1046 if(m_ggs_flag) {
1047 if(runFlag <3 ){
1048 ex = SaturCorrec( layid, costheta, ex );
1049 }
1050 else{ ex = ex;} // do not do the cos(theta) correction at hit level
1051 }
1052 //if(ntpFlag>0) cout<<"dE/dx value after costheta correction: "<<ex<<endl;
1053 return ex;
1054}
◆ StandardTrackCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 1058 of file DedxCorrecSvc.cxx.
1058 {
1059 if(m_debug) cout<<"DedxCorrecSvc::StandardTrackCorrec"<<endl;
1060 if( runNO<0 ) return ex;
1061 if( ntpFlag ==0) return ex;
1062 if( runFlag <3) return ex;
1063 if( calib_rec_Flag ==1){
1064 ex = CosthetaCorrec( costheta, ex );
1066 return ex;
1067 }
1068
1069 //if(ntpFlag>0) cout<<"trcak value before costheta correction: "<<ex<<endl;
1070 if( m_ggs_flag) {
1071 ex = CosthetaCorrec( costheta, ex );
1072 }
1073 //if(ntpFlag>0) cout<<"trcak value after costheta correction: "<<ex<<endl;
1074 if(m_t0_flag){
1076 else if(runFlag>3) {ex=ex;}
1077 }
1078 //if(ntpFlag>0) cout<<"trcak value after all correction: "<<ex<<endl;
1079 if( m_sat_flag) {
1080 ex = HadronCorrec( costheta, ex );
1081 }
1082
1083 if(m_rung_flag && calib_rec_Flag ==2){
1085 ex = ex/rungain_temp;
1086 }
1087
1088 //if(ntpFlag>0) cout<<"trcak value no run gain correction: "<<ex<<endl;
1089 return ex;
1090
1091}
◆ T0Correc()
| double DedxCorrecSvc::T0Correc | ( | double | t0, |
| double | ex | ||
| ) | const |
Definition at line 619 of file DedxCorrecSvc.cxx.
619 {
620 // cout<<"DedxCorrecSvc::T0Correc"<<endl;
621 double dedx_t0;
622 const int nbin=t0_k.size()+1;
623 if(nbin!=35)
624 cout<<"there should be 35 bins for t0 correction, check it!"<<endl;
625
627 if(t0<575)
628 n=0;
629 else if(t0<610)
630 n=1;
631 else if(t0>=610 && t0<1190){
632 n=(int)( (t0-610.0)/20.0 ) + 2;
633 }
634 else if(t0>=1190 && t0<1225)
635 n=31;
636 else if(t0>=1225 && t0<1275)
637 n=32;
638 else if(t0>=1275)
639 n=33;
640
642
643 if(dedx_t0>0){
644 dedx_t0 = dedx/dedx_t0;
645 return dedx_t0;
646 }
647 else return dedx;
648}
Referenced by StandardTrackCorrec().
◆ TrkCorrec()
|
virtual |
dEdx calib. for each track
Implements IDedxCorrecSvc.
Definition at line 904 of file DedxCorrecSvc.cxx.
904 {
905 // cout<<"DedxCorrecSvc::TrkCorrec"<<endl;
906 if( m_mdcg_flag == 0 ) return dedx;
907 ///dEdx calib. for each track
908 double calib_ex = dedx;
909
910 double run_const = m_dedx_gain;
911 if( run_const > 0 && m_mdcg_flag != 0 ) calib_ex /= run_const;
912
913 return calib_ex;
914
915}
◆ WireGainCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 533 of file DedxCorrecSvc.cxx.
533 {
534 if( m_wire_g[wireid] > 0 ){
535 double ch_dedx = (ex/m_wire_g[wireid])*m_valid[wireid];
536 return ch_dedx;
537 }
538 else if( m_wire_g[wireid] == 0 ){
539 return ex;
540 }
541 else return 0;
542}
Referenced by StandardCorrec(), and StandardHitCorrec().
◆ ZdepCorrec()
|
virtual |
Implements IDedxCorrecSvc.
Definition at line 766 of file DedxCorrecSvc.cxx.
766 {
767 // cout<<"DedxCorrecSvc::ZdepCorrec"<<endl;
768 double dedx_zdep;
769 if(m_par_flag == 1) {
770 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*z +
771 m_zdep[2][layer]*z*z + m_zdep[3][layer]*pow(z,3);
772 } else if(m_par_flag == 0) {
773 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*T1(z) +
774 m_zdep[2][layer]*T2(z) + m_zdep[3][layer]*T3(z);
775 }
776 dedx_zdep = dedx/dedx_zdep;
777 if (dedx_zdep>0.0) return dedx_zdep;
778 else return dedx;
779
780 //return dedx;
781}
Referenced by CellCorrec(), LayerCorrec(), and StandardCorrec().
The documentation for this class was generated from the following files:
