CGEM BOSS 6.6.5.f
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 CDKsList ksList(ksSelector);
199 dc_fill(ksList, recVeeVertexCol->begin(), recVeeVertexCol->end());
200
201 CDPi0List pi0List(pi0Selector);
202 dc_fill(pi0List, recPi0Col->begin(), recPi0Col->end());
203
204 CDEtaList etaList(etatoGGSelector);
205 dc_fill(etaList, recEtaToGGCol->begin(), recEtaToGGCol->end());
206
207
208 //pion/kaon list with PID
209 pionSelector.setpidtype(1);
210 kaonSelector.setpidtype(1);
211 CDChargedPionList pionList_tight(charged_begin, charged_end, pionSelector);
212 CDChargedKaonList kaonList_tight(charged_begin, charged_end, kaonSelector);
213
214 int run = eventHeader->runNumber();
215 m_ievt = eventHeader->eventNumber();
216 m_nChrg = recEvent->totalCharged();
217 m_nNeu = recEvent->totalNeutral();
218 m_nPion = pionList.size();
219 m_nKaon = kaonList.size();
220 m_nPi0 = pi0List.size();
221 m_nKs = ksList.size();
222
223 ///////////////////////
224 // get beam energy and beta
225 ///////////////////////
226
227 if(m_ReadBeamEFromDB && m_irun!=run){
228 m_irun=run;
229 if(m_usecalibBeamE)
230 m_readDb.setcalib(true);
231 m_beamE=m_readDb.getbeamE(m_irun, m_beamE);
232 if(run>0)
233 m_beta=m_readDb.getbeta();
234
235 //cout<<"use beam E from data base:"<<m_beamE<<endl;
236 }
237 double ebeam=m_beamE;
238
239 //////////////////////////////
240 //reconstruct decay lists
241 /////////////////////////////
242
243
244
245 for(int list=0;list<chanlist.size();list++){
246
247 string channel=chanlist[list];
248 vector<int> numchan;
249 neutralDSelector.setebeam(ebeam);
250 neutralDSelector.setbeta(m_beta);
251 CDDecayList decaylist(neutralDSelector);
252
253 bool isFlavorMode=false;
254
255 //K+/-: 1, Pi+/-:2, Pi0:3, Eta: 4, Ks:5, Eta'(pipiEta):6, Eta'(rhogam):7
256 //the fist element of the vector stands for decay mode,
257 //the rest will be particles, and size of the vector minus 1 will be number of daughers.
258 //the first particle in the vector will be used to get charm
259
260 if(channel=="D0toKPi") {
261 numchan.push_back( EvtRecDTag::kD0toKPi );
262 numchan.push_back(1);
263 numchan.push_back(2);
264 decaylist=kaonList.minus() * pionList.plus();
265 isFlavorMode=true;
266 }
267 else if(channel=="D0toKPiPi0"){
268 numchan.push_back( EvtRecDTag::kD0toKPiPi0);
269 numchan.push_back(1);
270 numchan.push_back(2);
271 numchan.push_back(3);
272 decaylist=kaonList.minus() * pionList.plus() * pi0List;
273 isFlavorMode=true;
274 }
275 else if(channel=="D0toKPiPi0Pi0"){
276 numchan.push_back( EvtRecDTag::kD0toKPiPi0Pi0);
277 numchan.push_back(1);
278 numchan.push_back(2);
279 numchan.push_back(3);
280 numchan.push_back(3);
281 decaylist=kaonList.minus() * pionList.plus() * pi0List * pi0List;
282 isFlavorMode=true;
283 }
284 else if(channel=="D0toKPiPiPi"){
285 numchan.push_back( EvtRecDTag::kD0toKPiPiPi);
286 numchan.push_back(1);
287 numchan.push_back(2);
288 numchan.push_back(2);
289 numchan.push_back(2);
290 decaylist=kaonList.minus()* pionList.plus()* pionList.plus() * pionList.minus();
291 isFlavorMode=true;
292 }
293 else if(channel=="D0toKPiPiPiPi0"){
294 numchan.push_back( EvtRecDTag::kD0toKPiPiPiPi0);
295 numchan.push_back(1);
296 numchan.push_back(2);
297 numchan.push_back(2);
298 numchan.push_back(2);
299 numchan.push_back(3);
300 decaylist=kaonList.minus()* pionList.plus()* pionList.plus()* pionList.minus()* pi0List;
301 isFlavorMode=true;
302 }
303 else if(channel=="D0toKPiEta"){
304 numchan.push_back( EvtRecDTag::kD0toKPiEta);
305 numchan.push_back(1);
306 numchan.push_back(2);
307 numchan.push_back(4);
308 decaylist=kaonList.minus() * pionList.plus() * etaList;
309 isFlavorMode=true;
310 }
311 else if(channel=="D0toKsKPi"){
312 numchan.push_back( EvtRecDTag::kD0toKsKPi);
313 numchan.push_back(5);
314 numchan.push_back(1);
315 numchan.push_back(2);
316 decaylist=ksList* kaonList.minus()* pionList.plus();
317 }
318 else if(channel=="D0toKsKPiPi0"){
319 numchan.push_back( EvtRecDTag::kD0toKsKPiPi0);
320 numchan.push_back(5);
321 numchan.push_back(1);
322 numchan.push_back(2);
323 numchan.push_back(3);
324 decaylist=ksList* kaonList.minus() * pionList.plus()* pi0List;
325 }
326 else if(channel=="D0toKsPiPi"){
327 numchan.push_back( EvtRecDTag::kD0toKsPiPi);
328 numchan.push_back(5);
329 numchan.push_back(2);
330 numchan.push_back(2);
331 decaylist=ksList* pionList.plus()* pionList.minus();
332 }
333 else if(channel=="D0toKsPiPiPi0"){
334 numchan.push_back( EvtRecDTag::kD0toKsPiPiPi0);
335 numchan.push_back(5);
336 numchan.push_back(2);
337 numchan.push_back(2);
338 numchan.push_back(3);
339 decaylist=ksList* pionList.plus()* pionList.minus()* pi0List;
340 }
341 else if(channel=="D0toKsPi0"){
342 numchan.push_back( EvtRecDTag::kD0toKsPi0);
343 numchan.push_back(5);
344 numchan.push_back(3);
345 decaylist=ksList* pi0List;
346 }
347 else if(channel=="D0toPiPiPi0"){
348 numchan.push_back( EvtRecDTag::kD0toPiPiPi0);
349 numchan.push_back(2);
350 numchan.push_back(2);
351 numchan.push_back(3);
352 decaylist=pionList.plus()* pionList.minus()* pi0List;
353 }
354 else if(channel=="D0toPiPi"){
355 numchan.push_back( EvtRecDTag::kD0toPiPi);
356 numchan.push_back(2);
357 numchan.push_back(2);
358 decaylist=pionList.plus()* pionList.minus();
359 }
360 else if(channel=="D0toKK"){
361 numchan.push_back( EvtRecDTag::kD0toKK);
362 numchan.push_back(1);
363 numchan.push_back(1);
364 decaylist=kaonList.minus()* kaonList.plus();
365 }
366 else if(channel=="D0toKKPi0"){
367 numchan.push_back( EvtRecDTag::kD0toKKPi0);
368 numchan.push_back(1);
369 numchan.push_back(1);
370 numchan.push_back(3);
371 decaylist=kaonList.minus()* kaonList.plus()* pi0List;
372 }
373 else if(channel=="D0toPi0Pi0"){
374 numchan.push_back( EvtRecDTag::kD0toPi0Pi0);
375 numchan.push_back(3);
376 numchan.push_back(3);
377 decaylist=pi0List* pi0List;
378 }
379 else if(channel=="D0toKsKs"){
380 numchan.push_back( EvtRecDTag::kD0toKsKs);
381 numchan.push_back(5);
382 numchan.push_back(5);
383 decaylist=ksList* ksList;
384 }
385 else if(channel=="D0toKsKsPi0"){
386 numchan.push_back( EvtRecDTag::kD0toKsKsPi0);
387 numchan.push_back(5);
388 numchan.push_back(5);
389 numchan.push_back(3);
390 decaylist=ksList* ksList* pi0List;
391 }
392 else if(channel=="D0toKsPi0Pi0"){
393 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0);
394 numchan.push_back(5);
395 numchan.push_back(3);
396 numchan.push_back(3);
397 decaylist=ksList* pi0List* pi0List;
398 }
399 else if(channel=="D0toKsKK"){
400 numchan.push_back( EvtRecDTag::kD0toKsKK);
401 numchan.push_back(5);
402 numchan.push_back(1);
403 numchan.push_back(1);
404 decaylist=ksList* kaonList.minus()* kaonList.plus();
405 }
406 else if(channel=="D0toKsEta"){
407 numchan.push_back( EvtRecDTag::kD0toKsEta);
408 numchan.push_back(5);
409 numchan.push_back(4);
410 decaylist=ksList* etaList;
411 }
412 else if(channel=="D0toPi0Pi0Pi0"){
413 numchan.push_back( EvtRecDTag::kD0toPi0Pi0Pi0);
414 numchan.push_back(3);
415 numchan.push_back(3);
416 numchan.push_back(3);
417 decaylist=pi0List* pi0List* pi0List;
418 }
419 else if(channel=="D0toKsKsKs"){
420 numchan.push_back( EvtRecDTag::kD0toKsKsKs);
421 numchan.push_back(5);
422 numchan.push_back(5);
423 numchan.push_back(5);
424 decaylist=ksList* ksList* ksList;
425 }
426 else if(channel=="D0toPiPiPiPi"){
427 numchan.push_back( EvtRecDTag::kD0toPiPiPiPi);
428 numchan.push_back(2);
429 numchan.push_back(2);
430 numchan.push_back(2);
431 numchan.push_back(2);
432 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
433 }
434 else if(channel=="D0toPiPiPi0Pi0"){
435 numchan.push_back( EvtRecDTag::kD0toPiPiPi0Pi0);
436 numchan.push_back(2);
437 numchan.push_back(2);
438 numchan.push_back(3);
439 numchan.push_back(3);
440 decaylist=pionList.plus()* pionList.minus()* pi0List* pi0List;
441 }
442 else if(channel=="D0toKKPiPi"){
443 numchan.push_back( EvtRecDTag::kD0toKKPiPi);
444 numchan.push_back(1);
445 numchan.push_back(1);
446 numchan.push_back(2);
447 numchan.push_back(2);
448 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus();
449 }
450 else if(channel=="D0toKKPi0Pi0"){
451 numchan.push_back( EvtRecDTag::kD0toKKPi0Pi0);
452 numchan.push_back(1);
453 numchan.push_back(1);
454 numchan.push_back(3);
455 numchan.push_back(3);
456 decaylist=kaonList.minus()* kaonList.plus()* pi0List* pi0List;
457 }
458 else if(channel=="D0toKsKsPiPi"){
459 numchan.push_back( EvtRecDTag::kD0toKsKsPiPi);
460 numchan.push_back(5);
461 numchan.push_back(5);
462 numchan.push_back(2);
463 numchan.push_back(2);
464 decaylist=ksList* ksList* pionList.plus()* pionList.minus();
465 }
466 else if(channel=="D0toPiPiPiPiPi0"){
467 numchan.push_back( EvtRecDTag::kD0toPiPiPiPiPi0);
468 numchan.push_back(2);
469 numchan.push_back(2);
470 numchan.push_back(2);
471 numchan.push_back(2);
472 numchan.push_back(3);
473 decaylist=pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus()* pi0List;
474 }
475 else if(channel=="D0toKsPiPiPiPi"){
476 numchan.push_back( EvtRecDTag::kD0toKsPiPiPiPi);
477 numchan.push_back(5);
478 numchan.push_back(2);
479 numchan.push_back(2);
480 numchan.push_back(2);
481 numchan.push_back(2);
482 decaylist=ksList* pionList.plus()* pionList.plus()* pionList.minus()* pionList.minus();
483 }
484 else if(channel=="D0toKKPiPiPi0"){
485 numchan.push_back( EvtRecDTag::kD0toKKPiPiPi0);
486 numchan.push_back(1);
487 numchan.push_back(1);
488 numchan.push_back(2);
489 numchan.push_back(2);
490 numchan.push_back(3);
491 decaylist=kaonList.minus()* kaonList.plus()* pionList.plus()* pionList.minus()* pi0List;
492 }
493 else if(channel=="D0toKsPi0Eta"){
494 numchan.push_back( EvtRecDTag::kD0toKsPi0Eta);
495 numchan.push_back(5);
496 numchan.push_back(3);
497 numchan.push_back(4);
498 decaylist=ksList* pi0List* etaList;
499 }
500 else if(channel=="D0toKsEPPiPiEta"){
501 numchan.push_back( EvtRecDTag::kD0toKsEPPiPiEta);
502 numchan.push_back(5);
503 numchan.push_back(6);
504 CDDecayList epList(eptoPiPiEtaSelector);
505 epList=pionList.plus()* pionList.minus()* etaList;
506 decaylist= ksList* epList;
507 }
508 else if(channel=="D0toKsEPRhoGam"){
509 numchan.push_back( EvtRecDTag::kD0toKsEPRhoGam);
510 numchan.push_back(5);
511 numchan.push_back(7);
512 CDDecayList rhoList(rhotoPiPiSelector);
513 rhoList=pionList.plus()* pionList.minus();
514 CDDecayList epList(eptoRhoGamSelector);
515 epList=rhoList* photonList;
516 decaylist= ksList* epList;
517 }
518 else if(channel=="D0toKsPi0Pi0Pi0"){
519 numchan.push_back( EvtRecDTag::kD0toKsPi0Pi0Pi0);
520 numchan.push_back(5);
521 numchan.push_back(3);
522 numchan.push_back(3);
523 numchan.push_back(3);
524 decaylist=ksList* pi0List* pi0List* pi0List;
525 }
526
527
528 CDDecayList::iterator D_begin =decaylist.particle_begin();
529 CDDecayList::iterator D_end =decaylist.particle_end();
530
531 for ( CDDecayList::iterator it = D_begin; it != D_end; it++ ) {
532
533 EvtRecDTag* recDTag = new EvtRecDTag;
534 recDTag->setdecayMode( (EvtRecDTag::DecayMode)numchan[0] );
535
536
537 vector<int> trackid, showerid;
538 vector<int> kaonid, pionid;
539 int charm=0;
540 int numofchildren=numchan.size()-1;
541
542 for(int i=0; i< numofchildren;i++){
543
544 const CDCandidate& daughter=(*it).particle().child(i);
545 if(isFlavorMode && i==0)
546 charm=-daughter.charge();
547
548 if(numchan[i+1]==1){
549 const EvtRecTrack* track=daughter.track();
550 trackid.push_back(track->trackId());
551 kaonid.push_back(track->trackId());
552 }
553 else if(numchan[i+1]==2){
554 const EvtRecTrack* track=daughter.track();
555 trackid.push_back(track->trackId());
556 pionid.push_back(track->trackId());
557 }
558 else if ( numchan[i+1]==3){
559 const EvtRecTrack* hiEnGamma=daughter.navPi0()->hiEnGamma();
560 const EvtRecTrack* loEnGamma=daughter.navPi0()->loEnGamma();
561 showerid.push_back(hiEnGamma->trackId());
562 showerid.push_back(loEnGamma->trackId());
563 }
564 else if ( numchan[i+1]==4){
565 const EvtRecTrack* hiEnGamma=daughter.navEta()->hiEnGamma();
566 const EvtRecTrack* loEnGamma=daughter.navEta()->loEnGamma();
567 showerid.push_back(hiEnGamma->trackId());
568 showerid.push_back(loEnGamma->trackId());
569 }
570 else if ( numchan[i+1]==5){
571 EvtRecVeeVertex* aKsCand = const_cast<EvtRecVeeVertex*>( daughter.navKshort() );
572 EvtRecTrack* pion1Trk = aKsCand->daughter(0);
573 EvtRecTrack* pion2Trk = aKsCand->daughter(1);
574 trackid.push_back(pion1Trk->trackId());
575 trackid.push_back(pion2Trk->trackId());
576 }
577 else if (numchan[i+1]==6){
578 const CDCandidate& apion = daughter.decay().child(0);
579 const CDCandidate& spion = daughter.decay().child(1);
580 const CDCandidate& eta = daughter.decay().child(2);
581 const EvtRecTrack* apiontrk = apion.track();
582 const EvtRecTrack* spiontrk = spion.track();
583 const EvtRecTrack* hiEnGamma=eta.navEta()->hiEnGamma();
584 const EvtRecTrack* loEnGamma=eta.navEta()->loEnGamma();
585
586 trackid.push_back(apiontrk->trackId());
587 trackid.push_back(spiontrk->trackId());
588 showerid.push_back(hiEnGamma->trackId());
589 showerid.push_back(loEnGamma->trackId());
590
591 }
592 else if (numchan[i+1]==7){
593 const CDCandidate& rho = daughter.decay().child(0);
594 const CDCandidate& gamma = daughter.decay().child(1);
595 const CDCandidate& apion = rho.decay().child(0);
596 const CDCandidate& spion = rho.decay().child(1);
597
598
599 const EvtRecTrack* apiontrk = apion.track();
600 const EvtRecTrack* spiontrk = spion.track();
601 const EvtRecTrack* gammatrk = gamma.photon();
602
603
604 trackid.push_back(apiontrk->trackId());
605 trackid.push_back(spiontrk->trackId());
606 showerid.push_back(gammatrk->trackId());
607
608 }
609
610
611 }//end of filling track and shower ids
612
613
614 saveD0Info(it, ebeam, charm, numofchildren, recDTag);
615 savetrack(trackid,showerid,charged_begin,charged_end,neutral_begin,neutral_end,recDTag);
616 pidtag(kaonid,pionid,kaonList_tight, pionList_tight,recDTag);
617
618 if(m_usevertexfit){
619 HepLorentzVector p4change_vfit=util.vfit(channel, kaonid, pionid, xorigin, charged_begin);
620 recDTag->setp4(recDTag->p4()+p4change_vfit);
621 }
622
623
624 trackid.clear();
625 showerid.clear();
626 kaonid.clear();
627 pionid.clear();
628
629
630 //write recdtag out
631 recDTagCol->push_back(recDTag);
632
633 } //end of decay list iterator
634
635 numchan.clear();
636
637 } //end of reconstrucing all D0 decay lists
638
639
640 return StatusCode::SUCCESS;
641}
642
643void NeutralDReconstruction::saveD0Info(CDDecayList::iterator it, double ebeam, int charm, int numofchildren, EvtRecDTag* recDTag){
644
645 double mass = (*it).particle().mass();
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 int tag=(*it).particle().userTag();
655 if(tag==1)
656 recDTag->settype( EvtRecDTag::Tight );
657 else
658 recDTag->settype( EvtRecDTag::Loose );
659 recDTag->setcharge(0);
660 recDTag->setcharm(charm);
661 recDTag->setnumOfChildren(numofchildren);
662 recDTag->setmass(mass);
663 recDTag->setmBC(mbc_CMS);
664 recDTag->setbeamE(ebeam);
665 recDTag->setdeltaE(deltaE_CMS);
666
667}
668
669void NeutralDReconstruction::savetrack(vector<int> trackid, vector<int> showerid, EvtRecTrackIterator charged_begin, EvtRecTrackIterator charged_end,
670 EvtRecTrackIterator neutral_begin, EvtRecTrackIterator neutral_end,EvtRecDTag* recDTag){
671
672 vector<EvtRecTrackIterator> trktemp;
673 vector<EvtRecTrackIterator> shrtemp;
674
675 //fill tracks
676 for(EvtRecTrackIterator trk=charged_begin; trk<charged_end;trk++){
677
678 bool isothertrack=true;
679 for(int i=0; i<trackid.size(); i++){
680 if( (*trk)->trackId()==trackid[i] ){
681 trktemp.push_back(trk);
682 isothertrack=false;
683 break;
684 }
685 }
686 if(isothertrack)
687 recDTag->addOtherTrack(*trk);
688 }
689 for(int i=0; i<trackid.size();i++){
690 for(int j=0; j<trktemp.size(); j++){
691 EvtRecTrackIterator trk=trktemp[j];
692 if( (*trk)->trackId()==trackid[i])
693 recDTag->addTrack(*trktemp[j]);
694 }
695 }
696
697
698 //fill showers
699 for(EvtRecTrackIterator shr=neutral_begin; shr<neutral_end;shr++){
700 bool isothershower=true;
701 for(int i=0; i<showerid.size(); i++){
702 if( (*shr)->trackId()==showerid[i] ){
703 shrtemp.push_back(shr);
704 isothershower=false;
705 break;
706 }
707 }
708 if(isothershower)
709 recDTag->addOtherShower(*shr);
710 }
711
712 for(int i=0; i<showerid.size();i++){
713 for(int j=0; j<shrtemp.size(); j++){
714 EvtRecTrackIterator shr=shrtemp[j];
715 if( (*shr)->trackId()==showerid[i])
716 recDTag->addShower(*shrtemp[j]);
717 }
718 }
719
720
721}
722
723
724void NeutralDReconstruction::pidtag(vector<int> kaonid, vector<int> pionid, CDChargedKaonList& kaonList, CDChargedPionList& pionList,EvtRecDTag* recDTag){
725
726 bool iskaon=false,ispion=false;
727
728
729 // save track ids which passed pion/kaon cuts
730
731 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
732 recDTag->addKaonId( (*kit).particle().track() );
733 }
734
735 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
736 recDTag->addPionId( (*pit).particle().track() );
737 }
738
739
740
741 /*
742 for(int i=0; i<kaonid.size(); i++){
743 bool ithkaon=false;
744 for (CDChargedKaonList::iterator kit = kaonList.particle_begin(); kit != kaonList.particle_end(); kit++) {
745 if((*kit).particle().track()->trackId()==kaonid[i]){
746 ithkaon=true;
747 break;
748 }
749 }
750 if(!ithkaon) break;
751 if(i==kaonid.size()-1)
752 iskaon=true;
753 }
754
755 for(int i=0; i<pionid.size(); i++){
756 bool ithpion=false;
757 for (CDChargedPionList::iterator pit = pionList.particle_begin(); pit != pionList.particle_end(); pit++) {
758 if((*pit).particle().track()->trackId()==pionid[i]){
759 ithpion=true;
760 break;
761 }
762 }
763 if(!ithpion) break;
764 if(i==pionid.size()-1)
765 ispion=true;
766 }
767
768
769 if( iskaon && ispion)
770 recDTag->settype( EvtRecDTag::Tight );
771 else if( (kaonid.size()==0 && ispion) || (pionid.size()==0 && iskaon))
772 recDTag->settype( EvtRecDTag::Tight );
773 else if( kaonid.size()==0 && pionid.size()==0 )
774 recDTag->settype( EvtRecDTag::Tight );
775 */
776
777}
778
779
780
781vector<string> NeutralDReconstruction::getlist(string& filename ){
782
783 string channel;
784 vector<string> temp;
785
786 ifstream inFile;
787
788 inFile.open(filename.c_str());
789 if (!inFile) {
790 cout << "Unable to open decay list file";
791 exit(1); // terminate with error
792 }
793
794 while (inFile >> channel) {
795 temp.push_back(channel);
796 }
797
798 inFile.close();
799
800 return temp;
801
802}
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:40
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-03/EvtRecEvent/EvtRecDTag.h:187
EvtRecDTag::addOtherShower
void addOtherShower(const SmartRef< EvtRecTrack > shower)
Definition: Event/EvtRecEvent/EvtRecEvent-00-02-03/EvtRecEvent/EvtRecDTag.h:189
EvtRecDTag::addKaonId
void addKaonId(const SmartRef< EvtRecTrack > kaonId)
Definition: Event/EvtRecEvent/EvtRecEvent-00-02-03/EvtRecEvent/EvtRecDTag.h:193
EvtRecDTag::addPionId
void addPionId(const SmartRef< EvtRecTrack > pionId)
Definition: Event/EvtRecEvent/EvtRecEvent-00-02-03/EvtRecEvent/EvtRecDTag.h:191
EvtRecDTag::addShower
void addShower(const SmartRef< EvtRecTrack > shower)
Definition: Event/EvtRecEvent/EvtRecEvent-00-02-03/EvtRecEvent/EvtRecDTag.h:185
EvtRecDTag::addTrack
void addTrack(const SmartRef< EvtRecTrack > track)
Definition: Event/EvtRecEvent/EvtRecEvent-00-02-03/EvtRecEvent/EvtRecDTag.h:183
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:669
NeutralDReconstruction::getlist
vector< string > getlist(string &filename)
Definition: NeutralDReconstruction.cxx:781
NeutralDReconstruction::pidtag
void pidtag(vector< int >, vector< int >, CDChargedKaonList &, CDChargedPionList &, EvtRecDTag *)
Definition: NeutralDReconstruction.cxx:724
NeutralDReconstruction::saveD0Info
void saveD0Info(CDDecayList::iterator, double, int, int, EvtRecDTag *)
Definition: NeutralDReconstruction.cxx:643
utility::vfit
HepLorentzVector vfit(string channel, vector< int > kaonid, vector< int > pionid, HepPoint3D vx, EvtRecTrackIterator charged_begin)
Definition: utility.cxx:57