BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
BhabhaPreSelect Class Reference

#include <BhabhaPreSelect.h>

+ Inheritance diagram for BhabhaPreSelect:

Public Member Functions

 BhabhaPreSelect (const std::string &name, ISvcLocator *pSvcLocator)
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
bool WhetherSector (double, double=0., double=CLHEP::twopi)
 
int selectedType () const
 
int selectedTrkID1 () const
 
int selectedTrkID2 () const
 

Detailed Description

Definition at line 11 of file BhabhaPreSelect.h.

Constructor & Destructor Documentation

◆ BhabhaPreSelect()

BhabhaPreSelect::BhabhaPreSelect ( const std::string & name,
ISvcLocator * pSvcLocator )

Definition at line 48 of file BhabhaPreSelect.cxx.

48 :
49 Algorithm(name, pSvcLocator) {
50
51 //Declare the properties
52
53 declareProperty ("BarrelOrEndcap", m_BarrelOrEndcap = 1);
54 declareProperty ("Output", m_output = false);
55 //cout<<" BarrelOrEndcap "<<m_BarrelOrEndcap<<endl;
56
57 }

Member Function Documentation

◆ execute()

StatusCode BhabhaPreSelect::execute ( )

Definition at line 129 of file BhabhaPreSelect.cxx.

