BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
EventPreSelect.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
8#include "EventModel/Event.h"
12
13#include "TMath.h"
14#include "GaudiKernel/INTupleSvc.h"
15#include "GaudiKernel/NTuple.h"
16#include "GaudiKernel/Bootstrap.h"
17#include "GaudiKernel/IHistogramSvc.h"
18#include "CLHEP/Vector/ThreeVector.h"
19#include "CLHEP/Vector/LorentzVector.h"
20#include "CLHEP/Vector/TwoVector.h"
21#include "EmcRawEvent/EmcDigi.h"
23#include "MdcRawEvent/MdcDigi.h"
24
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/ISvcLocator.h"
27using CLHEP::Hep3Vector;
28using CLHEP::Hep2Vector;
29using CLHEP::HepLorentzVector;
30#include "CLHEP/Geometry/Point3D.h"
31#ifndef ENABLE_BACKWARDS_COMPATIBILITY
33#endif
36
37#include <vector>
38
39typedef std::vector<int> Vint;
40typedef std::vector<HepLorentzVector> Vp4;
41
42const double mpi = 0.13957;
43
44int EventPreSelect::idmax[43]={40,44,48,56,64,72,80,80,76,76,
45 88,88,100,100,112,112,128,128,140,140,
46 160,160,160,160,176,176,176,176,208,208,
47 208,208,240,240,240,240,256,256,256,256,
48 288,288,288};
49/////////////////////////////////////////////////////////////////////////////
50
51EventPreSelect::EventPreSelect(const std::string& name, ISvcLocator* pSvcLocator) :
52 Algorithm(name, pSvcLocator) {
53 //Declare the properties
54
55 declareProperty("Output", m_output = false);
56 declareProperty("Ecm", m_ecm=3.686);
57 declareProperty("SelectBhabha", m_selectBhabha=false);
58 declareProperty("SelectDimu", m_selectDimu=false);
59 declareProperty("SelectHadron", m_selectHadron=false);
60 declareProperty("SelectDiphoton", m_selectDiphoton=false);
61 declareProperty("WriteDst", m_writeDst=false);
62 declareProperty("WriteRec", m_writeRec=false);
63 declareProperty("Vr0cut", m_vr0cut=1.0);
64 declareProperty("Vz0cut", m_vz0cut=5.0);
65 declareProperty("Pt0HighCut", m_pt0HighCut=5.0);
66 declareProperty("Pt0LowCut", m_pt0LowCut=0.05);
67 declareProperty("EnergyThreshold", m_energyThreshold=0.05);
68 declareProperty("GammaPhiCut", m_gammaPhiCut=20.0);
69 declareProperty("GammaThetaCut", m_gammaThetaCut=20.0);
70
71 declareProperty("BhabhaEmcECut", m_bhabhaEmcECut=0.7*m_ecm);
72 declareProperty("BhabhaMaxECut", m_bhabhaMaxECut=0.3*m_ecm);
73 declareProperty("BhabhaSecECut", m_bhabhaSecECut=0.1*m_ecm);
74 declareProperty("BhabhaDTheCut", m_bhabhaDTheCut=3.);
75 declareProperty("BhabhaDPhiCut1", m_bhabhaDPhiCut1=-25.);
76 declareProperty("BhabhaDPhiCut2", m_bhabhaDPhiCut2=-4.);
77 declareProperty("BhabhaDPhiCut3", m_bhabhaDPhiCut3=2.);
78 declareProperty("BhabhaDPhiCut4", m_bhabhaDPhiCut4=20.);
79 declareProperty("BhabhaMdcHitCutB", m_bhabhaMdcHitCutB=15);
80 declareProperty("BhabhaMdcHitCutE", m_bhabhaMdcHitCutE=5);
81
82 declareProperty("DimuEHighCut", m_dimuEHighCut=0.27*m_ecm);
83 declareProperty("DimuELowCut", m_dimuELowCut=0.027*m_ecm);
84 declareProperty("DimuDTheCut", m_dimuDTheCut=3.);
85 declareProperty("DimuDPhiCut", m_dimuDPhiCut=23.);
86
87 declareProperty("HadronChaECut", m_hadronChaECut=0.3*m_ecm);
88 declareProperty("HadronNeuECut", m_hadronNeuECut=0.3*m_ecm);
89
90 declareProperty("DiphotonEmcECut", m_diphotonEmcECut=0.7*m_ecm);
91 declareProperty("DiphotonSecECut", m_diphotonSecECut=0.3*m_ecm);
92 declareProperty("DiphotonDTheCut", m_diphotonDTheCut=3.);
93 declareProperty("DiphotonDPhiCut1", m_diphotonDPhiCut1=-4.);
94 declareProperty("DiphotonDPhiCut2", m_diphotonDPhiCut2=2.);
95}
96
97// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
99 MsgStream log(msgSvc(), name());
100
101 log << MSG::INFO << "in initialize()" << endmsg;
102
103 m_bhabhaEmcECut=0.7*m_ecm;
104 m_bhabhaMaxECut=0.3*m_ecm;
105 m_bhabhaSecECut=0.1*m_ecm;
106 m_dimuEHighCut=0.27*m_ecm;
107 m_dimuELowCut=0.027*m_ecm;
108 m_diphotonEmcECut=0.7*m_ecm;
109 m_diphotonSecECut=0.3*m_ecm;
110 m_hadronChaECut=0.3*m_ecm;
111 m_hadronNeuECut=0.3*m_ecm;
112
113
114 if (m_ecm==2.4) {
115 m_bhabhaEmcECut=0.7*m_ecm;
116 m_bhabhaMaxECut=0.15*m_ecm;
117 m_bhabhaSecECut=0.05*m_ecm;
118 m_diphotonEmcECut=0.7*m_ecm;
119 m_diphotonSecECut=0.2*m_ecm;
120
121 }
122
123
124
125 StatusCode sc;
126
127 log << MSG::INFO << "creating sub-algorithms...." << endreq;
128
129 if(m_selectBhabha) {
130 sc = createSubAlgorithm( "EventWriter", "SelectBarrelBhabha", m_subalg1);
131 if( sc.isFailure() ) {
132 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabhaBarrel" <<endreq;
133 return sc;
134 } else {
135 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabhaBarrel" <<endreq;
136 }
137 sc = createSubAlgorithm( "EventWriter", "SelectEndcapBhabha", m_subalg2);
138 if( sc.isFailure() ) {
139 log << MSG::ERROR << "Error creating Sub-Algorithm SelectBhabhaEndcap" <<endreq;
140 return sc;
141 } else {
142 log << MSG::INFO << "Success creating Sub-Algorithm SelectBhabhaEndcap" <<endreq;
143 }
144 }
145
146 if(m_selectDimu) {
147 m_dimuAlg = new DimuPreSelect;
148 sc = createSubAlgorithm( "EventWriter", "SelectBarrelDimu", m_subalg3);
149 if( sc.isFailure() ) {
150 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimuBarrel" <<endreq;
151 return sc;
152 } else {
153 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimuBarrel" <<endreq;
154 }
155 sc = createSubAlgorithm( "EventWriter", "SelectEndcapDimu", m_subalg4);
156 if( sc.isFailure() ) {
157 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDimuEndcap" <<endreq;
158 return sc;
159 } else {
160 log << MSG::INFO << "Success creating Sub-Algorithm SelectDimuEndap" <<endreq;
161 }
162 }
163
164 if(m_selectHadron) {
165 sc = createSubAlgorithm( "EventWriter", "SelectHadron", m_subalg5);
166 if( sc.isFailure() ) {
167 log << MSG::ERROR << "Error creating Sub-Algorithm SelectHadron" <<endreq;
168 return sc;
169 } else {
170 log << MSG::INFO << "Success creating Sub-Algorithm SelectHadron" <<endreq;
171 }
172 }
173
174 if(m_selectDiphoton) {
175 sc = createSubAlgorithm( "EventWriter", "SelectBarrelDiphoton", m_subalg6);
176 if( sc.isFailure() ) {
177 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDiphotonBarrel" <<endreq;
178 return sc;
179 } else {
180 log << MSG::INFO << "Success creating Sub-Algorithm SelectDiphotonBarrel" <<endreq;
181 }
182 sc = createSubAlgorithm( "EventWriter", "SelectEndcapDiphoton", m_subalg7);
183 if( sc.isFailure() ) {
184 log << MSG::ERROR << "Error creating Sub-Algorithm SelectDiphotonEndcap" <<endreq;
185 return sc;
186 } else {
187 log << MSG::INFO << "Success creating Sub-Algorithm SelectDiphotonEndcap" <<endreq;
188 }
189 }
190
191 if(m_writeDst) {
192 sc = createSubAlgorithm( "EventWriter", "WriteDst", m_subalg8);
193 if( sc.isFailure() ) {
194 log << MSG::ERROR << "Error creating Sub-Algorithm WriteDst" <<endreq;
195 return sc;
196 } else {
197 log << MSG::INFO << "Success creating Sub-Algorithm WriteDst" <<endreq;
198 }
199 }
200
201 if(m_writeRec) {
202 sc = createSubAlgorithm( "EventWriter", "WriteRec", m_subalg9);
203 if( sc.isFailure() ) {
204 log << MSG::ERROR << "Error creating Sub-Algorithm WriteRec" <<endreq;
205 return sc;
206 } else {
207 log << MSG::INFO << "Success creating Sub-Algorithm WriteRec" <<endreq;
208 }
209 }
210
211 if(m_output) {
212 StatusCode status;
213 NTuplePtr nt0(ntupleSvc(), "FILE1/hadron");
214 if ( nt0 ) m_tuple0 = nt0;
215 else {
216 m_tuple0 = ntupleSvc()->book ("FILE1/hadron", CLID_ColumnWiseTuple, "N-Tuple example");
217 if ( m_tuple0 ) {
218 status = m_tuple0->addItem ("esum", m_esum);
219 status = m_tuple0->addItem ("eemc", m_eemc);
220 status = m_tuple0->addItem ("etot", m_etot);
221 status = m_tuple0->addItem ("nGood", m_nGood);
222 status = m_tuple0->addItem ("nCharge", m_nCharge);
223 status = m_tuple0->addItem ("nGam", m_nGam);
224 status = m_tuple0->addItem ("ptot", m_ptot);
225 status = m_tuple0->addItem ("pp", m_pp);
226 status = m_tuple0->addItem ("pm", m_pm);
227 status = m_tuple0->addItem ("run", m_runnb);
228 status = m_tuple0->addItem ("event", m_evtnb);
229 status = m_tuple0->addItem ("maxE", m_maxE);
230 status = m_tuple0->addItem ("secE", m_secE);
231 status = m_tuple0->addItem ("dthe", m_dThe);
232 status = m_tuple0->addItem ("dphi", m_dPhi);
233 status = m_tuple0->addItem ("mdcHit1", m_mdcHit1);
234 status = m_tuple0->addItem ("mdcHit2", m_mdcHit2);
235 }
236 else {
237 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple0) << endmsg;
238 return StatusCode::FAILURE;
239 }
240 }
241
242 NTuplePtr nt1(ntupleSvc(), "FILE1/vxyz");
243 if ( nt1 ) m_tuple1 = nt1;
244 else {
245 m_tuple1 = ntupleSvc()->book ("FILE1/vxyz", CLID_ColumnWiseTuple, "ks N-Tuple example");
246 if ( m_tuple1 ) {
247 status = m_tuple1->addItem ("vx0", m_vx0);
248 status = m_tuple1->addItem ("vy0", m_vy0);
249 status = m_tuple1->addItem ("vz0", m_vz0);
250 status = m_tuple1->addItem ("vr0", m_vr0);
251 status = m_tuple1->addItem ("theta0", m_theta0);
252 status = m_tuple1->addItem ("p0", m_p0);
253 status = m_tuple1->addItem ("pt0", m_pt0);
254 }
255 else {
256 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
257 return StatusCode::FAILURE;
258 }
259 }
260
261 NTuplePtr nt2(ntupleSvc(), "FILE1/photon");
262 if ( nt2 ) m_tuple2 = nt2;
263 else {
264 m_tuple2 = ntupleSvc()->book ("FILE1/photon", CLID_ColumnWiseTuple, "ks N-Tuple example");
265 if ( m_tuple2 ) {
266 status = m_tuple2->addItem ("dthe", m_dthe);
267 status = m_tuple2->addItem ("dphi", m_dphi);
268 status = m_tuple2->addItem ("dang", m_dang);
269 status = m_tuple2->addItem ("eraw", m_eraw);
270 }
271 else {
272 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
273 return StatusCode::FAILURE;
274 }
275 }
276
277 if(m_selectDimu) {
278 NTuplePtr nt3(ntupleSvc(), "FILE1/dimu");
279 if ( nt3 ) m_tuple3 = nt3;
280 else {
281 m_tuple3 = ntupleSvc()->book ("FILE1/dimu", CLID_ColumnWiseTuple, "ks N-Tuple example");
282 if ( m_tuple3 ) {
283 m_dimuAlg->BookNtuple(m_tuple3);
284 }
285 else {
286 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple3) << endmsg;
287 return StatusCode::FAILURE;
288 }
289 }
290 }
291 }
292 //
293 //--------end of book--------
294 //
295
296 m_events=0;
297 m_barrelBhabhaNumber=0;
298 m_endcapBhabhaNumber=0;
299 m_barrelDimuNumber=0;
300 m_endcapDimuNumber=0;
301 m_hadronNumber=0;
302 m_barrelDiphotonNumber=0;
303 m_endcapDiphotonNumber=0;
304
305 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
306 return StatusCode::SUCCESS;
307
308}
309
310// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
312
313 //setFilterPassed(false);
314
315 MsgStream log(msgSvc(), name());
316 log << MSG::INFO << "in execute()" << endreq;
317
318 if( m_writeDst) m_subalg8->execute();
319 if( m_writeRec) m_subalg9->execute();
320
321 m_isBarrelBhabha = false;
322 m_isEndcapBhabha = false;
323 m_isBarrelDimu = false;
324 m_isEndcapDimu = false;
325 m_isHadron = false;
326 m_isBarrelDiphoton = false;
327 m_isEndcapDiphoton = false;
328
329 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
330 if(!eventHeader)
331 {
332 cout<<" eventHeader "<<endl;
333 return StatusCode::FAILURE;
334 }
335
336 int run=eventHeader->runNumber();
337 int event=eventHeader->eventNumber();
338
339 if(m_events%1000==0) cout<< m_events << " -------- run,event: "<<run<<","<<event<<endl;
340 m_events++;
341
342 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
343 if(!evtRecEvent ) {
344 cout<<" evtRecEvent "<<endl;
345 return StatusCode::FAILURE;
346 }
347
348 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
349 << evtRecEvent->totalCharged() << " , "
350 << evtRecEvent->totalNeutral() << " , "
351 << evtRecEvent->totalTracks() <<endreq;
352 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
353 if(!evtRecTrkCol){
354 cout<<" evtRecTrkCol "<<endl;
355 return StatusCode::FAILURE;
356 }
357
358 if(evtRecEvent->totalTracks()!=evtRecTrkCol->size()) return StatusCode::SUCCESS;
359
360 // -------- Good Charged Track Selection
361 Vint iGood;
362 iGood.clear();
363 int nCharge = 0;
364 for(int i = 0;i < evtRecEvent->totalCharged(); i++)
365 {
366 //if(i>=evtRecTrkCol->size()) break;
367 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
368 if(!(*itTrk)->isMdcTrackValid()) continue;
369 RecMdcTrack *mdcTrk =(*itTrk)->mdcTrack();
370 double vx0 = mdcTrk->x();
371 double vy0 = mdcTrk->y();
372 double vz0 = mdcTrk->z();
373 double vr0 = mdcTrk->r();
374 double theta0 = mdcTrk->theta();
375 double p0 = mdcTrk->p();
376 double pt0 = mdcTrk->pxy();
377
378 if(m_output) {
379 m_vx0 = vx0;
380 m_vy0 = vy0;
381 m_vz0 = vz0;
382 m_vr0 = vr0;
383 m_theta0 = theta0;
384 m_p0 = p0;
385 m_pt0 = pt0;
386 m_tuple1->write();
387 }
388
389 if(fabs(vz0) >= m_vz0cut) continue;
390 if(vr0 >= m_vr0cut) continue;
391 if(pt0 >= m_pt0HighCut) continue;
392 if(pt0 <= m_pt0LowCut) continue;
393
394 iGood.push_back((*itTrk)->trackId());
395 nCharge += mdcTrk->charge();
396 }
397 int nGood = iGood.size();
398
399 // -------- Good Photon Selection
400 Vint iGam;
401 iGam.clear();
402 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
403 //if(i>=evtRecTrkCol->size()) break;
404 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
405 if(!(*itTrk)->isEmcShowerValid()) continue;
406 RecEmcShower *emcTrk = (*itTrk)->emcShower();
407 double eraw = emcTrk->energy();
408 if(m_output) {
409 m_eraw = eraw;
410 m_tuple2->write();
411 }
412 if(eraw < m_energyThreshold) continue;
413 iGam.push_back((*itTrk)->trackId());
414 }
415 int nGam = iGam.size();
416
417 // -------- Assign 4-momentum to each charged track
418 Vint ipip, ipim;
419 ipip.clear();
420 ipim.clear();
421 Vp4 ppip, ppim;
422 ppip.clear();
423 ppim.clear();
424
425 //cout<<"charged track:"<<endl;
426 double echarge = 0.; //total energy of charged track
427 double ptot = 0.; //total momentum of charged track
428 double etot = 0.; //total energy in MDC and EMC
429 double eemc = 0.; //total energy in EMC
430 double pp = 0.;
431 double pm = 0.;
432
433 for(int i = 0; i < nGood; i++) {
434 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGood[i];
435
436 if((*itTrk)->isMdcTrackValid()) {
437
438 //RecMdcKalTrack* mdcKalTrk = (*itTrk)->mdcKalTrack();//After ParticleID, use RecMdcKalTrack substitute RecMdcTrack
439 //RecMdcKalTrack::setPidType(RecMdcKalTrack::pion);//PID can set to electron, muon, pion, kaon and proton;The default setting is pion
440 //If use this algorithm in reconstruction job, mdcKalTrk->p()=0!!! I don't know why.
441
442 RecMdcTrack* mdcTrk = (*itTrk)->mdcTrack();
443
444 ptot += mdcTrk->p();
445
446 HepLorentzVector ptrk;
447 ptrk.setPx(mdcTrk->px());
448 ptrk.setPy(mdcTrk->py());
449 ptrk.setPz(mdcTrk->pz());
450 double p3 = ptrk.mag();
451 ptrk.setE(sqrt(p3*p3+mpi*mpi));
452 ptrk = ptrk.boost(-0.011,0,0);//boost to cms
453
454 echarge += ptrk.e();
455 etot += ptrk.e();
456
457 if(mdcTrk->charge() >0) {
458 ppip.push_back(ptrk);
459 pp = ptrk.rho();
460 } else {
461 ppim.push_back(ptrk);
462 pm = ptrk.rho();
463 }
464
465 }
466
467 if((*itTrk)->isEmcShowerValid()) {
468
469 RecEmcShower* emcTrk = (*itTrk)->emcShower();
470 double eraw = emcTrk->energy();
471 double phi = emcTrk->phi();
472 double the = emcTrk->theta();
473
474 HepLorentzVector ptrk;
475 ptrk.setPx(eraw*sin(the)*cos(phi));
476 ptrk.setPy(eraw*sin(the)*sin(phi));
477 ptrk.setPz(eraw*cos(the));
478 ptrk.setE(eraw);
479 ptrk = ptrk.boost(-0.011,0,0);// boost to cms
480
481 eemc += ptrk.e();
482 etot += ptrk.e();
483
484 }
485
486 }
487
488
489 // -------- Assign 4-momentum to each photon
490 //cout<<"neutral:"<<endl;
491 Vp4 pGam;
492 pGam.clear();
493 double eneu=0; //total energy of neutral track
494 for(int i = 0; i < nGam; i++) {
495 EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + iGam[i];
496 RecEmcShower* emcTrk = (*itTrk)->emcShower();
497 double eraw = emcTrk->energy();
498 double phi = emcTrk->phi();
499 double the = emcTrk->theta();
500 HepLorentzVector ptrk;
501 ptrk.setPx(eraw*sin(the)*cos(phi));
502 ptrk.setPy(eraw*sin(the)*sin(phi));
503 ptrk.setPz(eraw*cos(the));
504 ptrk.setE(eraw);
505 ptrk = ptrk.boost(-0.011,0,0);// boost to cms
506 pGam.push_back(ptrk);
507 eneu += ptrk.e();
508 etot += ptrk.e();
509 eemc += ptrk.e();
510 }
511
512 double esum = echarge + eneu;
513
514 // -------- Use EMC shower information only
515 double maxE = 0.;
516 double secE = 0.;
517 double maxThe = 999.;
518 double maxPhi = 999.;
519 double secThe = 999.;
520 double secPhi = 999.;
521 int npart = 999.;
522 double dphi = 999.;
523 double dthe = 999.;
524 int mdcHit1 = 0.;
525 int mdcHit2 = 0.;
526
527 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
528 if (!aShowerCol) {
529 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
530 return( StatusCode::SUCCESS);
531 }
532
533 int ishower = 0;
534 RecEmcShowerCol::iterator iShowerCol;
535 for(iShowerCol=aShowerCol->begin();
536 iShowerCol!=aShowerCol->end();
537 iShowerCol++) {
538
539 if(ishower == 0) {
540 maxE = (*iShowerCol)->e5x5();
541 maxThe = (*iShowerCol)->theta();
542 maxPhi = (*iShowerCol)->phi();
543 npart = (*iShowerCol)->module();
544 } else if(ishower == 1) {
545 secE = (*iShowerCol)->e5x5();
546 secThe = (*iShowerCol)->theta();
547 secPhi = (*iShowerCol)->phi();
548 } else if(ishower == 2) {
549 break;
550 }
551
552 ishower++;
553
554 }
555
556 if(aShowerCol->size() >= 2) {
557
558 dphi = (fabs(maxPhi-secPhi)-CLHEP::pi)*180./CLHEP::pi;
559 dthe = (fabs(maxThe+secThe)-CLHEP::pi)*180./CLHEP::pi;
560
561 double phi1 = maxPhi<0 ? maxPhi+CLHEP::twopi : maxPhi;
562 double phi2 = secPhi<0 ? secPhi+CLHEP::twopi : secPhi;
563
564 //Define sector (phi11,phi12) and (phi21,phi22)
565 double phi11=min(phi1,phi2);
566 double phi22=max(phi1,phi2);
567 double phi12=(phi11+phi22-CLHEP::pi)*0.5;
568 double phi21=(phi11+phi22+CLHEP::pi)*0.5;
569 if(phi12<0.) phi12 += CLHEP::twopi;
570 if(phi21>CLHEP::twopi) phi21 -= CLHEP::twopi;
571
572 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(), "/Event/Digi/MdcDigiCol");
573 if (!mdcDigiCol) {
574 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
575 return StatusCode::FAILURE;
576 }
577 int hitnum = mdcDigiCol->size();
578 for (int i = 0;i< hitnum;i++ ) {
579 MdcDigi* digi= dynamic_cast<MdcDigi*>(mdcDigiCol->containedObject(i));
580 double time = digi->getTimeChannel();
581 double charge = digi->getChargeChannel();
582 if (time == 0x7FFFFFFF || charge == 0x7FFFFFFF) continue;
583 Identifier id(digi->identify());
584 unsigned int iphi=MdcID::wire(id);
585 unsigned int ilayer=MdcID::layer(id);
586 if(ilayer>=43)
587 log << MSG::ERROR << "MDC(" << ilayer <<","<<iphi <<")"<<endreq;
588 double phi=CLHEP::twopi*iphi/idmax[ilayer];
589 if(WhetherSector(phi,phi11,phi12)) mdcHit1++;
590 else if(WhetherSector(phi,phi21,phi22)) mdcHit2++;
591 //cout<<"phi="<<phi<<",phi11="<<phi11<<",phi12="<<phi12
592 //<<",phi21="<<phi21<<",phi22="<<phi22<<endl;
593 }
594 }
595
596 //If it's bhabha, return.
597 //if(nGood==2 && pp+pm>3.4) return StatusCode::SUCCESS;
598
599 // -------- Select each event
600 // Select bhabha
601 if( eemc>m_bhabhaEmcECut && maxE>m_bhabhaMaxECut && secE>m_bhabhaSecECut
602 && abs(dthe)<m_bhabhaDTheCut &&
603 ( (dphi>m_bhabhaDPhiCut1&&dphi<m_bhabhaDPhiCut2) ||
604 (dphi>m_bhabhaDPhiCut3&&dphi<m_bhabhaDPhiCut4) ) ) {
605 if( npart==1 && mdcHit1>m_bhabhaMdcHitCutB && mdcHit2>m_bhabhaMdcHitCutB ) {
606 m_isBarrelBhabha = true;
607 m_barrelBhabhaNumber++;
608 } else if( ( npart==0 || npart==2 )
609 && mdcHit1>m_bhabhaMdcHitCutE && mdcHit2>m_bhabhaMdcHitCutE ) {
610 m_isEndcapBhabha = true;
611 m_endcapBhabhaNumber++;
612 }
613 }
614
615 // Select dimu
616 /*if( maxE<m_dimuEHighCut && maxE>m_dimuELowCut
617 && secE<m_dimuEHighCut && secE>m_dimuELowCut
618 && fabs(dthe)<m_dimuDTheCut && fabs(dphi)<m_dimuDPhiCut ) {
619 if( npart==1 ) {
620 m_isBarrelDimu = true;
621 m_barrelDimuNumber++;
622 } else if( npart==0 || npart==2 ) {
623 m_isEndcapDimu = true;
624 m_endcapDimuNumber++;
625 }
626 }*/
627 if(m_selectDimu) {
628 if( m_dimuAlg->IsDimu() == 1 ) {
629 m_isBarrelDimu = true;
630 m_barrelDimuNumber++;
631 } else if(m_dimuAlg->IsDimu() == 2 ) {
632 m_isEndcapDimu = true;
633 m_endcapDimuNumber++;
634 }
635 }
636
637 // Select hadron
638 if( (nGood>=1 && esum>m_hadronChaECut)
639 || (nGood==0 && esum>m_hadronNeuECut) ) {
640 m_isHadron = true;
641 m_hadronNumber++;
642 }
643
644 // Select diphoton
645 if( nGood==0 && eemc>m_diphotonEmcECut && secE>m_diphotonSecECut
646 && fabs(dthe)<m_diphotonDTheCut
647 && dphi>m_diphotonDPhiCut1 && dphi<m_diphotonDPhiCut2 ) {
648 if( npart==1 ) {
649 m_isBarrelDiphoton = true;
650 m_barrelDiphotonNumber++;
651 } else if( npart==0 || npart==2 ) {
652 m_isEndcapDiphoton = true;
653 m_endcapDiphotonNumber++;
654 }
655 }
656
657 // -------- Write to root file
658 if( m_selectBhabha && m_isBarrelBhabha ) m_subalg1->execute();
659 if( m_selectBhabha && m_isEndcapBhabha ) m_subalg2->execute();
660 if( m_selectDimu && m_isBarrelDimu ) m_subalg3->execute();
661 if( m_selectDimu && m_isEndcapDimu ) m_subalg4->execute();
662 if( m_selectHadron && m_isHadron ) m_subalg5->execute();
663 if( m_selectDiphoton && m_isBarrelDiphoton ) m_subalg6->execute();
664 if( m_selectDiphoton && m_isEndcapDiphoton ) m_subalg7->execute();
665
666
667 if(m_output) {
668 m_runnb = run;
669 m_evtnb = event;
670 m_esum = esum;
671 m_eemc = eemc;
672 m_etot = etot;
673 m_nCharge = nCharge;
674 m_nGood = nGood;
675 m_nGam = nGam;
676 m_ptot = ptot;
677 m_pp = pp;
678 m_pm = pm;
679 m_maxE = maxE;
680 m_secE = secE;
681 m_dThe = dthe;
682 m_dPhi = dphi;
683 m_mdcHit1 = mdcHit1;
684 m_mdcHit2 = mdcHit2;
685 m_tuple0->write();
686 }
687
688 return StatusCode::SUCCESS;
689
690}
691
692// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
694
695 MsgStream log(msgSvc(), name());
696 log << MSG::INFO << "in finalize()" << endmsg;
697
698 if(m_selectDimu&&m_output) {
699 m_dimuAlg->Print();
700 }
701
702 cout<<"Total events: "<<m_events<<endl;
703
704 if(m_selectBhabha) {
705 cout << "Selected barrel bhabha: " << m_barrelBhabhaNumber << endl;
706 cout << "Selected endcap bhabha: " << m_endcapBhabhaNumber << endl;
707 }
708
709 if(m_selectDimu) {
710 delete m_dimuAlg;
711 cout << "Selected barrel dimu: " << m_barrelDimuNumber << endl;
712 cout << "Selected endcap dimu: " << m_endcapDimuNumber << endl;
713 }
714
715 if(m_selectHadron) {
716 cout << "Selected hadron: " << m_hadronNumber << endl;
717 }
718
719 if(m_selectDiphoton) {
720 cout << "Selected barrel diphoton: " << m_barrelDiphotonNumber << endl;
721 cout << "Selected endcap diphoton: " << m_endcapDiphotonNumber << endl;
722 }
723
724 return StatusCode::SUCCESS;
725}
726
727bool EventPreSelect::WhetherSector(double ph,double ph1,double ph2) {
728 double phi1=min(ph1,ph2);
729 double phi2=max(ph1,ph2);
730 double delta=0.0610865; //3.5*3.1415926/180.
731 if((phi2-phi1)<CLHEP::pi){
732 phi1 -=delta;
733 phi2 +=delta;
734 if(phi1<0.) phi1 += CLHEP::twopi;
735 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
736 double tmp1=min(phi1,phi2);
737 double tmp2=max(phi1,phi2);
738 phi1=tmp1;
739 phi2=tmp2;
740 }
741 else{
742 phi1 +=delta;
743 phi2 -=delta;
744 }
745 if((phi2-phi1)<CLHEP::pi){
746 if(ph<=phi2&&ph>=phi1) return true;
747 else return false;
748 }
749 else{
750 if(ph>=phi2||ph<=phi1) return true;
751 else return false;
752 }
753}
754
755
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
const double delta
HepGeom::Point3D< double > HepPoint3D
std::vector< HepLorentzVector > Vp4
const double mpi
std::vector< int > Vint
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition: EvtRecTrack.h:111
std::vector< HepLorentzVector > Vp4
Definition: Gam4pikp.cxx:53
const double mpi
Definition: Gam4pikp.cxx:47
std::vector< int > Vint
Definition: Gam4pikp.cxx:52
INTupleSvc * ntupleSvc()
IMessageSvc * msgSvc()
double theta() const
Definition: DstEmcShower.h:38
double phi() const
Definition: DstEmcShower.h:39
double energy() const
Definition: DstEmcShower.h:45
const double theta() const
Definition: DstMdcTrack.h:59
const double r() const
Definition: DstMdcTrack.h:64
const double py() const
Definition: DstMdcTrack.h:56
const int charge() const
Definition: DstMdcTrack.h:53
const double px() const
Definition: DstMdcTrack.h:55
const double pz() const
Definition: DstMdcTrack.h:57
const double pxy() const
Definition: DstMdcTrack.h:54
const double z() const
Definition: DstMdcTrack.h:63
const double p() const
Definition: DstMdcTrack.h:58
const double y() const
Definition: DstMdcTrack.h:62
const double x() const
Definition: DstMdcTrack.h:61
StatusCode execute()
EventPreSelect(const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize()
bool WhetherSector(double ph, double ph1, double ph2)
StatusCode finalize()
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
virtual Identifier identify() const
Definition: RawData.cxx:15
unsigned int getChargeChannel() const
Definition: RawData.cxx:45
unsigned int getTimeChannel() const
Definition: RawData.cxx:40
_EXTERN_ std::string EvtRecEvent
Definition: EventModel.h:116
_EXTERN_ std::string EvtRecTrackCol
Definition: EventModel.h:117
float charge
float ptrk