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