BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
NeutralDReconstruction.cxx
Go to the documentation of this file.
1//
2// All D0 decay modes Reconstruction
3//
4
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"
10#include "GaudiKernel/PropertyMgr.h"
11
12#include "EventModel/EventHeader.h"
13#include "EventModel/Event.h"
14#include "EventModel/EventModel.h"
15
16#include "EvtRecEvent/EvtRecEvent.h"
17#include "EvtRecEvent/EvtRecTrack.h"
18#include "EvtRecEvent/EvtRecDTag.h"
19#include "EvtRecEvent/EvtRecPi0.h"
20#include "EvtRecEvent/EvtRecEtaToGG.h"
21#include "EvtRecEvent/EvtRecVeeVertex.h"
22
23#include "BesDChain/CDChargedPionList.h"
24#include "BesDChain/CDChargedKaonList.h"
25#include "BesDChain/CDPhotonList.h"
26#include "BesDChain/CDPi0List.h"
27#include "BesDChain/CDEtaList.h"
28#include "BesDChain/CDKsList.h"
29#include "BesDChain/CDDecayList.h"
30
31#include "McTruth/McParticle.h"
32#include "ParticleID/ParticleID.h"
33#include "VertexFit/VertexFit.h"
34#include "VertexFit/SecondVertexFit.h"
35
36
37#include "DTagAlg/NeutralDReconstruction.h"
38#include "DTagAlg/LocalKaonSelector.h"
39#include "DTagAlg/LocalPionSelector.h"
40#include "DTagAlg/LocalKsSelector.h"
41#include "DTagAlg/LocalPi0Selector.h"
42#include "DTagAlg/LocalEtatoGGSelector.h"
43#include "DTagAlg/LocalPhotonSelector.h"
44#include "DTagAlg/LocalEptoPiPiEtaSelector.h"
45#include "DTagAlg/LocalEptoRhoGamSelector.h"
46#include "DTagAlg/LocalRhotoPiPiSelector.h"
47#include "DTagAlg/NeutralDSelector.h"
48
49#include <fstream>
50
51using namespace Event;
52
53
54//*******************************************************************************************
55NeutralDReconstruction::NeutralDReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
56Algorithm(name, pSvcLocator) {
57 //Declare the properties
58 declareProperty( "debug", m_debug = false );
59 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
60 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
61 declareProperty( "BeamE", m_beamE = 1.8865 );
62 declareProperty( "D0List", m_decaylist = "test.txt" );
63}
64
65//******************************************************************************************
66StatusCode NeutralDReconstruction::initialize() {
67 MsgStream log(msgSvc(), name());
68 log << MSG::INFO << "in initialize()" <<endreq;
69
70 m_irun=-100;
71 m_beta.setX(0.011);
72 m_beta.setY(0);
73 m_beta.setZ(0);
74
75 chanlist=getlist(m_decaylist);
76
77 return StatusCode::SUCCESS;
78}
79
80//********************************************************************************************
81StatusCode NeutralDReconstruction::finalize() {
82 MsgStream log(msgSvc(), name());
83 log << MSG::INFO << "in finalize()" << endreq;
84
85 chanlist.clear();
86
87 return StatusCode::SUCCESS;
88}
89
90StatusCode NeutralDReconstruction::registerEvtRecDTagCol(
91 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
92 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
93 aNewEvtRecDTagCol);
94 if (sc != StatusCode::SUCCESS) {
95 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
96 }
97 return sc;
98}
99
100//*********************************************************************************************
101StatusCode NeutralDReconstruction::execute() {
102 MsgStream log(msgSvc(), name());
103 log << MSG::INFO << "in execute()" << endreq;
104
105 StatusCode sc;
106
107 //////////////////
108 // Read REC data
109 /////////////////
110 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
111 int event= eventHeader->eventNumber();
112 // if ( m_debug || ( (event & 0x3FF) == 0 ) )
113 //std::cout << "event: " << event << std::endl;
114
115 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
116 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
117 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
118 << " " << eventHeader->eventNumber() << endreq;
119 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
120 << recEvent->totalCharged() << " , "
121 << recEvent->totalNeutral() << " , "
122 << recEvent->totalTracks() <<endreq;
123
124 EvtRecTrackIterator charged_begin = recTrackCol->begin();
125 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
126
127 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
128 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
129
130
131 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
132 if ( ! recPi0Col ) {
133 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
134 return StatusCode::FAILURE;
135 }
136
137
138
139 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
140 if ( ! recEtaToGGCol ) {
141 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
142 return StatusCode::FAILURE;
143 }
144
145
146 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
147 if ( ! recVeeVertexCol ) {
148 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
149 return StatusCode::FAILURE;
150 }
151
152
153 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
154 if (!recDTagCol) {
155 log << MSG::FATAL << "EvtRecDTagCol is not registered in TDS" << endreq;
156 return StatusCode::FAILURE;
157 }
158
159 //registered in the parent file DTag.cxx
160 /*
161 if (!recDTagCol) {
162 recDTagCol = new EvtRecDTagCol;
163 sc = registerEvtRecDTagCol(recDTagCol, log);
164 if (sc != StatusCode::SUCCESS) {
165 return sc;
166 }
167 }
168 */
169
170 /////////////////////////////
171 //reconstruct particle lists
172 /////////////////////////////
173
174 pionSelector.setpidtype(0);
175 kaonSelector.setpidtype(0);
176 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
177 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
178 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
179 CDKsList ksList(ksSelector);
180 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
181
182 CDPi0List pi0List(pi0Selector);
183 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
184
185 CDEtaList etaList(etatoGGSelector);
186 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
187
188
189 //pion/kaon list with PID
190 pionSelector.setpidtype(1);
191 kaonSelector.setpidtype(1);
192 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
193 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
194
195 int run = eventHeader->runNumber();
196 m_ievt = eventHeader->eventNumber();
197 m_nChrg = recEvent->totalCharged();
198 m_nNeu = recEvent->totalNeutral();
199 m_nPion = pionList.size();
200 m_nKaon = kaonList.size();
201 m_nPi0 = pi0List.size();
202 m_nKs = ksList.size();
203
204 ///////////////////////
205 // get beam energy and beta
206 ///////////////////////
207
208 if(m_ReadBeamEFromDB && m_irun!=run){
209 m_irun=run;
210 if(m_usecalibBeamE)
211 m_readDb.setcalib(true);
212 m_beamE=m_readDb.getbeamE(m_irun, m_beamE);
213 if(run>0)
214 m_beta=m_readDb.getbeta();
215
216 //cout<<"use beam E from data base:"<<m_beamE<<endl;
217 }
218 double ebeam=m_beamE;
219
220 //////////////////////////////
221 //reconstruct decay lists
222 /////////////////////////////
223
224
225
226 for(int list=0;list<chanlist.size();list++){
227
228 string channel=chanlist[list];
229 vector<int> numchan;
230 neutralDSelector.setebeam(ebeam);
231 neutralDSelector.setbeta(m_beta);
232 CDDecayList decaylist(neutralDSelector);
233
234 bool isFlavorMode=false;
235
236 //K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5, Eta'(pipiEta):6, Eta'(rhogam):7
237 //the fist element of the vector stands for decay mode,
238 //the rest will be particles, and size of the vector minus 1 will be number of daughers.
239 //the first particle in the vector will be used to get charm
240
241 if(channel=="D0toKPi") {
242 numchan.push_back( EvtRecDTag::kD0toKPi );
243 numchan.push_back(1);
244 numchan.push_back(2);
245 decaylist=kaonList.minus() * pionList.plus();
246 isFlavorMode=true;
247 }
248 else if(channel=="D0toKPiPi0"){
249 numchan.push_back( EvtRecDTag::kD0toKPiPi0);
250 numchan.push_back(1);
251 numchan.push_back(2);
252 numchan.push_back(3);
253 decaylist=kaonList.minus() * pionList.plus() * pi0List;
254 isFlavorMode=true;
255 }
256 else if(channel=="D0toKPiPi0Pi0"){
257 numchan.push_back( EvtRecDTag::kD0toKPiPi0Pi0);
258 numchan.push_back(1);
259 numchan.push_back(2);
260 numchan.push_back(3);
261 numchan.push_back(3);
262 decaylist=kaonList.minus() * pionList.plus() * pi0List * pi0List;
263 isFlavorMode=true;
264 }
265 else if(channel=="D0toKPiPiPi"){
266 numchan.push_back( EvtRecDTag::kD0toKPiPiPi);
267 numchan.push_back(1);
268 numchan.push_back(2);
269 numchan.push_back(2);
270 numchan.push_back(2);
271 decaylist=kaonList.minus()* pionList.plus()* pionList.plus() * pionList.minus();
272 isFlavorMode=true;
273 }
274 else if(channel=="D0toKPiPiPiPi0"){
275 numchan.push_back( EvtRecDTag::kD0toKPiPiPiPi0);
276 numchan.push_back(1);
277 numchan.push_back(2);
278 numchan.push_back(2);
279 numchan.push_back(2);
280 numchan.push_back(3);
281 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
282 isFlavorMode=true;
283 }
284 else if(channel=="D0toKPiEta"){
285 numchan.push_back( EvtRecDTag::kD0toKPiEta);
286 numchan.push_back(1);
287 numchan.push_back(2);
288 numchan.push_back(4);
289 decaylist=kaonList.minus() * pionList.plus() * etaList;
290 isFlavorMode=true;
291 }
292 else if(channel=="D0toKsKPi"){
293 numchan.push_back( EvtRecDTag::kD0toKsKPi);
294 numchan.push_back(5);
295 numchan.push_back(1);
296 numchan.push_back(2);
297 decaylist=ksList* kaonList.minus()* pionList.plus();
298 }
299 else if(channel=="D0toKsKPiPi0"){
300 numchan.push_back( EvtRecDTag::kD0toKsKPiPi0);
301 numchan.push_back(5);
302 numchan.push_back(1);
303 numchan.push_back(2);
304 numchan.push_back(3);
305 decaylist=ksList* kaonList.minus() * pionList.plus()* pi0List;
306 }
307 else if(channel=="D0toKsPiPi"){
308 numchan.push_back( EvtRecDTag::kD0toKsPiPi);
309 numchan.push_back(5);
310 numchan.push_back(2);
311 numchan.push_back(2);
312 decaylist=ksList* pionList.plus()* pionList.minus();
313 }
314 else if(channel=="D0toKsPiPiPi0"){
315 numchan.push_back( EvtRecDTag::kD0toKsPiPiPi0);
316 numchan.push_back(5);
317 numchan.push_back(2);
318 numchan.push_back(2);
319 numchan.push_back(3);
320 decaylist=ksList* pionList.plus()* pionList.minus()* pi0List;
321 }
322 else if(channel=="D0toKsPi0"){
323 numchan.push_back( EvtRecDTag::kD0toKsPi0);
324 numchan.push_back(5);
325 numchan.push_back(3);
326 decaylist=ksList* pi0List;
327 }
328 else if(channel=="D0toPiPiPi0"){
329 numchan.push_back( EvtRecDTag::kD0toPiPiPi0);
330 numchan.push_back(2);
331 numchan.push_back(2);
332 numchan.push_back(3);
333 decaylist=pionList.plus()* pionList.minus()* pi0List;
334 }
335 else if(channel=="D0toPiPi"){
336 numchan.push_back( EvtRecDTag::kD0toPiPi);
337 numchan.push_back(2);
338 numchan.push_back(2);
339 decaylist=pionList.plus()* pionList.minus();
340 }
341 else if(channel=="D0toKK"){
342 numchan.push_back( EvtRecDTag::kD0toKK);
343 numchan.push_back(1);
344 numchan.push_back(1);
345 decaylist=kaonList.minus()* kaonList.plus();
346 }
347 else if(channel=="D0toKKPi0"){
348 numchan.push_back( EvtRecDTag::kD0toKKPi0);
349 numchan.push_back(1);
350 numchan.push_back(1);
351 numchan.push_back(3);
352 decaylist=kaonList.minus()* kaonList.plus()* pi0List;
353 }
354 else if(channel=="D0toPi0Pi0"){
355 numchan.push_back( EvtRecDTag::kD0toPi0Pi0);
356 numchan.push_back(3);
357 numchan.push_back(3);
358 decaylist=pi0List* pi0List;
359 }
360 else if(channel=="D0toKsKs"){
361 numchan.push_back( EvtRecDTag::kD0toKsKs);
362 numchan.push_back(5);
363 numchan.push_back(5);
364 decaylist=ksList* ksList;
365 }
366 else if(channel=="D0toKsKsPi0"){
367 numchan.push_back( EvtRecDTag::kD0toKsKsPi0);
368 numchan.push_back(5);
369 numchan.push_back(5);
370 numchan.push_back(3);
371 decaylist=ksList* ksList* pi0List;
372 }
373 else if(channel=="D0toKsPi0Pi0"){
374 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0);
375 numchan.push_back(5);
376 numchan.push_back(3);
377 numchan.push_back(3);
378 decaylist=ksList* pi0List* pi0List;
379 }
380 else if(channel=="D0toKsKK"){
381 numchan.push_back( EvtRecDTag::kD0toKsKK);
382 numchan.push_back(5);
383 numchan.push_back(1);
384 numchan.push_back(1);
385 decaylist=ksList* kaonList.minus()* kaonList.plus();
386 }
387 else if(channel=="D0toKsEta"){
388 numchan.push_back( EvtRecDTag::kD0toKsEta);
389 numchan.push_back(5);
390 numchan.push_back(4);
391 decaylist=ksList* etaList;
392 }
393 else if(channel=="D0toPi0Pi0Pi0"){
394 numchan.push_back( EvtRecDTag::kD0toPi0Pi0Pi0);
395 numchan.push_back(3);
396 numchan.push_back(3);
397 numchan.push_back(3);
398 decaylist=pi0List* pi0List* pi0List;
399 }
400 else if(channel=="D0toKsKsKs"){
401 numchan.push_back( EvtRecDTag::kD0toKsKsKs);
402 numchan.push_back(5);
403 numchan.push_back(5);
404 numchan.push_back(5);
405 decaylist=ksList* ksList* ksList;
406 }
407 else if(channel=="D0toPiPiPiPi"){
408 numchan.push_back( EvtRecDTag::kD0toPiPiPiPi);
409 numchan.push_back(2);
410 numchan.push_back(2);
411 numchan.push_back(2);
412 numchan.push_back(2);
413 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
414 }
415 else if(channel=="D0toPiPiPi0Pi0"){
416 numchan.push_back( EvtRecDTag::kD0toPiPiPi0Pi0);
417 numchan.push_back(2);
418 numchan.push_back(2);
419 numchan.push_back(3);
420 numchan.push_back(3);
421 decaylist=pionList.plus()* pionList.minus()* pi0List* pi0List;
422 }
423 else if(channel=="D0toKKPiPi"){
424 numchan.push_back( EvtRecDTag::kD0toKKPiPi);
425 numchan.push_back(1);
426 numchan.push_back(1);
427 numchan.push_back(2);
428 numchan.push_back(2);
429 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus();
430 }
431 else if(channel=="D0toKKPi0Pi0"){
432 numchan.push_back( EvtRecDTag::kD0toKKPi0Pi0);
433 numchan.push_back(1);
434 numchan.push_back(1);
435 numchan.push_back(3);
436 numchan.push_back(3);
437 decaylist=kaonList.minus()* kaonList.plus()* pi0List* pi0List;
438 }
439 else if(channel=="D0toKsKsPiPi"){
440 numchan.push_back( EvtRecDTag::kD0toKsKsPiPi);
441 numchan.push_back(5);
442 numchan.push_back(5);
443 numchan.push_back(2);
444 numchan.push_back(2);
445 decaylist=ksList* ksList* pionList.plus()* pionList.minus();
446 }
447 else if(channel=="D0toPiPiPiPiPi0"){
448 numchan.push_back( EvtRecDTag::kD0toPiPiPiPiPi0);
449 numchan.push_back(2);
450 numchan.push_back(2);
451 numchan.push_back(2);
452 numchan.push_back(2);
453 numchan.push_back(3);
454 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
455 }
456 else if(channel=="D0toKsPiPiPiPi"){
457 numchan.push_back( EvtRecDTag::kD0toKsPiPiPiPi);
458 numchan.push_back(5);
459 numchan.push_back(2);
460 numchan.push_back(2);
461 numchan.push_back(2);
462 numchan.push_back(2);
463 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
464 }
465 else if(channel=="D0toKKPiPiPi0"){
466 numchan.push_back( EvtRecDTag::kD0toKKPiPiPi0);
467 numchan.push_back(1);
468 numchan.push_back(1);
469 numchan.push_back(2);
470 numchan.push_back(2);
471 numchan.push_back(3);
472 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
473 }
474 else if(channel=="D0toKsPi0Eta"){
475 numchan.push_back( EvtRecDTag::kD0toKsPi0Eta);
476 numchan.push_back(5);
477 numchan.push_back(3);
478 numchan.push_back(4);
479 decaylist=ksList* pi0List* etaList;
480 }
481 else if(channel=="D0toKsEPPiPiEta"){
482 numchan.push_back( EvtRecDTag::kD0toKsEPPiPiEta);
483 numchan.push_back(5);
484 numchan.push_back(6);
485 CDDecayList epList(eptoPiPiEtaSelector);
486 epList=pionList.plus()* pionList.minus()* etaList;
487 decaylist= ksList* epList;
488 }
489 else if(channel=="D0toKsEPRhoGam"){
490 numchan.push_back( EvtRecDTag::kD0toKsEPRhoGam);
491 numchan.push_back(5);
492 numchan.push_back(7);
493 CDDecayList rhoList(rhotoPiPiSelector);
494 rhoList=pionList.plus()* pionList.minus();
495 CDDecayList epList(eptoRhoGamSelector);
496 epList=rhoList* photonList;
497 decaylist= ksList* epList;
498 }
499
500
501 CDDecayList::iterator D_begin =decaylist.particle_begin();
502 CDDecayList::iterator D_end =decaylist.particle_end();
503
504 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
505
506 EvtRecDTag* recDTag = new EvtRecDTag;
507 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
508
509
510 vector<int> trackid, showerid;
511 vector<int> kaonid, pionid;
512 int charm=0;
513 int numofchildren=numchan.size()-1;
514
515 for(int i=0; i< numofchildren;i++){
516
517 const CDCandidate& daughter=(*it).particle().child(i);
518 if(isFlavorMode && i==0)
519 charm=-daughter.charge();
520
521 if(numchan[i+1]==1){
522 const EvtRecTrack* track=daughter.track();
523 trackid.push_back(track->trackId());
524 kaonid.push_back(track->trackId());
525 }
526 else if(numchan[i+1]==2){
527 const EvtRecTrack* track=daughter.track();
528 trackid.push_back(track->trackId());
529 pionid.push_back(track->trackId());
530 }
531 else if ( numchan[i+1]==3){
532 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
533 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
534 showerid.push_back(hiEnGamma->trackId());
535 showerid.push_back(loEnGamma->trackId());
536 }
537 else if ( numchan[i+1]==4){
538 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
539 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
540 showerid.push_back(hiEnGamma->trackId());
541 showerid.push_back(loEnGamma->trackId());
542 }
543 else if ( numchan[i+1]==5){
544 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
545 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
546 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
547 trackid.push_back(pion1Trk->trackId());
548 trackid.push_back(pion2Trk->trackId());
549 }
550 else if (numchan[i+1]==6){
551 const CDCandidate& apion = daughter.decay().child(0);
552 const CDCandidate& spion = daughter.decay().child(1);
553 const CDCandidate& eta = daughter.decay().child(2);
554 const EvtRecTrack* apiontrk = apion.track();
555 const EvtRecTrack* spiontrk = spion.track();
556 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
557 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
558
559 trackid.push_back(apiontrk->trackId());
560 trackid.push_back(spiontrk->trackId());
561 showerid.push_back(hiEnGamma->trackId());
562 showerid.push_back(loEnGamma->trackId());
563
564 }
565 else if (numchan[i+1]==7){
566 const CDCandidate& rho = daughter.decay().child(0);
567 const CDCandidate& gamma = daughter.decay().child(1);
568 const CDCandidate& apion = rho.decay().child(0);
569 const CDCandidate& spion = rho.decay().child(1);
570
571
572 const EvtRecTrack* apiontrk = apion.track();
573 const EvtRecTrack* spiontrk = spion.track();
574 const EvtRecTrack* gammatrk = gamma.photon();
575
576
577 trackid.push_back(apiontrk->trackId());
578 trackid.push_back(spiontrk->trackId());
579 showerid.push_back(gammatrk->trackId());
580
581 }
582
583
584 }//end of filling track and shower ids
585
586
587 saveD0Info(it, ebeam, charm, numofchildren, recDTag);
588 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
589 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
590
591
592 trackid.clear();
593 showerid.clear();
594 kaonid.clear();
595 pionid.clear();
596
597
598 //write recdtag out
599 recDTagCol->push_back(recDTag);
600
601 } //end of decay list iterator
602
603 numchan.clear();
604
605 } //end of reconstrucing all D0 decay lists
606
607
608 return StatusCode::SUCCESS;
609}
610
611void NeutralDReconstruction::saveD0Info(CDDecayList::iterator it, double ebeam, int charm, int numofchildren, EvtRecDTag* recDTag){
612
613 double mass = (*it).particle().mass();
614 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
615 recDTag->setp4(p4);
616
617 p4.boost(-m_beta);
618 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
619 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
620 double deltaE_CMS = p4.t() - ebeam;
621
622 int tag=(*it).particle().userTag();
623 if(tag==1)
624 recDTag->settype( EvtRecDTag::Tight );
625 else
626 recDTag->settype( EvtRecDTag::Loose );
627 recDTag->setcharge(0);
628 recDTag->setcharm(charm);
629 recDTag->setnumOfChildren(numofchildren);
630 recDTag->setmass(mass);
631 recDTag->setmBC(mbc_CMS);
632 recDTag->setbeamE(ebeam);
633 recDTag->setdeltaE(deltaE_CMS);
634
635}
636
637void NeutralDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
638 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
639
640 vector<EvtRecTrackIterator> trktemp;
641 vector<EvtRecTrackIterator> shrtemp;
642
643 //fill tracks
644 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
645
646 bool isothertrack=true;
647 for(int i=0; i<trackid.size(); i++){
648 if( (*trk)->trackId()==trackid[i] ){
649 trktemp.push_back(trk);
650 isothertrack=false;
651 break;
652 }
653 }
654 if(isothertrack)
655 recDTag->addOtherTrack(*trk);
656 }
657 for(int i=0; i<trackid.size();i++){
658 for(int j=0; j<trktemp.size(); j++){
659 EvtRecTrackIterator trk=trktemp[j];
660 if( (*trk)->trackId()==trackid[i])
661 recDTag->addTrack(*trktemp[j]);
662 }
663 }
664
665
666 //fill showers
667 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
668 bool isothershower=true;
669 for(int i=0; i<showerid.size(); i++){
670 if( (*shr)->trackId()==showerid[i] ){
671 shrtemp.push_back(shr);
672 isothershower=false;
673 break;
674 }
675 }
676 if(isothershower)
677 recDTag->addOtherShower(*shr);
678 }
679
680 for(int i=0; i<showerid.size();i++){
681 for(int j=0; j<shrtemp.size(); j++){
682 EvtRecTrackIterator shr=shrtemp[j];
683 if( (*shr)->trackId()==showerid[i])
684 recDTag->addShower(*shrtemp[j]);
685 }
686 }
687
688
689}
690
691
692void NeutralDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
693
694 bool iskaon=false,ispion=false;
695
696
697 // save track ids which passed pion/kaon cuts
698
699 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
700 recDTag->addKaonId( (*kit).particle().track() );
701 }
702
703 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
704 recDTag->addPionId( (*pit).particle().track() );
705 }
706
707
708
709 /*
710 for(int i=0; i<kaonid.size(); i++){
711 bool ithkaon=false;
712 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
713 if((*kit).particle().track()->trackId()==kaonid[i]){
714 ithkaon=true;
715 break;
716 }
717 }
718 if(!ithkaon) break;
719 if(i==kaonid.size()-1)
720 iskaon=true;
721 }
722
723 for(int i=0; i<pionid.size(); i++){
724 bool ithpion=false;
725 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
726 if((*pit).particle().track()->trackId()==pionid[i]){
727 ithpion=true;
728 break;
729 }
730 }
731 if(!ithpion) break;
732 if(i==pionid.size()-1)
733 ispion=true;
734 }
735
736
737 if( iskaon && ispion)
738 recDTag->settype( EvtRecDTag::Tight );
739 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
740 recDTag->settype( EvtRecDTag::Tight );
741 else if( kaonid.size()==0 && pionid.size()==0 )
742 recDTag->settype( EvtRecDTag::Tight );
743 */
744
745}
746
747
748vector<string> NeutralDReconstruction::getlist(string& filename ){
749
750 string channel;
751 vector<string> temp;
752
753 ifstream inFile;
754
755 inFile.open(filename.c_str());
756 if (!inFile) {
757 cout << "Unable to open decay list file";
758 exit(1); // terminate with error
759 }
760
761 while (inFile >> channel) {
762 temp.push_back(channel);
763 }
764
765 inFile.close();
766
767 return temp;
768
769}
mass
double mass
Definition: CosmicGenerator.cxx:138
dc_fill
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
Definition: DCFillableChargedList.h:152
EvtRecDTagCol
ObjectVector< EvtRecDTag > EvtRecDTagCol
Definition: EvtRecDTag.h:233
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
eptoPiPiEtaSelector
LocalEptoPiPiEtaSelector eptoPiPiEtaSelector
Definition: LocalEptoPiPiEtaSelector.cxx:39
eptoRhoGamSelector
LocalEptoRhoGamSelector eptoRhoGamSelector
Definition: LocalEptoRhoGamSelector.cxx:38
etatoGGSelector
LocalEtatoGGSelector etatoGGSelector
Definition: LocalEtatoGGSelector.cxx:37
kaonSelector
LocalKaonSelector kaonSelector
Definition: LocalKaonSelector.cxx:112
ksSelector
LocalKsSelector ksSelector
Definition: LocalKsSelector.cxx:40
photonSelector
LocalPhotonSelector photonSelector
Definition: LocalPhotonSelector.cxx:115
pi0Selector
LocalPi0Selector pi0Selector
Definition: LocalPi0Selector.cxx:39
pionSelector
LocalPionSelector pionSelector
Definition: LocalPionSelector.cxx:113
rhotoPiPiSelector
LocalRhotoPiPiSelector rhotoPiPiSelector
Definition: LocalRhotoPiPiSelector.cxx:36
neutralDSelector
NeutralDSelector neutralDSelector
Definition: NeutralDSelector.cxx:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
CDCandidate::decay
virtual const DecayEvidence & decay() const
Definition: CDCandidate.cxx:356
CDCandidate::photon
virtual const EvtRecTrack * photon() const
Definition: CDCandidate.cxx:256
CDCandidate::track
virtual const EvtRecTrack * track() const
Definition: CDCandidate.cxx:236
CDCandidate::charge
int charge() const
Definition: CDCandidate.cxx:154
CDCandidate::navKshort
virtual const EvtRecVeeVertex * navKshort() const
Definition: CDCandidate.cxx:316
CDCandidate::navPi0
virtual const EvtRecPi0 * navPi0() const
Definition: CDCandidate.cxx:276
CDCandidate::navEta
virtual const EvtRecEtaToGG * navEta() const
Definition: CDCandidate.cxx:296
CDDecay::child
const CDCandidate & child(unsigned int aPosition) const
Definition: CDDecay.cxx:247
EvtRecDTag::kD0toKsEPRhoGam
@ kD0toKsEPRhoGam
Definition: EvtRecDTag.h:88
EvtRecDTag::kD0toKKPiPiPi0
@ kD0toKKPiPiPi0
Definition: EvtRecDTag.h:84
EvtRecDTag::kD0toPiPiPiPiPi0
@ kD0toPiPiPiPiPi0
Definition: EvtRecDTag.h:82
EvtRecDTag::kD0toKsPiPiPiPi
@ kD0toKsPiPiPiPi
Definition: EvtRecDTag.h:83
EvtRecDTag::kD0toKPiPiPi
@ kD0toKPiPiPi
Definition: EvtRecDTag.h:52
EvtRecDTag::kD0toKKPi0Pi0
@ kD0toKKPi0Pi0
Definition: EvtRecDTag.h:80
EvtRecDTag::kD0toKsKsPi0
@ kD0toKsKsPi0
Definition: EvtRecDTag.h:70
EvtRecDTag::kD0toKsKsPiPi
@ kD0toKsKsPiPi
Definition: EvtRecDTag.h:81
EvtRecDTag::kD0toKsPiPiPi0
@ kD0toKsPiPiPi0
Definition: EvtRecDTag.h:61
EvtRecDTag::kD0toKsPi0Pi0
@ kD0toKsPi0Pi0
Definition: EvtRecDTag.h:71
EvtRecDTag::kD0toPiPiPiPi
@ kD0toPiPiPiPi
Definition: EvtRecDTag.h:77
EvtRecDTag::kD0toPiPiPi0Pi0
@ kD0toPiPiPi0Pi0
Definition: EvtRecDTag.h:78
EvtRecDTag::kD0toKsEPPiPiEta
@ kD0toKsEPPiPiEta
Definition: EvtRecDTag.h:87
EvtRecDTag::kD0toKsKPiPi0
@ kD0toKsKPiPi0
Definition: EvtRecDTag.h:58
EvtRecDTag::kD0toPi0Pi0Pi0
@ kD0toPi0Pi0Pi0
Definition: EvtRecDTag.h:74
EvtRecDTag::kD0toKPiPiPiPi0
@ kD0toKPiPiPiPi0
Definition: EvtRecDTag.h:54
EvtRecDTag::kD0toPiPiPi0
@ kD0toPiPiPi0
Definition: EvtRecDTag.h:64
EvtRecDTag::kD0toKPiPi0Pi0
@ kD0toKPiPi0Pi0
Definition: EvtRecDTag.h:51
EvtRecDTag::kD0toKsPi0Eta
@ kD0toKsPi0Eta
Definition: EvtRecDTag.h:85
EvtRecDTag::addOtherTrack
void addOtherTrack(const SmartRef< EvtRecTrack > track)
Definition: EvtRecDTag.h:186
EvtRecDTag::settype
void settype(SelType type)
Definition: EvtRecDTag.h:172
EvtRecDTag::setdecayMode
void setdecayMode(DecayMode decayMode)
Definition: EvtRecDTag.h:171
EvtRecDTag::setp4
void setp4(HepLorentzVector p4)
Definition: EvtRecDTag.h:180
EvtRecDTag::setmass
void setmass(double mass)
Definition: EvtRecDTag.h:174
EvtRecDTag::addOtherShower
void addOtherShower(const SmartRef< EvtRecTrack > shower)
Definition: EvtRecDTag.h:188
EvtRecDTag::addKaonId
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
Definition: EvtRecDTag.h:192
EvtRecDTag::setdeltaE
void setdeltaE(double deltaE)
Definition: EvtRecDTag.h:176
EvtRecDTag::setcharm
void setcharm(int charm)
Definition: EvtRecDTag.h:178
EvtRecDTag::setmBC
void setmBC(double mBC)
Definition: EvtRecDTag.h:175
EvtRecDTag::setcharge
void setcharge(int charge)
Definition: EvtRecDTag.h:177
EvtRecDTag::addPionId
void addPionId(const SmartRef< EvtRecTrack > pionId)
Definition: EvtRecDTag.h:190
EvtRecDTag::setbeamE
void setbeamE(double beamE)
Definition: EvtRecDTag.h:173
EvtRecDTag::addShower
void addShower(const SmartRef< EvtRecTrack > shower)
Definition: EvtRecDTag.h:184
EvtRecDTag::setnumOfChildren
void setnumOfChildren(int numOfChildren)
Definition: EvtRecDTag.h:179
EvtRecDTag::addTrack
void addTrack(const SmartRef< EvtRecTrack > track)
Definition: EvtRecDTag.h:182
EvtRecEtaToGG::hiEnGamma
const EvtRecTrack * hiEnGamma() const
Definition: EvtRecEtaToGG.h:30
EvtRecEtaToGG::loEnGamma
const EvtRecTrack * loEnGamma() const
Definition: EvtRecEtaToGG.h:31
EvtRecPi0::loEnGamma
const EvtRecTrack * loEnGamma() const
Definition: EvtRecPi0.h:31
EvtRecPi0::hiEnGamma
const EvtRecTrack * hiEnGamma() const
Definition: EvtRecPi0.h:30
EvtRecTrack::trackId
int trackId() const
Definition: EvtRecTrack.h:32
EvtRecVeeVertex::daughter
SmartRef< EvtRecTrack > & daughter(int i)
Definition: EvtRecVeeVertex.h:89
LocalKaonSelector::setpidtype
void setpidtype(int type)
Definition: LocalKaonSelector.h:15
LocalPionSelector::setpidtype
void setpidtype(int type)
Definition: LocalPionSelector.h:14
NeutralDReconstruction::NeutralDReconstruction
NeutralDReconstruction(const std::string &name, ISvcLocator *pSvcLocator)
Definition: NeutralDReconstruction.cxx:55
NeutralDReconstruction::savetrack
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
Definition: NeutralDReconstruction.cxx:637
NeutralDReconstruction::getlist
vector< string > getlist(string &filename)
Definition: NeutralDReconstruction.cxx:748
NeutralDReconstruction::pidtag
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
Definition: NeutralDReconstruction.cxx:692
NeutralDReconstruction::saveD0Info
void saveD0Info(CDDecayList::iterator, double, int, int, EvtRecDTag *)
Definition: NeutralDReconstruction.cxx:611
NeutralDSelector::setbeta
void setbeta(Hep3Vector beta)
Definition: NeutralDSelector.h:14
NeutralDSelector::setebeam
void setebeam(double ebeam)
Definition: NeutralDSelector.h:13
dchain::CandidateList::size
virtual int size() const
Definition: CandidateList.cc:143
dchain::ChargedVisibleList::plus
ChosenChargeList< Charged, CandidateClass > & plus() const
Definition: ChargedVisibleList.cc:176
dchain::ChargedVisibleList::minus
ChosenChargeList< Charged, CandidateClass > & minus() const
Definition: ChargedVisibleList.cc:186
dchain::DecayList::particle_begin
virtual iterator particle_begin()
Definition: DecayList.cc:169
dchain::DecayList::particle_end
virtual iterator particle_end()
Definition: DecayList.cc:176
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:111
EventModel::EvtRec::EvtRecDTagCol
_EXTERN_ std::string EvtRecDTagCol
Definition: EventModel.h:117
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:112
Event
Definition: Event.h:21