129 {
130
131 MsgStream log(msgSvc(), name());
132 log << MSG::INFO << "in execute()" << endreq;
133
134 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
135 //if(!eventHeader)
136 //{
137 // cout<<" eventHeader "<<endl;
138 // return StatusCode::FAILURE;
139 //}
140
141 int run=eventHeader->runNumber();
142 int event=eventHeader->eventNumber();
143 //cout<<" event "<<event<<endl;
144 if(event%1000==0) cout<<" run,event: "<<run<<","<<event<<endl;
145
146 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
147 if(!evtRecEvent ) {
148 cout<<" evtRecEvent "<<endl;
149 return StatusCode::FAILURE;
150 }
151
152 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
153 << evtRecEvent->totalCharged() << " , "
154 << evtRecEvent->totalNeutral() << " , "
155 << evtRecEvent->totalTracks() <<endreq;
156 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
157 if(!evtRecTrkCol){
158 cout<<" evtRecTrkCol "<<endl;
159 return StatusCode::FAILURE;
160 }
161
162 m_events++;
163 setFilterPassed(false);
164
165 Vint iGood;
166 iGood.clear();
167
168 double ene5x5,theta,phi,eseed;
169 double showerX,showerY,showerZ;
170 long int thetaIndex,phiIndex;
171
172 RecEmcID showerId;
173 unsigned int npart;
174
175 Vint ipip, ipim;
176 ipip.clear();
177 ipim.clear();
178 Vp4 ppip, ppim;
179 ppip.clear();
180 ppim.clear();
181
182 vector<RecEmcShower* > GoodShowers;
183 GoodShowers.clear();
184 double etot=0;
185 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
186 if(i>=evtRecTrkCol->size()) break;
187 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
188 if((*itTrk)->isEmcShowerValid()) {
189 RecEmcShower *theShower = (*itTrk)->emcShower();
190 etot+=theShower->e5x5();
191 showerId = theShower->getShowerId();
192 npart = EmcID::barrel_ec(showerId);
193 ene5x5=theShower->e5x5();
194 thetaIndex=EmcID::theta_module(showerId);
195 phiIndex=EmcID::phi_module(showerId);
196 if (ene5x5>0.4&&ene5x5<4&&npart==1&&(m_BarrelOrEndcap==1)){
197 GoodShowers.push_back(theShower);
198 iGood.push_back((*itTrk)->trackId());
199 }
200 else if (ene5x5>0.4&&ene5x5<4&&(npart==2||npart==0)&&(m_BarrelOrEndcap==2)){
201 GoodShowers.push_back(theShower);
202 iGood.push_back((*itTrk)->trackId());
203 }
204 else if (ene5x5>0.4&&ene5x5<4&&(m_BarrelOrEndcap==0)){
205 GoodShowers.push_back(theShower);
206 iGood.push_back((*itTrk)->trackId());
207 }
208 }
209 // good photon cut will be set here
210 }
211 // Finish Good Photon Selection
212 //
213
214 double MaxE(0), MaxPhi, MaxThe;
215 int MaxId;
216 double SecE(0), SecPhi, SecThe;
217 for(int i=0; i<GoodShowers.size(); i++) {
218 RecEmcShower *theShower = GoodShowers[i];
219 double eraw = theShower->energy();
220 if(eraw> MaxE) {
221 MaxId =i;
222 MaxE =eraw;
223 MaxPhi = theShower->phi();
224 MaxThe = theShower->theta();
225 }
226 }
227 for(int i=0; i<GoodShowers.size(); i++) {
228 RecEmcShower *theShower = GoodShowers[i];
229 double eraw = theShower->energy();
230 if(i!=MaxId&&eraw>SecE) {
231 SecE =eraw;
232 SecPhi = theShower->phi();
233 SecThe = theShower->theta();
234 }
235 }
236
237 double dphi = (fabs(MaxPhi-SecPhi)-PI)*180./PI;
238 double dthe = (fabs(MaxThe+SecThe)-PI)*180./PI;
239
240 int total1 = 0; //mdc hit
241 int total2 = 0; //mdc hit
242 if(GoodShowers.size()>=2 && MaxE>1.0 && abs(dthe)<3
243 && ((dphi>-25&&dphi<-4)||(dphi>2&&dphi<20)) && etot>2.7) {
244
245 double phi1 = MaxPhi<0 ? MaxPhi+CLHEP::twopi : MaxPhi;
246 double phi2 = SecPhi<0 ? SecPhi+CLHEP::twopi : SecPhi;
247
248 //Define sector (phi11,phi12) and (phi21,phi22)
249 double phi11=min(phi1,phi2);
250 double phi22=max(phi1,phi2);
251 double phi12=(phi11+phi22-CLHEP::pi)*0.5;
252 double phi21=(phi11+phi22+CLHEP::pi)*0.5;
253 if(phi12<0.) phi12 += CLHEP::twopi;
254 if(phi21>CLHEP::twopi) phi21 -= CLHEP::twopi;
255
256 SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(), "/Event/Digi/MdcDigiCol");
257 if (!mdcDigiCol) {
258 log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
259 return StatusCode::FAILURE;
260 }
261 int hitnum = mdcDigiCol->size();
262 for (int i = 0;i< hitnum;i++ ) {
263 MdcDigi* digi= dynamic_cast<MdcDigi*>(mdcDigiCol->containedObject(i));
264 double time = digi->getTimeChannel();
265 double charge = digi->getChargeChannel();
266 if (time == 0x7FFFFFFF || charge == 0x7FFFFFFF) continue;
267 Identifier id(digi->identify());
268 unsigned int iphi=MdcID::wire(id);
269 unsigned int ilayer=MdcID::layer(id);
270 if(ilayer>=43)
271 log << MSG::ERROR << "MDC(" << ilayer <<","<<iphi <<")"<<endreq;
272 double phi=CLHEP::twopi*iphi/idmax[ilayer];
273 if(WhetherSector(phi,phi11,phi12)) total1++;
274 else if(WhetherSector(phi,phi21,phi22)) total2++;
275 //cout<<"phi="<<phi<<",phi11="<<phi11<<",phi12="<<phi12
276 //<<",phi21="<<phi21<<",phi22="<<phi22<<endl;
277 }
278 //cout<<"total1="<<total1<<"total2="<<total2<<endl;
279 if ( (m_BarrelOrEndcap==1&&total1>15 &&total2>15)
280 || (m_BarrelOrEndcap!=1&&total1>5 &&total2>5) ) {
281 setFilterPassed(true);
282 selected++;
283 }
284 }
285
286 if(m_output) {
287 m_etot=etot;
288 m_mdc_hit1=total1;
289 m_mdc_hit2=total2;
290 m_sh1_ene = MaxE;
291 m_sh1_theta = MaxThe;
292 m_sh1_phi = MaxPhi;
293 m_sh2_ene = SecE;
294 m_sh2_theta = SecThe;
295 m_sh2_phi = SecPhi;
296 m_di_phi = dphi;
297 m_di_the = dthe;
298 m_tuple1->write();
299 }
300
301 return StatusCode::SUCCESS;
302
303}
selected
int selected
Definition BhabhaPreSelect.cxx:32
phi2
Double_t phi2
Definition DataBase/tau_mode.c:7
etot
Double_t etot
Definition DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition DataBase/tau_mode.c:7
time
Double_t time
Definition DataBase/tau_mode.c:7
PI
#define PI
Definition EvtD0ToKKpi0.cc:41
EvtRecTrackIterator
EvtRecTrackCol::iterator EvtRecTrackIterator
Definition EvtRecTrack.h:111
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:53
Vint
std::vector< int > Vint
Definition Gam4pikp.cxx:52
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
BhabhaPreSelect::WhetherSector
bool WhetherSector(double, double=0., double=CLHEP::twopi)
Definition BhabhaPreSelect.cxx:320
DstEmcShower::theta
double theta() const
Definition DstEmcShower.h:38
DstEmcShower::phi
double phi() const
Definition DstEmcShower.h:39
DstEmcShower::e5x5
double e5x5() const
Definition DstEmcShower.h:49
DstEmcShower::energy
double energy() const
Definition DstEmcShower.h:45
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition EmcID.cxx:38
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:43
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:48
Identifier::clear
void clear()
Reset to invalid state.
Definition Identifier.h:152
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:54
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition RawData.cxx:45
RawData::getTimeChannel
unsigned int getTimeChannel() const
Definition RawData.cxx:40
RecEmcShower::getShowerId
RecEmcID getShowerId() const
Definition RecEmcShower.h:55
EventModel::EvtRec::EvtRecEvent
_EXTERN_ std::string EvtRecEvent
Definition EventModel.h:116
EventModel::EvtRec::EvtRecTrackCol
_EXTERN_ std::string EvtRecTrackCol
Definition EventModel.h:117
charge
float charge
Definition prepare_data.cxx:28

