BOSS 7.1.2
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
13#include "EventModel/Event.h"
15
22
26#include "BesDChain/CDPi0List.h"
27#include "BesDChain/CDEtaList.h"
28#include "BesDChain/CDKsList.h"
30
31#include "McTruth/McParticle.h"
33#include "VertexFit/VertexFit.h"
35
36
48
49#include "DTagAlg/utility.h"
50
51#include <fstream>
52
53using namespace Event;
54
55DECLARE_COMPONENT(NeutralDReconstruction)
56//*******************************************************************************************
57NeutralDReconstruction::NeutralDReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
58Algorithm(name, pSvcLocator) {
59 //Declare the properties
60 declareProperty( "debug", m_debug = false );
61 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
62 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
63 declareProperty( "UseVertexfit", m_usevertexfit = false );
64 declareProperty( "BeamE", m_beamE = 1.8865 );
65 declareProperty( "D0List", m_decaylist = "test.txt" );
66 declareProperty("UseVFRefine", m_useVFrefine = true);
67 declareProperty( "UseBFieldCorr", m_useBFC = true);
68}
69
70//******************************************************************************************
72 MsgStream log(msgSvc(), name());
73 log << MSG::INFO << "in initialize()" <<endreq;
74
75 m_irun=-100;
76 m_beta.setX(0.011);
77 m_beta.setY(0);
78 m_beta.setZ(0);
79
80 chanlist=getlist(m_decaylist);
81
82 return StatusCode::SUCCESS;
83}
84
85//********************************************************************************************
87 MsgStream log(msgSvc(), name());
88 log << MSG::INFO << "in finalize()" << endreq;
89
90 chanlist.clear();
91
92 return StatusCode::SUCCESS;
93}
94
95StatusCode NeutralDReconstruction::registerEvtRecDTagCol(
96 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
97 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
98 aNewEvtRecDTagCol);
99 if (sc != StatusCode::SUCCESS) {
100 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
101 }
102 return sc;
103}
104
105//*********************************************************************************************
107 MsgStream log(msgSvc(), name());
108 log << MSG::INFO << "in execute()" << endreq;
109
110 StatusCode sc;
111
112 //////////////////
113 // Read REC data
114 /////////////////
115 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
116 int event= eventHeader->eventNumber();
117 // if ( m_debug || ( (event & 0x3FF) == 0 ) )
118 //std::cout << "event: " << event << std::endl;
119
120 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
121 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
122 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
123 << " " << eventHeader->eventNumber() << endreq;
124 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
125 << recEvent->totalCharged() << " , "
126 << recEvent->totalNeutral() << " , "
127 << recEvent->totalTracks() <<endreq;
128
129 EvtRecTrackIterator charged_begin = recTrackCol->begin();
130 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
131
132 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
133 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
134
135
136 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
137 if ( ! recPi0Col ) {
138 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
139 return StatusCode::FAILURE;
140 }
141
142
143
144 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
145 if ( ! recEtaToGGCol ) {
146 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
147 return StatusCode::FAILURE;
148 }
149
150
151 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
152 if ( ! recVeeVertexCol ) {
153 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
154 return StatusCode::FAILURE;
155 }
156
157
158 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
159 if (!recDTagCol) {
160 log << MSG::FATAL << "EvtRecDTagCol is not registered in TDS" << endreq;
161 return StatusCode::FAILURE;
162 }
163
164 //get primary vertex from db
165 Hep3Vector xorigin(0,0,0);
166 IVertexDbSvc* vtxsvc;
167 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
168 if (vtxsvc->isVertexValid()) {
169
170 //vertex[0] = vx; vertex[1]= vy; vertex[2] = vz;
171 double* vertex = vtxsvc->PrimaryVertex();
172 xorigin.setX(vertex[0]);
173 xorigin.setY(vertex[1]);
174 xorigin.setZ(vertex[2]);
175 }
176 utility util;
177
178
179
180 //registered in the parent file DTag.cxx
181 /*
182 if (!recDTagCol) {
183 recDTagCol = new EvtRecDTagCol;
184 sc = registerEvtRecDTagCol(recDTagCol, log);
185 if (sc != StatusCode::SUCCESS) {
186 return sc;
187 }
188 }
189 */
190
191 /////////////////////////////
192 //reconstruct particle lists
193 /////////////////////////////
194
197 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
198 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
199 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
200
201 CDKsList ksList(ksSelector);
202 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
203
204 // do a secondary vertex fit and cut on the results
205 map<EvtRecVeeVertex*, vector< double > > fitinfo;
206 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
207 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
208
209 if(m_useVFrefine){
210 fitinfo[ks] = util.SecondaryVFitref(ks, vtxsvc);
211 }else{
212 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
213 }
214
215 }
216
217 CDPi0List pi0List(pi0Selector);
218 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
219
220 CDEtaList etaList(etatoGGSelector);
221 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
222
223
224 //pion/kaon list with PID
227 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
228 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
229
230 int run = eventHeader->runNumber();
231 m_ievt = eventHeader->eventNumber();
232 m_nChrg = recEvent->totalCharged();
233 m_nNeu = recEvent->totalNeutral();
234 m_nPion = pionList.size();
235 m_nKaon = kaonList.size();
236 m_nPi0 = pi0List.size();
237 m_nKs = ksList.size();
238
239 ///////////////////////
240 // get beam energy and beta
241 ///////////////////////
242
243 if(m_ReadBeamEFromDB && m_irun!=run){
244 m_irun=run;
245 if(m_usecalibBeamE)
246 m_readDb.setcalib(true);
247 m_beamE=m_readDb.getbeamE(m_irun, m_beamE);
248 if(run>0)
249 m_beta=m_readDb.getbeta();
250
251 //cout<<"use beam E from data base:"<<m_beamE<<endl;
252 }
253 double ebeam=m_beamE;
254
255 //////////////////////////////
256 //reconstruct decay lists
257 /////////////////////////////
258
259
260
261 for(int list=0;list<chanlist.size();list++){
262
263 string channel=chanlist[list];
264 vector<int> numchan;
267 CDDecayList decaylist(neutralDSelector);
268
269 bool isFlavorMode=false;
270
271 //K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5, Eta'(pipiEta):6, Eta'(rhogam):7
272 //the fist element of the vector stands for decay mode,
273 //the rest will be particles, and size of the vector minus 1 will be number of daughers.
274 //the first particle in the vector will be used to get charm
275
276 if(channel=="D0toKPi") {
277 numchan.push_back( EvtRecDTag::kD0toKPi );
278 numchan.push_back(1);
279 numchan.push_back(2);
280 decaylist=kaonList.minus() * pionList.plus();
281 isFlavorMode=true;
282 }
283 else if(channel=="D0toKPiPi0"){
284 numchan.push_back( EvtRecDTag::kD0toKPiPi0);
285 numchan.push_back(1);
286 numchan.push_back(2);
287 numchan.push_back(3);
288 decaylist=kaonList.minus() * pionList.plus() * pi0List;
289 isFlavorMode=true;
290 }
291 else if(channel=="D0toKPiPi0Pi0"){
292 numchan.push_back( EvtRecDTag::kD0toKPiPi0Pi0);
293 numchan.push_back(1);
294 numchan.push_back(2);
295 numchan.push_back(3);
296 numchan.push_back(3);
297 decaylist=kaonList.minus() * pionList.plus() * pi0List * pi0List;
298 isFlavorMode=true;
299 }
300 else if(channel=="D0toKPiPiPi"){
301 numchan.push_back( EvtRecDTag::kD0toKPiPiPi);
302 numchan.push_back(1);
303 numchan.push_back(2);
304 numchan.push_back(2);
305 numchan.push_back(2);
306 decaylist=kaonList.minus()* pionList.plus()* pionList.plus() * pionList.minus();
307 isFlavorMode=true;
308 }
309 else if(channel=="D0toKPiPiPiPi0"){
310 numchan.push_back( EvtRecDTag::kD0toKPiPiPiPi0);
311 numchan.push_back(1);
312 numchan.push_back(2);
313 numchan.push_back(2);
314 numchan.push_back(2);
315 numchan.push_back(3);
316 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
317 isFlavorMode=true;
318 }
319 else if(channel=="D0toKPiEta"){
320 numchan.push_back( EvtRecDTag::kD0toKPiEta);
321 numchan.push_back(1);
322 numchan.push_back(2);
323 numchan.push_back(4);
324 decaylist=kaonList.minus() * pionList.plus() * etaList;
325 isFlavorMode=true;
326 }
327 else if(channel=="D0toKPiPi0Pi0Pi0"){
328 numchan.push_back( EvtRecDTag::kD0toKPiPi0Pi0Pi0);
329 numchan.push_back(1);
330 numchan.push_back(2);
331 numchan.push_back(3);
332 numchan.push_back(3);
333 numchan.push_back(3);
334 decaylist=kaonList.minus() * pionList.plus() * pi0List * pi0List * pi0List;
335 isFlavorMode=true;
336 }
337 else if(channel=="D0toKPiPi0Eta"){
338 numchan.push_back( EvtRecDTag::kD0toKPiPi0Eta);
339 numchan.push_back(1);
340 numchan.push_back(2);
341 numchan.push_back(3);
342 numchan.push_back(4);
343 decaylist=kaonList.minus() * pionList.plus() * pi0List * etaList;
344 isFlavorMode=true;
345 }
346 else if(channel=="D0toKPiEPPiPiEta"){
347 numchan.push_back( EvtRecDTag::kD0toKPiEPPiPiEta);
348 numchan.push_back(1);
349 numchan.push_back(2);
350 numchan.push_back(6);
352 epList=pionList.plus()* pionList.minus()* etaList;
353 decaylist=kaonList.minus() * pionList.plus() * epList;
354 isFlavorMode=true;
355 }
356 else if(channel=="D0toKPiEPRhoGam"){
357 numchan.push_back( EvtRecDTag::kD0toKPiEPRhoGam);
358 numchan.push_back(1);
359 numchan.push_back(2);
360 numchan.push_back(7);
362 rhoList=pionList.plus()* pionList.minus();
364 epList=rhoList* photonList;
365 decaylist=kaonList.minus() * pionList.plus() * epList;
366 isFlavorMode=true;
367 }
368 else if(channel=="D0toKKKPi"){
369 numchan.push_back( EvtRecDTag::kD0toKKKPi);
370 numchan.push_back(1);
371 numchan.push_back(1);
372 numchan.push_back(1);
373 numchan.push_back(2);
374 decaylist=kaonList.minus()* kaonList.plus()* kaonList.minus() * pionList.plus();
375 isFlavorMode=true;
376 }
377 else if(channel=="D0toKsKPi"){
378 numchan.push_back( EvtRecDTag::kD0toKsKPi);
379 numchan.push_back(5);
380 numchan.push_back(1);
381 numchan.push_back(2);
382 decaylist=ksList* kaonList.minus()* pionList.plus();
383 }
384 else if(channel=="D0toKsKPiPi0"){
385 numchan.push_back( EvtRecDTag::kD0toKsKPiPi0);
386 numchan.push_back(5);
387 numchan.push_back(1);
388 numchan.push_back(2);
389 numchan.push_back(3);
390 decaylist=ksList* kaonList.minus() * pionList.plus()* pi0List;
391 }
392 else if(channel=="D0toKsPiPi"){
393 numchan.push_back( EvtRecDTag::kD0toKsPiPi);
394 numchan.push_back(5);
395 numchan.push_back(2);
396 numchan.push_back(2);
397 decaylist=ksList* pionList.plus()* pionList.minus();
398 }
399 else if(channel=="D0toKsPiPiPi0"){
400 numchan.push_back( EvtRecDTag::kD0toKsPiPiPi0);
401 numchan.push_back(5);
402 numchan.push_back(2);
403 numchan.push_back(2);
404 numchan.push_back(3);
405 decaylist=ksList* pionList.plus()* pionList.minus()* pi0List;
406 }
407 else if(channel=="D0toKsPi0"){
408 numchan.push_back( EvtRecDTag::kD0toKsPi0);
409 numchan.push_back(5);
410 numchan.push_back(3);
411 decaylist=ksList* pi0List;
412 }
413 else if(channel=="D0toPiPiPi0"){
414 numchan.push_back( EvtRecDTag::kD0toPiPiPi0);
415 numchan.push_back(2);
416 numchan.push_back(2);
417 numchan.push_back(3);
418 decaylist=pionList.plus()* pionList.minus()* pi0List;
419 }
420 else if(channel=="D0toPiPi"){
421 numchan.push_back( EvtRecDTag::kD0toPiPi);
422 numchan.push_back(2);
423 numchan.push_back(2);
424 decaylist=pionList.plus()* pionList.minus();
425 }
426 else if(channel=="D0toKK"){
427 numchan.push_back( EvtRecDTag::kD0toKK);
428 numchan.push_back(1);
429 numchan.push_back(1);
430 decaylist=kaonList.minus()* kaonList.plus();
431 }
432 else if(channel=="D0toKKPi0"){
433 numchan.push_back( EvtRecDTag::kD0toKKPi0);
434 numchan.push_back(1);
435 numchan.push_back(1);
436 numchan.push_back(3);
437 decaylist=kaonList.minus()* kaonList.plus()* pi0List;
438 }
439 else if(channel=="D0toPi0Pi0"){
440 numchan.push_back( EvtRecDTag::kD0toPi0Pi0);
441 numchan.push_back(3);
442 numchan.push_back(3);
443 decaylist=pi0List* pi0List;
444 }
445 else if(channel=="D0toKsKs"){
446 numchan.push_back( EvtRecDTag::kD0toKsKs);
447 numchan.push_back(5);
448 numchan.push_back(5);
449 decaylist=ksList* ksList;
450 }
451 else if(channel=="D0toKsKsPi0"){
452 numchan.push_back( EvtRecDTag::kD0toKsKsPi0);
453 numchan.push_back(5);
454 numchan.push_back(5);
455 numchan.push_back(3);
456 decaylist=ksList* ksList* pi0List;
457 }
458 else if(channel=="D0toKsPi0Pi0"){
459 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0);
460 numchan.push_back(5);
461 numchan.push_back(3);
462 numchan.push_back(3);
463 decaylist=ksList* pi0List* pi0List;
464 }
465 else if(channel=="D0toKsKK"){
466 numchan.push_back( EvtRecDTag::kD0toKsKK);
467 numchan.push_back(5);
468 numchan.push_back(1);
469 numchan.push_back(1);
470 decaylist=ksList* kaonList.minus()* kaonList.plus();
471 }
472 else if(channel=="D0toKsEta"){
473 numchan.push_back( EvtRecDTag::kD0toKsEta);
474 numchan.push_back(5);
475 numchan.push_back(4);
476 decaylist=ksList* etaList;
477 }
478 else if(channel=="D0toPi0Pi0Pi0"){
479 numchan.push_back( EvtRecDTag::kD0toPi0Pi0Pi0);
480 numchan.push_back(3);
481 numchan.push_back(3);
482 numchan.push_back(3);
483 decaylist=pi0List* pi0List* pi0List;
484 }
485 else if(channel=="D0toKsKsKs"){
486 numchan.push_back( EvtRecDTag::kD0toKsKsKs);
487 numchan.push_back(5);
488 numchan.push_back(5);
489 numchan.push_back(5);
490 decaylist=ksList* ksList* ksList;
491 }
492 else if(channel=="D0toPiPiPiPi"){
493 numchan.push_back( EvtRecDTag::kD0toPiPiPiPi);
494 numchan.push_back(2);
495 numchan.push_back(2);
496 numchan.push_back(2);
497 numchan.push_back(2);
498 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
499 }
500 else if(channel=="D0toPiPiPi0Pi0"){
501 numchan.push_back( EvtRecDTag::kD0toPiPiPi0Pi0);
502 numchan.push_back(2);
503 numchan.push_back(2);
504 numchan.push_back(3);
505 numchan.push_back(3);
506 decaylist=pionList.plus()* pionList.minus()* pi0List* pi0List;
507 }
508 else if(channel=="D0toKKPiPi"){
509 numchan.push_back( EvtRecDTag::kD0toKKPiPi);
510 numchan.push_back(1);
511 numchan.push_back(1);
512 numchan.push_back(2);
513 numchan.push_back(2);
514 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus();
515 }
516 else if(channel=="D0toKKPi0Pi0"){
517 numchan.push_back( EvtRecDTag::kD0toKKPi0Pi0);
518 numchan.push_back(1);
519 numchan.push_back(1);
520 numchan.push_back(3);
521 numchan.push_back(3);
522 decaylist=kaonList.minus()* kaonList.plus()* pi0List* pi0List;
523 }
524 else if(channel=="D0toKsKsPiPi"){
525 numchan.push_back( EvtRecDTag::kD0toKsKsPiPi);
526 numchan.push_back(5);
527 numchan.push_back(5);
528 numchan.push_back(2);
529 numchan.push_back(2);
530 decaylist=ksList* ksList* pionList.plus()* pionList.minus();
531 }
532 else if(channel=="D0toPiPiPiPiPi0"){
533 numchan.push_back( EvtRecDTag::kD0toPiPiPiPiPi0);
534 numchan.push_back(2);
535 numchan.push_back(2);
536 numchan.push_back(2);
537 numchan.push_back(2);
538 numchan.push_back(3);
539 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
540 }
541 else if(channel=="D0toKsPiPiPiPi"){
542 numchan.push_back( EvtRecDTag::kD0toKsPiPiPiPi);
543 numchan.push_back(5);
544 numchan.push_back(2);
545 numchan.push_back(2);
546 numchan.push_back(2);
547 numchan.push_back(2);
548 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
549 }
550 else if(channel=="D0toKKPiPiPi0"){
551 numchan.push_back( EvtRecDTag::kD0toKKPiPiPi0);
552 numchan.push_back(1);
553 numchan.push_back(1);
554 numchan.push_back(2);
555 numchan.push_back(2);
556 numchan.push_back(3);
557 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
558 }
559 else if(channel=="D0toKsPi0Eta"){
560 numchan.push_back( EvtRecDTag::kD0toKsPi0Eta);
561 numchan.push_back(5);
562 numchan.push_back(3);
563 numchan.push_back(4);
564 decaylist=ksList* pi0List* etaList;
565 }
566 else if(channel=="D0toKsEPPiPiEta"){
567 numchan.push_back( EvtRecDTag::kD0toKsEPPiPiEta);
568 numchan.push_back(5);
569 numchan.push_back(6);
571 epList=pionList.plus()* pionList.minus()* etaList;
572 decaylist= ksList* epList;
573 }
574 else if(channel=="D0toKsEPRhoGam"){
575 numchan.push_back( EvtRecDTag::kD0toKsEPRhoGam);
576 numchan.push_back(5);
577 numchan.push_back(7);
579 rhoList=pionList.plus()* pionList.minus();
581 epList=rhoList* photonList;
582 decaylist= ksList* epList;
583 }
584 else if(channel=="D0toKsPi0Pi0Pi0"){
585 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0Pi0);
586 numchan.push_back(5);
587 numchan.push_back(3);
588 numchan.push_back(3);
589 numchan.push_back(3);
590 decaylist=ksList* pi0List* pi0List* pi0List;
591 }
592 else if(channel=="D0toKsPiPiPi0Pi0"){
593 numchan.push_back( EvtRecDTag::kD0toKsPiPiPi0Pi0);
594 numchan.push_back(5);
595 numchan.push_back(2);
596 numchan.push_back(2);
597 numchan.push_back(3);
598 numchan.push_back(3);
599 decaylist=ksList* pionList.plus()* pionList.minus()* pi0List* pi0List;
600 }
601 else if(channel=="D0toKsPiPiEta"){
602 numchan.push_back( EvtRecDTag::kD0toKsPiPiEta);
603 numchan.push_back(5);
604 numchan.push_back(2);
605 numchan.push_back(2);
606 numchan.push_back(4);
607 decaylist=ksList* pionList.plus()* pionList.minus()* etaList;
608 }
609 else if(channel=="D0toKsPi0EPPiPiEta"){
610 numchan.push_back( EvtRecDTag::kD0toKsPi0EPPiPiEta);
611 numchan.push_back(5);
612 numchan.push_back(3);
613 numchan.push_back(6);
615 epList=pionList.plus()* pionList.minus()* etaList;
616 decaylist= ksList* pi0List* epList;
617 }
618 else if(channel=="D0toKsPi0EPRhoGam"){
619 numchan.push_back( EvtRecDTag::kD0toKsPi0EPRhoGam);
620 numchan.push_back(5);
621 numchan.push_back(3);
622 numchan.push_back(7);
624 rhoList=pionList.plus()* pionList.minus();
626 epList=rhoList* photonList;
627 decaylist= ksList* pi0List* epList;
628 }
629 else if(channel=="D0toPiPiEta"){
630 numchan.push_back( EvtRecDTag::kD0toPiPiEta);
631 numchan.push_back(2);
632 numchan.push_back(2);
633 numchan.push_back(4);
634 decaylist=pionList.plus()* pionList.minus()* etaList;
635 }
636 else if(channel=="D0toPiPiPi0Eta"){
637 numchan.push_back( EvtRecDTag::kD0toPiPiPi0Eta);
638 numchan.push_back(2);
639 numchan.push_back(2);
640 numchan.push_back(3);
641 numchan.push_back(4);
642 decaylist=pionList.plus()* pionList.minus()* pi0List* etaList;
643 }
644 else if(channel=="D0toPiPiEPPiPiEta"){
645 numchan.push_back( EvtRecDTag::kD0toPiPiEPPiPiEta);
646 numchan.push_back(2);
647 numchan.push_back(2);
648 numchan.push_back(6);
650 epList=pionList.plus()* pionList.minus()* etaList;
651 decaylist=pionList.plus()* pionList.minus()* epList;
652 }
653 else if(channel=="D0toPiPiEPRhoGam"){
654 numchan.push_back( EvtRecDTag::kD0toPiPiEPRhoGam);
655 numchan.push_back(2);
656 numchan.push_back(2);
657 numchan.push_back(7);
659 rhoList=pionList.plus()* pionList.minus();
661 epList=rhoList* photonList;
662 decaylist=pionList.plus()* pionList.minus()* epList;
663 }
664 else if(channel=="D0toKKEta"){
665 numchan.push_back( EvtRecDTag::kD0toKKEta);
666 numchan.push_back(1);
667 numchan.push_back(1);
668 numchan.push_back(4);
669 decaylist=kaonList.minus()* kaonList.plus()* etaList;
670 }
671
672
673 CDDecayList::iterator D_begin =decaylist.particle_begin();
674 CDDecayList::iterator D_end =decaylist.particle_end();
675
676 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
677
678 EvtRecDTag* recDTag = new EvtRecDTag;
679 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
680
681
682 vector<int> trackid, showerid;
683 vector<int> kaonid, pionid;
684 int charm=0;
685 int numofchildren=numchan.size()-1;
686
687 for(int i=0; i< numofchildren;i++){
688
689 const CDCandidate& daughter=(*it).particle().child(i);
690 if(isFlavorMode && i==0)
691 charm=-daughter.charge();
692
693 if(numchan[i+1]==1){
694 const EvtRecTrack* track=daughter.track();
695 trackid.push_back(track->trackId());
696 kaonid.push_back(track->trackId());
697 }
698 else if(numchan[i+1]==2){
699 const EvtRecTrack* track=daughter.track();
700 trackid.push_back(track->trackId());
701 pionid.push_back(track->trackId());
702 }
703 else if ( numchan[i+1]==3){
704 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
705 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
706 showerid.push_back(hiEnGamma->trackId());
707 showerid.push_back(loEnGamma->trackId());
708 }
709 else if ( numchan[i+1]==4){
710 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
711 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
712 showerid.push_back(hiEnGamma->trackId());
713 showerid.push_back(loEnGamma->trackId());
714 }
715 else if ( numchan[i+1]==5){
716 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
717 // fill fit info
718 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
719 // fill tracks
720 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
721 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
722 trackid.push_back(pion1Trk->trackId());
723 trackid.push_back(pion2Trk->trackId());
724 }
725 else if (numchan[i+1]==6){
726 const CDCandidate& apion = daughter.decay().child(0);
727 const CDCandidate& spion = daughter.decay().child(1);
728 const CDCandidate& eta = daughter.decay().child(2);
729 const EvtRecTrack* apiontrk = apion.track();
730 const EvtRecTrack* spiontrk = spion.track();
731 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
732 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
733
734 trackid.push_back(apiontrk->trackId());
735 trackid.push_back(spiontrk->trackId());
736 showerid.push_back(hiEnGamma->trackId());
737 showerid.push_back(loEnGamma->trackId());
738
739 }
740 else if (numchan[i+1]==7){
741 const CDCandidate& rho = daughter.decay().child(0);
742 const CDCandidate& gamma = daughter.decay().child(1);
743 const CDCandidate& apion = rho.decay().child(0);
744 const CDCandidate& spion = rho.decay().child(1);
745
746
747 const EvtRecTrack* apiontrk = apion.track();
748 const EvtRecTrack* spiontrk = spion.track();
749 const EvtRecTrack* gammatrk = gamma.photon();
750
751
752 trackid.push_back(apiontrk->trackId());
753 trackid.push_back(spiontrk->trackId());
754 showerid.push_back(gammatrk->trackId());
755
756 }
757
758
759 }//end of filling track and shower ids
760
761
762 saveD0Info(it, ebeam, charm, numofchildren, recDTag);
763 if(m_useBFC){
764 //After use BFC package, we need to update information of Ks and Lambda to calculate D.
765 updateKsInfo(it, ebeam, charm, numofchildren, recDTag, numchan, vtxsvc, m_useVFrefine);
766 }
767 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
768 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
769
770 if(m_usevertexfit){
771 if(m_debug) cout<<"beforevfit:"<<endl;
772
773 HepLorentzVector p4change_vfit;
774
775 if(m_useVFrefine){
776 p4change_vfit=util.vfitref(channel, kaonid, pionid, xorigin, charged_begin);
777 }else{
778 p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
779 }
780
781 recDTag->setp4(recDTag->p4()+p4change_vfit);
782 }
783
784
785 trackid.clear();
786 showerid.clear();
787 kaonid.clear();
788 pionid.clear();
789
790
791 //write recdtag out
792 recDTagCol->push_back(recDTag);
793
794 } //end of decay list iterator
795
796 numchan.clear();
797
798 } //end of reconstrucing all D0 decay lists
799
800
801 return StatusCode::SUCCESS;
802}
803
804void NeutralDReconstruction::saveD0Info(CDDecayList::iterator it, double ebeam, int charm, int numofchildren, EvtRecDTag* recDTag){
805
806 double mass = (*it).particle().mass();
807 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
808 recDTag->setp4(p4);
809
810 p4.boost(-m_beta);
811 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
812 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
813 double deltaE_CMS = p4.t() - ebeam;
814
815 int tag=(*it).particle().userTag();
816 if(tag==1)
817 recDTag->settype( EvtRecDTag::Tight );
818 else
819 recDTag->settype( EvtRecDTag::Loose );
820 recDTag->setcharge(0);
821 recDTag->setcharm(charm);
822 recDTag->setnumOfChildren(numofchildren);
823 recDTag->setmass(mass);
824 recDTag->setmBC(mbc_CMS);
825 recDTag->setbeamE(ebeam);
826 recDTag->setdeltaE(deltaE_CMS);
827
828}
829
830void NeutralDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
831 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
832
833 vector<EvtRecTrackIterator> trktemp;
834 vector<EvtRecTrackIterator> shrtemp;
835
836 //fill tracks
837 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
838
839 bool isothertrack=true;
840 for(int i=0; i<trackid.size(); i++){
841 if( (*trk)->trackId()==trackid[i] ){
842 trktemp.push_back(trk);
843 isothertrack=false;
844 break;
845 }
846 }
847 if(isothertrack)
848 recDTag->addOtherTrack(*trk);
849 }
850 for(int i=0; i<trackid.size();i++){
851 for(int j=0; j<trktemp.size(); j++){
852 EvtRecTrackIterator trk=trktemp[j];
853 if( (*trk)->trackId()==trackid[i])
854 recDTag->addTrack(*trktemp[j]);
855 }
856 }
857
858
859 //fill showers
860 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
861 bool isothershower=true;
862 for(int i=0; i<showerid.size(); i++){
863 if( (*shr)->trackId()==showerid[i] ){
864 shrtemp.push_back(shr);
865 isothershower=false;
866 break;
867 }
868 }
869 if(isothershower)
870 recDTag->addOtherShower(*shr);
871 }
872
873 for(int i=0; i<showerid.size();i++){
874 for(int j=0; j<shrtemp.size(); j++){
875 EvtRecTrackIterator shr=shrtemp[j];
876 if( (*shr)->trackId()==showerid[i])
877 recDTag->addShower(*shrtemp[j]);
878 }
879 }
880
881
882}
883
884
885void NeutralDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
886
887 bool iskaon=false,ispion=false;
888
889
890 // save track ids which passed pion/kaon cuts
891
892 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
893 recDTag->addKaonId( (*kit).particle().track() );
894 }
895
896 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
897 recDTag->addPionId( (*pit).particle().track() );
898 }
899
900
901
902 /*
903 for(int i=0; i<kaonid.size(); i++){
904 bool ithkaon=false;
905 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
906 if((*kit).particle().track()->trackId()==kaonid[i]){
907 ithkaon=true;
908 break;
909 }
910 }
911 if(!ithkaon) break;
912 if(i==kaonid.size()-1)
913 iskaon=true;
914 }
915
916 for(int i=0; i<pionid.size(); i++){
917 bool ithpion=false;
918 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
919 if((*pit).particle().track()->trackId()==pionid[i]){
920 ithpion=true;
921 break;
922 }
923 }
924 if(!ithpion) break;
925 if(i==pionid.size()-1)
926 ispion=true;
927 }
928
929
930 if( iskaon && ispion)
931 recDTag->settype( EvtRecDTag::Tight );
932 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
933 recDTag->settype( EvtRecDTag::Tight );
934 else if( kaonid.size()==0 && pionid.size()==0 )
935 recDTag->settype( EvtRecDTag::Tight );
936 */
937
938}
939
940
941
942vector<string> NeutralDReconstruction::getlist(string& filename ){
943
944 string channel;
945 vector<string> temp;
946
947 ifstream inFile;
948
949 inFile.open(filename.c_str());
950 if (!inFile) {
951 cout << "Unable to open decay list file";
952 exit(1); // terminate with error
953 }
954
955 while (inFile >> channel) {
956 temp.push_back(channel);
957 }
958
959 inFile.close();
960
961 return temp;
962
963}
964
965
966void NeutralDReconstruction::updateKsInfo(CDDecayList::iterator it, double ebeam, int charm, int numofchildren, EvtRecDTag* recDTag, vector<int> numchan, IVertexDbSvc* vtxsvc, bool m_useVFrefine){
967
968//This function is used to update the information of Ks and Lambda. Because DTagAlg package uses the list from VeeVertexFit package which input the information of Ks and Lambda while reconstructing the dst, the information of them don't use the MBF correction and we need to update them here.
969
970 //If we don't use BF Correctoin, please close this function.
971
972 HepLorentzVector p4D0((*it).particle().momentum(), (*it).particle().energy());
973 HepLorentzVector p4D0_New;
974
975 for(int i = 0; i < numofchildren; i++){
976 const CDCandidate& daughter=(*it).particle().child(i);
977
978 if(numchan[i+1]==5){
979 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>(daughter.navKshort());
980 HepVector veewks = aKsCand->w(); //px, py, pz, E, x, y, z
981 HepLorentzVector p4KsCand;
982 p4KsCand.setX(veewks[0]);
983 p4KsCand.setY(veewks[1]);
984 p4KsCand.setZ(veewks[2]);
985 p4KsCand.setE(veewks[3]);
986 HepLorentzVector p4wks = p4KsCand;
987
988 utility util;
989
990 //by default: px, py, pz, E, chisq, ifok
991 vector<double> vp4ks_new = util.UpdatedKsIfo(aKsCand,vtxsvc,m_useVFrefine);
992 HepLorentzVector p4wks_new;
993 p4wks_new.setX(vp4ks_new[0]);
994 p4wks_new.setY(vp4ks_new[1]);
995 p4wks_new.setZ(vp4ks_new[2]);
996 p4wks_new.setE(vp4ks_new[3]);
997
998 p4D0 = p4D0 - p4wks + p4wks_new;
999
1000 }
1001
1002 }
1003
1004 p4D0_New = p4D0;
1005
1006 double mass = p4D0_New.m();
1007 recDTag->setp4(p4D0_New);
1008 p4D0_New.boost(-m_beta);
1009 double mbc2_CMS = ebeam*ebeam - p4D0_New.v().mag2();
1010 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
1011 double deltaE_CMS = p4D0_New.t() - ebeam;
1012
1013 int tag=(*it).particle().userTag();
1014 if(tag==1)
1015 recDTag->settype( EvtRecDTag::Tight );
1016 else
1017 recDTag->settype( EvtRecDTag::Loose );
1018 recDTag->setcharge(0);
1019 recDTag->setcharm(charm);
1020 recDTag->setnumOfChildren(numofchildren);
1021 recDTag->setmass(mass);
1022 recDTag->setmBC(mbc_CMS);
1023 recDTag->setbeamE(ebeam);
1024 recDTag->setdeltaE(deltaE_CMS);
1025
1026
1027}
1028
1029
1030
1031
1032
1033
double mass
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
double ebeam
ObjectVector< EvtRecDTag > EvtRecDTagCol
Definition EvtRecDTag.h:329
EvtRecTrackCol::iterator EvtRecTrackIterator
LocalEptoPiPiEtaSelector eptoPiPiEtaSelector
LocalEptoRhoGamSelector eptoRhoGamSelector
LocalEtatoGGSelector etatoGGSelector
LocalKaonSelector kaonSelector
LocalKsSelector ksSelector
LocalPhotonSelector photonSelector
LocalPi0Selector pi0Selector
LocalPionSelector pionSelector
LocalRhotoPiPiSelector rhotoPiPiSelector
NeutralDSelector neutralDSelector
IMessageSvc * msgSvc()
virtual const DecayEvidence & decay() const
virtual const EvtRecTrack * photon() const
virtual const EvtRecTrack * track() const
int charge() const
virtual const EvtRecVeeVertex * navKshort() const
virtual const EvtRecPi0 * navPi0() const
virtual const EvtRecEtaToGG * navEta() const
const CDCandidate & child(unsigned int aPosition) const
Definition CDDecay.cxx:247
@ kD0toKsEPRhoGam
Definition EvtRecDTag.h:97
@ kD0toKsPiPiPi0Pi0
Definition EvtRecDTag.h:100
@ kD0toPiPiPiPiPi0
Definition EvtRecDTag.h:91
@ kD0toKsPiPiPiPi
Definition EvtRecDTag.h:92
@ kD0toKPiEPRhoGam
Definition EvtRecDTag.h:63
@ kD0toPiPiEPPiPiEta
Definition EvtRecDTag.h:106
@ kD0toPiPiEPRhoGam
Definition EvtRecDTag.h:107
@ kD0toKsPi0EPRhoGam
Definition EvtRecDTag.h:103
@ kD0toKPiEPPiPiEta
Definition EvtRecDTag.h:62
@ kD0toKsPi0Pi0Pi0
Definition EvtRecDTag.h:98
@ kD0toKsPi0EPPiPiEta
Definition EvtRecDTag.h:102
@ kD0toPiPiPi0Pi0
Definition EvtRecDTag.h:87
@ kD0toKsEPPiPiEta
Definition EvtRecDTag.h:96
@ kD0toKPiPi0Pi0Pi0
Definition EvtRecDTag.h:60
@ kD0toKPiPiPiPi0
Definition EvtRecDTag.h:57
void addOtherTrack(const SmartRef< EvtRecTrack > track)
Definition EvtRecDTag.h:275
void settype(SelType type)
Definition EvtRecDTag.h:254
void setdecayMode(DecayMode decayMode)
Definition EvtRecDTag.h:253
HepLorentzVector p4() const
Definition EvtRecDTag.h:237
void setp4(HepLorentzVector p4)
Definition EvtRecDTag.h:262
void setmass(double mass)
Definition EvtRecDTag.h:256
void addToFitInfo(double ksmass, double chi2, double length, double error)
Definition EvtRecDTag.h:264
void addOtherShower(const SmartRef< EvtRecTrack > shower)
Definition EvtRecDTag.h:277
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
Definition EvtRecDTag.h:281
void setdeltaE(double deltaE)
Definition EvtRecDTag.h:258
void setcharm(int charm)
Definition EvtRecDTag.h:260
void setmBC(double mBC)
Definition EvtRecDTag.h:257
void setcharge(int charge)
Definition EvtRecDTag.h:259
void addPionId(const SmartRef< EvtRecTrack > pionId)
Definition EvtRecDTag.h:279
void setbeamE(double beamE)
Definition EvtRecDTag.h:255
void addShower(const SmartRef< EvtRecTrack > shower)
Definition EvtRecDTag.h:273
void setnumOfChildren(int numOfChildren)
Definition EvtRecDTag.h:261
void addTrack(const SmartRef< EvtRecTrack > track)
Definition EvtRecDTag.h:271
const EvtRecTrack * hiEnGamma() const
const EvtRecTrack * loEnGamma() const
const EvtRecTrack * loEnGamma() const
Definition EvtRecPi0.h:31
const EvtRecTrack * hiEnGamma() const
Definition EvtRecPi0.h:30
int trackId() const
Definition EvtRecTrack.h:32
SmartRef< EvtRecTrack > & daughter(int i)
double mass() const
const HepVector & w() const
virtual bool isVertexValid()=0
virtual double * PrimaryVertex()=0
void setpidtype(int type)
void setpidtype(int type)
void updateKsInfo(CDDecayList::iterator, double, int, int, EvtRecDTag *, vector< int >, IVertexDbSvc *, bool)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
vector< string > getlist(string &filename)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void saveD0Info(CDDecayList::iterator, double, int, int, EvtRecDTag *)
void setbeta(Hep3Vector beta)
void setebeam(double ebeam)
virtual int size() const
ChosenChargeList< Charged, CandidateClass > & plus() const
ChosenChargeList< Charged, CandidateClass > & minus() const
virtual iterator particle_begin()
Definition DecayList.cc:169
virtual iterator particle_end()
Definition DecayList.cc:176
HepLorentzVector vfitref(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:59
vector< double > SecondaryVFitref(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:199
vector< double > UpdatedKsIfo(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc, bool m_useVFrefine)
Definition utility.cxx:752
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:406
vector< double > SecondaryVFit(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:540
_EXTERN_ std::string EvtRecEvent
Definition EventModel.h:116
_EXTERN_ std::string EvtRecDTagCol
Definition EventModel.h:122
_EXTERN_ std::string EvtRecTrackCol
Definition EventModel.h:117
Definition Event.h:21
double double double * p4
Definition qcdloop1.h:77