#include <DedxCorrecSvc.h>
Inheritance diagram for DedxCorrecSvc:
Public Member Functions | |
double | CellCorrec (int ser, double adc, double dd, double enta, double z, double costheta) const |
double | CellCorrec (int ser, double adc, double dd, double enta, double z, double costheta) const |
double | CosthetaCorrec (double costheta, double ex) const |
double | CosthetaCorrec (double costheta, double ex) const |
double | D2I (const double &cosTheta, const double &D) const |
double | D2I (const double &cosTheta, const double &D) const |
DedxCorrecSvc (const std::string &name, ISvcLocator *svcloc) | |
DedxCorrecSvc (const std::string &name, ISvcLocator *svcloc) | |
double | DipAngCorrec (int layid, double costheta, double ex) const |
double | DipAngCorrec (int layid, double costheta, double ex) const |
double | DocaSinCorrec (int layid, double doca, double enta, double ex) const |
double | DocaSinCorrec (int layid, double doca, double enta, double ex) const |
double | DriftDistCorrec (int layid, double ddrift, double ex) const |
double | DriftDistCorrec (int layid, double ddrift, double ex) const |
double | EntaCorrec (int layid, double enta, double ex) const |
double | EntaCorrec (int layid, double enta, double ex) const |
virtual StatusCode | finalize () |
virtual StatusCode | finalize () |
double | GlobalCorrec (double dedx) const |
double | GlobalCorrec (double dedx) const |
double | HadronCorrec (double costheta, double dedx) const |
double | HadronCorrec (double costheta, double dedx) const |
void | handle (const Incident &) |
void | handle (const Incident &) |
double | I2D (const double &cosTheta, const double &I) const |
double | I2D (const double &cosTheta, const double &I) const |
virtual StatusCode | initialize () |
virtual StatusCode | initialize () |
double | LayerCorrec (int layer, double z, double costheta, double ex) const |
double | LayerCorrec (int layer, double z, double costheta, double ex) const |
double | LayerGainCorrec (int layid, double ex) const |
double | LayerGainCorrec (int layid, double ex) const |
double | PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const |
double | PathL (int ntpFlag, const Dedx_Helix &hel, int layer, int cellid, double z) const |
virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvUnknown) |
virtual StatusCode | queryInterface (const InterfaceID &riid, void **ppvUnknown) |
double | RungCorrec (int runNO, double ex) const |
double | RungCorrec (int runNO, double ex) const |
double | SaturCorrec (int layid, double costheta, double ex) const |
double | SaturCorrec (int layid, double costheta, double ex) const |
void | set_flag (int calib_F) |
void | set_flag (int calib_F) |
double | StandardCorrec (int typFlag, int ntpFlag, int runNO, Dedx_Helix &hel, Identifier mdcid, double adc, double dd, double enta, double z, double costheta) const |
double | StandardCorrec (int typFlag, int ntpFlag, int runNO, Dedx_Helix &hel, Identifier mdcid, double adc, double dd, double enta, double z, double costheta) const |
double | StandardHitCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, Dedx_Helix &hel, Identifier mdcid, double adc, double dd, double enta, double z, double costheta) const |
double | StandardHitCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, Dedx_Helix &hel, Identifier mdcid, double adc, double dd, double enta, double z, double costheta) const |
double | StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const |
double | StandardTrackCorrec (int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, double ex, double costheta, double t0) const |
double | T0Correc (double t0, double ex) const |
double | T0Correc (double t0, double ex) const |
double | TrkCorrec (double costheta, double dedx) const |
double | TrkCorrec (double costheta, double dedx) const |
double | WireGainCorrec (int wireid, double ex) const |
double | WireGainCorrec (int wireid, double ex) const |
double | ZdepCorrec (int layid, double zhit, double ex) const |
double | ZdepCorrec (int layid, double zhit, double ex) const |
~DedxCorrecSvc () | |
~DedxCorrecSvc () | |
Static Public Member Functions | |
const InterfaceID & | interfaceID () |
const InterfaceID & | interfaceID () |
Private Member Functions | |
double | cub (double &x) const |
double | cub (double &x) const |
void | init_param () |
void | init_param () |
void | init_param_svc () |
void | init_param_svc () |
double | sq (double &x) const |
double | sq (double &x) const |
double | sq10 (double &x) const |
double | sq10 (double &x) const |
double | sq11 (double &x) const |
double | sq11 (double &x) const |
double | sq4 (double &x) const |
double | sq4 (double &x) const |
double | sq5 (double &x) const |
double | sq5 (double &x) const |
double | sq6 (double &x) const |
double | sq6 (double &x) const |
double | sq7 (double &x) const |
double | sq7 (double &x) const |
double | sq8 (double &x) const |
double | sq8 (double &x) const |
double | sq9 (double &x) const |
double | sq9 (double &x) const |
double | T0 () |
double | T0 () |
double | T1 (double &x) const |
double | T1 (double &x) const |
double | T10 (double &x) const |
double | T10 (double &x) const |
double | T11 (double &x) const |
double | T11 (double &x) const |
double | T12 (double &x) const |
double | T12 (double &x) const |
double | T13 (double &x) const |
double | T13 (double &x) const |
double | T14 (double &x) const |
double | T14 (double &x) const |
double | T15 (double &x) const |
double | T15 (double &x) const |
double | T16 (double &x) const |
double | T16 (double &x) const |
double | T17 (double &x) const |
double | T17 (double &x) const |
double | T18 (double &x) const |
double | T18 (double &x) const |
double | T19 (double &x) const |
double | T19 (double &x) const |
double | T2 (double &x) const |
double | T2 (double &x) const |
double | T20 (double &x) const |
double | T20 (double &x) const |
double | T3 (double &x) const |
double | T3 (double &x) const |
double | T4 (double &x) const |
double | T4 (double &x) const |
double | T5 (double &x) const |
double | T5 (double &x) const |
double | T6 (double &x) const |
double | T6 (double &x) const |
double | T7 (double &x) const |
double | T7 (double &x) const |
double | T8 (double &x) const |
double | T8 (double &x) const |
double | T9 (double &x) const |
double | T9 (double &x) const |
Private Attributes | |
vector< double > | cos_b |
vector< double > | cos_b |
vector< double > | cos_k |
vector< double > | cos_k |
double | curve_par [5] |
IMdcGeomSvc * | geosvc |
mdc geometry | |
IMdcGeomSvc * | geosvc |
mdc geometry | |
double | Iner_Stepdoca |
double | m_alpha |
int | m_calib_flag |
int | m_dcasin_flag |
correction : doca and sinenta correction ( calib_F & 128 )? ON:OFF | |
double | m_ddg [4][43] |
int | m_ddg_flag |
calibration flag : drift distance ( calib_F & 8 ) ? ON : OFF | |
double | m_dedx_gain |
double | m_dedx_resol |
double | m_delta |
int | m_dip_flag |
correction : Qsaturation correction ( calib_F & 256 )? ON:OFF | |
double | m_docaeangle [40][40] |
double | m_enta [4][43] |
int | m_enta_flag |
calibration flag : entrance angle depsndence ( calib_F & 16 )? ON : OFF | |
double | m_gamma |
double | m_ggs [4][43] |
int | m_ggs_flag |
calibration flag : electron dip angle correction ( calib_F & 4 )? ON : OFF | |
int | m_init |
control flags | |
bool | m_initConstFlg |
double | m_layer_g [43] |
int | m_layerg_flag |
correction : mdc layer gain ( calib_F & 64 )? ON:OFF | |
int | m_mdcg_flag |
correction : mdc gain (run by run) ( calib_F & 512 )? ON:OFF | |
int | m_par_flag |
IMagneticFieldSvc * | m_pIMF |
IMagneticFieldSvc * | m_pIMF |
double | m_power |
double | m_ratio |
IntegerProperty | m_run |
double | m_rung [4][3000] |
int | m_rung_flag |
calibration flag : run by run correction ( calib_F & 1 )? ON : OFF | |
int | m_sat_flag |
calibration flag : hadron saturation correction ( calib_F & 32 ) ? ON : OFF | |
int | m_t0_flag |
calibration flag : t0 correction ( calib_F & 16 ) ? ON : OFF | |
double | m_valid [6796] |
vector< double > | m_venangle |
vector< double > | m_venangle |
double | m_wire_g [6796] |
int | m_wireg_flag |
calibration flag : wire by wire correction ( calib_F & 2 )? ON : OFF | |
double | m_zdep [4][43] |
int | m_zdep_flag |
calibration flag : z depsndence ( calib_F & 32 )? ON : OFF | |
int | N_run |
double | Out_Stepdoca |
double | sigma_par [6] |
vector< double > | t0_b |
vector< double > | t0_b |
vector< double > | t0_k |
vector< double > | t0_k |
|
00039 : 00040 geosvc(0),Service (name, svcloc) { 00041 declareProperty("Run",m_run=1); 00042 declareProperty("init",m_init=1); 00043 declareProperty("par_flag",m_par_flag=0); 00044 m_initConstFlg = false; 00045 // cout<<"DedxCorrecSvc::DedxCorrecSvc"<<endl; 00046 }
|
|
00048 { 00049 }
|
|
|
|
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00769 { 00770 00771 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0 00772 && m_layerg_flag == 0 && m_zdep_flag == 0 && 00773 m_ggs_flag == 0) return adc; 00774 adc = m_valid[ser]*m_wire_g[ser]*adc; 00775 // int lyr = Wire2Lyr(ser); 00776 int lyr = ser; 00777 double ex = adc; 00778 double correct1_ex , correct2_ex, correct3_ex, correct4_ex, correct5_ex; 00779 00780 if( m_ggs_flag) { 00781 correct1_ex = SaturCorrec( lyr, costheta, adc ); 00782 ex = correct1_ex; 00783 } else { 00784 correct1_ex = adc; 00785 } 00786 00787 if( m_ddg_flag) { 00788 correct2_ex = DriftDistCorrec( lyr, dd, correct1_ex ); 00789 ex = correct2_ex; 00790 } else { 00791 correct2_ex =correct1_ex; 00792 } 00793 if( m_enta_flag) { 00794 correct3_ex = DriftDistCorrec( lyr, sinenta, correct2_ex ); 00795 ex = correct3_ex; 00796 } else { 00797 correct3_ex =correct2_ex; 00798 } 00799 00800 if( m_zdep_flag) { 00801 correct4_ex = ZdepCorrec( lyr, z, correct3_ex ); 00802 ex = correct4_ex; 00803 } else { 00804 correct4_ex = correct3_ex; 00805 } 00806 00807 if( m_layerg_flag) { 00808 correct5_ex = LayerGainCorrec( lyr, correct4_ex ); 00809 ex = correct5_ex; 00810 } else { 00811 correct5_ex = correct4_ex; 00812 } 00813 return ex; 00814 00815 }
|
|
|
|
00507 { 00508 // cout<<"DedxCorrecSvc::CosthetaCorrec"<<endl; 00509 double dedx_cos; 00510 //cout<<"costheta vaule is = "<<costheta<<endl; 00511 if(fabs(costheta)>1.0) 00512 return dedx; 00513 00514 const int nbin=cos_k.size()+1; 00515 const double step=2.00/nbin; 00516 int n= (int)((costheta+1.00-0.5*step)/step); 00517 if(costheta>(1.00-0.5*step)) 00518 n=nbin-2; 00519 00520 if(costheta>-0.5*step && costheta<=0) 00521 dedx_cos=cos_k[n-1]*costheta+cos_b[n-1]; 00522 else if(costheta>0 && costheta<0.5*step) 00523 dedx_cos=cos_k[n+1]*costheta+cos_b[n+1]; 00524 else 00525 dedx_cos=cos_k[n]*costheta+cos_b[n]; 00526 00527 //conside physical edge around 0.93 00528 /* 00529 if(costheta>0.9125 && costheta<0.95){ 00530 n= (int)((0.9124+1.00-0.5*step)/step); 00531 dedx_cos=cos_k[n]*costheta+cos_b[n]; 00532 } 00533 if(costheta>-0.95 && costheta<-0.9125){ 00534 n= (int)((-0.9124+1.00-0.5*step)/step); 00535 dedx_cos=cos_k[n]*costheta+cos_b[n]; 00536 } 00537 */ 00538 00539 if(dedx_cos>0){ 00540 dedx_cos = dedx/dedx_cos; 00541 return dedx_cos; 00542 } 00543 else return dedx; 00544 }
|
|
00067 {return (x*x*x);}
|
|
00067 {return (x*x*x);}
|
|
|
|
02072 { 02073 //cout<<"DedxCorrecSvc::D2I"<<endl; 02074 double absCosTheta = fabs(cosTheta); 02075 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire 02076 double chargeDensity = D/projection; 02077 double numerator = 1 + m_alpha*chargeDensity; 02078 double denominator = 1 + m_gamma*chargeDensity; 02079 double I = D; 02080 02081 //if(denominator > 0.1) 02082 02083 I = D * m_ratio* numerator/denominator; 02084 02085 return I; 02086 }
|
|
|
|
00754 { 00755 }
|
|
|
|
00678 { 00679 // cout<<"DedxCorrecSvc::DocaSinCorrec"<<endl; 00680 //cout<<"doca = "<<doca<<" eangle = "<<eangle<<endl; 00681 if(eangle>0.25) eangle = eangle -0.5; 00682 else if(eangle < -0.25) eangle = eangle +0.5; 00683 double dedx_docaeangle; 00684 int iDoca; 00685 int iEAng = 0; 00686 //int iEAng; 00687 doca = doca/doca_norm[layer]; 00688 iDoca =(Int_t)floor((doca-Out_DocaMin)/Out_Stepdoca); 00689 if(doca<Out_DocaMin) iDoca=0; 00690 //if(iDoca < 0 ) return dedx; 00691 if(iDoca >=(Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca) ) 00692 iDoca = (Int_t)floor((Out_DocaMax-Out_DocaMin)/Out_Stepdoca)-1; 00693 //cout<<"iDoca : "<<iDoca<<" doca : "<<doca<<" Out_DocaMin : "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<endl; 00694 00695 //iEAng = (Int_t)floor((eangle-Phi_Min)/IEangle_step); 00696 for(int i =0; i<14; i++){ 00697 //iEAng =0; 00698 if(eangle <Eangle_value_cut[i] || eangle > Eangle_value_cut[i+1]) continue; 00699 else if(eangle>= Eangle_value_cut[i] && eangle < Eangle_value_cut[i+1]) { 00700 for(int k =0; k<i; k++){ 00701 iEAng+= Eangle_cut_bin[k]; 00702 } 00703 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i]; 00704 int temp_bin = int((eangle-Eangle_value_cut[i])/eangle_bin_cut_step); 00705 iEAng+= temp_bin; 00706 } 00707 } 00708 //cout<<iDoca <<" "<<iEAng<<endl; 00709 if(m_docaeangle[iDoca][iEAng]>0) { 00710 dedx_docaeangle = dedx/m_docaeangle[iDoca][iEAng]; 00711 } 00712 else{ 00713 return dedx; 00714 } 00715 00716 //auxillary 1d entrance angle correction 00717 if(eangle>-0.25 && eangle<0.25){ 00718 double stepsize= 0.5/m_venangle.size(); 00719 int nsize= m_venangle.size(); 00720 double slope=0; 00721 double offset=1; 00722 double y1=0,y2=0,x1=0,x2=0; 00723 00724 if(eangle>=-0.25+0.5*stepsize && eangle<=0.25-0.5*stepsize){ 00725 int bin = (int)( (eangle-(-0.25+0.5*stepsize))/stepsize); 00726 y1=m_venangle[bin]; 00727 x1=-0.25+0.5*stepsize+bin*stepsize; 00728 y2=m_venangle[bin+1]; 00729 x2=-0.25+1.5*stepsize+bin*stepsize; 00730 } 00731 else if(eangle<=-0.25+0.5*stepsize){ 00732 y1=m_venangle[0]; 00733 x1=-0.25+0.5*stepsize; 00734 y2=m_venangle[1]; 00735 x2=-0.25+1.5*stepsize; 00736 } 00737 else if( eangle>=0.25-0.5*stepsize){ 00738 y1=m_venangle[nsize-2]; 00739 x1=0.25-1.5*stepsize; 00740 y2=m_venangle[nsize-1]; 00741 x2=0.25-0.5*stepsize; 00742 } 00743 double kcof= (y2-y1)/(x2-x1); 00744 double bcof= y1-kcof*x1; 00745 double ratio = kcof*eangle+bcof; 00746 dedx_docaeangle=dedx_docaeangle/ratio; 00747 } 00748 00749 00750 return dedx_docaeangle; 00751 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00619 { 00620 // cout<<"DedxCorrecSvc::DriftDistCorrec"<<endl; 00621 double dedx_ddg; 00622 //cout<<"m_par_flag = "<<m_par_flag<<endl; 00623 //cout<<"dd vaule is = "<<dd<<endl; 00624 dd = fabs(dd); 00625 if(m_par_flag == 1) { 00626 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*dd + 00627 m_ddg[2][layer]*dd*dd + m_ddg[3][layer]*pow(dd,3); 00628 } else if(m_par_flag == 0) { 00629 dedx_ddg = m_ddg[0][layer] + m_ddg[1][layer]*T1(dd) + 00630 m_ddg[2][layer]*T2(dd) + m_ddg[3][layer]*T3(dd); 00631 } 00632 //cout<<"dedx_ddg is = "<<dedx_ddg<<endl; 00633 dedx_ddg = dedx/dedx_ddg; 00634 if (dedx_ddg>0.0) return dedx_ddg; 00635 else return dedx; 00636 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00640 { 00641 // cout<<"DedxCorrecSvc::EntaCorrec"<<endl; 00642 double dedx_enta; 00643 if(m_par_flag == 1) { 00644 dedx_enta = m_enta[0][layer] + m_enta[1][layer]*sinenta + 00645 m_enta[2][layer]*sinenta*sinenta + m_enta[3][layer]*pow(sinenta,3); 00646 } else if(m_par_flag == 0) { 00647 dedx_enta = m_enta[0][layer] + m_enta[1][layer]*T1(sinenta) + 00648 m_enta[2][layer]*T2(sinenta) + m_enta[3][layer]*T3(sinenta); 00649 } 00650 dedx_enta = dedx/dedx_enta; 00651 if (dedx_enta>0.0) return dedx_enta; 00652 else return dedx; 00653 00654 //return dedx; 00655 }
|
|
|
|
00091 { 00092 // cout<<"DedxCorrecSvc::finalize()"<<endl; 00093 MsgStream log(messageService(), name()); 00094 log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq; 00095 return StatusCode::SUCCESS; 00096 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00758 { 00759 if( m_mdcg_flag == 0 ) return dedx; 00760 double calib_ex = dedx; 00761 if( m_dedx_gain > 0 ) calib_ex /= m_dedx_gain; 00762 return calib_ex; 00763 }
|
|
|
|
00580 { 00581 // cout<<"DedxCorrecSvc::HadronCorrec"<<endl; 00582 //constant 1.00 in the following function is the mean dedx of normalized electrons. 00583 dedx=dedx/550.00; 00584 return D2I(costheta, I2D(costheta,1.00)/1.00*dedx)*550; 00585 }
|
|
|
|
00098 { 00099 // cout<<"DedxCorrecSvc::handle"<<endl; 00100 MsgStream log( messageService(), name() ); 00101 log << MSG::DEBUG << "handle: " << inc.type() << endreq; 00102 00103 if ( inc.type() == "NewRun" ){ 00104 log << MSG::DEBUG << "New Run" << endreq; 00105 00106 if( ! m_initConstFlg ) 00107 { 00108 if( m_init == 0 ) { init_param(); } 00109 else if( m_init ==1 ) { init_param_svc(); } 00110 /* if( ! init_param() ){ 00111 log << MSG::ERROR 00112 << "can not initilize Mdc Calib Constants" << endreq; 00113 }*/ 00114 } 00115 } 00116 }
|
|
|
|
02090 { 02091 // cout<<" DedxCorrecSvc::I2D"<<endl; 02092 double absCosTheta = fabs(cosTheta); 02093 double projection = pow(absCosTheta,m_power) + m_delta; // Projected length on wire 02094 02095 // Do the quadratic to compute d of the electron 02096 double a = m_alpha / projection; 02097 double b = 1 - m_gamma / projection*(I/m_ratio); 02098 double c = -I/m_ratio; 02099 02100 if(b==0 && a==0){ 02101 cout<<"both a and b coefficiants for hadron correction are 0"<<endl; 02102 return I; 02103 } 02104 02105 double D = a != 0? (-b + sqrt(b*b - 4.0*a*c))/(2.0*a) : -c/b; 02106 if(D<0){ 02107 cout<<"D is less 0! will try another solution"<<endl; 02108 D=a != 0? (-b - sqrt(b*b + 4.0*a*c))/(2.0*a) : -c/b; 02109 if(D<0){ 02110 cout<<"D is still less 0! just return uncorrectecd value"<<endl; 02111 return I; 02112 } 02113 } 02114 02115 return D; 02116 }
|
|
|
|
01016 { 01017 01018 //cout<<"DedxCorrecSvc::init_param()"<<endl; 01019 01020 for( int i = 0; i<6796 ; i++) { 01021 m_valid[i] = 1.0; 01022 m_wire_g[i] = 1.0; 01023 } 01024 m_dedx_gain = 1.0; 01025 m_dedx_resol = 1.0; 01026 for(int j = 0; j<100; j++){ 01027 m_rung[0][j] =1; 01028 } 01029 for(int j = 0; j<43; j++) { 01030 m_layer_g[j] = 1.0; 01031 m_ggs[0][j] = 1.0; 01032 m_ddg[0][j] = 1.0; 01033 m_enta[0][j] = 1.0; 01034 m_zdep[0][j] = 1.0; 01035 for(int k = 1; k<4; k++ ) { 01036 m_ggs[k][j] = 0.0; 01037 m_ddg[k][j] = 0.0; 01038 m_enta[k][j] = 0.0; 01039 m_zdep[k][j] = 0.0; 01040 } 01041 } 01042 01043 std::cout<<"DedxCorrecSvc::init_param()!!!"<<std::endl; 01044 std::cout<<"hello,init_param"<<endl; 01045 m_initConstFlg =true; 01046 01047 }
|
|
|
|
00119 { 00120 // cout<<"DedxCorrecSvc::init_param_svc()"<<endl; 00121 MsgStream log(messageService(), name()); 00122 IDataProviderSvc* pCalibDataSvc; 00123 StatusCode sc = service("CalibDataSvc", pCalibDataSvc, true); 00124 if ( !sc.isSuccess() ) { 00125 log << MSG::ERROR 00126 << "Could not get IDataProviderSvc interface of CalibXmlCnvSvc" 00127 << endreq; 00128 return ; 00129 } else { 00130 log << MSG::DEBUG 00131 << "Retrieved IDataProviderSvc interface of CalibXmlCnvSvc" 00132 << endreq; 00133 } 00134 00135 Iner_Stepdoca = (Iner_DocaMax-Iner_DocaMin)/NumSlicesDoca; 00136 Out_Stepdoca = (Out_DocaMax-Out_DocaMin)/NumSlicesDoca; 00137 //cout<<"Out_DocaMax: "<<Out_DocaMax<<" Out_DocaMin: "<<Out_DocaMin<<" Out_Stepdoca : "<<Out_Stepdoca<<" NumSlicesDoca : "<<NumSlicesDoca<<endl; 00138 std::string fullPath = "/Calib/DedxCal"; 00139 00140 SmartDataPtr<CalibData::DedxCalibData> test(pCalibDataSvc, fullPath); 00141 00142 //rung par---------------------------------------------------------- 00143 N_run = test->getrunNO(); 00144 cout<<"N_run = "<<N_run<<endl; 00145 for(int i=0;i<4;i++) { 00146 for(int j=0;j<N_run;j++) { 00147 m_rung[i][j] = test->getrung(i,j); 00148 //std::cout<<"rung["<<i<<"]["<<j<<"]="<<m_rung[i][j]<<endl; 00149 //std::cout<<"\n"; 00150 } 00151 } 00152 00153 //TFile* fruncalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/rung_calib/first_calib/DedxConst.root"); 00154 /*TFile* fruncalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/rung_calib/fourth_calib/temp/DedxConst.root"); 00155 00156 double rungain, runmean, runresol; 00157 int runno; 00158 TTree *rungtree= (TTree*) fruncalib-> Get("runcalib"); 00159 rungtree -> SetBranchAddress("rungain", &rungain); 00160 rungtree -> SetBranchAddress("runmean", &runmean); 00161 rungtree -> SetBranchAddress("runno", &runno); 00162 rungtree -> SetBranchAddress("runresol", &runresol); 00163 rungtree -> GetEntry(0); 00164 int NRUN = rungtree -> GetEntries(); 00165 N_run = rungtree -> GetEntries(); 00166 //cout<<"N_run = "<<N_run<<endl; 00167 //cout<<"NRUN = "<<NRUN<<endl; 00168 for(int runid =0; runid<NRUN; runid++) { 00169 rungtree->GetEntry(runid); 00170 m_rung[0][runid] = rungain; 00171 m_rung[1][runid] = runmean; 00172 m_rung[2][runid] = runno; 00173 m_rung[3][runid] = runresol; 00174 //cout<<"runid :"<<runid<<" run_no = "<<m_rung[2][runid]<<" run_mean = "<<m_rung[1][runid]<<" run_gain = "<<m_rung[0][runid]<<" runresol = "<<m_rung[3][runid]<<endl; 00175 } 00176 00177 //cout<<"run by run const ------------------------------------------------- "<<endl; 00178 for(int i=0;i<4;i++) { 00179 for(int j=0;j<NRUN;j++) { 00180 //std::cout<<"rung["<<i<<"]["<<j<<"]="<<m_rung[i][j]<<endl; 00181 //std::cout<<"\n"; 00182 } 00183 }*/ 00184 00185 for(int i=0;i<4;i++) { 00186 for(int j=0;j<43;j++) { 00187 m_ddg[i][j] = test->getddg(i,j); 00188 m_ggs[i][j] = test->getggs(i,j); 00189 m_zdep[i][j] = test->getzdep(i,j); 00190 // m_enta[i][j] = test->getenta(i,j); 00191 //std::cout<<"ddg["<<i<<"]["<<j<<"]="<<m_ddg[i][j]; 00192 //std::cout<<" ggs["<<i<<"]["<<j<<"]="<<m_ggs[i][j]; 00193 //std::cout<<" zdep["<<i<<"]["<<j<<"]="<<m_zdep[i][j]; 00194 //std::cout<<"\n"; 00195 } 00196 } 00197 00198 m_dedx_gain = test->getgain(); 00199 std::cout<<"gain="<<m_dedx_gain<<"\n"; 00200 00201 m_dedx_resol = test->getresol(); 00202 std::cout<<"resol="<<m_dedx_resol<<"\n"; 00203 00204 for(int i=0;i<43;i++) { 00205 m_layer_g[i] = test -> getlayerg(i); 00206 if(m_layer_g[i]>2.0 || m_layer_g[i]<0.5) m_layer_g[i] =1; 00207 //std::cout<<"layerg["<<i<<"]="<<m_layer_g[i]<<endl; 00208 } 00209 00210 /*const int n = 6796; 00211 double wireg[n]; 00212 TFile* fwiregcalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/wireg_calib/wiregcalib.root"); 00213 TTree *wiregtree= (TTree*)fwiregcalib -> Get("wiregcalib"); 00214 wiregtree -> SetBranchAddress("wireg", &wireg); 00215 wiregtree -> GetEntry(0); 00216 //cout<<"wire gain ------------------------------------------------- "<<endl; 00217 for(int i=0; i<n; i++){ 00218 m_wire_g[i] = wireg[i]; 00219 //std::cout<<"wireg["<<i<<"] = "<<m_wire_g[i]<<"\n"; 00220 }*/ 00221 00222 for(int i=0;i<6796;i++) { 00223 m_wire_g[i] = test -> getwireg(i); 00224 //std::cout<<"wireg["<<i<<"] = "<<m_wire_g[i]<<"\n"; 00225 } 00226 00227 for(int i=0;i<6796;i++) { 00228 m_valid[i] = 1; 00229 /*if(i<=483) {m_valid[i] =0;}*/ 00230 for(int j=0; j<12; j++){ 00231 if( i == dead_wire[j] ) m_valid[i] = 0; 00232 } 00233 //std::cout<<"valid["<<i<<"]="<<m_valid[i]<<"\n"; 00234 } 00235 00236 //tempoary way to get costheta constants 00237 cos_k.clear(); 00238 cos_b.clear(); 00239 //double cos1[80]; 00240 if(true){ 00241 const int nbin=80; 00242 vector<double> cos; 00243 /*TFile* fcosthecalib=new TFile("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/costhe_calib/costhetacalib2.root"); 00244 TTree *costhetree= (TTree*)fcosthecalib -> Get("costhetacalib"); 00245 TBranch* b_costheta; 00246 double costheta; 00247 costhetree->SetBranchAddress("costheta",&costheta,&b_costheta); 00248 //const int nbin=tc->GetEntries(); 00249 //cout<<"costheta gain ------------------------------------------------- "<<endl; 00250 for(int i=0; i<nbin; i++){ 00251 costhetree->GetEntry(i); 00252 cos.push_back(costheta); 00253 //cout<<"i : "<<i<<" costheta gain : "<<costheta<<endl; 00254 }*/ 00255 00256 for(int i=0; i<80; i++){ 00257 cos.push_back(test -> get_costheta(i)); 00258 //cout<<"i : "<<i<<" costheta gain : "<<test -> get_costheta(i)<<endl; 00259 } 00260 for(int i=0;i<nbin-1;i++){ 00261 double k=0; 00262 double b=0; 00263 if(cos[i]>0){ 00264 if(cos[i+1]>0){ 00265 double x1=-1.00+(0.5+i)*(2.00/nbin); 00266 double y1=cos[i]; 00267 double x2=-1.00+(0.5+i+1)*(2.00/nbin); 00268 double y2=cos[i+1]; 00269 k=(y2-y1)/(x2-x1); 00270 b=(y1*x2-y2*x1)/(x2-x1); 00271 } 00272 else if(i>0){ 00273 if(cos[i-1]>0){ 00274 double x1=-1.00+(0.5+i-1)*(2.00/nbin); 00275 double y1=cos[i-1]; 00276 double x2=-1.00+(0.5+i)*(2.00/nbin); 00277 double y2=cos[i]; 00278 k=(y2-y1)/(x2-x1); 00279 b=(y1*x2-y2*x1)/(x2-x1); 00280 } 00281 } 00282 } 00283 else { 00284 if(i>0){ 00285 if(cos[i-1]>0 && cos[i+1]>0){ 00286 double x1=-1.00+(0.5+i-1)*(2.00/nbin); 00287 double y1=cos[i-1]; 00288 double x2=-1.00+(0.5+i+1)*(2.00/nbin); 00289 double y2=cos[i+1]; 00290 k=(y2-y1)/(x2-x1); 00291 b=(y1*x2-y2*x1)/(x2-x1); 00292 } 00293 } 00294 } 00295 cos_k.push_back(k); 00296 cos_b.push_back(b); 00297 } 00298 } 00299 00300 t0_k.clear(); 00301 t0_b.clear(); 00302 if(true){ 00303 //const int nbin=35; 00304 /*vector<double> xbin; 00305 vector<double> ybin; 00306 TFile* ft0calib=new TFile("/ihepbatch/besd13/chunlei/calib/t0calib.root"); 00307 TTree *tree_t0=(TTree*)ft0calib->Get("t0calib"); 00308 TBranch* b_t0; 00309 TBranch* b_dedx; 00310 double t0, dedx; 00311 tree_t0->SetBranchAddress("t0",&t0,&b_t0); 00312 tree_t0->SetBranchAddress("dedx",&dedx,&b_dedx); 00313 const int nbin=tree_t0->GetEntries(); 00314 //cout<<"costheta t0 ------------------------------------------------- "<<endl; 00315 for(int i=0; i<tree_t0->GetEntries(); i++){ 00316 tree_t0->GetEntry(i); 00317 xbin.push_back(t0); 00318 ybin.push_back(dedx); 00319 //cout<<"xbin = "<<t0<<" ybin = "<<dedx<<endl; 00320 }*/ 00321 00322 const int nbin=35; 00323 vector<double> xbin; 00324 vector<double> ybin; 00325 for(int i=0; i<35; i++){ 00326 xbin.push_back(test ->get_t0(i)); 00327 ybin.push_back(test ->get_dedx(i)); 00328 //cout<<"xbin = "<<test ->get_t0(i)<<" ybin = "<<test ->get_dedx(i)<<endl; 00329 } 00330 for(int i=0;i<nbin-1;i++){ 00331 double k=0; 00332 double b=0; 00333 if(ybin[i]>0){ 00334 if(ybin[i+1]>0){ 00335 double x1=xbin[i]; 00336 double y1=ybin[i]; 00337 double x2=xbin[i+1]; 00338 double y2=ybin[i+1]; 00339 k=(y2-y1)/(x2-x1); 00340 b=(y1*x2-y2*x1)/(x2-x1); 00341 } 00342 else if(i>0){ 00343 if(ybin[i-1]>0){ 00344 double x1=xbin[i-1]; 00345 double y1=ybin[i-1]; 00346 double x2=xbin[i]; 00347 double y2=ybin[i]; 00348 k=(y2-y1)/(x2-x1); 00349 b=(y1*x2-y2*x1)/(x2-x1); 00350 } 00351 } 00352 } 00353 else { 00354 if(i>0){ 00355 if(ybin[i-1]>0 && ybin[i+1]>0){ 00356 double x1=xbin[i-1]; 00357 double y1=ybin[i-1]; 00358 double x2=xbin[i+1]; 00359 double y2=ybin[i+1]; 00360 k=(y2-y1)/(x2-x1); 00361 b=(y1*x2-y2*x1)/(x2-x1); 00362 } 00363 } 00364 } 00365 t0_k.push_back(k); 00366 t0_b.push_back(b); 00367 } 00368 } 00369 00370 if(true){ 00371 /*TFile fconst9("/home/bes/xcao/besd26/Calib/651_Calib/Psai_cal/doca_eangle_calib/check/doca_eangle_kal_doca110.root"); 00372 const int n = 1600; 00373 double Iner_gain[n], Iner_chi[n], Iner_hits[n]; 00374 double Out_gain[n], Out_chi[n], Out_hits[n]; 00375 double Id_doca[n], Ip_eangle[n]; 00376 00377 TTree *tree_docasin= (TTree*)fconst9.Get("docasincalib"); 00378 tree_docasin -> SetBranchAddress("Iner_gain", &Iner_gain); 00379 tree_docasin -> SetBranchAddress("Iner_chi", &Iner_chi); 00380 tree_docasin -> SetBranchAddress("Iner_hits", &Iner_hits); 00381 tree_docasin -> SetBranchAddress("Out_gain", &Out_gain); 00382 tree_docasin -> SetBranchAddress("Out_chi", &Out_chi); 00383 tree_docasin -> SetBranchAddress("Out_hits", &Out_hits); 00384 tree_docasin -> SetBranchAddress("Id_doca", &Id_doca); 00385 tree_docasin -> SetBranchAddress("Ip_eangle", &Ip_eangle); 00386 tree_docasin -> GetEntry(0); 00387 double docaeangle_gain[n]; 00388 double docaeangle_chisq[n]; 00389 double docaeangle_entries[n]; 00390 double cell[n]; 00391 //cout<<"doca eangle gain ------------------------------------------------- "<<endl; 00392 for(int i=0; i<n; i++){ 00393 tree_docasin->GetEntry(i); 00394 cell[i] = i; 00395 docaeangle_gain[i] = Out_gain[i]; 00396 docaeangle_chisq[i] = Out_chi[i]; 00397 docaeangle_entries[i] = Out_hits[i]; 00398 int Id_Doca_temp = Id_doca[i]; 00399 int Ip_Eangle_temp = Ip_eangle[i]; 00400 m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = Out_gain[i]; 00401 //cout<<i<<" "<<Id_Doca_temp<<" "<<Ip_Eangle_temp<<" "<<docaeangle_gain[i]<<endl; 00402 //m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = docaeangle_gain[i]; 00403 //cout<<i<<" "<<Id_doca[i]<<" "<<Ip_eangle[i]<<" Entries : "<<docaeangle_entries[i]<<" Gain value : "<<docaeangle_gain[i]<<" chisq : "<<docaeangle_chisq[i] <<endl; 00404 00405 } 00406 m_docaeangle[38][28]=0.72805;*/ 00407 00408 const int n = 1600; 00409 double docaeangle_gain[n]; 00410 double docaeangle_chisq[n]; 00411 double docaeangle_entries[n]; 00412 double cell[n]; 00413 for(int i=0; i<n; i++){ 00414 cell[i] = i; 00415 docaeangle_gain[i] = test -> get_out_gain(i); 00416 docaeangle_chisq[i] = test -> get_out_chi(i); 00417 docaeangle_entries[i] = test -> get_out_hits(i); 00418 int Id_Doca_temp = test -> get_id_doca(i); 00419 int Ip_Eangle_temp = test -> get_ip_eangle(i); 00420 m_docaeangle[Id_Doca_temp][Ip_Eangle_temp] = docaeangle_gain[i]; 00421 //cout<<i<<" "<<Id_Doca_temp<<" "<<Ip_Eangle_temp<<" "<<docaeangle_gain[i]<<endl; 00422 } 00423 } 00424 00425 00426 //get 1d entrance angle constants 00427 int onedsize=test->get_enanglesize(); 00428 m_venangle.clear(); 00429 for(int i=0; i< onedsize; i++){ 00430 m_venangle.push_back(test->get_enangle(i)); 00431 } 00432 00433 00434 //temporary way to get hadron saturation constants 00435 //for Psai hadron saturation 00436 00437 //only valide for jpsi data now !!!!! 00438 //have to find ways pass constans for different data 00439 /* m_alpha= 1.35630e-02; 00440 m_gamma= -6.78907e-04; 00441 m_delta= 1.18037e-02; 00442 m_power= 1.66268e+00; 00443 m_ratio= 9.94775e-01;*/ 00444 00445 const int hadronNo=test -> get_hadronNo(); 00446 // cout<<"@@@hadronNO===="<<hadronNo<<endl; 00447 m_alpha=test -> get_hadron(0); 00448 // cout<<"@@@@m_alpha===="<<m_alpha<<endl; 00449 m_gamma=test -> get_hadron(1); 00450 // cout<<"@@@@m_gamma===="<<m_gamma<<endl; 00451 m_delta=test -> get_hadron(2); 00452 // cout<<"@@@@m_delta====="<<m_delta<<endl; 00453 m_power=test -> get_hadron(3); 00454 // cout<<"@@@@m_power====="<<m_power<<endl; 00455 m_ratio=test -> get_hadron(4); 00456 // cout<<"@@@m_ratio======"<<m_ratio<<endl; 00457 00458 00459 00460 /* 00461 m_delta=0.1; 00462 m_alpha=0.0282; 00463 m_gamma=-0.030; 00464 m_power=1.0; 00465 m_ratio=1.0; 00466 //m_delta=0.1; 00467 //m_alpha=0.0382; 00468 //m_gamma=-0.0438; 00469 */ 00470 00471 //m_initConstFlg =true; 00472 }
|
|
|
|
00061 { 00062 // cout<<"DedxCorrecSvc::initialize"<<endl; 00063 MsgStream log(messageService(), name()); 00064 log << MSG::INFO << name() << "DedxCorrecSvc::initialize()" << endreq; 00065 00066 StatusCode sc = Service::initialize(); 00067 if( sc.isFailure() ) return sc; 00068 00069 IIncidentSvc* incsvc; 00070 sc = service("IncidentSvc", incsvc); 00071 int priority = 100; 00072 if( sc.isSuccess() ){ 00073 //incsvc -> addListener(this, "BeginEvent", priority); 00074 incsvc -> addListener(this, "NewRun", priority); 00075 } 00076 00077 StatusCode scgeo = service("MdcGeomSvc", geosvc); 00078 if(scgeo.isFailure() ) { 00079 log << MSG::ERROR << "GeoSvc failing!!!!!!!" << endreq; 00080 return scgeo; 00081 } 00082 00083 StatusCode scmgn = service ("MagneticFieldSvc",m_pIMF); 00084 if(scmgn!=StatusCode::SUCCESS) { 00085 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq; 00086 } 00087 00088 return StatusCode::SUCCESS; 00089 }
|
|
00015 { return IID_IDedxCorrecSvc; }
|
|
00015 { return IID_IDedxCorrecSvc; }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00818 { 00819 //cout<<"DedxCorrecSvc::LayerCorrec"<<endl; 00820 00821 if( m_zdep_flag == 0 && m_ggs_flag == 0 ) return ex; 00822 00823 double calib_ex = ZdepCorrec( layer, z, ex ); 00824 if( m_ggs_flag != 0 ) calib_ex = SaturCorrec( layer, costheta, calib_ex ); 00825 return calib_ex; 00826 00827 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00587 { 00588 // cout<<"DedxCorrecSvc::LayerGainCorrec"<<endl; 00589 if( m_layer_g[layid] > 0 ){ 00590 double ch_dedx = dedx/m_layer_g[layid]; 00591 return ch_dedx; 00592 } 00593 else if( m_layer_g[layid] == 0 ){ 00594 return dedx; 00595 } 00596 else return 0; 00597 }
|
|
Implements IDedxCorrecSvc. |
|
***********-------------------***************------------------****************// ***********-------------------***************------------------****************// assumed the first half circle ***********-------------------***************------------------****************// ***********-------------------***************------------------****************// in the Local Helix case, phi must be small Implements IDedxCorrecSvc. 01407 { 01408 01409 double length = geosvc->Layer(layer)->Length(); 01410 int genlay = geosvc->Layer(layer)->Gen(); 01411 double East_lay_X = geosvc->GeneralLayer(genlay)->SxEast(); 01412 double East_lay_Y = geosvc->GeneralLayer(genlay)->SyEast(); 01413 double East_lay_Z = geosvc->GeneralLayer(genlay)->SzEast(); 01414 double West_lay_X = geosvc->GeneralLayer(genlay)->SxWest(); 01415 double West_lay_Y = geosvc->GeneralLayer(genlay)->SyWest(); 01416 double West_lay_Z = geosvc->GeneralLayer(genlay)->SzWest(); 01417 01418 HepPoint3D East_origin(East_lay_X, East_lay_Y, East_lay_Z); 01419 HepPoint3D West_origin(West_lay_X, West_lay_Y, West_lay_Z); 01420 Hep3Vector wire = (CLHEP::Hep3Vector)East_origin - (CLHEP::Hep3Vector)West_origin; 01421 HepPoint3D piovt_z =(z*10+length/2 )/length * wire + West_origin; 01422 piovt_z = piovt_z*0.1; // conversion of the units(mm=>cm) 01423 01424 01425 //-------------------------------temp -------------------------------// 01426 HepPoint3D piv(hel.pivot()); 01427 //-------------------------------temp -------------------------------// 01428 01429 double dr0 = hel.a()[0]; 01430 double phi0 = hel.a()[1]; 01431 double kappa = hel.a()[2]; 01432 double dz0 = hel.a()[3]; 01433 double tanl = hel.a()[4]; 01434 01435 // Choose the local field !! 01436 double Bz = 1000*(m_pIMF->getReferField()); 01437 double ALPHA_loc = 1000/(2.99792458*Bz); 01438 01439 // Choose the local field !! 01440 //double Bz = 1.0; // will be obtained from magnetic field service 01441 //double ALPHA_loc = 1000/(2.99792458*Bz); 01442 //double ALPHA_loc = 333.564095; 01443 int charge = ( kappa >= 0 )? 1 : -1; 01444 double rho = ALPHA_loc/kappa; 01445 double pt = fabs( 1.0/kappa ); 01446 double lambda = atan( tanl ); 01447 double theta = M_PI_2- lambda; 01448 double sintheta = sin(M_PI_2-atan(tanl)); 01449 // theta = 2.0*M_PI*theta/360.; 01450 double phi = fmod(phi0 + M_PI*4, M_PI*2); 01451 double csf0 = cos(phi); 01452 double snf0 = (1. - csf0) * (1. + csf0); 01453 snf0 = sqrt((snf0 > 0.) ? snf0 : 0.); 01454 if(phi > M_PI) snf0 = - snf0; 01455 //if(ntpFlag>0) 01456 //cout<<"rho = "<<rho<<" hel.dr() + rho = "<<hel.dr() + rho<<endl; 01457 01458 //-------------------------------temp -------------------------------// 01459 double x_c = piv.x() + ( hel.dr() + rho )*csf0; 01460 double y_c = piv.y() + ( hel.dr() + rho )*snf0; 01461 double z_c = piv.z() + hel.dz(); 01462 HepPoint3D ccenter(x_c, y_c, 0.0); 01463 double m_c_perp(ccenter.perp()); 01464 Hep3Vector m_c_unit((HepPoint3D)ccenter.unit()); 01465 //-------------------------------temp -------------------------------// 01466 01468 double x_c_boost = x_c - piovt_z.x(); 01469 double y_c_boost = y_c - piovt_z.y(); 01470 double z_c_boost = z_c - piovt_z.z(); 01471 HepPoint3D ccenter_boost(x_c_boost, y_c_boost, 0.0); 01472 double m_c_perp_boost(ccenter_boost.perp()); 01473 //if(ntpFlag>0) 01474 //cout<<" ccenter = "<<ccenter<<" m_c_perp ="<<m_c_perp<<endl; 01475 Hep3Vector m_c_unit_boost((HepPoint3D)ccenter_boost.unit()); 01476 01477 //phi region estimation 01478 double phi_io[2]; 01479 HepPoint3D IO = (-1)*charge*m_c_unit; 01480 double dphi0 = fmod( IO.phi()+4*M_PI,2*M_PI ) - phi; 01481 double IO_phi = fmod( IO.phi()+4*M_PI,2*M_PI ); 01482 //double dphi0_bak = IO_phi - phi; 01483 //if(dphi0 != 0) 01484 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl; 01485 if( dphi0 > M_PI ) dphi0 -= 2*M_PI; 01486 else if( dphi0 < -M_PI ) dphi0 += 2*M_PI; 01487 //cout<<"dphi value is : "<<dphi0<<" phi:"<<phi<<" IO.phi:"<<IO_phi<<endl; 01488 //cout<<"charge is = "<<charge<<endl; 01489 phi_io[0] = -(1+charge)*M_PI_2 - charge*dphi0; 01490 //phi_io[0] = -(1+charge)*M_PI_2 + charge*dphi0; 01491 phi_io[1] = phi_io[0]+1.5*M_PI; 01492 //cout<<"phi_io[0] is : "<<phi_io[0]<<" phi_io[1]:"<<phi_io[1]<<endl; 01493 double m_crio[2]; 01494 double m_zb, m_zf, Calpha; 01495 01496 // retrieve Mdc geometry information 01497 double rcsiz1= geosvc->Layer(layer)->RCSiz1(); 01498 double rcsiz2= geosvc->Layer(layer)->RCSiz2(); 01499 double rlay= geosvc->Layer(layer)->Radius(); 01500 int ncell= geosvc->Layer(layer)->NCell(); 01501 double phioffset= geosvc->Layer(layer)->Offset(); 01502 float shift= geosvc->Layer(layer)->Shift(); 01503 double slant= geosvc->Layer(layer)->Slant(); 01504 //double length = geosvc->Layer(layer)->Length(); 01505 int type = geosvc->Layer(layer)->Sup()->Type(); 01506 01508 //check the value if same // 01510 int w0id = geosvc->Layer(layer)->Wirst(); 01511 int wid = w0id + cellid; 01512 HepPoint3D backkward = geosvc->Wire(wid)->BWirePos(); 01513 HepPoint3D forward = geosvc->Wire(wid)->FWirePos(); 01514 double x_lay_backward = geosvc->Wire(layer, cellid)->Backward().x(); 01515 double y_lay_backward = geosvc->Wire(layer, cellid)->Backward().y(); 01516 double x_lay_forward = geosvc->Wire(layer, cellid)->Forward().x(); 01517 double y_lay_forward = geosvc->Wire(layer, cellid)->Forward().y(); 01518 double r_lay_backward = sqrt(x_lay_backward*x_lay_backward+y_lay_backward*y_lay_backward); 01519 double r_lay_forward = sqrt(x_lay_forward*x_lay_forward+y_lay_forward*y_lay_forward); 01520 double r_lay_use = ((z*10+length/2)/length)*(r_lay_backward-r_lay_forward) + r_lay_forward; 01521 /*for(int i=0; i<43; i++){ 01522 double r_up= geosvc->Layer(i)->RCSiz1(); 01523 double r_down= geosvc->Layer(i)->RCSiz2(); 01524 cout<<i<<" "<<r_up<<" "<<r_down<<endl; 01525 }*/ 01526 // shift = shift*type; 01527 // cout<< "type "<< type <<endl; 01528 r_lay_use = 0.1*r_lay_use; 01529 rcsiz1 = 0.1*rcsiz1; 01530 rcsiz2 = 0.1*rcsiz2; 01531 rlay = 0.1*rlay; 01532 length = 0.1*length; 01533 m_zb = 0.5*length; 01534 m_zf = -0.5*length; 01535 m_crio[0] = rlay - rcsiz1; 01536 m_crio[1] = rlay + rcsiz2; 01537 //conversion of the units(mm=>cm) 01538 int sign = -1; 01539 int epflag[2]; 01540 Hep3Vector iocand[2]; 01541 Hep3Vector cell_IO[2]; 01542 double dphi, downin; 01543 Hep3Vector zvector; 01544 //if (ntpFlag>0) cout<<"z = "<<z<<" , m_zb = "<<m_zb<<" , m_zf = "<<m_zf<<endl; 01545 if( type ) { 01546 downin = (z*z-m_zb*m_zb)*pow(tan(slant),2); 01547 m_crio[0] = sqrt(m_crio[0]*m_crio[0]+downin); 01548 m_crio[1] = sqrt(m_crio[1]*m_crio[1]+downin); 01549 } 01550 01551 //int stced[2]; 01552 01553 for( int i = 0; i < 2; i++ ) { 01554 double cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[i]*m_crio[i] - rho*rho; 01555 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[i] ); 01556 if(fabs(cos_alpha)>1&&i==0) return(-1.0); 01557 if(fabs(cos_alpha)>1&&i==1) { 01558 cos_alpha = m_c_perp_boost*m_c_perp_boost + m_crio[0]*m_crio[0] - rho*rho; 01559 cos_alpha = 0.5*cos_alpha/( m_c_perp_boost*m_crio[0] ); 01560 Calpha = 2.0*M_PI-acos( cos_alpha ); 01561 } else { 01562 Calpha = acos( cos_alpha ); 01563 } 01564 epflag[i] = 0; 01565 iocand[i] = m_c_unit_boost; 01566 iocand[i].rotateZ( charge*sign*Calpha ); 01567 iocand[i]*= m_crio[i]; 01568 //if (ntpFlag>0) cout<<"iocand corridate is : "<<iocand[i]<<endl; 01569 01571 //compare with standard coordinate system results // 01573 //-------------------------------temp-------------------------// 01574 iocand[i] = iocand[i]+ piovt_z; //change to standard coordinate system 01575 //iocand[i].y() = iocand[i].y() + piovt_z.y(); 01576 //if (ntpFlag>0) cout<<i<<setw(10)<<iocand[i].x()<<setw(10)<<iocand[i].y()<<endl; 01577 //------------------------------temp -------------------------// 01578 01579 double xx = iocand[i].x() - x_c; 01580 double yy = iocand[i].y() - y_c; 01581 //dphi = atan2(yy, xx) - phi0 - M_PI_2*(1+charge); 01582 dphi = atan2(yy, xx) - phi0 - M_PI_2*(1-charge); 01583 dphi = fmod( dphi + 8.0*M_PI, 2*M_PI ); 01584 01585 if( dphi < phi_io[0] ) { 01586 dphi += 2*M_PI; 01587 } 01588 else if( phi_io[1] < dphi ) { 01589 dphi -= 2*M_PI; 01590 } 01592 01593 //Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl); 01594 Hep3Vector zvector( 0., 0., z_c-rho*dphi*tanl-piovt_z.z()); 01595 //if (ntpFlag>0) cout<<"z_c-rho*dphi*tanl : "<<z_c-rho*dphi*tanl<<endl; 01596 cell_IO[i] = iocand[i]; 01597 cell_IO[i] += zvector; 01598 //---------------------------------temp ---------------------------------// 01599 //cell_IO[i] = hel.x(dphi);//compare with above results 01600 //---------------------------------temp ---------------------------------// 01601 01602 double xcio, ycio, phip; 01604 /* 01605 float delrin = 2.0; 01606 if( m_zf-delrin > zvector.z() ){ 01607 phip = z_c - m_zb + delrin; 01608 phip = phip/( rho*tanl ); 01609 phip = phip + phi0; 01610 xcio = x_c - rho*cos( phip ); 01611 ycio = y_c - rho*sin( phip ); 01612 cell_IO[i].setX( xcio ); 01613 cell_IO[i].setY( ycio ); 01614 cell_IO[i].setZ( m_zb+delrin ); 01615 epflag[i] = 1; 01616 } 01617 else if( m_zb+delrin < zvector.z() ){ 01618 phip = z_c - m_zf-delrin; 01619 phip = phip/( rho*tanl ); 01620 phip = phip + phi0; 01621 xcio = x_c - rho*cos( phip ); 01622 ycio = y_c - rho*sin( phip ); 01623 cell_IO[i].setX( xcio ); 01624 cell_IO[i].setY( ycio ); 01625 cell_IO[i].setZ( m_zf-delrin ); 01626 epflag[i] = 1; 01627 } 01628 else{ 01629 */ 01630 // cell_IO[i] = iocand; 01631 // cell_IO[i] += zvector; 01632 // } 01633 //dphi = dphi -M_PI; 01634 cell_IO[i] = hel.x(dphi); 01635 //if (ntpFlag>0) { cout<<"cell_IO corridate : "<<cell_IO[i]<<endl;} 01636 // if(i==0) cout<<"zhit value is = "<<z<<endl; 01637 } 01638 01639 // path length estimation 01640 // phi calculation 01641 Hep3Vector cl = cell_IO[1] - cell_IO[0]; 01642 01643 //float ch_tphi, ch_tdphi; 01644 double ch_theta; 01645 double ch_dphi; 01646 double ch_ltrk = 0; 01647 double ch_ltrk_rp = 0; 01648 ch_dphi = cl.perp()*0.5/(ALPHA_loc*pt); 01649 ch_dphi = 2.0 * asin( ch_dphi ); 01650 ch_ltrk_rp = ch_dphi*ALPHA_loc*pt; 01651 double rpi_path = sqrt(cl.x()*cl.x()+cl.y()*cl.y()); 01652 ch_ltrk = sqrt( ch_ltrk_rp*ch_ltrk_rp + cl.z()*cl.z() ); 01653 double path = ch_ltrk_rp/ sintheta; 01654 ch_theta = cl.theta(); 01655 /* if (ntpFlag>0) 01656 { 01657 // if((ch_ltrk_rp-rpi_path)>0.001 || (ch_ltrk-path)>0.001) 01658 cout<<"ch_ltrk_rp : "<<ch_ltrk_rp<<" rpi_path: "<<rpi_path<<endl; 01659 cout<<"ch_ltrk : "<<ch_ltrk<<" path:"<<path<<endl; 01660 } 01661 */ 01662 /* 01663 if( ch_theta < theta*0.85 || 1.15*theta < ch_theta) 01664 ch_ltrk *= -1; /// miss calculation 01665 01666 if( epflag[0] == 1 || epflag [1] == 1 ) 01667 ch_ltrk *= -1; /// invalid region for dE/dx or outside of Mdc 01668 */ 01669 // judge how many cells are traversed and deal with different situations 01670 double phibin, phi_in, phi_out, inlow, inup, outlow, outup, gap, phi1, phi2, phi_mid, phi_midin, phi_midout; 01671 int cid_in, cid_out; 01672 double inlow_z, inup_z, outlow_z, outup_z, gap_z, phi1_z, phi2_z, phi_mid_z, phi_midin_z, phi_midout_z; 01673 //cache sampl in each cell for this layer 01674 01675 std::vector<double> sampl; 01676 sampl.resize(ncell); 01677 for(int k=0; k<ncell; k++) { 01678 sampl[k] = -1.0; 01679 } 01680 01681 cid_in = cid_out = -1; 01682 phi_in = cell_IO[0].phi(); 01683 phi_out = cell_IO[1].phi(); 01684 //phi_in = iocand[0].phi(); 01685 //phi_out = iocand[1].phi(); 01686 phi_in = fmod( phi_in+4*M_PI,2*M_PI ); 01687 phi_out = fmod( phi_out+4*M_PI,2*M_PI ); 01688 phibin = 2.0*M_PI/ncell; 01689 // gap = fabs(phi_out-phi_in); 01690 01691 //determine the in/out cell id 01692 Hep3Vector cell_mid=0.5*(cell_IO[0]+cell_IO[1]); 01693 //Hep3Vector cell_mid=0.5*(iocand[0]+iocand[0]); 01694 //cout<<cell_mid<<endl; 01695 //double stereophi = shift*phibin*(0.5-cell_mid.z()/length); 01696 //phioffset = phioffset+stereophi; 01697 double stphi[2], phioff[2]; 01698 stphi[0] = shift*phibin*(0.5-cell_IO[0].z()/length); 01699 stphi[1] = shift*phibin*(0.5-cell_IO[1].z()/length); 01700 //stphi[0] = shift*phibin*(0.5-z/length); 01701 //stphi[1] = shift*phibin*(0.5-z/length); 01702 phioff[0] = phioffset+stphi[0]; 01703 phioff[1] = phioffset+stphi[1]; 01704 01705 for(int i=0; i<ncell; i++) { 01706 01707 double x_lay_backward_cell = geosvc->Wire(layer, i)->Backward().x()*0.1; 01708 double y_lay_backward_cell = geosvc->Wire(layer, i)->Backward().y()*0.1; 01709 double x_lay_forward_cell = geosvc->Wire(layer, i)->Forward().x()*0.1; 01710 double y_lay_forward_cell = geosvc->Wire(layer, i)->Forward().y()*0.1; 01711 //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; 01712 //Hep3Vector lay_backward, lay_forward; 01713 Hep3Vector lay_backward(x_lay_backward_cell, y_lay_backward_cell, 0); 01714 Hep3Vector lay_forward(x_lay_forward_cell, y_lay_forward_cell, 0); 01715 Hep3Vector Cell_z[2]; 01716 Cell_z[0] = ((cell_IO[0].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward; 01717 Cell_z[1] = ((cell_IO[1].z()+length/2)/length)*(lay_backward - lay_forward) + lay_forward; 01718 double z_phi[2]; 01719 z_phi[0] = Cell_z[0].phi(); 01720 z_phi[1] = Cell_z[1].phi(); 01721 z_phi[0] = fmod( z_phi[0]+4*M_PI,2*M_PI ); 01722 z_phi[1] = fmod( z_phi[1]+4*M_PI,2*M_PI ); 01723 //double backward_phi = lay_backward.phi(); 01724 //double forward_phi = lay_forward.phi(); 01725 //backward_phi = fmod( backward_phi+4*M_PI,2*M_PI ); 01726 //forward_phi = fmod( forward_phi+4*M_PI,2*M_PI ); 01727 //if(ntpFlag>0) cout<<"backward_phi cal : "<<atan2(y_lay_backward,x_lay_backward)<<" forward_phi : "<<atan2(y_lay_forward,x_lay_forward)<<endl; 01728 //if(ntpFlag>0) cout<<layer<<" backward_phi : "<<backward_phi<<" forward_phi : "<<forward_phi<<endl; 01729 01730 //double z_phi[2]; 01731 //z_phi[0] = ((cell_IO[0].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi; 01732 //z_phi[1] = ((cell_IO[1].z()+length/2)/length)*(backward_phi - forward_phi) + forward_phi; 01733 //if(ntpFlag>0) cout<<"phi z : "<<z_phi[0]<<" "<<z_phi[1]<<endl; 01734 inlow_z = z_phi[0] - phibin*0.5; 01735 inup_z = z_phi[0] + phibin*0.5; 01736 outlow_z = z_phi[1] - phibin*0.5; 01737 outup_z = z_phi[1] + phibin*0.5; 01738 inlow_z = fmod( inlow_z+4*M_PI,2*M_PI ); 01739 inup_z = fmod( inup_z+4*M_PI,2*M_PI ); 01740 outlow_z = fmod( outlow_z+4*M_PI,2*M_PI ); 01741 outup_z = fmod( outup_z+4*M_PI,2*M_PI ); 01742 01743 01744 // for stereo layer 01745 inlow = phioff[0]+phibin*i-phibin*0.5; 01746 inup = phioff[0]+phibin*i+phibin*0.5; 01747 outlow = phioff[1]+phibin*i-phibin*0.5; 01748 outup = phioff[1]+phibin*i+phibin*0.5; 01749 inlow = fmod( inlow+4*M_PI,2*M_PI ); 01750 inup = fmod( inup+4*M_PI,2*M_PI ); 01751 outlow = fmod( outlow+4*M_PI,2*M_PI ); 01752 outup = fmod( outup+4*M_PI,2*M_PI ); 01753 01754 #ifdef YDEBUG 01755 if(ntpFlag > 0) cout << "shift " << shift 01756 <<" phi_in " << phi_in << " phi_out " << phi_out 01757 << " inup "<< inup << " inlow " << inlow 01758 << " outup "<< outup << " outlow " << outlow << endl; 01759 #endif 01760 01761 /*if(phi_in>=inlow&&phi_in<inup) cid_in = i; 01762 if(phi_out>=outlow&&phi_out<outup) cid_out = i; 01763 if ( inlow>inup) { 01764 if((phi_in>=inlow&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup)) cid_in = i; 01765 } 01766 if ( outlow>outup) { 01767 if((phi_out>=outlow&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup)) cid_out = i; 01768 }*/ 01769 01770 if(phi_in>=inlow_z&&phi_in<inup_z) cid_in = i; 01771 if(phi_out>=outlow_z&&phi_out<outup_z) cid_out = i; 01772 if ( inlow_z>inup_z) { 01773 if((phi_in>=inlow_z&&phi_in<2.0*M_PI)||(phi_in>=0.0&&phi_in<inup_z)) cid_in = i; 01774 } 01775 if ( outlow_z>outup_z) { 01776 if((phi_out>=outlow_z&&phi_out<2.0*M_PI)||(phi_out>=0.0&&phi_out<outup_z)) cid_out = i; 01777 } 01778 } 01779 01780 phi_midin = phi_midout = phi1 = phi2 = -999.0; 01781 gap = -999.0; 01782 //only one cell traversed 01783 //if(ntpFlag > 0) cout<<"cid_in = "<<cid_in<<" cid_out = "<<cid_out<<endl; 01784 if(cid_in == -1 || cid_out == -1) return -1; 01785 01786 if( cid_in == cid_out) { 01787 sampl[cid_in]= ch_ltrk; 01788 //return ch_ltrk; 01789 return sampl[cellid]; 01790 } else if(cid_in < cid_out) { 01791 //in cell id less than out cell id 01792 //deal with the special case crossing x axis 01793 if( cid_out-cid_in>ncell/2 ) { 01794 01795 // if(gap>=M_PI) gap = 2.0*M_PI-gap; 01796 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1; 01797 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1; 01798 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1; 01799 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1; 01800 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0); 01801 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0); 01802 Hep3Vector Cell_z[2]; 01803 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin; 01804 double phi_cin_z = Cell_z[0].phi(); 01805 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI ); 01806 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1; 01807 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1; 01808 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1; 01809 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1; 01810 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0); 01811 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0); 01812 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout; 01813 double phi_cout_z = Cell_z[1].phi(); 01814 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI ); 01815 01816 phi_midin_z = phi_cin_z-phibin*0.5; 01817 phi_midout_z = phi_cout_z+phibin*0.5; 01818 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI ); 01819 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI ); 01820 phi1_z = phi_midout_z-phi_out; 01821 phi2_z = phi_in-phi_midin_z; 01822 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI); 01823 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI); 01824 gap_z = phi1_z+phi2_z+(ncell-1-cid_out+cid_in)*phibin; 01825 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI); 01826 sampl[cid_in]=phi2_z/gap_z*ch_ltrk; 01827 sampl[cid_out]=phi1_z/gap_z*ch_ltrk; 01828 for( int jj = cid_out+1; jj<ncell; jj++) { 01829 sampl[jj]=phibin/gap_z*ch_ltrk; 01830 } 01831 for( int jj = 0; jj<cid_in; jj++) { 01832 sampl[jj]=phibin/gap_z*ch_ltrk; 01833 } 01834 01835 /* 01836 // if(gap>=M_PI) gap = 2.0*M_PI-gap; 01837 phi_midin = phioff[0]+phibin*cid_in-phibin*0.5; 01838 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5; 01839 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI ); 01840 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI ); 01841 phi1 = phi_midout-phi_out; 01842 phi2 = phi_in-phi_midin; 01843 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI); 01844 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI); 01845 gap = phi1+phi2+(ncell-1-cid_out+cid_in)*phibin; 01846 gap = fmod(gap+2.0*M_PI,2.0*M_PI); 01847 sampl[cid_in]=phi2/gap*ch_ltrk; 01848 sampl[cid_out]=phi1/gap*ch_ltrk; 01849 for( int jj = cid_out+1; jj<ncell; jj++) { 01850 sampl[jj]=phibin/gap*ch_ltrk; 01851 } 01852 for( int jj = 0; jj<cid_in; jj++) { 01853 sampl[jj]=phibin/gap*ch_ltrk; 01854 }*/ 01855 /* 01856 cout<<"LLLLLLLLLLLLL"<<endl; 01857 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift 01858 << " cid_in " << cid_in << " cid_out " << cid_out << endl; 01859 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out " 01860 << phi_out << " phi_midout " << phi_midout <<endl; 01861 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 " 01862 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl; 01863 */ 01864 } else { 01865 //normal case 01866 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1; 01867 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1; 01868 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1; 01869 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1; 01870 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0); 01871 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0); 01872 Hep3Vector Cell_z[2]; 01873 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin; 01874 double phi_cin_z = Cell_z[0].phi(); 01875 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI ); 01876 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1; 01877 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1; 01878 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1; 01879 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1; 01880 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0); 01881 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0); 01882 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout; 01883 double phi_cout_z = Cell_z[1].phi(); 01884 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI ); 01885 01886 phi_midin_z = phi_cin_z+phibin*0.5; 01887 phi_midout_z = phi_cout_z-phibin*0.5; 01888 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI ); 01889 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI ); 01890 phi1_z = phi_midin_z-phi_in; 01891 phi2_z = phi_out-phi_midout_z; 01892 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI); 01893 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI); 01894 gap_z = phi1_z+phi2_z+(cid_out-cid_in-1)*phibin; 01895 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI); 01896 sampl[cid_in]=phi1_z/gap_z*ch_ltrk; 01897 sampl[cid_out]=phi2_z/gap_z*ch_ltrk; 01898 for( int jj = cid_in+1; jj<cid_out; jj++) { 01899 sampl[jj]=phibin/gap_z*ch_ltrk; 01900 } 01901 //normal case 01902 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5; 01903 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5; 01904 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI ); 01905 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI ); 01906 phi1 = phi_midin-phi_in; 01907 phi2 = phi_out-phi_midout; 01908 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI); 01909 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI); 01910 gap = phi1+phi2+(cid_out-cid_in-1)*phibin; 01911 gap = fmod(gap+2.0*M_PI,2.0*M_PI); 01912 sampl[cid_in]=phi1/gap*ch_ltrk; 01913 sampl[cid_out]=phi2/gap*ch_ltrk; 01914 for( int jj = cid_in+1; jj<cid_out; jj++) { 01915 sampl[jj]=phibin/gap*ch_ltrk; 01916 }*/ 01917 } 01918 01919 } else if(cid_in > cid_out) { 01920 //in cell id greater than out cell id 01921 //deal with the special case crossing x axis 01922 if( cid_in-cid_out>ncell/2 ) { 01923 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1; 01924 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1; 01925 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1; 01926 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1; 01927 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0); 01928 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0); 01929 Hep3Vector Cell_z[2]; 01930 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin; 01931 double phi_cin_z = Cell_z[0].phi(); 01932 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI ); 01933 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1; 01934 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1; 01935 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1; 01936 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1; 01937 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0); 01938 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0); 01939 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout; 01940 double phi_cout_z = Cell_z[1].phi(); 01941 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI ); 01942 01943 phi_midin_z = phi_cin_z+phibin*0.5; 01944 phi_midout_z = phi_cout_z-phibin*0.5; 01945 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI ); 01946 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI ); 01947 phi1_z = phi_midin_z-phi_in; 01948 phi2_z = phi_out-phi_midout_z; 01949 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI); 01950 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI); 01951 gap_z = phi1_z+phi2_z+(ncell-1-cid_in+cid_out)*phibin; 01952 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI); 01953 sampl[cid_out]=phi2_z/gap_z*ch_ltrk; 01954 sampl[cid_in]=phi1_z/gap_z*ch_ltrk; 01955 for( int jj = cid_in+1; jj<ncell; jj++) { 01956 sampl[jj]=phibin/gap_z*ch_ltrk; 01957 } 01958 for( int jj = 0; jj<cid_out; jj++) { 01959 sampl[jj]=phibin/gap_z*ch_ltrk; 01960 } 01961 01962 // if(gap>=M_PI) gap = 2.0*M_PI-gap; 01963 /*phi_midin = phioff[0]+phibin*cid_in+phibin*0.5; 01964 phi_midout = phioff[1]+phibin*cid_out-phibin*0.5; 01965 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI ); 01966 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI ); 01967 phi1 = phi_midin-phi_in; 01968 phi2 = phi_out-phi_midout; 01969 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI); 01970 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI); 01971 gap = phi1+phi2+(ncell-1-cid_in+cid_out)*phibin; 01972 gap = fmod(gap+2.0*M_PI,2.0*M_PI); 01973 sampl[cid_out]=phi2/gap*ch_ltrk; 01974 sampl[cid_in]=phi1/gap*ch_ltrk; 01975 for( int jj = cid_in+1; jj<ncell; jj++) { 01976 sampl[jj]=phibin/gap*ch_ltrk; 01977 } 01978 for( int jj = 0; jj<cid_out; jj++) { 01979 sampl[jj]=phibin/gap*ch_ltrk; 01980 }*/ 01981 /* 01982 cout<<"LLLLLLLLLLLLL"<<endl; 01983 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift 01984 << " cid_in " << cid_in << " cid_out " << cid_out << endl; 01985 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out " 01986 << phi_out << " phi_midout " << phi_midout <<endl; 01987 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 " 01988 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl; 01989 */ 01990 } else { 01991 double x_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().x()*0.1; 01992 double y_lay_backward_cin = geosvc->Wire(layer, cid_in)->Backward().y()*0.1; 01993 double x_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().x()*0.1; 01994 double y_lay_forward_cin = geosvc->Wire(layer, cid_in)->Forward().y()*0.1; 01995 Hep3Vector lay_backward_cin(x_lay_backward_cin, y_lay_backward_cin, 0); 01996 Hep3Vector lay_forward_cin(x_lay_forward_cin, y_lay_forward_cin, 0); 01997 Hep3Vector Cell_z[2]; 01998 Cell_z[0]=((cell_IO[0].z()+length/2)/length)*(lay_backward_cin-lay_forward_cin)+lay_forward_cin; 01999 double phi_cin_z = Cell_z[0].phi(); 02000 phi_cin_z = fmod( phi_cin_z+4*M_PI,2*M_PI ); 02001 double x_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().x()*0.1; 02002 double y_lay_backward_cout = geosvc->Wire(layer, cid_out)->Backward().y()*0.1; 02003 double x_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().x()*0.1; 02004 double y_lay_forward_cout = geosvc->Wire(layer, cid_out)->Forward().y()*0.1; 02005 Hep3Vector lay_backward_cout(x_lay_backward_cout, y_lay_backward_cout, 0); 02006 Hep3Vector lay_forward_cout(x_lay_forward_cout, y_lay_forward_cout, 0); 02007 Cell_z[1]=((cell_IO[1].z()+length/2)/length)*(lay_backward_cout-lay_forward_cout)+lay_forward_cout; 02008 double phi_cout_z = Cell_z[1].phi(); 02009 phi_cout_z = fmod( phi_cout_z+4*M_PI,2*M_PI ); 02010 02011 phi_midin_z = phi_cin_z-phibin*0.5; 02012 phi_midout_z = phi_cout_z+phibin*0.5; 02013 phi_midin_z = fmod( phi_midin_z+4*M_PI,2*M_PI ); 02014 phi_midout_z = fmod( phi_midout_z+4*M_PI,2*M_PI ); 02015 phi1_z = phi_midout_z-phi_out; 02016 phi2_z = phi_in-phi_midin_z; 02017 phi1_z = fmod(phi1_z+2.0*M_PI,2.0*M_PI); 02018 phi2_z = fmod(phi2_z+2.0*M_PI,2.0*M_PI); 02019 gap_z = phi1_z+phi2_z+(cid_in-cid_out-1)*phibin; 02020 gap_z = fmod(gap_z+2.0*M_PI,2.0*M_PI); 02021 sampl[cid_in]=phi2_z/gap_z*ch_ltrk; 02022 sampl[cid_out]=phi1_z/gap_z*ch_ltrk; 02023 for( int jj = cid_out+1; jj<cid_in; jj++) { 02024 sampl[jj]=phibin/gap_z*ch_ltrk; 02025 } 02026 02027 //normal case 02028 /*phi_midin = phioff[0]+phibin*cid_in-phibin*0.5; 02029 phi_midout = phioff[1]+phibin*cid_out+phibin*0.5; 02030 phi_midin = fmod( phi_midin+4*M_PI,2*M_PI ); 02031 phi_midout = fmod( phi_midout+4*M_PI,2*M_PI ); 02032 phi1 = phi_midout-phi_out; 02033 phi2 = phi_in-phi_midin; 02034 phi1 = fmod(phi1+2.0*M_PI,2.0*M_PI); 02035 phi2 = fmod(phi2+2.0*M_PI,2.0*M_PI); 02036 gap = phi1+phi2+(cid_in-cid_out-1)*phibin; 02037 gap = fmod(gap+2.0*M_PI,2.0*M_PI); 02038 sampl[cid_in]=phi2/gap*ch_ltrk; 02039 sampl[cid_out]=phi1/gap*ch_ltrk; 02040 for( int jj = cid_out+1; jj<cid_in; jj++) { 02041 sampl[jj]=phibin/gap*ch_ltrk; 02042 }*/ 02043 } 02044 } 02045 02046 #ifdef YDEBUG 02047 if(sampl[cellid]<0.0) { 02048 if(cid_in!=cid_out) cout<<"?????????"<<endl; 02049 cout<< "layerId " << layer <<" cell id "<< cellid <<" shift " << shift 02050 << " cid_in " << cid_in << " cid_out " << cid_out << endl; 02051 02052 cout <<" phi_in " << phi_in <<" phi_midin " << phi_midin << " phi_out " 02053 << phi_out << " phi_midout " << phi_midout <<endl; 02054 cout<<"total sampl " << ch_ltrk << " gap "<< gap << " phi1 " 02055 << phi1 << " phi2 " << phi2 << " phibin " << phibin << endl; 02056 02057 02058 for(int l=0; l<ncell; l++) { 02059 if(sampl[l]>=0.0) 02060 cout<<"picked cellid " << l << " sampling length "<< sampl[l]<< endl; 02061 } 02062 } 02063 #endif 02064 return sampl[cellid]; 02065 }
|
|
|
|
00051 { 00052 //cout<<"DedxCorrecSvc::queryInterface"<<endl; 00053 if( IID_IDedxCorrecSvc.versionMatch(riid) ){ 00054 *ppvInterface = static_cast<IDedxCorrecSvc*> (this); 00055 } else{ 00056 return Service::queryInterface(riid, ppvInterface); 00057 } 00058 return StatusCode::SUCCESS; 00059 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00601 { 00602 // cout<<"DedxCorrecSvc::RungCorrec"<<endl; 00603 double dedx_rung; 00604 int run_ture =0; 00605 //cout<<"N_run : "<<N_run<<" runNO : "<<runNO<<endl; 00606 00607 for(int j=0;j<N_run;j++) { 00608 if(runNO == m_rung[2][j]) { 00609 dedx_rung = dedx/m_rung[0][j]; 00610 run_ture =1; 00611 return dedx_rung; 00612 } 00613 } 00614 if(run_ture ==0) return dedx; 00615 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00486 { 00487 //cout<<"DedxCorrecSvc::SaturCorrec"<<endl; 00488 double dedx_ggs; 00489 //cout<<"costheta vaule is = "<<costheta<<endl; 00490 costheta = fabs(costheta); 00491 if(m_par_flag == 1) { 00492 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*costheta + 00493 m_ggs[2][layer]*pow(costheta,2) + m_ggs[3][layer]*pow(costheta,3); 00494 } else if(m_par_flag == 0) { 00495 dedx_ggs = m_ggs[0][layer] + m_ggs[1][layer]*T1(costheta) + 00496 m_ggs[2][layer]*T2(costheta) + m_ggs[3][layer]*T3(costheta); 00497 } 00498 //cout<<"dedx_ggs is = "<<dedx_ggs<<endl; 00499 dedx_ggs = dedx/dedx_ggs; 00500 if(dedx_ggs>0.0) return dedx_ggs; 00501 else return dedx; 00502 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 01051 { 01052 // cout<<"DedxCorrecSvc::set_flag"<<endl; 01053 cout<<"calib_F is = "<<calib_F<<endl; 01054 m_calib_flag = calib_F; 01055 m_rung_flag = ( calib_F & 1 )? 1 : 0; 01056 m_wireg_flag = ( calib_F & 2 )? 1 : 0; 01057 m_dcasin_flag = ( calib_F & 4 )? 1 : 0; 01058 m_ggs_flag = ( calib_F & 8 )? 1 : 0; 01059 //m_ddg_flag = ( calib_F & 8 )? 1 : 0; 01060 m_t0_flag = ( calib_F & 16 )? 1 : 0; 01061 m_sat_flag = ( calib_F & 32 )? 1 : 0; 01062 m_layerg_flag = ( calib_F & 64 )? 1 : 0; 01063 //m_dcasin_flag = ( calib_F & 128 )? 1 : 0; 01064 m_dip_flag = ( calib_F & 256 )? 1 : 0; 01065 m_mdcg_flag = ( calib_F & 512 )? 1 : 0; 01066 }
|
|
00066 {return (x*x);}
|
|
00066 {return (x*x);}
|
|
|
|
|
|
|
|
|
|
00068 {return (x*x*x*x);}
|
|
00068 {return (x*x*x*x);}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00847 { 00848 // cout<<"DedxCorrecSvc::StandardCorrec"<<endl; 00849 int layid = MdcID::layer(mdcid); 00850 int localwid = MdcID::wire(mdcid); 00851 int w0id = geosvc->Layer(layid)->Wirst(); 00852 int wid = w0id + localwid; 00853 double pathl = PathL(ntpFlag, hel, layid, localwid, z); 00854 //pathl= PathLCosmic(hel, layid, localwid, z, sigmaz ); 00855 double ex = adc; 00856 if( runNO<0 ) return ex; 00857 ex = ex*1.5/pathl; 00858 //double ex = adc/pathl; 00859 //if( runNO<0 ) return ex; 00860 if( ntpFlag ==0) return ex; 00861 //double ex = adc/pathl; 00862 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0 00863 && m_layerg_flag == 0 && m_zdep_flag == 0 && 00864 m_ggs_flag == 0) return ex; 00865 00866 if(m_rung_flag) { 00867 ex = RungCorrec( runNO, ex ); 00868 } 00869 00870 if( m_wireg_flag) { 00871 ex = WireGainCorrec(wid, ex); 00872 } 00873 00874 if( m_dcasin_flag) { 00875 if(runFlag <3) { 00876 ex = DriftDistCorrec( layid, dd, ex ); 00877 } 00878 else if(runFlag >=3){ ex = DocaSinCorrec(layid, dd, eangle, ex);} 00879 } 00880 00881 if( m_ddg_flag) { 00882 ex = DriftDistCorrec( layid, dd, ex ); 00883 } 00884 00885 if(m_ggs_flag) { 00886 if(runFlag <3 ){ 00887 ex = SaturCorrec( layid, costheta, ex ); 00888 } 00889 else if(runFlag >=3){ ex = CosthetaCorrec( costheta, ex );} 00890 //ex = CosthetaCorrec( costheta, ex ); 00891 } 00892 00893 00894 if( m_sat_flag) { 00895 ex = HadronCorrec( costheta, ex ); 00896 } 00897 00898 //if( m_enta_flag) { 00899 //ex = EntaCorrec( layid, sinenta, ex ); 00900 //} 00901 00902 if( m_zdep_flag) { 00903 ex = ZdepCorrec( layid, z, ex ); 00904 } 00905 00906 if( m_layerg_flag) { 00907 ex = LayerGainCorrec( layid, ex ); 00908 } 00909 00910 if( m_dip_flag){ 00911 ex = DipAngCorrec(layid, costheta, ex); 00912 } 00913 00914 if( m_mdcg_flag) { 00915 ex = GlobalCorrec( ex ); 00916 } 00917 return ex; 00918 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00921 { 00922 // cout<<"DedxCorrecSvc::StandardHitCorrec"<<endl; 00923 int layid = MdcID::layer(mdcid); 00924 int localwid = MdcID::wire(mdcid); 00925 int w0id = geosvc->Layer(layid)->Wirst(); 00926 int wid = w0id + localwid; 00927 double pathl = PathL(ntpFlag, hel, layid, localwid, z); 00928 //pathl= PathLCosmic(hel, layid, localwid, z, sigmaz ); 00929 double ex = adc; 00930 if( runNO<0 ) return ex; 00931 ex = ex*1.5/pathl; 00932 //if( runNO<0 ) return ex; 00933 //double ex = adc/pathl; 00934 if( ntpFlag ==0) return ex; 00935 //if(ntpFlag>0) cout<<"dE/dx value after path correction: "<<ex<<endl; 00936 //double ex = adc/pathl; 00937 //cout<<m_rung_flag << " "<<m_wireg_flag<<" "<<m_ddg_flag<<" "<<m_ggs_flag<<endl; 00938 if( m_rung_flag == 0 && m_wireg_flag == 0 && m_ddg_flag == 0 00939 && m_layerg_flag == 0 && m_zdep_flag == 0 && 00940 m_ggs_flag == 0) return ex; 00941 if(m_rung_flag) { 00942 ex = RungCorrec( runNO, ex ); 00943 } 00944 //if(ntpFlag>0) cout<<"dE/dx value after run by run correction: "<<ex<<endl; 00945 00946 if( m_wireg_flag) { 00947 ex = WireGainCorrec(wid, ex); 00948 } 00949 //if(ntpFlag>0) cout<<"dE/dx value after wire gain correction: "<<ex<<endl; 00950 00951 if( m_dcasin_flag) { 00952 if(runFlag <3) { 00953 ex = DriftDistCorrec( layid, dd, ex ); 00954 } 00955 else if(runFlag >=3){ 00956 //cout<<layid<<" "<<dd<<" "<<eangle<<" "<<ex<<endl; 00957 ex = DocaSinCorrec(layid, dd, eangle, ex); 00958 } 00959 } 00960 00961 if( m_ddg_flag) { 00962 ex = DriftDistCorrec( layid, dd, ex ); 00963 } 00964 //if(ntpFlag>0) cout<<"dE/dx value after ddg correction: "<<ex<<endl; 00965 00966 if(m_ggs_flag) { 00967 if(runFlag <3 ){ 00968 ex = SaturCorrec( layid, costheta, ex ); 00969 } 00970 else if(runFlag >=3){ ex = ex;} 00971 //ex = CosthetaCorrec( costheta, ex ); 00972 } 00973 //if(ntpFlag>0) cout<<"dE/dx value after costheta correction: "<<ex<<endl; 00974 return ex; 00975 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00979 { 00980 // cout<<"DedxCorrecSvc::StandardTrackCorrec"<<endl; 00981 if( runNO<0 ) return ex; 00982 if( ntpFlag ==0) return ex; 00983 if( runFlag <3) return ex; 00984 if( calib_rec_Flag ==1){ 00985 ex = CosthetaCorrec( costheta, ex ); 00986 return ex; 00987 } 00988 00989 //if(ntpFlag>0) cout<<"trcak value before costheta correction: "<<ex<<endl; 00990 if( m_ggs_flag) { 00991 ex = CosthetaCorrec( costheta, ex ); 00992 } 00993 //if(ntpFlag>0) cout<<"trcak value after costheta correction: "<<ex<<endl; 00994 if(m_t0_flag){ 00995 if(runFlag==3) {ex= T0Correc(t0, ex);} 00996 else if(runFlag>3) {ex=ex;} 00997 } 00998 if( m_sat_flag) { 00999 ex = HadronCorrec( costheta, ex ); 01000 } 01001 //if(ntpFlag>0) cout<<"trcak value after all correction: "<<ex<<endl; 01002 01003 if( calib_rec_Flag ==2){ 01004 double rungain_temp =RungCorrec( runNO, ex )/ex; 01005 ex = ex/rungain_temp; 01006 } 01007 01008 //if(ntpFlag>0) cout<<"trcak value no run gain correction: "<<ex<<endl; 01009 return ex; 01010 01011 }
|
|
00077 {return 1;}
|
|
00077 {return 1;}
|
|
|
|
00548 { 00549 // cout<<"DedxCorrecSvc::T0Correc"<<endl; 00550 double dedx_t0; 00551 const int nbin=t0_k.size()+1; 00552 if(nbin!=35) 00553 cout<<"there should be 35 bins for t0 correction, check it!"<<endl; 00554 00555 int n=0; 00556 if(t0<575) 00557 n=0; 00558 else if(t0<610) 00559 n=1; 00560 else if(t0>=610 && t0<1190){ 00561 n=(int)( (t0-610.0)/20.0 ) + 2; 00562 } 00563 else if(t0>=1190 && t0<1225) 00564 n=31; 00565 else if(t0>=1225 && t0<1275) 00566 n=32; 00567 else if(t0>=1275) 00568 n=33; 00569 00570 dedx_t0=t0_k[n]*t0+t0_b[n]; 00571 00572 if(dedx_t0>0){ 00573 dedx_t0 = dedx/dedx_t0; 00574 return dedx_t0; 00575 } 00576 else return dedx; 00577 }
|
|
00078 {return x;}
|
|
00078 {return x;}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
00079 {return (2*sq(x)-1);}
|
|
00079 {return (2*sq(x)-1);}
|
|
|
|
|
|
00080 {return (4*cub(x)-3*x);}
|
|
00080 {return (4*cub(x)-3*x);}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Implements IDedxCorrecSvc. |
|
dEdx calib. for each track Implements IDedxCorrecSvc. 00830 { 00831 // cout<<"DedxCorrecSvc::TrkCorrec"<<endl; 00832 if( m_mdcg_flag == 0 ) return dedx; 00834 double calib_ex = dedx; 00835 00836 double run_const = m_dedx_gain; 00837 if( run_const > 0 && m_mdcg_flag != 0 ) calib_ex /= run_const; 00838 00839 return calib_ex; 00840 00841 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00474 { 00475 if( m_wire_g[wireid] > 0 ){ 00476 double ch_dedx = (ex/m_wire_g[wireid])*m_valid[wireid]; 00477 return ch_dedx; 00478 } 00479 else if( m_wire_g[wireid] == 0 ){ 00480 return ex; 00481 } 00482 else return 0; 00483 }
|
|
Implements IDedxCorrecSvc. |
|
Implements IDedxCorrecSvc. 00659 { 00660 // cout<<"DedxCorrecSvc::ZdepCorrec"<<endl; 00661 double dedx_zdep; 00662 if(m_par_flag == 1) { 00663 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*z + 00664 m_zdep[2][layer]*z*z + m_zdep[3][layer]*pow(z,3); 00665 } else if(m_par_flag == 0) { 00666 dedx_zdep = m_zdep[0][layer] + m_zdep[1][layer]*T1(z) + 00667 m_zdep[2][layer]*T2(z) + m_zdep[3][layer]*T3(z); 00668 } 00669 dedx_zdep = dedx/dedx_zdep; 00670 if (dedx_zdep>0.0) return dedx_zdep; 00671 else return dedx; 00672 00673 //return dedx; 00674 }
|
|
|
|
|
|
|
|
|
|
|
|
mdc geometry
|
|
mdc geometry
|
|
|
|
|
|
|
|
correction : doca and sinenta correction ( calib_F & 128 )? ON:OFF
|
|
|
|
calibration flag : drift distance ( calib_F & 8 ) ? ON : OFF
|
|
|
|
|
|
|
|
correction : Qsaturation correction ( calib_F & 256 )? ON:OFF
|
|
|
|
|
|
calibration flag : entrance angle depsndence ( calib_F & 16 )? ON : OFF
|
|
|
|
|
|
calibration flag : electron dip angle correction ( calib_F & 4 )? ON : OFF
|
|
control flags
|
|
|
|
|
|
correction : mdc layer gain ( calib_F & 64 )? ON:OFF
|
|
correction : mdc gain (run by run) ( calib_F & 512 )? ON:OFF
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
calibration flag : run by run correction ( calib_F & 1 )? ON : OFF
|
|
calibration flag : hadron saturation correction ( calib_F & 32 ) ? ON : OFF
|
|
calibration flag : t0 correction ( calib_F & 16 ) ? ON : OFF
|
|
|
|
|
|
|
|
|
|
calibration flag : wire by wire correction ( calib_F & 2 )? ON : OFF
|
|
|
|
calibration flag : z depsndence ( calib_F & 32 )? ON : OFF
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|