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