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"
51Algorithm(name, pSvcLocator) {
53 declareProperty(
"debug", m_debug =
false );
54 declareProperty(
"ReadBeamEFromDB", m_ReadBeamEFromDB =
false );
55 declareProperty(
"UseCalibBeamE", m_usecalibBeamE =
false );
56 declareProperty(
"BeamE", m_beamE = 1.8865 );
57 declareProperty(
"DpList", m_decaylist =
"test.txt" );
62 MsgStream log(
msgSvc(), name());
63 log << MSG::INFO <<
"in initialize()" <<endreq;
71 return StatusCode::SUCCESS;
76 MsgStream log(
msgSvc(), name());
77 log << MSG::INFO <<
"in finalize()" << endreq;
81 return StatusCode::SUCCESS;
84StatusCode ChargedDReconstruction::registerEvtRecDTagCol(
86 StatusCode sc = eventSvc()->registerObject(
"/Event/EvtRec/EvtRecDTagCol",
88 if (sc != StatusCode::SUCCESS) {
89 log << MSG::FATAL <<
"Could not register EvtRecDTagCol in TDS!" << endreq;
96 MsgStream log(
msgSvc(), name());
97 log << MSG::INFO <<
"in execute()" << endreq;
104 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
105 int event= eventHeader->eventNumber();
111 log << MSG::DEBUG <<
"run and event = " << eventHeader->runNumber()
112 <<
" " << eventHeader->eventNumber() << endreq;
113 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
114 << recEvent->totalCharged() <<
" , "
115 << recEvent->totalNeutral() <<
" , "
116 << recEvent->totalTracks() <<endreq;
125 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
127 log << MSG::FATAL <<
"Could not find EvtRecPi0Col" << endreq;
128 return StatusCode::FAILURE;
132 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(),
"/Event/EvtRec/EvtRecEtaToGGCol");
133 if ( ! recEtaToGGCol ) {
134 log << MSG::FATAL <<
"Could not find EvtRecEtaToGGCol" << endreq;
135 return StatusCode::FAILURE;
139 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(),
"/Event/EvtRec/EvtRecVeeVertexCol");
140 if ( ! recVeeVertexCol ) {
141 log << MSG::FATAL <<
"Could not find EvtRecVeeVertexCol" << endreq;
142 return StatusCode::FAILURE;
148 log << MSG::FATAL <<
"EvtRecDTagCol is not registered yet" << endreq;
149 return StatusCode::FAILURE;
172 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
175 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
178 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
188 int run = eventHeader->runNumber();
189 m_ievt = eventHeader->eventNumber();
190 m_nChrg = recEvent->totalCharged();
191 m_nNeu = recEvent->totalNeutral();
192 m_nPion = pionList.
size();
193 m_nKaon = kaonList.
size();
194 m_nPi0 = pi0List.size();
195 m_nKs = ksList.size();
203 if(m_ReadBeamEFromDB && m_irun!=run){
207 m_beamE=m_readDb.
getbeamE(m_irun,m_beamE);
212 double ebeam=m_beamE;
220 for(
int list=0;list<chanlist.size();list++){
222 string channel=chanlist[list];
232 if(channel==
"DptoKPiPi") {
234 numchan.push_back(1);
235 numchan.push_back(2);
236 numchan.push_back(2);
237 decaylist=kaonList.
minus()* pionList.
plus()* pionList.
plus();
239 else if(channel==
"DptoKPiPiPi0") {
241 numchan.push_back(1);
242 numchan.push_back(2);
243 numchan.push_back(2);
244 numchan.push_back(3);
245 decaylist=kaonList.
minus()* pionList.
plus()* pionList.
plus()* pi0List;
247 else if(channel==
"DptoKsPi") {
249 numchan.push_back(5);
250 numchan.push_back(2);
251 decaylist=ksList* pionList.
plus();
253 else if(channel==
"DptoKsPiPi0") {
255 numchan.push_back(5);
256 numchan.push_back(2);
257 numchan.push_back(3);
258 decaylist=ksList* pionList.
plus()* pi0List;
260 else if(channel==
"DptoKsPiPiPi") {
262 numchan.push_back(5);
263 numchan.push_back(2);
264 numchan.push_back(2);
265 numchan.push_back(2);
266 decaylist=ksList* pionList.
plus()* pionList.
plus()* pionList.
minus();
268 else if(channel==
"DptoKKPi") {
270 numchan.push_back(1);
271 numchan.push_back(1);
272 numchan.push_back(2);
273 decaylist=kaonList.
minus()* kaonList.
plus()* pionList.
plus();
275 else if(channel==
"DptoPiPi0") {
277 numchan.push_back(2);
278 numchan.push_back(3);
279 decaylist=pionList.
plus()* pi0List;
281 else if(channel==
"DptoKPi0") {
283 numchan.push_back(1);
284 numchan.push_back(3);
285 decaylist=kaonList.
plus()* pi0List;
287 else if(channel==
"DptoKsK") {
289 numchan.push_back(5);
290 numchan.push_back(1);
291 decaylist=ksList* kaonList.
plus();
293 else if(channel==
"DptoPiPiPi") {
295 numchan.push_back(2);
296 numchan.push_back(2);
297 numchan.push_back(2);
298 decaylist=pionList.
plus()* pionList.
plus()* pionList.
minus() ;
300 else if(channel==
"DptoPiPi0Pi0") {
302 numchan.push_back(2);
303 numchan.push_back(3);
304 numchan.push_back(3);
305 decaylist=pionList.
plus()* pi0List* pi0List;
307 else if(channel==
"DptoKsKsPi") {
309 numchan.push_back(5);
310 numchan.push_back(5);
311 numchan.push_back(2);
312 decaylist=ksList* ksList* pionList.
plus();
314 else if(channel==
"DptoKsKPi0") {
316 numchan.push_back(5);
317 numchan.push_back(1);
318 numchan.push_back(3);
319 decaylist=ksList* kaonList.
plus()* pi0List;
321 else if(channel==
"DptoKsKsK") {
323 numchan.push_back(5);
324 numchan.push_back(5);
325 numchan.push_back(1);
326 decaylist=ksList* ksList* kaonList.
plus();
328 else if(channel==
"DptoPiPiPiPi0") {
330 numchan.push_back(2);
331 numchan.push_back(2);
332 numchan.push_back(2);
333 numchan.push_back(3);
334 decaylist=pionList.
plus()* pionList.
plus()* pionList.
minus()* pi0List;
336 else if(channel==
"DptoKsPiPi0Pi0") {
338 numchan.push_back(5);
339 numchan.push_back(2);
340 numchan.push_back(3);
341 numchan.push_back(3);
342 decaylist=ksList* pionList.
plus()* pi0List* pi0List;
344 else if(channel==
"DptoKsKplusPiPi") {
346 numchan.push_back(5);
347 numchan.push_back(1);
348 numchan.push_back(2);
349 numchan.push_back(2);
350 decaylist=ksList* kaonList.
plus()* pionList.
plus()* pionList.
minus();
352 else if(channel==
"DptoKsKminusPiPi") {
354 numchan.push_back(5);
355 numchan.push_back(1);
356 numchan.push_back(2);
357 numchan.push_back(2);
358 decaylist=ksList* kaonList.
minus()* pionList.
plus()* pionList.
plus();
360 else if(channel==
"DptoKKPiPi0") {
362 numchan.push_back(1);
363 numchan.push_back(1);
364 numchan.push_back(2);
365 numchan.push_back(3);
366 decaylist=kaonList.
minus()* kaonList.
plus()* pionList.
plus()* pi0List;
368 else if(channel==
"DptoPiPiPiPiPi") {
370 numchan.push_back(2);
371 numchan.push_back(2);
372 numchan.push_back(2);
373 numchan.push_back(2);
374 numchan.push_back(2);
377 else if(channel==
"DptoKPiPiPiPi") {
379 numchan.push_back(1);
380 numchan.push_back(2);
381 numchan.push_back(2);
382 numchan.push_back(2);
383 numchan.push_back(2);
386 else if(channel==
"DptoPiEta") {
388 numchan.push_back(2);
389 numchan.push_back(4);
390 decaylist=pionList.
plus()* etaList;
392 else if(channel==
"DptoKsPiEta") {
394 numchan.push_back(5);
395 numchan.push_back(2);
396 numchan.push_back(4);
397 decaylist=ksList* pionList.
plus()* etaList;
408 vector<int> trackid, showerid;
409 vector<int> kaonid, pionid;
410 int numofchildren=numchan.size()-1;
412 for(
int i=0; i< numofchildren;i++){
414 const CDCandidate& daughter=(*it).particle().child(i);
418 trackid.push_back(track->
trackId());
419 kaonid.push_back(track->
trackId());
421 else if(numchan[i+1]==2){
423 trackid.push_back(track->
trackId());
424 pionid.push_back(track->
trackId());
426 else if ( numchan[i+1]==3){
429 showerid.push_back(hiEnGamma->
trackId());
430 showerid.push_back(loEnGamma->
trackId());
432 else if ( numchan[i+1]==4){
435 showerid.push_back(hiEnGamma->
trackId());
436 showerid.push_back(loEnGamma->
trackId());
438 else if ( numchan[i+1]==5){
442 trackid.push_back(pion1Trk->
trackId());
443 trackid.push_back(pion2Trk->
trackId());
449 saveDpInfo(it, ebeam, numofchildren, recDTag);
450 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
451 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
460 recDTagCol->push_back(recDTag);
470 return StatusCode::SUCCESS;
476 double mass = (*it).particle().mass();
477 int charge= (*it).particle().charge();
478 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
482 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
483 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
484 double deltaE_CMS = p4.t() - ebeam;
486 if((*it).particle().userTag()==1)
503 vector<EvtRecTrackIterator> trktemp;
504 vector<EvtRecTrackIterator> shrtemp;
509 bool isothertrack=
true;
510 for(
int i=0; i<trackid.size(); i++){
511 if( (*trk)->trackId()==trackid[i] ){
512 trktemp.push_back(trk);
520 for(
int i=0; i<trackid.size();i++){
521 for(
int j=0; j<trktemp.size(); j++){
523 if( (*trk)->trackId()==trackid[i])
531 bool isothershower=
true;
532 for(
int i=0; i<showerid.size(); i++){
533 if( (*shr)->trackId()==showerid[i] ){
534 shrtemp.push_back(shr);
543 for(
int i=0; i<showerid.size();i++){
544 for(
int j=0; j<shrtemp.size(); j++){
546 if( (*shr)->trackId()==showerid[i])
556 bool iskaon=
false,ispion=
false;
562 recDTag->
addKaonId( (*kit).particle().track() );
566 recDTag->
addPionId( (*pit).particle().track() );
617 inFile.open(filename.c_str());
619 cout <<
"Unable to open decay list file";
623 while (inFile >> channel) {
624 temp.push_back(channel);
ChargedDSelector chargedDSelector
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
ObjectVector< EvtRecDTag > EvtRecDTagCol
EvtRecTrackCol::iterator EvtRecTrackIterator
LocalEtatoGGSelector etatoGGSelector
LocalKaonSelector kaonSelector
LocalKsSelector ksSelector
LocalPi0Selector pi0Selector
LocalPionSelector pionSelector
virtual const EvtRecTrack * track() const
virtual const EvtRecVeeVertex * navKshort() const
virtual const EvtRecPi0 * navPi0() const
virtual const EvtRecEtaToGG * navEta() const
void saveDpInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
vector< string > getlist(string &filename)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
ChargedDReconstruction(const std::string &name, ISvcLocator *pSvcLocator)
void setbeta(Hep3Vector beta)
void setebeam(double ebeam)
void addOtherTrack(const SmartRef< EvtRecTrack > track)
void settype(SelType type)
void setdecayMode(DecayMode decayMode)
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)
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()
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecDTagCol
_EXTERN_ std::string EvtRecTrackCol