◆ finalize()

StatusCode BhabhaPreSelect::finalize ( )

Definition at line 306 of file BhabhaPreSelect.cxx.

306 {
307
308 MsgStream log(msgSvc(), name());
309 log << MSG::INFO << "in finalize()" << endmsg;
310 cout<<"m_events "<<m_events<<endl;
311 cout<<" selected "<<selected<<endl;
312 cout<<"m_rejected "<<m_rejected<<endl;
313 cout<<"m_oneProngsSelected "<<m_oneProngsSelected<<endl;
314 cout<<"m_twoProngsMatchedSelected "<<m_twoProngsMatchedSelected<<endl;
315 cout<<"m_twoProngsOneMatchedSelected "<<m_twoProngsOneMatchedSelected<<endl;
316
317 return StatusCode::SUCCESS;
318}

◆ initialize()

StatusCode BhabhaPreSelect::initialize ( )

Definition at line 60 of file BhabhaPreSelect.cxx.

60 {
61 MsgStream log(msgSvc(), name());
62
63 log << MSG::INFO << "in initialize()" << endmsg;
64 cout<<" BarrelOrEndcap "<<m_BarrelOrEndcap<<endl;
65
66 if(m_output) {
67 StatusCode status;
68 NTuplePtr nt1(ntupleSvc(), "FILE1/bhabha");
69 if ( nt1 ) m_tuple1 = nt1;
70 else {
71 m_tuple1 = ntupleSvc()->book ("FILE1/bhabha", CLID_ColumnWiseTuple, "N-Tuple example");
72 if ( m_tuple1 ) {
73 status = m_tuple1->addItem ("sh1_ene",m_sh1_ene);
74 status = m_tuple1->addItem ("sh1_theta", m_sh1_theta);
75 status = m_tuple1->addItem ("sh1_phi", m_sh1_phi);
76 status = m_tuple1->addItem ("sh2_ene",m_sh2_ene);
77 status = m_tuple1->addItem ("sh2_theta", m_sh2_theta);
78 status = m_tuple1->addItem ("sh2_phi", m_sh2_phi);
79 status = m_tuple1->addItem ("di_phi", m_di_phi);
80 status = m_tuple1->addItem ("di_the", m_di_the);
81 status = m_tuple1->addItem ("acolli", m_acolli);
82 status = m_tuple1->addItem ("etot", m_etot);
83 status = m_tuple1->addItem ("mdc_hit1", m_mdc_hit1);
84 status = m_tuple1->addItem ("mdc_hit2", m_mdc_hit2);
85
86 }
87 else {
88 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
89 return StatusCode::FAILURE;
90 }
91 }
92
93
94 NTuplePtr nt2(ntupleSvc(), "FILE1/bha1");
95 if ( nt2 ) m_tuple2 = nt2;
96 else {
97 m_tuple2 = ntupleSvc()->book ("FILE1/bha1", CLID_ColumnWiseTuple, "N-Tuple example");
98 if ( m_tuple2 ) {
99 status = m_tuple2->addItem ("sh_ene",m_sh_ene);
100 status = m_tuple2->addItem ("sh_theta", m_sh_theta);
101 status = m_tuple2->addItem ("sh_phi", m_sh_phi);
102 }
103 else {
104 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple2) << endmsg;
105 return StatusCode::FAILURE;
106 }
107 }
108 }
109
110
111 //
112 //--------end of book--------
113 //
114
115 m_rejected=0;
116 m_events=0;
117 m_oneProngsSelected=0;
118 m_twoProngsMatchedSelected=0;
119 m_twoProngsOneMatchedSelected=0;
120 m_selectedTrkID1=999;
121 m_selectedTrkID2=999;
122
123 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
124 return StatusCode::SUCCESS;
125
126}
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55

