BOSS 7.0.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 /////////////////////////////////////////////////////////////////////////////
120
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( 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 double time1=-99,depth1=-99;
898 double time2=-99,depth2=-99;
899 if( (*itp)->isTofTrackValid() ){
900 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
901 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
902
903 for(;iter_tof!=tofTrkCol.end();iter_tof++){
904 TofHitStatus* status =new TofHitStatus;
905 status->setStatus( (*iter_tof)->status() );
906 if(status->is_cluster()){
907 time1=(*iter_tof)->tof();
908 }
909 delete status;
910 }
911 }
912
913 if( (*itp)->isMucTrackValid() ){
914 RecMucTrack* mucTrk=(*itp)->mucTrack();
915 depth1=mucTrk->depth();
916 }
917
918 if( (*itn)->isTofTrackValid() ){
919 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
920 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
921
922 for(;iter_tof!=tofTrkCol.end();iter_tof++){
923 TofHitStatus* status =new TofHitStatus;
924 status->setStatus( (*iter_tof)->status() );
925 if(status->is_cluster()){
926 time2=(*iter_tof)->tof();
927 }
928 delete status;
929 }
930 }
931
932 if( (*itn)->isMucTrackValid() ){
933 RecMucTrack* mucTrk=(*itn)->mucTrack();
934 depth2=mucTrk->depth();
935 }
936
937
938 //dimu
939 //if( m_selectDimu && time1!=-99 && time2!=-99 && fabs(time1-time2)<5 && depth1>5 && depth2>5){ //tight, no endcap
940 // 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
941 // if( eemc<0.3*m_ecm && (fabs(pmax)>0.45*m_ecm || fabs(psec)>0.45*m_ecm) ){
942 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)
943 && (depth1>=3 || depth2>=3)){
944 m_isDimu=true;
945 m_dimuNumber++;
946 }
947 }//end of dimu
948
949
950
951 }//end of bhabha, dimu
952
953
954
955 if(m_selectCosmic && nCosmicGood == 2 && eemc/m_ecm<0.3){
956
957 EvtRecTrackIterator itp = evtRecTrkCol->begin() + iCosmicGood[0];
958 EvtRecTrackIterator itn = evtRecTrkCol->begin() + iCosmicGood[1];
959
960 double time1=-99,depth1=-99;
961 double time2=-99,depth2=-99;
962 if( (*itp)->isTofTrackValid() ){
963 SmartRefVector<RecTofTrack> tofTrkCol= (*itp)->tofTrack();
964 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
965
966 for(;iter_tof!=tofTrkCol.end();iter_tof++){
967 TofHitStatus* status =new TofHitStatus;
968 status->setStatus( (*iter_tof)->status() );
969 if(status->is_cluster()){
970 time1=(*iter_tof)->tof();
971 }
972 delete status;
973 }
974 }
975
976 if( (*itp)->isMucTrackValid() ){
977 RecMucTrack* mucTrk=(*itp)->mucTrack();
978 depth1=mucTrk->depth();
979 }
980
981 if( (*itn)->isTofTrackValid() ){
982 SmartRefVector<RecTofTrack> tofTrkCol= (*itn)->tofTrack();
983 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
984
985 for(;iter_tof!=tofTrkCol.end();iter_tof++){
986 TofHitStatus* status =new TofHitStatus;
987 status->setStatus( (*iter_tof)->status() );
988 if(status->is_cluster()){
989 time2=(*iter_tof)->tof();
990 }
991 delete status;
992 }
993 }
994
995 if( (*itn)->isMucTrackValid() ){
996 RecMucTrack* mucTrk=(*itn)->mucTrack();
997 depth2=mucTrk->depth();
998 }
999
1000 //cosmic
1001 //if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 && depth1>5 && depth2>5){
1002 if( m_selectCosmic && time1!=-99 && time2!=-99 && fabs(time1-time2)>5 ){
1003 m_isCosmic=true;
1004 m_cosmicNumber++;
1005 }
1006
1007
1008 }//end of cosmic
1009
1010
1011 //select generic protons
1012
1013 if(m_events%m_proton_scale ==0 ){
1014
1015 bool protoncand=false;
1016 for(int i=0; i<nGood; i++){
1017
1018 EvtRecTrackIterator iter = evtRecTrkCol->begin() + iGood[i];
1019 RecMdcKalTrack *mdcTrk = (*iter)->mdcKalTrack();
1020 mdcTrk->setPidType(RecMdcKalTrack::proton);
1021
1022 double pp = P(mdcTrk);
1023 double chiP=-99;
1024 if((*iter)->isMdcDedxValid()){
1025 RecMdcDedx* dedxTrk=(*iter)->mdcDedx();
1026 chiP=dedxTrk->chiP();
1027 }
1028
1029 if( fabs(pp)<0.6 && fabs(chiP)<5){
1030 protoncand=true;
1031 break;
1032 }
1033 }
1034
1035 if( protoncand ){
1036 m_isGenProton=true;
1037 m_genProtonNumber++;
1038 }
1039
1040 }//end of generic proton
1041
1042
1043 //fill monitoring histograms for rad bhabha
1044
1045
1046 if(m_printmonitor){
1047 h_nGood->fill(nGood);
1048 h_nCharge->fill(nCharge);
1049 h_pmaxobeam->fill(pmax/(m_ecm/2.0));
1050 h_eopmax->fill(eopmax);
1051 h_eopsec->fill(eopsec);
1052 h_deltap->fill(pmax-psec);
1053 h_esumoecm->fill(esum/m_ecm);
1054 h_egmax->fill(egmax);
1055 h_deltaphi->fill(phid);
1056 h_Pt->fill(Pt.mag());
1057 }
1058
1059
1060
1061 //select radbhabha
1062 if(nGood>1 && pmax/(m_ecm/2.0)>0.4 && eopmax>0.5 && esum/m_ecm>0.75 &&
1063 egmax>0.08 && fabs( fabs(phid)-CLHEP::pi)*180.0/CLHEP::pi>2.86 )
1064 {
1065 m_isRadBhabha=true;
1066 m_radBhabhaNumber++;
1067 }
1068 //select radbhabha tight
1069 if( m_isRadBhabha )
1070 {
1071 //prescale high momentum tracks
1072 if(nGood==2 && nCharge==0 && eemc/m_ecm>=0.7 && eopmax>=0.85 && eopmax<=1.15 && eopsec>=0.85 && eopsec<=1.15 ){
1073
1074 if( fabs(fabs(pmax)-m_ecm/2.0)<0.1 && fabs(fabs(psec)-m_ecm/2.0)<0.1 ){
1075 if(m_events%1000==0){
1076 m_isRadBhabhaT=true;
1077 m_radBhabhaTNumber++;
1078 }
1079 }
1080 else{
1081 m_isRadBhabhaT=true;
1082 m_radBhabhaTNumber++;
1083 }
1084
1085 }
1086
1087
1088
1089 }
1090
1091 //select ggee events, (exclude radee tight)
1092 //if(!m_isRadBhabhaT && nGood==2 && nCharge==0 && eopmax>=0.85 && eopsec>=0.85 && eemc/m_ecm<=0.8 && Pt.mag()<=0.2)
1093 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)
1094 {
1095 m_isGGEE=true;
1096 m_GGEENumber++;
1097 }
1098
1099 //select gg4pi events
1100 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)
1101 {
1102 m_isGG4Pi=true;
1103 m_GG4PiNumber++;
1104 }
1105
1106 //select rhopi/kstark events
1107 if( (m_selectRhoPi || m_selectKstarK) && nGood == 2 && nCharge==0 && nPi0 == 1){
1108
1109 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iGood[0];
1110 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iGood[1];
1111
1112
1113 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1114
1115 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1116 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1117
1118 HepLorentzVector p4trk1;
1119 p4trk1.setPx(Px(trk1));
1120 p4trk1.setPy(Py(trk1));
1121 p4trk1.setPz(Pz(trk1));
1122 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
1123
1124 HepLorentzVector p4trk2;
1125 p4trk2.setPx(Px(trk2));
1126 p4trk2.setPy(Py(trk2));
1127 p4trk2.setPz(Pz(trk2));
1128 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
1129
1130
1131 HepLorentzVector p4trk3;
1132 p4trk3.setPx(Px(trk1));
1133 p4trk3.setPy(Py(trk1));
1134 p4trk3.setPz(Pz(trk1));
1135 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1136
1137 HepLorentzVector p4trk4;
1138 p4trk4.setPx(Px(trk2));
1139 p4trk4.setPy(Py(trk2));
1140 p4trk4.setPz(Pz(trk2));
1141 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1142
1143
1144 //EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1145 itpi0 = recPi0Col->begin();
1146 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1147 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1148 HepLorentzVector p4pi0 = p4gam1+p4gam2;
1149
1150 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0; //kstark
1151 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0; // rhopi
1152
1153 double rhopimass=p4total2.m();
1154 double kstarkmass=p4total.m();
1155 if( (kstarkmass > 3.0 && kstarkmass < 3.15) || (rhopimass > 3.0 && rhopimass < 3.15) ){
1156
1157 //tight cuts
1158
1159 //remove bhabha background
1160 double eop1=0.0,eop2=0.0;
1161 if( (*itone)->isEmcShowerValid() ){
1162 RecEmcShower *emcTrk = (*itone)->emcShower();
1163 double etrk=emcTrk->energy();
1164 //cout<<"trk1 p="<<P(trk1)<<endl;
1165 if(P(trk1)!=0){
1166 eop1 = etrk/P(trk1);
1167 //if( fabs(eop1)> 0.8 ) return StatusCode::SUCCESS;
1168 }
1169 }
1170
1171 if( (*ittwo)->isEmcShowerValid() ){
1172 RecEmcShower *emcTrk = (*ittwo)->emcShower();
1173 double etrk=emcTrk->energy();
1174 //cout<<"trk2 p="<<P(trk2)<<endl;
1175 if(P(trk2)!=0){
1176 eop2 = etrk/P(trk2);
1177 //if( fabs(eop2)> 0.8 ) return StatusCode::SUCCESS;
1178 }
1179 }
1180
1181 if(eop1<0.8 && eop2<0.8){
1182
1183 if(rhopimass>3.0 && rhopimass<3.15){
1184 //kinematic fit
1185 //HepLorentzVector ecms(0.034,0,0,3.097);
1186 HepLorentzVector ecms(0.034,0,0,m_ecm);
1187
1188
1189 WTrackParameter wvpiTrk1, wvpiTrk2;
1190 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1191 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1192
1193 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1194 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1195 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1196 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1197
1198 KinematicFit * kmfit = KinematicFit::instance();
1199 kmfit->init();
1200 kmfit->AddTrack(0, wvpiTrk1);
1201 kmfit->AddTrack(1, wvpiTrk2);
1202 kmfit->AddTrack(2, 0.0, gShr1);
1203 kmfit->AddTrack(3, 0.0, gShr2);
1204 kmfit->AddFourMomentum(0, ecms);
1205
1206 bool oksq1 = kmfit->Fit();
1207 double chi1 = kmfit->chisq();
1208
1209 //
1210 if(oksq1 && chi1<=60){
1211 m_isRhoPi = true;
1212 m_rhoPiNumber++;
1213 }
1214 } //end of selecting rhopi
1215
1216
1217 if(kstarkmass>3.0 && kstarkmass<3.15){
1218
1219 //kstark resonances
1220 double mkstarp=0, mkstarm=0;
1221 if( trk1->charge() >0 ){
1222 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1223 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1224 mkstarp = p4kstarp.m();
1225 mkstarm = p4kstarm.m();
1226 }
1227 else{
1228 HepLorentzVector p4kstarm = p4trk1 + p4pi0;
1229 HepLorentzVector p4kstarp = p4trk2 + p4pi0;
1230 mkstarp = p4kstarp.m();
1231 mkstarm = p4kstarm.m();
1232 }
1233
1234 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 ) ){
1235 //kinematic fit
1236 //HepLorentzVector ecms(0.034,0,0,3.097);
1237 HepLorentzVector ecms(0.034,0,0,m_ecm);
1238
1239 WTrackParameter wvpiTrk1, wvpiTrk2;
1240 wvpiTrk1 = WTrackParameter(mkaon, trk1->getZHelix(), trk1->getZError());
1241 wvpiTrk2 = WTrackParameter(mkaon, trk2->getZHelix(), trk2->getZError());
1242
1243 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1244 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1245 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1246 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1247
1248 KinematicFit * kmfit = KinematicFit::instance();
1249 kmfit->init();
1250 kmfit->AddTrack(0, wvpiTrk1);
1251 kmfit->AddTrack(1, wvpiTrk2);
1252 kmfit->AddTrack(2, 0.0, gShr1);
1253 kmfit->AddTrack(3, 0.0, gShr2);
1254 kmfit->AddFourMomentum(0, ecms);
1255
1256 bool oksq1 = kmfit->Fit();
1257 double chi1 = kmfit->chisq();
1258 //
1259
1260 if(oksq1 && chi1<=60){
1261 m_isKstarK = true;
1262 m_kstarKNumber++;
1263 }
1264
1265 }
1266
1267 } //end of selecting kstark
1268
1269 }//end of non di-electron
1270
1271 //end of tight cuts
1272
1273 }
1274 }
1275
1276
1277
1278 } //end of selecting rhopi/kstark events
1279
1280 //select ppPiPi events
1281 if(m_selectPPPiPi && nGood ==4 && nCharge == 0 && nPi0==0){
1282
1283 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iCP[0];
1284 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iCP[1];
1285 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iCN[0];
1286 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iCN[1];
1287 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1288 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1289 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
1290 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
1291
1292 HepLorentzVector p4trkpp;
1293 HepLorentzVector p4trkpm;
1294 HepLorentzVector p4trkpip;
1295 HepLorentzVector p4trkpim;
1296
1297 //hypothesis 1
1298
1299 trk1->setPidType(RecMdcKalTrack::proton);
1300 trk2->setPidType(RecMdcKalTrack::pion);
1301 trk3->setPidType(RecMdcKalTrack::proton);
1302 trk4->setPidType(RecMdcKalTrack::pion);
1303
1304
1305 p4trkpp.setPx(Px(trk1));
1306 p4trkpp.setPy(Py(trk1));
1307 p4trkpp.setPz(Pz(trk1));
1308 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1309
1310 p4trkpm.setPx(Px(trk3));
1311 p4trkpm.setPy(Py(trk3));
1312 p4trkpm.setPz(Pz(trk3));
1313 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1314
1315
1316 p4trkpip.setPx(Px(trk2));
1317 p4trkpip.setPy(Py(trk2));
1318 p4trkpip.setPz(Pz(trk2));
1319 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1320
1321 p4trkpim.setPx(Px(trk4));
1322 p4trkpim.setPy(Py(trk4));
1323 p4trkpim.setPz(Pz(trk4));
1324 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1325
1326 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1327
1328
1329 //hypothesis 2
1330
1331 trk1->setPidType(RecMdcKalTrack::proton);
1332 trk2->setPidType(RecMdcKalTrack::pion);
1333 trk3->setPidType(RecMdcKalTrack::pion);
1334 trk4->setPidType(RecMdcKalTrack::proton);
1335
1336
1337 p4trkpp.setPx(Px(trk1));
1338 p4trkpp.setPy(Py(trk1));
1339 p4trkpp.setPz(Pz(trk1));
1340 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1341
1342 p4trkpm.setPx(Px(trk4));
1343 p4trkpm.setPy(Py(trk4));
1344 p4trkpm.setPz(Pz(trk4));
1345 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1346
1347
1348 p4trkpip.setPx(Px(trk2));
1349 p4trkpip.setPy(Py(trk2));
1350 p4trkpip.setPz(Pz(trk2));
1351 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1352
1353 p4trkpim.setPx(Px(trk3));
1354 p4trkpim.setPy(Py(trk3));
1355 p4trkpim.setPz(Pz(trk3));
1356 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1357
1358 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1359
1360
1361 //hypothesis 3
1362
1363 trk1->setPidType(RecMdcKalTrack::pion);
1364 trk2->setPidType(RecMdcKalTrack::proton);
1365 trk3->setPidType(RecMdcKalTrack::proton);
1366 trk4->setPidType(RecMdcKalTrack::pion);
1367
1368
1369 p4trkpp.setPx(Px(trk2));
1370 p4trkpp.setPy(Py(trk2));
1371 p4trkpp.setPz(Pz(trk2));
1372 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1373
1374 p4trkpm.setPx(Px(trk3));
1375 p4trkpm.setPy(Py(trk3));
1376 p4trkpm.setPz(Pz(trk3));
1377 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1378
1379
1380 p4trkpip.setPx(Px(trk1));
1381 p4trkpip.setPy(Py(trk1));
1382 p4trkpip.setPz(Pz(trk1));
1383 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1384
1385 p4trkpim.setPx(Px(trk4));
1386 p4trkpim.setPy(Py(trk4));
1387 p4trkpim.setPz(Pz(trk4));
1388 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1389
1390 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1391
1392
1393 //hypothesis 4
1394
1395 trk1->setPidType(RecMdcKalTrack::pion);
1396 trk2->setPidType(RecMdcKalTrack::proton);
1397 trk3->setPidType(RecMdcKalTrack::pion);
1398 trk4->setPidType(RecMdcKalTrack::proton);
1399
1400
1401 p4trkpp.setPx(Px(trk2));
1402 p4trkpp.setPy(Py(trk2));
1403 p4trkpp.setPz(Pz(trk2));
1404 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1405
1406 p4trkpm.setPx(Px(trk4));
1407 p4trkpm.setPy(Py(trk4));
1408 p4trkpm.setPz(Pz(trk4));
1409 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1410
1411
1412 p4trkpip.setPx(Px(trk1));
1413 p4trkpip.setPy(Py(trk1));
1414 p4trkpip.setPz(Pz(trk1));
1415 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1416
1417 p4trkpim.setPx(Px(trk3));
1418 p4trkpim.setPy(Py(trk3));
1419 p4trkpim.setPz(Pz(trk3));
1420 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1421
1422 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1423
1424
1425
1426
1427 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1428 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1429
1430 //tight cuts
1431
1432 //kinematic fits
1433 double chi1, chi2, chi3, chi4;
1434 int type=0;
1435 double chimin=-999;
1436 HepLorentzVector ecms(0.034,0,0,m_ecm);
1437
1438 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 ;
1439
1440 {
1441 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1442 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1443 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1444 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1445
1446
1447 KinematicFit * kmfit = KinematicFit::instance();
1448 kmfit->init();
1449 kmfit->AddTrack(0, wvpiTrk1);
1450 kmfit->AddTrack(1, wvpiTrk2);
1451 kmfit->AddTrack(2, wvpiTrk3);
1452 kmfit->AddTrack(3, wvpiTrk4);
1453 kmfit->AddFourMomentum(0, ecms);
1454
1455 bool oksq1 = kmfit->Fit();
1456 chi1 = kmfit->chisq();
1457 if(oksq1) {
1458 chimin=chi1;
1459 type=1;
1460 }
1461 //
1462 }
1463
1464
1465 {
1466 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1467 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1468 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1469 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1470
1471
1472 KinematicFit * kmfit = KinematicFit::instance();
1473 kmfit->init();
1474 kmfit->AddTrack(0, wvpiTrk1);
1475 kmfit->AddTrack(1, wvpiTrk2);
1476 kmfit->AddTrack(2, wvpiTrk3);
1477 kmfit->AddTrack(3, wvpiTrk4);
1478 kmfit->AddFourMomentum(0, ecms);
1479
1480 bool oksq1 = kmfit->Fit();
1481 chi2 = kmfit->chisq();
1482 if(oksq1){
1483 if(type==0){
1484 chimin=chi2;
1485 type=2;
1486 }
1487 else if(chi2<chimin){
1488 chimin=chi2;
1489 type=2;
1490 }
1491 }
1492 //
1493 }
1494
1495 {
1496 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1497 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
1498 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1499 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1500
1501 KinematicFit * kmfit = KinematicFit::instance();
1502 kmfit->init();
1503 kmfit->AddTrack(0, wvpiTrk1);
1504 kmfit->AddTrack(1, wvpiTrk2);
1505 kmfit->AddTrack(2, wvpiTrk3);
1506 kmfit->AddTrack(3, wvpiTrk4);
1507
1508 kmfit->AddFourMomentum(0, ecms);
1509
1510 bool oksq1 = kmfit->Fit();
1511 chi3 = kmfit->chisq();
1512 if(oksq1){
1513 if(type==0){
1514 chimin=chi3;
1515 type=3;
1516 }
1517 else if(chi3<chimin){
1518 chimin=chi3;
1519 type=3;
1520 }
1521 }
1522
1523 //
1524 }
1525
1526
1527 {
1528 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1529 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
1530 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1531 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1532
1533 KinematicFit * kmfit = KinematicFit::instance();
1534 kmfit->init();
1535 kmfit->AddTrack(0, wvpiTrk1);
1536 kmfit->AddTrack(1, wvpiTrk2);
1537 kmfit->AddTrack(2, wvpiTrk3);
1538 kmfit->AddTrack(3, wvpiTrk4);
1539
1540 kmfit->AddFourMomentum(0, ecms);
1541
1542 bool oksq1 = kmfit->Fit();
1543 chi4 = kmfit->chisq();
1544
1545 if(oksq1){
1546 if(type==0){
1547 chimin=chi4;
1548 type=4;
1549 }
1550 else if(chi4<chimin){
1551 chimin=chi4;
1552 type=4;
1553 }
1554 }
1555
1556 //
1557 }
1558
1559 if(type!=0 && chimin<100){
1560 m_isPPPiPi = true;
1561 m_ppPiPiNumber++;
1562 }
1563
1564 //end of tight cuts
1565
1566
1567 } //end of loose cuts
1568
1569 }//end of selecting pppipi
1570
1571 //select recoil events
1572 //if( m_selectRecoJpsi && (nGood==4 || nGood==6) && nCharge==0 ){
1573 if( (m_selectPsipRhoPi || m_selectPsipKstarK || m_selectPsipPPPiPi) && (nGood==4 || nGood==6) && nCharge==0 ){
1574 EvtRecTrackIterator pione, pitwo;
1575 RecMdcKalTrack *pitrk1;
1576 RecMdcKalTrack *pitrk2;
1577
1578 double bestmass=1.0;
1579 int pindex,nindex;
1580 vector<int> iJPsiP,iJPsiN;
1581 for(int ip=0; ip<iCP.size(); ip++){
1582 for(int in=0; in<iCN.size();in++){
1583 pione = evtRecTrkCol->begin() + iCP[ip];
1584 pitwo = evtRecTrkCol->begin() + iCN[in];
1585 pitrk1 = (*pione)->mdcKalTrack();
1586 pitrk2 = (*pitwo)->mdcKalTrack();
1587 Hep3Vector p1(Px(pitrk1), Py(pitrk1), Pz(pitrk1));
1588 Hep3Vector p2(Px(pitrk2), Py(pitrk2), Pz(pitrk2));
1589 double E1=sqrt( pow(P(pitrk1),2)+mpi*mpi);
1590 double E2=sqrt( pow(P(pitrk2),2)+mpi*mpi);
1591 double recomass = sqrt(pow(3.686-E1-E2,2)- (-(p1+p2)).mag2() );
1592 //if(fabs(recomass-3.686)< fabs(bestmass-3.686)){
1593 if(fabs(recomass-3.096)< fabs(bestmass-3.096)){
1594 bestmass=recomass;
1595 pindex=ip;
1596 nindex=in;
1597 }
1598 }
1599 }
1600
1601
1602 //soft pions
1603 pione = evtRecTrkCol->begin() + iCP[pindex];
1604 pitwo = evtRecTrkCol->begin() + iCN[nindex];
1605 pitrk1 = (*pione)->mdcKalTrack();
1606 pitrk2 = (*pitwo)->mdcKalTrack();
1607
1608
1609 //tracks from jpsi
1610 double jpsicharge=0;
1611 for(int ip=0; ip<iCP.size(); ip++){
1612 if(ip!=pindex){
1613 iJPsiP.push_back(iCP[ip]);
1614 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCP[ip];
1615 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
1616 jpsicharge+=trktmp->charge();
1617 }
1618 }
1619
1620 for(int in=0; in<iCN.size(); in++){
1621 if(in!=nindex){
1622 iJPsiN.push_back(iCN[in]);
1623 EvtRecTrackIterator itertmp=evtRecTrkCol->begin() + iCN[in];
1624 RecMdcKalTrack* trktmp=(*itertmp)->mdcKalTrack();
1625 jpsicharge+=trktmp->charge();
1626 }
1627 }
1628
1629
1630 HepLorentzVector ecms(0.034,0,0,m_ecm);
1631
1632 //jpsi to 2 charged track and 1 pi0
1633 if( (m_selectPsipRhoPi || m_selectPsipKstarK) && (bestmass>3.075 && bestmass<3.120) && nGood==4 && jpsicharge==0 && nPi0==1){
1634
1635 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1636 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiN[0];
1637
1638
1639 if( (*itone)->isMdcKalTrackValid() && (*ittwo)->isMdcKalTrackValid() ) {
1640
1641 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1642 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1643
1644 HepLorentzVector p4trk1;
1645 p4trk1.setPx(Px(trk1));
1646 p4trk1.setPy(Py(trk1));
1647 p4trk1.setPz(Pz(trk1));
1648 p4trk1.setE(sqrt( pow(P(trk1),2)+ mkaon*mkaon) );
1649
1650 HepLorentzVector p4trk2;
1651 p4trk2.setPx(Px(trk2));
1652 p4trk2.setPy(Py(trk2));
1653 p4trk2.setPz(Pz(trk2));
1654 p4trk2.setE(sqrt( pow(P(trk2),2)+ mkaon*mkaon) );
1655
1656
1657 HepLorentzVector p4trk3;
1658 p4trk3.setPx(Px(trk1));
1659 p4trk3.setPy(Py(trk1));
1660 p4trk3.setPz(Pz(trk1));
1661 p4trk3.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1662
1663 HepLorentzVector p4trk4;
1664 p4trk4.setPx(Px(trk2));
1665 p4trk4.setPy(Py(trk2));
1666 p4trk4.setPz(Pz(trk2));
1667 p4trk4.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1668
1669
1670 EvtRecPi0Col::iterator itpi0 = recPi0Col->begin();
1671 const EvtRecTrack* gTrk1 = (*itpi0)->hiEnGamma();
1672 const EvtRecTrack* gTrk2 = (*itpi0)->loEnGamma();
1673 RecEmcShower *gShr1 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1674 RecEmcShower *gShr2 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1675 RecEmcShower *gShr3 = const_cast<EvtRecTrack*>(gTrk1)->emcShower();
1676 RecEmcShower *gShr4 = const_cast<EvtRecTrack*>(gTrk2)->emcShower();
1677
1678
1679 HepLorentzVector p4gam1 = (*itpi0)->hiPfit();
1680 HepLorentzVector p4gam2 = (*itpi0)->loPfit();
1681 HepLorentzVector p4pi0 = p4gam1 + p4gam2;
1682 HepLorentzVector p4total = p4trk1 + p4trk2 + p4pi0;
1683 HepLorentzVector p4total2 = p4trk3 + p4trk4 + p4pi0;
1684
1685
1686 HepLorentzVector p4kstarp = p4trk1 + p4pi0;
1687 HepLorentzVector p4kstarm = p4trk2 + p4pi0;
1688 double mkstarp = p4kstarp.m();
1689 double mkstarm = p4kstarm.m();
1690
1691 if( (p4total.m() > 3.0 && p4total.m() < 3.15) || (p4total2.m() > 3.0 && p4total2.m() < 3.15) ){
1692 //m_isRecoJpsi = true;
1693 //m_recoJpsiNumber++;
1694
1695
1696 //tighten cuts for rhopi and kstark
1697
1698
1699 WTrackParameter wvpiTrk1, wvpiTrk2;
1700 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1701 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1702
1703 WTrackParameter wvkTrk1, wvkTrk2;
1704 wvkTrk1 = WTrackParameter(mkaon, trk1->getZHelixK(), trk1->getZErrorK());//kaon
1705 wvkTrk2 = WTrackParameter(mkaon, trk2->getZHelixK(), trk2->getZErrorK());//kaon
1706
1707 //soft pions
1708 WTrackParameter wvpiTrk3, wvpiTrk4;
1709 wvpiTrk3 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1710 wvpiTrk4 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1711
1712 if(m_selectPsipRhoPi){
1713 KinematicFit * kmfit = KinematicFit::instance();
1714 kmfit->init();
1715 kmfit->AddTrack(0, wvpiTrk1);
1716 kmfit->AddTrack(1, wvpiTrk2);
1717 kmfit->AddTrack(2, 0.0, gShr1);
1718 kmfit->AddTrack(3, 0.0, gShr2);
1719 kmfit->AddTrack(4, wvpiTrk3);
1720 kmfit->AddTrack(5, wvpiTrk4);
1721 kmfit->AddFourMomentum(0, ecms);
1722
1723 bool oksq1 = kmfit->Fit();
1724 double chi1 = kmfit->chisq();
1725 //
1726
1727 if(oksq1 && chi1>0 && chi1<100){
1728 m_isPsipRhoPi = true;
1729 m_psipRhoPiNumber++;
1730 }
1731 }
1732 if(m_selectPsipKstarK){
1733 KinematicFit * kmfit2 = KinematicFit::instance();
1734 kmfit2->init();
1735 kmfit2->AddTrack(0, wvkTrk1);
1736 kmfit2->AddTrack(1, wvkTrk2);
1737 kmfit2->AddTrack(2, 0.0, gShr3);
1738 kmfit2->AddTrack(3, 0.0, gShr4);
1739 kmfit2->AddTrack(4, wvpiTrk3);
1740 kmfit2->AddTrack(5, wvpiTrk4);
1741 kmfit2->AddFourMomentum(0, ecms);
1742
1743
1744 bool oksq2 = kmfit2->Fit();
1745 double chi2 = kmfit2->chisq();
1746
1747 if(oksq2 && chi2>0 && chi2<200 &&
1748 ( (fabs(mkstarp-0.89166)<0.1 && fabs(mkstarm-0.89166)>0.1)
1749 || (fabs(mkstarm-0.89166)<0.1 && fabs(mkstarp-0.89166)>0.1 ) )){
1750 m_isPsipKstarK = true;
1751 m_psipKstarKNumber++;
1752 }
1753
1754 }
1755
1756
1757 }//end of loose cuts
1758
1759 }
1760
1761 } //recoil jpsi to 2tracks and 1 pi0
1762
1763
1764
1765 //jpsi to pppipi
1766 if(m_selectPsipPPPiPi && (bestmass>3.075 && bestmass<3.120) && nGood==6 && jpsicharge==0 && nPi0==0){
1767
1768
1769 EvtRecTrackIterator itone = evtRecTrkCol->begin() + iJPsiP[0];
1770 EvtRecTrackIterator ittwo = evtRecTrkCol->begin() + iJPsiP[1];
1771 EvtRecTrackIterator itthr = evtRecTrkCol->begin() + iJPsiN[0];
1772 EvtRecTrackIterator itfor = evtRecTrkCol->begin() + iJPsiN[1];
1773 RecMdcKalTrack *trk1 = (*itone)->mdcKalTrack();
1774 RecMdcKalTrack *trk2 = (*ittwo)->mdcKalTrack();
1775 RecMdcKalTrack *trk3 = (*itthr)->mdcKalTrack();
1776 RecMdcKalTrack *trk4 = (*itfor)->mdcKalTrack();
1777
1778 HepLorentzVector p4trkpp;
1779 HepLorentzVector p4trkpm;
1780 HepLorentzVector p4trkpip;
1781 HepLorentzVector p4trkpim;
1782
1783 //hypothesis 1
1784
1785 trk1->setPidType(RecMdcKalTrack::proton);
1786 trk2->setPidType(RecMdcKalTrack::pion);
1787 trk3->setPidType(RecMdcKalTrack::proton);
1788 trk4->setPidType(RecMdcKalTrack::pion);
1789
1790
1791 p4trkpp.setPx(Px(trk1));
1792 p4trkpp.setPy(Py(trk1));
1793 p4trkpp.setPz(Pz(trk1));
1794 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1795
1796 p4trkpm.setPx(Px(trk3));
1797 p4trkpm.setPy(Py(trk3));
1798 p4trkpm.setPz(Pz(trk3));
1799 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1800
1801
1802 p4trkpip.setPx(Px(trk2));
1803 p4trkpip.setPy(Py(trk2));
1804 p4trkpip.setPz(Pz(trk2));
1805 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1806
1807 p4trkpim.setPx(Px(trk4));
1808 p4trkpim.setPy(Py(trk4));
1809 p4trkpim.setPz(Pz(trk4));
1810 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1811
1812 double jpsimass1= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1813
1814
1815 //hypothesis 2
1816
1817 trk1->setPidType(RecMdcKalTrack::proton);
1818 trk2->setPidType(RecMdcKalTrack::pion);
1819 trk3->setPidType(RecMdcKalTrack::pion);
1820 trk4->setPidType(RecMdcKalTrack::proton);
1821
1822
1823 p4trkpp.setPx(Px(trk1));
1824 p4trkpp.setPy(Py(trk1));
1825 p4trkpp.setPz(Pz(trk1));
1826 p4trkpp.setE(sqrt( pow(P(trk1),2)+ mproton*mproton) );
1827
1828 p4trkpm.setPx(Px(trk4));
1829 p4trkpm.setPy(Py(trk4));
1830 p4trkpm.setPz(Pz(trk4));
1831 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1832
1833
1834 p4trkpip.setPx(Px(trk2));
1835 p4trkpip.setPy(Py(trk2));
1836 p4trkpip.setPz(Pz(trk2));
1837 p4trkpip.setE(sqrt( pow(P(trk2),2)+ mpi*mpi) );
1838
1839 p4trkpim.setPx(Px(trk3));
1840 p4trkpim.setPy(Py(trk3));
1841 p4trkpim.setPz(Pz(trk3));
1842 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1843
1844 double jpsimass2= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1845
1846
1847 //hypothesis 3
1848
1849 trk1->setPidType(RecMdcKalTrack::pion);
1850 trk2->setPidType(RecMdcKalTrack::proton);
1851 trk3->setPidType(RecMdcKalTrack::proton);
1852 trk4->setPidType(RecMdcKalTrack::pion);
1853
1854
1855 p4trkpp.setPx(Px(trk2));
1856 p4trkpp.setPy(Py(trk2));
1857 p4trkpp.setPz(Pz(trk2));
1858 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1859
1860 p4trkpm.setPx(Px(trk3));
1861 p4trkpm.setPy(Py(trk3));
1862 p4trkpm.setPz(Pz(trk3));
1863 p4trkpm.setE(sqrt( pow(P(trk3),2)+ mproton*mproton) );
1864
1865
1866 p4trkpip.setPx(Px(trk1));
1867 p4trkpip.setPy(Py(trk1));
1868 p4trkpip.setPz(Pz(trk1));
1869 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1870
1871 p4trkpim.setPx(Px(trk4));
1872 p4trkpim.setPy(Py(trk4));
1873 p4trkpim.setPz(Pz(trk4));
1874 p4trkpim.setE(sqrt( pow(P(trk4),2)+ mpi*mpi) );
1875
1876 double jpsimass3= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1877
1878
1879 //hypothesis 4
1880
1881 trk1->setPidType(RecMdcKalTrack::pion);
1882 trk2->setPidType(RecMdcKalTrack::proton);
1883 trk3->setPidType(RecMdcKalTrack::pion);
1884 trk4->setPidType(RecMdcKalTrack::proton);
1885
1886
1887 p4trkpp.setPx(Px(trk2));
1888 p4trkpp.setPy(Py(trk2));
1889 p4trkpp.setPz(Pz(trk2));
1890 p4trkpp.setE(sqrt( pow(P(trk2),2)+ mproton*mproton) );
1891
1892 p4trkpm.setPx(Px(trk4));
1893 p4trkpm.setPy(Py(trk4));
1894 p4trkpm.setPz(Pz(trk4));
1895 p4trkpm.setE(sqrt( pow(P(trk4),2)+ mproton*mproton) );
1896
1897
1898 p4trkpip.setPx(Px(trk1));
1899 p4trkpip.setPy(Py(trk1));
1900 p4trkpip.setPz(Pz(trk1));
1901 p4trkpip.setE(sqrt( pow(P(trk1),2)+ mpi*mpi) );
1902
1903 p4trkpim.setPx(Px(trk3));
1904 p4trkpim.setPy(Py(trk3));
1905 p4trkpim.setPz(Pz(trk3));
1906 p4trkpim.setE(sqrt( pow(P(trk3),2)+ mpi*mpi) );
1907
1908 double jpsimass4= (p4trkpp+p4trkpm+p4trkpip+p4trkpim).m();
1909
1910 if( fabs(jpsimass1-3.075)<=0.075 || fabs(jpsimass2-3.075)<=0.075 ||
1911 fabs(jpsimass3-3.075)<=0.075 || fabs(jpsimass4-3.075)<=0.075){
1912
1913
1914 //tighten cuts
1915 double chi1, chi2, chi3, chi4;
1916 int type=0;
1917 double chimin=-999;
1918
1919
1920 WTrackParameter wvpiTrk1, wvpiTrk2, wvpiTrk3, wvpiTrk4 , wvpiTrk5, wvpiTrk6 ;
1921
1922 {
1923 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1924 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1925 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
1926 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
1927 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1928 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1929
1930
1931
1932 KinematicFit * kmfit = KinematicFit::instance();
1933 kmfit->init();
1934 kmfit->AddTrack(0, wvpiTrk1);
1935 kmfit->AddTrack(1, wvpiTrk2);
1936 kmfit->AddTrack(2, wvpiTrk3);
1937 kmfit->AddTrack(3, wvpiTrk4);
1938 kmfit->AddTrack(4, wvpiTrk5);
1939 kmfit->AddTrack(5, wvpiTrk6);
1940 kmfit->AddFourMomentum(0, ecms);
1941
1942 bool oksq1 = kmfit->Fit();
1943 chi1 = kmfit->chisq();
1944 if(oksq1){
1945 chimin=chi1;
1946 type=1;
1947 }
1948
1949 }
1950
1951
1952 {
1953 wvpiTrk1 = WTrackParameter(mproton, trk1->getZHelix(), trk1->getZError());
1954 wvpiTrk2 = WTrackParameter(mpi, trk2->getZHelix(), trk2->getZError());
1955 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1956 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1957 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1958 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1959
1960
1961
1962 KinematicFit * kmfit = KinematicFit::instance();
1963 kmfit->init();
1964 kmfit->AddTrack(0, wvpiTrk1);
1965 kmfit->AddTrack(1, wvpiTrk2);
1966 kmfit->AddTrack(2, wvpiTrk3);
1967 kmfit->AddTrack(3, wvpiTrk4);
1968 kmfit->AddTrack(4, wvpiTrk5);
1969 kmfit->AddTrack(5, wvpiTrk6);
1970 kmfit->AddFourMomentum(0, ecms);
1971
1972 bool oksq1 = kmfit->Fit();
1973 chi2 = kmfit->chisq();
1974 if(oksq1){
1975 if(type==0){
1976 chimin=chi2;
1977 type=2;
1978 }
1979 else if(chi2<chimin){
1980 chimin=chi2;
1981 type=2;
1982 }
1983
1984 }
1985 //
1986 }
1987
1988 {
1989 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
1990 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
1991 wvpiTrk3 = WTrackParameter(mpi, trk3->getZHelix(), trk3->getZError());
1992 wvpiTrk4 = WTrackParameter(mproton, trk4->getZHelix(), trk4->getZError());
1993 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
1994 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
1995
1996
1997
1998 KinematicFit * kmfit = KinematicFit::instance();
1999 kmfit->init();
2000 kmfit->AddTrack(0, wvpiTrk1);
2001 kmfit->AddTrack(1, wvpiTrk2);
2002 kmfit->AddTrack(2, wvpiTrk3);
2003 kmfit->AddTrack(3, wvpiTrk4);
2004 kmfit->AddTrack(4, wvpiTrk5);
2005 kmfit->AddTrack(5, wvpiTrk6);
2006 kmfit->AddFourMomentum(0, ecms);
2007
2008 bool oksq1 = kmfit->Fit();
2009 chi3 = kmfit->chisq();
2010 if(oksq1){
2011 if(type==0){
2012 chimin=chi3;
2013 type=3;
2014 }
2015 else if(chi3<chimin){
2016 chimin=chi3;
2017 type=3;
2018 }
2019 }
2020 //delete kmfit;
2021 }
2022
2023 {
2024 wvpiTrk1 = WTrackParameter(mpi, trk1->getZHelix(), trk1->getZError());
2025 wvpiTrk2 = WTrackParameter(mproton, trk2->getZHelix(), trk2->getZError());
2026 wvpiTrk3 = WTrackParameter(mproton, trk3->getZHelix(), trk3->getZError());
2027 wvpiTrk4 = WTrackParameter(mpi, trk4->getZHelix(), trk4->getZError());
2028 wvpiTrk5 = WTrackParameter(mpi, pitrk1->getZHelix(), pitrk1->getZError());
2029 wvpiTrk6 = WTrackParameter(mpi, pitrk2->getZHelix(), pitrk2->getZError());
2030
2031
2032 KinematicFit * kmfit = KinematicFit::instance();
2033 kmfit->init();
2034 kmfit->AddTrack(0, wvpiTrk1);
2035 kmfit->AddTrack(1, wvpiTrk2);
2036 kmfit->AddTrack(2, wvpiTrk3);
2037 kmfit->AddTrack(3, wvpiTrk4);
2038 kmfit->AddTrack(4, wvpiTrk5);
2039 kmfit->AddTrack(5, wvpiTrk6);
2040 kmfit->AddFourMomentum(0, ecms);
2041
2042 bool oksq1 = kmfit->Fit();
2043 chi4 = kmfit->chisq();
2044 if(oksq1){
2045 if(type==0){
2046 chimin=chi4;
2047 type=4;
2048 }
2049 else if(chi4<chimin){
2050 chimin=chi4;
2051 type=4;
2052 }
2053 }
2054
2055 }
2056
2057
2058 if(chimin>0 && chimin<200){
2059 m_isPsipPPPiPi = true;
2060 m_psipPPPiPiNumber++;
2061 }
2062
2063 }//end of loose cuts
2064
2065 }//end of selecting recol jpsi to pppipi
2066
2067
2068
2069
2070 }//end of recoil jpsi selection
2071
2072
2073
2074 //select psi'' events using dtaging
2075
2076 if(m_selectPsippCand){
2077
2078 SmartDataPtr<EvtRecDTagCol> evtRecDTagCol(eventSvc(), EventModel::EvtRec::EvtRecDTagCol);
2079 if ( ! evtRecDTagCol ) {
2080 log << MSG::FATAL << "Could not find EvtRecDTagCol" << endreq;
2081 return StatusCode::FAILURE;
2082 }
2083 if(evtRecDTagCol->size()>0){
2084
2085 EvtRecDTagCol::iterator m_iterbegin=evtRecDTagCol->begin();
2086 EvtRecDTagCol::iterator m_iterend=evtRecDTagCol->end();
2087 for(EvtRecDTagCol::iterator iter=m_iterbegin; iter != m_iterend; iter++){
2088
2089 if( ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2090 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPi0Pi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2091 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2092 ((*iter)->decayMode()==EvtRecDTag::kD0toKPiPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2093 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPi && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2094 ((*iter)->decayMode()==EvtRecDTag::kD0toKsPiPiPi0 && fabs((*iter)->mBC()-1.865)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2095 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPi && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ||
2096 ((*iter)->decayMode()==EvtRecDTag::kDptoKPiPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03)||
2097 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPi0 && fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.05 && (*iter)->deltaE()<0.03) ||
2098 ((*iter)->decayMode()==EvtRecDTag::kDptoKsPiPiPi&& fabs((*iter)->mBC()-1.87)<0.006 && (*iter)->deltaE()>-0.03 && (*iter)->deltaE()<0.03) ) {
2099 m_isPsippCand = true;
2100 m_psippCandNumber++;
2101 }
2102
2103
2104 }//end of looping D modes
2105
2106
2107 }//end of non-zero dtag list
2108
2109 }//end of selecting psi'' events
2110
2111 // -------- Write to root file
2112
2113 if( m_selectRadBhabha && m_isRadBhabha ) m_subalg3->execute();
2114 if( m_selectGGEE && m_isGGEE ) m_subalg4->execute();
2115 if( m_selectGG4Pi && m_isGG4Pi ) m_subalg5->execute();
2116 if( m_selectRadBhabhaT && m_isRadBhabhaT ) m_subalg6->execute();
2117 if( m_selectRhoPi && m_isRhoPi ) m_subalg7->execute();
2118 if( m_selectKstarK && m_isKstarK ) m_subalg8->execute();
2119 if( m_selectPPPiPi && m_isPPPiPi ) m_subalg9->execute();
2120 if( m_selectRecoJpsi && m_isRecoJpsi ) m_subalg10->execute();
2121 if( m_selectBhabha && m_isBhabha ) m_subalg11->execute();
2122 if( m_selectDimu && m_isDimu ) m_subalg12->execute();
2123 if( m_selectCosmic && m_isCosmic ) m_subalg13->execute();
2124 if( m_selectGenProton && m_isGenProton ) m_subalg14->execute();
2125 if( m_selectPsipRhoPi && m_isPsipRhoPi ) m_subalg15->execute();
2126 if( m_selectPsipKstarK && m_isPsipKstarK ) m_subalg16->execute();
2127 if( m_selectPsipPPPiPi && m_isPsipPPPiPi ) m_subalg17->execute();
2128 if( m_selectPsippCand && m_isPsippCand ) m_subalg18->execute();
2129
2130
2131 //if(m_output) {
2132 // m_runnb = run;
2133 // m_evtnb = event;
2134 // m_esum = esum;
2135 // m_eemc = eemc;
2136 // m_etot = etot;
2137 // m_nCharge = nCharge;
2138 // m_nGood = nGood;
2139 // m_nGam = nGam;
2140 // m_ptot = ptot;
2141 // m_pp = pp;
2142 // m_pm = pm;
2143 // m_maxE = maxE;
2144 // m_secE = secE;
2145 // m_dThe = dthe;
2146 // m_dPhi = dphi;
2147 // m_mdcHit1 = mdcHit1;
2148 // m_mdcHit2 = mdcHit2;
2149 // m_tuple0->write();
2150 //}
2151
2152 return StatusCode::SUCCESS;
2153
2154}
2155
2156// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2157StatusCode CalibEventSelect::finalize() {
2158
2159 MsgStream log(msgSvc(), name());
2160 log << MSG::INFO << "in finalize()" << endmsg;
2161
2162 cout<<"Total events: "<<m_events<<endl;
2163
2164
2165 if(m_selectRadBhabha) {
2166 cout << "Selected rad bhabha: " << m_radBhabhaNumber << endl;
2167 }
2168
2169
2170 if(m_selectGGEE) {
2171 cout << "Selected ggee events: " << m_GGEENumber << endl;
2172 }
2173
2174 if(m_selectGG4Pi) {
2175 cout << "Selected gg4pi events: " << m_GG4PiNumber << endl;
2176 }
2177
2178 if(m_selectRadBhabhaT) {
2179 cout << "Selected rad bhabha tight: " << m_radBhabhaTNumber << endl;
2180 }
2181
2182 if(m_selectRhoPi) {
2183 cout << "Selected rhopi events: " << m_rhoPiNumber << endl;
2184 }
2185
2186 if(m_selectKstarK) {
2187 cout << "Selected kstark events: " << m_kstarKNumber << endl;
2188 }
2189
2190 if(m_selectPPPiPi) {
2191 cout << "Selected pppipi events: " << m_ppPiPiNumber << endl;
2192 }
2193
2194 if(m_selectRecoJpsi) {
2195 cout << "Selected recoil jsi events: " << m_recoJpsiNumber << endl;
2196 }
2197
2198
2199 if(m_selectBhabha) {
2200 cout << "Selected bhabha events: " << m_bhabhaNumber << endl;
2201 }
2202
2203
2204 if(m_selectDimu) {
2205 cout << "Selected dimu events: " << m_dimuNumber << endl;
2206 }
2207
2208
2209 if(m_selectCosmic) {
2210 cout << "Selected cosmic events: " << m_cosmicNumber << endl;
2211 }
2212
2213 if(m_selectGenProton) {
2214 cout << "Selected generic proton events: " << m_genProtonNumber << endl;
2215 }
2216
2217 if(m_selectPsipRhoPi) {
2218 cout << "Selected recoil rhopi events: " << m_psipRhoPiNumber << endl;
2219 }
2220
2221 if(m_selectPsipKstarK) {
2222 cout << "Selected recoil kstark events: " << m_psipKstarKNumber << endl;
2223 }
2224
2225 if(m_selectPsipPPPiPi) {
2226 cout << "Selected recoil pppipi events: " << m_psipPPPiPiNumber << endl;
2227 }
2228
2229 if(m_selectPsippCand) {
2230 cout << "Selected psi'' candi events: " << m_psippCandNumber << endl;
2231 }
2232
2233 return StatusCode::SUCCESS;
2234}
2235
2236bool CalibEventSelect::WhetherSector(double ph,double ph1,double ph2) {
2237 double phi1=min(ph1,ph2);
2238 double phi2=max(ph1,ph2);
2239 double delta=0.0610865; //3.5*3.1415926/180.
2240 if((phi2-phi1)<CLHEP::pi){
2241 phi1 -=delta;
2242 phi2 +=delta;
2243 if(phi1<0.) phi1 += CLHEP::twopi;
2244 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
2245 double tmp1=min(phi1,phi2);
2246 double tmp2=max(phi1,phi2);
2247 phi1=tmp1;
2248 phi2=tmp2;
2249 }
2250 else{
2251 phi1 +=delta;
2252 phi2 -=delta;
2253 }
2254 if((phi2-phi1)<CLHEP::pi){
2255 if(ph<=phi2&&ph>=phi1) return true;
2256 else return false;
2257 }
2258 else{
2259 if(ph>=phi2||ph<=phi1) return true;
2260 else return false;
2261 }
2262}
2263
2264
2265 //*************************************************************************
2266void CalibEventSelect::readbeamEfromDb(int runNo, double & beamE){
2267
2268 const char host[] = "202.122.37.69";
2269 const char user[] = "guest";
2270 const char passwd[] = "guestpass";
2271 const char db[] = "run";
2272 unsigned int port_num = 3306;
2273
2274
2275 MYSQL* mysql = mysql_init(NULL);
2276 mysql = mysql_real_connect(mysql, host, user, passwd, db, port_num,
2277 NULL, // socket
2278 0); // client_flag
2279
2280 if (mysql == NULL) {
2281 fprintf(stderr, "can not open database: RunInfo for run %d lum\n", runNo);
2282 }
2283
2284 char stmt[1024];
2285
2286 snprintf(stmt, 1024,
2287 "select BER_PRB, BPR_PRB "
2288 "from RunParams where run_number = %d", runNo);
2289 if (mysql_real_query(mysql, stmt, strlen(stmt))) {
2290 fprintf(stderr, "query error\n");
2291 }
2292
2293 MYSQL_RES* result_set = mysql_store_result(mysql);
2294 MYSQL_ROW row = mysql_fetch_row(result_set);
2295 if (!row) {
2296 fprintf(stderr, "cannot find data for RunNo %d\n", runNo);
2297 return ;
2298 }
2299
2300 double E_E=0, E_P=0;
2301 sscanf(row[0], "%lf", &E_E);
2302 sscanf(row[1], "%lf", &E_P);
2303
2304 beamE=(E_E+E_P)/2.0;
2305}
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
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
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
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
mproton
const double mproton
Definition: PipiJpsi.cxx:50
CalibEventSelect::WhetherSector
bool WhetherSector(double ph, double ph1, double ph2)
Definition: CalibEventSelect.cxx:2236
CalibEventSelect::readbeamEfromDb
void readbeamEfromDb(int runNo, double &beamE)
Definition: CalibEventSelect.cxx:2266
CalibEventSelect::CalibEventSelect
CalibEventSelect(const std::string &name, ISvcLocator *pSvcLocator)
Definition: CalibEventSelect.cxx:121
CalibEventSelect::initialize
StatusCode initialize()
Definition: CalibEventSelect.cxx:168
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::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition: KinematicFit.cxx:354
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-02-77/VertexFit/Helix.h:186
VFHelix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Analysis/VertexFit/VertexFit-00-02-77/VertexFit/Helix.h:240
ecms
const double ecms
Definition: inclkstar.cxx:42