5#include "GaudiKernel/MsgStream.h"
6#include "GaudiKernel/AlgFactory.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/SmartDataPtr.h"
9#include "GaudiKernel/IDataProviderSvc.h"
10#include "GaudiKernel/PropertyMgr.h"
12#include "EventModel/EventHeader.h"
13#include "EventModel/Event.h"
14#include "EventModel/EventModel.h"
16#include "EvtRecEvent/EvtRecEvent.h"
17#include "EvtRecEvent/EvtRecTrack.h"
18#include "EvtRecEvent/EvtRecDTag.h"
19#include "EvtRecEvent/EvtRecPi0.h"
20#include "EvtRecEvent/EvtRecEtaToGG.h"
21#include "EvtRecEvent/EvtRecVeeVertex.h"
23#include "BesDChain/CDChargedPionList.h"
24#include "BesDChain/CDChargedKaonList.h"
25#include "BesDChain/CDPhotonList.h"
26#include "BesDChain/CDPi0List.h"
27#include "BesDChain/CDEtaList.h"
28#include "BesDChain/CDKsList.h"
29#include "BesDChain/CDDecayList.h"
31#include "McTruth/McParticle.h"
32#include "ParticleID/ParticleID.h"
33#include "VertexFit/VertexFit.h"
34#include "VertexFit/SecondVertexFit.h"
37#include "DTagAlg/DsReconstruction.h"
38#include "DTagAlg/LocalKaonSelector.h"
39#include "DTagAlg/LocalPionSelector.h"
40#include "DTagAlg/LocalPhotonSelector.h"
41#include "DTagAlg/LocalKsSelector.h"
42#include "DTagAlg/LocalPi0Selector.h"
43#include "DTagAlg/LocalEtatoGGSelector.h"
44#include "DTagAlg/LocalEptoPiPiEtaSelector.h"
45#include "DTagAlg/LocalEptoRhoGamSelector.h"
46#include "DTagAlg/LocalRhotoPiPiSelector.h"
47#include "DTagAlg/DsSelector.h"
49#include "DTagAlg/utility.h"
55Algorithm(name, pSvcLocator) {
57 declareProperty(
"debug", m_debug =
false );
58 declareProperty(
"ReadBeamEFromDB", m_ReadBeamEFromDB =
false );
59 declareProperty(
"UseCalibBeamE", m_usecalibBeamE =
false );
60 declareProperty(
"UseVertexfit", m_usevertexfit =
false );
61 declareProperty(
"BeamE", m_beamE=2.015 );
62 declareProperty(
"DsList", m_decaylist =
"test.txt" );
67 MsgStream log(
msgSvc(), name());
68 log << MSG::INFO <<
"in initialize()" <<endreq;
76 return StatusCode::SUCCESS;
81 MsgStream log(
msgSvc(), name());
82 log << MSG::INFO <<
"in finalize()" << endreq;
86 return StatusCode::SUCCESS;
89StatusCode DsReconstruction::registerEvtRecDTagCol(
91 StatusCode sc = eventSvc()->registerObject(
"/Event/EvtRec/EvtRecDTagCol",
93 if (sc != StatusCode::SUCCESS) {
94 log << MSG::FATAL <<
"Could not register EvtRecDTagCol in TDS!" << endreq;
101 MsgStream log(
msgSvc(), name());
102 log << MSG::INFO <<
"in execute()" << endreq;
109 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
110 int event= eventHeader->eventNumber();
116 log << MSG::DEBUG <<
"run and event = " << eventHeader->runNumber()
117 <<
" " << eventHeader->eventNumber() << endreq;
118 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
119 << recEvent->totalCharged() <<
" , "
120 << recEvent->totalNeutral() <<
" , "
121 << recEvent->totalTracks() <<endreq;
130 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
132 log << MSG::FATAL <<
"Could not find EvtRecPi0Col" << endreq;
133 return StatusCode::FAILURE;
136 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(),
"/Event/EvtRec/EvtRecEtaToGGCol");
137 if ( ! recEtaToGGCol ) {
138 log << MSG::FATAL <<
"Could not find EvtRecEtaToGGCol" << endreq;
139 return StatusCode::FAILURE;
143 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(),
"/Event/EvtRec/EvtRecVeeVertexCol");
144 if ( ! recVeeVertexCol ) {
145 log << MSG::FATAL <<
"Could not find EvtRecVeeVertexCol" << endreq;
146 return StatusCode::FAILURE;
152 log << MSG::FATAL <<
"EvtRecDTagCol is not registered yet" << endreq;
153 return StatusCode::FAILURE;
158 Hep3Vector xorigin(0,0,0);
160 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
165 xorigin.setX(vertex[0]);
166 xorigin.setY(vertex[1]);
167 xorigin.setZ(vertex[2]);
194 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
197 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
200 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
209 int run = eventHeader->runNumber();
210 m_ievt = eventHeader->eventNumber();
211 m_nChrg = recEvent->totalCharged();
212 m_nNeu = recEvent->totalNeutral();
213 m_nPion = pionList.
size();
214 m_nKaon = kaonList.
size();
215 m_nPi0 = pi0List.size();
216 m_nKs = ksList.size();
223 if(m_ReadBeamEFromDB && m_irun!=run){
227 m_beamE=m_readDb.
getbeamE(m_irun,m_beamE);
232 double ebeam=m_beamE;
239 for(
int list=0;list<chanlist.size();list++){
241 string channel=chanlist[list];
251 if(channel==
"DstoKsK") {
253 numchan.push_back(5);
254 numchan.push_back(1);
255 decaylist=ksList* kaonList.
plus();
257 else if(channel==
"DstoKKPi") {
259 numchan.push_back(1);
260 numchan.push_back(1);
261 numchan.push_back(2);
262 decaylist=kaonList.
minus()* kaonList.
plus()* pionList.
plus();
264 else if(channel==
"DstoKsKPi0") {
266 numchan.push_back(5);
267 numchan.push_back(1);
268 numchan.push_back(3);
269 decaylist=ksList* kaonList.
plus()* pi0List;
271 else if(channel==
"DstoKsKsPi") {
273 numchan.push_back(5);
274 numchan.push_back(5);
275 numchan.push_back(2);
276 decaylist=ksList* ksList* pionList.
plus();
278 else if(channel==
"DstoKKPiPi0") {
280 numchan.push_back(1);
281 numchan.push_back(1);
282 numchan.push_back(2);
283 numchan.push_back(3);
284 decaylist=kaonList.
minus()* kaonList.
plus()* pionList.
plus()* pi0List;
286 else if(channel==
"DstoKsKplusPiPi") {
288 numchan.push_back(5);
289 numchan.push_back(1);
290 numchan.push_back(2);
291 numchan.push_back(2);
292 decaylist=ksList* kaonList.
plus()* pionList.
plus()* pionList.
minus();
294 else if(channel==
"DstoKsKminusPiPi") {
296 numchan.push_back(5);
297 numchan.push_back(1);
298 numchan.push_back(2);
299 numchan.push_back(2);
300 decaylist=ksList* kaonList.
minus()* pionList.
plus()* pionList.
plus();
302 else if(channel==
"DstoKKPiPiPi") {
304 numchan.push_back(1);
305 numchan.push_back(1);
306 numchan.push_back(2);
307 numchan.push_back(2);
308 numchan.push_back(2);
311 else if(channel==
"DstoPiPi0") {
313 numchan.push_back(2);
314 numchan.push_back(3);
315 decaylist=pionList.
plus()* pi0List;
317 else if(channel==
"DstoPiPiPi") {
319 numchan.push_back(2);
320 numchan.push_back(2);
321 numchan.push_back(2);
322 decaylist=pionList.
plus()* pionList.
plus()* pionList.
minus();
324 else if(channel==
"DstoPiPiPiPi0") {
326 numchan.push_back(2);
327 numchan.push_back(2);
328 numchan.push_back(2);
329 numchan.push_back(3);
330 decaylist=pionList.
plus()* pionList.
plus()* pionList.
minus()* pi0List;
332 else if(channel==
"DstoPiPiPiPiPi") {
334 numchan.push_back(2);
335 numchan.push_back(2);
336 numchan.push_back(2);
337 numchan.push_back(2);
338 numchan.push_back(2);
341 else if(channel==
"DstoPiPiPiPiPiPi0") {
343 numchan.push_back(2);
344 numchan.push_back(2);
345 numchan.push_back(2);
346 numchan.push_back(2);
347 numchan.push_back(2);
348 numchan.push_back(3);
351 else if(channel==
"DstoPiEta") {
353 numchan.push_back(2);
354 numchan.push_back(4);
355 decaylist=pionList.
plus()* etaList;
357 else if(channel==
"DstoPiPi0Eta") {
359 numchan.push_back(2);
360 numchan.push_back(3);
361 numchan.push_back(4);
362 decaylist=pionList.
plus()* pi0List* etaList;
364 else if(channel==
"DstoPiPiPiEta") {
366 numchan.push_back(2);
367 numchan.push_back(2);
368 numchan.push_back(2);
369 numchan.push_back(4);
370 decaylist=pionList.
plus()* pionList.
plus()* pionList.
minus()* etaList;
372 else if(channel==
"DstoPiEPPiPiEta") {
374 numchan.push_back(2);
375 numchan.push_back(6);
377 epList=pionList.
plus()* pionList.
minus()* etaList;
378 decaylist=pionList.
plus()* epList;
380 else if(channel==
"DstoPiPi0EPPiPiEta") {
382 numchan.push_back(2);
383 numchan.push_back(3);
384 numchan.push_back(6);
386 epList=pionList.
plus()* pionList.
minus()* etaList;
387 decaylist=pionList.
plus()* pi0List* epList;
389 else if(channel==
"DstoPiEPRhoGam") {
391 numchan.push_back(2);
392 numchan.push_back(7);
394 rhoList=pionList.
plus()* pionList.
minus();
396 epList=rhoList* photonList;
397 decaylist=pionList.
plus()* epList;
399 else if(channel==
"DstoPiPi0EPRhoGam") {
401 numchan.push_back(2);
402 numchan.push_back(3);
403 numchan.push_back(7);
405 rhoList=pionList.
plus()* pionList.
minus();
407 epList=rhoList* photonList;
408 decaylist=pionList.
plus()* pi0List* epList;
410 else if(channel==
"DstoKsPi") {
412 numchan.push_back(5);
413 numchan.push_back(2);
414 decaylist=ksList* pionList.
plus();
416 else if(channel==
"DstoKsPiPi0") {
418 numchan.push_back(5);
419 numchan.push_back(2);
420 numchan.push_back(3);
421 decaylist=ksList* pionList.
plus()* pi0List;
423 else if(channel==
"DstoKPiPi") {
425 numchan.push_back(1);
426 numchan.push_back(2);
427 numchan.push_back(2);
428 decaylist=kaonList.
plus()* pionList.
plus()* pionList.
minus();
430 else if(channel==
"DstoKPiPiPi0") {
432 numchan.push_back(1);
433 numchan.push_back(2);
434 numchan.push_back(2);
435 numchan.push_back(3);
436 decaylist=kaonList.
plus()* pionList.
plus()* pionList.
minus()* pi0List;
438 else if(channel==
"DstoKKK") {
440 numchan.push_back(1);
441 numchan.push_back(1);
442 numchan.push_back(1);
443 decaylist=kaonList.
minus()* kaonList.
plus()* kaonList.
plus();
454 vector<int> trackid, showerid;
455 vector<int> kaonid, pionid;
456 int numofchildren=numchan.size()-1;
458 for(
int i=0; i< numofchildren;i++){
460 const CDCandidate& daughter=(*it).particle().child(i);
464 trackid.push_back(track->
trackId());
465 kaonid.push_back(track->
trackId());
467 else if(numchan[i+1]==2){
469 trackid.push_back(track->
trackId());
470 pionid.push_back(track->
trackId());
472 else if ( numchan[i+1]==3){
475 showerid.push_back(hiEnGamma->
trackId());
476 showerid.push_back(loEnGamma->
trackId());
478 else if ( numchan[i+1]==4){
481 showerid.push_back(hiEnGamma->
trackId());
482 showerid.push_back(loEnGamma->
trackId());
484 else if ( numchan[i+1]==5){
488 trackid.push_back(pion1Trk->
trackId());
489 trackid.push_back(pion2Trk->
trackId());
491 else if (numchan[i+1]==6){
500 trackid.push_back(apiontrk->
trackId());
501 trackid.push_back(spiontrk->
trackId());
502 showerid.push_back(hiEnGamma->
trackId());
503 showerid.push_back(loEnGamma->
trackId());
506 else if (numchan[i+1]==7){
517 trackid.push_back(apiontrk->
trackId());
518 trackid.push_back(spiontrk->
trackId());
519 showerid.push_back(gammatrk->
trackId());
527 saveDsInfo(it, ebeam, numofchildren, recDTag);
528 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
529 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
533 HepLorentzVector p4change_vfit=util.
vfit(channel, kaonid, pionid, xorigin, charged_begin);
534 recDTag->
setp4(recDTag->
p4()+p4change_vfit);
546 recDTagCol->push_back(recDTag);
556 return StatusCode::SUCCESS;
562 double mass = (*it).particle().mass();
563 int charge= (*it).particle().charge();
564 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
568 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
569 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
570 double deltaE_CMS = p4.t() - ebeam;
572 if((*it).particle().userTag()==1)
589 vector<EvtRecTrackIterator> trktemp;
590 vector<EvtRecTrackIterator> shrtemp;
595 bool isothertrack=
true;
596 for(
int i=0; i<trackid.size(); i++){
597 if( (*trk)->trackId()==trackid[i] ){
598 trktemp.push_back(trk);
606 for(
int i=0; i<trackid.size();i++){
607 for(
int j=0; j<trktemp.size(); j++){
609 if( (*trk)->trackId()==trackid[i])
617 bool isothershower=
true;
618 for(
int i=0; i<showerid.size(); i++){
619 if( (*shr)->trackId()==showerid[i] ){
620 shrtemp.push_back(shr);
629 for(
int i=0; i<showerid.size();i++){
630 for(
int j=0; j<shrtemp.size(); j++){
632 if( (*shr)->trackId()==showerid[i])
642 bool iskaon=
false,ispion=
false;
648 recDTag->
addKaonId( (*kit).particle().track() );
652 recDTag->
addPionId( (*pit).particle().track() );
702 inFile.open(filename.c_str());
704 cout <<
"Unable to open decay list file";
708 while (inFile >> channel) {
709 temp.push_back(channel);
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
ObjectVector< EvtRecDTag > EvtRecDTagCol
EvtRecTrackCol::iterator EvtRecTrackIterator
LocalEptoPiPiEtaSelector eptoPiPiEtaSelector
LocalEptoRhoGamSelector eptoRhoGamSelector
LocalEtatoGGSelector etatoGGSelector
LocalKaonSelector kaonSelector
LocalKsSelector ksSelector
LocalPhotonSelector photonSelector
LocalPi0Selector pi0Selector
LocalPionSelector pionSelector
LocalRhotoPiPiSelector rhotoPiPiSelector
virtual const DecayEvidence & decay() const
virtual const EvtRecTrack * photon() const
virtual const EvtRecTrack * track() const
virtual const EvtRecVeeVertex * navKshort() const
virtual const EvtRecPi0 * navPi0() const
virtual const EvtRecEtaToGG * navEta() const
const CDCandidate & child(unsigned int aPosition) const
vector< string > getlist(string &filename)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
DsReconstruction(const std::string &name, ISvcLocator *pSvcLocator)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void saveDsInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
void setbeta(Hep3Vector beta)
void setebeam(double ebeam)
void addOtherTrack(const SmartRef< EvtRecTrack > track)
void settype(SelType type)
void setdecayMode(DecayMode decayMode)
HepLorentzVector p4() const
void setp4(HepLorentzVector p4)
void setmass(double mass)
void addOtherShower(const SmartRef< EvtRecTrack > shower)
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
void setdeltaE(double deltaE)
void setcharge(int charge)
void addPionId(const SmartRef< EvtRecTrack > pionId)
void setbeamE(double beamE)
void addShower(const SmartRef< EvtRecTrack > shower)
void setnumOfChildren(int numOfChildren)
void addTrack(const SmartRef< EvtRecTrack > track)
const EvtRecTrack * hiEnGamma() const
const EvtRecTrack * loEnGamma() const
const EvtRecTrack * loEnGamma() const
const EvtRecTrack * hiEnGamma() const
SmartRef< EvtRecTrack > & daughter(int i)
virtual bool isVertexValid()=0
virtual double * PrimaryVertex()=0
void setpidtype(int type)
void setpidtype(int type)
bool setcalib(bool calib)
CLHEP::Hep3Vector getbeta()
ChosenChargeList< Charged, CandidateClass > & plus() const
ChosenChargeList< Charged, CandidateClass > & minus() const
iterator particle_begin()
virtual iterator particle_begin()
virtual iterator particle_end()
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecDTagCol
_EXTERN_ std::string EvtRecTrackCol