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) |
| Public Member Functions inherited from IDedxCorrecSvc | |
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;
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);
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;
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];
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
double SzWest(void) const
double SxEast(void) const
double SyWest(void) const
double SxWest(void) const
double SzEast(void) const
double SyEast(void) const
double Radius(void) const
double Slant(void) const
double RCSiz2(void) const
double Length(void) const
MdcGeoSuper * Sup(void) const
double Shift(void) const
double RCSiz1(void) const
int NCell(void) const
int Wirst(void) const
double Offset(void) const
HepPoint3D FWirePos(void) const
HepPoint3D BWirePos(void) const
HepPoint3D Forward(void) const
HepPoint3D Backward(void) const
◆ 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:
