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