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