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