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