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
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 #include "VertexFit/IVertexDbSvc.h"
00036
00037 #include "DTagAlg/ChargedDReconstruction.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/LocalEptoPiPiEtaSelector.h"
00045 #include "DTagAlg/LocalEptoRhoGamSelector.h"
00046 #include "DTagAlg/LocalRhotoPiPiSelector.h"
00047 #include "DTagAlg/ChargedDSelector.h"
00048
00049
00050 #include "DTagAlg/utility.h"
00051
00052 using namespace Event;
00053
00054
00055 ChargedDReconstruction::ChargedDReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
00056 Algorithm(name, pSvcLocator) {
00057
00058 declareProperty( "debug", m_debug = false );
00059 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
00060 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
00061 declareProperty( "UseVertexfit", m_usevertexfit = false );
00062 declareProperty( "BeamE", m_beamE = 1.8865 );
00063 declareProperty( "DpList", m_decaylist = "test.txt" );
00064 declareProperty("ReadMeasuredEcms", m_RdMeasuredEcms = false);
00065 }
00066
00067
00068 StatusCode ChargedDReconstruction::initialize() {
00069 MsgStream log(msgSvc(), name());
00070 log << MSG::INFO << "in initialize()" <<endreq;
00071
00072 m_irun=-100;
00073 m_beta.setX(0.011);
00074 m_beta.setY(0);
00075 m_beta.setZ(0);
00076 chanlist=getlist(m_decaylist);
00077
00078
00079 if(m_RdMeasuredEcms){
00080 StatusCode status=serviceLocator()->service("MeasuredEcmsSvc", ecmsSvc, true);
00081 if(!status.isSuccess()){
00082 std::cout<<"ERROR: Can not initial the IMeasuredEcmsSvc right"<<std::endl;
00083 return status;
00084 }
00085 }
00086
00087
00088 return StatusCode::SUCCESS;
00089 }
00090
00091
00092 StatusCode ChargedDReconstruction::finalize() {
00093 MsgStream log(msgSvc(), name());
00094 log << MSG::INFO << "in finalize()" << endreq;
00095
00096 chanlist.clear();
00097
00098 return StatusCode::SUCCESS;
00099 }
00100
00101 StatusCode ChargedDReconstruction::registerEvtRecDTagCol(
00102 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
00103 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
00104 aNewEvtRecDTagCol);
00105 if (sc != StatusCode::SUCCESS) {
00106 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
00107 }
00108 return sc;
00109 }
00110
00111
00112 StatusCode ChargedDReconstruction::execute() {
00113 MsgStream log(msgSvc(), name());
00114 log << MSG::INFO << "in execute()" << endreq;
00115
00116 StatusCode sc;
00117
00119
00121 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00122 int event= eventHeader->eventNumber();
00123
00124
00125
00126 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00127 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00128 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
00129 << " " << eventHeader->eventNumber() << endreq;
00130 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
00131 << recEvent->totalCharged() << " , "
00132 << recEvent->totalNeutral() << " , "
00133 << recEvent->totalTracks() <<endreq;
00134
00135 EvtRecTrackIterator charged_begin = recTrackCol->begin();
00136 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
00137
00138 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
00139 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
00140
00141
00142 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
00143 if ( ! recPi0Col ) {
00144 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
00145 return StatusCode::FAILURE;
00146 }
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
00172 Hep3Vector xorigin(0,0,0);
00173
00174 IVertexDbSvc* vtxsvc;
00175 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
00176 if (vtxsvc->isVertexValid()) {
00177
00178
00179 double* vertex = vtxsvc->PrimaryVertex();
00180 xorigin.setX(vertex[0]);
00181 xorigin.setY(vertex[1]);
00182 xorigin.setZ(vertex[2]);
00183 }
00184
00185 utility util;
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00200
00202 pionSelector.setpidtype(0);
00203 kaonSelector.setpidtype(0);
00204 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
00205 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
00206
00207 CDKsList ksList(ksSelector);
00208 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
00209
00210
00211 map<EvtRecVeeVertex*, vector< double > > fitinfo;
00212 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
00213 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
00214 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
00215 }
00216
00217 CDPi0List pi0List(pi0Selector);
00218 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
00219
00220 CDEtaList etaList(etatoGGSelector);
00221 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
00222
00223
00224 pionSelector.setpidtype(1);
00225 kaonSelector.setpidtype(1);
00226 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
00227 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
00228
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
00246
00247 if(m_RdMeasuredEcms&&m_irun!=run){
00248 if(ecmsSvc->isReadDBValid(run)) m_beamE=ecmsSvc->getInfo(run);
00249 if(ecmsSvc->isRunNoValid(run)){
00250 m_beta.setX(ecmsSvc->getPx(run));
00251 m_beta.setY(ecmsSvc->getPy(run));
00252 m_beta.setZ(ecmsSvc->getPz(run));
00253 }
00254 }
00255 log << MSG::INFO << "MeasuredEcmsSvc: " <<m_beamE<< endreq;
00256
00257
00258 if(m_ReadBeamEFromDB && m_irun!=run){
00259 m_irun=run;
00260 if(m_usecalibBeamE)
00261 m_readDb.setcalib(true);
00262 m_beamE=m_readDb.getbeamE(m_irun,m_beamE);
00263 if(run>0)
00264 m_beta=m_readDb.getbeta();
00265
00266 }
00267 double ebeam=m_beamE;
00268
00269
00271
00273
00274
00275 for(int list=0;list<chanlist.size();list++){
00276
00277 string channel=chanlist[list];
00278 vector<int> numchan;
00279 chargedDSelector.setebeam(ebeam);
00280 chargedDSelector.setbeta(m_beta);
00281 CDDecayList decaylist(chargedDSelector);
00282
00283
00284
00285
00286
00287 if(channel=="DptoKPiPi") {
00288 numchan.push_back( EvtRecDTag::kDptoKPiPi );
00289 numchan.push_back(1);
00290 numchan.push_back(2);
00291 numchan.push_back(2);
00292 decaylist=kaonList.minus()* pionList.plus()* pionList.plus();
00293 }
00294 else if(channel=="DptoKPiPiPi0") {
00295 numchan.push_back( EvtRecDTag::kDptoKPiPiPi0 );
00296 numchan.push_back(1);
00297 numchan.push_back(2);
00298 numchan.push_back(2);
00299 numchan.push_back(3);
00300 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pi0List;
00301 }
00302 else if(channel=="DptoKsPi") {
00303 numchan.push_back( EvtRecDTag::kDptoKsPi );
00304 numchan.push_back(5);
00305 numchan.push_back(2);
00306 decaylist=ksList* pionList.plus();
00307 }
00308 else if(channel=="DptoKsPiPi0") {
00309 numchan.push_back( EvtRecDTag::kDptoKsPiPi0 );
00310 numchan.push_back(5);
00311 numchan.push_back(2);
00312 numchan.push_back(3);
00313 decaylist=ksList* pionList.plus()* pi0List;
00314 }
00315 else if(channel=="DptoKsPiPiPi") {
00316 numchan.push_back( EvtRecDTag::kDptoKsPiPiPi );
00317 numchan.push_back(5);
00318 numchan.push_back(2);
00319 numchan.push_back(2);
00320 numchan.push_back(2);
00321 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus();
00322 }
00323 else if(channel=="DptoKKPi") {
00324 numchan.push_back( EvtRecDTag::kDptoKKPi );
00325 numchan.push_back(1);
00326 numchan.push_back(1);
00327 numchan.push_back(2);
00328 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus();
00329 }
00330 else if(channel=="DptoPiPi0") {
00331 numchan.push_back( EvtRecDTag::kDptoPiPi0 );
00332 numchan.push_back(2);
00333 numchan.push_back(3);
00334 decaylist=pionList.plus()* pi0List;
00335 }
00336 else if(channel=="DptoKPi0") {
00337 numchan.push_back( EvtRecDTag::kDptoKPi0 );
00338 numchan.push_back(1);
00339 numchan.push_back(3);
00340 decaylist=kaonList.plus()* pi0List;
00341 }
00342 else if(channel=="DptoKsK") {
00343 numchan.push_back( EvtRecDTag::kDptoKsK );
00344 numchan.push_back(5);
00345 numchan.push_back(1);
00346 decaylist=ksList* kaonList.plus();
00347 }
00348 else if(channel=="DptoPiPiPi") {
00349 numchan.push_back( EvtRecDTag::kDptoPiPiPi );
00350 numchan.push_back(2);
00351 numchan.push_back(2);
00352 numchan.push_back(2);
00353 decaylist=pionList.plus()* pionList.plus()* pionList.minus() ;
00354 }
00355 else if(channel=="DptoPiPi0Pi0") {
00356 numchan.push_back( EvtRecDTag::kDptoPiPi0Pi0 );
00357 numchan.push_back(2);
00358 numchan.push_back(3);
00359 numchan.push_back(3);
00360 decaylist=pionList.plus()* pi0List* pi0List;
00361 }
00362 else if(channel=="DptoKsKsPi") {
00363 numchan.push_back( EvtRecDTag::kDptoKsKsPi );
00364 numchan.push_back(5);
00365 numchan.push_back(5);
00366 numchan.push_back(2);
00367 decaylist=ksList* ksList* pionList.plus();
00368 }
00369 else if(channel=="DptoKsKPi0") {
00370 numchan.push_back( EvtRecDTag::kDptoKsKPi0 );
00371 numchan.push_back(5);
00372 numchan.push_back(1);
00373 numchan.push_back(3);
00374 decaylist=ksList* kaonList.plus()* pi0List;
00375 }
00376 else if(channel=="DptoKsKsK") {
00377 numchan.push_back( EvtRecDTag::kDptoKsKsK );
00378 numchan.push_back(5);
00379 numchan.push_back(5);
00380 numchan.push_back(1);
00381 decaylist=ksList* ksList* kaonList.plus();
00382 }
00383 else if(channel=="DptoPiPiPiPi0") {
00384 numchan.push_back( EvtRecDTag::kDptoPiPiPiPi0 );
00385 numchan.push_back(2);
00386 numchan.push_back(2);
00387 numchan.push_back(2);
00388 numchan.push_back(3);
00389 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
00390 }
00391 else if(channel=="DptoKsPiPi0Pi0") {
00392 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Pi0 );
00393 numchan.push_back(5);
00394 numchan.push_back(2);
00395 numchan.push_back(3);
00396 numchan.push_back(3);
00397 decaylist=ksList* pionList.plus()* pi0List* pi0List;
00398 }
00399 else if(channel=="DptoKsKplusPiPi") {
00400 numchan.push_back( EvtRecDTag::kDptoKsKplusPiPi );
00401 numchan.push_back(5);
00402 numchan.push_back(1);
00403 numchan.push_back(2);
00404 numchan.push_back(2);
00405 decaylist=ksList* kaonList.plus()* pionList.plus()* pionList.minus();
00406 }
00407 else if(channel=="DptoKsKminusPiPi") {
00408 numchan.push_back( EvtRecDTag::kDptoKsKminusPiPi );
00409 numchan.push_back(5);
00410 numchan.push_back(1);
00411 numchan.push_back(2);
00412 numchan.push_back(2);
00413 decaylist=ksList* kaonList.minus()* pionList.plus()* pionList.plus();
00414 }
00415 else if(channel=="DptoKKPiPi0") {
00416 numchan.push_back( EvtRecDTag::kDptoKKPiPi0 );
00417 numchan.push_back(1);
00418 numchan.push_back(1);
00419 numchan.push_back(2);
00420 numchan.push_back(3);
00421 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pi0List;
00422 }
00423 else if(channel=="DptoPiPiPiPiPi") {
00424 numchan.push_back( EvtRecDTag::kDptoPiPiPiPiPi );
00425 numchan.push_back(2);
00426 numchan.push_back(2);
00427 numchan.push_back(2);
00428 numchan.push_back(2);
00429 numchan.push_back(2);
00430 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
00431 }
00432 else if(channel=="DptoKPiPiPiPi") {
00433 numchan.push_back( EvtRecDTag::kDptoKPiPiPiPi );
00434 numchan.push_back(1);
00435 numchan.push_back(2);
00436 numchan.push_back(2);
00437 numchan.push_back(2);
00438 numchan.push_back(2);
00439 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus();
00440 }
00441 else if(channel=="DptoPiEta") {
00442 numchan.push_back( EvtRecDTag::kDptoPiEta );
00443 numchan.push_back(2);
00444 numchan.push_back(4);
00445 decaylist=pionList.plus()* etaList;
00446 }
00447 else if(channel=="DptoKsPiEta") {
00448 numchan.push_back( EvtRecDTag::kDptoKsPiEta );
00449 numchan.push_back(5);
00450 numchan.push_back(2);
00451 numchan.push_back(4);
00452 decaylist=ksList* pionList.plus()* etaList;
00453 }
00454
00455 CDDecayList::iterator D_begin =decaylist.particle_begin();
00456 CDDecayList::iterator D_end =decaylist.particle_end();
00457
00458 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
00459
00460 EvtRecDTag* recDTag = new EvtRecDTag;
00461 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
00462
00463 vector<int> trackid, showerid;
00464 vector<int> kaonid, pionid;
00465 int numofchildren=numchan.size()-1;
00466
00467 for(int i=0; i< numofchildren;i++){
00468
00469 const CDCandidate& daughter=(*it).particle().child(i);
00470
00471 if(numchan[i+1]==1){
00472 const EvtRecTrack* track=daughter.track();
00473 trackid.push_back(track->trackId());
00474 kaonid.push_back(track->trackId());
00475 }
00476 else if(numchan[i+1]==2){
00477 const EvtRecTrack* track=daughter.track();
00478 trackid.push_back(track->trackId());
00479 pionid.push_back(track->trackId());
00480 }
00481 else if ( numchan[i+1]==3){
00482 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
00483 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
00484 showerid.push_back(hiEnGamma->trackId());
00485 showerid.push_back(loEnGamma->trackId());
00486 }
00487 else if ( numchan[i+1]==4){
00488 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
00489 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
00490 showerid.push_back(hiEnGamma->trackId());
00491 showerid.push_back(loEnGamma->trackId());
00492 }
00493 else if ( numchan[i+1]==5){
00494 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
00495
00496 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
00497
00498 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
00499 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
00500 trackid.push_back(pion1Trk->trackId());
00501 trackid.push_back(pion2Trk->trackId());
00502 }
00503
00504 }
00505
00506
00507 saveDpInfo(it, ebeam, numofchildren, recDTag);
00508 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
00509 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
00510
00511 if(m_usevertexfit){
00512 HepLorentzVector p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
00513 recDTag->setp4(recDTag->p4()+p4change_vfit);
00514 }
00515
00516
00517 trackid.clear();
00518 showerid.clear();
00519 kaonid.clear();
00520 pionid.clear();
00521
00522
00523 recDTagCol->push_back(recDTag);
00524
00525 }
00526
00527 numchan.clear();
00528
00529 }
00530
00531 return StatusCode::SUCCESS;
00532 }
00533
00534
00535 void ChargedDReconstruction::saveDpInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag){
00536
00537 double mass = (*it).particle().mass();
00538 int charge= (*it).particle().charge();
00539 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
00540 recDTag->setp4(p4);
00541
00542 p4.boost(-m_beta);
00543 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
00544 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
00545 double deltaE_CMS = p4.t() - ebeam;
00546
00547 if((*it).particle().userTag()==1)
00548 recDTag->settype( EvtRecDTag::Tight );
00549 else
00550 recDTag->settype( EvtRecDTag::Loose );
00551 recDTag->setcharge(charge);
00552 recDTag->setcharm(charge);
00553 recDTag->setnumOfChildren(numofchildren);
00554 recDTag->setmass(mass);
00555 recDTag->setmBC(mbc_CMS);
00556 recDTag->setbeamE(ebeam);
00557 recDTag->setdeltaE(deltaE_CMS);
00558
00559 }
00560
00561 void ChargedDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
00562 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
00563
00564 vector<EvtRecTrackIterator> trktemp;
00565 vector<EvtRecTrackIterator> shrtemp;
00566
00567
00568 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
00569
00570 bool isothertrack=true;
00571 for(int i=0; i<trackid.size(); i++){
00572 if( (*trk)->trackId()==trackid[i] ){
00573 trktemp.push_back(trk);
00574 isothertrack=false;
00575 break;
00576 }
00577 }
00578 if(isothertrack)
00579 recDTag->addOtherTrack(*trk);
00580 }
00581 for(int i=0; i<trackid.size();i++){
00582 for(int j=0; j<trktemp.size(); j++){
00583 EvtRecTrackIterator trk=trktemp[j];
00584 if( (*trk)->trackId()==trackid[i])
00585 recDTag->addTrack(*trktemp[j]);
00586 }
00587 }
00588
00589
00590
00591 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
00592 bool isothershower=true;
00593 for(int i=0; i<showerid.size(); i++){
00594 if( (*shr)->trackId()==showerid[i] ){
00595 shrtemp.push_back(shr);
00596 isothershower=false;
00597 break;
00598 }
00599 }
00600 if(isothershower)
00601 recDTag->addOtherShower(*shr);
00602 }
00603
00604 for(int i=0; i<showerid.size();i++){
00605 for(int j=0; j<shrtemp.size(); j++){
00606 EvtRecTrackIterator shr=shrtemp[j];
00607 if( (*shr)->trackId()==showerid[i])
00608 recDTag->addShower(*shrtemp[j]);
00609 }
00610 }
00611
00612
00613 }
00614
00615 void ChargedDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
00616
00617 bool iskaon=false,ispion=false;
00618
00619
00620
00621
00622 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
00623 recDTag->addKaonId( (*kit).particle().track() );
00624 }
00625
00626 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
00627 recDTag->addPionId( (*pit).particle().track() );
00628 }
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667 }
00668
00669
00670
00671 vector<string> ChargedDReconstruction::getlist(string& filename ){
00672
00673 string channel;
00674 vector<string> temp;
00675
00676 ifstream inFile;
00677
00678 inFile.open(filename.c_str());
00679 if (!inFile) {
00680 cout << "Unable to open decay list file";
00681 exit(1);
00682 }
00683
00684 while (inFile >> channel) {
00685 temp.push_back(channel);
00686 }
00687
00688 inFile.close();
00689
00690 return temp;
00691
00692 }