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