BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
ChargedDReconstruction.cxx
Go to the documentation of this file.
1//
2// All D+ 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
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"
36
48
49
50#include "DTagAlg/utility.h"
51
52using namespace Event;
53
54DECLARE_COMPONENT(ChargedDReconstruction)
55//*******************************************************************************************
56ChargedDReconstruction::ChargedDReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
57Algorithm(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 = 1.8865 );
64 declareProperty( "DpList", 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 ChargedDReconstruction::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
141 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
142 if ( ! recEtaToGGCol ) {
143 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
144 return StatusCode::FAILURE;
145 }
146
147
148 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
149 if ( ! recVeeVertexCol ) {
150 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
151 return StatusCode::FAILURE;
152 }
153
154
155 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
156 if (!recDTagCol) {
157 log << MSG::FATAL << "EvtRecDTagCol is not registered yet" << endreq;
158 return StatusCode::FAILURE;
159 }
160
161
162
163 //get primary vertex from db
164 Hep3Vector xorigin(0,0,0);
165
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
177 utility util;
178
179
180 //registered in 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 /////////////////////////////
196 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
197 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
198 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
199
200 CDKsList ksList(ksSelector);
201 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
202
203 // do a secondary vertex fit and cut on the results
204 map<EvtRecVeeVertex*, vector< double > > fitinfo;
205 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
206 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
207
208 if(m_useVFrefine){
209 fitinfo[ks] = util.SecondaryVFitref(ks, vtxsvc);
210 }else{
211 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
212 }
213
214 }
215
216 CDPi0List pi0List(pi0Selector);
217 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
218
219 CDEtaList etaList(etatoGGSelector);
220 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
221
222 //pion/kaon list with PID
225 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
226 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
227
228
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 ///////////////////////
241 // get beam energy and beta
242 ///////////////////////
243
244
245 if(m_ReadBeamEFromDB && m_irun!=run){
246 m_irun=run;
247 if(m_usecalibBeamE)
248 m_readDb.setcalib(true);
249 m_beamE=m_readDb.getbeamE(m_irun,m_beamE);
250 if(run>0)
251 m_beta=m_readDb.getbeta();
252 // cout<<"use beam E from data base:"<<m_beamE<<endl;
253 }
254 double ebeam=m_beamE;
255
256
257 //////////////////////////////
258 //reconstruct decay lists
259 /////////////////////////////
260
261
262 for(int list=0;list<chanlist.size();list++){
263
264 string channel=chanlist[list];
265 vector<int> numchan;
268 CDDecayList decaylist(chargedDSelector);
269
270 //K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5
271 //the fist element of the vector stands for decay mode,
272 //the rest will be particles, and size of the vector minus 1 will be number of daughers.
273
274 if(channel=="DptoKPiPi") {
275 numchan.push_back( EvtRecDTag::kDptoKPiPi );
276 numchan.push_back(1);
277 numchan.push_back(2);
278 numchan.push_back(2);
279 decaylist=kaonList.minus()* pionList.plus()* pionList.plus();
280 }
281 else if(channel=="DptoKPiPiPi0") {
282 numchan.push_back( EvtRecDTag::kDptoKPiPiPi0 );
283 numchan.push_back(1);
284 numchan.push_back(2);
285 numchan.push_back(2);
286 numchan.push_back(3);
287 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pi0List;
288 }
289 else if(channel=="DptoKsPi") {
290 numchan.push_back( EvtRecDTag::kDptoKsPi );
291 numchan.push_back(5);
292 numchan.push_back(2);
293 decaylist=ksList* pionList.plus();
294 }
295 else if(channel=="DptoKsPiPi0") {
296 numchan.push_back( EvtRecDTag::kDptoKsPiPi0 );
297 numchan.push_back(5);
298 numchan.push_back(2);
299 numchan.push_back(3);
300 decaylist=ksList* pionList.plus()* pi0List;
301 }
302 else if(channel=="DptoKsPiPiPi") {
303 numchan.push_back( EvtRecDTag::kDptoKsPiPiPi );
304 numchan.push_back(5);
305 numchan.push_back(2);
306 numchan.push_back(2);
307 numchan.push_back(2);
308 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus();
309 }
310 else if(channel=="DptoKKPi") {
311 numchan.push_back( EvtRecDTag::kDptoKKPi );
312 numchan.push_back(1);
313 numchan.push_back(1);
314 numchan.push_back(2);
315 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus();
316 }
317 else if(channel=="DptoPiPi0") {
318 numchan.push_back( EvtRecDTag::kDptoPiPi0 );
319 numchan.push_back(2);
320 numchan.push_back(3);
321 decaylist=pionList.plus()* pi0List;
322 }
323 else if(channel=="DptoKPi0") {
324 numchan.push_back( EvtRecDTag::kDptoKPi0 );
325 numchan.push_back(1);
326 numchan.push_back(3);
327 decaylist=kaonList.plus()* pi0List;
328 }
329 else if(channel=="DptoKsK") {
330 numchan.push_back( EvtRecDTag::kDptoKsK );
331 numchan.push_back(5);
332 numchan.push_back(1);
333 decaylist=ksList* kaonList.plus();
334 }
335 else if(channel=="DptoPiPiPi") {
336 numchan.push_back( EvtRecDTag::kDptoPiPiPi );
337 numchan.push_back(2);
338 numchan.push_back(2);
339 numchan.push_back(2);
340 decaylist=pionList.plus()* pionList.plus()* pionList.minus() ;
341 }
342 else if(channel=="DptoPiPi0Pi0") {
343 numchan.push_back( EvtRecDTag::kDptoPiPi0Pi0 );
344 numchan.push_back(2);
345 numchan.push_back(3);
346 numchan.push_back(3);
347 decaylist=pionList.plus()* pi0List* pi0List;
348 }
349 else if(channel=="DptoKsKsPi") {
350 numchan.push_back( EvtRecDTag::kDptoKsKsPi );
351 numchan.push_back(5);
352 numchan.push_back(5);
353 numchan.push_back(2);
354 decaylist=ksList* ksList* pionList.plus();
355 }
356 else if(channel=="DptoKsKPi0") {
357 numchan.push_back( EvtRecDTag::kDptoKsKPi0 );
358 numchan.push_back(5);
359 numchan.push_back(1);
360 numchan.push_back(3);
361 decaylist=ksList* kaonList.plus()* pi0List;
362 }
363 else if(channel=="DptoKsKsK") {
364 numchan.push_back( EvtRecDTag::kDptoKsKsK );
365 numchan.push_back(5);
366 numchan.push_back(5);
367 numchan.push_back(1);
368 decaylist=ksList* ksList* kaonList.plus();
369 }
370 else if(channel=="DptoPiPiPiPi0") {
371 numchan.push_back( EvtRecDTag::kDptoPiPiPiPi0 );
372 numchan.push_back(2);
373 numchan.push_back(2);
374 numchan.push_back(2);
375 numchan.push_back(3);
376 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
377 }
378 else if(channel=="DptoKsPiPi0Pi0") {
379 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Pi0 );
380 numchan.push_back(5);
381 numchan.push_back(2);
382 numchan.push_back(3);
383 numchan.push_back(3);
384 decaylist=ksList* pionList.plus()* pi0List* pi0List;
385 }
386 else if(channel=="DptoKsKplusPiPi") {
387 numchan.push_back( EvtRecDTag::kDptoKsKplusPiPi );
388 numchan.push_back(5);
389 numchan.push_back(1);
390 numchan.push_back(2);
391 numchan.push_back(2);
392 decaylist=ksList* kaonList.plus()* pionList.plus()* pionList.minus();
393 }
394 else if(channel=="DptoKsKminusPiPi") {
395 numchan.push_back( EvtRecDTag::kDptoKsKminusPiPi );
396 numchan.push_back(5);
397 numchan.push_back(1);
398 numchan.push_back(2);
399 numchan.push_back(2);
400 decaylist=ksList* kaonList.minus()* pionList.plus()* pionList.plus();
401 }
402 else if(channel=="DptoKKPiPi0") {
403 numchan.push_back( EvtRecDTag::kDptoKKPiPi0 );
404 numchan.push_back(1);
405 numchan.push_back(1);
406 numchan.push_back(2);
407 numchan.push_back(3);
408 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pi0List;
409 }
410 else if(channel=="DptoPiPiPiPiPi") {
411 numchan.push_back( EvtRecDTag::kDptoPiPiPiPiPi );
412 numchan.push_back(2);
413 numchan.push_back(2);
414 numchan.push_back(2);
415 numchan.push_back(2);
416 numchan.push_back(2);
417 decaylist=pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
418 }
419 else if(channel=="DptoKPiPiPiPi") {
420 numchan.push_back( EvtRecDTag::kDptoKPiPiPiPi );
421 numchan.push_back(1);
422 numchan.push_back(2);
423 numchan.push_back(2);
424 numchan.push_back(2);
425 numchan.push_back(2);
426 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.plus()* pionList.minus();
427 }
428 else if(channel=="DptoPiEta") {
429 numchan.push_back( EvtRecDTag::kDptoPiEta );
430 numchan.push_back(2);
431 numchan.push_back(4);
432 decaylist=pionList.plus()* etaList;
433 }
434 else if(channel=="DptoKsPiEta") {
435 numchan.push_back( EvtRecDTag::kDptoKsPiEta );
436 numchan.push_back(5);
437 numchan.push_back(2);
438 numchan.push_back(4);
439 decaylist=ksList* pionList.plus()* etaList;
440 }
441 else if(channel=="DptoKPiPiPi0Pi0") {
442 numchan.push_back( EvtRecDTag::kDptoKPiPiPi0Pi0 );
443 numchan.push_back(1);
444 numchan.push_back(2);
445 numchan.push_back(2);
446 numchan.push_back(3);
447 numchan.push_back(3);
448 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pi0List* pi0List;
449 }
450 else if(channel=="DptoKsPiPiPiPi0") {
451 numchan.push_back( EvtRecDTag::kDptoKsPiPiPiPi0 );
452 numchan.push_back(5);
453 numchan.push_back(2);
454 numchan.push_back(2);
455 numchan.push_back(2);
456 numchan.push_back(3);
457 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
458 }
459 else if(channel=="DptoKsPiPi0Pi0Pi0") {
460 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Pi0Pi0 );
461 numchan.push_back(5);
462 numchan.push_back(2);
463 numchan.push_back(3);
464 numchan.push_back(3);
465 numchan.push_back(3);
466 decaylist=ksList* pionList.plus()* pi0List* pi0List* pi0List;
467 }
468 else if(channel=="DptoKsPiEPPiPiEta"){
469 numchan.push_back( EvtRecDTag::kDptoKsPiEPPiPiEta);
470 numchan.push_back(5);
471 numchan.push_back(2);
472 numchan.push_back(6);
474 epList=pionList.plus()* pionList.minus()* etaList;
475 decaylist= ksList* pionList.plus()* epList;
476 }
477 else if(channel=="DptoKsPiEPRhoGam"){
478 numchan.push_back( EvtRecDTag::kDptoKsPiEPRhoGam);
479 numchan.push_back(5);
480 numchan.push_back(2);
481 numchan.push_back(7);
483 rhoList=pionList.plus()* pionList.minus();
485 epList=rhoList* photonList;
486 decaylist= ksList* pionList.plus()* epList;
487 }
488 else if(channel=="DptoKPiPiEta") {
489 numchan.push_back( EvtRecDTag::kDptoKPiPiEta );
490 numchan.push_back(1);
491 numchan.push_back(2);
492 numchan.push_back(2);
493 numchan.push_back(4);
494 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* etaList;
495 }
496 else if(channel=="DptoKsPiPi0Eta") {
497 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Eta );
498 numchan.push_back(5);
499 numchan.push_back(2);
500 numchan.push_back(3);
501 numchan.push_back(4);
502 decaylist=ksList* pionList.plus()* pi0List* etaList;
503 }
504 else if(channel=="DptoKsKKPi") {
505 numchan.push_back( EvtRecDTag::kDptoKsKKPi );
506 numchan.push_back(5);
507 numchan.push_back(1);
508 numchan.push_back(1);
509 numchan.push_back(2);
510 decaylist=ksList* kaonList.minus()* kaonList.plus()* pionList.plus();
511 }
512 else if(channel=="DptoPiPiPiPi0Pi0") {
513 numchan.push_back( EvtRecDTag::kDptoPiPiPiPi0Pi0 );
514 numchan.push_back(2);
515 numchan.push_back(2);
516 numchan.push_back(2);
517 numchan.push_back(3);
518 numchan.push_back(3);
519 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pi0List* pi0List;
520 }
521 else if(channel=="DptoPiPi0Eta") {
522 numchan.push_back( EvtRecDTag::kDptoPiPi0Eta );
523 numchan.push_back(2);
524 numchan.push_back(3);
525 numchan.push_back(4);
526 decaylist=pionList.plus()* pi0List* etaList;
527 }
528 else if(channel=="DptoPiPiPiEta") {
529 numchan.push_back( EvtRecDTag::kDptoPiPiPiEta );
530 numchan.push_back(2);
531 numchan.push_back(2);
532 numchan.push_back(2);
533 numchan.push_back(4);
534 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* etaList ;
535 }
536 else if(channel=="DptoPiEtaEta") {
537 numchan.push_back( EvtRecDTag::kDptoPiEtaEta );
538 numchan.push_back(2);
539 numchan.push_back(4);
540 numchan.push_back(4);
541 decaylist=pionList.plus()* etaList* etaList;
542 }
543 else if(channel=="DptoPiEPPiPiEta"){
544 numchan.push_back( EvtRecDTag::kDptoPiEPPiPiEta);
545 numchan.push_back(2);
546 numchan.push_back(6);
548 epList=pionList.plus()* pionList.minus()* etaList;
549 decaylist= pionList.plus()* epList;
550 }
551 else if(channel=="DptoPiEPRhoGam"){
552 numchan.push_back( EvtRecDTag::kDptoPiEPRhoGam);
553 numchan.push_back(2);
554 numchan.push_back(7);
556 rhoList=pionList.plus()* pionList.minus();
558 epList=rhoList* photonList;
559 decaylist= pionList.plus()* epList;
560 }
561 else if(channel=="DptoPiPi0EPPiPiEta"){
562 numchan.push_back( EvtRecDTag::kDptoPiPi0EPPiPiEta);
563 numchan.push_back(2);
564 numchan.push_back(3);
565 numchan.push_back(6);
567 epList=pionList.plus()* pionList.minus()* etaList;
568 decaylist= pionList.plus()* pi0List* epList;
569 }
570 else if(channel=="DptoPiPi0EPRhoGam"){
571 numchan.push_back( EvtRecDTag::kDptoPiPi0EPRhoGam);
572 numchan.push_back(2);
573 numchan.push_back(3);
574 numchan.push_back(7);
576 rhoList=pionList.plus()* pionList.minus();
578 epList=rhoList* photonList;
579 decaylist= pionList.plus()* pi0List* epList;
580 }
581 else if(channel=="DptoKsKsPiPi0") {
582 numchan.push_back( EvtRecDTag::kDptoKsKsPiPi0 );
583 numchan.push_back(5);
584 numchan.push_back(5);
585 numchan.push_back(2);
586 numchan.push_back(3);
587 decaylist=ksList* ksList* pionList.plus()* pi0List;
588 }
589 else if(channel=="DptoKsKPi0Pi0") {
590 numchan.push_back( EvtRecDTag::kDptoKsKPi0Pi0 );
591 numchan.push_back(5);
592 numchan.push_back(1);
593 numchan.push_back(3);
594 numchan.push_back(3);
595 decaylist=ksList* kaonList.plus()* pi0List* pi0List;
596 }
597 else if(channel=="DptoKsKEta") {
598 numchan.push_back( EvtRecDTag::kDptoKsKEta );
599 numchan.push_back(5);
600 numchan.push_back(1);
601 numchan.push_back(4);
602 decaylist=ksList* kaonList.plus()* etaList;
603 }
604 else if(channel=="DptoKKPiPiPi") {
605 numchan.push_back( EvtRecDTag::kDptoKKPiPiPi );
606 numchan.push_back(1);
607 numchan.push_back(1);
608 numchan.push_back(2);
609 numchan.push_back(2);
610 numchan.push_back(2);
611 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus()* pionList.plus();
612 }
613 else if(channel=="DptoKpPiPi") {
614 numchan.push_back( EvtRecDTag::kDptoKpPiPi );
615 numchan.push_back(1);
616 numchan.push_back(2);
617 numchan.push_back(2);
618 decaylist=kaonList.plus()* pionList.plus()* pionList.minus();
619 }
620 else if(channel=="DptoKpPi0Pi0") {
621 numchan.push_back( EvtRecDTag::kDptoKpPi0Pi0 );
622 numchan.push_back(1);
623 numchan.push_back(3);
624 numchan.push_back(3);
625 decaylist=kaonList.plus()* pi0List* pi0List;
626 }
627 else if(channel=="DptoKpPi0Eta") {
628 numchan.push_back( EvtRecDTag::kDptoKpPi0Eta );
629 numchan.push_back(1);
630 numchan.push_back(3);
631 numchan.push_back(4);
632 decaylist=kaonList.plus()* pi0List* etaList;
633 }
634 else if(channel=="DptoKpPiPiPi0") {
635 numchan.push_back( EvtRecDTag::kDptoKpPiPiPi0 );
636 numchan.push_back(1);
637 numchan.push_back(2);
638 numchan.push_back(2);
639 numchan.push_back(3);
640 decaylist=kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
641 }
642 else if(channel=="DptoKpPiPiEta") {
643 numchan.push_back( EvtRecDTag::kDptoKpPiPiEta );
644 numchan.push_back(1);
645 numchan.push_back(2);
646 numchan.push_back(2);
647 numchan.push_back(4);
648 decaylist=kaonList.plus()* pionList.plus()* pionList.minus()* etaList;
649 }
650
651
652 CDDecayList::iterator D_begin =decaylist.particle_begin();
653 CDDecayList::iterator D_end =decaylist.particle_end();
654
655 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
656
657 EvtRecDTag* recDTag = new EvtRecDTag;
658 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
659
660 vector<int> trackid, showerid;
661 vector<int> kaonid, pionid;
662 int numofchildren=numchan.size()-1;
663
664 for(int i=0; i< numofchildren;i++){
665
666 const CDCandidate& daughter=(*it).particle().child(i);
667
668 if(numchan[i+1]==1){
669 const EvtRecTrack* track=daughter.track();
670 trackid.push_back(track->trackId());
671 kaonid.push_back(track->trackId());
672 }
673 else if(numchan[i+1]==2){
674 const EvtRecTrack* track=daughter.track();
675 trackid.push_back(track->trackId());
676 pionid.push_back(track->trackId());
677 }
678 else if ( numchan[i+1]==3){
679 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
680 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
681 showerid.push_back(hiEnGamma->trackId());
682 showerid.push_back(loEnGamma->trackId());
683 }
684 else if ( numchan[i+1]==4){
685 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
686 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
687 showerid.push_back(hiEnGamma->trackId());
688 showerid.push_back(loEnGamma->trackId());
689 }
690 else if ( numchan[i+1]==5){
691 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
692 // fill fit info
693 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
694 // fill tracks
695 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
696 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
697 trackid.push_back(pion1Trk->trackId());
698 trackid.push_back(pion2Trk->trackId());
699 }
700 else if (numchan[i+1]==6){
701 const CDCandidate& apion = daughter.decay().child(0);
702 const CDCandidate& spion = daughter.decay().child(1);
703 const CDCandidate& eta = daughter.decay().child(2);
704 const EvtRecTrack* apiontrk = apion.track();
705 const EvtRecTrack* spiontrk = spion.track();
706 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
707 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
708
709 trackid.push_back(apiontrk->trackId());
710 trackid.push_back(spiontrk->trackId());
711 showerid.push_back(hiEnGamma->trackId());
712 showerid.push_back(loEnGamma->trackId());
713
714 }
715 else if (numchan[i+1]==7){
716 const CDCandidate& rho = daughter.decay().child(0);
717 const CDCandidate& gamma = daughter.decay().child(1);
718 const CDCandidate& apion = rho.decay().child(0);
719 const CDCandidate& spion = rho.decay().child(1);
720
721
722 const EvtRecTrack* apiontrk = apion.track();
723 const EvtRecTrack* spiontrk = spion.track();
724 const EvtRecTrack* gammatrk = gamma.photon();
725
726
727 trackid.push_back(apiontrk->trackId());
728 trackid.push_back(spiontrk->trackId());
729 showerid.push_back(gammatrk->trackId());
730
731 }
732
733 }//end of filling track and shower ids
734
735
736 saveDpInfo(it, ebeam, numofchildren, recDTag);
737 if(m_useBFC){
738 //After use BFC package, we need to update information of Ks and Lambda to calculate D.
739 updateKsInfo(it, ebeam, numofchildren, recDTag, numchan, vtxsvc, m_useVFrefine);
740 }
741 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
742 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
743
744 if(m_usevertexfit){
745 if(m_debug) cout<<"beforevfit:"<<endl;
746
747 HepLorentzVector p4change_vfit;
748
749 if(m_useVFrefine){
750 p4change_vfit=util.vfitref(channel, kaonid, pionid, xorigin, charged_begin);
751 }else{
752 p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
753 }
754
755 recDTag->setp4(recDTag->p4()+p4change_vfit);
756
757 }
758
759
760 trackid.clear();
761 showerid.clear();
762 kaonid.clear();
763 pionid.clear();
764
765 //write dtag object out
766 recDTagCol->push_back(recDTag);
767
768 }//end of decaylist iterator
769
770 numchan.clear();
771
772 }//end of reconstrucing all D+ decay lists
773
774 return StatusCode::SUCCESS;
775}
776
777
778void ChargedDReconstruction::saveDpInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag){
779
780 double mass = (*it).particle().mass();
781 int charge= (*it).particle().charge();
782 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
783 recDTag->setp4(p4);
784
785 p4.boost(-m_beta);
786 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
787 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
788 double deltaE_CMS = p4.t() - ebeam;
789
790 if((*it).particle().userTag()==1)
791 recDTag->settype( EvtRecDTag::Tight );
792 else
793 recDTag->settype( EvtRecDTag::Loose );
794 recDTag->setcharge(charge);
795 recDTag->setcharm(charge);
796 recDTag->setnumOfChildren(numofchildren);
797 recDTag->setmass(mass);
798 recDTag->setmBC(mbc_CMS);
799 recDTag->setbeamE(ebeam);
800 recDTag->setdeltaE(deltaE_CMS);
801
802}
803
804void ChargedDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
805 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
806
807 vector<EvtRecTrackIterator> trktemp;
808 vector<EvtRecTrackIterator> shrtemp;
809
810 //fill tracks
811 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
812
813 bool isothertrack=true;
814 for(int i=0; i<trackid.size(); i++){
815 if( (*trk)->trackId()==trackid[i] ){
816 trktemp.push_back(trk);
817 isothertrack=false;
818 break;
819 }
820 }
821 if(isothertrack)
822 recDTag->addOtherTrack(*trk);
823 }
824 for(int i=0; i<trackid.size();i++){
825 for(int j=0; j<trktemp.size(); j++){
826 EvtRecTrackIterator trk=trktemp[j];
827 if( (*trk)->trackId()==trackid[i])
828 recDTag->addTrack(*trktemp[j]);
829 }
830 }
831
832
833 //fill showers
834 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
835 bool isothershower=true;
836 for(int i=0; i<showerid.size(); i++){
837 if( (*shr)->trackId()==showerid[i] ){
838 shrtemp.push_back(shr);
839 isothershower=false;
840 break;
841 }
842 }
843 if(isothershower)
844 recDTag->addOtherShower(*shr);
845 }
846
847 for(int i=0; i<showerid.size();i++){
848 for(int j=0; j<shrtemp.size(); j++){
849 EvtRecTrackIterator shr=shrtemp[j];
850 if( (*shr)->trackId()==showerid[i])
851 recDTag->addShower(*shrtemp[j]);
852 }
853 }
854
855
856}
857
858void ChargedDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
859
860 bool iskaon=false,ispion=false;
861
862
863 // save track ids which passed pion/kaon cuts
864
865 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
866 recDTag->addKaonId( (*kit).particle().track() );
867 }
868
869 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
870 recDTag->addPionId( (*pit).particle().track() );
871 }
872
873
874 /*
875 for(int i=0; i<kaonid.size(); i++){
876 bool ithkaon=false;
877 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
878 if((*kit).particle().track()->trackId()==kaonid[i]){
879 ithkaon=true;
880 break;
881 }
882 }
883 if(!ithkaon) break;
884 if(i==kaonid.size()-1)
885 iskaon=true;
886 }
887
888 for(int i=0; i<pionid.size(); i++){
889 bool ithpion=false;
890 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
891 if((*pit).particle().track()->trackId()==pionid[i]){
892 ithpion=true;
893 break;
894 }
895 }
896 if(!ithpion) break;
897 if(i==pionid.size()-1)
898 ispion=true;
899 }
900
901
902 if( iskaon && ispion)
903 recDTag->settype( EvtRecDTag::Tight );
904 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
905 recDTag->settype( EvtRecDTag::Tight );
906 else if( kaonid.size()==0 && pionid.size()==0 )
907 recDTag->settype( EvtRecDTag::Tight );
908 */
909
910}
911
912
913
914vector<string> ChargedDReconstruction::getlist(string& filename ){
915
916 string channel;
917 vector<string> temp;
918
919 ifstream inFile;
920
921 inFile.open(filename.c_str());
922 if (!inFile) {
923 cout << "Unable to open decay list file";
924 exit(1); // terminate with error
925 }
926
927 while (inFile >> channel) {
928 temp.push_back(channel);
929 }
930
931 inFile.close();
932
933 return temp;
934
935}
936
937
938void ChargedDReconstruction::updateKsInfo(CDDecayList::iterator it, double ebeam, int numofchildren, EvtRecDTag* recDTag, vector<int> numchan, IVertexDbSvc* vtxsvc, bool m_useVFrefine){
939
940//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.
941
942 //If we don't use BF Correctoin, please close this function.
943
944 HepLorentzVector p4Dp((*it).particle().momentum(), (*it).particle().energy());
945 HepLorentzVector p4Dp_New;
946
947 for(int i = 0; i < numofchildren; i++){
948 const CDCandidate& daughter=(*it).particle().child(i);
949
950 if(numchan[i+1]==5){
951 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>(daughter.navKshort());
952 HepVector veewks = aKsCand->w(); //px, py, pz, E, x, y, z
953 HepLorentzVector p4KsCand;
954 p4KsCand.setX(veewks[0]);
955 p4KsCand.setY(veewks[1]);
956 p4KsCand.setZ(veewks[2]);
957 p4KsCand.setE(veewks[3]);
958 HepLorentzVector p4wks = p4KsCand;
959
960 utility util;
961
962 //by default: px, py, pz, E, chisq, ifok
963 vector<double> vp4ks_new = util.UpdatedKsIfo(aKsCand,vtxsvc,m_useVFrefine);
964 HepLorentzVector p4wks_new;
965 p4wks_new.setX(vp4ks_new[0]);
966 p4wks_new.setY(vp4ks_new[1]);
967 p4wks_new.setZ(vp4ks_new[2]);
968 p4wks_new.setE(vp4ks_new[3]);
969
970 p4Dp = p4Dp - p4wks + p4wks_new;
971
972 }
973
974 }
975
976 p4Dp_New = p4Dp;
977
978 double mass = p4Dp_New.m();
979 int charge = (*it).particle().charge();
980 recDTag->setp4(p4Dp_New);
981 p4Dp_New.boost(-m_beta);
982 double mbc2_CMS = ebeam*ebeam - p4Dp_New.v().mag2();
983 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
984 double deltaE_CMS = p4Dp_New.t() - ebeam;
985
986 int tag=(*it).particle().userTag();
987 if(tag==1)
988 recDTag->settype( EvtRecDTag::Tight );
989 else
990 recDTag->settype( EvtRecDTag::Loose );
991 recDTag->setcharge(charge);
992 recDTag->setcharm(charge);
993 recDTag->setnumOfChildren(numofchildren);
994 recDTag->setmass(mass);
995 recDTag->setmBC(mbc_CMS);
996 recDTag->setbeamE(ebeam);
997 recDTag->setdeltaE(deltaE_CMS);
998
999}
1000
1001
1002
1003
ChargedDSelector chargedDSelector
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
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
void updateKsInfo(CDDecayList::iterator, double, int, EvtRecDTag *, vector< int >, IVertexDbSvc *, bool)
void saveDpInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
vector< string > getlist(string &filename)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
void setbeta(Hep3Vector beta)
void setebeam(double ebeam)
@ kDptoKsKminusPiPi
Definition EvtRecDTag.h:128
@ kDptoKsPiEPPiPiEta
Definition EvtRecDTag.h:138
@ kDptoKsPiEPRhoGam
Definition EvtRecDTag.h:139
@ kDptoPiPi0EPPiPiEta
Definition EvtRecDTag.h:149
@ kDptoPiPi0EPRhoGam
Definition EvtRecDTag.h:150
@ kDptoKsPiPi0Pi0Pi0
Definition EvtRecDTag.h:137
@ kDptoPiPiPiPi0Pi0
Definition EvtRecDTag.h:143
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