BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
NeutralDReconstruction.cxx
Go to the documentation of this file.
1//
2// All D0 decay modes Reconstruction
3//
4
5#include "GaudiKernel/MsgStream.h"
6#include "GaudiKernel/AlgFactory.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/SmartDataPtr.h"
9#include "GaudiKernel/IDataProviderSvc.h"
10#include "GaudiKernel/PropertyMgr.h"
11
12#include "EventModel/EventHeader.h"
13#include "EventModel/Event.h"
14#include "EventModel/EventModel.h"
15
16#include "EvtRecEvent/EvtRecEvent.h"
17#include "EvtRecEvent/EvtRecTrack.h"
18#include "EvtRecEvent/EvtRecDTag.h"
19#include "EvtRecEvent/EvtRecPi0.h"
20#include "EvtRecEvent/EvtRecEtaToGG.h"
21#include "EvtRecEvent/EvtRecVeeVertex.h"
22
23#include "BesDChain/CDChargedPionList.h"
24#include "BesDChain/CDChargedKaonList.h"
25#include "BesDChain/CDPhotonList.h"
26#include "BesDChain/CDPi0List.h"
27#include "BesDChain/CDEtaList.h"
28#include "BesDChain/CDKsList.h"
29#include "BesDChain/CDDecayList.h"
30
31#include "McTruth/McParticle.h"
32#include "ParticleID/ParticleID.h"
33#include "VertexFit/VertexFit.h"
34#include "VertexFit/SecondVertexFit.h"
35
36
37#include "DTagAlg/NeutralDReconstruction.h"
38#include "DTagAlg/LocalKaonSelector.h"
39#include "DTagAlg/LocalPionSelector.h"
40#include "DTagAlg/LocalKsSelector.h"
41#include "DTagAlg/LocalPi0Selector.h"
42#include "DTagAlg/LocalEtatoGGSelector.h"
43#include "DTagAlg/LocalPhotonSelector.h"
44#include "DTagAlg/LocalEptoPiPiEtaSelector.h"
45#include "DTagAlg/LocalEptoRhoGamSelector.h"
46#include "DTagAlg/LocalRhotoPiPiSelector.h"
47#include "DTagAlg/NeutralDSelector.h"
48
49#include "DTagAlg/utility.h"
50
51#include <fstream>
52
53using namespace Event;
54
55DECLARE_COMPONENT(NeutralDReconstruction)
56//*******************************************************************************************
57NeutralDReconstruction::NeutralDReconstruction(const std::string& name, ISvcLocator* pSvcLocator):
58Algorithm(name, pSvcLocator) {
59 //Declare the properties
60 declareProperty( "debug", m_debug = false );
61 declareProperty( "ReadBeamEFromDB", m_ReadBeamEFromDB = false );
62 declareProperty( "UseCalibBeamE", m_usecalibBeamE = false );
63 declareProperty( "UseVertexfit", m_usevertexfit = false );
64 declareProperty( "BeamE", m_beamE = 1.8865 );
65 declareProperty( "D0List", m_decaylist = "test.txt" );
66 declareProperty("UseVFRefine", m_useVFrefine = true);
67}
68
69//******************************************************************************************
70StatusCode NeutralDReconstruction::initialize() {
71 MsgStream log(msgSvc(), name());
72 log << MSG::INFO << "in initialize()" <<endreq;
73
74 m_irun=-100;
75 m_beta.setX(0.011);
76 m_beta.setY(0);
77 m_beta.setZ(0);
78
79 chanlist=getlist(m_decaylist);
80
81 return StatusCode::SUCCESS;
82}
83
84//********************************************************************************************
85StatusCode NeutralDReconstruction::finalize() {
86 MsgStream log(msgSvc(), name());
87 log << MSG::INFO << "in finalize()" << endreq;
88
89 chanlist.clear();
90
91 return StatusCode::SUCCESS;
92}
93
94StatusCode NeutralDReconstruction::registerEvtRecDTagCol(
95 EvtRecDTagCol* aNewEvtRecDTagCol, MsgStream& log) {
96 StatusCode sc = eventSvc()->registerObject("/Event/EvtRec/EvtRecDTagCol",
97 aNewEvtRecDTagCol);
98 if (sc != StatusCode::SUCCESS) {
99 log << MSG::FATAL << "Could not register EvtRecDTagCol in TDS!" << endreq;
100 }
101 return sc;
102}
103
104//*********************************************************************************************
105StatusCode NeutralDReconstruction::execute() {
106 MsgStream log(msgSvc(), name());
107 log << MSG::INFO << "in execute()" << endreq;
108
109 StatusCode sc;
110
111 //////////////////
112 // Read REC data
113 /////////////////
114 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
115 int event= eventHeader->eventNumber();
116 // if ( m_debug || ( (event & 0x3FF) == 0 ) )
117 //std::cout << "event: " << event << std::endl;
118
119 SmartDataPtr<EvtRecEvent> recEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
120 SmartDataPtr<EvtRecTrackCol> recTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
121 log << MSG::DEBUG << "run and event = " << eventHeader->runNumber()
122 << " " << eventHeader->eventNumber() << endreq;
123 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
124 << recEvent->totalCharged() << " , "
125 << recEvent->totalNeutral() << " , "
126 << recEvent->totalTracks() <<endreq;
127
128 EvtRecTrackIterator charged_begin = recTrackCol->begin();
129 EvtRecTrackIterator charged_end = charged_begin + recEvent->totalCharged();
130
131 EvtRecTrackIterator neutral_begin = recTrackCol->begin()+recEvent->totalCharged();
132 EvtRecTrackIterator neutral_end = recTrackCol->begin()+recEvent->totalTracks();
133
134
135 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
136 if ( ! recPi0Col ) {
137 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
138 return StatusCode::FAILURE;
139 }
140
141
142
143 SmartDataPtr<EvtRecEtaToGGCol> recEtaToGGCol(eventSvc(), "/Event/EvtRec/EvtRecEtaToGGCol");
144 if ( ! recEtaToGGCol ) {
145 log << MSG::FATAL << "Could not find EvtRecEtaToGGCol" << endreq;
146 return StatusCode::FAILURE;
147 }
148
149
150 SmartDataPtr<EvtRecVeeVertexCol> recVeeVertexCol(eventSvc(), "/Event/EvtRec/EvtRecVeeVertexCol");
151 if ( ! recVeeVertexCol ) {
152 log << MSG::FATAL << "Could not find EvtRecVeeVertexCol" << endreq;
153 return StatusCode::FAILURE;
154 }
155
156
157 SmartDataPtr<EvtRecDTagCol> recDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
158 if (!recDTagCol) {
159 log << MSG::FATAL << "EvtRecDTagCol is not registered in TDS" << endreq;
160 return StatusCode::FAILURE;
161 }
162
163 //get primary vertex from db
164 Hep3Vector xorigin(0,0,0);
165 IVertexDbSvc* vtxsvc;
166 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
167 if (vtxsvc->isVertexValid()) {
168
169 //vertex[0] = vx; vertex[1]= vy; vertex[2] = vz;
170 double* vertex = vtxsvc->PrimaryVertex();
171 xorigin.setX(vertex[0]);
172 xorigin.setY(vertex[1]);
173 xorigin.setZ(vertex[2]);
174 }
175 utility util;
176
177
178
179 //registered in the parent file DTag.cxx
180 /*
181 if (!recDTagCol) {
182 recDTagCol = new EvtRecDTagCol;
183 sc = registerEvtRecDTagCol(recDTagCol, log);
184 if (sc != StatusCode::SUCCESS) {
185 return sc;
186 }
187 }
188 */
189
190 /////////////////////////////
191 //reconstruct particle lists
192 /////////////////////////////
193
194 pionSelector.setpidtype(0);
195 kaonSelector.setpidtype(0);
196 CDChargedPionList pionList(charged_begin, charged_end, pionSelector);
197 CDChargedKaonList kaonList(charged_begin, charged_end, kaonSelector);
198 CDPhotonList photonList(neutral_begin, neutral_end, photonSelector);
199
200 CDKsList ksList(ksSelector);
201 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
202
203 // do a secondary vertex fit and cut on the results
204 map<EvtRecVeeVertex*, vector< double > > fitinfo;
205 for( CDKsList::iterator ksit = ksList.particle_begin(); ksit != ksList.particle_end(); ++ksit ){
206 EvtRecVeeVertex* ks = const_cast<EvtRecVeeVertex*>( (*ksit).particle().navKshort() );
207
208 if(m_useVFrefine){
209 fitinfo[ks] = util.SecondaryVFitref(ks, vtxsvc);
210 }else{
211 fitinfo[ks] = util.SecondaryVFit(ks, vtxsvc);
212 }
213
214 }
215
216 CDPi0List pi0List(pi0Selector);
217 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
218
219 CDEtaList etaList(etatoGGSelector);
220 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
221
222
223 //pion/kaon list with PID
224 pionSelector.setpidtype(1);
225 kaonSelector.setpidtype(1);
226 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
227 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
228
229 int run = eventHeader->runNumber();
230 m_ievt = eventHeader->eventNumber();
231 m_nChrg = recEvent->totalCharged();
232 m_nNeu = recEvent->totalNeutral();
233 m_nPion = pionList.size();
234 m_nKaon = kaonList.size();
235 m_nPi0 = pi0List.size();
236 m_nKs = ksList.size();
237
238 ///////////////////////
239 // get beam energy and beta
240 ///////////////////////
241
242 if(m_ReadBeamEFromDB && m_irun!=run){
243 m_irun=run;
244 if(m_usecalibBeamE)
245 m_readDb.setcalib(true);
246 m_beamE=m_readDb.getbeamE(m_irun, m_beamE);
247 if(run>0)
248 m_beta=m_readDb.getbeta();
249
250 //cout<<"use beam E from data base:"<<m_beamE<<endl;
251 }
252 double ebeam=m_beamE;
253
254 //////////////////////////////
255 //reconstruct decay lists
256 /////////////////////////////
257
258
259
260 for(int list=0;list<chanlist.size();list++){
261
262 string channel=chanlist[list];
263 vector<int> numchan;
264 neutralDSelector.setebeam(ebeam);
265 neutralDSelector.setbeta(m_beta);
266 CDDecayList decaylist(neutralDSelector);
267
268 bool isFlavorMode=false;
269
270 //K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5, Eta'(pipiEta):6, Eta'(rhogam):7
271 //the fist element of the vector stands for decay mode,
272 //the rest will be particles, and size of the vector minus 1 will be number of daughers.
273 //the first particle in the vector will be used to get charm
274
275 if(channel=="D0toKPi") {
276 numchan.push_back( EvtRecDTag::kD0toKPi );
277 numchan.push_back(1);
278 numchan.push_back(2);
279 decaylist=kaonList.minus() * pionList.plus();
280 isFlavorMode=true;
281 }
282 else if(channel=="D0toKPiPi0"){
283 numchan.push_back( EvtRecDTag::kD0toKPiPi0);
284 numchan.push_back(1);
285 numchan.push_back(2);
286 numchan.push_back(3);
287 decaylist=kaonList.minus() * pionList.plus() * pi0List;
288 isFlavorMode=true;
289 }
290 else if(channel=="D0toKPiPi0Pi0"){
291 numchan.push_back( EvtRecDTag::kD0toKPiPi0Pi0);
292 numchan.push_back(1);
293 numchan.push_back(2);
294 numchan.push_back(3);
295 numchan.push_back(3);
296 decaylist=kaonList.minus() * pionList.plus() * pi0List * pi0List;
297 isFlavorMode=true;
298 }
299 else if(channel=="D0toKPiPiPi"){
300 numchan.push_back( EvtRecDTag::kD0toKPiPiPi);
301 numchan.push_back(1);
302 numchan.push_back(2);
303 numchan.push_back(2);
304 numchan.push_back(2);
305 decaylist=kaonList.minus()* pionList.plus()* pionList.plus() * pionList.minus();
306 isFlavorMode=true;
307 }
308 else if(channel=="D0toKPiPiPiPi0"){
309 numchan.push_back( EvtRecDTag::kD0toKPiPiPiPi0);
310 numchan.push_back(1);
311 numchan.push_back(2);
312 numchan.push_back(2);
313 numchan.push_back(2);
314 numchan.push_back(3);
315 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
316 isFlavorMode=true;
317 }
318 else if(channel=="D0toKPiEta"){
319 numchan.push_back( EvtRecDTag::kD0toKPiEta);
320 numchan.push_back(1);
321 numchan.push_back(2);
322 numchan.push_back(4);
323 decaylist=kaonList.minus() * pionList.plus() * etaList;
324 isFlavorMode=true;
325 }
326 else if(channel=="D0toKsKPi"){
327 numchan.push_back( EvtRecDTag::kD0toKsKPi);
328 numchan.push_back(5);
329 numchan.push_back(1);
330 numchan.push_back(2);
331 decaylist=ksList* kaonList.minus()* pionList.plus();
332 }
333 else if(channel=="D0toKsKPiPi0"){
334 numchan.push_back( EvtRecDTag::kD0toKsKPiPi0);
335 numchan.push_back(5);
336 numchan.push_back(1);
337 numchan.push_back(2);
338 numchan.push_back(3);
339 decaylist=ksList* kaonList.minus() * pionList.plus()* pi0List;
340 }
341 else if(channel=="D0toKsPiPi"){
342 numchan.push_back( EvtRecDTag::kD0toKsPiPi);
343 numchan.push_back(5);
344 numchan.push_back(2);
345 numchan.push_back(2);
346 decaylist=ksList* pionList.plus()* pionList.minus();
347 }
348 else if(channel=="D0toKsPiPiPi0"){
349 numchan.push_back( EvtRecDTag::kD0toKsPiPiPi0);
350 numchan.push_back(5);
351 numchan.push_back(2);
352 numchan.push_back(2);
353 numchan.push_back(3);
354 decaylist=ksList* pionList.plus()* pionList.minus()* pi0List;
355 }
356 else if(channel=="D0toKsPi0"){
357 numchan.push_back( EvtRecDTag::kD0toKsPi0);
358 numchan.push_back(5);
359 numchan.push_back(3);
360 decaylist=ksList* pi0List;
361 }
362 else if(channel=="D0toPiPiPi0"){
363 numchan.push_back( EvtRecDTag::kD0toPiPiPi0);
364 numchan.push_back(2);
365 numchan.push_back(2);
366 numchan.push_back(3);
367 decaylist=pionList.plus()* pionList.minus()* pi0List;
368 }
369 else if(channel=="D0toPiPi"){
370 numchan.push_back( EvtRecDTag::kD0toPiPi);
371 numchan.push_back(2);
372 numchan.push_back(2);
373 decaylist=pionList.plus()* pionList.minus();
374 }
375 else if(channel=="D0toKK"){
376 numchan.push_back( EvtRecDTag::kD0toKK);
377 numchan.push_back(1);
378 numchan.push_back(1);
379 decaylist=kaonList.minus()* kaonList.plus();
380 }
381 else if(channel=="D0toKKPi0"){
382 numchan.push_back( EvtRecDTag::kD0toKKPi0);
383 numchan.push_back(1);
384 numchan.push_back(1);
385 numchan.push_back(3);
386 decaylist=kaonList.minus()* kaonList.plus()* pi0List;
387 }
388 else if(channel=="D0toPi0Pi0"){
389 numchan.push_back( EvtRecDTag::kD0toPi0Pi0);
390 numchan.push_back(3);
391 numchan.push_back(3);
392 decaylist=pi0List* pi0List;
393 }
394 else if(channel=="D0toKsKs"){
395 numchan.push_back( EvtRecDTag::kD0toKsKs);
396 numchan.push_back(5);
397 numchan.push_back(5);
398 decaylist=ksList* ksList;
399 }
400 else if(channel=="D0toKsKsPi0"){
401 numchan.push_back( EvtRecDTag::kD0toKsKsPi0);
402 numchan.push_back(5);
403 numchan.push_back(5);
404 numchan.push_back(3);
405 decaylist=ksList* ksList* pi0List;
406 }
407 else if(channel=="D0toKsPi0Pi0"){
408 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0);
409 numchan.push_back(5);
410 numchan.push_back(3);
411 numchan.push_back(3);
412 decaylist=ksList* pi0List* pi0List;
413 }
414 else if(channel=="D0toKsKK"){
415 numchan.push_back( EvtRecDTag::kD0toKsKK);
416 numchan.push_back(5);
417 numchan.push_back(1);
418 numchan.push_back(1);
419 decaylist=ksList* kaonList.minus()* kaonList.plus();
420 }
421 else if(channel=="D0toKsEta"){
422 numchan.push_back( EvtRecDTag::kD0toKsEta);
423 numchan.push_back(5);
424 numchan.push_back(4);
425 decaylist=ksList* etaList;
426 }
427 else if(channel=="D0toPi0Pi0Pi0"){
428 numchan.push_back( EvtRecDTag::kD0toPi0Pi0Pi0);
429 numchan.push_back(3);
430 numchan.push_back(3);
431 numchan.push_back(3);
432 decaylist=pi0List* pi0List* pi0List;
433 }
434 else if(channel=="D0toKsKsKs"){
435 numchan.push_back( EvtRecDTag::kD0toKsKsKs);
436 numchan.push_back(5);
437 numchan.push_back(5);
438 numchan.push_back(5);
439 decaylist=ksList* ksList* ksList;
440 }
441 else if(channel=="D0toPiPiPiPi"){
442 numchan.push_back( EvtRecDTag::kD0toPiPiPiPi);
443 numchan.push_back(2);
444 numchan.push_back(2);
445 numchan.push_back(2);
446 numchan.push_back(2);
447 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
448 }
449 else if(channel=="D0toPiPiPi0Pi0"){
450 numchan.push_back( EvtRecDTag::kD0toPiPiPi0Pi0);
451 numchan.push_back(2);
452 numchan.push_back(2);
453 numchan.push_back(3);
454 numchan.push_back(3);
455 decaylist=pionList.plus()* pionList.minus()* pi0List* pi0List;
456 }
457 else if(channel=="D0toKKPiPi"){
458 numchan.push_back( EvtRecDTag::kD0toKKPiPi);
459 numchan.push_back(1);
460 numchan.push_back(1);
461 numchan.push_back(2);
462 numchan.push_back(2);
463 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus();
464 }
465 else if(channel=="D0toKKPi0Pi0"){
466 numchan.push_back( EvtRecDTag::kD0toKKPi0Pi0);
467 numchan.push_back(1);
468 numchan.push_back(1);
469 numchan.push_back(3);
470 numchan.push_back(3);
471 decaylist=kaonList.minus()* kaonList.plus()* pi0List* pi0List;
472 }
473 else if(channel=="D0toKsKsPiPi"){
474 numchan.push_back( EvtRecDTag::kD0toKsKsPiPi);
475 numchan.push_back(5);
476 numchan.push_back(5);
477 numchan.push_back(2);
478 numchan.push_back(2);
479 decaylist=ksList* ksList* pionList.plus()* pionList.minus();
480 }
481 else if(channel=="D0toPiPiPiPiPi0"){
482 numchan.push_back( EvtRecDTag::kD0toPiPiPiPiPi0);
483 numchan.push_back(2);
484 numchan.push_back(2);
485 numchan.push_back(2);
486 numchan.push_back(2);
487 numchan.push_back(3);
488 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
489 }
490 else if(channel=="D0toKsPiPiPiPi"){
491 numchan.push_back( EvtRecDTag::kD0toKsPiPiPiPi);
492 numchan.push_back(5);
493 numchan.push_back(2);
494 numchan.push_back(2);
495 numchan.push_back(2);
496 numchan.push_back(2);
497 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
498 }
499 else if(channel=="D0toKKPiPiPi0"){
500 numchan.push_back( EvtRecDTag::kD0toKKPiPiPi0);
501 numchan.push_back(1);
502 numchan.push_back(1);
503 numchan.push_back(2);
504 numchan.push_back(2);
505 numchan.push_back(3);
506 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
507 }
508 else if(channel=="D0toKsPi0Eta"){
509 numchan.push_back( EvtRecDTag::kD0toKsPi0Eta);
510 numchan.push_back(5);
511 numchan.push_back(3);
512 numchan.push_back(4);
513 decaylist=ksList* pi0List* etaList;
514 }
515 else if(channel=="D0toKsEPPiPiEta"){
516 numchan.push_back( EvtRecDTag::kD0toKsEPPiPiEta);
517 numchan.push_back(5);
518 numchan.push_back(6);
519 CDDecayList epList(eptoPiPiEtaSelector);
520 epList=pionList.plus()* pionList.minus()* etaList;
521 decaylist= ksList* epList;
522 }
523 else if(channel=="D0toKsEPRhoGam"){
524 numchan.push_back( EvtRecDTag::kD0toKsEPRhoGam);
525 numchan.push_back(5);
526 numchan.push_back(7);
527 CDDecayList rhoList(rhotoPiPiSelector);
528 rhoList=pionList.plus()* pionList.minus();
529 CDDecayList epList(eptoRhoGamSelector);
530 epList=rhoList* photonList;
531 decaylist= ksList* epList;
532 }
533 else if(channel=="D0toKsPi0Pi0Pi0"){
534 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0Pi0);
535 numchan.push_back(5);
536 numchan.push_back(3);
537 numchan.push_back(3);
538 numchan.push_back(3);
539 decaylist=ksList* pi0List* pi0List* pi0List;
540 }
541
542
543 CDDecayList::iterator D_begin =decaylist.particle_begin();
544 CDDecayList::iterator D_end =decaylist.particle_end();
545
546 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
547
548 EvtRecDTag* recDTag = new EvtRecDTag;
549 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
550
551
552 vector<int> trackid, showerid;
553 vector<int> kaonid, pionid;
554 int charm=0;
555 int numofchildren=numchan.size()-1;
556
557 for(int i=0; i< numofchildren;i++){
558
559 const CDCandidate& daughter=(*it).particle().child(i);
560 if(isFlavorMode && i==0)
561 charm=-daughter.charge();
562
563 if(numchan[i+1]==1){
564 const EvtRecTrack* track=daughter.track();
565 trackid.push_back(track->trackId());
566 kaonid.push_back(track->trackId());
567 }
568 else if(numchan[i+1]==2){
569 const EvtRecTrack* track=daughter.track();
570 trackid.push_back(track->trackId());
571 pionid.push_back(track->trackId());
572 }
573 else if ( numchan[i+1]==3){
574 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
575 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
576 showerid.push_back(hiEnGamma->trackId());
577 showerid.push_back(loEnGamma->trackId());
578 }
579 else if ( numchan[i+1]==4){
580 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
581 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
582 showerid.push_back(hiEnGamma->trackId());
583 showerid.push_back(loEnGamma->trackId());
584 }
585 else if ( numchan[i+1]==5){
586 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
587 // fill fit info
588 recDTag->addToFitInfo(aKsCand->mass(),fitinfo[aKsCand][0],fitinfo[aKsCand][1],fitinfo[aKsCand][2]);
589 // fill tracks
590 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
591 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
592 trackid.push_back(pion1Trk->trackId());
593 trackid.push_back(pion2Trk->trackId());
594 }
595 else if (numchan[i+1]==6){
596 const CDCandidate& apion = daughter.decay().child(0);
597 const CDCandidate& spion = daughter.decay().child(1);
598 const CDCandidate& eta = daughter.decay().child(2);
599 const EvtRecTrack* apiontrk = apion.track();
600 const EvtRecTrack* spiontrk = spion.track();
601 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
602 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
603
604 trackid.push_back(apiontrk->trackId());
605 trackid.push_back(spiontrk->trackId());
606 showerid.push_back(hiEnGamma->trackId());
607 showerid.push_back(loEnGamma->trackId());
608
609 }
610 else if (numchan[i+1]==7){
611 const CDCandidate& rho = daughter.decay().child(0);
612 const CDCandidate& gamma = daughter.decay().child(1);
613 const CDCandidate& apion = rho.decay().child(0);
614 const CDCandidate& spion = rho.decay().child(1);
615
616
617 const EvtRecTrack* apiontrk = apion.track();
618 const EvtRecTrack* spiontrk = spion.track();
619 const EvtRecTrack* gammatrk = gamma.photon();
620
621
622 trackid.push_back(apiontrk->trackId());
623 trackid.push_back(spiontrk->trackId());
624 showerid.push_back(gammatrk->trackId());
625
626 }
627
628
629 }//end of filling track and shower ids
630
631
632 saveD0Info(it, ebeam, charm, numofchildren, recDTag);
633 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
634 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
635
636 if(m_usevertexfit){
637 if(m_debug) cout<<"beforevfit:"<<endl;
638
639 HepLorentzVector p4change_vfit;
640
641 if(m_useVFrefine){
642 p4change_vfit=util.vfitref(channel, kaonid, pionid, xorigin, charged_begin);
643 }else{
644 p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
645 }
646
647 recDTag->setp4(recDTag->p4()+p4change_vfit);
648 }
649
650
651 trackid.clear();
652 showerid.clear();
653 kaonid.clear();
654 pionid.clear();
655
656
657 //write recdtag out
658 recDTagCol->push_back(recDTag);
659
660 } //end of decay list iterator
661
662 numchan.clear();
663
664 } //end of reconstrucing all D0 decay lists
665
666
667 return StatusCode::SUCCESS;
668}
669
670void NeutralDReconstruction::saveD0Info(CDDecayList::iterator it, double ebeam, int charm, int numofchildren, EvtRecDTag* recDTag){
671
672 double mass = (*it).particle().mass();
673 HepLorentzVector p4((*it).particle().momentum(), (*it).particle().energy());
674 recDTag->setp4(p4);
675
676 p4.boost(-m_beta);
677 double mbc2_CMS = ebeam*ebeam - p4.v().mag2();
678 double mbc_CMS = mbc2_CMS > 0 ? sqrt( mbc2_CMS ) : -10;
679 double deltaE_CMS = p4.t() - ebeam;
680
681 int tag=(*it).particle().userTag();
682 if(tag==1)
683 recDTag->settype( EvtRecDTag::Tight );
684 else
685 recDTag->settype( EvtRecDTag::Loose );
686 recDTag->setcharge(0);
687 recDTag->setcharm(charm);
688 recDTag->setnumOfChildren(numofchildren);
689 recDTag->setmass(mass);
690 recDTag->setmBC(mbc_CMS);
691 recDTag->setbeamE(ebeam);
692 recDTag->setdeltaE(deltaE_CMS);
693
694}
695
696void NeutralDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
697 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
698
699 vector<EvtRecTrackIterator> trktemp;
700 vector<EvtRecTrackIterator> shrtemp;
701
702 //fill tracks
703 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
704
705 bool isothertrack=true;
706 for(int i=0; i<trackid.size(); i++){
707 if( (*trk)->trackId()==trackid[i] ){
708 trktemp.push_back(trk);
709 isothertrack=false;
710 break;
711 }
712 }
713 if(isothertrack)
714 recDTag->addOtherTrack(*trk);
715 }
716 for(int i=0; i<trackid.size();i++){
717 for(int j=0; j<trktemp.size(); j++){
718 EvtRecTrackIterator trk=trktemp[j];
719 if( (*trk)->trackId()==trackid[i])
720 recDTag->addTrack(*trktemp[j]);
721 }
722 }
723
724
725 //fill showers
726 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
727 bool isothershower=true;
728 for(int i=0; i<showerid.size(); i++){
729 if( (*shr)->trackId()==showerid[i] ){
730 shrtemp.push_back(shr);
731 isothershower=false;
732 break;
733 }
734 }
735 if(isothershower)
736 recDTag->addOtherShower(*shr);
737 }
738
739 for(int i=0; i<showerid.size();i++){
740 for(int j=0; j<shrtemp.size(); j++){
741 EvtRecTrackIterator shr=shrtemp[j];
742 if( (*shr)->trackId()==showerid[i])
743 recDTag->addShower(*shrtemp[j]);
744 }
745 }
746
747
748}
749
750
751void NeutralDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
752
753 bool iskaon=false,ispion=false;
754
755
756 // save track ids which passed pion/kaon cuts
757
758 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
759 recDTag->addKaonId( (*kit).particle().track() );
760 }
761
762 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
763 recDTag->addPionId( (*pit).particle().track() );
764 }
765
766
767
768 /*
769 for(int i=0; i<kaonid.size(); i++){
770 bool ithkaon=false;
771 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
772 if((*kit).particle().track()->trackId()==kaonid[i]){
773 ithkaon=true;
774 break;
775 }
776 }
777 if(!ithkaon) break;
778 if(i==kaonid.size()-1)
779 iskaon=true;
780 }
781
782 for(int i=0; i<pionid.size(); i++){
783 bool ithpion=false;
784 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
785 if((*pit).particle().track()->trackId()==pionid[i]){
786 ithpion=true;
787 break;
788 }
789 }
790 if(!ithpion) break;
791 if(i==pionid.size()-1)
792 ispion=true;
793 }
794
795
796 if( iskaon && ispion)
797 recDTag->settype( EvtRecDTag::Tight );
798 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
799 recDTag->settype( EvtRecDTag::Tight );
800 else if( kaonid.size()==0 && pionid.size()==0 )
801 recDTag->settype( EvtRecDTag::Tight );
802 */
803
804}
805
806
807
808vector<string> NeutralDReconstruction::getlist(string& filename ){
809
810 string channel;
811 vector<string> temp;
812
813 ifstream inFile;
814
815 inFile.open(filename.c_str());
816 if (!inFile) {
817 cout << "Unable to open decay list file";
818 exit(1); // terminate with error
819 }
820
821 while (inFile >> channel) {
822 temp.push_back(channel);
823 }
824
825 inFile.close();
826
827 return temp;
828
829}
mass
double mass
Definition CosmicGenerator.cxx:138
dc_fill
void dc_fill(DCFillableChargedList< Charged > &aFillableList, WitnessIterator first, WitnessIterator last)
Definition DCFillableChargedList.h:152
ebeam
double ebeam
Definition EvtPhokharaDef.hh:17
EvtRecDTagCol
ObjectVector< EvtRecDTag > EvtRecDTagCol
Definition EvtRecDTag.h:286
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition EvtRecTrack.h:111
eptoPiPiEtaSelector
LocalEptoPiPiEtaSelector eptoPiPiEtaSelector
Definition LocalEptoPiPiEtaSelector.cxx:45
eptoRhoGamSelector
LocalEptoRhoGamSelector eptoRhoGamSelector
Definition LocalEptoRhoGamSelector.cxx:40
etatoGGSelector
LocalEtatoGGSelector etatoGGSelector
Definition LocalEtatoGGSelector.cxx:37
kaonSelector
LocalKaonSelector kaonSelector
Definition LocalKaonSelector.cxx:135
ksSelector
LocalKsSelector ksSelector
Definition LocalKsSelector.cxx:207
photonSelector
LocalPhotonSelector photonSelector
Definition LocalPhotonSelector.cxx:115
pi0Selector
LocalPi0Selector pi0Selector
Definition LocalPi0Selector.cxx:37
pionSelector
LocalPionSelector pionSelector
Definition LocalPionSelector.cxx:134
rhotoPiPiSelector
LocalRhotoPiPiSelector rhotoPiPiSelector
Definition LocalRhotoPiPiSelector.cxx:42
neutralDSelector
NeutralDSelector neutralDSelector
Definition NeutralDSelector.cxx:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
CDCandidate::decay
virtual const DecayEvidence & decay() const
Definition CDCandidate.cxx:356
CDCandidate::photon
virtual const EvtRecTrack * photon() const
Definition CDCandidate.cxx:256
CDCandidate::track
virtual const EvtRecTrack * track() const
Definition CDCandidate.cxx:236
CDCandidate::charge
int charge() const
Definition CDCandidate.cxx:154
CDCandidate::navKshort
virtual const EvtRecVeeVertex * navKshort() const
Definition CDCandidate.cxx:316
CDCandidate::navPi0
virtual const EvtRecPi0 * navPi0() const
Definition CDCandidate.cxx:276
CDCandidate::navEta
virtual const EvtRecEtaToGG * navEta() const
Definition CDCandidate.cxx:296
CDDecay::child
const CDCandidate & child(unsigned int aPosition) const
Definition CDDecay.cxx:247
EvtRecDTag::kD0toKsEPRhoGam
@ kD0toKsEPRhoGam
Definition EvtRecDTag.h:91
EvtRecDTag::kD0toPiPiPiPiPi0
@ kD0toPiPiPiPiPi0
Definition EvtRecDTag.h:85
EvtRecDTag::kD0toKsPiPiPiPi
@ kD0toKsPiPiPiPi
Definition EvtRecDTag.h:86
EvtRecDTag::kD0toKsPi0Pi0Pi0
@ kD0toKsPi0Pi0Pi0
Definition EvtRecDTag.h:92
EvtRecDTag::kD0toPiPiPi0Pi0
@ kD0toPiPiPi0Pi0
Definition EvtRecDTag.h:81
EvtRecDTag::kD0toKsEPPiPiEta
@ kD0toKsEPPiPiEta
Definition EvtRecDTag.h:90
EvtRecDTag::kD0toKPiPiPiPi0
@ kD0toKPiPiPiPi0
Definition EvtRecDTag.h:57
EvtRecDTag::addOtherTrack
void addOtherTrack(const SmartRef< EvtRecTrack > track)
Definition EvtRecDTag.h:232
EvtRecDTag::settype
void settype(SelType type)
Definition EvtRecDTag.h:211
EvtRecDTag::setdecayMode
void setdecayMode(DecayMode decayMode)
Definition EvtRecDTag.h:210
EvtRecDTag::p4
HepLorentzVector p4() const
Definition EvtRecDTag.h:194
EvtRecDTag::setp4
void setp4(HepLorentzVector p4)
Definition EvtRecDTag.h:219
EvtRecDTag::setmass
void setmass(double mass)
Definition EvtRecDTag.h:213
EvtRecDTag::addToFitInfo
void addToFitInfo(double ksmass, double chi2, double length, double error)
Definition EvtRecDTag.h:221
EvtRecDTag::addOtherShower
void addOtherShower(const SmartRef< EvtRecTrack > shower)
Definition EvtRecDTag.h:234
EvtRecDTag::addKaonId
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
Definition EvtRecDTag.h:238
EvtRecDTag::setdeltaE
void setdeltaE(double deltaE)
Definition EvtRecDTag.h:215
EvtRecDTag::setcharm
void setcharm(int charm)
Definition EvtRecDTag.h:217
EvtRecDTag::setmBC
void setmBC(double mBC)
Definition EvtRecDTag.h:214
EvtRecDTag::setcharge
void setcharge(int charge)
Definition EvtRecDTag.h:216
EvtRecDTag::addPionId
void addPionId(const SmartRef< EvtRecTrack > pionId)
Definition EvtRecDTag.h:236
EvtRecDTag::setbeamE
void setbeamE(double beamE)
Definition EvtRecDTag.h:212
EvtRecDTag::addShower
void addShower(const SmartRef< EvtRecTrack > shower)
Definition EvtRecDTag.h:230
EvtRecDTag::setnumOfChildren
void setnumOfChildren(int numOfChildren)
Definition EvtRecDTag.h:218
EvtRecDTag::addTrack
void addTrack(const SmartRef< EvtRecTrack > track)
Definition EvtRecDTag.h:228
EvtRecEtaToGG::hiEnGamma
const EvtRecTrack * hiEnGamma() const
Definition EvtRecEtaToGG.h:30
EvtRecEtaToGG::loEnGamma
const EvtRecTrack * loEnGamma() const
Definition EvtRecEtaToGG.h:31
EvtRecPi0::loEnGamma
const EvtRecTrack * loEnGamma() const
Definition EvtRecPi0.h:31
EvtRecPi0::hiEnGamma
const EvtRecTrack * hiEnGamma() const
Definition EvtRecPi0.h:30
EvtRecTrack::trackId
int trackId() const
Definition EvtRecTrack.h:32
EvtRecVeeVertex::daughter
SmartRef< EvtRecTrack > & daughter(int i)
Definition EvtRecVeeVertex.h:89
EvtRecVeeVertex::mass
double mass() const
Definition EvtRecVeeVertex.h:39
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
LocalKaonSelector::setpidtype
void setpidtype(int type)
Definition LocalKaonSelector.h:15
LocalPionSelector::setpidtype
void setpidtype(int type)
Definition LocalPionSelector.h:14
NeutralDReconstruction::savetrack
void savetrack(vector< int >, vector< int >, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecTrackIterator, EvtRecDTag *)
Definition NeutralDReconstruction.cxx:696
NeutralDReconstruction::getlist
vector< string > getlist(string &filename)
Definition NeutralDReconstruction.cxx:808
NeutralDReconstruction::pidtag
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
Definition NeutralDReconstruction.cxx:751
NeutralDReconstruction::saveD0Info
void saveD0Info(CDDecayList::iterator, double, int, int, EvtRecDTag *)
Definition NeutralDReconstruction.cxx:670
NeutralDSelector::setbeta
void setbeta(Hep3Vector beta)
Definition NeutralDSelector.h:14
NeutralDSelector::setebeam
void setebeam(double ebeam)
Definition NeutralDSelector.h:13
dchain::CandidateList::size
virtual int size() const
Definition CandidateList.cc:143
dchain::ChargedVisibleList::plus
ChosenChargeList< Charged, CandidateClass > & plus() const
Definition ChargedVisibleList.cc:176
dchain::ChargedVisibleList::minus
ChosenChargeList< Charged, CandidateClass > & minus() const
Definition ChargedVisibleList.cc:186
dchain::DecayList::particle_begin
virtual iterator particle_begin()
Definition DecayList.cc:169
dchain::DecayList::particle_end
virtual iterator particle_end()
Definition DecayList.cc:176
utility::vfitref
HepLorentzVector vfitref(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:59
utility::SecondaryVFitref
vector< double > SecondaryVFitref(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:199
utility::vfit
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition utility.cxx:406
utility::SecondaryVFit
vector< double > SecondaryVFit(EvtRecVeeVertex *ks, IVertexDbSvc *vtxsvc)
Definition utility.cxx:540
m_debug
bool m_debug
Definition draw_charge_space_coarse_II.cxx:3
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition EventModel.h:116
EventModel::EvtRec::EvtRecDTagCol
_EXTERN_ std::string EvtRecDTagCol
Definition EventModel.h:122
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition EventModel.h:117
Event
Definition Event.h:21
p4
double double double * p4
Definition qcdloop1.h:77