BOSS 7.1.1
BESIII Offline Software System
Loading...
Searching...
No Matches
CalibEventSelect.cxx
Go to the documentation of this file.
1 #include "GaudiKernel/MsgStream.h"
2 #include "GaudiKernel/AlgFactory.h"
3 #include "GaudiKernel/SmartDataPtr.h"
4 #include "GaudiKernel/IDataProviderSvc.h"
5 #include "GaudiKernel/PropertyMgr.h"
6
7 #include "EventModel/EventModel.h"
8 #include "EventModel/Event.h"
9 #include "EventModel/EventHeader.h"
10 #include "EvtRecEvent/EvtRecEvent.h"
11 #include "EvtRecEvent/EvtRecTrack.h"
12 #include "EvtRecEvent/EvtRecDTag.h"
13
14 #include "EvtRecEvent/EvtRecPi0.h"
15
16 #include "TMath.h"
17 #include "GaudiKernel/INTupleSvc.h"
18 #include "GaudiKernel/NTuple.h"
19 #include "GaudiKernel/Bootstrap.h"
20 #include "GaudiKernel/IHistogramSvc.h"
21 #include "CLHEP/Vector/ThreeVector.h"
22 #include "CLHEP/Vector/LorentzVector.h"
23 #include "CLHEP/Vector/TwoVector.h"
24 #include "EmcRawEvent/EmcDigi.h"
25 #include "RawEvent/RawDataUtil.h"
26 #include "MdcRawEvent/MdcDigi.h"
27
28 #include "VertexFit/IVertexDbSvc.h"
29 #include "VertexFit/KinematicFit.h"
30 #include "VertexFit/VertexFit.h"
31 #include "VertexFit/Helix.h"
32
33 #include "DstEvent/TofHitStatus.h"
34
35
36 #include "GaudiKernel/Bootstrap.h"
37 #include "GaudiKernel/ISvcLocator.h"
38 using CLHEP::Hep3Vector;
39 using CLHEP::Hep2Vector;
40 using CLHEP::HepLorentzVector;
41 #include "CLHEP/Geometry/Point3D.h"
42 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
43 typedef HepGeom::Point3D<double> HepPoint3D;
44 #endif
45 #include "CalibEventSelect/CalibEventSelect.h"
46
47 #include <vector>
48 #include "mysql.h"
49
50 typedef std::vector<int> Vint;
51 typedef std::vector<HepLorentzVector> Vp4;
52
53 const double mpi = 0.13957;
54 const double mkaon = 0.49367;
55 const double mproton = 0.93827;
56
57
58 double Px(RecMdcKalTrack *trk){
59 double phi= trk->fi0();
60 double kappa= trk->kappa();
61 double tanl = trk->tanl();
62
63 double px=-sin(phi)/fabs(kappa);
64 return px;
65 }
66
67 double Py(RecMdcKalTrack *trk){
68 double phi= trk->fi0();
69 double kappa= trk->kappa();
70 double tanl = trk->tanl();
71
72 double py=cos(phi)/fabs(kappa);
73
74 return py;
75 }
76
77
78 double Pz(RecMdcKalTrack *trk){
79 double phi= trk->fi0();
80 double kappa= trk->kappa();
81 double tanl = trk->tanl();
82
83 double pz=tanl/fabs(kappa);
84 return pz;
85 }
86
87 double P(RecMdcKalTrack *trk){
88 double phi= trk->fi0();
89 double kappa= trk->kappa();
90 double tanl = trk->tanl();
91
92 double px=-sin(phi)/fabs(kappa);
93 double py=cos(phi)/fabs(kappa);
94 double pz=tanl/fabs(kappa);
95
96 return sqrt(px*px+py*py+pz*pz);
97 }
98
99double Phi(RecMdcKalTrack *trk){
100 double phi0=trk->fi0();
101 double kappa = trk->kappa();
102 double pxy=0;
103 if(kappa!=0) pxy = 1.0/fabs(kappa);
104
105
106 double px = pxy * (-sin(phi0));
107 double py = pxy * cos(phi0);
108
109 return atan2(py,px);
110
111}
112
113
114 int CalibEventSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
115 88,88,100,100,112,112,128,128,140,140,
116 160,160,160,160,176,176,176,176,208,208,
117 208,208,240,240,240,240,256,256,256,256,
118 288,288,288};
119 /////////////////////////////////////////////////////////////////////////////
120DECLARE_COMPONENT(CalibEventSelect)
121 CalibEventSelect::CalibEventSelect(const std::string& name, ISvcLocator* pSvcLocator) :
122 Algorithm(name, pSvcLocator) {
123 //Declare the properties
124
125 declareProperty("Output", m_output = false);
126 declareProperty("Display", m_display = false);
127 declareProperty("PrintMonitor", m_printmonitor=false);
128 declareProperty("SelectRadBhabha", m_selectRadBhabha=false);
129 declareProperty("SelectGGEE", m_selectGGEE=false);
130 declareProperty("SelectGG4Pi", m_selectGG4Pi=false);
131 declareProperty("SelectRadBhabhaT", m_selectRadBhabhaT=false);
132 declareProperty("SelectRhoPi", m_selectRhoPi=false);
133 declareProperty("SelectKstarK", m_selectKstarK=false);
134 declareProperty("SelectPPPiPi", m_selectPPPiPi=false);
135 declareProperty("SelectRecoJpsi", m_selectRecoJpsi=false);
136 declareProperty("SelectBhabha", m_selectBhabha=false);
137 declareProperty("SelectDimu", m_selectDimu=false);
138 declareProperty("SelectCosmic", m_selectCosmic=false);
139 declareProperty("SelectGenProton", m_selectGenProton=false);
140 declareProperty("SelectPsipRhoPi", m_selectPsipRhoPi=false);
141 declareProperty("SelectPsipKstarK", m_selectPsipKstarK=false);
142 declareProperty("SelectPsipPPPiPi", m_selectPsipPPPiPi=false);
143 declareProperty("SelectPsippCand", m_selectPsippCand=false);
144
145 declareProperty("BhabhaScale", m_radbhabha_scale=1000);
146 declareProperty("RadBhabhaScale", m_bhabha_scale=1000);
147 declareProperty("DimuScale", m_dimu_scale=10);
148 declareProperty("CosmicScale", m_cosmic_scale=3);
149 declareProperty("ProtonScale", m_proton_scale=100);
150
151 declareProperty("CosmicLowp", m_cosmic_lowp=1.0);
152
153 declareProperty("WriteDst", m_writeDst=false);
154 declareProperty("WriteRec", m_writeRec=false);
155 declareProperty("Ecm", m_ecm=3.1);
156 declareProperty("Vr0cut", m_vr0cut=1.0);
157 declareProperty("Vz0cut", m_vz0cut=10.0);
158 declareProperty("Pt0HighCut", m_pt0HighCut=5.0);
159 declareProperty("Pt0LowCut", m_pt0LowCut=0.05);
160 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
161 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
162 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
163
164
165 }
166
167 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
168 StatusCode CalibEventSelect::initialize() {
169 MsgStream log(msgSvc(), name());
170
171 log << MSG::INFO << "in initialize()" << endmsg;
172
173 m_run=0;
174 //m_ecm=3.1;
175
176
177 h_nGood= histoSvc()->book("1D/nGoodtracks", 1, "num of good chaged tracks", 20, 0, 20);
178 h_nCharge= histoSvc()->book("1D/nCharge", 1, "net charge", 20, -10, 10);
179 h_pmaxobeam= histoSvc()->book("1D/pmax", 1, "pmax over beam ratio", 100, 0, 3);
180 h_eopmax= histoSvc()->book("1D/eopmax", 1, "e over pmax ratio", 100, 0, 3);
181 h_eopsec= histoSvc()->book("1D/eopsec", 1, "e over psec ratio", 100, 0, 3);
182 h_deltap= histoSvc()->book("1D/deltap", 1, "pmax minus psec", 100, -3, 3);
183 h_esumoecm= histoSvc()->book("1D/esumoverecm", 1, "total energy over ecm ratio", 100, 0, 3);
184 h_egmax= histoSvc()->book("1D/egmax", 1, "most energetic photon energy", 100, 0, 3);
185 h_deltaphi= histoSvc()->book("1D/deltaphi", 1, "phi_pmax minus phi_sec", 150, -4, 4);
186 h_Pt= histoSvc()->book("1D/Pt", 1, "total Transverse Momentum", 200,-1, 4);
187
188 StatusCode sc;
189
190 log << MSG::INFO << "creating sub-algorithms...." << endreq;
191
192
193
194
195
196 if(m_writeDst) {
197 sc = createSubAlgorithm( "EventWriter", "WriteDst", m_subalg1);
198 if( sc.isFailure() ) {
199 log << MSG::ERROR << "Error creating Sub-Algorithm WriteDst" <<endreq;
200 return sc;
201 } else {
202 log << MSG::INFO << "Success creating Sub-Algorithm WriteDst" <<endreq;
203 }
204 }
205
206 if(m_writeRec) {
207 sc = createSubAlgorithm( "EventWriter", "WriteRec", m_subalg2);
208 if( sc.isFailure() ) {
209 log << MSG::ERROR << "Error creating Sub-Algorithm WriteRec" <<endreq;
210 return sc;
211 } else {
212 log << MSG::INFO << "Success creating Sub-Algorithm WriteRec" <<endreq;
213 }
214 }
215
216
217 if(m_selectRadBhabha) {
218 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabha", m_subalg3);
219 if( sc.isFailure() ) {
220 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabha" <<endreq;
221 return sc;
222 } else {
223 log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabha" <<endreq;
224 }
225 }
226
227 if(m_selectGGEE) {
228 sc = createSubAlgorithm( "EventWriter", "SelectGGEE", m_subalg4);
229 if( sc.isFailure() ) {
230 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGGEE" <<endreq;
231 return sc;
232 } else {
233 log << MSG::INFO << "Success creating Sub-Algorithm SelectGGEE" <<endreq;
234 }
235 }
236
237 if(m_selectGG4Pi) {
238 sc = createSubAlgorithm( "EventWriter", "SelectGG4Pi", m_subalg5);
239 if( sc.isFailure() ) {
240 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGG4Pi" <<endreq;
241 return sc;
242 } else {
243 log << MSG::INFO << "Success creating Sub-Algorithm SelectGG4Pi" <<endreq;
244 }
245 }
246
247
248 if(m_selectRadBhabhaT) {
249 sc = createSubAlgorithm( "EventWriter", "SelectRadBhabhaT", m_subalg6);
250 if( sc.isFailure() ) {
251 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
252 return sc;
253 } else {
254 log << MSG::INFO << "Success creating Sub-Algorithm SelectRadBhabhaT" <<endreq;
255 }
256 }
257
258
259 if(m_selectRhoPi) {
260 sc = createSubAlgorithm( "EventWriter", "SelectRhoPi", m_subalg7);
261 if( sc.isFailure() ) {
262 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRhoPi" <<endreq;
263 return sc;
264 } else {
265 log << MSG::INFO << "Success creating Sub-Algorithm SelectRhoPi" <<endreq;
266 }
267 }
268
269
270 if(m_selectKstarK) {
271 sc = createSubAlgorithm( "EventWriter", "SelectKstarK", m_subalg8);
272 if( sc.isFailure() ) {
273 log << MSG::ERROR << "Error creating Sub-Algorithm SelectKstarK" <<endreq;
274 return sc;
275 } else {
276 log << MSG::INFO << "Success creating Sub-Algorithm SelectKstarK" <<endreq;
277 }
278 }
279
280
281
282 if(m_selectPPPiPi) {
283 sc = createSubAlgorithm( "EventWriter", "SelectPPPiPi", m_subalg9);
284 if( sc.isFailure() ) {
285 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPPPiPi" <<endreq;
286 return sc;
287 } else {
288 log << MSG::INFO << "Success creating Sub-Algorithm SelectPPPiPi" <<endreq;
289 }
290 }
291
292
293 if(m_selectRecoJpsi) {
294 sc = createSubAlgorithm( "EventWriter", "SelectRecoJpsi", m_subalg10);
295 if( sc.isFailure() ) {
296 log << MSG::ERROR << "Error creating Sub-Algorithm SelectRecoJpsi" <<endreq;
297 return sc;
298 } else {
299 log << MSG::INFO << "Success creating Sub-Algorithm SelectRecoJpsi" <<endreq;
300 }
301 }
302
303
304 if(m_selectBhabha) {
305 sc = createSubAlgorithm( "EventWriter", "SelectBhabha", m_subalg11);
306 if( sc.isFailure() ) {
307 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabha" <<endreq;
308 return sc;
309 } else {
310 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabha" <<endreq;
311 }
312 }
313
314 if(m_selectDimu) {
315 sc = createSubAlgorithm( "EventWriter", "SelectDimu", m_subalg12);
316 if( sc.isFailure() ) {
317 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimu" <<endreq;
318 return sc;
319 } else {
320 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimu" <<endreq;
321 }
322 }
323
324 if(m_selectCosmic) {
325 sc = createSubAlgorithm( "EventWriter", "SelectCosmic", m_subalg13);
326 if( sc.isFailure() ) {
327 log << MSG::ERROR << "Error creating Sub-Algorithm SelectCosmic" <<endreq;
328 return sc;
329 } else {
330 log << MSG::INFO << "Success creating Sub-Algorithm SelectCosmic" <<endreq;
331 }
332 }
333
334 if(m_selectGenProton) {
335 sc = createSubAlgorithm( "EventWriter", "SelectGenProton", m_subalg14);
336 if( sc.isFailure() ) {
337 log << MSG::ERROR << "Error creating Sub-Algorithm SelectGenProton" <<endreq;
338 return sc;
339 } else {
340 log << MSG::INFO << "Success creating Sub-Algorithm SelectGenProton" <<endreq;
341 }
342 }
343
344
345 if(m_selectPsipRhoPi) {
346 sc = createSubAlgorithm( "EventWriter", "SelectPsipRhoPi", m_subalg15);
347 if( sc.isFailure() ) {
348 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
349 return sc;
350 } else {
351 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipRhoPi" <<endreq;
352 }
353 }
354
355
356 if(m_selectPsipKstarK) {
357 sc = createSubAlgorithm( "EventWriter", "SelectPsipKstarK", m_subalg16);
358 if( sc.isFailure() ) {
359 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipKstarK" <<endreq;
360 return sc;
361 } else {
362 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipKstarK" <<endreq;
363 }
364 }
365
366
367 if(m_selectPsipPPPiPi) {
368 sc = createSubAlgorithm( "EventWriter", "SelectPsipPPPiPi", m_subalg17);
369 if( sc.isFailure() ) {
370 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
371 return sc;
372 } else {
373 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsipPPPiPi" <<endreq;
374 }
375 }
376
377
378 if(m_selectPsippCand) {
379 sc = createSubAlgorithm( "EventWriter", "SelectPsippCand", m_subalg18);
380 if( sc.isFailure() ) {
381 log << MSG::ERROR << "Error creating Sub-Algorithm SelectPsippCand" <<endreq;
382 return sc;
383 } else {
384 log << MSG::INFO << "Success creating Sub-Algorithm SelectPsippCand" <<endreq;
385 }
386 }
387
388
389
390
391 if(m_output) {
392 StatusCode status;
393 NTuplePtr nt0(ntupleSvc(), "FILE1/hadron");
394 if ( nt0 ) m_tuple0 = nt0;
395 else {
396 m_tuple0 = ntupleSvc()->book ("FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example");
397 if ( m_tuple0 ) {
398 status = m_tuple0->addItem ("esum", m_esum);
399 status = m_tuple0->addItem ("eemc", m_eemc);
400 status = m_tuple0->addItem ("etot", m_etot);
401 status = m_tuple0->addItem ("nGood", m_nGood);
402 status = m_tuple0->addItem ("nCharge", m_nCharge);
403 status = m_tuple0->addItem ("nGam", m_nGam);
404 status = m_tuple0->addItem ("ptot", m_ptot);
405 status = m_tuple0->addItem ("pp", m_pp);
406 status = m_tuple0->addItem ("pm", m_pm);
407 status = m_tuple0->addItem ("run", m_runnb);
408 status = m_tuple0->addItem ("event", m_evtnb);
409 status = m_tuple0->addItem ("maxE", m_maxE);
410 status = m_tuple0->addItem ("secE", m_secE);
411 status = m_tuple0->addItem ("dthe", m_dthe);
412 status = m_tuple0->addItem ("dphi", m_dphi);
413 status = m_tuple0->addItem ("mdcHit1", m_mdcHit1);
414 status = m_tuple0->addItem ("mdcHit2", m_mdcHit2);
415 }
416 else {
417 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple0) << endmsg;
418 return StatusCode::FAILURE;
419 }
420 }
421
422 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
423 if ( nt1 ) m_tuple1 = nt1;
424 else {
425 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
426 if ( m_tuple1 ) {
427 status = m_tuple1->addItem ("vx0", m_vx0);
428 status = m_tuple1->addItem ("vy0", m_vy0);
429 status = m_tuple1->addItem ("vz0", m_vz0);
430 status = m_tuple1->addItem ("vr0", m_vr0);
431 status = m_tuple1->addItem ("theta0", m_theta0);
432 status = m_tuple1->addItem ("p0", m_p0);
433 status = m_tuple1->addItem ("pt0", m_pt0);
434 }
435 else {
436 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
437 return StatusCode::FAILURE;
438 }
439 }
440
441 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
442 if ( nt2 ) m_tuple2 = nt2;
443 else {
444 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
445 if ( m_tuple2 ) {
446 status = m_tuple2->addItem ("dthe", m_dthe);
447 status = m_tuple2->addItem ("dphi", m_dphi);
448 status = m_tuple2->addItem ("dang", m_dang);
449 status = m_tuple2->addItem ("eraw", m_eraw);
450 }
451 else {
452 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
453 return StatusCode::FAILURE;
454 }
455 }
456 }
457 //
458 //--------end of book--------
459 //
460
461 m_events=0;
462 m_radBhabhaNumber=0;
463 m_GGEENumber=0;
464 m_GG4PiNumber=0;
465 m_radBhabhaTNumber=0;
466 m_rhoPiNumber=0;
467 m_kstarKNumber=0;
468 m_ppPiPiNumber=0;
469 m_recoJpsiNumber=0;
470 m_bhabhaNumber=0;
471 m_dimuNumber=0;
472 m_cosmicNumber=0;
473 m_genProtonNumber=0;
474 m_psipRhoPiNumber=0;
475 m_psipKstarKNumber=0;
476 m_psipPPPiPiNumber=0;
477
478 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
479 return StatusCode::SUCCESS;
480
481 }
482
483 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
484 StatusCode CalibEventSelect::execute() {
485
486 //setFilterPassed(false);
487
488 MsgStream log(msgSvc(), name());
489 log << MSG::INFO << "in execute()" << endreq;
490
491 if( m_writeDst) m_subalg1->execute();
492 if( m_writeRec) m_subalg2->execute();
493
494
495 m_isRadBhabha = false;
496 m_isGGEE = false;
497 m_isGG4Pi = false;
498 m_isRadBhabhaT = false;
499 m_isRhoPi = false;
500 m_isKstarK = false;
501 m_isRecoJpsi = false;
502 m_isPPPiPi = false;
503 m_isBhabha = false;
504 m_isDimu = false;
505 m_isCosmic = false;
506 m_isGenProton = false;
507 m_isPsipRhoPi = false;
508 m_isPsipKstarK = false;
509 m_isPsipPPPiPi = false;
510 m_isPsippCand = false;
511
512 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
513 if(!eventHeader)
514 {
515 cout<<" eventHeader "<<endl;
516 return StatusCode::FAILURE;
517 }
518
519 int run=eventHeader->runNumber();
520 int event=eventHeader->eventNumber();
521
522 //get run by run ebeam
523 if(m_run !=run){
524 m_run=run;
525 double beamE=0;
526 readbeamEfromDb(run,beamE);
527 cout<<"new run is:"<<m_run<<endl;
528 cout<<"beamE is:"<<beamE<<endl;
529 if(beamE>0 && beamE<3)
530 m_ecm=2*beamE;
531 }
532
533
534
535 if(m_display && m_events%1000==0){
536 cout<< m_events << " -------- run,event: "<<run<<","<<event<<endl;
537 cout<<"m_ecm="<<m_ecm<<endl;
538 }
539 m_events++;
540
541 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
542 if(!evtRecEvent ) {
543 cout<<" evtRecEvent "<<endl;
544 return StatusCode::FAILURE;
545 }
546
547 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
548 << evtRecEvent->totalCharged() << " , "
549 << evtRecEvent->totalNeutral() << " , "
550 << evtRecEvent->totalTracks() <<endreq;
551 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
552 if(!evtRecTrkCol){
553 cout<<" evtRecTrkCol "<<endl;
554 return StatusCode::FAILURE;
555 }
556
557 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size()) return StatusCode::SUCCESS;
558
559
560 //get pi0 reconstruction
561 SmartDataPtr<EvtRecPi0Col> recPi0Col(eventSvc(), "/Event/EvtRec/EvtRecPi0Col");
562 if ( ! recPi0Col ) {
563 log << MSG::FATAL << "Could not find EvtRecPi0Col" << endreq;
564 return StatusCode::FAILURE;
565 }
566
567
568 int nPi0 = recPi0Col->size();
569 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
570 if(nPi0==1){
571 double mpi0 = (*itpi0)->unconMass();
572 if ( fabs( mpi0 - 0.135 )> 0.015 )
573 nPi0=0;
574 }
575
576 /*
577 int nPi0=0;
578 for(EvtRecPi0Col::iterator it=itpi0; it!= recPi0Col->end(); it++){
579 double mpi0 = (*itpi0)->unconMass();
580 if ( fabs( mpi0 - 0.135 )<= 0.015 )
581 nPi0++;
582 }
583 */
584
585
586 // -------- Good Charged Track Selection
587 Vint iGood;
588 iGood.clear();
589 vector<int> iCP, iCN;
590 iCP.clear();
591 iCN.clear();
592
593 Vint iCosmicGood;
594 iCosmicGood.clear();
595
596 int nCharge = 0;
597 int nCosmicCharge =0;
598
599 for(int i = 0;i < evtRecEvent->totalCharged(); i++)
600 {
601 //if(i>=evtRecTrkCol->size()) break;
602 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
603
604 //if(!(*itTrk)->isMdcTrackValid()) continue;
605 //RecMdcTrack *mdcTrk =(*itTrk)->mdcTrack();
606 if(!(*itTrk)->isMdcKalTrackValid()) continue;
607 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
608 // mdcTrk->setPidType(RecMdcKalTrack::electron);
609
610
611 double vx0 = mdcTrk->x();
612 double vy0 = mdcTrk->y();
613 double vz0 = mdcTrk->z();
614 double vr0 = mdcTrk->r();
615 double theta0 = mdcTrk->theta();
616 double p0 = P(mdcTrk);
617 double pt0 = sqrt( pow(Px(mdcTrk),2)+pow(Py(mdcTrk),2));
618
619
620 if(m_output) {
621 m_vx0 = vx0;
622 m_vy0 = vy0;
623 m_vz0 = vz0;
624 m_vr0 = vr0;
625 m_theta0 = theta0;
626 m_p0 = p0;
627 m_pt0 = pt0;
628 m_tuple1->write();
629 }
630
631 //db
632
633 Hep3Vector xorigin(0,0,0);
634 IVertexDbSvc* vtxsvc;
635 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
636 if(vtxsvc->isVertexValid()){
637 double* dbv = vtxsvc->PrimaryVertex();
638 double* vv = vtxsvc->SigmaPrimaryVertex();
639 xorigin.setX(dbv[0]);
640 xorigin.setY(dbv[1]);
641 xorigin.setZ(dbv[2]);
642 }
643 HepVector a = mdcTrk->getZHelix();
644 HepSymMatrix Ea = mdcTrk->getZError();
645 HepPoint3D point0(0.,0.,0.); // the initial point for MDC recosntruction
646 HepPoint3D IP(xorigin[0],xorigin[1],xorigin[2]);
647 VFHelix helixip(point0,a,Ea);
648 helixip.pivot(IP);
649 HepVector vecipa = helixip.a();
650 double db=fabs(vecipa[0]);
651 double dz=vecipa[3];
652
653
654
655
656 if(fabs(dz) >= m_vz0cut) continue;
657 if(db >= m_vr0cut) continue;
658
659 //cosmic tracks
660 if(p0>m_cosmic_lowp && p0<20){
661 iCosmicGood.push_back((*itTrk)->trackId());
662 nCosmicCharge += mdcTrk->charge();
663 }
664
665
666
667 if(pt0 >= m_pt0HighCut) continue;
668 if(pt0 <= m_pt0LowCut) continue;
669
670 iGood.push_back((*itTrk)->trackId());
671 nCharge += mdcTrk->charge();
672 if(mdcTrk->charge()>0)
673 iCP.push_back((*itTrk)->trackId());
674 else if(mdcTrk->charge()<0)
675 iCN.push_back((*itTrk)->trackId());
676
677
678 }
679 int nGood = iGood.size();
680 int nCosmicGood = iCosmicGood.size();
681
682 // -------- Good Photon Selection
683 Vint iGam;
684 iGam.clear();
685 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
686 //if(i>=evtRecTrkCol->size()) break;
687 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
688 if(!(*itTrk)->isEmcShowerValid()) continue;
689 RecEmcShower *emcTrk = (*itTrk)->emcShower();
690 double eraw = emcTrk->energy();
691 double time = emcTrk->time();
692 double costh = cos(emcTrk->theta());
693 if(fabs(costh)<0.83 && eraw<0.025) continue;
694 if(fabs(costh)>=0.83 && eraw<0.05) continue;
695 if(time<0 || time>14) continue;
696
697 iGam.push_back((*itTrk)->trackId());
698 }
699 int nGam = iGam.size();
700
701 // -------- Assign 4-momentum to each charged track
702 Vint ipip, ipim;
703 ipip.clear();
704 ipim.clear();
705 Vp4 ppip, ppim;
706 ppip.clear();
707 ppim.clear();
708
709 //cout<<"charged track:"<<endl;
710 double echarge = 0.; //total energy of charged track
711 double ptot = 0.; //total momentum of charged track
712 double etot = 0.; //total energy in MDC and EMC
713 double eemc = 0.; //total energy in EMC
714 double pp = 0.;
715 double pm = 0.;
716 double pmax=0.0;
717 double psec=0.0;
718 double eccmax=0.0;
719 double eccsec=0.0;
720 double phimax=0.0;
721 double phisec=0.0;
722 double eopmax=0.0;
723 double eopsec=0.0;
724 Hep3Vector Pt_charge(0,0,0);
725
726 for(int i = 0; i < nGood; i++) {
727 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
728
729 double p3=0;
730 //if((*itTrk)->isMdcTrackValid()) {
731 //RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
732
733 if((*itTrk)->isMdcKalTrackValid()) {
734 RecMdcKalTrack *mdcTrk = (*itTrk)->mdcKalTrack();
735 mdcTrk->setPidType(RecMdcKalTrack::electron);
736
737 ptot += P(mdcTrk);
738
739 //double phi= mdcTrk->phi();
740 double phi= Phi(mdcTrk);
741
742
743 HepLorentzVector ptrk;
744 ptrk.setPx(Px(mdcTrk));
745 ptrk.setPy(Py(mdcTrk));
746 ptrk.setPz(Pz(mdcTrk));
747 p3 = fabs(ptrk.mag());
748
749
750
751 //cout<<"p3 before compa="<<p3<<endl;
752 //cout<<"pmax before compa ="<<pmax<<endl;
753 //cout<<"psec before compa ="<<psec<<endl;
754
755 Hep3Vector ptemp(Px(mdcTrk),Py(mdcTrk),0);
756 Pt_charge+=ptemp;
757
758 double ecc=0.0;
759 if((*itTrk)->isEmcShowerValid()) {
760 RecEmcShower* emcTrk = (*itTrk)->emcShower();
761 ecc = emcTrk->energy();
762 }
763
764 if( p3 > pmax){
765 psec=pmax;
766 eccsec=eccmax;
767 phisec=phimax;
768 pmax=p3;
769 eccmax=ecc;
770 phimax=phi;
771 }
772 else if( p3 < pmax && p3> psec){
773 psec=p3;
774 eccsec=ecc;
775 phisec=phi;
776 }
777
778 // cout<<"p3 after compa="<<p3<<endl;
779 // cout<<"pmax after compa ="<<pmax<<endl;
780 //cout<<"psec after compa ="<<psec<<endl;
781
782 ptrk.setE(sqrt(p3*p3+mpi*mpi));
783 ptrk = ptrk.boost(-0.011,0,0);//boost to cms
784
785
786 echarge += ptrk.e();
787 etot += ptrk.e();
788
789 if(mdcTrk->charge() >0) {
790 ppip.push_back(ptrk);
791 pp = ptrk.rho();
792 } else {
793 ppim.push_back(ptrk);
794 pm = ptrk.rho();
795 }
796
797 }
798
799 if((*itTrk)->isEmcShowerValid()) {
800
801 RecEmcShower* emcTrk = (*itTrk)->emcShower();
802 double eraw = emcTrk->energy();
803 double phiemc = emcTrk->phi();
804 double the = emcTrk->theta();
805
806 HepLorentzVector pemc;
807 pemc.setPx(eraw*sin(the)*cos(phiemc));
808 pemc.setPy(eraw*sin(the)*sin(phiemc));
809 pemc.setPz(eraw*cos(the));
810 pemc.setE(eraw);
811 pemc = pemc.boost(-0.011,0,0);// boost to cms
812
813 // eemc += eraw;
814 etot += pemc.e();
815
816 }
817
818
819
820
821 } //end of looping over good charged track
822
823 if(pmax!=0) eopmax=eccmax/pmax;
824 if(psec!=0) eopsec=eccsec/psec;
825
826 eemc=eccmax+eccsec;
827
828 /*
829 if(nGood>1){
830 cout<<"pmax="<<pmax<<endl;
831 cout<<"psec="<<psec<<endl;
832 cout<<"eopmax="<<eopmax<<endl;
833 cout<<"dphi-180="<< (fabs(phimax-phisec)-CLHEP::pi)*180/CLHEP::pi <<endl;
834 }
835 */
836
837 // -------- Assign 4-momentum to each photon
838 //cout<<"neutral:"<<endl;
839 Vp4 pGam;
840 pGam.clear();
841 double eneu=0; //total energy of neutral track
842 double egmax=0;
843 Hep3Vector Pt_neu(0,0,0);
844
845 for(int i = 0; i < nGam; i++) {
846 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
847 RecEmcShower* emcTrk = (*itTrk)->emcShower();
848 double eraw = emcTrk->energy();
849 double phi = emcTrk->phi();
850 double the = emcTrk->theta();
851 HepLorentzVector ptrk;
852 ptrk.setPx(eraw*sin(the)*cos(phi));
853 ptrk.setPy(eraw*sin(the)*sin(phi));
854 ptrk.setPz(eraw*cos(the));
855 ptrk.setE(eraw);
856 ptrk = ptrk.boost(-0.011,0,0);// boost to cms
857 pGam.push_back(ptrk);
858 eneu += ptrk.e();
859 etot += ptrk.e();
860 eemc += eraw;
861
862 Hep3Vector ptemp(eraw*sin(the)*cos(phi), eraw*sin(the)*sin(phi),0);
863 Pt_neu+=ptemp;
864
865 if(ptrk.e()>egmax)
866 egmax=ptrk.e();
867
868 }
869
870 double esum = echarge + eneu;
871 Hep3Vector Pt=Pt_charge+Pt_neu;
872
873
874 double phid=phimax-phisec;
875 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+pi, CLHEP::twopi) - CLHEP::pi;
876
877 // -------- Select each event
878
879
880 //select bhabha/dimu/cosmic events, need prescale
881
882 if( nGood == 2 && nCharge==0 && (m_selectBhabha || m_selectDimu) ){
883
884 //bhabha
885 if( m_events%m_bhabha_scale == 0 && m_selectBhabha && echarge/m_ecm>0.8 && eopmax>0.8 && eopsec>0.8){
886 m_isBhabha=true;
887 m_bhabhaNumber++;
888 }
889
890
891 //dimu
892 if( iCP.size()==1 && iCN.size()==1 && m_events%m_dimu_scale == 0 && m_selectDimu && eemc/m_ecm<0.3){
893
894 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCP[0];
895 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCN[0];
896
897 /*
898 double time1=-99,depth1=-99;
899 double time2=-99,depth2=-99;
900 if( (*itp)->isTofTrackValid() ){
901 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
902 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
903
904 for(;iter_tof!=tofTrkCol.end();iter_tof++){
905 TofHitStatus* status =new TofHitStatus;
906 status->setStatus( (*iter_tof)->status() );
907 if(status->is_cluster()){
908 time1=(*iter_tof)->tof();
909 }
910 delete status;
911 }
912 }
913
914 if( (*itp)->isMucTrackValid() ){
915 RecMucTrack* mucTrk=(*itp)->mucTrack();
916 depth1=mucTrk->depth();
917 }
918
919 if( (*itn)->isTofTrackValid() ){
920 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
921 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
922
923 for(;iter_tof!=tofTrkCol.end();iter_tof++){
924 TofHitStatus* status =new TofHitStatus;
925 status->setStatus( (*iter_tof)->status() );
926 if(status->is_cluster()){
927 time2=(*iter_tof)->tof();
928 }
929 delete status;
930 }
931 }
932
933 if( (*itn)->isMucTrackValid() ){
934 RecMucTrack* mucTrk=(*itn)->mucTrack();
935 depth2=mucTrk->depth();
936 }
937
938
939 //dimu
940 //if( m_selectDimu && time1!=-99 && time2!=-99 && fabs(time1-time2)<5 && depth1>5 && depth2>5){ //tight, no endcap
941 // if( eopmax<0.3 && eopsec<0.3 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi<6 && fabs(psec)>m_ecm/4.0 ){ //tight, no rad dimu
942 // if( eemc<0.3*m_ecm && (fabs(pmax)>0.45*m_ecm || fabs(psec)>0.45*m_ecm) ){
943 if( ((fabs(pmax)/0.5/m_ecm>0.15 && fabs(pmax)/0.5/m_ecm<.75) || (fabs(psec)/0.5/m_ecm>0.15 && fabs(psec)/0.5/m_ecm<.75)) && (eopmax<0.4 || eopsec<0.4)
944 && (depth1>=3 || depth2>=3)){
945 m_isDimu=true;
946 m_dimuNumber++;
947 }
948 */
949
950 double time1=-99,depth1=-99;
951 double time2=-99,depth2=-99;
952 double p1=-99, p2=-99;
953 double emc1=0, emc2=0;
954 if( (*itp)->isTofTrackValid() ){
955 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
956 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
957
958 for(;iter_tof!=tofTrkCol.end();iter_tof++){
959 TofHitStatus* status =new TofHitStatus;
960 status->setStatus( (*iter_tof)->status() );
961 if(status->is_cluster()){
962 time1=(*iter_tof)->tof();
963 }
964 delete status;
965 }
966 }
967
968 if( (*itp)->isMucTrackValid() ){
969 RecMucTrack* mucTrk=(*itp)->mucTrack();
970 depth1=mucTrk->depth();
971 }
972
973
974 if((*itp)->isMdcKalTrackValid()) {
975 RecMdcKalTrack *mdcTrk = (*itp)->mdcKalTrack();
976 mdcTrk->setPidType(RecMdcKalTrack::muon);
977 p1= P(mdcTrk);
978 }
979
980 if((*itp)->isEmcShowerValid()) {
981 RecEmcShower* emcTrk = (*itp)->emcShower();
982 emc1= emcTrk->energy();
983 }
984
985 if( (*itn)->isTofTrackValid() ){
986 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
987 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
988
989 for(;iter_tof!=tofTrkCol.end();iter_tof++){
990 TofHitStatus* status =new TofHitStatus;
991 status->setStatus( (*iter_tof)->status() );
992 if(status->is_cluster()){
993 time2=(*iter_tof)->tof();
994 }
995 delete status;
996 }
997 }
998
999 if( (*itn)->isMucTrackValid() ){
1000 RecMucTrack* mucTrk=(*itn)->mucTrack();
1001 depth2=mucTrk->depth();
1002 }
1003
1004 if((*itn)->isMdcKalTrackValid()) {
1005 RecMdcKalTrack *mdcTrk = (*itn)->mdcKalTrack();
1006 mdcTrk->setPidType(RecMdcKalTrack::muon);
1007 p2= P(mdcTrk);
1008 }
1009
1010
1011 if((*itn)->isEmcShowerValid()) {
1012 RecEmcShower* emcTrk = (*itn)->emcShower();
1013 emc2= emcTrk->energy();
1014 }
1015
1016 double ebeam=m_ecm/2.0;
1017 if( fabs(time1-time2)<5 && ((p1/ebeam>0.15 && p1/ebeam<0.75) || (p2/ebeam>0.15 && p2/ebeam<0.75))
1018 && (emc1/p1<0.4 || emc2/p2<0.4) && (depth1>3 || depth2>3) &&
1019 emc1>0.15 && emc1<0.3 && emc2>0.15 && emc2<0.3 && emc1/p1<0.8 && emc2/p2<0.8 &&
1020 ((depth1>80*p1-60 || depth1>40) || (depth2>80*p2-60 || depth2>40)) ){
1021 m_isDimu=true;
1022 m_dimuNumber++;
1023 }
1024
1025 }//end of dimu
1026
1027
1028
1029 }//end of bhabha, dimu
1030
1031
1032
1033 if(m_selectCosmic && nCosmicGood == 2 && eemc/m_ecm<0.3){
1034
1035 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCosmicGood[0];
1036 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCosmicGood[1];
1037
1038 double time1=-99,depth1=-99;
1039 double time2=-99,depth2=-99;
1040 if( (*itp)->isTofTrackValid() ){
1041 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
1042 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1043
1044 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1045 TofHitStatus* status =new TofHitStatus;
1046 status->setStatus( (*iter_tof)->status() );
1047 if(status->is_cluster()){
1048 time1=(*iter_tof)->tof();
1049 }
1050 delete status;
1051 }
1052 }
1053
1054 if( (*itp)->isMucTrackValid() ){
1055 RecMucTrack* mucTrk=(*itp)->mucTrack();
1056 depth1=mucTrk->depth();
1057 }
1058
1059 if( (*itn)->isTofTrackValid() ){
1060 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
1061 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
1062
1063 for(;iter_tof!=tofTrkCol.end();iter_tof++){
1064 TofHitStatus* status =new TofHitStatus;
1065 status->setStatus( (*iter_tof)->status() );
1066 if(status->is_cluster()){
1067 time2=(*iter_tof)->tof();
1068 }
1069 delete status;
1070 }
1071 }
1072
1073 if( (*itn)->isMucTrackValid() ){
1074 RecMucTrack* mucTrk=(*itn)->mucTrack();
1075 depth2=mucTrk->depth();
1076 }
1077
1078 //cosmic
1079 //if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 && depth1>5 && depth2>5){
1080 if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 ){
1081 m_isCosmic=true;
1082 m_cosmicNumber++;
1083 }
1084
1085
1086 }//end of cosmic
1087
1088
1089 //select generic protons
1090
1091 if(m_events%m_proton_scale ==0 ){
1092
1093 bool protoncand=false;
1094 int ncharge=0;
1095 int nproton=0;
1096 for(int i=0; i<nGood; i++){
1097
1098 EvtRecTrackIterator iter = evtRecTrkCol->begin() + iGood[i];
1099 RecMdcKalTrack *mdcTrk = (*iter)->mdcKalTrack();
1100 mdcTrk->setPidType(RecMdcKalTrack::proton);
1101
1102 double pp = P(mdcTrk);
1103 double dedx=-99;
1104 double chiP=-99;
1105
1106 if((*iter)->isMdcDedxValid()){
1107 RecMdcDedx* dedxTrk=(*iter)->mdcDedx();
1108 dedx=dedxTrk->normPH();
1109 chiP=dedxTrk->chiP();
1110 }
1111
1112 if( fabs(pp)<0.6 && dedx>1.2 && fabs(chiP)<6){
1113 ncharge+=mdcTrk->charge();
1114 nproton++;
1115 //protoncand=true;
1116 //break;
1117 }
1118 }
1119 if( (nproton==1 && ncharge==-1) || (nproton>=2 && ncharge<=nproton-2))
1120 protoncand=true;
1121 if( protoncand ){
1122 m_isGenProton=true;
1123 m_genProtonNumber++;
1124 }
1125
1126 }//end of generic proton
1127
1128
1129 //fill monitoring histograms for rad bhabha
1130
1131
1132 if(m_printmonitor){
1133 h_nGood->fill(nGood);
1134 h_nCharge->fill(nCharge);
1135 h_pmaxobeam->fill(pmax/(m_ecm/2.0));
1136 h_eopmax->fill(eopmax);
1137 h_eopsec->fill(eopsec);
1138 h_deltap->fill(pmax-psec);
1139 h_esumoecm->fill(esum/m_ecm);
1140 h_egmax->fill(egmax);
1141 h_deltaphi->fill(phid);
1142 h_Pt->fill(Pt.mag());
1143 }
1144
1145
1146
1147 //select radbhabha
1148 if(nGood>1 && pmax/(m_ecm/2.0)>0.4 && eopmax>0.5 && esum/m_ecm>0.75 &&
1149 egmax>0.08 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi>2.86 )
1150 {
1151 m_isRadBhabha=true;
1152 m_radBhabhaNumber++;
1153 }
1154 //select radbhabha tight
1155 if( m_isRadBhabha )
1156 {
1157 //prescale high momentum tracks
1158 if(nGood==2 && nCharge==0 && eemc/m_ecm>=0.7 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 ){
1159
1160 if( fabs(fabs(pmax)-m_ecm/2.0)<0.1 && fabs(fabs(psec)-m_ecm/2.0)<0.1 ){
1161 if(m_events%1000==0){
1162 m_isRadBhabhaT=true;
1163 m_radBhabhaTNumber++;
1164 }
1165 }
1166 else{
1167 m_isRadBhabhaT=true;
1168 m_radBhabhaTNumber++;
1169 }
1170
1171 }
1172
1173
1174
1175 }
1176
1177 //select ggee events, (exclude radee tight)
1178 //if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopsec>=0.85 && eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
1179 if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 && eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
1180 {
1181 m_isGGEE=true;
1182 m_GGEENumber++;
1183 }
1184
1185 //select gg4pi events
1186 if(m_selectGG4Pi && nGood==4 && nCharge==0 && pmax/(m_ecm/2.0)>0.04 && pmax/(m_ecm/2.0)<0.9 && esum/m_ecm>0.04 && esum/m_ecm<0.75 && Pt.mag()<=0.2)
1187 {
1188 m_isGG4Pi=true;
1189 m_GG4PiNumber++;
1190 }
1191
1192 //select rhopi/kstark events
1193 if( (m_selectRhoPi || m_selectKstarK) && nGood == 2 && nCharge==0 && nPi0 == 1){
1194
1195 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iGood[0];
1196 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iGood[1];
1197
1198
1199 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1200
1201 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1202 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1203
1204 HepLorentzVector p4trk1;
1205 p4trk1.setPx(Px(trk1));
1206 p4trk1.setPy(Py(trk1));
1207 p4trk1.setPz(Pz(trk1));
1208 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
1209
1210 HepLorentzVector p4trk2;
1211 p4trk2.setPx(Px(trk2));
1212 p4trk2.setPy(Py(trk2));
1213 p4trk2.setPz(Pz(trk2));
1214 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
1215
1216
1217 HepLorentzVector p4trk3;
1218 p4trk3.setPx(Px(trk1));
1219 p4trk3.setPy(Py(trk1));
1220 p4trk3.setPz(Pz(trk1));
1221 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1222
1223 HepLorentzVector p4trk4;
1224 p4trk4.setPx(Px(trk2));
1225 p4trk4.setPy(Py(trk2));
1226 p4trk4.setPz(Pz(trk2));
1227 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1228
1229
1230 //EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1231 itpi0 = recPi0Col->begin();
1232 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1233 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1234 HepLorentzVector p4pi0 = p4gam1+p4gam2;
1235
1236 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0; //kstark
1237 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0; // rhopi
1238
1239 double rhopimass=p4total2.m();
1240 double kstarkmass=p4total.m();
1241 if( (kstarkmass > 3.0 && kstarkmass < 3.15) || (rhopimass > 3.0 && rhopimass < 3.15) ){
1242
1243 //tight cuts
1244
1245 //remove bhabha background
1246 double eop1=0.0,eop2=0.0;
1247 if( (*itone)->isEmcShowerValid() ){
1248 RecEmcShower *emcTrk = (*itone)->emcShower();
1249 double etrk=emcTrk->energy();
1250 //cout<<"trk1 p="<<P(trk1)<<endl;
1251 if(P(trk1)!=0){
1252 eop1 = etrk/P(trk1);
1253 //if( fabs(eop1)> 0.8 ) return StatusCode::SUCCESS;
1254 }
1255 }
1256
1257 if( (*ittwo)->isEmcShowerValid() ){
1258 RecEmcShower *emcTrk = (*ittwo)->emcShower();
1259 double etrk=emcTrk->energy();
1260 //cout<<"trk2 p="<<P(trk2)<<endl;
1261 if(P(trk2)!=0){
1262 eop2 = etrk/P(trk2);
1263 //if( fabs(eop2)> 0.8 ) return StatusCode::SUCCESS;
1264 }
1265 }
1266
1267 if(eop1<0.8 && eop2<0.8){
1268
1269 if(rhopimass>3.0 && rhopimass<3.15){
1270 //kinematic fit
1271 //HepLorentzVector ecms(0.034,0,0,3.097);
1272 HepLorentzVector ecms(0.034,0,0,m_ecm);
1273
1274
1275 WTrackParameter wvpiTrk1, wvpiTrk2;
1276 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1277 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1278
1279 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1280 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1281 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1282 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1283
1284 KinematicFit * kmfit = KinematicFit::instance();
1285 kmfit->init();
1286 kmfit->AddTrack(0, wvpiTrk1);
1287 kmfit->AddTrack(1, wvpiTrk2);
1288 kmfit->AddTrack(2, 0.0, gShr1);
1289 kmfit->AddTrack(3, 0.0, gShr2);
1290 kmfit->AddFourMomentum(0, ecms);
1291
1292 bool oksq1 = kmfit->Fit();
1293 double chi1 = kmfit->chisq();
1294
1295 //
1296 if(oksq1 && chi1<=60){
1297 m_isRhoPi = true;
1298 m_rhoPiNumber++;
1299 }
1300 } //end of selecting rhopi
1301
1302
1303 if(kstarkmass>3.0 && kstarkmass<3.15){
1304
1305 //kstark resonances
1306 double mkstarp=0, mkstarm=0;
1307 if( trk1->charge() >0 ){
1308 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1309 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1310 mkstarp = p4kstarp.m();
1311 mkstarm = p4kstarm.m();
1312 }
1313 else{
1314 HepLorentzVector p4kstarm = p4trk1 + p4pi0;
1315 HepLorentzVector p4kstarp = p4trk2 + p4pi0;
1316 mkstarp = p4kstarp.m();
1317 mkstarm = p4kstarm.m();
1318 }
1319
1320 if ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1) || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) ){
1321 //kinematic fit
1322 //HepLorentzVector ecms(0.034,0,0,3.097);
1323 HepLorentzVector ecms(0.034,0,0,m_ecm);
1324
1325 WTrackParameter wvpiTrk1, wvpiTrk2;
1326 wvpiTrk1 = WTrackParameter(mkaon, trk1->getZHelix(), trk1->getZError());
1327 wvpiTrk2 = WTrackParameter(mkaon, trk2->getZHelix(), trk2->getZError());
1328
1329 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1330 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1331 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1332 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1333
1334 KinematicFit * kmfit = KinematicFit::instance();
1335 kmfit->init();
1336 kmfit->AddTrack(0, wvpiTrk1);
1337 kmfit->AddTrack(1, wvpiTrk2);
1338 kmfit->AddTrack(2, 0.0, gShr1);
1339 kmfit->AddTrack(3, 0.0, gShr2);
1340 kmfit->AddFourMomentum(0, ecms);
1341
1342 bool oksq1 = kmfit->Fit();
1343 double chi1 = kmfit->chisq();
1344 //
1345
1346 if(oksq1 && chi1<=60){
1347 m_isKstarK = true;
1348 m_kstarKNumber++;
1349 }
1350
1351 }
1352
1353 } //end of selecting kstark
1354
1355 }//end of non di-electron
1356
1357 //end of tight cuts
1358
1359 }
1360 }
1361
1362
1363
1364 } //end of selecting rhopi/kstark events
1365
1366 //select ppPiPi events
1367 if(m_selectPPPiPi && nGood ==4 && nCharge == 0 && nPi0==0){
1368
1369 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iCP[0];
1370 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iCP[1];
1371 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iCN[0];
1372 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iCN[1];
1373 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1374 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1375 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
1376 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
1377
1378 HepLorentzVector p4trkpp;
1379 HepLorentzVector p4trkpm;
1380 HepLorentzVector p4trkpip;
1381 HepLorentzVector p4trkpim;
1382
1383 //hypothesis 1
1384
1385 trk1->setPidType(RecMdcKalTrack::proton);
1386 p4trkpp.setPx(Px(trk1));
1387 p4trkpp.setPy(Py(trk1));
1388 p4trkpp.setPz(Pz(trk1));
1389 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1390
1391 trk3->setPidType(RecMdcKalTrack::proton);
1392 p4trkpm.setPx(Px(trk3));
1393 p4trkpm.setPy(Py(trk3));
1394 p4trkpm.setPz(Pz(trk3));
1395 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1396
1397 trk2->setPidType(RecMdcKalTrack::pion);
1398 p4trkpip.setPx(Px(trk2));
1399 p4trkpip.setPy(Py(trk2));
1400 p4trkpip.setPz(Pz(trk2));
1401 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1402
1403 trk4->setPidType(RecMdcKalTrack::pion);
1404 p4trkpim.setPx(Px(trk4));
1405 p4trkpim.setPy(Py(trk4));
1406 p4trkpim.setPz(Pz(trk4));
1407 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1408
1409 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1410
1411
1412 //hypothesis 2
1413
1414 trk1->setPidType(RecMdcKalTrack::proton);
1415 p4trkpp.setPx(Px(trk1));
1416 p4trkpp.setPy(Py(trk1));
1417 p4trkpp.setPz(Pz(trk1));
1418 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1419
1420 trk4->setPidType(RecMdcKalTrack::proton);
1421 p4trkpm.setPx(Px(trk4));
1422 p4trkpm.setPy(Py(trk4));
1423 p4trkpm.setPz(Pz(trk4));
1424 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1425
1426 trk2->setPidType(RecMdcKalTrack::pion);
1427 p4trkpip.setPx(Px(trk2));
1428 p4trkpip.setPy(Py(trk2));
1429 p4trkpip.setPz(Pz(trk2));
1430 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1431
1432 trk3->setPidType(RecMdcKalTrack::pion);
1433 p4trkpim.setPx(Px(trk3));
1434 p4trkpim.setPy(Py(trk3));
1435 p4trkpim.setPz(Pz(trk3));
1436 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1437
1438 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1439
1440
1441 //hypothesis 3
1442
1443
1444 trk2->setPidType(RecMdcKalTrack::proton);
1445 p4trkpp.setPx(Px(trk2));
1446 p4trkpp.setPy(Py(trk2));
1447 p4trkpp.setPz(Pz(trk2));
1448 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1449
1450 trk3->setPidType(RecMdcKalTrack::proton);
1451 p4trkpm.setPx(Px(trk3));
1452 p4trkpm.setPy(Py(trk3));
1453 p4trkpm.setPz(Pz(trk3));
1454 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1455
1456 trk1->setPidType(RecMdcKalTrack::pion);
1457 p4trkpip.setPx(Px(trk1));
1458 p4trkpip.setPy(Py(trk1));
1459 p4trkpip.setPz(Pz(trk1));
1460 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1461
1462 trk4->setPidType(RecMdcKalTrack::pion);
1463 p4trkpim.setPx(Px(trk4));
1464 p4trkpim.setPy(Py(trk4));
1465 p4trkpim.setPz(Pz(trk4));
1466 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1467
1468 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1469
1470
1471 //hypothesis 4
1472
1473
1474 trk2->setPidType(RecMdcKalTrack::proton);
1475 p4trkpp.setPx(Px(trk2));
1476 p4trkpp.setPy(Py(trk2));
1477 p4trkpp.setPz(Pz(trk2));
1478 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1479
1480 trk4->setPidType(RecMdcKalTrack::proton);
1481 p4trkpm.setPx(Px(trk4));
1482 p4trkpm.setPy(Py(trk4));
1483 p4trkpm.setPz(Pz(trk4));
1484 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1485
1486 trk1->setPidType(RecMdcKalTrack::pion);
1487 p4trkpip.setPx(Px(trk1));
1488 p4trkpip.setPy(Py(trk1));
1489 p4trkpip.setPz(Pz(trk1));
1490 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1491
1492 trk3->setPidType(RecMdcKalTrack::pion);
1493 p4trkpim.setPx(Px(trk3));
1494 p4trkpim.setPy(Py(trk3));
1495 p4trkpim.setPz(Pz(trk3));
1496 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1497
1498 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1499
1500
1501
1502
1503 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1504 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1505
1506 //tight cuts
1507
1508 //kinematic fits
1509 double chi1, chi2, chi3, chi4;
1510 int type=0;
1511 double chimin=-999;
1512 HepLorentzVector ecms(0.034,0,0,m_ecm);
1513
1514 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 ;
1515
1516 {
1517 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1518 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1519 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1520 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1521
1522
1523 KinematicFit * kmfit = KinematicFit::instance();
1524 kmfit->init();
1525 kmfit->AddTrack(0, wvpiTrk1);
1526 kmfit->AddTrack(1, wvpiTrk2);
1527 kmfit->AddTrack(2, wvpiTrk3);
1528 kmfit->AddTrack(3, wvpiTrk4);
1529 kmfit->AddFourMomentum(0, ecms);
1530
1531 bool oksq1 = kmfit->Fit();
1532 chi1 = kmfit->chisq();
1533 if(oksq1) {
1534 chimin=chi1;
1535 type=1;
1536 }
1537 //
1538 }
1539
1540
1541 {
1542 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1543 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1544 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1545 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1546
1547
1548 KinematicFit * kmfit = KinematicFit::instance();
1549 kmfit->init();
1550 kmfit->AddTrack(0, wvpiTrk1);
1551 kmfit->AddTrack(1, wvpiTrk2);
1552 kmfit->AddTrack(2, wvpiTrk3);
1553 kmfit->AddTrack(3, wvpiTrk4);
1554 kmfit->AddFourMomentum(0, ecms);
1555
1556 bool oksq1 = kmfit->Fit();
1557 chi2 = kmfit->chisq();
1558 if(oksq1){
1559 if(type==0){
1560 chimin=chi2;
1561 type=2;
1562 }
1563 else if(chi2<chimin){
1564 chimin=chi2;
1565 type=2;
1566 }
1567 }
1568 //
1569 }
1570
1571 {
1572 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1573 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
1574 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1575 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1576
1577 KinematicFit * kmfit = KinematicFit::instance();
1578 kmfit->init();
1579 kmfit->AddTrack(0, wvpiTrk1);
1580 kmfit->AddTrack(1, wvpiTrk2);
1581 kmfit->AddTrack(2, wvpiTrk3);
1582 kmfit->AddTrack(3, wvpiTrk4);
1583
1584 kmfit->AddFourMomentum(0, ecms);
1585
1586 bool oksq1 = kmfit->Fit();
1587 chi3 = kmfit->chisq();
1588 if(oksq1){
1589 if(type==0){
1590 chimin=chi3;
1591 type=3;
1592 }
1593 else if(chi3<chimin){
1594 chimin=chi3;
1595 type=3;
1596 }
1597 }
1598
1599 //
1600 }
1601
1602
1603 {
1604 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1605 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
1606 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1607 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1608
1609 KinematicFit * kmfit = KinematicFit::instance();
1610 kmfit->init();
1611 kmfit->AddTrack(0, wvpiTrk1);
1612 kmfit->AddTrack(1, wvpiTrk2);
1613 kmfit->AddTrack(2, wvpiTrk3);
1614 kmfit->AddTrack(3, wvpiTrk4);
1615
1616 kmfit->AddFourMomentum(0, ecms);
1617
1618 bool oksq1 = kmfit->Fit();
1619 chi4 = kmfit->chisq();
1620
1621 if(oksq1){
1622 if(type==0){
1623 chimin=chi4;
1624 type=4;
1625 }
1626 else if(chi4<chimin){
1627 chimin=chi4;
1628 type=4;
1629 }
1630 }
1631
1632 //
1633 }
1634
1635 if(type!=0 && chimin<100){
1636 m_isPPPiPi = true;
1637 m_ppPiPiNumber++;
1638 }
1639
1640 //end of tight cuts
1641
1642
1643 } //end of loose cuts
1644
1645 }//end of selecting pppipi
1646
1647 //select recoil events
1648 //if( m_selectRecoJpsi && (nGood==4 || nGood==6) && nCharge==0 ){
1649 if( (m_selectPsipRhoPi || m_selectPsipKstarK || m_selectPsipPPPiPi) && (nGood==4 || nGood==6) && nCharge==0 ){
1650 EvtRecTrackIterator pione, pitwo;
1651 RecMdcKalTrack *pitrk1;
1652 RecMdcKalTrack *pitrk2;
1653
1654 double bestmass=1.0;
1655 int pindex,nindex;
1656 vector<int> iJPsiP,iJPsiN;
1657 for(int ip=0; ip<iCP.size(); ip++){
1658 for(int in=0; in<iCN.size();in++){
1659 pione = evtRecTrkCol->begin() + iCP[ip];
1660 pitwo = evtRecTrkCol->begin() + iCN[in];
1661 pitrk1 = (*pione)->mdcKalTrack();
1662 pitrk2 = (*pitwo)->mdcKalTrack();
1663 Hep3Vector p1(Px(pitrk1), Py(pitrk1), Pz(pitrk1));
1664 Hep3Vector p2(Px(pitrk2), Py(pitrk2), Pz(pitrk2));
1665 double E1=sqrt( pow(P(pitrk1),2)+mpi*mpi);
1666 double E2=sqrt( pow(P(pitrk2),2)+mpi*mpi);
1667 double recomass = sqrt(pow(3.686-E1-E2,2)- (-(p1+p2)).mag2() );
1668 //if(fabs(recomass-3.686)< fabs(bestmass-3.686)){
1669 if(fabs(recomass-3.096)< fabs(bestmass-3.096)){
1670 bestmass=recomass;
1671 pindex=ip;
1672 nindex=in;
1673 }
1674 }
1675 }
1676
1677
1678 //soft pions
1679 pione = evtRecTrkCol->begin() + iCP[pindex];
1680 pitwo = evtRecTrkCol->begin() + iCN[nindex];
1681 pitrk1 = (*pione)->mdcKalTrack();
1682 pitrk2 = (*pitwo)->mdcKalTrack();
1683
1684
1685 //tracks from jpsi
1686 double jpsicharge=0;
1687 for(int ip=0; ip<iCP.size(); ip++){
1688 if(ip!=pindex){
1689 iJPsiP.push_back(iCP[ip]);
1690 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCP[ip];
1691 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
1692 jpsicharge+=trktmp->charge();
1693 }
1694 }
1695
1696 for(int in=0; in<iCN.size(); in++){
1697 if(in!=nindex){
1698 iJPsiN.push_back(iCN[in]);
1699 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCN[in];
1700 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
1701 jpsicharge+=trktmp->charge();
1702 }
1703 }
1704
1705
1706 HepLorentzVector ecms(0.034,0,0,m_ecm);
1707
1708 //jpsi to 2 charged track and 1 pi0
1709 if( (m_selectPsipRhoPi || m_selectPsipKstarK) && (bestmass>3.075 && bestmass<3.120) && nGood==4 && jpsicharge==0 && nPi0==1){
1710
1711 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1712 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiN[0];
1713
1714
1715 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1716
1717 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1718 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1719
1720 HepLorentzVector p4trk1;
1721 p4trk1.setPx(Px(trk1));
1722 p4trk1.setPy(Py(trk1));
1723 p4trk1.setPz(Pz(trk1));
1724 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
1725
1726 HepLorentzVector p4trk2;
1727 p4trk2.setPx(Px(trk2));
1728 p4trk2.setPy(Py(trk2));
1729 p4trk2.setPz(Pz(trk2));
1730 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
1731
1732
1733 HepLorentzVector p4trk3;
1734 p4trk3.setPx(Px(trk1));
1735 p4trk3.setPy(Py(trk1));
1736 p4trk3.setPz(Pz(trk1));
1737 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1738
1739 HepLorentzVector p4trk4;
1740 p4trk4.setPx(Px(trk2));
1741 p4trk4.setPy(Py(trk2));
1742 p4trk4.setPz(Pz(trk2));
1743 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1744
1745
1746 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1747 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1748 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1749 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1750 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1751 RecEmcShower *gShr3 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1752 RecEmcShower *gShr4 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1753
1754
1755 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1756 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1757 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
1758 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
1759 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
1760
1761
1762 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1763 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1764 double mkstarp = p4kstarp.m();
1765 double mkstarm = p4kstarm.m();
1766
1767 if( (p4total.m() > 3.0 && p4total.m() < 3.15) || (p4total2.m() > 3.0 && p4total2.m() < 3.15) ){
1768 //m_isRecoJpsi = true;
1769 //m_recoJpsiNumber++;
1770
1771
1772 //tighten cuts for rhopi and kstark
1773
1774
1775 WTrackParameter wvpiTrk1, wvpiTrk2;
1776 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1777 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1778
1779 WTrackParameter wvkTrk1, wvkTrk2;
1780 wvkTrk1 = WTrackParameter(mkaon, trk1->getZHelixK(), trk1->getZErrorK());//kaon
1781 wvkTrk2 = WTrackParameter(mkaon, trk2->getZHelixK(), trk2->getZErrorK());//kaon
1782
1783 //soft pions
1784 WTrackParameter wvpiTrk3, wvpiTrk4;
1785 wvpiTrk3 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1786 wvpiTrk4 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1787
1788 if(m_selectPsipRhoPi){
1789 KinematicFit * kmfit = KinematicFit::instance();
1790 kmfit->init();
1791 kmfit->AddTrack(0, wvpiTrk1);
1792 kmfit->AddTrack(1, wvpiTrk2);
1793 kmfit->AddTrack(2, 0.0, gShr1);
1794 kmfit->AddTrack(3, 0.0, gShr2);
1795 kmfit->AddTrack(4, wvpiTrk3);
1796 kmfit->AddTrack(5, wvpiTrk4);
1797 kmfit->AddFourMomentum(0, ecms);
1798
1799 bool oksq1 = kmfit->Fit();
1800 double chi1 = kmfit->chisq();
1801 //
1802
1803 if(oksq1 && chi1>0 && chi1<100){
1804 m_isPsipRhoPi = true;
1805 m_psipRhoPiNumber++;
1806 }
1807 }
1808 if(m_selectPsipKstarK){
1809 KinematicFit * kmfit2 = KinematicFit::instance();
1810 kmfit2->init();
1811 kmfit2->AddTrack(0, wvkTrk1);
1812 kmfit2->AddTrack(1, wvkTrk2);
1813 kmfit2->AddTrack(2, 0.0, gShr3);
1814 kmfit2->AddTrack(3, 0.0, gShr4);
1815 kmfit2->AddTrack(4, wvpiTrk3);
1816 kmfit2->AddTrack(5, wvpiTrk4);
1817 kmfit2->AddFourMomentum(0, ecms);
1818
1819
1820 bool oksq2 = kmfit2->Fit();
1821 double chi2 = kmfit2->chisq();
1822
1823 if(oksq2 && chi2>0 && chi2<200 &&
1824 ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1)
1825 || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) )){
1826 m_isPsipKstarK = true;
1827 m_psipKstarKNumber++;
1828 }
1829
1830 }
1831
1832
1833 }//end of loose cuts
1834
1835 }
1836
1837 } //recoil jpsi to 2tracks and 1 pi0
1838
1839
1840
1841 //jpsi to pppipi
1842 if(m_selectPsipPPPiPi && (bestmass>3.075 && bestmass<3.120) && nGood==6 && jpsicharge==0 && nPi0==0){
1843
1844
1845 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1846 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiP[1];
1847 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iJPsiN[0];
1848 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iJPsiN[1];
1849 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1850 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1851 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
1852 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
1853
1854 HepLorentzVector p4trkpp;
1855 HepLorentzVector p4trkpm;
1856 HepLorentzVector p4trkpip;
1857 HepLorentzVector p4trkpim;
1858
1859 //hypothesis 1
1860
1861 trk1->setPidType(RecMdcKalTrack::proton);
1862 p4trkpp.setPx(Px(trk1));
1863 p4trkpp.setPy(Py(trk1));
1864 p4trkpp.setPz(Pz(trk1));
1865 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1866
1867 trk3->setPidType(RecMdcKalTrack::proton);
1868 p4trkpm.setPx(Px(trk3));
1869 p4trkpm.setPy(Py(trk3));
1870 p4trkpm.setPz(Pz(trk3));
1871 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1872
1873
1874 trk2->setPidType(RecMdcKalTrack::pion);
1875 p4trkpip.setPx(Px(trk2));
1876 p4trkpip.setPy(Py(trk2));
1877 p4trkpip.setPz(Pz(trk2));
1878 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1879
1880 trk4->setPidType(RecMdcKalTrack::pion);
1881 p4trkpim.setPx(Px(trk4));
1882 p4trkpim.setPy(Py(trk4));
1883 p4trkpim.setPz(Pz(trk4));
1884 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1885
1886 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1887
1888
1889 //hypothesis 2
1890
1891 trk1->setPidType(RecMdcKalTrack::proton);
1892 p4trkpp.setPx(Px(trk1));
1893 p4trkpp.setPy(Py(trk1));
1894 p4trkpp.setPz(Pz(trk1));
1895 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1896
1897 trk4->setPidType(RecMdcKalTrack::proton);
1898 p4trkpm.setPx(Px(trk4));
1899 p4trkpm.setPy(Py(trk4));
1900 p4trkpm.setPz(Pz(trk4));
1901 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1902
1903 trk2->setPidType(RecMdcKalTrack::pion);
1904 p4trkpip.setPx(Px(trk2));
1905 p4trkpip.setPy(Py(trk2));
1906 p4trkpip.setPz(Pz(trk2));
1907 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1908
1909 trk3->setPidType(RecMdcKalTrack::pion);
1910 p4trkpim.setPx(Px(trk3));
1911 p4trkpim.setPy(Py(trk3));
1912 p4trkpim.setPz(Pz(trk3));
1913 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1914
1915 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1916
1917
1918 //hypothesis 3
1919
1920
1921 trk2->setPidType(RecMdcKalTrack::proton);
1922 p4trkpp.setPx(Px(trk2));
1923 p4trkpp.setPy(Py(trk2));
1924 p4trkpp.setPz(Pz(trk2));
1925 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1926
1927 trk3->setPidType(RecMdcKalTrack::proton);
1928 p4trkpm.setPx(Px(trk3));
1929 p4trkpm.setPy(Py(trk3));
1930 p4trkpm.setPz(Pz(trk3));
1931 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1932
1933 trk1->setPidType(RecMdcKalTrack::pion);
1934 p4trkpip.setPx(Px(trk1));
1935 p4trkpip.setPy(Py(trk1));
1936 p4trkpip.setPz(Pz(trk1));
1937 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1938
1939 trk4->setPidType(RecMdcKalTrack::pion);
1940 p4trkpim.setPx(Px(trk4));
1941 p4trkpim.setPy(Py(trk4));
1942 p4trkpim.setPz(Pz(trk4));
1943 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1944
1945 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1946
1947
1948 //hypothesis 4
1949
1950
1951 trk2->setPidType(RecMdcKalTrack::proton);
1952 p4trkpp.setPx(Px(trk2));
1953 p4trkpp.setPy(Py(trk2));
1954 p4trkpp.setPz(Pz(trk2));
1955 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1956
1957 trk4->setPidType(RecMdcKalTrack::proton);
1958 p4trkpm.setPx(Px(trk4));
1959 p4trkpm.setPy(Py(trk4));
1960 p4trkpm.setPz(Pz(trk4));
1961 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1962
1963 trk1->setPidType(RecMdcKalTrack::pion);
1964 p4trkpip.setPx(Px(trk1));
1965 p4trkpip.setPy(Py(trk1));
1966 p4trkpip.setPz(Pz(trk1));
1967 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1968
1969 trk3->setPidType(RecMdcKalTrack::pion);
1970 p4trkpim.setPx(Px(trk3));
1971 p4trkpim.setPy(Py(trk3));
1972 p4trkpim.setPz(Pz(trk3));
1973 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1974
1975 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1976
1977 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1978 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1979
1980
1981 //tighten cuts
1982 double chi1, chi2, chi3, chi4;
1983 int type=0;
1984 double chimin=-999;
1985
1986
1987 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 , wvpiTrk5, wvpiTrk6 ;
1988
1989 {
1990 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1991 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1992 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1993 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1994 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1995 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1996
1997
1998
1999 KinematicFit * kmfit = KinematicFit::instance();
2000 kmfit->init();
2001 kmfit->AddTrack(0, wvpiTrk1);
2002 kmfit->AddTrack(1, wvpiTrk2);
2003 kmfit->AddTrack(2, wvpiTrk3);
2004 kmfit->AddTrack(3, wvpiTrk4);
2005 kmfit->AddTrack(4, wvpiTrk5);
2006 kmfit->AddTrack(5, wvpiTrk6);
2007 kmfit->AddFourMomentum(0, ecms);
2008
2009 bool oksq1 = kmfit->Fit();
2010 chi1 = kmfit->chisq();
2011 if(oksq1){
2012 chimin=chi1;
2013 type=1;
2014 }
2015
2016 }
2017
2018
2019 {
2020 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
2021 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
2022 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
2023 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
2024 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
2025 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
2026
2027
2028
2029 KinematicFit * kmfit = KinematicFit::instance();
2030 kmfit->init();
2031 kmfit->AddTrack(0, wvpiTrk1);
2032 kmfit->AddTrack(1, wvpiTrk2);
2033 kmfit->AddTrack(2, wvpiTrk3);
2034 kmfit->AddTrack(3, wvpiTrk4);
2035 kmfit->AddTrack(4, wvpiTrk5);
2036 kmfit->AddTrack(5, wvpiTrk6);
2037 kmfit->AddFourMomentum(0, ecms);
2038
2039 bool oksq1 = kmfit->Fit();
2040 chi2 = kmfit->chisq();
2041 if(oksq1){
2042 if(type==0){
2043 chimin=chi2;
2044 type=2;
2045 }
2046 else if(chi2<chimin){
2047 chimin=chi2;
2048 type=2;
2049 }
2050
2051 }
2052 //
2053 }
2054
2055 {
2056 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
2057 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
2058 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
2059 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
2060 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
2061 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
2062
2063
2064
2065 KinematicFit * kmfit = KinematicFit::instance();
2066 kmfit->init();
2067 kmfit->AddTrack(0, wvpiTrk1);
2068 kmfit->AddTrack(1, wvpiTrk2);
2069 kmfit->AddTrack(2, wvpiTrk3);
2070 kmfit->AddTrack(3, wvpiTrk4);
2071 kmfit->AddTrack(4, wvpiTrk5);
2072 kmfit->AddTrack(5, wvpiTrk6);
2073 kmfit->AddFourMomentum(0, ecms);
2074
2075 bool oksq1 = kmfit->Fit();
2076 chi3 = kmfit->chisq();
2077 if(oksq1){
2078 if(type==0){
2079 chimin=chi3;
2080 type=3;
2081 }
2082 else if(chi3<chimin){
2083 chimin=chi3;
2084 type=3;
2085 }
2086 }
2087 //delete kmfit;
2088 }
2089
2090 {
2091 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
2092 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
2093 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
2094 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
2095 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
2096 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
2097
2098
2099 KinematicFit * kmfit = KinematicFit::instance();
2100 kmfit->init();
2101 kmfit->AddTrack(0, wvpiTrk1);
2102 kmfit->AddTrack(1, wvpiTrk2);
2103 kmfit->AddTrack(2, wvpiTrk3);
2104 kmfit->AddTrack(3, wvpiTrk4);
2105 kmfit->AddTrack(4, wvpiTrk5);
2106 kmfit->AddTrack(5, wvpiTrk6);
2107 kmfit->AddFourMomentum(0, ecms);
2108
2109 bool oksq1 = kmfit->Fit();
2110 chi4 = kmfit->chisq();
2111 if(oksq1){
2112 if(type==0){
2113 chimin=chi4;
2114 type=4;
2115 }
2116 else if(chi4<chimin){
2117 chimin=chi4;
2118 type=4;
2119 }
2120 }
2121
2122 }
2123
2124
2125 if(chimin>0 && chimin<200){
2126 m_isPsipPPPiPi = true;
2127 m_psipPPPiPiNumber++;
2128 }
2129
2130 }//end of loose cuts
2131
2132 }//end of selecting recol jpsi to pppipi
2133
2134
2135
2136
2137 }//end of recoil jpsi selection
2138
2139
2140
2141 //select psi'' events using dtaging
2142
2143 if(m_selectPsippCand){
2144
2145 SmartDataPtr<EvtRecDTagCol> evtRecDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
2146 if ( ! evtRecDTagCol ) {
2147 log << MSG::FATAL << "Could not find EvtRecDTagCol" << endreq;
2148 return StatusCode::FAILURE;
2149 }
2150 if(evtRecDTagCol->size()>0){
2151
2152 EvtRecDTagCol::iterator m_iterbegin=evtRecDTagCol->begin();
2153 EvtRecDTagCol::iterator m_iterend=evtRecDTagCol->end();
2154 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
2155
2156 if( ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2157 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0Pi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2158 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2159 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2160 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2161 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2162 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPi && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2163 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2164 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2165 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPiPi&& fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ) {
2166 m_isPsippCand = true;
2167 m_psippCandNumber++;
2168 break;
2169 }
2170
2171
2172 }//end of looping D modes
2173
2174
2175 }//end of non-zero dtag list
2176
2177 }//end of selecting psi'' events
2178
2179 // -------- Write to root file
2180
2181 if( m_selectRadBhabha && m_isRadBhabha ) m_subalg3->execute();
2182 if( m_selectGGEE && m_isGGEE ) m_subalg4->execute();
2183 if( m_selectGG4Pi && m_isGG4Pi ) m_subalg5->execute();
2184 if( m_selectRadBhabhaT && m_isRadBhabhaT ) m_subalg6->execute();
2185 if( m_selectRhoPi && m_isRhoPi ) m_subalg7->execute();
2186 if( m_selectKstarK && m_isKstarK ) m_subalg8->execute();
2187 if( m_selectPPPiPi && m_isPPPiPi ) m_subalg9->execute();
2188 if( m_selectRecoJpsi && m_isRecoJpsi ) m_subalg10->execute();
2189 if( m_selectBhabha && m_isBhabha ) m_subalg11->execute();
2190 if( m_selectDimu && m_isDimu ) m_subalg12->execute();
2191 if( m_selectCosmic && m_isCosmic ) m_subalg13->execute();
2192 if( m_selectGenProton && m_isGenProton ) m_subalg14->execute();
2193 if( m_selectPsipRhoPi && m_isPsipRhoPi ) m_subalg15->execute();
2194 if( m_selectPsipKstarK && m_isPsipKstarK ) m_subalg16->execute();
2195 if( m_selectPsipPPPiPi && m_isPsipPPPiPi ) m_subalg17->execute();
2196 if( m_selectPsippCand && m_isPsippCand ) m_subalg18->execute();
2197
2198
2199 //if(m_output) {
2200 // m_runnb = run;
2201 // m_evtnb = event;
2202 // m_esum = esum;
2203 // m_eemc = eemc;
2204 // m_etot = etot;
2205 // m_nCharge = nCharge;
2206 // m_nGood = nGood;
2207 // m_nGam = nGam;
2208 // m_ptot = ptot;
2209 // m_pp = pp;
2210 // m_pm = pm;
2211 // m_maxE = maxE;
2212 // m_secE = secE;
2213 // m_dThe = dthe;
2214 // m_dPhi = dphi;
2215 // m_mdcHit1 = mdcHit1;
2216 // m_mdcHit2 = mdcHit2;
2217 // m_tuple0->write();
2218 //}
2219
2220 return StatusCode::SUCCESS;
2221
2222}
2223
2224// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2225StatusCode CalibEventSelect::finalize() {
2226
2227 MsgStream log(msgSvc(), name());
2228 log << MSG::INFO << "in finalize()" << endmsg;
2229
2230 cout<<"Total events: "<<m_events<<endl;
2231
2232
2233 if(m_selectRadBhabha) {
2234 cout << "Selected rad bhabha: " << m_radBhabhaNumber << endl;
2235 }
2236
2237
2238 if(m_selectGGEE) {
2239 cout << "Selected ggee events: " << m_GGEENumber << endl;
2240 }
2241
2242 if(m_selectGG4Pi) {
2243 cout << "Selected gg4pi events: " << m_GG4PiNumber << endl;
2244 }
2245
2246 if(m_selectRadBhabhaT) {
2247 cout << "Selected rad bhabha tight: " << m_radBhabhaTNumber << endl;
2248 }
2249
2250 if(m_selectRhoPi) {
2251 cout << "Selected rhopi events: " << m_rhoPiNumber << endl;
2252 }
2253
2254 if(m_selectKstarK) {
2255 cout << "Selected kstark events: " << m_kstarKNumber << endl;
2256 }
2257
2258 if(m_selectPPPiPi) {
2259 cout << "Selected pppipi events: " << m_ppPiPiNumber << endl;
2260 }
2261
2262 if(m_selectRecoJpsi) {
2263 cout << "Selected recoil jsi events: " << m_recoJpsiNumber << endl;
2264 }
2265
2266
2267 if(m_selectBhabha) {
2268 cout << "Selected bhabha events: " << m_bhabhaNumber << endl;
2269 }
2270
2271
2272 if(m_selectDimu) {
2273 cout << "Selected dimu events: " << m_dimuNumber << endl;
2274 }
2275
2276
2277 if(m_selectCosmic) {
2278 cout << "Selected cosmic events: " << m_cosmicNumber << endl;
2279 }
2280
2281 if(m_selectGenProton) {
2282 cout << "Selected generic proton events: " << m_genProtonNumber << endl;
2283 }
2284
2285 if(m_selectPsipRhoPi) {
2286 cout << "Selected recoil rhopi events: " << m_psipRhoPiNumber << endl;
2287 }
2288
2289 if(m_selectPsipKstarK) {
2290 cout << "Selected recoil kstark events: " << m_psipKstarKNumber << endl;
2291 }
2292
2293 if(m_selectPsipPPPiPi) {
2294 cout << "Selected recoil pppipi events: " << m_psipPPPiPiNumber << endl;
2295 }
2296
2297 if(m_selectPsippCand) {
2298 cout << "Selected psi'' candi events: " << m_psippCandNumber << endl;
2299 }
2300
2301 return StatusCode::SUCCESS;
2302}
2303
2304bool CalibEventSelect::WhetherSector(double ph,double ph1,double ph2) {
2305 double phi1=min(ph1,ph2);
2306 double phi2=max(ph1,ph2);
2307 double delta=0.0610865; //3.5*3.1415926/180.
2308 if((phi2-phi1)<CLHEP::pi){
2309 phi1 -=delta;
2310 phi2 +=delta;
2311 if(phi1<0.) phi1 += CLHEP::twopi;
2312 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
2313 double tmp1=min(phi1,phi2);
2314 double tmp2=max(phi1,phi2);
2315 phi1=tmp1;
2316 phi2=tmp2;
2317 }
2318 else{
2319 phi1 +=delta;
2320 phi2 -=delta;
2321 }
2322 if((phi2-phi1)<CLHEP::pi){
2323 if(ph<=phi2&&ph>=phi1) return true;
2324 else return false;
2325 }
2326 else{
2327 if(ph>=phi2||ph<=phi1) return true;
2328 else return false;
2329 }
2330}
2331
2332
2333 //*************************************************************************
2334void CalibEventSelect::readbeamEfromDb(int runNo, double & beamE){
2335
2336 const char host[] = "202.122.37.69";
2337 const char user[] = "guest";
2338 const char passwd[] = "guestpass";
2339 const char db[] = "run";
2340 unsigned int port_num = 3306;
2341
2342
2343 MYSQL* mysql = mysql_init(NULL);
2344 mysql = mysql_real_connect(mysql, host, user, passwd, db, port_num,
2345 NULL, // socket
2346 0); // client_flag
2347
2348 if (mysql == NULL) {
2349 fprintf(stderr, "can not open database: RunInfo for run %d lum\n", runNo);
2350 }
2351
2352 char stmt[1024];
2353
2354 snprintf(stmt, 1024,
2355 "select BER_PRB, BPR_PRB "
2356 "from RunParams where run_number = %d", runNo);
2357 if (mysql_real_query(mysql, stmt, strlen(stmt))) {
2358 fprintf(stderr, "query error\n");
2359 }
2360
2361 MYSQL_RES* result_set = mysql_store_result(mysql);
2362 MYSQL_ROW row = mysql_fetch_row(result_set);
2363 if (!row) {
2364 fprintf(stderr, "cannot find data for RunNo %d\n", runNo);
2365 return ;
2366 }
2367
2368 double E_E=0, E_P=0;
2369 sscanf(row[0], "%lf", &E_E);
2370 sscanf(row[1], "%lf", &E_P);
2371
2372 beamE=(E_E+E_P)/2.0;
2373}
p1
double p1[4]
Definition BabayagaNLO.cxx:43
p2
double p2[4]
Definition BabayagaNLO.cxx:43
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
Px
double Px(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:58
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CalibEventSelect.cxx:43
Pz
double Pz(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:78
mproton
const double mproton
Definition CalibEventSelect.cxx:55
Vp4
std::vector< HepLorentzVector > Vp4
Definition CalibEventSelect.cxx:51
Phi
double Phi(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:99
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:87
Py
double Py(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:67
mpi
const double mpi
Definition CalibEventSelect.cxx:53
mkaon
const double mkaon
Definition CalibEventSelect.cxx:54
Vint
std::vector< int > Vint
Definition CalibEventSelect.cxx:50
delta
const double delta
Definition CheckRung/check/gainVsDoca.cxx:32
mpi0
const double mpi0
Definition DQAJpsi2PPbarAlg.cxx:52
runNo
int runNo
Definition DQA_TO_DB.cxx:12
mkaon
const double mkaon
Definition DSemilepAlg.cxx:52
phi2
Double_t phi2
Definition DataBase/tau_mode.c:7
etot
Double_t etot
Definition DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
time
Double_t time
Definition DataBase/tau_mode.c:7
ebeam
double ebeam
Definition EvtPhokharaDef.hh:17
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition EvtRecTrack.h:111
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:100
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:53
mpi
const double mpi
Definition Gam4pikp.cxx:47
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:52
MYSQL_RES
struct st_mysql_res MYSQL_RES
Definition MysqlConnection.h:12
MYSQL
struct st_mysql MYSQL
Definition MysqlConnection.h:11
mproton
const double mproton
Definition PipiJpsi.cxx:50
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
histoSvc
IHistogramSvc * histoSvc()
Definition ServiceAccessor.h:25
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
NULL
#define NULL
Definition VertexExtrapolate.h:14
CalibEventSelect::WhetherSector
bool WhetherSector(double ph, double ph1, double ph2)
Definition CalibEventSelect.cxx:2304
CalibEventSelect::readbeamEfromDb
void readbeamEfromDb(int runNo, double &beamE)
Definition CalibEventSelect.cxx:2334
DstEmcShower::theta
double theta() const
Definition DstEmcShower.h:38
DstEmcShower::phi
double phi() const
Definition DstEmcShower.h:39
DstEmcShower::time
double time() const
Definition DstEmcShower.h:50
DstEmcShower::energy
double energy() const
Definition DstEmcShower.h:45
DstMdcDedx::normPH
double normPH() const
Definition DstMdcDedx.h:67
DstMdcDedx::chiP
double chiP() const
Definition DstMdcDedx.h:63
DstMdcKalTrack::y
const double y() const
Definition DstMdcKalTrack.h:74
DstMdcKalTrack::theta
const double theta() const
Definition DstMdcKalTrack.h:70
DstMdcKalTrack::z
const double z() const
Definition DstMdcKalTrack.h:75
DstMdcKalTrack::tanl
const double tanl(void) const
Definition DstMdcKalTrack.h:95
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition DstMdcKalTrack.h:59
DstMdcKalTrack::r
const double r() const
Definition DstMdcKalTrack.h:81
DstMdcKalTrack::kappa
const double kappa(void) const
Definition DstMdcKalTrack.h:93
DstMdcKalTrack::charge
const int charge() const
Definition DstMdcKalTrack.h:65
DstMdcKalTrack::fi0
const double fi0(void) const
Definition DstMdcKalTrack.h:92
DstMdcKalTrack::x
const double x() const
Definition DstMdcKalTrack.h:73
DstMucTrack::depth
double depth() const
Definition DstMucTrack.h:45
EvtRecDTag::kD0toKPiPiPiPi0
@ kD0toKPiPiPiPi0
Definition EvtRecDTag.h:57
IVertexDbSvc::isVertexValid
virtual bool isVertexValid()=0
IVertexDbSvc::SigmaPrimaryVertex
virtual double * SigmaPrimaryVertex()=0
IVertexDbSvc::PrimaryVertex
virtual double * PrimaryVertex()=0
KinematicFit::instance
static KinematicFit * instance()
Definition KinematicFit.cxx:20
KinematicFit::chisq
double chisq() const
Definition KinematicFit.h:150
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition KinematicFit.cxx:354
RecMdcKalTrack::getZHelix
const HepVector & getZHelix() const
Definition RecMdcKalTrack.h:66
RecMdcKalTrack::getZHelixK
HepVector & getZHelixK()
Definition RecMdcKalTrack.h:91
RecMdcKalTrack::getZError
const HepSymMatrix & getZError() const
Definition RecMdcKalTrack.h:67
RecMdcKalTrack::getZErrorK
HepSymMatrix & getZErrorK()
Definition RecMdcKalTrack.h:92
TofHitStatus::is_cluster
bool is_cluster() const
Definition TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition TofHitStatus.cxx:134
TrackPool::AddTrack
void AddTrack(const int number, const double mass, const RecMdcTrack *trk)
Definition TrackPool.cxx:22
VFHelix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Analysis/VertexFit/VertexFit-00-03-05/VertexFit/Helix.h:186
VFHelix::a
const HepVector & a(void) const
returns helix parameters.
Definition Analysis/VertexFit/VertexFit-00-03-05/VertexFit/Helix.h:240
ecms
const double ecms
Definition inclkstar.cxx:42
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
ptrk
float ptrk
Definition prepare_data.cxx:28
p3
double double * p3
Definition qcdloop1.h:76
pi
const float pi
Definition vector3.h:133