◆ selectedTrkID1()

int BhabhaPreSelect::selectedTrkID1 ( ) const
inline

Definition at line 25 of file BhabhaPreSelect.h.

26 {
27 return m_selectedTrkID1;
28 }

◆ selectedTrkID2()

int BhabhaPreSelect::selectedTrkID2 ( ) const
inline

Definition at line 29 of file BhabhaPreSelect.h.

30 {
31 return m_selectedTrkID2;
32 }

◆ selectedType()

int BhabhaPreSelect::selectedType ( ) const
inline

Definition at line 21 of file BhabhaPreSelect.h.

22 {
23 return m_selectedType;
24 }

◆ WhetherSector()

bool BhabhaPreSelect::WhetherSector ( double ph,
double ph1 = 0.,
double ph2 = CLHEP::twopi )

Definition at line 320 of file BhabhaPreSelect.cxx.

320 {
321 double phi1=min(ph1,ph2);
322 double phi2=max(ph1,ph2);
323 double delta=0.0610865; //3.5*3.1415926/180.
324 if((phi2-phi1)<CLHEP::pi){
325 phi1 -=delta;
326 phi2 +=delta;
327 if(phi1<0.) phi1 += CLHEP::twopi;
328 if(phi2>CLHEP::twopi) phi2 -= CLHEP::twopi;
329 double tmp1=min(phi1,phi2);
330 double tmp2=max(phi1,phi2);
331 phi1=tmp1;
332 phi2=tmp2;
333 }
334 else{
335 phi1 +=delta;
336 phi2 -=delta;
337 }
338 if((phi2-phi1)<CLHEP::pi){
339 if(ph<=phi2&&ph>=phi1) return true;
340 else return false;
341 }
342 else{
343 if(ph>=phi2||ph<=phi1) return true;
344 else return false;
345 }
346}
delta
const double delta
Definition CheckRung/check/gainVsDoca.cxx:32

Referenced by execute().

The documentation for this class was generated from the following files: