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