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"
63 Algorithm(name, pSvcLocator) {
65 declareProperty(
"debug", m_debug =
false );
66 declareProperty(
"ReadBeamEFromDB", m_ReadBeamEFromDB =
false );
67 declareProperty(
"UseCalibBeamE", m_usecalibBeamE =
false );
68 declareProperty(
"UseVertexfit", m_usevertexfit =
false );
69 declareProperty(
"BeamE", m_beamE = 2.3);
70 declareProperty(
"LcList", m_decaylist =
"test.txt" );
75 MsgStream log(
msgSvc(), name());
76 log << MSG::INFO <<
"in initialize()" <<endreq;
84 return StatusCode::SUCCESS;
89 MsgStream log(
msgSvc(), name());
90 log << MSG::INFO <<
"in finalize()" << endreq;
94 return StatusCode::SUCCESS;
97StatusCode LambdaCReconstruction::registerEvtRecDTagCol(
99 StatusCode sc = eventSvc()->registerObject(
"/Event/EvtRec/EvtRecDTagCol",
101 if (sc != StatusCode::SUCCESS) {
102 log << MSG::FATAL <<
"Could not register EvtRecDTagCol in TDS!" << endreq;
109 MsgStream log(
msgSvc(), name());
110 log << MSG::INFO <<
"in execute()" << endreq;
117 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
118 int event= eventHeader->eventNumber();
124 log << MSG::DEBUG <<
"run and event = " << eventHeader->runNumber()
125 <<
" " << eventHeader->eventNumber() << endreq;
126 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
127 << recEvent->totalCharged() <<
" , "
128 << recEvent->totalNeutral() <<
" , "
129 << recEvent->totalTracks() <<endreq;
138 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
140 log << MSG::FATAL <<
"Could not find EvtRecPi0Col" << endreq;
141 return StatusCode::FAILURE;
145 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(),
"/Event/EvtRec/EvtRecEtaToGGCol");
146 if ( ! recEtaToGGCol ) {
147 log << MSG::FATAL <<
"Could not find EvtRecEtaToGGCol" << endreq;
148 return StatusCode::FAILURE;
152 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(),
"/Event/EvtRec/EvtRecVeeVertexCol");
153 if ( ! recVeeVertexCol ) {
154 log << MSG::FATAL <<
"Could not find EvtRecVeeVertexCol" << endreq;
155 return StatusCode::FAILURE;
161 log << MSG::FATAL <<
"EvtRecDTagCol is not registered yet" << endreq;
162 return StatusCode::FAILURE;
168 Hep3Vector xorigin(0,0,0);
171 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
176 xorigin.setX(vertex[0]);
177 xorigin.setY(vertex[1]);
178 xorigin.setZ(vertex[2]);
197 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
200 map<EvtRecVeeVertex*, vector< double > > fitinfo_Ks;
201 for(
CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
209 map<EvtRecVeeVertex*, vector< double > > fitinfo_Lambda;
216 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
219 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
228 int run = eventHeader->runNumber();
229 m_ievt = eventHeader->eventNumber();
230 m_nChrg = recEvent->totalCharged();
231 m_nNeu = recEvent->totalNeutral();
232 m_nPion = pionList.
size();
233 m_nKaon = kaonList.
size();
234 m_nProton = protonList.
size();
235 m_nPi0 = pi0List.size();
236 m_nKs = ksList.size();
237 m_nLambda = lambdaList.
size();
245 if(m_ReadBeamEFromDB && m_irun!=run){
249 m_beamE=m_readDb.
getbeamE(m_irun,m_beamE);
254 double ebeam=m_beamE;
262 for(
int list=0;list<chanlist.size();list++){
264 string channel=chanlist[list];
274 if(channel==
"LambdacPtoKsP") {
276 numchan.push_back(8);
277 numchan.push_back(5);
278 decaylist=protonList.
plus()* ksList;
280 else if(channel==
"LambdacPtoKPiP") {
282 numchan.push_back(8);
283 numchan.push_back(1);
284 numchan.push_back(2);
285 decaylist=protonList.
plus()* kaonList.
minus()* pionList.
plus();
287 else if(channel==
"LambdacPtoKsPi0P") {
289 numchan.push_back(8);
290 numchan.push_back(5);
291 numchan.push_back(3);
292 decaylist=protonList.
plus()* ksList* pi0List;
294 else if(channel==
"LambdacPtoKsPiPiP") {
296 numchan.push_back(8);
297 numchan.push_back(5);
298 numchan.push_back(2);
299 numchan.push_back(2);
300 decaylist=protonList.
plus()* ksList* pionList.
plus()* pionList.
minus();
302 else if(channel==
"LambdacPtoKPiPi0P") {
304 numchan.push_back(8);
305 numchan.push_back(1);
306 numchan.push_back(2);
307 numchan.push_back(3);
308 decaylist=protonList.
plus()* kaonList.
minus()* pionList.
plus()* pi0List;
310 else if(channel==
"LambdacPtoPiPiP") {
312 numchan.push_back(8);
313 numchan.push_back(2);
314 numchan.push_back(2);
315 decaylist=protonList.
plus()* pionList.
plus()* pionList.
minus();
317 else if(channel==
"LambdacPtoLambdaPi") {
319 numchan.push_back(9);
320 numchan.push_back(2);
321 decaylist=lambdaList()* pionList.
plus();
323 else if(channel==
"LambdacPtoLambdaPiPi0") {
325 numchan.push_back(9);
326 numchan.push_back(2);
327 numchan.push_back(3);
328 decaylist=lambdaList()* pionList.
plus()* pi0List;
330 else if(channel==
"LambdacPtoLambdaPiEta") {
332 numchan.push_back(9);
333 numchan.push_back(2);
334 numchan.push_back(4);
335 decaylist=lambdaList()* pionList.
plus()* etaList;;
337 else if(channel==
"LambdacPtoLambdaPiPiPi") {
339 numchan.push_back(9);
340 numchan.push_back(2);
341 numchan.push_back(2);
342 numchan.push_back(2);
343 decaylist=lambdaList()* pionList.
plus()* pionList.
minus() *pionList.
plus();
345 else if(channel==
"LambdacPtoLambdaPiOmega") {
347 numchan.push_back(9);
348 numchan.push_back(2);
349 numchan.push_back(12);
351 omgList=pionList.
plus()* pionList.
minus() *pi0List;
352 decaylist=lambdaList()* pionList.
plus()* omgList;
354 else if(channel==
"LambdacPtoPiSIGMA0LambdaGam") {
356 numchan.push_back(2);
357 numchan.push_back(10);
359 sgmList=lambdaList()* photonList;
360 decaylist=pionList.
plus()* sgmList;
362 else if(channel==
"LambdacPtoPiPi0SIGMA0LambdaGam") {
364 numchan.push_back(2);
365 numchan.push_back(3);
366 numchan.push_back(10);
368 sgmList=lambdaList()* photonList;
369 decaylist=pionList.
plus()* pi0List* sgmList;
371 else if(channel==
"LambdacPtoPi0SIGMAPi0P") {
373 numchan.push_back(3);
374 numchan.push_back(11);
376 sgmList=pi0List* protonList.
plus();
377 decaylist=pi0List* sgmList;
379 else if(channel==
"LambdacPtoPiPiSIGMAPi0P") {
381 numchan.push_back(2);
382 numchan.push_back(2);
383 numchan.push_back(11);
385 sgmList=pi0List* protonList.
plus();
386 decaylist=pionList.
plus()* pionList.
minus()* sgmList;
388 else if(channel==
"LambdacPtoOmegaSIGMAPi0P") {
390 numchan.push_back(12);
391 numchan.push_back(11);
393 omgList=pionList.
plus()* pionList.
minus()* pi0List;
395 sgmList=pi0List* protonList.
plus();
396 decaylist=omgList* sgmList;
398 else if(channel==
"LambdacPtoPiPiPi0SIGMAPi0P") {
400 numchan.push_back(2);
401 numchan.push_back(2);
402 numchan.push_back(3);
403 numchan.push_back(11);
405 sgmList=pi0List* protonList.
plus();
406 decaylist=pionList.
plus()* pionList.
minus()*pi0List* sgmList;
418 vector<int> trackid, showerid;
419 vector<int> kaonid, pionid, protonid;
420 int numofchildren=numchan.size()-1;
422 for(
int i=0; i< numofchildren;i++){
424 const CDCandidate& daughter=(*it).particle().child(i);
428 trackid.push_back(track->
trackId());
429 kaonid.push_back(track->
trackId());
431 else if(numchan[i+1]==2){
433 trackid.push_back(track->
trackId());
434 pionid.push_back(track->
trackId());
436 else if(numchan[i+1]==8){
438 trackid.push_back(track->
trackId());
439 protonid.push_back(track->
trackId());
441 else if ( numchan[i+1]==3){
444 showerid.push_back(hiEnGamma->
trackId());
445 showerid.push_back(loEnGamma->
trackId());
447 else if ( numchan[i+1]==4){
450 showerid.push_back(hiEnGamma->
trackId());
451 showerid.push_back(loEnGamma->
trackId());
453 else if ( numchan[i+1]==5){
455 recDTag->
addToFitInfo(aKsCand->
mass(),fitinfo_Ks[aKsCand][0],fitinfo_Ks[aKsCand][1],fitinfo_Ks[aKsCand][2]);
458 trackid.push_back(pion1Trk->
trackId());
459 trackid.push_back(pion2Trk->
trackId());
461 else if ( numchan[i+1]==9){
463 recDTag->
addToFitInfo(aLambdaCand->
mass(),fitinfo_Lambda[aLambdaCand][0],fitinfo_Lambda[aLambdaCand][1],fitinfo_Lambda[aLambdaCand][2]);
465 int index[2] = {0, 1};
473 trackid.push_back(protonTrk->
trackId());
474 trackid.push_back(pionTrk->
trackId());
476 else if (numchan[i+1]==6){
485 trackid.push_back(apiontrk->
trackId());
486 trackid.push_back(spiontrk->
trackId());
487 showerid.push_back(hiEnGamma->
trackId());
488 showerid.push_back(loEnGamma->
trackId());
491 else if (numchan[i+1]==7){
502 trackid.push_back(apiontrk->
trackId());
503 trackid.push_back(spiontrk->
trackId());
504 showerid.push_back(gammatrk->
trackId());
507 else if (numchan[i+1]==10){
510 showerid.push_back(gammatrk->
trackId());
513 recDTag->
addToFitInfo(aLambdaCand->
mass(),fitinfo_Lambda[aLambdaCand][0],fitinfo_Lambda[aLambdaCand][1],fitinfo_Lambda[aLambdaCand][2]);
515 int index[2] = {0, 1};
523 trackid.push_back(protonTrk->
trackId());
524 trackid.push_back(pionTrk->
trackId());
527 else if (numchan[i+1]==11){
538 showerid.push_back(agamtrk->
trackId());
539 showerid.push_back(sgamtrk->
trackId());
540 trackid.push_back(Ptrk->
trackId());
542 else if (numchan[i+1]==12){
551 trackid.push_back(apiontrk->
trackId());
552 trackid.push_back(spiontrk->
trackId());
553 showerid.push_back(hiEnGamma->
trackId());
554 showerid.push_back(loEnGamma->
trackId());
562 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
563 pidtag(kaonid,pionid,protonid, kaonList_tight, pionList_tight, protonList_tight, recDTag);
566 HepLorentzVector p4change_vfit=util.
vfit(channel, kaonid, pionid, protonid, xorigin, charged_begin);
567 recDTag->
setp4(recDTag->
p4()+p4change_vfit);
578 recDTagCol->push_back(recDTag);
586 return StatusCode::SUCCESS;
592 double mass = (*it).particle().mass();
593 int charge= (*it).particle().charge();
594 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
598 double mbc2_CMS =
ebeam*
ebeam - p4.v().mag2();
599 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
600 double deltaE_CMS = p4.t() -
ebeam;
602 if((*it).particle().userTag()==1)
619 vector<EvtRecTrackIterator> trktemp;
620 vector<EvtRecTrackIterator> shrtemp;
625 bool isothertrack=
true;
626 for(
int i=0; i<trackid.size(); i++){
627 if( (*trk)->trackId()==trackid[i] ){
628 trktemp.push_back(trk);
636 for(
int i=0; i<trackid.size();i++){
637 for(
int j=0; j<trktemp.size(); j++){
639 if( (*trk)->trackId()==trackid[i])
647 bool isothershower=
true;
648 for(
int i=0; i<showerid.size(); i++){
649 if( (*shr)->trackId()==showerid[i] ){
650 shrtemp.push_back(shr);
659 for(
int i=0; i<showerid.size();i++){
660 for(
int j=0; j<shrtemp.size(); j++){
662 if( (*shr)->trackId()==showerid[i])
672 bool iskaon=
false,ispion=
false,isproton=
false;
678 recDTag->
addKaonId( (*kit).particle().track() );
682 recDTag->
addPionId( (*pit).particle().track() );
736 inFile.open(filename.c_str());
738 cout <<
"Unable to open decay list file";
742 while (inFile >> channel) {
743 temp.push_back(channel);
double P(RecMdcKalTrack *trk)
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
ObjectVector< EvtRecDTag > EvtRecDTagCol
EvtRecTrackCol::iterator EvtRecTrackIterator
LambdaCSelector lambdaCSelector
LocalChargedSigmaSelector chargedsigmaSelector
LocalEtatoGGSelector etatoGGSelector
LocalKaonSelector kaonSelector
LocalKsSelector ksSelector
LocalLambdaSelector lambdaSelector
LocalPhotonSelector photonSelector
LocalPi0Selector pi0Selector
LocalPionSelector pionSelector
LocalProtonSelector protonSelector
LocalSigma0Selector sigma0Selector
LocalomegatoPiPiPi0Selector omegatoPiPiPi0Selector
virtual const DecayEvidence & decay() const
virtual const EvtRecTrack * photon() const
virtual const EvtRecVeeVertex * navLambda() 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
@ kLambdacPtoLambdaPiOmega
@ kLambdacPtoPiSIGMA0LambdaGam
@ kLambdacPtoOmegaSIGMAPi0P
@ kLambdacPtoPi0SIGMAPi0P
@ kLambdacPtoPiPiSIGMAPi0P
@ kLambdacPtoPiPiPi0SIGMAPi0P
@ kLambdacPtoPiPi0SIGMA0LambdaGam
@ kLambdacPtoLambdaPiPiPi
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 addToFitInfo(double ksmass, double chi2, double length, double error)
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 addProtonId(const SmartRef< EvtRecTrack > protonId)
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 savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
void pidtag(vector< int >, vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, CDProtonList &, EvtRecDTag *)
void saveLcInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
vector< string > getlist(string &filename)
LambdaCReconstruction(const std::string &name, ISvcLocator *pSvcLocator)
void setebeam(double ebeam)
void setbeta(Hep3Vector beta)
void setpidtype(int type)
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()
iterator particle_begin()
vector< double > SecondaryVFit_Lambda(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc)
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
vector< double > SecondaryVFit(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecDTagCol
_EXTERN_ std::string EvtRecTrackCol