#include <SD0Tag.h>
Public Member Functions | |
SD0Tag (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
Private Attributes | |
int | m_MC_sample |
int | PID_flag |
int | Seperate_Charge |
int | Charge_default |
int | p_run [3467] |
double | p_Ecm [3467] |
NTuple::Tuple * | m_tuple1 |
NTuple::Item< int > | m_tagmode |
NTuple::Item< double > | m_mass_bc |
NTuple::Item< double > | m_delE_tag |
NTuple::Item< double > | m_EGam_max_0 |
NTuple::Item< int > | m_nGood |
NTuple::Item< int > | m_nGam |
NTuple::Item< int > | m_runNo |
NTuple::Item< int > | m_event |
NTuple::Item< int > | m_cosmic_ok |
Definition at line 13 of file SD0Tag.h.
SD0Tag::SD0Tag | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 63 of file SD0Tag.cxx.
References Charge_default, m_MC_sample, PID_flag, and Seperate_Charge.
00063 : 00064 Algorithm(name, pSvcLocator) { 00065 //Declare the properties 00066 declareProperty("PID_FLAG", PID_flag = 1); 00067 declareProperty("MC_sample", m_MC_sample = 1); 00068 00069 declareProperty("m_Seperate_Charge", Seperate_Charge = 1); 00070 declareProperty("m_Charge_default", Charge_default = 1); 00071 00072 00073 }
StatusCode SD0Tag::execute | ( | ) |
Definition at line 127 of file SD0Tag.cxx.
References VFHelix::a(), abs, DstEmcShower::cellId(), Charge_default, cos(), DTagTool::cosmicandleptonVeto(), ebeam, DstExtTrack::emcPosition(), DstExtTrack::emcVolumeNumber(), DstEmcShower::energy(), DstMdcTrack::err(), EventModel::EvtRec::EvtRecEvent, EventModel::EvtRec::EvtRecTrackCol, Sing::GetdelE_tag(), Sing::Getmass_bc(), Sing::Getoktg(), Sing::Gettagmd(), Sing::Gettagtrk1(), DstMdcTrack::helix(), genRecEmupikp::i, IVertexDbSvc::isVertexValid(), ganga-rec::j, m_cosmic_ok, m_delE_tag, m_EGam_max_0, m_event, m_mass_bc, m_nGam, m_nGood, m_runNo, m_tagmode, m_tuple1, Sing::Mdset(), msgSvc(), nD0, p_Ecm, p_run, pi, PID_flag, VFHelix::pivot(), IVertexDbSvc::PrimaryVertex(), runNo, Seperate_Charge, IVertexDbSvc::SigmaPrimaryVertex(), DstEmcShower::theta(), DstMdcTrack::theta(), DstEmcShower::time(), DstEmcShower::x(), DstEmcShower::y(), and DstEmcShower::z().
00127 { 00128 00129 MsgStream log(msgSvc(), name()); 00130 00131 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00132 int runNo=eventHeader->runNumber(); 00133 int event=eventHeader->eventNumber(); 00134 00135 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent); 00136 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol); 00137 //cout<<"//////////////////////////////////////"<<endl; 00138 nD0++; 00139 00140 // 00141 // The default energy is 3.773 GeV for psi(3770) data. 00142 // Alayizer can add calibrated energy here. 00143 double Ebeam = 3.773/2.0; 00144 00145 if(runNo >=11414 && runNo <= 23500) { 00146 for(int i = 0; i < 3467; i++) { 00147 if(runNo == p_run[i]) {Ebeam = p_Ecm[i]/2.0;} 00148 } 00149 } 00150 00151 if(runNo < 0) 00152 { 00153 int irun = abs(runNo); 00154 // Ecmset* ECMSET = Ecmset::instance(); 00155 // ECMSET->EcmSet(irun); 00156 // Ebeam = ECMSET->ECM()/2.0; 00157 } 00158 00159 if(nD0%1000==0) cout<<"SD0TagAlg-13-11-15pp = "<<nD0<<" "<<runNo<<" "<<event<<" "<<Ebeam*2<<endl; 00160 00161 Hep3Vector xorigin(0,0,0); 00162 IVertexDbSvc* vtxsvc; 00163 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc); 00164 if(vtxsvc->isVertexValid()){ 00165 double* dbv = vtxsvc->PrimaryVertex(); 00166 double* vv = vtxsvc->SigmaPrimaryVertex(); 00167 xorigin.setX(dbv[0]); 00168 xorigin.setY(dbv[1]); 00169 xorigin.setZ(dbv[2]); 00170 } 00171 00173 Vint iCharge_good; iCharge_good.clear(); 00174 for(int i = 0; i < evtRecEvent->totalCharged(); i++){ 00175 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i; 00176 if(!(*itTrk)->isMdcTrackValid()) continue; 00177 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack(); 00178 //%%%%%%%%%%%%%%%%%%%%%%%% 00179 HepVector a = mdcTrk->helix(); 00180 HepSymMatrix Ea = mdcTrk->err(); 00181 HepPoint3D point0(0.,0.,0.); // the initial point for MDC recosntruction 00182 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]); 00183 VFHelix helixip(point0,a,Ea); 00184 helixip.pivot(IP); 00185 HepVector vecipa = helixip.a(); 00186 double Rvxy0=fabs(vecipa[0]); //the nearest distance to IP in xy plane 00187 double Rvz0=vecipa[3]; //the nearest distance to IP in z direction 00188 double costheta = cos(mdcTrk->theta()); 00189 //%%%%%%%%%%%%%%%%%%%%%%%% 00190 if(fabs(Rvxy0) >= 1.0) continue; 00191 if(fabs(Rvz0) >= 10.0) continue; 00192 if(fabs(costheta) >= 0.930) continue; 00193 iCharge_good.push_back(i); 00194 } 00195 00197 Vint iGam; iGam.clear(); 00198 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) { 00199 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i; 00200 if(!(*itTrk)->isEmcShowerValid()) continue; 00201 RecEmcShower *emcTrk = (*itTrk)->emcShower(); 00202 // 00203 // We here remove the hot channels of EMC 00204 // 00205 if(abs(runNo)>=11615&&abs(runNo)<=11617&&emcTrk->cellId()==805438481) continue; 00206 if(abs(runNo)>=11615&&abs(runNo)<=11655&&emcTrk->cellId()==805376008) continue; 00207 if(abs(runNo)>=13629&&abs(runNo)<=13669&&emcTrk->cellId()==805374754) continue; 00208 if(abs(runNo)>=21190&&abs(runNo)<=21231&&emcTrk->cellId()==805375262) continue; 00209 00210 int emctdc = emcTrk->time(); 00211 if(emctdc>14||emctdc<0) continue; 00212 00213 double eraw = emcTrk->energy(); 00214 double theta = emcTrk->theta(); 00215 int Module = 0; 00216 if(fabs(cos(theta)) < 0.80 && eraw > 0.025){ Module = 1; } 00217 if(fabs(cos(theta)) > 0.86 && fabs(cos(theta)) < 0.92 && eraw > 0.05) { Module = 2; } 00218 if(Module == 0) continue; 00219 iGam.push_back((*itTrk)->trackId()); 00220 } 00222 DTagTool trk; 00223 bool cosmic_ok = trk.cosmicandleptonVeto(); 00225 00226 int ntk; 00227 double tagmass,ebeam,tagmode,ksmass,dlte; 00228 00229 Vint tagtrk; tagtrk.clear(); 00230 Vint tagGam; tagGam.clear(); 00231 00232 HepLorentzVector tagp; 00233 00234 double mass_bc_temp, mass_kf_temp, kf_chi2_temp, mks_temp, mpi0_temp, ptag_temp; 00235 int Charge_candidate_D = 0; 00236 double EGam_max_0 = 0; 00237 00238 for(int i1 = 0; i1 < 2; i1++) { 00239 if(Seperate_Charge == 2 ) { Charge_candidate_D = Charge_default; i1 = 2;} 00240 if(Seperate_Charge == 1 ) { Charge_candidate_D = pow(-1,i1); } 00241 if(Seperate_Charge == 0 ) { Charge_candidate_D = 0; i1 = 2; } 00242 00243 for(int i = 0; i < 15; i++) { 00244 EGam_max_0 = 0; 00245 int mdslct=pow(2.,i); 00246 Sing sing; 00247 sing.Mdset(event,evtRecTrkCol,iCharge_good,iGam,mdslct,Ebeam, PID_flag,Charge_candidate_D); 00248 bool oktag=sing.Getoktg(); 00249 if(oktag) { 00250 // 00251 // Here analysizer should pick up variables from anti-D0 tags 00252 // to do physics analysis 00253 // 00254 tagtrk=sing.Gettagtrk1(); 00255 tagmode = sing.Gettagmd(); 00256 //========================================================================== 00257 if((abs(tagmode) == 11 || abs(tagmode) == 15 || abs(tagmode) == 19)) { 00258 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00259 for(int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++) { 00260 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i; 00261 int itrk = (*itTrk)->trackId(); 00262 if(!(*itTrk)->isEmcShowerValid()) continue; 00263 RecEmcShower *emcTrk = (*itTrk)->emcShower(); 00264 00265 int emctdc = emcTrk->time(); 00266 if(emctdc>14||emctdc<0) continue; 00267 00268 Hep3Vector emcpos(emcTrk->x(), emcTrk->y(), emcTrk->z()); 00269 double dang = 200.; 00270 for(int j = 0; j < tagtrk.size(); j++) { 00271 EvtRecTrackIterator jtTrk = evtRecTrkCol->begin() + tagtrk[j]; 00272 if(!(*jtTrk)->isExtTrackValid()) continue; 00273 RecExtTrack *extTrk = (*jtTrk)->extTrack(); 00274 if(extTrk->emcVolumeNumber() == -1) continue; 00275 Hep3Vector extpos = extTrk->emcPosition(); 00276 double angd = extpos.angle(emcpos); 00277 if(angd < dang) dang = angd; 00278 } 00279 dang = dang * 180 / (CLHEP::pi); 00280 if(fabs(dang)<20.) continue; 00281 00282 double eraw = emcTrk->energy(); 00283 if(eraw > EGam_max_0) EGam_max_0 = eraw; 00284 } 00285 } 00286 //========================================================================== 00287 m_cosmic_ok = cosmic_ok; 00288 m_EGam_max_0 = EGam_max_0; 00289 m_nGood = iCharge_good.size(); 00290 m_nGam = iGam.size(); 00291 m_runNo = runNo; 00292 m_event = event; 00293 00294 m_tagmode = sing.Gettagmd(); 00295 m_mass_bc = sing.Getmass_bc(); 00296 m_delE_tag = sing.GetdelE_tag(); 00297 00298 m_tuple1->write(); 00299 00300 // Here one can do double tag analysis based on the variables 00301 // from single anti-D0 tag analysis. 00302 // To do so, analyzer should make his/her codes 00303 // 00304 } 00305 00306 } 00307 } // The end of loopping over mdslct 00308 00309 return StatusCode::SUCCESS; 00310 00311 }
StatusCode SD0Tag::finalize | ( | ) |
Definition at line 318 of file SD0Tag.cxx.
References Bes_Common::INFO, msgSvc(), n1, n2, ND0, nD0, and NDp.
00318 { 00319 MsgStream log(msgSvc(), name()); 00320 cout<<"nD0 = "<<nD0<<endl; 00321 cout<<"n1 = "<<n1<<endl; 00322 cout<<"n2 = "<<n2<<endl; 00323 cout<<"ND0 ==== "<<ND0<<endl; 00324 cout<<"NDp ==== "<<NDp<<endl; 00325 log << MSG::INFO << "in finalize()" << endmsg; 00326 return StatusCode::SUCCESS; 00327 }
StatusCode SD0Tag::initialize | ( | ) |
Definition at line 77 of file SD0Tag.cxx.
References calibUtil::ERROR, genRecEmupikp::i, Bes_Common::INFO, m_cosmic_ok, m_delE_tag, m_EGam_max_0, m_event, m_mass_bc, m_nGam, m_nGood, m_runNo, m_tagmode, m_tuple1, msgSvc(), ntupleSvc(), p_Ecm, and p_run.
00077 { 00078 00079 00080 MsgStream log(msgSvc(), name()); 00081 00082 log << MSG::INFO << "in initialize()" << endmsg; 00083 00084 StatusCode status; 00085 00086 ifstream readrunEcm; 00087 // 00088 // Read in the energy calibration constant for run-by-run data 00089 readrunEcm.open("../share/psipp_ecms_calrunbyrun_runno_3648runs.dat"); 00090 cout<<"readrunEcm = "<<readrunEcm.is_open()<<endl; 00091 for(int i = 0; i < 3467; i++) { 00092 readrunEcm>>p_run[i]; readrunEcm>>p_Ecm[i]; 00093 } 00094 00095 00096 NTuplePtr nt1(ntupleSvc(), "FILE1/tag"); 00097 if ( nt1 ) m_tuple1 = nt1; 00098 else { 00099 m_tuple1 = ntupleSvc()->book ("FILE1/tag", CLID_ColumnWiseTuple, "tag N-Tuple example"); 00100 if ( m_tuple1 ) { 00101 status = m_tuple1->addItem ("tagmode", m_tagmode); 00102 status = m_tuple1->addItem ("mass_bc", m_mass_bc); 00103 status = m_tuple1->addItem ("delE_tag", m_delE_tag); 00104 00105 status = m_tuple1->addItem ("cosmic_ok", m_cosmic_ok); 00106 status = m_tuple1->addItem ("EGam_max_0", m_EGam_max_0); 00107 status = m_tuple1->addItem ("nGood", m_nGood); 00108 status = m_tuple1->addItem ("nGam", m_nGam); 00109 status = m_tuple1->addItem ("runNo", m_runNo); 00110 status = m_tuple1->addItem ("event", m_event); 00111 } 00112 else { 00113 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) <<endmsg; 00114 return StatusCode::FAILURE; 00115 } 00116 } 00117 00118 log << MSG::INFO << "successfully return from initialize()" <<endmsg; 00119 return StatusCode::SUCCESS; 00120 }
int SD0Tag::Charge_default [private] |
NTuple::Item<int> SD0Tag::m_cosmic_ok [private] |
NTuple::Item<double> SD0Tag::m_delE_tag [private] |
NTuple::Item<double> SD0Tag::m_EGam_max_0 [private] |
NTuple::Item<int> SD0Tag::m_event [private] |
NTuple::Item<double> SD0Tag::m_mass_bc [private] |
int SD0Tag::m_MC_sample [private] |
NTuple::Item<int> SD0Tag::m_nGam [private] |
NTuple::Item<int> SD0Tag::m_nGood [private] |
NTuple::Item<int> SD0Tag::m_runNo [private] |
NTuple::Item<int> SD0Tag::m_tagmode [private] |
NTuple::Tuple* SD0Tag::m_tuple1 [private] |
double SD0Tag::p_Ecm[3467] [private] |
int SD0Tag::p_run[3467] [private] |
int SD0Tag::PID_flag [private] |
int SD0Tag::Seperate_Charge [private] |