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