BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Bhlumi Class Reference

#include <Bhlumi.h>

Inheritance diagram for Bhlumi:

Public Member Functions

 Bhlumi (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode initialize ()
StatusCode execute ()
StatusCode finalize ()

Detailed Description

Definition at line 22 of file Bhlumi.h.

Constructor & Destructor Documentation

◆ Bhlumi()

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

Definition at line 99 of file Bhlumi.cxx.

100 : Algorithm( name, pSvcLocator ) {
101 declareProperty( "CMEnergy", m_cmEnergy = 3.097 ); // 2*Ebeam [GeV]
102 declareProperty( "AngleMode",
103 m_angleMode = 0 ); // 0: rad; 1: degree; 2: cos(theta);
104 // if 1 or 2, angle will be controlled absolutely by user
105 declareProperty( "MinThetaAngle", m_minThetaAngle = 0.24 ); // [rad]
106 declareProperty( "MaxThetaAngle", m_maxThetaAngle = 0.58 ); // [rad]
107 declareProperty( "SoftPhotonCut",
108 m_infraredCut = 1E-4 ); // dimensionless, Ephoton > m_infraredCut*sqrt(s)/2
109 m_initSeed.clear();
110 // m_initSeed.push_back(54217137); m_initSeed.push_back(0); m_initSeed.push_back(0);
111 // declareProperty("InitializedSeed", m_initSeed);
112}

Referenced by Bhlumi().

Member Function Documentation

◆ execute()

StatusCode Bhlumi::execute ( )

Definition at line 224 of file Bhlumi.cxx.

224 {
225 MsgStream log( msgSvc(), name() );
226 log << MSG::INFO << "Bhlumi executing" << endmsg;
227
228 BHLUMI( 0, xpar, npar );
229
230 if ( log.level() < MSG::INFO )
231 {
232 DUMPS( 6 );
233 // dump output to file
234 // DUMPS(16);
235 }
236
237 int npart = 0;
238
239 // Fill event information
240 GenEvent* evt = new GenEvent( 1, 1 );
241
242 GenVertex* prod_vtx = new GenVertex();
243 // prod_vtx->add_particle_out( p );
244 evt->add_vertex( prod_vtx );
245
246 // incoming beam e+
247 GenParticle* p = new GenParticle(
248 CLHEP::HepLorentzVector( MOMSET.p1[0], MOMSET.p1[1], MOMSET.p1[2], MOMSET.p1[3] ), -11,
249 3 );
250 p->suggest_barcode( ++npart );
251 prod_vtx->add_particle_in( p );
252
253 // incoming beam e-
254 p = new GenParticle(
255 CLHEP::HepLorentzVector( MOMSET.q1[0], MOMSET.q1[1], MOMSET.q1[2], MOMSET.q1[3] ), 11,
256 3 );
257 p->suggest_barcode( ++npart );
258 prod_vtx->add_particle_in( p );
259
260 // scattered e+
261 p = new GenParticle(
262 CLHEP::HepLorentzVector( MOMSET.p2[0], MOMSET.p2[1], MOMSET.p2[2], MOMSET.p2[3] ), -11,
263 1 );
264 p->suggest_barcode( ++npart );
265 prod_vtx->add_particle_out( p );
266
267 // scattered e-
268 p = new GenParticle(
269 CLHEP::HepLorentzVector( MOMSET.q2[0], MOMSET.q2[1], MOMSET.q2[2], MOMSET.q2[3] ), 11,
270 1 );
271 p->suggest_barcode( ++npart );
272 prod_vtx->add_particle_out( p );
273
274 int iphot = 0;
275 for ( iphot = 0; iphot < MOMSET.nphot; iphot++ )
276 {
277 // gamma
278 p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.phot[0][iphot], MOMSET.phot[1][iphot],
279 MOMSET.phot[2][iphot],
280 MOMSET.phot[3][iphot] ),
281 22, 1 );
282 p->suggest_barcode( ++npart );
283 prod_vtx->add_particle_out( p );
284 }
285
286 if ( log.level() < MSG::INFO ) { evt->print(); }
287
288 // Check if the McCollection already exists
289 SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
290 if ( anMcCol != 0 )
291 {
292 // Add event to existing collection
293 MsgStream log( msgSvc(), name() );
294 log << MSG::INFO << "Add McGenEvent to existing collection" << endmsg;
295 McGenEvent* mcEvent = new McGenEvent( evt );
296 anMcCol->push_back( mcEvent );
297 }
298 else
299 {
300 // Create Collection and add to the transient store
301 McGenEventCol* mcColl = new McGenEventCol;
302 McGenEvent* mcEvent = new McGenEvent( evt );
303 mcColl->push_back( mcEvent );
304 StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
305 if ( sc != StatusCode::SUCCESS )
306 {
307 log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
308 delete mcColl;
309 delete evt;
310 delete mcEvent;
311 return StatusCode::FAILURE;
312 }
313 else
314 {
315 log << MSG::INFO << "McGenEventCol created and " << npart
316 << " particles stored in McGenEvent" << endmsg;
317 }
318 }
319
320 // retrieve event from Transient Store (Storegate)
321 /* SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen");
322 if ( McEvtColl == 0 )
323 {
324 log << MSG::ERROR << "Could not retrieve McGenEventCollection" << endmsg;
325 return StatusCode::FAILURE;
326 };
327
328 McGenEventCol::iterator mcItr;
329 for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ )
330 {
331 HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt();
332 // MeVToGeV( hEvt );
333 // callEvtGen( hEvt );
334 // GeVToMeV( hEvt );
335 };
336 */
337
338 return StatusCode::SUCCESS;
339}
MOMSET
#define MOMSET
Definition Babayaga.cxx:43
BHLUMI
#define BHLUMI(MODE, XPAR, NPAR)
Definition Bhlumi.cxx:58
DUMPS
#define DUMPS(NOUT)
Definition Bhlumi.cxx:55
McGenEventCol
ObjectVector< McGenEvent > McGenEventCol
Definition Event/GeneratorObject/include/GeneratorObject/McGenEvent.h:39
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12

