BOSS 7.1.2
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 25 of file Bhlumi.h.

Constructor & Destructor Documentation

◆ Bhlumi()

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

Definition at line 106 of file Bhlumi.cxx.

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

Member Function Documentation

◆ execute()

StatusCode Bhlumi::execute ( )

Definition at line 221 of file Bhlumi.cxx.

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

◆ finalize()

StatusCode Bhlumi::finalize ( )

Definition at line 339 of file Bhlumi.cxx.

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

◆ 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 119 of file Bhlumi.cxx.

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

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