BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
DsReconstruction.cxx
Go to the documentation of this file.
1//
2// All Ds decay modes Reconstruction
3//
4#include <fstream>
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
49#include "DTagAlg/DsSelector.h"
50
51#include "DTagAlg/utility.h"
52
53using namespace Event;
54DECLARE_COMPONENT(DsReconstruction)
55//*******************************************************************************************
56DsReconstruction::DsReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
57 Algorithm(name, pSvcLocator) {
58 //Declare the properties
59 declareProperty( "debug", m_debug = false );
60 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
61 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
62 declareProperty( "UseVertexfit", m_usevertexfit = false );
63 declareProperty( "BeamE", m_beamE=2.015 );
64 declareProperty( "DsList", m_decaylist = "test.txt" );
65 declareProperty("UseVFRefine", m_useVFrefine = true);
66 declareProperty( "UseBFieldCorr", m_useBFC = true);
67}
68
69//******************************************************************************************
71 MsgStream log(msgSvc(), name());
72 log << MSG::INFO << "in initialize()" <<endreq;
73
74 m_irun=-100;
75 m_beta.setX(0.011);
76 m_beta.setY(0);
77 m_beta.setZ(0);
78 chanlist=getlist(m_decaylist);
79
80 return StatusCode::SUCCESS;
81}
82
83//********************************************************************************************
85 MsgStream log(msgSvc(), name());
86 log << MSG::INFO << "in finalize()" << endreq;
87
88 chanlist.clear();
89
90 return StatusCode::SUCCESS;
91}
92
93StatusCode DsReconstruction::registerEvtRecDTagCol(
94 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
95 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
96 aNewEvtRecDTagCol);
97 if (sc != StatusCode::SUCCESS) {
98 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
99 }
100 return sc;
101}
102
103//*********************************************************************************************
105 MsgStream log(msgSvc(), name());
106 log << MSG::INFO << "in execute()" << endreq;
107
108 StatusCode sc;
109
110 //////////////////
111 // Read REC data
112 /////////////////
113 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
114 int event= eventHeader->eventNumber();
115 // if ( m_debug || ( (event & 0x3FF) == 0 ) )
116 //std::cout << "event: " << event << std::endl;
117
118 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
119 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
120 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
121 << " " << eventHeader->eventNumber() << endreq;
122 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
123 << recEvent->totalCharged() << " , "
124 << recEvent->totalNeutral() << " , "
125 << recEvent->totalTracks() <<endreq;
126
127 EvtRecTrackIterator charged_begin = recTrackCol->begin();
128 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
129
130 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
131 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
132
133
134 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
135 if ( ! recPi0Col ) {
136 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
137 return StatusCode::FAILURE;
138 }
139
140 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
141 if ( ! recEtaToGGCol ) {
142 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
143 return StatusCode::FAILURE;
144 }
145
146
147 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
148 if ( ! recVeeVertexCol ) {
149 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
150 return StatusCode::FAILURE;
151 }
152
153
154 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
155 if (!recDTagCol) {
156 log << MSG::FATAL << "EvtRecDTagCol is not registered yet" << endreq;
157 return StatusCode::FAILURE;
158 }
159
160
161 //get primary vertex from db
162 Hep3Vector xorigin(0,0,0);
163 IVertexDbSvc* vtxsvc;
164 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
165 if (vtxsvc->isVertexValid()) {
166
167 //vertex[0] = vx; vertex[1]= vy; vertex[2] = vz;
168 double* vertex = vtxsvc->PrimaryVertex();
169 xorigin.setX(vertex[0]);
170 xorigin.setY(vertex[1]);
171 xorigin.setZ(vertex[2]);
172 }
173 utility util;
174
175
176 //registered in DTag.cxx
177 /*
178 if (!recDTagCol) {
179 recDTagCol = new EvtRecDTagCol;
180 sc = registerEvtRecDTagCol(recDTagCol, log);
181 if (sc != StatusCode::SUCCESS) {
182 return sc;
183 }
184 }
185 */
186
187
188 /////////////////////////////
189 //reconstruct particle lists
190 /////////////////////////////
191
194 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
195 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
196 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
197
198 CDKsList ksList(ksSelector);
199 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
200
201
202 // do a secondary vertex fit and cut on the results
203 map<EvtRecVeeVertex*, vector< double > > fitinfo;
204 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
205 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
206
207 if(m_useVFrefine){
208 fitinfo[ks] = util.SecondaryVFitref(ks, vtxsvc);
209 }else{
210 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
211 }
212
213 }
214
215 CDPi0List pi0List(pi0Selector);
216 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
217
218 CDEtaList etaList(etatoGGSelector);
219 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
220
221 //pion/kaon list with PID
224 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
225 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
226
227
228 int run = eventHeader->runNumber();
229 m_ievt = eventHeader->eventNumber();
230 m_nChrg = recEvent->totalCharged();
231 m_nNeu = recEvent->totalNeutral();
232 m_nPion = pionList.size();
233 m_nKaon = kaonList.size();
234 m_nPi0 = pi0List.size();
235 m_nKs = ksList.size();
236
237
238 ///////////////////////
239 // get beam energy and beta
240 ///////////////////////
241
242 if(m_ReadBeamEFromDB && m_irun!=run){
243 m_irun=run;
244 if(m_usecalibBeamE)
245 m_readDb.setcalib(true);
246 m_beamE=m_readDb.getbeamE(m_irun,m_beamE);
247 if(run>0)
248 m_beta=m_readDb.getbeta();
249 //cout<<"use beam E from data base:"<<m_beamE<<endl;
250 }
251 double ebeam=m_beamE;
252
253 //////////////////////////////
254 //reconstruct decay lists
255 /////////////////////////////
256
257
258 for(int list=0;list<chanlist.size();list++){
259
260 string channel=chanlist[list];
261 vector<int> numchan;
263 dsSelector.setbeta(m_beta);
264 CDDecayList decaylist(dsSelector);
265
266 //K+/-: 1, Pi+/-:2, Pi0:3,
267 //Eta: 4, Ks:5,
268 //eta'(pipieta): 6,
269 //eta'(rhogamma): 7
270 //eta'(pipieta(pipipi0)): 8
271 //eta(pipipi0): 9,
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
275 if(channel=="DstoKsK") {
276 numchan.push_back( EvtRecDTag::kDstoKsK );
277 numchan.push_back(5);
278 numchan.push_back(1);
279 decaylist=ksList* kaonList.plus();
280 }
281 else if(channel=="DstoKKPi") {
282 numchan.push_back( EvtRecDTag::kDstoKKPi );
283 numchan.push_back(1);
284 numchan.push_back(1);
285 numchan.push_back(2);
286 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus();
287 }
288 else if(channel=="DstoKsKPi0") {
289 numchan.push_back( EvtRecDTag::kDstoKsKPi0 );
290 numchan.push_back(5);
291 numchan.push_back(1);
292 numchan.push_back(3);
293 decaylist=ksList* kaonList.plus()* pi0List;
294 }
295 else if(channel=="DstoKsKsPi") {
296 numchan.push_back( EvtRecDTag::kDstoKsKsPi );
297 numchan.push_back(5);
298 numchan.push_back(5);
299 numchan.push_back(2);
300 decaylist=ksList* ksList* pionList.plus();
301 }
302 else if(channel=="DstoKKPiPi0") {
303 numchan.push_back( EvtRecDTag::kDstoKKPiPi0 );
304 numchan.push_back(1);
305 numchan.push_back(1);
306 numchan.push_back(2);
307 numchan.push_back(3);
308 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pi0List;
309 }
310 else if(channel=="DstoKsKplusPiPi") {
311 numchan.push_back( EvtRecDTag::kDstoKsKplusPiPi );
312 numchan.push_back(5);
313 numchan.push_back(1);
314 numchan.push_back(2);
315 numchan.push_back(2);
316 decaylist=ksList* kaonList.plus()* pionList.plus()* pionList.minus();
317 }
318 else if(channel=="DstoKsKminusPiPi") {
319 numchan.push_back( EvtRecDTag::kDstoKsKminusPiPi );
320 numchan.push_back(5);
321 numchan.push_back(1);
322 numchan.push_back(2);
323 numchan.push_back(2);
324 decaylist=ksList* kaonList.minus()* pionList.plus()* pionList.plus();
325 }
326 else if(channel=="DstoKKPiPiPi") {
327 numchan.push_back( EvtRecDTag::kDstoKKPiPiPi );
328 numchan.push_back(1);
329 numchan.push_back(1);
330 numchan.push_back(2);
331 numchan.push_back(2);
332 numchan.push_back(2);
333 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.plus()* pionList.minus();
334 }
335 else if(channel=="DstoPiPi0") {
336 numchan.push_back( EvtRecDTag::kDstoPiPi0 );
337 numchan.push_back(2);
338 numchan.push_back(3);
339 decaylist=pionList.plus()* pi0List;
340 }
341 else if(channel=="DstoPiPiPi") {
342 numchan.push_back( EvtRecDTag::kDstoPiPiPi );
343 numchan.push_back(2);
344 numchan.push_back(2);
345 numchan.push_back(2);
346 decaylist=pionList.plus()* pionList.plus()* pionList.minus();
347 }
348 else if(channel=="DstoPiPiPiPi0") {
349 numchan.push_back( EvtRecDTag::kDstoPiPiPiPi0 );
350 numchan.push_back(2);
351 numchan.push_back(2);
352 numchan.push_back(2);
353 numchan.push_back(3);
354 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
355 }
356 else if(channel=="DstoPiPiPiPiPi") {
357 numchan.push_back( EvtRecDTag::kDstoPiPiPiPiPi );
358 numchan.push_back(2);
359 numchan.push_back(2);
360 numchan.push_back(2);
361 numchan.push_back(2);
362 numchan.push_back(2);
363 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
364 }
365 else if(channel=="DstoPiPiPiPiPiPi0") {
366 numchan.push_back( EvtRecDTag::kDstoPiPiPiPiPiPi0 );
367 numchan.push_back(2);
368 numchan.push_back(2);
369 numchan.push_back(2);
370 numchan.push_back(2);
371 numchan.push_back(2);
372 numchan.push_back(3);
373 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
374 }
375 else if(channel=="DstoPiPiPiPi0Pi0") { // <------------- NEW MODE
376 numchan.push_back( EvtRecDTag::kDstoPiPiPiPi0Pi0 );
377 numchan.push_back(2);
378 numchan.push_back(2);
379 numchan.push_back(2);
380 numchan.push_back(3);
381 numchan.push_back(3);
382 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List* pi0List;
383 }
384 else if(channel=="DstoPiEta") {
385 numchan.push_back( EvtRecDTag::kDstoPiEta );
386 numchan.push_back(2);
387 numchan.push_back(4);
388 decaylist=pionList.plus()* etaList;
389 }
390 else if(channel=="DstoPiEtaPiPiPi0") { // <------------- NEW MODE
391 numchan.push_back( EvtRecDTag::kDstoPiEtaPiPiPi0 );
392 numchan.push_back(2);
393 numchan.push_back(9);
395 EtaList=pionList.plus()* pionList.minus()* pi0List;
396 decaylist=pionList.plus()* EtaList;
397 }
398 else if(channel=="DstoPiPi0Eta") {
399 numchan.push_back( EvtRecDTag::kDstoPiPi0Eta );
400 numchan.push_back(2);
401 numchan.push_back(3);
402 numchan.push_back(4);
403 decaylist=pionList.plus()* pi0List* etaList;
404 }
405 else if(channel=="DstoPiPi0EtaPiPiPi0") { // <---------- NEW MODE
406 numchan.push_back( EvtRecDTag::kDstoPiPi0EtaPiPiPi0 );
407 numchan.push_back(2);
408 numchan.push_back(3);
409 numchan.push_back(9);
411 EtaList=pionList.plus()* pionList.minus()* pi0List;
412 decaylist=pionList.plus()* pi0List* EtaList;
413 }
414 else if(channel=="DstoPiPiPiEta") {
415 numchan.push_back( EvtRecDTag::kDstoPiPiPiEta );
416 numchan.push_back(2);
417 numchan.push_back(2);
418 numchan.push_back(2);
419 numchan.push_back(4);
420 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* etaList;
421 }
422 else if(channel=="DstoPiPiPiEtaPiPiPi0") { // <---------- NEW MODE
423 numchan.push_back( EvtRecDTag::kDstoPiPiPiEtaPiPiPi0 );
424 numchan.push_back(2);
425 numchan.push_back(2);
426 numchan.push_back(2);
427 numchan.push_back(9);
429 EtaList=pionList.plus()* pionList.minus()* pi0List;
430 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* EtaList;
431 }
432 else if(channel=="DstoPiEPPiPiEta") {
433 numchan.push_back( EvtRecDTag::kDstoPiEPPiPiEta );
434 numchan.push_back(2);
435 numchan.push_back(6);
437 epList=pionList.plus()* pionList.minus()* etaList;
438 decaylist=pionList.plus()* epList;
439 }
440 else if(channel=="DstoPiPi0EPPiPiEta") {
441 numchan.push_back( EvtRecDTag::kDstoPiPi0EPPiPiEta );
442 numchan.push_back(2);
443 numchan.push_back(3);
444 numchan.push_back(6);
446 epList=pionList.plus()* pionList.minus()* etaList;
447 decaylist=pionList.plus()* pi0List* epList;
448 }
449 else if(channel=="DstoPiEPPiPiEtaPiPiPi0") {// <------------- New mode: 470 eta'(pipieta, eta->pipipi0)
450 numchan.push_back( EvtRecDTag::kDstoPiEPPiPiEtaPiPiPi0 );
451 numchan.push_back(2);
452 numchan.push_back(8);
454 EtaList=pionList.plus()* pionList.minus()* pi0List;
456 EtapList=pionList.plus()* pionList.minus()* EtaList;
457 decaylist=pionList.plus()* EtapList;
458 }
459 else if(channel=="DstoPiPi0EPPiPiEtaPiPiPi0") { //<------------ New mode: 471 eta'(pipieta, eta->pipipi0)
460 numchan.push_back( EvtRecDTag::kDstoPiPi0EPPiPiEtaPiPiPi0 );
461 numchan.push_back(2);
462 numchan.push_back(3);
463 numchan.push_back(8);
465 EtaList=pionList.plus()* pionList.minus()* pi0List;
467 EtapList=pionList.plus()* pionList.minus()* EtaList;
468 decaylist=pionList.plus()* pi0List* EtapList;
469 }
470
471 else if(channel=="DstoPiEPRhoGam") {
472 numchan.push_back( EvtRecDTag::kDstoPiEPRhoGam );
473 numchan.push_back(2);
474 numchan.push_back(7);
476 rhoList=pionList.plus()* pionList.minus();
478 epList=rhoList* photonList;
479 decaylist=pionList.plus()* epList;
480 }
481 else if(channel=="DstoPiPi0EPRhoGam") {
482 numchan.push_back( EvtRecDTag::kDstoPiPi0EPRhoGam );
483 numchan.push_back(2);
484 numchan.push_back(3);
485 numchan.push_back(7);
487 rhoList=pionList.plus()* pionList.minus();
489 epList=rhoList* photonList;
490 decaylist=pionList.plus()* pi0List* epList;
491 }
492 else if(channel=="DstoKsPi") {
493 numchan.push_back( EvtRecDTag::kDstoKsPi );
494 numchan.push_back(5);
495 numchan.push_back(2);
496 decaylist=ksList* pionList.plus();
497 }
498 else if(channel=="DstoKsPiPi0") {
499 numchan.push_back( EvtRecDTag::kDstoKsPiPi0 );
500 numchan.push_back(5);
501 numchan.push_back(2);
502 numchan.push_back(3);
503 decaylist=ksList* pionList.plus()* pi0List;
504 }
505 else if(channel=="DstoKPiPi") {
506 numchan.push_back( EvtRecDTag::kDstoKPiPi );
507 numchan.push_back(1);
508 numchan.push_back(2);
509 numchan.push_back(2);
510 decaylist=kaonList.plus()* pionList.plus()* pionList.minus();
511 }
512 else if(channel=="DstoKPiPiPi0") {
513 numchan.push_back( EvtRecDTag::kDstoKPiPiPi0 );
514 numchan.push_back(1);
515 numchan.push_back(2);
516 numchan.push_back(2);
517 numchan.push_back(3);
518 decaylist=kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
519 }
520 else if(channel=="DstoKKK") {
521 numchan.push_back( EvtRecDTag::kDstoKKK );
522 numchan.push_back(1);
523 numchan.push_back(1);
524 numchan.push_back(1);
525 decaylist=kaonList.minus()* kaonList.plus()* kaonList.plus();
526 }
527
528 CDDecayList::iterator D_begin =decaylist.particle_begin();
529 CDDecayList::iterator D_end =decaylist.particle_end();
530
531 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
532
533 EvtRecDTag* recDTag = new EvtRecDTag;
534 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
535
536 vector<int> trackid, showerid;
537 vector<int> kaonid, pionid;
538 int numofchildren=numchan.size()-1;
539
540 for(int i=0; i< numofchildren;i++){
541
542 const CDCandidate& daughter=(*it).particle().child(i);
543
544 if(numchan[i+1]==1){
545 const EvtRecTrack* track=daughter.track();
546 trackid.push_back(track->trackId());
547 kaonid.push_back(track->trackId());
548 }
549 else if(numchan[i+1]==2){
550 const EvtRecTrack* track=daughter.track();
551 trackid.push_back(track->trackId());
552 pionid.push_back(track->trackId());
553 }
554 else if ( numchan[i+1]==3){
555 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
556 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
557 showerid.push_back(hiEnGamma->trackId());
558 showerid.push_back(loEnGamma->trackId());
559 }
560 else if ( numchan[i+1]==4){
561 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
562 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
563 showerid.push_back(hiEnGamma->trackId());
564 showerid.push_back(loEnGamma->trackId());
565 }
566 else if ( numchan[i+1]==5){
567 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
568
569 // fill fit info
570 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
571 // fill tracks
572 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
573 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
574 trackid.push_back(pion1Trk->trackId());
575 trackid.push_back(pion2Trk->trackId());
576 }
577 else if (numchan[i+1]==6){
578 const CDCandidate& apion = daughter.decay().child(0);
579 const CDCandidate& spion = daughter.decay().child(1);
580 const CDCandidate& eta = daughter.decay().child(2);
581 const EvtRecTrack* apiontrk = apion.track();
582 const EvtRecTrack* spiontrk = spion.track();
583 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
584 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
585
586 trackid.push_back(apiontrk->trackId());
587 trackid.push_back(spiontrk->trackId());
588 showerid.push_back(hiEnGamma->trackId());
589 showerid.push_back(loEnGamma->trackId());
590
591 }
592 else if (numchan[i+1]==7){
593 const CDCandidate& rho = daughter.decay().child(0);
594 const CDCandidate& apion = rho.decay().child(0);
595 const CDCandidate& spion = rho.decay().child(1);
596 const CDCandidate& gamma = daughter.decay().child(1);
597
598 const EvtRecTrack* apiontrk = apion.track();
599 const EvtRecTrack* spiontrk = spion.track();
600 const EvtRecTrack* gammatrk = gamma.photon();
601
602
603 trackid.push_back(apiontrk->trackId());
604 trackid.push_back(spiontrk->trackId());
605 showerid.push_back(gammatrk->trackId());
606
607 }
608 else if (numchan[i+1]==8){ // <------------ NEW PART FOR ETAP TO PIPIETA TO PIPIPI0
609 const CDCandidate& apion = daughter.decay().child(0);
610 const CDCandidate& spion = daughter.decay().child(1);
611 const CDCandidate& eta = daughter.decay().child(2);
612 const CDCandidate& e0pion = eta.decay().child(0);
613 const CDCandidate& e1pion = eta.decay().child(1);
614 const CDCandidate& pi0 = eta.decay().child(2);
615
616 const EvtRecTrack* apiontrk = apion.track();
617 const EvtRecTrack* spiontrk = spion.track();
618 const EvtRecTrack* e0piontrk = e0pion.track();
619 const EvtRecTrack* e1piontrk = e1pion.track();
620
621 const EvtRecTrack* hiEnGamma=pi0.navPi0()->hiEnGamma();
622 const EvtRecTrack* loEnGamma=pi0.navPi0()->loEnGamma();
623
624 trackid.push_back(apiontrk->trackId());
625 trackid.push_back(spiontrk->trackId());
626 trackid.push_back(e0piontrk->trackId());
627 trackid.push_back(e1piontrk->trackId());
628 showerid.push_back(hiEnGamma->trackId());
629 showerid.push_back(loEnGamma->trackId());
630
631 }
632 else if (numchan[i+1]==9){ // <------- NEW PART FOR ETA TO PIPIPI0
633 const CDCandidate& apion = daughter.decay().child(0);
634 const CDCandidate& spion = daughter.decay().child(1);
635 const CDCandidate& pi0 = daughter.decay().child(2);
636 const EvtRecTrack* apiontrk = apion.track();
637 const EvtRecTrack* spiontrk = spion.track();
638 const EvtRecTrack* hiEnGamma=pi0.navPi0()->hiEnGamma();
639 const EvtRecTrack* loEnGamma=pi0.navPi0()->loEnGamma();
640
641 trackid.push_back(apiontrk->trackId());
642 trackid.push_back(spiontrk->trackId());
643 showerid.push_back(hiEnGamma->trackId());
644 showerid.push_back(loEnGamma->trackId());
645
646 }
647
648
649 }//end of filling track and shower ids
650
651
652 saveDsInfo(it, ebeam, numofchildren, recDTag);
653 if(m_useBFC){
654 //After use BFC package, we need to update information of Ks and Lambda to calculate D.
655 updateKsInfo(it, ebeam, numofchildren, recDTag, numchan, vtxsvc, m_useVFrefine);
656 }
657 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
658 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
659
660
661 if(m_usevertexfit){
662 if(m_debug) cout<<"beforevfit:"<<endl;
663
664 HepLorentzVector p4change_vfit;
665
666 if(m_useVFrefine){
667 p4change_vfit=util.vfitref(channel, kaonid, pionid, xorigin, charged_begin);
668 }else{
669 p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
670 }
671
672 recDTag->setp4(recDTag->p4()+p4change_vfit);
673 }
674
675
676 trackid.clear();
677 showerid.clear();
678 kaonid.clear();
679 pionid.clear();
680
681
682
683 //write dtag object out
684 recDTagCol->push_back(recDTag);
685
686 }//end of decaylist iterator
687
688 numchan.clear();
689
690 }//end of reconstrucing all D+ decay lists
691
692
693
694 return StatusCode::SUCCESS;
695}
696
697
698void DsReconstruction::saveDsInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag){
699
700 double mass = (*it).particle().mass();
701 int charge= (*it).particle().charge();
702 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
703 recDTag->setp4(p4);
704
705 p4.boost(-m_beta);
706 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
707 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
708 double deltaE_CMS = p4.t() - ebeam;
709
710 if((*it).particle().userTag()==1)
711 recDTag->settype( EvtRecDTag::Tight );
712 else
713 recDTag->settype( EvtRecDTag::Loose );
714 recDTag->setcharge(charge);
715 recDTag->setcharm(charge);
716 recDTag->setnumOfChildren(numofchildren);
717 recDTag->setmass(mass);
718 recDTag->setmBC(mbc_CMS);
719 recDTag->setbeamE(ebeam);
720 recDTag->setdeltaE(deltaE_CMS);
721
722}
723
724void DsReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
725 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
726
727 vector<EvtRecTrackIterator> trktemp;
728 vector<EvtRecTrackIterator> shrtemp;
729
730 //fill tracks
731 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
732
733 bool isothertrack=true;
734 for(int i=0; i<trackid.size(); i++){
735 if( (*trk)->trackId()==trackid[i] ){
736 trktemp.push_back(trk);
737 isothertrack=false;
738 break;
739 }
740 }
741 if(isothertrack)
742 recDTag->addOtherTrack(*trk);
743 }
744 for(int i=0; i<trackid.size();i++){
745 for(int j=0; j<trktemp.size(); j++){
746 EvtRecTrackIterator trk=trktemp[j];
747 if( (*trk)->trackId()==trackid[i])
748 recDTag->addTrack(*trktemp[j]);
749 }
750 }
751
752
753 //fill showers
754 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
755 bool isothershower=true;
756 for(int i=0; i<showerid.size(); i++){
757 if( (*shr)->trackId()==showerid[i] ){
758 shrtemp.push_back(shr);
759 isothershower=false;
760 break;
761 }
762 }
763 if(isothershower)
764 recDTag->addOtherShower(*shr);
765 }
766
767 for(int i=0; i<showerid.size();i++){
768 for(int j=0; j<shrtemp.size(); j++){
769 EvtRecTrackIterator shr=shrtemp[j];
770 if( (*shr)->trackId()==showerid[i])
771 recDTag->addShower(*shrtemp[j]);
772 }
773 }
774
775
776}
777
778void DsReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
779
780 bool iskaon=false,ispion=false;
781
782
783 // save track ids which passed pion/kaon cuts
784
785 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
786 recDTag->addKaonId( (*kit).particle().track() );
787 }
788
789 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
790 recDTag->addPionId( (*pit).particle().track() );
791 }
792
793
794 /*
795 for(int i=0; i<kaonid.size(); i++){
796 bool ithkaon=false;
797 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
798 if((*kit).particle().track()->trackId()==kaonid[i]){
799 ithkaon=true;
800 break;
801 }
802 }
803 if(!ithkaon) break;
804 if(i==kaonid.size()-1)
805 iskaon=true;
806 }
807
808 for(int i=0; i<pionid.size(); i++){
809 bool ithpion=false;
810 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
811 if((*pit).particle().track()->trackId()==pionid[i]){
812 ithpion=true;
813 break;
814 }
815 }
816 if(!ithpion) break;
817 if(i==pionid.size()-1)
818 ispion=true;
819 }
820
821
822 if( iskaon && ispion)
823 recDTag->settype( EvtRecDTag::Tight );
824 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
825 recDTag->settype( EvtRecDTag::Tight );
826 else if( kaonid.size()==0 && pionid.size()==0 )
827 recDTag->settype( EvtRecDTag::Tight );
828 */
829
830}
831
832
833vector<string> DsReconstruction::getlist(string& filename ){
834
835 string channel;
836 vector<string> temp;
837
838 ifstream inFile;
839
840 inFile.open(filename.c_str());
841 if (!inFile) {
842 cout << "Unable to open decay list file";
843 exit(1); // terminate with error
844 }
845
846 while (inFile >> channel) {
847 temp.push_back(channel);
848 }
849
850 inFile.close();
851
852 return temp;
853
854}
855
856
857void DsReconstruction::updateKsInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag, vector<int> numchan, IVertexDbSvc* vtxsvc, bool m_useVFrefine){
858
859//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.
860
861 //If we don't use BF Correctoin, please close this function.
862
863 HepLorentzVector p4Ds((*it).particle().momentum(), (*it).particle().energy());
864 HepLorentzVector p4Ds_New;
865
866 for(int i = 0; i < numofchildren; i++){
867 const CDCandidate& daughter=(*it).particle().child(i);
868
869 if(numchan[i+1]==5){
870 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>(daughter.navKshort());
871 HepVector veewks = aKsCand->w(); //px, py, pz, E, x, y, z
872 HepLorentzVector p4KsCand;
873 p4KsCand.setX(veewks[0]);
874 p4KsCand.setY(veewks[1]);
875 p4KsCand.setZ(veewks[2]);
876 p4KsCand.setE(veewks[3]);
877 HepLorentzVector p4wks = p4KsCand;
878
879 utility util;
880
881 //by default: px, py, pz, E, chisq, ifok
882 vector<double> vp4ks_new = util.UpdatedKsIfo(aKsCand,vtxsvc,m_useVFrefine);
883 HepLorentzVector p4wks_new;
884 p4wks_new.setX(vp4ks_new[0]);
885 p4wks_new.setY(vp4ks_new[1]);
886 p4wks_new.setZ(vp4ks_new[2]);
887 p4wks_new.setE(vp4ks_new[3]);
888
889 p4Ds = p4Ds - p4wks + p4wks_new;
890
891 }
892
893 }
894
895 p4Ds_New = p4Ds;
896
897 double mass = p4Ds_New.m();
898 int charge = (*it).particle().charge();
899 recDTag->setp4(p4Ds_New);
900 p4Ds_New.boost(-m_beta);
901 double mbc2_CMS = ebeam*ebeam - p4Ds_New.v().mag2();
902 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
903 double deltaE_CMS = p4Ds_New.t() - ebeam;
904
905 int tag=(*it).particle().userTag();
906 if(tag==1)
907 recDTag->settype( EvtRecDTag::Tight );
908 else
909 recDTag->settype( EvtRecDTag::Loose );
910 recDTag->setcharge(charge);
911 recDTag->setcharm(charge);
912 recDTag->setnumOfChildren(numofchildren);
913 recDTag->setmass(mass);
914 recDTag->setmBC(mbc_CMS);
915 recDTag->setbeamE(ebeam);
916 recDTag->setdeltaE(deltaE_CMS);
917
918
919}
920
921
922
923
924
925
double mass
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
DsSelector dsSelector
double ebeam
ObjectVector< EvtRecDTag > EvtRecDTagCol
Definition EvtRecDTag.h:329
EvtRecTrackCol::iterator EvtRecTrackIterator
LocalEptoPiPiEta3PiSelector eptoPiPiEta3PiSelector
LocalEptoPiPiEtaSelector eptoPiPiEtaSelector
LocalEptoRhoGamSelector eptoRhoGamSelector
LocalEtatoGGSelector etatoGGSelector
LocalEtatoPiPiPi0Selector etatoPiPiPi0Selector
LocalKaonSelector kaonSelector
LocalKsSelector ksSelector
LocalPhotonSelector photonSelector
LocalPi0Selector pi0Selector
LocalPionSelector pionSelector
LocalRhotoPiPiSelector rhotoPiPiSelector
IMessageSvc * msgSvc()
virtual const DecayEvidence & decay() const
virtual const EvtRecTrack * photon() const
virtual const EvtRecTrack * track() 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
StatusCode initialize()
vector< string > getlist(string &filename)
void updateKsInfo(CDDecayList::iterator, double, int, EvtRecDTag *, vector< int >, IVertexDbSvc *, bool)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void saveDsInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
void setbeta(Hep3Vector beta)
Definition DsSelector.h:14
void setebeam(double ebeam)
Definition DsSelector.h:13
@ kDstoPiPiPiPi0Pi0
Definition EvtRecDTag.h:176
@ kDstoPiPi0EPPiPiEta
Definition EvtRecDTag.h:187
@ kDstoPiEPPiPiEtaPiPiPi0
Definition EvtRecDTag.h:189
@ kDstoPiEtaPiPiPi0
Definition EvtRecDTag.h:182
@ kDstoPiPiPiEtaPiPiPi0
Definition EvtRecDTag.h:184
@ kDstoKsKminusPiPi
Definition EvtRecDTag.h:168
@ kDstoPiPi0EtaPiPiPi0
Definition EvtRecDTag.h:183
@ kDstoPiPi0EPPiPiEtaPiPiPi0
Definition EvtRecDTag.h:190
@ kDstoPiPiPiPiPiPi0
Definition EvtRecDTag.h:175
@ kDstoPiPi0EPRhoGam
Definition EvtRecDTag.h:193
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)
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
float charge
double double double * p4
Definition qcdloop1.h:77