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" );
71 declareProperty(
"UseVFRefine", m_useVFrefine =
true);
76 MsgStream log(
msgSvc(), name());
77 log << MSG::INFO <<
"in initialize()" <<endreq;
85 return StatusCode::SUCCESS;
90 MsgStream log(
msgSvc(), name());
91 log << MSG::INFO <<
"in finalize()" << endreq;
95 return StatusCode::SUCCESS;
98StatusCode LambdaCReconstruction::registerEvtRecDTagCol(
100 StatusCode sc = eventSvc()->registerObject(
"/Event/EvtRec/EvtRecDTagCol",
102 if (sc != StatusCode::SUCCESS) {
103 log << MSG::FATAL <<
"Could not register EvtRecDTagCol in TDS!" << endreq;
110 MsgStream log(
msgSvc(), name());
111 log << MSG::INFO <<
"in execute()" << endreq;
118 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
119 int event= eventHeader->eventNumber();
125 log << MSG::DEBUG <<
"run and event = " << eventHeader->runNumber()
126 <<
" " << eventHeader->eventNumber() << endreq;
127 log << MSG::DEBUG <<
"ncharg, nneu, tottks = "
128 << recEvent->totalCharged() <<
" , "
129 << recEvent->totalNeutral() <<
" , "
130 << recEvent->totalTracks() <<endreq;
139 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(),
"/Event/EvtRec/EvtRecPi0Col");
141 log << MSG::FATAL <<
"Could not find EvtRecPi0Col" << endreq;
142 return StatusCode::FAILURE;
146 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(),
"/Event/EvtRec/EvtRecEtaToGGCol");
147 if ( ! recEtaToGGCol ) {
148 log << MSG::FATAL <<
"Could not find EvtRecEtaToGGCol" << endreq;
149 return StatusCode::FAILURE;
153 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(),
"/Event/EvtRec/EvtRecVeeVertexCol");
154 if ( ! recVeeVertexCol ) {
155 log << MSG::FATAL <<
"Could not find EvtRecVeeVertexCol" << endreq;
156 return StatusCode::FAILURE;
162 log << MSG::FATAL <<
"EvtRecDTagCol is not registered yet" << endreq;
163 return StatusCode::FAILURE;
169 Hep3Vector xorigin(0,0,0);
172 Gaudi::svcLocator()->service(
"VertexDbSvc", vtxsvc);
177 xorigin.setX(vertex[0]);
178 xorigin.setY(vertex[1]);
179 xorigin.setZ(vertex[2]);
198 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
201 map<EvtRecVeeVertex*, vector< double > > fitinfo_Ks;
202 for(
CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
215 map<EvtRecVeeVertex*, vector< double > > fitinfo_Lambda;
227 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
230 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
239 int run = eventHeader->runNumber();
240 m_ievt = eventHeader->eventNumber();
241 m_nChrg = recEvent->totalCharged();
242 m_nNeu = recEvent->totalNeutral();
243 m_nPion = pionList.
size();
244 m_nKaon = kaonList.
size();
245 m_nProton = protonList.
size();
246 m_nPi0 = pi0List.size();
247 m_nKs = ksList.size();
248 m_nLambda = lambdaList.
size();
256 if(m_ReadBeamEFromDB && m_irun!=run){
260 m_beamE=m_readDb.
getbeamE(m_irun,m_beamE);
265 double ebeam=m_beamE;
273 for(
int list=0;list<chanlist.size();list++){
275 string channel=chanlist[list];
285 if(channel==
"LambdacPtoKsP") {
287 numchan.push_back(8);
288 numchan.push_back(5);
289 decaylist=protonList.
plus()* ksList;
291 else if(channel==
"LambdacPtoKPiP") {
293 numchan.push_back(8);
294 numchan.push_back(1);
295 numchan.push_back(2);
296 decaylist=protonList.
plus()* kaonList.
minus()* pionList.
plus();
298 else if(channel==
"LambdacPtoKsPi0P") {
300 numchan.push_back(8);
301 numchan.push_back(5);
302 numchan.push_back(3);
303 decaylist=protonList.
plus()* ksList* pi0List;
305 else if(channel==
"LambdacPtoKsPiPiP") {
307 numchan.push_back(8);
308 numchan.push_back(5);
309 numchan.push_back(2);
310 numchan.push_back(2);
311 decaylist=protonList.
plus()* ksList* pionList.
plus()* pionList.
minus();
313 else if(channel==
"LambdacPtoKPiPi0P") {
315 numchan.push_back(8);
316 numchan.push_back(1);
317 numchan.push_back(2);
318 numchan.push_back(3);
319 decaylist=protonList.
plus()* kaonList.
minus()* pionList.
plus()* pi0List;
321 else if(channel==
"LambdacPtoPiPiP") {
323 numchan.push_back(8);
324 numchan.push_back(2);
325 numchan.push_back(2);
326 decaylist=protonList.
plus()* pionList.
plus()* pionList.
minus();
328 else if(channel==
"LambdacPtoLambdaPi") {
330 numchan.push_back(9);
331 numchan.push_back(2);
332 decaylist=lambdaList()* pionList.
plus();
334 else if(channel==
"LambdacPtoLambdaPiPi0") {
336 numchan.push_back(9);
337 numchan.push_back(2);
338 numchan.push_back(3);
339 decaylist=lambdaList()* pionList.
plus()* pi0List;
341 else if(channel==
"LambdacPtoLambdaPiEta") {
343 numchan.push_back(9);
344 numchan.push_back(2);
345 numchan.push_back(4);
346 decaylist=lambdaList()* pionList.
plus()* etaList;;
348 else if(channel==
"LambdacPtoLambdaPiPiPi") {
350 numchan.push_back(9);
351 numchan.push_back(2);
352 numchan.push_back(2);
353 numchan.push_back(2);
354 decaylist=lambdaList()* pionList.
plus()* pionList.
minus() *pionList.
plus();
356 else if(channel==
"LambdacPtoLambdaPiOmega") {
358 numchan.push_back(9);
359 numchan.push_back(2);
360 numchan.push_back(12);
362 omgList=pionList.
plus()* pionList.
minus() *pi0List;
363 decaylist=lambdaList()* pionList.
plus()* omgList;
365 else if(channel==
"LambdacPtoPiSIGMA0LambdaGam") {
367 numchan.push_back(2);
368 numchan.push_back(10);
370 sgmList=lambdaList()* photonList;
371 decaylist=pionList.
plus()* sgmList;
373 else if(channel==
"LambdacPtoPiPi0SIGMA0LambdaGam") {
375 numchan.push_back(2);
376 numchan.push_back(3);
377 numchan.push_back(10);
379 sgmList=lambdaList()* photonList;
380 decaylist=pionList.
plus()* pi0List* sgmList;
382 else if(channel==
"LambdacPtoPi0SIGMAPi0P") {
384 numchan.push_back(3);
385 numchan.push_back(11);
387 sgmList=pi0List* protonList.
plus();
388 decaylist=pi0List* sgmList;
390 else if(channel==
"LambdacPtoPiPiSIGMAPi0P") {
392 numchan.push_back(2);
393 numchan.push_back(2);
394 numchan.push_back(11);
396 sgmList=pi0List* protonList.
plus();
397 decaylist=pionList.
plus()* pionList.
minus()* sgmList;
399 else if(channel==
"LambdacPtoOmegaSIGMAPi0P") {
401 numchan.push_back(12);
402 numchan.push_back(11);
404 omgList=pionList.
plus()* pionList.
minus()* pi0List;
406 sgmList=pi0List* protonList.
plus();
407 decaylist=omgList* sgmList;
409 else if(channel==
"LambdacPtoPiPiPi0SIGMAPi0P") {
411 numchan.push_back(2);
412 numchan.push_back(2);
413 numchan.push_back(3);
414 numchan.push_back(11);
416 sgmList=pi0List* protonList.
plus();
417 decaylist=pionList.
plus()* pionList.
minus()*pi0List* sgmList;
429 vector<int> trackid, showerid;
430 vector<int> kaonid, pionid, protonid;
431 int numofchildren=numchan.size()-1;
433 for(
int i=0; i< numofchildren;i++){
435 const CDCandidate& daughter=(*it).particle().child(i);
439 trackid.push_back(track->
trackId());
440 kaonid.push_back(track->
trackId());
442 else if(numchan[i+1]==2){
444 trackid.push_back(track->
trackId());
445 pionid.push_back(track->
trackId());
447 else if(numchan[i+1]==8){
449 trackid.push_back(track->
trackId());
450 protonid.push_back(track->
trackId());
452 else if ( numchan[i+1]==3){
455 showerid.push_back(hiEnGamma->
trackId());
456 showerid.push_back(loEnGamma->
trackId());
458 else if ( numchan[i+1]==4){
461 showerid.push_back(hiEnGamma->
trackId());
462 showerid.push_back(loEnGamma->
trackId());
464 else if ( numchan[i+1]==5){
466 recDTag->
addToFitInfo(aKsCand->
mass(),fitinfo_Ks[aKsCand][0],fitinfo_Ks[aKsCand][1],fitinfo_Ks[aKsCand][2]);
469 trackid.push_back(pion1Trk->
trackId());
470 trackid.push_back(pion2Trk->
trackId());
472 else if ( numchan[i+1]==9){
474 recDTag->
addToFitInfo(aLambdaCand->
mass(),fitinfo_Lambda[aLambdaCand][0],fitinfo_Lambda[aLambdaCand][1],fitinfo_Lambda[aLambdaCand][2]);
476 int index[2] = {0, 1};
484 trackid.push_back(protonTrk->
trackId());
485 trackid.push_back(pionTrk->
trackId());
487 else if (numchan[i+1]==6){
496 trackid.push_back(apiontrk->
trackId());
497 trackid.push_back(spiontrk->
trackId());
498 showerid.push_back(hiEnGamma->
trackId());
499 showerid.push_back(loEnGamma->
trackId());
502 else if (numchan[i+1]==7){
513 trackid.push_back(apiontrk->
trackId());
514 trackid.push_back(spiontrk->
trackId());
515 showerid.push_back(gammatrk->
trackId());
518 else if (numchan[i+1]==10){
521 showerid.push_back(gammatrk->
trackId());
524 recDTag->
addToFitInfo(aLambdaCand->
mass(),fitinfo_Lambda[aLambdaCand][0],fitinfo_Lambda[aLambdaCand][1],fitinfo_Lambda[aLambdaCand][2]);
526 int index[2] = {0, 1};
534 trackid.push_back(protonTrk->
trackId());
535 trackid.push_back(pionTrk->
trackId());
538 else if (numchan[i+1]==11){
549 showerid.push_back(agamtrk->
trackId());
550 showerid.push_back(sgamtrk->
trackId());
551 trackid.push_back(Ptrk->
trackId());
553 else if (numchan[i+1]==12){
562 trackid.push_back(apiontrk->
trackId());
563 trackid.push_back(spiontrk->
trackId());
564 showerid.push_back(hiEnGamma->
trackId());
565 showerid.push_back(loEnGamma->
trackId());
573 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
574 pidtag(kaonid,pionid,protonid, kaonList_tight, pionList_tight, protonList_tight, recDTag);
577 if(m_debug) cout<<
"beforevfit:"<<endl;
579 HepLorentzVector p4change_vfit;
582 p4change_vfit=util.
vfitref(channel, kaonid, pionid, xorigin, charged_begin);
584 p4change_vfit=util.
vfit(channel, kaonid, pionid, xorigin, charged_begin);
587 recDTag->
setp4(recDTag->
p4()+p4change_vfit);
598 recDTagCol->push_back(recDTag);
606 return StatusCode::SUCCESS;
612 double mass = (*it).particle().mass();
613 int charge= (*it).particle().charge();
614 HepLorentzVector
p4((*it).particle().momentum(), (*it).particle().energy());
619 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
620 double deltaE_CMS =
p4.t() -
ebeam;
622 if((*it).particle().userTag()==1)
639 vector<EvtRecTrackIterator> trktemp;
640 vector<EvtRecTrackIterator> shrtemp;
645 bool isothertrack=
true;
646 for(
int i=0; i<trackid.size(); i++){
647 if( (*trk)->trackId()==trackid[i] ){
648 trktemp.push_back(trk);
656 for(
int i=0; i<trackid.size();i++){
657 for(
int j=0; j<trktemp.size(); j++){
659 if( (*trk)->trackId()==trackid[i])
667 bool isothershower=
true;
668 for(
int i=0; i<showerid.size(); i++){
669 if( (*shr)->trackId()==showerid[i] ){
670 shrtemp.push_back(shr);
679 for(
int i=0; i<showerid.size();i++){
680 for(
int j=0; j<shrtemp.size(); j++){
682 if( (*shr)->trackId()==showerid[i])
692 bool iskaon=
false,ispion=
false,isproton=
false;
698 recDTag->
addKaonId( (*kit).particle().track() );
702 recDTag->
addPionId( (*pit).particle().track() );
756 inFile.open(filename.c_str());
758 cout <<
"Unable to open decay list file";
762 while (inFile >> channel) {
763 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)
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()
HepLorentzVector vfitref(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
vector< double > SecondaryVFit_Lambdaref(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc)
vector< double > SecondaryVFit_Lambda(EvtRecVeeVertex *lambda, IVertexDbSvc *vtxsvc)
vector< double > SecondaryVFitref(EvtRecVeeVertex *ks, 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
double double double * p4