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