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