◆ finalize()

StatusCode Bhlumi::finalize ( )

Definition at line 341 of file Bhlumi.cxx.

341 {
342 MsgStream log( msgSvc(), name() );
343
344 BHLUMI( 2, xpar, npar );
345
346 log << MSG::INFO << "Bhlumi finalized" << endmsg;
347
348 return StatusCode::SUCCESS;
349}

◆ initialize()

StatusCode Bhlumi::initialize ( )

||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! User should cross-check the folowing two output cross sections ! which are calculated and printed at the very end of the output: ! Workshop95, Table14, BARE1 WW for zmin=0.5: KeyGen=3, KeyPia=0, KeyZet=0 ! Workshop95, Table18, CALO2 WW for zmin=0.5: KeyGen=3, KeyPia=2, KeyZet=1 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

Input parameters for Bhlumi

Technical parameters, nothing should depend on them:

Try both options for KeyWgt, result should be the same

Physics parameters:

2*Ebeam [GeV], as in Workshop95

Definition at line 114 of file Bhlumi.cxx.

114 {
115
116 MsgStream log( msgSvc(), name() );
117
118 log << MSG::INFO << "Bhlumi initialize" << endmsg;
119
120 static const bool CREATEIFNOTTHERE( true );
121 StatusCode RndmStatus = service( "BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE );
122 if ( !RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc )
123 {
124 log << MSG::ERROR << " Could not initialize Random Number Service" << endmsg;
125 return RndmStatus;
126 }
127 CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine( "Bhlumi" );
128 engine->showStatus();
129 BhlumiRandom::setRandomEngine( engine );
130
131 GLIMIT( 50000 );
132
133 /*!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
134 ! User should cross-check the folowing two output cross sections
135 ! which are calculated and printed at the very end of the output:
136 ! Workshop95, Table14, BARE1 WW for zmin=0.5: KeyGen=3, KeyPia=0,
137 KeyZet=0
138 ! Workshop95, Table18, CALO2 WW for zmin=0.5: KeyGen=3, KeyPia=2,
139 KeyZet=1
140 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
141 */
142 //!---------------------------------------------------------------------------
143 //! Input parameters for Bhlumi
144 //!----------------------------------------------------------------------
145 //! Technical parameters, nothing should depend on them:
146 int KeyGen = 3; // ! Multiphoton Bhlumi
147 int KeyRem = 1; // ! No remooval of photons below epsCM
148 KeyRem = 0; // ! Remooval of photons below epsCM, Necessary for Z!!!
149 //! Try both options for KeyWgt, result should be the same
150 int KeyWgt = 2; // ! weighted events, with t generation down to zero
151 KeyWgt = 0; // ! WT=1 events, for detector simulation
152 int KeyRnd = 1; // ! RANMAR random numbers
153 int KeyOpt = 1000 * KeyGen + 100 * KeyRem + 10 * KeyWgt + KeyRnd;
154 //!--------------------------------------------------
155 //! Physics parameters:
156 int KeyPia = 0; // ! Vacuum polarization OFF
157 int KeyMod = 2; // ! Matrix element default version 4.x
158 KeyPia = 2; // ! Vacuum polarization ON
159 int KeyZet = 0; // ! Z contribution OFF
160 KeyZet = 1; // ! Z contribution ON
161 int KeyRad = 1000 * KeyZet + 10 * KeyMod + KeyPia;
162 //!--------------------------------------
163 npar[0] = KeyOpt;
164 npar[1] = KeyRad;
165 double CmsEne = m_cmEnergy; //! 2*Ebeam [GeV], as in Workshop95
166 xpar[0] = CmsEne;
167 double th1, th2, thmin, thmax;
168 if ( m_angleMode == 0 )
169 {
170 th1 = m_minThetaAngle; // ! Detector range ThetaMin [rad]
171 th2 = m_maxThetaAngle; // ! Detector range ThetaMax [rad]
172 thmin = 0.7 * th1; // ! thmin has to be lower than th1
173 thmax = 2.0 * th2; // ! thmax has to be higher than th2
174 if ( thmin < 0. || thmax > 3.1415926 )
175 {
176 log << MSG::WARNING
177 << "input angles exceed range (0~pi), so effect will be unprospectable" << endmsg;
178 return StatusCode::FAILURE;
179 }
180 }
181 else if ( m_angleMode == 1 )
182 {
183 th1 = m_minThetaAngle * 3.1415926 / 180.;
184 th2 = m_maxThetaAngle * 3.1415926 / 180.;
185 // not multiply 0.7(2.0) coefficient while large angle
186 thmin = th1;
187 thmax = th2;
188 }
189 else if ( m_angleMode == 2 )
190 {
191 th1 = acos( max( m_minThetaAngle, m_maxThetaAngle ) );
192 th2 = acos( min( m_minThetaAngle, m_maxThetaAngle ) );
193 // not multiply 0.7(2.0) coefficient while large angle
194 thmin = th1;
195 thmax = th2;
196 }
197 else
198 {
199 log << MSG::FATAL << "unknown angle mode!" << endmsg;
200 return StatusCode::FAILURE;
201 }
202 if ( thmin < 0. || thmax > 3.1415926 )
203 {
204 log << MSG::FATAL << "input angles exceed range (0~pi), unprospectable" << endmsg;
205 return StatusCode::FAILURE;
206 }
207 else if ( thmin > thmax )
208 {
209 log << MSG::FATAL << "thmin>thmax, unprospectable" << endmsg;
210 return StatusCode::FAILURE;
211 }
212 if ( KeyWgt == 2 ) thmin = th1; // ! Because generation below th1 is on!!!
213 xpar[1] = CmsEne * CmsEne * ( 1 - cos( thmin ) ) / 2; // ! TransMin [GeV**2]
214 xpar[2] = CmsEne * CmsEne * ( 1 - cos( thmax ) ) / 2; // ! TransMax [GeV**2]
215 xpar[3] = m_infraredCut; // ! Infrared cut on photon energy
216
217 // MARINI(m_initSeed[0], m_initSeed[1], m_initSeed[2]);
218
219 BHLUMI( -1, xpar, npar );
220
221 return StatusCode::SUCCESS;
222}
GLIMIT
#define GLIMIT(LENMX)
Definition Bhlumi.cxx:52
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
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
BhlumiRandom::setRandomEngine
static void setRandomEngine(CLHEP::HepRandomEngine *randomEngine)
Definition BhlumiRandom.cxx:9
The documentation for this class was generated from the following files: