Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

DedxCorrecSvc Class Reference

#include <DedxCorrecSvc.h>

Inheritance diagram for DedxCorrecSvc:

IDedxCorrecSvc IDedxCorrecSvc List of all members.

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]
IMdcGeomSvcgeosvc
 mdc geometry
IMdcGeomSvcgeosvc
 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
IMagneticFieldSvcm_pIMF
IMagneticFieldSvcm_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

Constructor & Destructor Documentation

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

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   }

DedxCorrecSvc::~DedxCorrecSvc  ) 
 

00048                              {
00049 }

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

DedxCorrecSvc::~DedxCorrecSvc  ) 
 


Member Function Documentation

double DedxCorrecSvc::CellCorrec int  ser,
double  adc,
double  dd,
double  enta,
double  z,
double  costheta
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::CellCorrec int  ser,
double  adc,
double  dd,
double  enta,
double  z,
double  costheta
const [virtual]
 

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 }

double DedxCorrecSvc::CosthetaCorrec double  costheta,
double  ex
const
 

double DedxCorrecSvc::CosthetaCorrec double  costheta,
double  ex
const
 

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 }

double DedxCorrecSvc::cub double &  x  )  const [inline, private]
 

00067 {return (x*x*x);}

double DedxCorrecSvc::cub double &  x  )  const [inline, private]
 

00067 {return (x*x*x);}

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

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

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 }

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

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

00754                                                                          {
00755 }

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

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

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 }

double DedxCorrecSvc::DriftDistCorrec int  layid,
double  ddrift,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::DriftDistCorrec int  layid,
double  ddrift,
double  ex
const [virtual]
 

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 }

double DedxCorrecSvc::EntaCorrec int  layid,
double  enta,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::EntaCorrec int  layid,
double  enta,
double  ex
const [virtual]
 

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 }

virtual StatusCode DedxCorrecSvc::finalize  )  [virtual]
 

StatusCode DedxCorrecSvc::finalize  )  [virtual]
 

00091                                   {
00092   // cout<<"DedxCorrecSvc::finalize()"<<endl;
00093   MsgStream log(messageService(), name());
00094   log << MSG::INFO << name() << "DedxCorrecSvc::finalize()" << endreq;
00095   return StatusCode::SUCCESS;
00096 }

double DedxCorrecSvc::GlobalCorrec double  dedx  )  const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::GlobalCorrec double  dedx  )  const [virtual]
 

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 }

double DedxCorrecSvc::HadronCorrec double  costheta,
double  dedx
const
 

double DedxCorrecSvc::HadronCorrec double  costheta,
double  dedx
const
 

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 }

void DedxCorrecSvc::handle const Incident &   ) 
 

void DedxCorrecSvc::handle const Incident &   ) 
 

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 }

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

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

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 }

void DedxCorrecSvc::init_param  )  [private]
 

void DedxCorrecSvc::init_param  )  [private]
 

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 }

void DedxCorrecSvc::init_param_svc  )  [private]
 

void DedxCorrecSvc::init_param_svc  )  [private]
 

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 }

virtual StatusCode DedxCorrecSvc::initialize  )  [virtual]
 

StatusCode DedxCorrecSvc::initialize  )  [virtual]
 

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 }

const InterfaceID& IDedxCorrecSvc::interfaceID  )  [inline, static, inherited]
 

00015 { return IID_IDedxCorrecSvc; }

const InterfaceID& IDedxCorrecSvc::interfaceID  )  [inline, static, inherited]
 

00015 { return IID_IDedxCorrecSvc; }

double DedxCorrecSvc::LayerCorrec int  layer,
double  z,
double  costheta,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::LayerCorrec int  layer,
double  z,
double  costheta,
double  ex
const [virtual]
 

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 }

double DedxCorrecSvc::LayerGainCorrec int  layid,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::LayerGainCorrec int  layid,
double  ex
const [virtual]
 

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 }

