00001
00002
00003
00004 #include <fstream>
00005 #include "GaudiKernel/MsgStream.h"
00006 #include "GaudiKernel/AlgFactory.h"
00007 #include "GaudiKernel/ISvcLocator.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009 #include "GaudiKernel/IDataProviderSvc.h"
00010 #include "GaudiKernel/PropertyMgr.h"
00011
00012 #include "EventModel/EventHeader.h"
00013 #include "EventModel/Event.h"
00014 #include "EventModel/EventModel.h"
00015
00016 #include "EvtRecEvent/EvtRecEvent.h"
00017 #include "EvtRecEvent/EvtRecTrack.h"
00018 #include "EvtRecEvent/EvtRecDTag.h"
00019 #include "EvtRecEvent/EvtRecPi0.h"
00020 #include "EvtRecEvent/EvtRecEtaToGG.h"
00021 #include "EvtRecEvent/EvtRecVeeVertex.h"
00022
00023 #include "BesDChain/CDChargedPionList.h"
00024 #include "BesDChain/CDChargedKaonList.h"
00025 #include "BesDChain/CDPhotonList.h"
00026 #include "BesDChain/CDPi0List.h"
00027 #include "BesDChain/CDEtaList.h"
00028 #include "BesDChain/CDKsList.h"
00029 #include "BesDChain/CDDecayList.h"
00030
00031 #include "McTruth/McParticle.h"
00032 #include "ParticleID/ParticleID.h"
00033 #include "VertexFit/VertexFit.h"
00034 #include "VertexFit/SecondVertexFit.h"
00035
00036
00037 #include "DTagAlg/DsReconstruction.h"
00038 #include "DTagAlg/LocalKaonSelector.h"
00039 #include "DTagAlg/LocalPionSelector.h"
00040 #include "DTagAlg/LocalPhotonSelector.h"
00041 #include "DTagAlg/LocalKsSelector.h"
00042 #include "DTagAlg/LocalPi0Selector.h"
00043 #include "DTagAlg/LocalEtatoGGSelector.h"
00044 #include "DTagAlg/LocalEtatoPiPiPi0Selector.h"
00045 #include "DTagAlg/LocalEptoPiPiEtaSelector.h"
00046 #include "DTagAlg/LocalEptoPiPiEta3PiSelector.h"
00047 #include "DTagAlg/LocalEptoRhoGamSelector.h"
00048 #include "DTagAlg/LocalRhotoPiPiSelector.h"
00049 #include "DTagAlg/DsSelector.h"
00050
00051 #include "DTagAlg/utility.h"
00052
00053 using namespace Event;
00054
00055
00056 DsReconstruction::DsReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
00057 Algorithm(name, pSvcLocator) {
00058
00059 declareProperty( "debug", m_debug = false );
00060 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
00061 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
00062 declareProperty( "UseVertexfit", m_usevertexfit = false );
00063 declareProperty( "BeamE", m_beamE=2.015 );
00064 declareProperty( "DsList", m_decaylist = "test.txt" );
00065 declareProperty("ReadMeasuredEcms", m_RdMeasuredEcms = false);
00066 }
00067
00068
00069 StatusCode DsReconstruction::initialize() {
00070 MsgStream log(msgSvc(), name());
00071 log << MSG::INFO << "in initialize()" <<endreq;
00072
00073 m_irun=-100;
00074 m_beta.setX(0.011);
00075 m_beta.setY(0);
00076 m_beta.setZ(0);
00077 chanlist=getlist(m_decaylist);
00078
00079
00080 if(m_RdMeasuredEcms){
00081 StatusCode status=serviceLocator()->service("MeasuredEcmsSvc", ecmsSvc, true);
00082 if(!status.isSuccess()){
00083 std::cout<<"ERROR: Can not initial the IMeasuredEcmsSvc right"<<std::endl;
00084 return status;
00085 }
00086 }
00087
00088
00089 return StatusCode::SUCCESS;
00090 }
00091
00092
00093 StatusCode DsReconstruction::finalize() {
00094 MsgStream log(msgSvc(), name());
00095 log << MSG::INFO << "in finalize()" << endreq;
00096
00097 chanlist.clear();
00098
00099 return StatusCode::SUCCESS;
00100 }
00101
00102 StatusCode DsReconstruction::registerEvtRecDTagCol(
00103 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
00104 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
00105 aNewEvtRecDTagCol);
00106 if (sc != StatusCode::SUCCESS) {
00107 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
00108 }
00109 return sc;
00110 }
00111
00112
00113 StatusCode DsReconstruction::execute() {
00114 MsgStream log(msgSvc(), name());
00115 log << MSG::INFO << "in execute()" << endreq;
00116
00117 StatusCode sc;
00118
00120
00122 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00123 int event= eventHeader->eventNumber();
00124
00125
00126
00127 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00128 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00129 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00130 << " " << eventHeader->eventNumber() << endreq;
00131 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00132 << recEvent->totalCharged() << " , "
00133 << recEvent->totalNeutral() << " , "
00134 << recEvent->totalTracks() <<endreq;
00135
00136 EvtRecTrackIterator charged_begin = recTrackCol->begin();
00137 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
00138
00139 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
00140 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
00141
00142
00143 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
00144 if ( ! recPi0Col ) {
00145 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
00146 return StatusCode::FAILURE;
00147 }
00148
00149 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
00150 if ( ! recEtaToGGCol ) {
00151 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
00152 return StatusCode::FAILURE;
00153 }
00154
00155
00156 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
00157 if ( ! recVeeVertexCol ) {
00158 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
00159 return StatusCode::FAILURE;
00160 }
00161
00162
00163 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
00164 if (!recDTagCol) {
00165 log << MSG::FATAL << "EvtRecDTagCol is not registered yet" << endreq;
00166 return StatusCode::FAILURE;
00167 }
00168
00169
00170
00171 Hep3Vector xorigin(0,0,0);
00172 IVertexDbSvc* vtxsvc;
00173 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00174 if (vtxsvc->isVertexValid()) {
00175
00176
00177 double* vertex = vtxsvc->PrimaryVertex();
00178 xorigin.setX(vertex[0]);
00179 xorigin.setY(vertex[1]);
00180 xorigin.setZ(vertex[2]);
00181 }
00182 utility util;
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00198
00200
00201 pionSelector.setpidtype(0);
00202 kaonSelector.setpidtype(0);
00203 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
00204 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
00205 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
00206
00207 CDKsList ksList(ksSelector);
00208 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
00209
00210
00211
00212 map<EvtRecVeeVertex*, vector< double > > fitinfo;
00213 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
00214 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
00215 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
00216 }
00217
00218 CDPi0List pi0List(pi0Selector);
00219 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
00220
00221 CDEtaList etaList(etatoGGSelector);
00222 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
00223
00224
00225 pionSelector.setpidtype(1);
00226 kaonSelector.setpidtype(1);
00227 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
00228 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
00229
00230
00231 int run = eventHeader->runNumber();
00232 m_ievt = eventHeader->eventNumber();
00233 m_nChrg = recEvent->totalCharged();
00234 m_nNeu = recEvent->totalNeutral();
00235 m_nPion = pionList.size();
00236 m_nKaon = kaonList.size();
00237 m_nPi0 = pi0List.size();
00238 m_nKs = ksList.size();
00239
00240
00242
00244
00245 if(m_RdMeasuredEcms&&m_irun!=run){
00246 if(ecmsSvc->isReadDBValid(run)) m_beamE=ecmsSvc->getInfo(run);
00247 if(ecmsSvc->isRunNoValid(run)){
00248 m_beta.setX(ecmsSvc->getPx(run));
00249 m_beta.setY(ecmsSvc->getPy(run));
00250 m_beta.setZ(ecmsSvc->getPz(run));
00251 }
00252 }
00253 log << MSG::INFO << "MeasuredEcmsSvc: " <<m_beamE<< endreq;
00254
00255
00256 if(m_ReadBeamEFromDB && m_irun!=run){
00257 m_irun=run;
00258 if(m_usecalibBeamE)
00259 m_readDb.setcalib(true);
00260 m_beamE=m_readDb.getbeamE(m_irun,m_beamE);
00261 if(run>0)
00262 m_beta=m_readDb.getbeta();
00263
00264 }
00265 double ebeam=m_beamE;
00266
00268
00270
00271
00272 for(int list=0;list<chanlist.size();list++){
00273
00274 string channel=chanlist[list];
00275 vector<int> numchan;
00276 dsSelector.setebeam(ebeam);
00277 dsSelector.setbeta(m_beta);
00278 CDDecayList decaylist(dsSelector);
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 if(channel=="DstoKsK") {
00290 numchan.push_back( EvtRecDTag::kDstoKsK );
00291 numchan.push_back(5);
00292 numchan.push_back(1);
00293 decaylist=ksList* kaonList.plus();
00294 }
00295 else if(channel=="DstoKKPi") {
00296 numchan.push_back( EvtRecDTag::kDstoKKPi );
00297 numchan.push_back(1);
00298 numchan.push_back(1);
00299 numchan.push_back(2);
00300 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus();
00301 }
00302 else if(channel=="DstoKsKPi0") {
00303 numchan.push_back( EvtRecDTag::kDstoKsKPi0 );
00304 numchan.push_back(5);
00305 numchan.push_back(1);
00306 numchan.push_back(3);
00307 decaylist=ksList* kaonList.plus()* pi0List;
00308 }
00309 else if(channel=="DstoKsKsPi") {
00310 numchan.push_back( EvtRecDTag::kDstoKsKsPi );
00311 numchan.push_back(5);
00312 numchan.push_back(5);
00313 numchan.push_back(2);
00314 decaylist=ksList* ksList* pionList.plus();
00315 }
00316 else if(channel=="DstoKKPiPi0") {
00317 numchan.push_back( EvtRecDTag::kDstoKKPiPi0 );
00318 numchan.push_back(1);
00319 numchan.push_back(1);
00320 numchan.push_back(2);
00321 numchan.push_back(3);
00322 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pi0List;
00323 }
00324 else if(channel=="DstoKsKplusPiPi") {
00325 numchan.push_back( EvtRecDTag::kDstoKsKplusPiPi );
00326 numchan.push_back(5);
00327 numchan.push_back(1);
00328 numchan.push_back(2);
00329 numchan.push_back(2);
00330 decaylist=ksList* kaonList.plus()* pionList.plus()* pionList.minus();
00331 }
00332 else if(channel=="DstoKsKminusPiPi") {
00333 numchan.push_back( EvtRecDTag::kDstoKsKminusPiPi );
00334 numchan.push_back(5);
00335 numchan.push_back(1);
00336 numchan.push_back(2);
00337 numchan.push_back(2);
00338 decaylist=ksList* kaonList.minus()* pionList.plus()* pionList.plus();
00339 }
00340 else if(channel=="DstoKKPiPiPi") {
00341 numchan.push_back( EvtRecDTag::kDstoKKPiPiPi );
00342 numchan.push_back(1);
00343 numchan.push_back(1);
00344 numchan.push_back(2);
00345 numchan.push_back(2);
00346 numchan.push_back(2);
00347 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.plus()* pionList.minus();
00348 }
00349 else if(channel=="DstoPiPi0") {
00350 numchan.push_back( EvtRecDTag::kDstoPiPi0 );
00351 numchan.push_back(2);
00352 numchan.push_back(3);
00353 decaylist=pionList.plus()* pi0List;
00354 }
00355 else if(channel=="DstoPiPiPi") {
00356 numchan.push_back( EvtRecDTag::kDstoPiPiPi );
00357 numchan.push_back(2);
00358 numchan.push_back(2);
00359 numchan.push_back(2);
00360 decaylist=pionList.plus()* pionList.plus()* pionList.minus();
00361 }
00362 else if(channel=="DstoPiPiPiPi0") {
00363 numchan.push_back( EvtRecDTag::kDstoPiPiPiPi0 );
00364 numchan.push_back(2);
00365 numchan.push_back(2);
00366 numchan.push_back(2);
00367 numchan.push_back(3);
00368 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
00369 }
00370 else if(channel=="DstoPiPiPiPiPi") {
00371 numchan.push_back( EvtRecDTag::kDstoPiPiPiPiPi );
00372 numchan.push_back(2);
00373 numchan.push_back(2);
00374 numchan.push_back(2);
00375 numchan.push_back(2);
00376 numchan.push_back(2);
00377 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
00378 }
00379 else if(channel=="DstoPiPiPiPiPiPi0") {
00380 numchan.push_back( EvtRecDTag::kDstoPiPiPiPiPiPi0 );
00381 numchan.push_back(2);
00382 numchan.push_back(2);
00383 numchan.push_back(2);
00384 numchan.push_back(2);
00385 numchan.push_back(2);
00386 numchan.push_back(3);
00387 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
00388 }
00389 else if(channel=="DstoPiPiPiPi0Pi0") {
00390 numchan.push_back( EvtRecDTag::kDstoPiPiPiPi0Pi0 );
00391 numchan.push_back(2);
00392 numchan.push_back(2);
00393 numchan.push_back(2);
00394 numchan.push_back(3);
00395 numchan.push_back(3);
00396 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List* pi0List;
00397 }
00398 else if(channel=="DstoPiEta") {
00399 numchan.push_back( EvtRecDTag::kDstoPiEta );
00400 numchan.push_back(2);
00401 numchan.push_back(4);
00402 decaylist=pionList.plus()* etaList;
00403 }
00404 else if(channel=="DstoPiEtaPiPiPi0") {
00405 numchan.push_back( EvtRecDTag::kDstoPiEtaPiPiPi0 );
00406 numchan.push_back(2);
00407 numchan.push_back(9);
00408 CDDecayList EtaList(etatoPiPiPi0Selector);
00409 EtaList=pionList.plus()* pionList.minus()* pi0List;
00410 decaylist=pionList.plus()* EtaList;
00411 }
00412 else if(channel=="DstoPiPi0Eta") {
00413 numchan.push_back( EvtRecDTag::kDstoPiPi0Eta );
00414 numchan.push_back(2);
00415 numchan.push_back(3);
00416 numchan.push_back(4);
00417 decaylist=pionList.plus()* pi0List* etaList;
00418 }
00419 else if(channel=="DstoPiPi0EtaPiPiPi0") {
00420 numchan.push_back( EvtRecDTag::kDstoPiPi0EtaPiPiPi0 );
00421 numchan.push_back(2);
00422 numchan.push_back(3);
00423 numchan.push_back(9);
00424 CDDecayList EtaList(etatoPiPiPi0Selector);
00425 EtaList=pionList.plus()* pionList.minus()* pi0List;
00426 decaylist=pionList.plus()* pi0List* EtaList;
00427 }
00428 else if(channel=="DstoPiPiPiEta") {
00429 numchan.push_back( EvtRecDTag::kDstoPiPiPiEta );
00430 numchan.push_back(2);
00431 numchan.push_back(2);
00432 numchan.push_back(2);
00433 numchan.push_back(4);
00434 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* etaList;
00435 }
00436 else if(channel=="DstoPiPiPiEtaPiPiPi0") {
00437 numchan.push_back( EvtRecDTag::kDstoPiPiPiEtaPiPiPi0 );
00438 numchan.push_back(2);
00439 numchan.push_back(2);
00440 numchan.push_back(2);
00441 numchan.push_back(9);
00442 CDDecayList EtaList(etatoPiPiPi0Selector);
00443 EtaList=pionList.plus()* pionList.minus()* pi0List;
00444 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* EtaList;
00445 }
00446 else if(channel=="DstoPiEPPiPiEta") {
00447 numchan.push_back( EvtRecDTag::kDstoPiEPPiPiEta );
00448 numchan.push_back(2);
00449 numchan.push_back(6);
00450 CDDecayList epList(eptoPiPiEtaSelector);
00451 epList=pionList.plus()* pionList.minus()* etaList;
00452 decaylist=pionList.plus()* epList;
00453 }
00454 else if(channel=="DstoPiPi0EPPiPiEta") {
00455 numchan.push_back( EvtRecDTag::kDstoPiPi0EPPiPiEta );
00456 numchan.push_back(2);
00457 numchan.push_back(3);
00458 numchan.push_back(6);
00459 CDDecayList epList(eptoPiPiEtaSelector);
00460 epList=pionList.plus()* pionList.minus()* etaList;
00461 decaylist=pionList.plus()* pi0List* epList;
00462 }
00463 else if(channel=="DstoPiEPPiPiEtaPiPiPi0") {
00464 numchan.push_back( EvtRecDTag::kDstoPiEPPiPiEtaPiPiPi0 );
00465 numchan.push_back(2);
00466 numchan.push_back(8);
00467 CDDecayList EtaList(etatoPiPiPi0Selector);
00468 EtaList=pionList.plus()* pionList.minus()* pi0List;
00469 CDDecayList EtapList(eptoPiPiEta3PiSelector);
00470 EtapList=pionList.plus()* pionList.minus()* EtaList;
00471 decaylist=pionList.plus()* EtapList;
00472 }
00473 else if(channel=="DstoPiPi0EPPiPiEtaPiPiPi0") {
00474 numchan.push_back( EvtRecDTag::kDstoPiPi0EPPiPiEtaPiPiPi0 );
00475 numchan.push_back(2);
00476 numchan.push_back(3);
00477 numchan.push_back(8);
00478 CDDecayList EtaList(etatoPiPiPi0Selector);
00479 EtaList=pionList.plus()* pionList.minus()* pi0List;
00480 CDDecayList EtapList(eptoPiPiEta3PiSelector);
00481 EtapList=pionList.plus()* pionList.minus()* EtaList;
00482 decaylist=pionList.plus()* pi0List* EtapList;
00483 }
00484
00485 else if(channel=="DstoPiEPRhoGam") {
00486 numchan.push_back( EvtRecDTag::kDstoPiEPRhoGam );
00487 numchan.push_back(2);
00488 numchan.push_back(7);
00489 CDDecayList rhoList(rhotoPiPiSelector);
00490 rhoList=pionList.plus()* pionList.minus();
00491 CDDecayList epList(eptoRhoGamSelector);
00492 epList=rhoList* photonList;
00493 decaylist=pionList.plus()* epList;
00494 }
00495 else if(channel=="DstoPiPi0EPRhoGam") {
00496 numchan.push_back( EvtRecDTag::kDstoPiPi0EPRhoGam );
00497 numchan.push_back(2);
00498 numchan.push_back(3);
00499 numchan.push_back(7);
00500 CDDecayList rhoList(rhotoPiPiSelector);
00501 rhoList=pionList.plus()* pionList.minus();
00502 CDDecayList epList(eptoRhoGamSelector);
00503 epList=rhoList* photonList;
00504 decaylist=pionList.plus()* pi0List* epList;
00505 }
00506 else if(channel=="DstoKsPi") {
00507 numchan.push_back( EvtRecDTag::kDstoKsPi );
00508 numchan.push_back(5);
00509 numchan.push_back(2);
00510 decaylist=ksList* pionList.plus();
00511 }
00512 else if(channel=="DstoKsPiPi0") {
00513 numchan.push_back( EvtRecDTag::kDstoKsPiPi0 );
00514 numchan.push_back(5);
00515 numchan.push_back(2);
00516 numchan.push_back(3);
00517 decaylist=ksList* pionList.plus()* pi0List;
00518 }
00519 else if(channel=="DstoKPiPi") {
00520 numchan.push_back( EvtRecDTag::kDstoKPiPi );
00521 numchan.push_back(1);
00522 numchan.push_back(2);
00523 numchan.push_back(2);
00524 decaylist=kaonList.plus()* pionList.plus()* pionList.minus();
00525 }
00526 else if(channel=="DstoKPiPiPi0") {
00527 numchan.push_back( EvtRecDTag::kDstoKPiPiPi0 );
00528 numchan.push_back(1);
00529 numchan.push_back(2);
00530 numchan.push_back(2);
00531 numchan.push_back(3);
00532 decaylist=kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
00533 }
00534 else if(channel=="DstoKKK") {
00535 numchan.push_back( EvtRecDTag::kDstoKKK );
00536 numchan.push_back(1);
00537 numchan.push_back(1);
00538 numchan.push_back(1);
00539 decaylist=kaonList.minus()* kaonList.plus()* kaonList.plus();
00540 }
00541
00542 CDDecayList::iterator D_begin =decaylist.particle_begin();
00543 CDDecayList::iterator D_end =decaylist.particle_end();
00544
00545 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
00546
00547 EvtRecDTag* recDTag = new EvtRecDTag;
00548 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
00549
00550 vector<int> trackid, showerid;
00551 vector<int> kaonid, pionid;
00552 int numofchildren=numchan.size()-1;
00553
00554 for(int i=0; i< numofchildren;i++){
00555
00556 const CDCandidate& daughter=(*it).particle().child(i);
00557
00558 if(numchan[i+1]==1){
00559 const EvtRecTrack* track=daughter.track();
00560 trackid.push_back(track->trackId());
00561 kaonid.push_back(track->trackId());
00562 }
00563 else if(numchan[i+1]==2){
00564 const EvtRecTrack* track=daughter.track();
00565 trackid.push_back(track->trackId());
00566 pionid.push_back(track->trackId());
00567 }
00568 else if ( numchan[i+1]==3){
00569 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
00570 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
00571 showerid.push_back(hiEnGamma->trackId());
00572 showerid.push_back(loEnGamma->trackId());
00573 }
00574 else if ( numchan[i+1]==4){
00575 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
00576 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
00577 showerid.push_back(hiEnGamma->trackId());
00578 showerid.push_back(loEnGamma->trackId());
00579 }
00580 else if ( numchan[i+1]==5){
00581 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
00582
00583
00584 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
00585
00586 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
00587 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
00588 trackid.push_back(pion1Trk->trackId());
00589 trackid.push_back(pion2Trk->trackId());
00590 }
00591 else if (numchan[i+1]==6){
00592 const CDCandidate& apion = daughter.decay().child(0);
00593 const CDCandidate& spion = daughter.decay().child(1);
00594 const CDCandidate& eta = daughter.decay().child(2);
00595 const EvtRecTrack* apiontrk = apion.track();
00596 const EvtRecTrack* spiontrk = spion.track();
00597 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
00598 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
00599
00600 trackid.push_back(apiontrk->trackId());
00601 trackid.push_back(spiontrk->trackId());
00602 showerid.push_back(hiEnGamma->trackId());
00603 showerid.push_back(loEnGamma->trackId());
00604
00605 }
00606 else if (numchan[i+1]==7){
00607 const CDCandidate& rho = daughter.decay().child(0);
00608 const CDCandidate& apion = rho.decay().child(0);
00609 const CDCandidate& spion = rho.decay().child(1);
00610 const CDCandidate& gamma = daughter.decay().child(1);
00611
00612 const EvtRecTrack* apiontrk = apion.track();
00613 const EvtRecTrack* spiontrk = spion.track();
00614 const EvtRecTrack* gammatrk = gamma.photon();
00615
00616
00617 trackid.push_back(apiontrk->trackId());
00618 trackid.push_back(spiontrk->trackId());
00619 showerid.push_back(gammatrk->trackId());
00620
00621 }
00622 else if (numchan[i+1]==8){
00623 const CDCandidate& apion = daughter.decay().child(0);
00624 const CDCandidate& spion = daughter.decay().child(1);
00625 const CDCandidate& eta = daughter.decay().child(2);
00626 const CDCandidate& e0pion = eta.decay().child(0);
00627 const CDCandidate& e1pion = eta.decay().child(1);
00628 const CDCandidate& pi0 = eta.decay().child(2);
00629
00630 const EvtRecTrack* apiontrk = apion.track();
00631 const EvtRecTrack* spiontrk = spion.track();
00632 const EvtRecTrack* e0piontrk = e0pion.track();
00633 const EvtRecTrack* e1piontrk = e1pion.track();
00634
00635 const EvtRecTrack* hiEnGamma=pi0.navPi0()->hiEnGamma();
00636 const EvtRecTrack* loEnGamma=pi0.navPi0()->loEnGamma();
00637
00638 trackid.push_back(apiontrk->trackId());
00639 trackid.push_back(spiontrk->trackId());
00640 trackid.push_back(e0piontrk->trackId());
00641 trackid.push_back(e1piontrk->trackId());
00642 showerid.push_back(hiEnGamma->trackId());
00643 showerid.push_back(loEnGamma->trackId());
00644
00645 }
00646 else if (numchan[i+1]==9){
00647 const CDCandidate& apion = daughter.decay().child(0);
00648 const CDCandidate& spion = daughter.decay().child(1);
00649 const CDCandidate& pi0 = daughter.decay().child(2);
00650 const EvtRecTrack* apiontrk = apion.track();
00651 const EvtRecTrack* spiontrk = spion.track();
00652 const EvtRecTrack* hiEnGamma=pi0.navPi0()->hiEnGamma();
00653 const EvtRecTrack* loEnGamma=pi0.navPi0()->loEnGamma();
00654
00655 trackid.push_back(apiontrk->trackId());
00656 trackid.push_back(spiontrk->trackId());
00657 showerid.push_back(hiEnGamma->trackId());
00658 showerid.push_back(loEnGamma->trackId());
00659
00660 }
00661
00662
00663 }
00664
00665
00666 saveDsInfo(it, ebeam, numofchildren, recDTag);
00667 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
00668 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
00669
00670
00671 if(m_usevertexfit){
00672 HepLorentzVector p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
00673 recDTag->setp4(recDTag->p4()+p4change_vfit);
00674 }
00675
00676
00677 trackid.clear();
00678 showerid.clear();
00679 kaonid.clear();
00680 pionid.clear();
00681
00682
00683
00684
00685 recDTagCol->push_back(recDTag);
00686
00687 }
00688
00689 numchan.clear();
00690
00691 }
00692
00693
00694
00695 return StatusCode::SUCCESS;
00696 }
00697
00698
00699 void DsReconstruction::saveDsInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag){
00700
00701 double mass = (*it).particle().mass();
00702 int charge= (*it).particle().charge();
00703 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
00704 recDTag->setp4(p4);
00705
00706 p4.boost(-m_beta);
00707 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
00708 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
00709 double deltaE_CMS = p4.t() - ebeam;
00710
00711 if((*it).particle().userTag()==1)
00712 recDTag->settype( EvtRecDTag::Tight );
00713 else
00714 recDTag->settype( EvtRecDTag::Loose );
00715 recDTag->setcharge(charge);
00716 recDTag->setcharm(charge);
00717 recDTag->setnumOfChildren(numofchildren);
00718 recDTag->setmass(mass);
00719 recDTag->setmBC(mbc_CMS);
00720 recDTag->setbeamE(ebeam);
00721 recDTag->setdeltaE(deltaE_CMS);
00722
00723 }
00724
00725 void DsReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
00726 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
00727
00728 vector<EvtRecTrackIterator> trktemp;
00729 vector<EvtRecTrackIterator> shrtemp;
00730
00731
00732 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
00733
00734 bool isothertrack=true;
00735 for(int i=0; i<trackid.size(); i++){
00736 if( (*trk)->trackId()==trackid[i] ){
00737 trktemp.push_back(trk);
00738 isothertrack=false;
00739 break;
00740 }
00741 }
00742 if(isothertrack)
00743 recDTag->addOtherTrack(*trk);
00744 }
00745 for(int i=0; i<trackid.size();i++){
00746 for(int j=0; j<trktemp.size(); j++){
00747 EvtRecTrackIterator trk=trktemp[j];
00748 if( (*trk)->trackId()==trackid[i])
00749 recDTag->addTrack(*trktemp[j]);
00750 }
00751 }
00752
00753
00754
00755 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
00756 bool isothershower=true;
00757 for(int i=0; i<showerid.size(); i++){
00758 if( (*shr)->trackId()==showerid[i] ){
00759 shrtemp.push_back(shr);
00760 isothershower=false;
00761 break;
00762 }
00763 }
00764 if(isothershower)
00765 recDTag->addOtherShower(*shr);
00766 }
00767
00768 for(int i=0; i<showerid.size();i++){
00769 for(int j=0; j<shrtemp.size(); j++){
00770 EvtRecTrackIterator shr=shrtemp[j];
00771 if( (*shr)->trackId()==showerid[i])
00772 recDTag->addShower(*shrtemp[j]);
00773 }
00774 }
00775
00776
00777 }
00778
00779 void DsReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
00780
00781 bool iskaon=false,ispion=false;
00782
00783
00784
00785
00786 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
00787 recDTag->addKaonId( (*kit).particle().track() );
00788 }
00789
00790 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
00791 recDTag->addPionId( (*pit).particle().track() );
00792 }
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831 }
00832
00833
00834 vector<string> DsReconstruction::getlist(string& filename ){
00835
00836 string channel;
00837 vector<string> temp;
00838
00839 ifstream inFile;
00840
00841 inFile.open(filename.c_str());
00842 if (!inFile) {
00843 cout << "Unable to open decay list file";
00844 exit(1);
00845 }
00846
00847 while (inFile >> channel) {
00848 temp.push_back(channel);
00849 }
00850
00851 inFile.close();
00852
00853 return temp;
00854
00855 }