1#include "GaudiKernel/IDataProviderSvc.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/SmartDataPtr.h"
16#include "TofCorrection.C"
19 : Service(name, svcLoc)
21 declareProperty(
"mindedxchi", m_dedxminchi=4);
22 declareProperty(
"mintofchi", m_tofminchi=2);
23 declareProperty(
"tofcorrection", m_tofcorrec=
true);
32 MsgStream log(messageService(), name());
33 log << MSG::INFO <<
"@initialize()" << endreq;
35 StatusCode sc = Service::initialize();
37 sc = serviceLocator()->service(
"EventDataSvc", eventSvc_,
true);
45 filename = string(getenv(
"SIMPLEPIDSVCROOT"))+
"/share/eop.root";
46 TFile* file=
new TFile(filename.c_str(),
"read");
47 h_ebarlp=(TH1F*)file->Get(
"barlp");
48 h_ebarhp=(TH1F*)file->Get(
"barhp");
49 h_eendlp=(TH1F*)file->Get(
"endlp");
50 h_eendhp=(TH1F*)file->Get(
"endhp");
51 h_kbar=(TH1F*)file->Get(
"barkaon");
52 h_kend=(TH1F*)file->Get(
"endkaon");
53 h_pibar=(TH1F*)file->Get(
"barpion");
54 h_piend=(TH1F*)file->Get(
"endpion");
61 MsgStream log(messageService(), name());
62 log << MSG::INFO <<
"@initialize()" << endreq;
64 StatusCode sc = Service::finalize();
85 return Service::queryInterface(riid, ppvIF);
88 return StatusCode::SUCCESS;
94 filename = string(getenv(
"SIMPLEPIDSVCROOT"))+
"/share/constant/const.txt";
100 std::ifstream fs(filename.c_str());
102 int run_begin, run_end, method;
103 bool isdata=
false, ismc=
false;
104 while(fs.good() && !fs.eof() ){
114 if(tag==99 && run>=run_begin && run<=run_end){
118 if(tag==-99 && run>=run_begin && run<=run_end){
126 m_datatof.push_back(
num);
128 m_mctof.push_back(
num);
136 SmartDataPtr<Event::EventHeader> eventHeaderpid(eventSvc_,
"/Event/EventHeader");
137 m_run=eventHeaderpid->runNumber();
143 if(m_datatof.size()==0 || m_mctof.size()==0)
151 for(
int pid=0; pid<5; ++pid){
174 double kappa=zhelix[2];
175 double theta=CLHEP::pi/2.0-atan(zhelix[4]);
176 m_p[pid]=1.0/fabs(kappa)/
sin(theta);
177 m_costh[pid]=
cos(theta);
181 for(
int i=0; i<5; ++i){
197 for(
int i=0; i<5; ++i){
200 m_tofscale1[i]=toflayer1scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_datatof);
201 m_tofscale2[i]=toflayer2scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_datatof);
204 m_tofscale1[i]=mctoflayer1scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_mctof);
205 m_tofscale2[i]=mctoflayer2scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_mctof);
221 m_eop=
energy/fabs(m_p[0]);
229 m_dedxchi[0]= dedxTrk->
chi(0);
230 m_dedxchi[1]= dedxTrk->
chi(1);
231 m_dedxchi[2]= dedxTrk->
chi(2);
232 m_dedxchi[3]= dedxTrk->
chi(3);
233 m_dedxchi[4]= dedxTrk->
chi(4);
276 SmartRefVector<RecTofTrack> tofTrkCol=track->
tofTrack();
277 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
279 for(;iter_tof!=tofTrkCol.end();iter_tof++){
281 status->
setStatus( (*iter_tof)->status() );
314 for(
int i=0; i<5; i++)
315 m_tofdt1[i]=(*iter_tof)->tof() - m_tofscale1[i]*(*iter_tof)->texp(i);
316 m_sigma1=(*iter_tof)->sigma(0);;
321 if(status->
layer()==1){
324 for(
int i=0; i<5; i++){
325 m_tofdt1[i]=(*iter_tof)->tof() - m_tofscale1[i]*(*iter_tof)->texp(i);
327 m_sigma1=(*iter_tof)->sigma(0);;
329 else if(status->
layer()==2){
332 for(
int i=0; i<5; i++){
333 m_tofdt2[i]=(*iter_tof)->tof() - m_tofscale2[i]*(*iter_tof)->texp(i);
335 m_sigma2=(*iter_tof)->sigma(0);;
375 for(
int i=0; i<5 ; ++i){
377 double dt1= m_tofdt1[i];
378 double dt2= m_tofdt2[i];
379 double sigma1= m_sigma1;
380 double sigma2= m_sigma2;
382 double sqrtcov=0.041;
384 sqrtcov=tofsqrcov(m_p[i]/m_mass[i],m_datatof);
386 sqrtcov=mctofsqrcov(m_p[i]/m_mass[i],m_mctof);
389 double weight1=(sigma2*sigma2-sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov);
390 double weight2= (sigma1*sigma1-sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov);
393 double sigma_comb=-99;
397 if( dt1!=-99 && dt2!=-99){
399 dt_comb= weight1*dt1+weight2*dt2;
400 sigma_comb= sqrt((sigma1*sigma1*sigma2*sigma2 - sqrtcov*sqrtcov*sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov));
402 else if(dt1!=-99 && dt2==-99){
406 else if ( dt1==-99 && dt2!=-99){
412 double sigmascale=1.0;
414 sigmascale=tofsigmascale(m_p[i]/m_mass[i],m_datatof);
416 sigmascale=mctofsigmascale(m_p[i]/m_mass[i],m_mctof);
419 m_tofchi[i]= dt_comb/sigma_comb/sigmascale;
435 for(
int i=0; i<5 ;i++){
437 if(!usededx && fabs(m_dedxchi[i])<m_dedxminchi)
439 if(!usetof && fabs(m_tofchi[i])<m_tofminchi)
447 for(
int i=0; i<5; i++)
453 for(
int i=0; i<5; i++)
462 for(
int i=0; i<5; ++i){
468 if( usededx && usetof ){
469 chi2= pow(m_dedxchi[i],2)+pow(m_tofchi[i],2);
472 else if( usededx && !usetof ){
473 chi2= pow(m_dedxchi[i],2);
477 else if( !usededx && usetof ){
478 chi2= pow(m_tofchi[i],2);
499 m_prob[i]=TMath::Prob(chi2, ndf);
511 if(m_prob[0]>0 && m_prob[0]>m_prob[2] && m_prob[0]>m_prob[3]){
517 if( fabs(m_costh[0])<0.7 && m_eop>0 && m_eop<0.8)
520 if( fabs(m_costh[0])>=0.7 && fabs(m_costh[0])<0.8 && m_eop>0 && m_eop<-7.5*fabs(m_costh[0])+6.05)
523 if( fabs(m_costh[0])>0.85 && m_eop>0 && m_eop<0.6)
620 if(m_prob[2]>=0.00 && m_prob[2]>=m_prob[3])
633 if(m_prob[3]>=0.00 && m_prob[3]>=m_prob[2]){
double sin(const BesAngle a)
double cos(const BesAngle a)
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
static void setPidType(PidType pidType)
bool isMdcKalTrackValid()
SmartRefVector< RecTofTrack > tofTrack()
RecEmcShower * emcShower()
RecMdcKalTrack * mdcKalTrack()
static const InterfaceID & interfaceID()
const HepVector & getZHelix() const
HepVector & getZHelixMu()
SimplePIDSvc(const std::string &name, ISvcLocator *svcLoc)
bool iselectron(bool eop=false)
void preparePID(EvtRecTrack *track)
virtual StatusCode initialize()
virtual StatusCode finalize()
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvIF)
unsigned int layer() const
void setStatus(unsigned int status)