double DedxCorrecSvc::PathL int  ntpFlag,
const Dedx_Helix hel,
int  layer,
int  cellid,
double  z
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::PathL int  ntpFlag,
const Dedx_Helix hel,
int  layer,
int  cellid,
double  z
const [virtual]
 

***********-------------------***************------------------****************//

***********-------------------***************------------------****************//

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 }

virtual StatusCode DedxCorrecSvc::queryInterface const InterfaceID &  riid,
void **  ppvUnknown
[virtual]
 

StatusCode DedxCorrecSvc::queryInterface const InterfaceID &  riid,
void **  ppvUnknown
[virtual]
 

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 }

double DedxCorrecSvc::RungCorrec int  runNO,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::RungCorrec int  runNO,
double  ex
const [virtual]
 

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 }

double DedxCorrecSvc::SaturCorrec int  layid,
double  costheta,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::SaturCorrec int  layid,
double  costheta,
double  ex
const [virtual]
 

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   }

void DedxCorrecSvc::set_flag int  calib_F  )  [virtual]
 

Implements IDedxCorrecSvc.

void DedxCorrecSvc::set_flag int  calib_F  )  [virtual]
 

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 }

double DedxCorrecSvc::sq double &  x  )  const [inline, private]
 

00066 {return (x*x);}

double DedxCorrecSvc::sq double &  x  )  const [inline, private]
 

00066 {return (x*x);}

double DedxCorrecSvc::sq10 double &  x  )  const [inline, private]
 

00074 {return sq5(x)*sq5(x);}

double DedxCorrecSvc::sq10 double &  x  )  const [inline, private]
 

00074 {return sq5(x)*sq5(x);}

double DedxCorrecSvc::sq11 double &  x  )  const [inline, private]
 

00075 {return sq5(x)*sq6(x);} 

double DedxCorrecSvc::sq11 double &  x  )  const [inline, private]
 

00075 {return sq5(x)*sq6(x);} 

double DedxCorrecSvc::sq4 double &  x  )  const [inline, private]
 

00068 {return (x*x*x*x);}

double DedxCorrecSvc::sq4 double &  x  )  const [inline, private]
 

00068 {return (x*x*x*x);}

double DedxCorrecSvc::sq5 double &  x  )  const [inline, private]
 

00069 {return sq(x)*cub(x);}

double DedxCorrecSvc::sq5 double &  x  )  const [inline, private]
 

00069 {return sq(x)*cub(x);}

double DedxCorrecSvc::sq6 double &  x  )  const [inline, private]
 

00070 {return cub(x)*cub(x);}

double DedxCorrecSvc::sq6 double &  x  )  const [inline, private]
 

00070 {return cub(x)*cub(x);}

double DedxCorrecSvc::sq7 double &  x  )  const [inline, private]
 

00071 {return sq4(x)*cub(x);}

double DedxCorrecSvc::sq7 double &  x  )  const [inline, private]
 

00071 {return sq4(x)*cub(x);}

double DedxCorrecSvc::sq8 double &  x  )  const [inline, private]
 

00072 {return sq4(x)*sq4(x);}

double DedxCorrecSvc::sq8 double &  x  )  const [inline, private]
 

00072 {return sq4(x)*sq4(x);}

double DedxCorrecSvc::sq9 double &  x  )  const [inline, private]
 

00073 {return sq4(x)*sq5(x);}

double DedxCorrecSvc::sq9 double &  x  )  const [inline, private]
 

00073 {return sq4(x)*sq5(x);}

double DedxCorrecSvc::StandardCorrec int  typFlag,
int  ntpFlag,
int  runNO,
Dedx_Helix hel,
Identifier  mdcid,
double  adc,
double  dd,
double  enta,
double  z,
double  costheta
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::StandardCorrec int  typFlag,
int  ntpFlag,
int  runNO,
Dedx_Helix hel,
Identifier  mdcid,
double  adc,
double  dd,
double  enta,
double  z,
double  costheta
const [virtual]
 

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 }

double DedxCorrecSvc::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 [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::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 [virtual]
 

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 }

double DedxCorrecSvc::StandardTrackCorrec int  calib_rec_Flag,
int  typFlag,
int  ntpFlag,
int  runNO,
double  ex,
double  costheta,
double  t0
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::StandardTrackCorrec int  calib_rec_Flag,
int  typFlag,
int  ntpFlag,
int  runNO,
double  ex,
double  costheta,
double  t0
const [virtual]
 

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 }

double DedxCorrecSvc::T0  )  [inline, private]
 

00077 {return 1;}

double DedxCorrecSvc::T0  )  [inline, private]
 

00077 {return 1;}

double DedxCorrecSvc::T0Correc double  t0,
double  ex
const
 

double DedxCorrecSvc::T0Correc double  t0,
double  ex
const
 

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 }

double DedxCorrecSvc::T1 double &  x  )  const [inline, private]
 

00078 {return x;}

double DedxCorrecSvc::T1 double &  x  )  const [inline, private]
 

00078 {return x;}

double DedxCorrecSvc::T10 double &  x  )  const [inline, private]
 

00087 {return (512*sq10(x)-1280*sq8(x)+1120*sq6(x)-400*sq4(x)+50*sq(x)-1);}

double DedxCorrecSvc::T10 double &  x  )  const [inline, private]
 

00087 {return (512*sq10(x)-1280*sq8(x)+1120*sq6(x)-400*sq4(x)+50*sq(x)-1);}

double DedxCorrecSvc::T11 double &  x  )  const [inline, private]
 

00088 {return (1024*sq11(x)-2816*sq9(x)+2816*sq7(x)-1232*sq5(x)+220*cub(x)-11*x);}

double DedxCorrecSvc::T11 double &  x  )  const [inline, private]
 

00088 {return (1024*sq11(x)-2816*sq9(x)+2816*sq7(x)-1232*sq5(x)+220*cub(x)-11*x);}

double DedxCorrecSvc::T12 double &  x  )  const [inline, private]
 

00089 {return (2*x*T11(x)-T10(x));}

double DedxCorrecSvc::T12 double &  x  )  const [inline, private]
 

00089 {return (2*x*T11(x)-T10(x));}

double DedxCorrecSvc::T13 double &  x  )  const [inline, private]
 

00090 {return (2*x*T12(x)-T11(x));}

double DedxCorrecSvc::T13 double &  x  )  const [inline, private]
 

00090 {return (2*x*T12(x)-T11(x));}

double DedxCorrecSvc::T14 double &  x  )  const [inline, private]
 

00091 {return (2*x*T13(x)-T12(x));}

double DedxCorrecSvc::T14 double &  x  )  const [inline, private]
 

00091 {return (2*x*T13(x)-T12(x));}

double DedxCorrecSvc::T15 double &  x  )  const [inline, private]
 

00092 {return (2*x*T14(x)-T13(x));}

double DedxCorrecSvc::T15 double &  x  )  const [inline, private]
 

00092 {return (2*x*T14(x)-T13(x));}

double DedxCorrecSvc::T16 double &  x  )  const [inline, private]
 

00093 {return (2*x*T15(x)-T14(x));}

double DedxCorrecSvc::T16 double &  x  )  const [inline, private]
 

00093 {return (2*x*T15(x)-T14(x));}

double DedxCorrecSvc::T17 double &  x  )  const [inline, private]
 

00094 {return (2*x*T16(x)-T15(x));}

double DedxCorrecSvc::T17 double &  x  )  const [inline, private]
 

00094 {return (2*x*T16(x)-T15(x));}

double DedxCorrecSvc::T18 double &  x  )  const [inline, private]
 

00095 {return (2*x*T17(x)-T16(x));}

double DedxCorrecSvc::T18 double &  x  )  const [inline, private]
 

00095 {return (2*x*T17(x)-T16(x));}

double DedxCorrecSvc::T19 double &  x  )  const [inline, private]
 

00096 {return (2*x*T18(x)-T17(x));}

double DedxCorrecSvc::T19 double &  x  )  const [inline, private]
 

00096 {return (2*x*T18(x)-T17(x));}

double DedxCorrecSvc::T2 double &  x  )  const [inline, private]
 

00079 {return (2*sq(x)-1);}

double DedxCorrecSvc::T2 double &  x  )  const [inline, private]
 

00079 {return (2*sq(x)-1);}

double DedxCorrecSvc::T20 double &  x  )  const [inline, private]
 

00097 {return (2*x*T19(x)-T18(x));}      

double DedxCorrecSvc::T20 double &  x  )  const [inline, private]
 

00097 {return (2*x*T19(x)-T18(x));}      

double DedxCorrecSvc::T3 double &  x  )  const [inline, private]
 

00080 {return (4*cub(x)-3*x);}

double DedxCorrecSvc::T3 double &  x  )  const [inline, private]
 

00080 {return (4*cub(x)-3*x);}

double DedxCorrecSvc::T4 double &  x  )  const [inline, private]
 

00081 {return (8*sq4(x)-8*sq(x)+1);}

double DedxCorrecSvc::T4 double &  x  )  const [inline, private]
 

00081 {return (8*sq4(x)-8*sq(x)+1);}

double DedxCorrecSvc::T5 double &  x  )  const [inline, private]
 

00082 {return (16*sq5(x)-20*cub(x)+5*x);}

double DedxCorrecSvc::T5 double &  x  )  const [inline, private]
 

00082 {return (16*sq5(x)-20*cub(x)+5*x);}

double DedxCorrecSvc::T6 double &  x  )  const [inline, private]
 

00083 {return (32*sq6(x)-48*sq4(x)+18*sq(x)-1);}

double DedxCorrecSvc::T6 double &  x  )  const [inline, private]
 

00083 {return (32*sq6(x)-48*sq4(x)+18*sq(x)-1);}

double DedxCorrecSvc::T7 double &  x  )  const [inline, private]
 

00084 {return (64*sq7(x)-112*sq5(x)+56*cub(x)-7*x);}

double DedxCorrecSvc::T7 double &  x  )  const [inline, private]
 

00084 {return (64*sq7(x)-112*sq5(x)+56*cub(x)-7*x);}

double DedxCorrecSvc::T8 double &  x  )  const [inline, private]
 

00085 {return (128*sq8(x)-256*sq6(x)+160*sq4(x)-32*sq(x)+1);}

double DedxCorrecSvc::T8 double &  x  )  const [inline, private]
 

00085 {return (128*sq8(x)-256*sq6(x)+160*sq4(x)-32*sq(x)+1);}

double DedxCorrecSvc::T9 double &  x  )  const [inline, private]
 

00086 {return (256*sq9(x)-576*sq7(x)+432*sq5(x)-120*cub(x)+9*x);}

double DedxCorrecSvc::T9 double &  x  )  const [inline, private]
 

00086 {return (256*sq9(x)-576*sq7(x)+432*sq5(x)-120*cub(x)+9*x);}

double DedxCorrecSvc::TrkCorrec double  costheta,
double  dedx
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::TrkCorrec double  costheta,
double  dedx
const [virtual]
 

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 }

double DedxCorrecSvc::WireGainCorrec int  wireid,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::WireGainCorrec int  wireid,
double  ex
const [virtual]
 

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 }

double DedxCorrecSvc::ZdepCorrec int  layid,
double  zhit,
double  ex
const [virtual]
 

Implements IDedxCorrecSvc.

double DedxCorrecSvc::ZdepCorrec int  layid,
double  zhit,
double  ex
const [virtual]
 

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 }


Member Data Documentation

vector<double> DedxCorrecSvc::cos_b [private]
 

vector<double> DedxCorrecSvc::cos_b [private]
 

vector<double> DedxCorrecSvc::cos_k [private]
 

vector<double> DedxCorrecSvc::cos_k [private]
 

double DedxCorrecSvc::curve_par [private]
 

IMdcGeomSvc* DedxCorrecSvc::geosvc [private]
 

mdc geometry

IMdcGeomSvc* DedxCorrecSvc::geosvc [private]
 

mdc geometry

double DedxCorrecSvc::Iner_Stepdoca [private]
 

double DedxCorrecSvc::m_alpha [private]
 

int DedxCorrecSvc::m_calib_flag [private]
 

int DedxCorrecSvc::m_dcasin_flag [private]
 

correction : doca and sinenta correction ( calib_F & 128 )? ON:OFF

double DedxCorrecSvc::m_ddg [private]
 

int DedxCorrecSvc::m_ddg_flag [private]
 

calibration flag : drift distance ( calib_F & 8 ) ? ON : OFF

double DedxCorrecSvc::m_dedx_gain [private]
 

double DedxCorrecSvc::m_dedx_resol [private]
 

double DedxCorrecSvc::m_delta [private]
 

int DedxCorrecSvc::m_dip_flag [private]
 

correction : Qsaturation correction ( calib_F & 256 )? ON:OFF

double DedxCorrecSvc::m_docaeangle [private]
 

double DedxCorrecSvc::m_enta [private]
 

int DedxCorrecSvc::m_enta_flag [private]
 

calibration flag : entrance angle depsndence ( calib_F & 16 )? ON : OFF

double DedxCorrecSvc::m_gamma [private]
 

double DedxCorrecSvc::m_ggs [private]
 

int DedxCorrecSvc::m_ggs_flag [private]
 

calibration flag : electron dip angle correction ( calib_F & 4 )? ON : OFF

int DedxCorrecSvc::m_init [private]
 

control flags

bool DedxCorrecSvc::m_initConstFlg [private]
 

double DedxCorrecSvc::m_layer_g [private]
 

int DedxCorrecSvc::m_layerg_flag [private]
 

correction : mdc layer gain ( calib_F & 64 )? ON:OFF

int DedxCorrecSvc::m_mdcg_flag [private]
 

correction : mdc gain (run by run) ( calib_F & 512 )? ON:OFF

int DedxCorrecSvc::m_par_flag [private]
 

IMagneticFieldSvc* DedxCorrecSvc::m_pIMF [private]
 

IMagneticFieldSvc* DedxCorrecSvc::m_pIMF [private]
 

double DedxCorrecSvc::m_power [private]
 

double DedxCorrecSvc::m_ratio [private]
 

IntegerProperty DedxCorrecSvc::m_run [private]
 

double DedxCorrecSvc::m_rung [private]
 

int DedxCorrecSvc::m_rung_flag [private]
 

calibration flag : run by run correction ( calib_F & 1 )? ON : OFF

int DedxCorrecSvc::m_sat_flag [private]
 

calibration flag : hadron saturation correction ( calib_F & 32 ) ? ON : OFF

int DedxCorrecSvc::m_t0_flag [private]
 

calibration flag : t0 correction ( calib_F & 16 ) ? ON : OFF

double DedxCorrecSvc::m_valid [private]
 

vector<double> DedxCorrecSvc::m_venangle [private]
 

vector<double> DedxCorrecSvc::m_venangle [private]
 

double DedxCorrecSvc::m_wire_g [private]
 

int DedxCorrecSvc::m_wireg_flag [private]
 

calibration flag : wire by wire correction ( calib_F & 2 )? ON : OFF

double DedxCorrecSvc::m_zdep [private]
 

int DedxCorrecSvc::m_zdep_flag [private]
 

calibration flag : z depsndence ( calib_F & 32 )? ON : OFF

int DedxCorrecSvc::N_run [private]
 

double DedxCorrecSvc::Out_Stepdoca [private]
 

double DedxCorrecSvc::sigma_par [private]
 

vector<double> DedxCorrecSvc::t0_b [private]
 

vector<double> DedxCorrecSvc::t0_b [private]
 

vector<double> DedxCorrecSvc::t0_k [private]
 

vector<double> DedxCorrecSvc::t0_k [private]
 


The documentation for this class was generated from the following files:
Generated on Wed Feb 2 15:57:32 2011 for BOSS6.5.5 by  doxygen 1.3.9.1