BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
ChargedDReconstruction.cxx
Go to the documentation of this file.
1//
2// All D+ decay modes Reconstruction
3//
4#include "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 "VertexDbSvc/IVertexDbSvc.h"
29
31
32#include "utility.h"
33
34using namespace Event;
35
37//*******************************************************************************************
39 ISvcLocator* pSvcLocator )
40 : Algorithm( name, pSvcLocator ) {
41 // Declare the properties
42 declareProperty( "debug", m_debug = false );
43 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
44 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
45 declareProperty( "UseVertexfit", m_usevertexfit = false );
46 declareProperty( "BeamE", m_beamE = 1.8865 );
47 declareProperty( "DpList", m_decaylist = "test.txt" );
48 declareProperty( "UseVFRefine", m_useVFrefine = true );
49 declareProperty( "UseBFieldCorr", m_useBFC = true );
50}
51
52//******************************************************************************************
54 MsgStream log( msgSvc(), name() );
55 log << MSG::INFO << "in initialize()" << endmsg;
56
57 m_irun = -100;
58 m_beta.setX( 0.011 );
59 m_beta.setY( 0 );
60 m_beta.setZ( 0 );
61 chanlist = getlist( m_decaylist );
62
63 auto sc = toolSvc()->retrieveTool( "LocalPionSelector", m_pionSelector );
64 sc &= toolSvc()->retrieveTool( "LocalKaonSelector", m_kaonSelector );
65 sc &= toolSvc()->retrieveTool( "LocalPhotonSelector", m_photonSelector );
66 sc &= toolSvc()->retrieveTool( "LocalKsSelector", m_ksSelector );
67 sc &= toolSvc()->retrieveTool( "LocalPi0Selector", m_pi0Selector );
68 sc &= toolSvc()->retrieveTool( "LocalEtatoGGSelector", m_etatoGGSelector );
69 sc &= toolSvc()->retrieveTool( "ChargedDSelector", m_chargedDSelector );
70 sc &= toolSvc()->retrieveTool( "LocalEptoPiPiEtaSelector", m_eptoPiPiEtaSelector );
71 sc &= toolSvc()->retrieveTool( "LocalRhotoPiPiSelector", m_rhotoPiPiSelector );
72 sc &= toolSvc()->retrieveTool( "LocalEptoRhoGamSelector", m_eptoRhoGamSelector );
73
74 if ( !sc.isSuccess() )
75 {
76 log << MSG::FATAL << "Could not retrieve all selector tools!" << endmsg;
77 return sc;
78 }
79
80 return StatusCode::SUCCESS;
81}
82
83//********************************************************************************************
85 MsgStream log( msgSvc(), name() );
86 log << MSG::INFO << "in finalize()" << endmsg;
87
88 chanlist.clear();
89
90 return StatusCode::SUCCESS;
91}
92
93StatusCode ChargedDReconstruction::registerEvtRecDTagCol( EvtRecDTagCol* aNewEvtRecDTagCol,
94 MsgStream& log ) {
95 StatusCode sc =
96 eventSvc()->registerObject( "/Event/EvtRec/EvtRecDTagCol", aNewEvtRecDTagCol );
97 if ( sc != StatusCode::SUCCESS )
98 { log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endmsg; }
99 return sc;
100}
101
102//*********************************************************************************************
104 MsgStream log( msgSvc(), name() );
105 log << MSG::INFO << "in execute()" << endmsg;
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() << endmsg;
121 log << MSG::DEBUG << "ncharg, nneu, tottks = " << recEvent->totalCharged() << " , "
122 << recEvent->totalNeutral() << " , " << recEvent->totalTracks() << endmsg;
123
124 EvtRecTrackIterator charged_begin = recTrackCol->begin();
125 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
126
127 EvtRecTrackIterator neutral_begin = recTrackCol->begin() + recEvent->totalCharged();
128 EvtRecTrackIterator neutral_end = recTrackCol->begin() + recEvent->totalTracks();
129
130 SmartDataPtr<EvtRecPi0Col> recPi0Col( eventSvc(), "/Event/EvtRec/EvtRecPi0Col" );
131 if ( !recPi0Col )
132 {
133 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endmsg;
134 return StatusCode::FAILURE;
135 }
136
137 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol( eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol" );
138 if ( !recEtaToGGCol )
139 {
140 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endmsg;
141 return StatusCode::FAILURE;
142 }
143
144 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol( eventSvc(),
145 "/Event/EvtRec/EvtRecVeeVertexCol" );
146 if ( !recVeeVertexCol )
147 {
148 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endmsg;
149 return StatusCode::FAILURE;
150 }
151
152 SmartDataPtr<EvtRecDTagCol> recDTagCol( eventSvc(), EventModel::EvtRec::EvtRecDTagCol );
153 if ( !recDTagCol )
154 {
155 log << MSG::FATAL << "EvtRecDTagCol is not registered yet" << endmsg;
156 return StatusCode::FAILURE;
157 }
158
159 // get primary vertex from db
160 Hep3Vector xorigin( 0, 0, 0 );
161
162 IVertexDbSvc* vtxsvc;
163 sc = serviceLocator()->service( "VertexDbSvc", vtxsvc );
164 if ( vtxsvc->isVertexValid() )
165 {
166
167 // vertex[0] = vx; vertex[1]= vy; vertex[2] = vz;
168 double* vertex = vtxsvc->PrimaryVertex();
169 xorigin.setX( vertex[0] );
170 xorigin.setY( vertex[1] );
171 xorigin.setZ( vertex[2] );
172 }
173
174 utility util;
175
176 // registered in DTag.cxx
177 /*
178 if (!recDTagCol) {
179 recDTagCol = new EvtRecDTagCol;
180 sc = registerEvtRecDTagCol(recDTagCol, log);
181 if (sc != StatusCode::SUCCESS) {
182 return sc;
183 }
184 }
185 */
186
187 /////////////////////////////
188 // reconstruct particle lists
189 /////////////////////////////
190 m_pionSelector->setpidtype( 0 );
191 m_kaonSelector->setpidtype( 0 );
192 CDChargedPionList pionList( charged_begin, charged_end, *m_pionSelector );
193 CDChargedKaonList kaonList( charged_begin, charged_end, *m_kaonSelector );
194 CDPhotonList photonList( neutral_begin, neutral_end, *m_photonSelector );
195
196 CDKsList ksList( *m_ksSelector );
197 dc_fill( ksList, recVeeVertexCol->begin(), recVeeVertexCol->end() );
198
199 // do a secondary vertex fit and cut on the results
201 for ( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end();
202 ++ksit )
203 {
204 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( ( *ksit ).particle().navKshort() );
205
206 if ( m_useVFrefine ) { fitinfo[ks] = util.SecondaryVFitref( ks, vtxsvc ); }
207 else { fitinfo[ks] = util.SecondaryVFit( ks, vtxsvc ); }
208 }
209
210 CDPi0List pi0List( *m_pi0Selector );
211 dc_fill( pi0List, recPi0Col->begin(), recPi0Col->end() );
212
213 CDEtaList etaList( *m_etatoGGSelector );
214 dc_fill( etaList, recEtaToGGCol->begin(), recEtaToGGCol->end() );
215
216 // pion/kaon list with PID
217 m_pionSelector->setpidtype( 1 );
218 m_kaonSelector->setpidtype( 1 );
219 CDChargedPionList pionList_tight( charged_begin, charged_end, *m_pionSelector );
220 CDChargedKaonList kaonList_tight( charged_begin, charged_end, *m_kaonSelector );
221
222 int run = eventHeader->runNumber();
223 m_ievt = eventHeader->eventNumber();
224 m_nChrg = recEvent->totalCharged();
225 m_nNeu = recEvent->totalNeutral();
226 m_nPion = pionList.size();
227 m_nKaon = kaonList.size();
228 m_nPi0 = pi0List.size();
229 m_nKs = ksList.size();
230
231 ///////////////////////
232 // get beam energy and beta
233 ///////////////////////
234
235 if ( m_ReadBeamEFromDB && m_irun != run )
236 {
237 m_irun = run;
238 if ( m_usecalibBeamE ) m_readDb.setcalib( true );
239 m_beamE = m_readDb.getbeamE( m_irun, m_beamE );
240 if ( run > 0 ) m_beta = m_readDb.getbeta();
241 // cout<<"use beam E from data base:"<<m_beamE<<endl;
242 }
243 double ebeam = m_beamE;
244
245 //////////////////////////////
246 // reconstruct decay lists
247 /////////////////////////////
248
249 for ( int list = 0; list < chanlist.size(); list++ )
250 {
251
252 string channel = chanlist[list];
253 vector<int> numchan;
254 m_chargedDSelector->setebeam( ebeam );
255 m_chargedDSelector->setbeta( m_beta );
256 CDDecayList decaylist( *m_chargedDSelector );
257
258 // K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5
259 // the fist element of the vector stands for decay mode,
260 // the rest will be particles, and size of the vector minus 1 will be number of daughers.
261
262 if ( channel == "DptoKPiPi" )
263 {
264 numchan.push_back( EvtRecDTag::kDptoKPiPi );
265 numchan.push_back( 1 );
266 numchan.push_back( 2 );
267 numchan.push_back( 2 );
268 decaylist = kaonList.minus() * pionList.plus() * pionList.plus();
269 }
270 else if ( channel == "DptoKPiPiPi0" )
271 {
272 numchan.push_back( EvtRecDTag::kDptoKPiPiPi0 );
273 numchan.push_back( 1 );
274 numchan.push_back( 2 );
275 numchan.push_back( 2 );
276 numchan.push_back( 3 );
277 decaylist = kaonList.minus() * pionList.plus() * pionList.plus() * pi0List;
278 }
279 else if ( channel == "DptoKsPi" )
280 {
281 numchan.push_back( EvtRecDTag::kDptoKsPi );
282 numchan.push_back( 5 );
283 numchan.push_back( 2 );
284 decaylist = ksList * pionList.plus();
285 }
286 else if ( channel == "DptoKsPiPi0" )
287 {
288 numchan.push_back( EvtRecDTag::kDptoKsPiPi0 );
289 numchan.push_back( 5 );
290 numchan.push_back( 2 );
291 numchan.push_back( 3 );
292 decaylist = ksList * pionList.plus() * pi0List;
293 }
294 else if ( channel == "DptoKsPiPiPi" )
295 {
296 numchan.push_back( EvtRecDTag::kDptoKsPiPiPi );
297 numchan.push_back( 5 );
298 numchan.push_back( 2 );
299 numchan.push_back( 2 );
300 numchan.push_back( 2 );
301 decaylist = ksList * pionList.plus() * pionList.plus() * pionList.minus();
302 }
303 else if ( channel == "DptoKKPi" )
304 {
305 numchan.push_back( EvtRecDTag::kDptoKKPi );
306 numchan.push_back( 1 );
307 numchan.push_back( 1 );
308 numchan.push_back( 2 );
309 decaylist = kaonList.minus() * kaonList.plus() * pionList.plus();
310 }
311 else if ( channel == "DptoPiPi0" )
312 {
313 numchan.push_back( EvtRecDTag::kDptoPiPi0 );
314 numchan.push_back( 2 );
315 numchan.push_back( 3 );
316 decaylist = pionList.plus() * pi0List;
317 }
318 else if ( channel == "DptoKPi0" )
319 {
320 numchan.push_back( EvtRecDTag::kDptoKPi0 );
321 numchan.push_back( 1 );
322 numchan.push_back( 3 );
323 decaylist = kaonList.plus() * pi0List;
324 }
325 else if ( channel == "DptoKsK" )
326 {
327 numchan.push_back( EvtRecDTag::kDptoKsK );
328 numchan.push_back( 5 );
329 numchan.push_back( 1 );
330 decaylist = ksList * kaonList.plus();
331 }
332 else if ( channel == "DptoPiPiPi" )
333 {
334 numchan.push_back( EvtRecDTag::kDptoPiPiPi );
335 numchan.push_back( 2 );
336 numchan.push_back( 2 );
337 numchan.push_back( 2 );
338 decaylist = pionList.plus() * pionList.plus() * pionList.minus();
339 }
340 else if ( channel == "DptoPiPi0Pi0" )
341 {
342 numchan.push_back( EvtRecDTag::kDptoPiPi0Pi0 );
343 numchan.push_back( 2 );
344 numchan.push_back( 3 );
345 numchan.push_back( 3 );
346 decaylist = pionList.plus() * pi0List * pi0List;
347 }
348 else if ( channel == "DptoKsKsPi" )
349 {
350 numchan.push_back( EvtRecDTag::kDptoKsKsPi );
351 numchan.push_back( 5 );
352 numchan.push_back( 5 );
353 numchan.push_back( 2 );
354 decaylist = ksList * ksList * pionList.plus();
355 }
356 else if ( channel == "DptoKsKPi0" )
357 {
358 numchan.push_back( EvtRecDTag::kDptoKsKPi0 );
359 numchan.push_back( 5 );
360 numchan.push_back( 1 );
361 numchan.push_back( 3 );
362 decaylist = ksList * kaonList.plus() * pi0List;
363 }
364 else if ( channel == "DptoKsKsK" )
365 {
366 numchan.push_back( EvtRecDTag::kDptoKsKsK );
367 numchan.push_back( 5 );
368 numchan.push_back( 5 );
369 numchan.push_back( 1 );
370 decaylist = ksList * ksList * kaonList.plus();
371 }
372 else if ( channel == "DptoPiPiPiPi0" )
373 {
374 numchan.push_back( EvtRecDTag::kDptoPiPiPiPi0 );
375 numchan.push_back( 2 );
376 numchan.push_back( 2 );
377 numchan.push_back( 2 );
378 numchan.push_back( 3 );
379 decaylist = pionList.plus() * pionList.plus() * pionList.minus() * pi0List;
380 }
381 else if ( channel == "DptoKsPiPi0Pi0" )
382 {
383 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Pi0 );
384 numchan.push_back( 5 );
385 numchan.push_back( 2 );
386 numchan.push_back( 3 );
387 numchan.push_back( 3 );
388 decaylist = ksList * pionList.plus() * pi0List * pi0List;
389 }
390 else if ( channel == "DptoKsKplusPiPi" )
391 {
392 numchan.push_back( EvtRecDTag::kDptoKsKplusPiPi );
393 numchan.push_back( 5 );
394 numchan.push_back( 1 );
395 numchan.push_back( 2 );
396 numchan.push_back( 2 );
397 decaylist = ksList * kaonList.plus() * pionList.plus() * pionList.minus();
398 }
399 else if ( channel == "DptoKsKminusPiPi" )
400 {
401 numchan.push_back( EvtRecDTag::kDptoKsKminusPiPi );
402 numchan.push_back( 5 );
403 numchan.push_back( 1 );
404 numchan.push_back( 2 );
405 numchan.push_back( 2 );
406 decaylist = ksList * kaonList.minus() * pionList.plus() * pionList.plus();
407 }
408 else if ( channel == "DptoKKPiPi0" )
409 {
410 numchan.push_back( EvtRecDTag::kDptoKKPiPi0 );
411 numchan.push_back( 1 );
412 numchan.push_back( 1 );
413 numchan.push_back( 2 );
414 numchan.push_back( 3 );
415 decaylist = kaonList.minus() * kaonList.plus() * pionList.plus() * pi0List;
416 }
417 else if ( channel == "DptoPiPiPiPiPi" )
418 {
419 numchan.push_back( EvtRecDTag::kDptoPiPiPiPiPi );
420 numchan.push_back( 2 );
421 numchan.push_back( 2 );
422 numchan.push_back( 2 );
423 numchan.push_back( 2 );
424 numchan.push_back( 2 );
425 decaylist = pionList.plus() * pionList.plus() * pionList.plus() * pionList.minus() *
426 pionList.minus();
427 }
428 else if ( channel == "DptoKPiPiPiPi" )
429 {
430 numchan.push_back( EvtRecDTag::kDptoKPiPiPiPi );
431 numchan.push_back( 1 );
432 numchan.push_back( 2 );
433 numchan.push_back( 2 );
434 numchan.push_back( 2 );
435 numchan.push_back( 2 );
436 decaylist = kaonList.minus() * pionList.plus() * pionList.plus() * pionList.plus() *
437 pionList.minus();
438 }
439 else if ( channel == "DptoPiEta" )
440 {
441 numchan.push_back( EvtRecDTag::kDptoPiEta );
442 numchan.push_back( 2 );
443 numchan.push_back( 4 );
444 decaylist = pionList.plus() * etaList;
445 }
446 else if ( channel == "DptoKsPiEta" )
447 {
448 numchan.push_back( EvtRecDTag::kDptoKsPiEta );
449 numchan.push_back( 5 );
450 numchan.push_back( 2 );
451 numchan.push_back( 4 );
452 decaylist = ksList * pionList.plus() * etaList;
453 }
454 else if ( channel == "DptoKPiPiPi0Pi0" )
455 {
456 numchan.push_back( EvtRecDTag::kDptoKPiPiPi0Pi0 );
457 numchan.push_back( 1 );
458 numchan.push_back( 2 );
459 numchan.push_back( 2 );
460 numchan.push_back( 3 );
461 numchan.push_back( 3 );
462 decaylist = kaonList.minus() * pionList.plus() * pionList.plus() * pi0List * pi0List;
463 }
464 else if ( channel == "DptoKsPiPiPiPi0" )
465 {
466 numchan.push_back( EvtRecDTag::kDptoKsPiPiPiPi0 );
467 numchan.push_back( 5 );
468 numchan.push_back( 2 );
469 numchan.push_back( 2 );
470 numchan.push_back( 2 );
471 numchan.push_back( 3 );
472 decaylist = ksList * pionList.plus() * pionList.plus() * pionList.minus() * pi0List;
473 }
474 else if ( channel == "DptoKsPiPi0Pi0Pi0" )
475 {
476 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Pi0Pi0 );
477 numchan.push_back( 5 );
478 numchan.push_back( 2 );
479 numchan.push_back( 3 );
480 numchan.push_back( 3 );
481 numchan.push_back( 3 );
482 decaylist = ksList * pionList.plus() * pi0List * pi0List * pi0List;
483 }
484 else if ( channel == "DptoKsPiEPPiPiEta" )
485 {
486 numchan.push_back( EvtRecDTag::kDptoKsPiEPPiPiEta );
487 numchan.push_back( 5 );
488 numchan.push_back( 2 );
489 numchan.push_back( 6 );
490 CDDecayList epList( *m_eptoPiPiEtaSelector );
491 epList = pionList.plus() * pionList.minus() * etaList;
492 decaylist = ksList * pionList.plus() * epList;
493 }
494 else if ( channel == "DptoKsPiEPRhoGam" )
495 {
496 numchan.push_back( EvtRecDTag::kDptoKsPiEPRhoGam );
497 numchan.push_back( 5 );
498 numchan.push_back( 2 );
499 numchan.push_back( 7 );
500 CDDecayList rhoList( *m_rhotoPiPiSelector );
501 rhoList = pionList.plus() * pionList.minus();
502 CDDecayList epList( *m_eptoRhoGamSelector );
503 epList = rhoList * photonList;
504 decaylist = ksList * pionList.plus() * epList;
505 }
506 else if ( channel == "DptoKPiPiEta" )
507 {
508 numchan.push_back( EvtRecDTag::kDptoKPiPiEta );
509 numchan.push_back( 1 );
510 numchan.push_back( 2 );
511 numchan.push_back( 2 );
512 numchan.push_back( 4 );
513 decaylist = kaonList.minus() * pionList.plus() * pionList.plus() * etaList;
514 }
515 else if ( channel == "DptoKsPiPi0Eta" )
516 {
517 numchan.push_back( EvtRecDTag::kDptoKsPiPi0Eta );
518 numchan.push_back( 5 );
519 numchan.push_back( 2 );
520 numchan.push_back( 3 );
521 numchan.push_back( 4 );
522 decaylist = ksList * pionList.plus() * pi0List * etaList;
523 }
524 else if ( channel == "DptoKsKKPi" )
525 {
526 numchan.push_back( EvtRecDTag::kDptoKsKKPi );
527 numchan.push_back( 5 );
528 numchan.push_back( 1 );
529 numchan.push_back( 1 );
530 numchan.push_back( 2 );
531 decaylist = ksList * kaonList.minus() * kaonList.plus() * pionList.plus();
532 }
533 else if ( channel == "DptoPiPiPiPi0Pi0" )
534 {
535 numchan.push_back( EvtRecDTag::kDptoPiPiPiPi0Pi0 );
536 numchan.push_back( 2 );
537 numchan.push_back( 2 );
538 numchan.push_back( 2 );
539 numchan.push_back( 3 );
540 numchan.push_back( 3 );
541 decaylist = pionList.plus() * pionList.plus() * pionList.minus() * pi0List * pi0List;
542 }
543 else if ( channel == "DptoPiPi0Eta" )
544 {
545 numchan.push_back( EvtRecDTag::kDptoPiPi0Eta );
546 numchan.push_back( 2 );
547 numchan.push_back( 3 );
548 numchan.push_back( 4 );
549 decaylist = pionList.plus() * pi0List * etaList;
550 }
551 else if ( channel == "DptoPiPiPiEta" )
552 {
553 numchan.push_back( EvtRecDTag::kDptoPiPiPiEta );
554 numchan.push_back( 2 );
555 numchan.push_back( 2 );
556 numchan.push_back( 2 );
557 numchan.push_back( 4 );
558 decaylist = pionList.plus() * pionList.plus() * pionList.minus() * etaList;
559 }
560 else if ( channel == "DptoPiEtaEta" )
561 {
562 numchan.push_back( EvtRecDTag::kDptoPiEtaEta );
563 numchan.push_back( 2 );
564 numchan.push_back( 4 );
565 numchan.push_back( 4 );
566 decaylist = pionList.plus() * etaList * etaList;
567 }
568 else if ( channel == "DptoPiEPPiPiEta" )
569 {
570 numchan.push_back( EvtRecDTag::kDptoPiEPPiPiEta );
571 numchan.push_back( 2 );
572 numchan.push_back( 6 );
573 CDDecayList epList( *m_eptoPiPiEtaSelector );
574 epList = pionList.plus() * pionList.minus() * etaList;
575 decaylist = pionList.plus() * epList;
576 }
577 else if ( channel == "DptoPiEPRhoGam" )
578 {
579 numchan.push_back( EvtRecDTag::kDptoPiEPRhoGam );
580 numchan.push_back( 2 );
581 numchan.push_back( 7 );
582 CDDecayList rhoList( *m_rhotoPiPiSelector );
583 rhoList = pionList.plus() * pionList.minus();
584 CDDecayList epList( *m_eptoRhoGamSelector );
585 epList = rhoList * photonList;
586 decaylist = pionList.plus() * epList;
587 }
588 else if ( channel == "DptoPiPi0EPPiPiEta" )
589 {
590 numchan.push_back( EvtRecDTag::kDptoPiPi0EPPiPiEta );
591 numchan.push_back( 2 );
592 numchan.push_back( 3 );
593 numchan.push_back( 6 );
594 CDDecayList epList( *m_eptoPiPiEtaSelector );
595 epList = pionList.plus() * pionList.minus() * etaList;
596 decaylist = pionList.plus() * pi0List * epList;
597 }
598 else if ( channel == "DptoPiPi0EPRhoGam" )
599 {
600 numchan.push_back( EvtRecDTag::kDptoPiPi0EPRhoGam );
601 numchan.push_back( 2 );
602 numchan.push_back( 3 );
603 numchan.push_back( 7 );
604 CDDecayList rhoList( *m_rhotoPiPiSelector );
605 rhoList = pionList.plus() * pionList.minus();
606 CDDecayList epList( *m_eptoRhoGamSelector );
607 epList = rhoList * photonList;
608 decaylist = pionList.plus() * pi0List * epList;
609 }
610 else if ( channel == "DptoKsKsPiPi0" )
611 {
612 numchan.push_back( EvtRecDTag::kDptoKsKsPiPi0 );
613 numchan.push_back( 5 );
614 numchan.push_back( 5 );
615 numchan.push_back( 2 );
616 numchan.push_back( 3 );
617 decaylist = ksList * ksList * pionList.plus() * pi0List;
618 }
619 else if ( channel == "DptoKsKPi0Pi0" )
620 {
621 numchan.push_back( EvtRecDTag::kDptoKsKPi0Pi0 );
622 numchan.push_back( 5 );
623 numchan.push_back( 1 );
624 numchan.push_back( 3 );
625 numchan.push_back( 3 );
626 decaylist = ksList * kaonList.plus() * pi0List * pi0List;
627 }
628 else if ( channel == "DptoKsKEta" )
629 {
630 numchan.push_back( EvtRecDTag::kDptoKsKEta );
631 numchan.push_back( 5 );
632 numchan.push_back( 1 );
633 numchan.push_back( 4 );
634 decaylist = ksList * kaonList.plus() * etaList;
635 }
636 else if ( channel == "DptoKKPiPiPi" )
637 {
638 numchan.push_back( EvtRecDTag::kDptoKKPiPiPi );
639 numchan.push_back( 1 );
640 numchan.push_back( 1 );
641 numchan.push_back( 2 );
642 numchan.push_back( 2 );
643 numchan.push_back( 2 );
644 decaylist = kaonList.minus() * kaonList.plus() * pionList.plus() * pionList.minus() *
645 pionList.plus();
646 }
647 else if ( channel == "DptoKpPiPi" )
648 {
649 numchan.push_back( EvtRecDTag::kDptoKpPiPi );
650 numchan.push_back( 1 );
651 numchan.push_back( 2 );
652 numchan.push_back( 2 );
653 decaylist = kaonList.plus() * pionList.plus() * pionList.minus();
654 }
655 else if ( channel == "DptoKpPi0Pi0" )
656 {
657 numchan.push_back( EvtRecDTag::kDptoKpPi0Pi0 );
658 numchan.push_back( 1 );
659 numchan.push_back( 3 );
660 numchan.push_back( 3 );
661 decaylist = kaonList.plus() * pi0List * pi0List;
662 }
663 else if ( channel == "DptoKpPi0Eta" )
664 {
665 numchan.push_back( EvtRecDTag::kDptoKpPi0Eta );
666 numchan.push_back( 1 );
667 numchan.push_back( 3 );
668 numchan.push_back( 4 );
669 decaylist = kaonList.plus() * pi0List * etaList;
670 }
671 else if ( channel == "DptoKpPiPiPi0" )
672 {
673 numchan.push_back( EvtRecDTag::kDptoKpPiPiPi0 );
674 numchan.push_back( 1 );
675 numchan.push_back( 2 );
676 numchan.push_back( 2 );
677 numchan.push_back( 3 );
678 decaylist = kaonList.plus() * pionList.plus() * pionList.minus() * pi0List;
679 }
680 else if ( channel == "DptoKpPiPiEta" )
681 {
682 numchan.push_back( EvtRecDTag::kDptoKpPiPiEta );
683 numchan.push_back( 1 );
684 numchan.push_back( 2 );
685 numchan.push_back( 2 );
686 numchan.push_back( 4 );
687 decaylist = kaonList.plus() * pionList.plus() * pionList.minus() * etaList;
688 }
689
690 CDDecayList::iterator D_begin = decaylist.particle_begin();
691 CDDecayList::iterator D_end = decaylist.particle_end();
692
693 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ )
694 {
695
696 EvtRecDTag* recDTag = new EvtRecDTag;
697 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
698
699 vector<int> trackid, showerid;
700 vector<int> kaonid, pionid;
701 int numofchildren = numchan.size() - 1;
702
703 for ( int i = 0; i < numofchildren; i++ )
704 {
705
706 const CDCandidate& daughter = ( *it ).particle().child( i );
707
708 if ( numchan[i + 1] == 1 )
709 {
710 const EvtRecTrack* track = daughter.track();
711 trackid.push_back( track->trackId() );
712 kaonid.push_back( track->trackId() );
713 }
714 else if ( numchan[i + 1] == 2 )
715 {
716 const EvtRecTrack* track = daughter.track();
717 trackid.push_back( track->trackId() );
718 pionid.push_back( track->trackId() );
719 }
720 else if ( numchan[i + 1] == 3 )
721 {
722 const EvtRecTrack* hiEnGamma = daughter.navPi0()->hiEnGamma();
723 const EvtRecTrack* loEnGamma = daughter.navPi0()->loEnGamma();
724 showerid.push_back( hiEnGamma->trackId() );
725 showerid.push_back( loEnGamma->trackId() );
726 }
727 else if ( numchan[i + 1] == 4 )
728 {
729 const EvtRecTrack* hiEnGamma = daughter.navEta()->hiEnGamma();
730 const EvtRecTrack* loEnGamma = daughter.navEta()->loEnGamma();
731 showerid.push_back( hiEnGamma->trackId() );
732 showerid.push_back( loEnGamma->trackId() );
733 }
734 else if ( numchan[i + 1] == 5 )
735 {
736 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
737 // fill fit info
738 recDTag->addToFitInfo( aKsCand->mass(), fitinfo[aKsCand][0], fitinfo[aKsCand][1],
739 fitinfo[aKsCand][2] );
740 // fill tracks
741 EvtRecTrack* pion1Trk = aKsCand->daughter( 0 );
742 EvtRecTrack* pion2Trk = aKsCand->daughter( 1 );
743 trackid.push_back( pion1Trk->trackId() );
744 trackid.push_back( pion2Trk->trackId() );
745 }
746 else if ( numchan[i + 1] == 6 )
747 {
748 const CDCandidate& apion = daughter.decay().child( 0 );
749 const CDCandidate& spion = daughter.decay().child( 1 );
750 const CDCandidate& eta = daughter.decay().child( 2 );
751 const EvtRecTrack* apiontrk = apion.track();
752 const EvtRecTrack* spiontrk = spion.track();
753 const EvtRecTrack* hiEnGamma = eta.navEta()->hiEnGamma();
754 const EvtRecTrack* loEnGamma = eta.navEta()->loEnGamma();
755
756 trackid.push_back( apiontrk->trackId() );
757 trackid.push_back( spiontrk->trackId() );
758 showerid.push_back( hiEnGamma->trackId() );
759 showerid.push_back( loEnGamma->trackId() );
760 }
761 else if ( numchan[i + 1] == 7 )
762 {
763 const CDCandidate& rho = daughter.decay().child( 0 );
764 const CDCandidate& gamma = daughter.decay().child( 1 );
765 const CDCandidate& apion = rho.decay().child( 0 );
766 const CDCandidate& spion = rho.decay().child( 1 );
767
768 const EvtRecTrack* apiontrk = apion.track();
769 const EvtRecTrack* spiontrk = spion.track();
770 const EvtRecTrack* gammatrk = gamma.photon();
771
772 trackid.push_back( apiontrk->trackId() );
773 trackid.push_back( spiontrk->trackId() );
774 showerid.push_back( gammatrk->trackId() );
775 }
776
777 } // end of filling track and shower ids
778
779 saveDpInfo( it, ebeam, numofchildren, recDTag );
780 if ( m_useBFC )
781 {
782 // After use BFC package, we need to update information of Ks and Lambda to calculate
783 // D.
784 updateKsInfo( it, ebeam, numofchildren, recDTag, numchan, vtxsvc, m_useVFrefine );
785 }
786 savetrack( trackid, showerid, charged_begin, charged_end, neutral_begin, neutral_end,
787 recDTag );
788 pidtag( kaonid, pionid, kaonList_tight, pionList_tight, recDTag );
789
790 if ( m_usevertexfit )
791 {
792 if ( m_debug ) cout << "beforevfit:" << endl;
793
794 HepLorentzVector p4change_vfit;
795
796 if ( m_useVFrefine )
797 { p4change_vfit = util.vfitref( channel, kaonid, pionid, xorigin, charged_begin ); }
798 else { p4change_vfit = util.vfit( channel, kaonid, pionid, xorigin, charged_begin ); }
799
800 recDTag->setp4( recDTag->p4() + p4change_vfit );
801 }
802
803 trackid.clear();
804 showerid.clear();
805 kaonid.clear();
806 pionid.clear();
807
808 // write dtag object out
809 recDTagCol->push_back( recDTag );
810
811 } // end of decaylist iterator
812
813 numchan.clear();
814
815 } // end of reconstrucing all D+ decay lists
816
817 return StatusCode::SUCCESS;
818}
819
821 int numofchildren, EvtRecDTag* recDTag ) {
822
823 double mass = ( *it ).particle().mass();
824 int charge = ( *it ).particle().charge();
825 HepLorentzVector p4( ( *it ).particle().momentum(), ( *it ).particle().energy() );
826 recDTag->setp4( p4 );
827
828 p4.boost( -m_beta );
829 double mbc2_CMS = ebeam * ebeam - p4.v().mag2();
830 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
831 double deltaE_CMS = p4.t() - ebeam;
832
833 if ( ( *it ).particle().userTag() == 1 ) recDTag->settype( EvtRecDTag::Tight );
834 else recDTag->settype( EvtRecDTag::Loose );
835 recDTag->setcharge( charge );
836 recDTag->setcharm( charge );
837 recDTag->setnumOfChildren( numofchildren );
838 recDTag->setmass( mass );
839 recDTag->setmBC( mbc_CMS );
840 recDTag->setbeamE( ebeam );
841 recDTag->setdeltaE( deltaE_CMS );
842}
843
844void ChargedDReconstruction::savetrack( vector<int> trackid, vector<int> showerid,
845 EvtRecTrackIterator charged_begin,
846 EvtRecTrackIterator charged_end,
847 EvtRecTrackIterator neutral_begin,
848 EvtRecTrackIterator neutral_end,
849 EvtRecDTag* recDTag ) {
850
851 vector<EvtRecTrackIterator> trktemp;
852 vector<EvtRecTrackIterator> shrtemp;
853
854 // fill tracks
855 for ( EvtRecTrackIterator trk = charged_begin; trk < charged_end; trk++ )
856 {
857
858 bool isothertrack = true;
859 for ( int i = 0; i < trackid.size(); i++ )
860 {
861 if ( ( *trk )->trackId() == trackid[i] )
862 {
863 trktemp.push_back( trk );
864 isothertrack = false;
865 break;
866 }
867 }
868 if ( isothertrack ) recDTag->addOtherTrack( *trk );
869 }
870 for ( int i = 0; i < trackid.size(); i++ )
871 {
872 for ( int j = 0; j < trktemp.size(); j++ )
873 {
874 EvtRecTrackIterator trk = trktemp[j];
875 if ( ( *trk )->trackId() == trackid[i] ) recDTag->addTrack( *trktemp[j] );
876 }
877 }
878
879 // fill showers
880 for ( EvtRecTrackIterator shr = neutral_begin; shr < neutral_end; shr++ )
881 {
882 bool isothershower = true;
883 for ( int i = 0; i < showerid.size(); i++ )
884 {
885 if ( ( *shr )->trackId() == showerid[i] )
886 {
887 shrtemp.push_back( shr );
888 isothershower = false;
889 break;
890 }
891 }
892 if ( isothershower ) recDTag->addOtherShower( *shr );
893 }
894
895 for ( int i = 0; i < showerid.size(); i++ )
896 {
897 for ( int j = 0; j < shrtemp.size(); j++ )
898 {
899 EvtRecTrackIterator shr = shrtemp[j];
900 if ( ( *shr )->trackId() == showerid[i] ) recDTag->addShower( *shrtemp[j] );
901 }
902 }
903}
904
905void ChargedDReconstruction::pidtag( vector<int> kaonid, vector<int> pionid,
906 CDChargedKaonList& kaonList, CDChargedPionList& pionList,
907 EvtRecDTag* recDTag ) {
908
909 bool iskaon = false, ispion = false;
910
911 // save track ids which passed pion/kaon cuts
912
913 for ( CDChargedKaonList::iterator kit = kaonList.particle_begin();
914 kit != kaonList.particle_end(); kit++ )
915 { recDTag->addKaonId( ( *kit ).particle().track() ); }
916
917 for ( CDChargedPionList::iterator pit = pionList.particle_begin();
918 pit != pionList.particle_end(); pit++ )
919 { recDTag->addPionId( ( *pit ).particle().track() ); }
920
921 /*
922 for(int i=0; i<kaonid.size(); i++){
923 bool ithkaon=false;
924 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit !=
925 kaonList.particle_end(); kit++) { if((*kit).particle().track()->trackId()==kaonid[i]){
926 ithkaon=true;
927 break;
928 }
929 }
930 if(!ithkaon) break;
931 if(i==kaonid.size()-1)
932 iskaon=true;
933 }
934
935 for(int i=0; i<pionid.size(); i++){
936 bool ithpion=false;
937 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit !=
938 pionList.particle_end(); pit++) { if((*pit).particle().track()->trackId()==pionid[i]){
939 ithpion=true;
940 break;
941 }
942 }
943 if(!ithpion) break;
944 if(i==pionid.size()-1)
945 ispion=true;
946 }
947
948
949 if( iskaon && ispion)
950 recDTag->settype( EvtRecDTag::Tight );
951 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
952 recDTag->settype( EvtRecDTag::Tight );
953 else if( kaonid.size()==0 && pionid.size()==0 )
954 recDTag->settype( EvtRecDTag::Tight );
955 */
956}
957
958vector<string> ChargedDReconstruction::getlist( string& filename ) {
959
960 string channel;
961 vector<string> temp;
962
963 ifstream inFile;
964
965 inFile.open( filename.c_str() );
966 if ( !inFile )
967 {
968 cout << "Unable to open decay list file";
969 exit( 1 ); // terminate with error
970 }
971
972 while ( inFile >> channel ) { temp.push_back( channel ); }
973
974 inFile.close();
975
976 return temp;
977}
978
980 int numofchildren, EvtRecDTag* recDTag,
981 vector<int> numchan, IVertexDbSvc* vtxsvc,
982 bool m_useVFrefine ) {
983
984 // This function is used to update the information of Ks and Lambda. Because DTagAlg package
985 // uses the list from VeeVertexFit package which input the information of Ks and Lambda while
986 // reconstructing the dst, the information of them don't use the MBF correction and we need
987 // to update them here.
988
989 // If we don't use BF Correctoin, please close this function.
990
991 HepLorentzVector p4Dp( ( *it ).particle().momentum(), ( *it ).particle().energy() );
992 HepLorentzVector p4Dp_New;
993
994 for ( int i = 0; i < numofchildren; i++ )
995 {
996 const CDCandidate& daughter = ( *it ).particle().child( i );
997
998 if ( numchan[i + 1] == 5 )
999 {
1000 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
1001 HepVector veewks = aKsCand->w(); // px, py, pz, E, x, y, z
1002 HepLorentzVector p4KsCand;
1003 p4KsCand.setX( veewks[0] );
1004 p4KsCand.setY( veewks[1] );
1005 p4KsCand.setZ( veewks[2] );
1006 p4KsCand.setE( veewks[3] );
1007 HepLorentzVector p4wks = p4KsCand;
1008
1009 utility util;
1010
1011 // by default: px, py, pz, E, chisq, ifok
1012 vector<double> vp4ks_new = util.UpdatedKsIfo( aKsCand, vtxsvc, m_useVFrefine );
1013 HepLorentzVector p4wks_new;
1014 p4wks_new.setX( vp4ks_new[0] );
1015 p4wks_new.setY( vp4ks_new[1] );
1016 p4wks_new.setZ( vp4ks_new[2] );
1017 p4wks_new.setE( vp4ks_new[3] );
1018
1019 p4Dp = p4Dp - p4wks + p4wks_new;
1020 }
1021 }
1022
1023 p4Dp_New = p4Dp;
1024
1025 double mass = p4Dp_New.m();
1026 int charge = ( *it ).particle().charge();
1027 recDTag->setp4( p4Dp_New );
1028 p4Dp_New.boost( -m_beta );
1029 double mbc2_CMS = ebeam * ebeam - p4Dp_New.v().mag2();
1030 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
1031 double deltaE_CMS = p4Dp_New.t() - ebeam;
1032
1033 int tag = ( *it ).particle().userTag();
1034 if ( tag == 1 ) recDTag->settype( EvtRecDTag::Tight );
1035 else recDTag->settype( EvtRecDTag::Loose );
1036 recDTag->setcharge( charge );
1037 recDTag->setcharm( charge );
1038 recDTag->setnumOfChildren( numofchildren );
1039 recDTag->setmass( mass );
1040 recDTag->setmBC( mbc_CMS );
1041 recDTag->setbeamE( ebeam );
1042 recDTag->setdeltaE( deltaE_CMS );
1043}
DECLARE_COMPONENT(BesBdkRc)
double mass
DCFillableChargedList< CDChargedKaon > CDChargedKaonList
DCFillableChargedList< CDChargedPion > CDChargedPionList
DCDecayList< CDDecay, CDDecay::CandidateClass > CDDecayList
DCFillableNeutralList< CDPhoton > CDPhotonList
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
ObjectVector< EvtRecDTag > EvtRecDTagCol
EvtRecTrackCol::iterator EvtRecTrackIterator
double ebeam
IMessageSvc * msgSvc()
virtual const DecayEvidence & decay() const
virtual const EvtRecTrack * photon() const
virtual const EvtRecTrack * track() const
virtual const EvtRecVeeVertex * navKshort() const
virtual const EvtRecPi0 * navPi0() const
virtual const EvtRecEtaToGG * navEta() const
const CDCandidate & child(unsigned int aPosition) const
Definition CDDecay.cxx:231
void updateKsInfo(CDDecayList::iterator, double, int, EvtRecDTag *, vector< int >, IVertexDbSvc *, bool)
void saveDpInfo(CDDecayList::iterator, double, int, EvtRecDTag *)
vector< string > getlist(string &filename)
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
ChargedDReconstruction(const std::string &name, ISvcLocator *pSvcLocator)
void addOtherTrack(const SmartRef< EvtRecTrack > track)
void addToFitInfo(double ksmass, double chi2, double length, double error)
void addOtherShower(const SmartRef< EvtRecTrack > shower)
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
void addPionId(const SmartRef< EvtRecTrack > pionId)
void addShower(const SmartRef< EvtRecTrack > shower)
void addTrack(const SmartRef< EvtRecTrack > track)
const EvtRecTrack * loEnGamma() const
const EvtRecTrack * hiEnGamma() const
virtual bool isVertexValid()=0
virtual double * PrimaryVertex()=0
dchain::MuteWholeCandidateItr< CandidateClass > iterator
HepLorentzVector vfitref(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:59
vector< double > SecondaryVFitref(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:203
vector< double > UpdatedKsIfo(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc, bool m_useVFrefine)
Definition utility.cxx:793
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:425
vector< double > SecondaryVFit(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:567