BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
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 10 of file BhabhaPreSelect.h.

Constructor & Destructor Documentation

◆ BhabhaPreSelect()

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

Definition at line 47 of file BhabhaPreSelect.cxx.

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

Referenced by BhabhaPreSelect().

Member Function Documentation

◆ execute()

StatusCode BhabhaPreSelect::execute ( )

Definition at line 130 of file BhabhaPreSelect.cxx.

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

◆ finalize()

StatusCode BhabhaPreSelect::finalize ( )

Definition at line 321 of file BhabhaPreSelect.cxx.

321 {
322
323 MsgStream log( msgSvc(), name() );
324 log << MSG::INFO << "in finalize()" << endmsg;
325 cout << "m_events " << m_events << endl;
326 cout << " selected " << selected << endl;
327 cout << "m_rejected " << m_rejected << endl;
328 cout << "m_oneProngsSelected " << m_oneProngsSelected << endl;
329 cout << "m_twoProngsMatchedSelected " << m_twoProngsMatchedSelected << endl;
330 cout << "m_twoProngsOneMatchedSelected " << m_twoProngsOneMatchedSelected << endl;
331
332 return StatusCode::SUCCESS;
333}

◆ initialize()

StatusCode BhabhaPreSelect::initialize ( )

Definition at line 58 of file BhabhaPreSelect.cxx.

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

◆ selectedTrkID1()

int BhabhaPreSelect::selectedTrkID1 ( ) const
inline

Definition at line 21 of file BhabhaPreSelect.h.

21{ return m_selectedTrkID1; }

◆ selectedTrkID2()

int BhabhaPreSelect::selectedTrkID2 ( ) const
inline

Definition at line 22 of file BhabhaPreSelect.h.

22{ return m_selectedTrkID2; }

◆ selectedType()

int BhabhaPreSelect::selectedType ( ) const
inline

Definition at line 20 of file BhabhaPreSelect.h.

20{ return m_selectedType; }

◆ WhetherSector()

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

Definition at line 335 of file BhabhaPreSelect.cxx.

335 {
336 double phi1 = min( ph1, ph2 );
337 double phi2 = max( ph1, ph2 );
338 double delta = 0.0610865; // 3.5*3.1415926/180.
339 if ( ( phi2 - phi1 ) < CLHEP::pi )
340 {
341 phi1 -= delta;
342 phi2 += delta;
343 if ( phi1 < 0. ) phi1 += CLHEP::twopi;
344 if ( phi2 > CLHEP::twopi ) phi2 -= CLHEP::twopi;
345 double tmp1 = min( phi1, phi2 );
346 double tmp2 = max( phi1, phi2 );
347 phi1 = tmp1;
348 phi2 = tmp2;
349 }
350 else
351 {
352 phi1 += delta;
353 phi2 -= delta;
354 }
355 if ( ( phi2 - phi1 ) < CLHEP::pi )
356 {
357 if ( ph <= phi2 && ph >= phi1 ) return true;
358 else return false;
359 }
360 else
361 {
362 if ( ph >= phi2 || ph <= phi1 ) return true;
363 else return false;
364 }
365}

Referenced by execute().

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