BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Ekhara.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2#include "GaudiKernel/SmartDataPtr.h"
3#include <cstdio>
4#include <unistd.h>
5
6#include "BesRndmGenSvc/IBesRndmGenSvc.h"
7#include "GeneratorObject/McGenEvent.h"
8
9#include "HepMC/GenEvent.h"
10
11#include "Ekhara.h"
12#include "EkharaDef.h"
13#include "EkharaRandom.h"
14
15#include "CLHEP/Vector/LorentzVector.h"
16
17#include "TMath.h"
18
19#ifdef DECLARE_COMPONENT
20DECLARE_COMPONENT( Ekhara )
21#endif
22
23using namespace HepMC;
24using namespace CLHEP;
25
26Ekhara::Ekhara( const std::string& name, ISvcLocator* pSvcLocator )
27 : Algorithm( name, pSvcLocator ) {
28 declareProperty( "Ecm", m_cmsEnergy = 3.773 ); //!- CMS-energy (GeV)
29 declareProperty( "InitialEvents",
30 m_numUpperBound = 50000 ); //!- # of events to determine the max xs
31 declareProperty( "IgnoreUpperBound",
32 m_ignoreUpperBound = 0.0 ); //!- handle UB underestim <0> = halt,
33 //!- <positive> = use new UB if overshoot,
34 //!- <negative> = ignore
35 declareProperty( "UpperBoundEnlarge",
36 m_upperBoundEnlarge = 1.2 ); //!- Enlarge factor for upper bound
37 declareProperty( "FinalStateID", m_finalState = 2 ); //!- 0:pi0pi0, 1:pi+pi-, 2:pi0, 3:eta,
38 //! 4:eta', 5:Chi_c0, 6:Chi_c1, 7:Chi_c2
39
40 declareProperty( "MesonProdAmplitudes",
41 m_mesonAmplitudes = 2 ); //!- pi0,eta,eta': s(1); t(2); s+t(3); t:signal
42 declareProperty( "IncludeVacuumPolarization",
43 m_switchVP = false ); //!- include vacuum polarization (true or false)
44 declareProperty( "ConsiderNLO",
45 m_switchNLO = false ); //!- NLO radiative corrections (true or false)
46 declareProperty( "OutputWeightedNLOEvents",
47 m_nloWithWeights =
48 false ); //!- return weighted events as output for NLO (true or false)
49 declareProperty( "PhaseSpaceRegulator",
50 m_eps_ph = 0.01 ); //!- phase space regulator eps_ph: mgamma=me*eps_ph
51 declareProperty( "MesonFormfactor", m_TFF_ID = 8 ); //!- pi0: WZWconst(1); rho pole(2);
52 //!- LMD(3); LMD+V(4); LMD+Vnew(5);
53 //!- 1-Octet(7); 2-Octet(8);
54 //!- eta/eta': WZWconst(1); VMD(3);
55 //!- 1-Octet(7); 2-Octet(8);
56 declareProperty( "PositronThetaMin", m_posiThetaMin = 0.0 ); //!- e+ minimum angle
57 declareProperty( "PositronThetaMax", m_posiThetaMax = 180.0 ); //!- e+ maximum angle
58 declareProperty( "PositronEnergyMin", m_posiEnergyMin = 0.0 ); //!- e+ minimum Energy
59 declareProperty( "PositronEnergyMax", m_posiEnergyMax = 110.0 ); //!- e+ maximum Energy
60 declareProperty( "ElectronThetaMin", m_elecThetaMin = 0.0 ); //!- e- minimum angle
61 declareProperty( "ElectronThetaMax", m_elecThetaMax = 180.0 ); //!- e- maximum angle
62 declareProperty( "ElectronEnergyMin", m_elecEnergyMin = 0.0 ); //!- e- minimum Energy
63 declareProperty( "ElectronEnergyMax", m_elecEnergyMax = 110.0 ); //!- e- maximum Energy
64
65 declareProperty( "MesonThetaMin", m_mesonThetaMin = 0.0 ); //!- pi0, eta, eta' minimum angle
66 declareProperty( "MesonThetaMax",
67 m_mesonThetaMax = 180.0 ); //!- pi0, eta, eta' maximum angle
68 declareProperty( "MesonEnergyMin",
69 m_mesonEnergyMin = 0.0 ); //!- pi0, eta, eta' minimum Energy
70 declareProperty( "MesonEnergyMax",
71 m_mesonEnergyMax = 100.0 ); //!- pi0, eta, eta' maximum Energy
72
73 declareProperty( "TwoPionPhaseSpaceAlg",
74 m_twoPionPhspAlg =
75 1 ); //!- the choice of generation algorithm: 0-v.1.0; 1-v.3.1
76 declareProperty( "TwoPionProdAmplitudes",
77 m_twoPionAmplitudes =
78 4 ); //!- pi+pi-: s(1); s+t(2); s+t+2g_Born(3); s+t+2g_Full(4)
79 declareProperty( "TwoPionFormFactor",
80 m_twoPionFormFactor = 2 ); //!- the choice of the pion form factor to be
81 //! used: 0-KS; 1-GS; 2-GS new
82 declareProperty( "TwoPionThetaMin", m_twoPionThetaMin = 0.0 ); //!- pion angle cut
83 declareProperty( "TwoPionThetaMax", m_twoPionThetaMax = 180.0 );
84 declareProperty( "TwoPionMissThetaMin",
85 m_twoPionMissThetaMin = 0.0 ); //!- missing momentum angle cut
86 declareProperty( "TwoPionMissThetaMax", m_twoPionMissThetaMax = 180.0 );
87
88 declareProperty( "Chi_cjThetaMin", m_chicjThetaMin = 0.0 ); //!- Chi_cj minimum angle
89 declareProperty( "Chi_cjThetaMax", m_chicjThetaMax = 180.0 ); //!- Chi_cj maximum angle
90 declareProperty( "Chi_cjEnergyMin", m_chicjEnergyMin = 0.0 ); //!- Chi_cj minimum energy
91 declareProperty( "Chi_cjEnergyMax", m_chicjEnergyMax = 100.0 ); //!- Chi_cj minimum energy
92
93 declareProperty( "TaggingAngle",
94 m_taggingAngle = 20. ); //!- Tagging angle (symmetric) in degrees
95 declareProperty( "TaggingQsquare",
96 m_taggingQsquare = 0.3 ); //!- minimum Q^2 of tagged lepton in GeV^2
97 declareProperty( "TaggingMode",
98 m_taggingMode =
99 0 ); //!- Tagging mode 0:no tagging; 1/11:untagged
100 //!- 2/12: single tagging; 3/13:double tagging
101}
102
103StatusCode Ekhara::initialize() {
104 MsgStream log( msgSvc(), name() );
105 log << MSG::DEBUG << "Ekhara in initialize()" << endmsg;
106
107 //!- Set Bes unified random engine
108 static const bool CREATEIFNOTTHERE( true );
109 StatusCode RndmStatus = service( "BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE );
110 if ( !RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc )
111 {
112 log << MSG::ERROR << " Could not initialize Random Number Service" << endmsg;
113 return RndmStatus;
114 }
115 HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine( "Ekhara" );
116 EkharaRandom::setRandomEngine( engine );
117
118 //!- Perform consitency checks of jobOptions and print screen messages
119 if ( !checkAndPrintParameters() ) return StatusCode::FAILURE;
120
121 //!- Transfer integer parameters to FORTRAN common blocks
122 CHANNELSEL.channel_id = m_finalState;
123 CHANNELSEL.NLO1P = 0;
124 if ( m_switchNLO ) CHANNELSEL.NLO1P = 1;
125 CHANNELSEL.VPSW = 0;
126 if ( m_switchVP ) CHANNELSEL.VPSW = 1;
127 CHANNELSEL.phsp_2pi = m_twoPionPhspAlg;
128 SWDIAG.sw_2pi = m_twoPionAmplitudes;
129 SWDIAG.sw_1pi = m_mesonAmplitudes;
130 PIONFFSW.piggFFsw = m_TFF_ID;
131 TAGGINGMODE.tagmode = m_taggingMode;
132 FFPARAMSET.FF_pion = m_twoPionFormFactor;
133 NLOTYPE.weighted = 0;
134 if ( m_nloWithWeights ) NLOTYPE.weighted = 1;
135
136 //!- Prepare floating point parameters to be transfered to quadruple precision
137 double xpar[27];
138 xpar[0] = m_cmsEnergy;
139 xpar[1] = m_numUpperBound;
140 xpar[2] = m_ignoreUpperBound;
141 xpar[3] = m_upperBoundEnlarge;
142 xpar[4] = m_eps_ph;
143 xpar[5] = m_posiThetaMin;
144 xpar[6] = m_posiThetaMax;
145 xpar[7] = m_posiEnergyMin;
146 xpar[8] = m_posiEnergyMax;
147 xpar[9] = m_elecThetaMin;
148 xpar[10] = m_elecThetaMax;
149 xpar[11] = m_elecEnergyMin;
150 xpar[12] = m_elecEnergyMax;
151 xpar[13] = m_mesonThetaMin;
152 xpar[14] = m_mesonThetaMax;
153 xpar[15] = m_mesonEnergyMin;
154 xpar[16] = m_mesonEnergyMax;
155 xpar[17] = m_twoPionThetaMin;
156 xpar[18] = m_twoPionThetaMax;
157 xpar[19] = m_twoPionMissThetaMin;
158 xpar[20] = m_twoPionMissThetaMax;
159 xpar[21] = m_chicjThetaMin;
160 xpar[22] = m_chicjThetaMax;
161 xpar[23] = m_chicjEnergyMin;
162 xpar[24] = m_chicjEnergyMax;
163 xpar[25] = m_taggingAngle;
164 xpar[26] = m_taggingQsquare;
165
166 //!- Transfer floating point parameters to FORTAN routines
167 //!- and initialize FORTAN routines
168 BOSS_INIT_EKHARA( xpar );
169
170 if ( log.level() < MSG::INFO ) { DIAGNOSE(); }
171
172 //!- Initialize event counter
173 // should be replaced by request to ApplicationManager for EvtMax
174 m_totalEvents = 0;
175
176 return StatusCode::SUCCESS;
177}
178
179StatusCode Ekhara::execute() {
180 MsgStream log( msgSvc(), name() );
181 log << MSG::DEBUG << "Ekhara in execute()" << endmsg;
182
183 //!- Generate an event
184 if ( m_finalState >= 2 && m_finalState <= 4 && m_switchNLO && m_nloWithWeights )
185 { EKHARA( 2 ); }
186 else { EKHARA( 0 ); }
187
188 double p1[4], p2[4], q1[4], q2[4];
189 double pi1[4], pi2[4], qpion[4], qcj[4], kphp[4];
190
191 for ( int i = 0; i < 4; i++ )
192 {
193 p1[i] = 0.0; // four vector of incoming positron
194 p2[i] = 0.0; // four vector of incoming electron
195 q1[i] = 0.0; // four vector of outgoung positron
196 q2[i] = 0.0; // four vector of outgoing electron
197 pi1[i] = 0.0; // four vector of produced pi0,pi+ (m_finalState == 0,1)
198 pi2[i] = 0.0; // four vector of produced pi0,pi- (m_finalState == 0,1)
199 qpion[i] = 0.0; // four vector of produced pi0,eta,eta' (m_finalState == 2,3,4)
200 qcj[i] = 0.0; // four vector of produced chi_cj (m_finalState == 5,6,7)
201 kphp[i] = 0.0; // four vector of radiative photon
202 }
203
204 // Fill event information
205 GenEvent* evt = new GenEvent( 1, 1 );
206
207 GenVertex* prod_vtx = new GenVertex();
208 evt->add_vertex( prod_vtx );
209
210 GET_FOURMOMENTA_LEPTONS( p1, p2, q1, q2 );
211
212 // incoming beam e+
213 GenParticle* posi =
214 new GenParticle( CLHEP::HepLorentzVector( p1[1], p1[2], p1[3], p1[0] ), -11, 3 );
215 posi->suggest_barcode( 1 );
216 prod_vtx->add_particle_in( posi );
217
218 // incoming beam e-
219 GenParticle* elec =
220 new GenParticle( CLHEP::HepLorentzVector( p2[1], p2[2], p2[3], p2[0] ), 11, 3 );
221 elec->suggest_barcode( 2 );
222 prod_vtx->add_particle_in( elec );
223
224 evt->set_beam_particles( posi, elec );
225
226 // outcoming beam e+
227 GenParticle* p =
228 new GenParticle( CLHEP::HepLorentzVector( q1[1], q1[2], q1[3], q1[0] ), -11, 1 );
229 p->suggest_barcode( 3 );
230 prod_vtx->add_particle_out( p );
231
232 // outcoming beam e-
233 p = new GenParticle( CLHEP::HepLorentzVector( q2[1], q2[2], q2[3], q2[0] ), 11, 1 );
234 p->suggest_barcode( 4 );
235 prod_vtx->add_particle_out( p );
236
237 switch ( m_finalState )
238 {
239 case 0: // e+e-pi0 pi0
240 GET_FOURMOMENTA_TWOPI( pi1, pi2 );
241 p = new GenParticle( CLHEP::HepLorentzVector( pi1[1], pi1[2], pi1[3], pi1[0] ), 111, 1 );
242 p->suggest_barcode( 5 );
243 prod_vtx->add_particle_out( p );
244 p = new GenParticle( CLHEP::HepLorentzVector( pi2[1], pi2[2], pi2[3], pi2[0] ), 111, 1 );
245 p->suggest_barcode( 6 );
246 prod_vtx->add_particle_out( p );
247 break;
248 case 1: // e+e-pi+ pi-
249 GET_FOURMOMENTA_TWOPI( pi1, pi2 );
250 p = new GenParticle( CLHEP::HepLorentzVector( pi1[1], pi1[2], pi1[3], pi1[0] ), 211, 1 );
251 p->suggest_barcode( 5 );
252 prod_vtx->add_particle_out( p );
253 p = new GenParticle( CLHEP::HepLorentzVector( pi2[1], pi2[2], pi2[3], pi2[0] ), -211, 1 );
254 p->suggest_barcode( 6 );
255 prod_vtx->add_particle_out( p );
256 double weights[6]; // event weight + five weights for the various contributions to e+e- -->
257 // e+e-pi+pi-
258 for ( int i = 0; i < 6; i++ ) { weights[i] = 0.0; }
259 GET_TWOPI_WEIGHTS( weights );
260 for ( int ai = 0; ai < ( sizeof( weights ) / sizeof( double ) ); ai++ )
261 { evt->weights().push_back( weights[ai] ); }
262 p = new GenParticle( CLHEP::HepLorentzVector( 0, 0, 0, weights[0] ), 99, 1 );
263 p->suggest_barcode( 7 );
264 prod_vtx->add_particle_out( p );
265
266 break;
267 case 2: // e+e- pi0
268 GET_FOURMOMENTA_PION( qpion );
269 p = new GenParticle( CLHEP::HepLorentzVector( qpion[1], qpion[2], qpion[3], qpion[0] ),
270 111, 1 );
271 p->suggest_barcode( 5 );
272 prod_vtx->add_particle_out( p );
273 break;
274 case 3: // e+e- eta
275 GET_FOURMOMENTA_PION( qpion );
276 p = new GenParticle( CLHEP::HepLorentzVector( qpion[1], qpion[2], qpion[3], qpion[0] ),
277 221, 1 );
278 p->suggest_barcode( 5 );
279 prod_vtx->add_particle_out( p );
280 break;
281 case 4: // e+e- etaPrime
282 GET_FOURMOMENTA_PION( qpion );
283 p = new GenParticle( CLHEP::HepLorentzVector( qpion[1], qpion[2], qpion[3], qpion[0] ),
284 331, 1 );
285 p->suggest_barcode( 5 );
286 prod_vtx->add_particle_out( p );
287 break;
288 case 5: // e+e- chi_c0
289 GET_FOURMOMENTA_CHICJ( qcj );
290 p = new GenParticle( CLHEP::HepLorentzVector( qcj[1], qcj[2], qcj[3], qcj[0] ), 10441, 1 );
291 p->suggest_barcode( 5 );
292 prod_vtx->add_particle_out( p );
293 break;
294 case 6: // e+e- chi_c1
295 GET_FOURMOMENTA_CHICJ( qcj );
296 p = new GenParticle( CLHEP::HepLorentzVector( qcj[1], qcj[2], qcj[3], qcj[0] ), 20443, 1 );
297 p->suggest_barcode( 5 );
298 prod_vtx->add_particle_out( p );
299 break;
300 case 7: // e+e- chi_c2
301 GET_FOURMOMENTA_CHICJ( qcj );
302 p = new GenParticle( CLHEP::HepLorentzVector( qcj[1], qcj[2], qcj[3], qcj[0] ), 445, 1 );
303 p->suggest_barcode( 5 );
304 prod_vtx->add_particle_out( p );
305 break;
306 default: break;
307 }
308
309 if ( m_finalState >= 2 && m_finalState <= 4 && m_switchNLO )
310 {
311 if ( NLOTYPE.EvtPhTyp == 1 )
312 {
313 GET_FOURMOMENTA_PHOTON( kphp );
314 p = new GenParticle( CLHEP::HepLorentzVector( kphp[1], kphp[2], kphp[3], kphp[0] ), 22,
315 1 );
316 p->suggest_barcode( 6 );
317 prod_vtx->add_particle_out( p );
318 }
319 if ( m_nloWithWeights )
320 {
321 double weight = GET_WEIGHT();
322 evt->weights().push_back( weight );
323 //!-Store weigth as particle at rest:
324 if ( NLOTYPE.EvtPhTyp == 1 )
325 {
326 p = new GenParticle( CLHEP::HepLorentzVector( 0, 0, 0, weight ), 99, 1 );
327 p->suggest_barcode( 7 );
328 prod_vtx->add_particle_out( p );
329 }
330 else if ( NLOTYPE.EvtPhTyp == 0 )
331 {
332 p = new GenParticle( CLHEP::HepLorentzVector( 0, 0, 0, weight ), 98, 1 );
333 p->suggest_barcode( 6 );
334 prod_vtx->add_particle_out( p );
335 }
336 }
337 }
338
339 if ( log.level() < MSG::INFO ) { evt->print(); }
340
341 // Check if the McCollection already exists
342 SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
343 if ( anMcCol != 0 )
344 {
345 log << MSG::DEBUG << "Adding McGenEvent to existing collection" << endmsg;
346 McGenEvent* mcEvent = new McGenEvent( evt );
347 anMcCol->push_back( mcEvent );
348 }
349 else
350 {
351 McGenEventCol* mcColl = new McGenEventCol;
352 McGenEvent* mcEvent = new McGenEvent( evt );
353 mcColl->push_back( mcEvent );
354 StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
355 if ( sc != StatusCode::SUCCESS )
356 {
357 log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
358 delete mcColl;
359 delete evt;
360 delete mcEvent;
361 return StatusCode::FAILURE;
362 }
363 else { log << MSG::INFO << "McGenEventCol created " << endmsg; }
364 }
365
366 m_totalEvents++;
367
368 log << MSG::DEBUG << "before execute() return" << endmsg;
369 return StatusCode::SUCCESS;
370}
371
372StatusCode Ekhara::finalize() {
373 MsgStream log( msgSvc(), name() );
374 log << MSG::INFO << "Ekhara in finalize()" << endmsg;
375
376 EKHARA( 1 );
377
378 char xsect[100];
379 if ( m_finalState == 0 )
380 { //! e+e- --> e+e-pi0pi0
381 double tnpfinpar[10];
382 for ( int i = 0; i < 10; i++ ) { tnpfinpar[i] = 0.0; }
383
384 GET_FINAL_MESON_INFO( 0, tnpfinpar );
385
386 log << MSG::INFO << "Total number of iterations = " << tnpfinpar[0] << endmsg;
387 log << MSG::INFO << "Number of MC trials = " << tnpfinpar[1] << endmsg;
388 log << MSG::INFO << "ContribMax_MCloop = " << tnpfinpar[2] << endmsg;
389 log << MSG::INFO << "MC efficiency (evt/eff.iter) = " << m_totalEvents / tnpfinpar[1]
390 << endmsg;
391
392 log << MSG::INFO << "================================" << endmsg;
393 log << MSG::INFO << "From all weights:" << endmsg;
394 double CS = tnpfinpar[3] / tnpfinpar[0];
395 double CSerr =
396 sqrt( fabs( tnpfinpar[4] / tnpfinpar[0] - CS * CS ) / ( tnpfinpar[0] - 1. ) );
397 log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [nb]"
398 << endmsg;
399 log << MSG::INFO << endmsg;
400 log << MSG::INFO << "From UNWEIGHTED events:" << endmsg;
401 if ( tnpfinpar[1] <= 0. )
402 {
403 CS = 0.;
404 CSerr = 0.;
405 }
406 else
407 {
408 CS = tnpfinpar[5] * m_totalEvents / tnpfinpar[1];
409 CSerr = sqrt( CS * ( tnpfinpar[5] - CS ) ) / sqrt( tnpfinpar[1] );
410 }
411 log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [nb]"
412 << endmsg;
413 log << MSG::INFO << "================================" << endmsg;
414 if ( tnpfinpar[7] > 0. )
415 {
416 log << MSG::WARNING << "There were overshootings." << endmsg;
417 log << MSG::WARNING << "The events, corresponding to them were ignored." << endmsg;
418 log << MSG::WARNING << "The unweighted sample may be not reliable!" << endmsg;
419 log << MSG::WARNING << "\nContribution from overshooted events:" << endmsg;
420 sprintf( xsect, " Integral Cross Section = %1.5E +- %1.5E [GeV^-2]", tnpfinpar[8],
421 tnpfinpar[9] );
422 log << MSG::WARNING << xsect << endmsg;
423 sprintf( xsect, " = %1.5E +- %1.5E [nb]",
424 tnpfinpar[8] * tnpfinpar[6], tnpfinpar[9] * tnpfinpar[6] );
425 log << MSG::WARNING << xsect << endmsg;
426 log << MSG::INFO << "================================" << endmsg;
427 }
428 }
429 else if ( m_finalState == 1 )
430 { //! e+e- --> e+e-pi+pi-
431 double pipifinpar[2];
432 for ( int i = 0; i < 2; i++ ) { pipifinpar[i] = 0.0; }
433 GET_FINAL_TWOPI_INFO( pipifinpar );
434
435 double CS = pipifinpar[0] / m_totalEvents;
436 double CSerr =
437 sqrt( ( pipifinpar[1] / m_totalEvents - CS * CS ) / ( m_totalEvents - 1. ) );
438 log << MSG::INFO << "================================" << endmsg;
439 log << MSG::INFO << "Cross Section = " << CS << " +- " << CSerr << " [nb]" << endmsg;
440 log << MSG::INFO << "================================" << endmsg;
441 }
442 else if ( m_finalState > 1 && m_finalState < 5 )
443 {
444 double mfp[10];
445 for ( int i = 0; i < 10; i++ ) { mfp[i] = 0.0; }
446 GET_FINAL_MESON_INFO( 0, mfp );
447 if ( m_switchNLO )
448 {
449 log << MSG::INFO << "================================" << endmsg;
450 log << MSG::INFO << "Events 0ph:" << endmsg;
451 log << MSG::INFO << "-------------------------------- " << endmsg;
452 }
453 log << MSG::INFO << "Total number of iterations = " << mfp[0] << endmsg;
454 if ( !m_switchNLO )
455 {
456 log << MSG::INFO << "Number of MC trials = " << mfp[1] << endmsg;
457 log << MSG::INFO << "ContribMax_MCloop = " << mfp[2] << endmsg;
458 log << MSG::INFO << "MC efficiency (evt/eff.iter) = " << m_totalEvents / mfp[1]
459 << endmsg;
460 }
461 log << MSG::INFO << "================================" << endmsg;
462 log << MSG::INFO << "From all weights:" << endmsg;
463 double CS = mfp[3] / mfp[0];
464 double CSerr = sqrt( fabs( mfp[4] / mfp[0] - CS * CS ) / ( mfp[0] - 1. ) );
465 double sumCS = CS;
466 double sumCSerr = CSerr * CSerr;
467 log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [GeV^-2]"
468 << endmsg;
469 log << MSG::INFO << " = " << CS * mfp[6] << " +- " << CSerr * mfp[6]
470 << " [nb]" << endmsg;
471 log << MSG::INFO << endmsg;
472 if ( m_switchNLO )
473 {
474 for ( int i = 0; i < 10; i++ ) { mfp[i] = 0.0; }
475 GET_FINAL_MESON_INFO( 1, mfp );
476 log << MSG::INFO << "================================" << endmsg;
477 log << MSG::INFO << "Events 1ph:" << endmsg;
478 log << MSG::INFO << "-------------------------------- " << endmsg;
479 log << MSG::INFO << "Total number of iterations = " << mfp[0] << endmsg;
480 log << MSG::INFO << "================================" << endmsg;
481 log << MSG::INFO << "From all weights:" << endmsg;
482 CS = mfp[3] / mfp[0];
483 CSerr = sqrt( fabs( mfp[4] / mfp[0] - CS * CS ) / ( mfp[0] - 1. ) );
484 sumCS += CS;
485 sumCSerr += CSerr * CSerr;
486 log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [GeV^-2]"
487 << endmsg;
488 log << MSG::INFO << " = " << CS * mfp[6] << " +- "
489 << CSerr * mfp[6] << " [nb]" << endmsg;
490 log << MSG::INFO << endmsg;
491 log << MSG::INFO << " -------------------------------- " << endmsg;
492 log << MSG::INFO << " sigma (0ph+1ph events) = " << sumCS << " +- " << sqrt( sumCSerr )
493 << " [GeV^-2]" << endmsg;
494 log << MSG::INFO << " = " << sumCS * mfp[6] << " +- "
495 << sqrt( sumCSerr ) * mfp[6] << " [nb]" << endmsg;
496 log << MSG::INFO << "================================" << endmsg;
497 }
498 else
499 {
500 log << MSG::INFO << "From UNWEIGHTED events:" << endmsg;
501 if ( mfp[1] <= 0. )
502 {
503 CS = 0.;
504 CSerr = 0.;
505 }
506 else
507 {
508 CS = mfp[5] * m_totalEvents / mfp[1];
509 CSerr = sqrt( CS * ( mfp[5] - CS ) ) / sqrt( mfp[1] );
510 }
511 sprintf( xsect, " Integral Cross Section = %1.5E +- %1.5E [GeV^-2]", CS, CSerr );
512 log << MSG::INFO << xsect << endmsg;
513 sprintf( xsect, " = %1.5E +- %1.5E [nb]", CS * mfp[6],
514 CSerr * mfp[6] );
515 log << MSG::INFO << xsect << endmsg;
516 log << MSG::INFO << "================================" << endmsg;
517
518 if ( mfp[7] > 0. )
519 {
520 log << MSG::WARNING << "There were overshootings." << endmsg;
521 log << MSG::WARNING << "The events, corresponding to them were ignored." << endmsg;
522 log << MSG::WARNING << "The unweighted sample may be not reliable!" << endmsg;
523 log << MSG::WARNING << "\nContribution from overshooted events:" << endmsg;
524 sprintf( xsect, " Integral Cross Section = %1.5E +- %1.5E [GeV^-2]", mfp[8], mfp[9] );
525 log << MSG::WARNING << xsect << endmsg;
526 sprintf( xsect, " = %1.5E +- %1.5E [nb]", mfp[8] * mfp[6],
527 mfp[9] * mfp[6] );
528 log << MSG::WARNING << xsect << endmsg;
529 log << MSG::INFO << "================================" << endmsg;
530 }
531 }
532 }
533 else if ( m_finalState > 4 && m_finalState < 8 )
534 {
535 double chicjfinpar[3];
536 for ( int i = 0; i < 3; i++ ) { chicjfinpar[i] = 0.0; }
537 GET_FINAL_CHICJ_INFO( chicjfinpar );
538
539 double CS = chicjfinpar[1] / chicjfinpar[0];
540 double CSerr =
541 sqrt( ( chicjfinpar[2] / chicjfinpar[0] - CS * CS ) / ( chicjfinpar[0] - 1. ) );
542 log << MSG::INFO << "================================" << endmsg;
543 sprintf( xsect, " Cross Section = %1.5E +- %1.5E [nb]", CS, CSerr );
544 log << MSG::INFO << xsect << endmsg;
545 log << MSG::INFO << "================================" << endmsg;
546 }
547
548 log << MSG::DEBUG << "Ekhara has terminated successfully" << endmsg;
549
550 return StatusCode::SUCCESS;
551}
552
553StatusCode Ekhara::checkAndPrintParameters() {
554 MsgStream log( msgSvc(), name() );
555
556 log << MSG::INFO
557 << "=============================================================" << endmsg;
558 log << MSG::INFO << "This is EKHARA, Version 3.1" << endmsg;
559 switch ( m_finalState )
560 {
561 case 0:
562 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi0pi0" << endmsg;
563 break;
564 case 1:
565 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi+pi-" << endmsg;
566 break;
567 case 2: log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi0" << endmsg; break;
568 case 3: log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- eta" << endmsg; break;
569 case 4:
570 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- etaPRIME" << endmsg;
571 break;
572 case 5:
573 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c0" << endmsg;
574 break;
575 case 6:
576 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c1" << endmsg;
577 break;
578 case 7:
579 log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c2" << endmsg;
580 break;
581 default:
582 log << MSG::ERROR << "WRONG channel ID: Process not implemented!" << endmsg;
583 return StatusCode::FAILURE;
584 break;
585 }
586
587 if ( m_finalState == 0 )
588 { //!- pi0pi0
589 if ( m_switchNLO )
590 { //!- Check switches for NLO
591 log << MSG::ERROR << "NLO not implemented for this Process!" << endmsg;
592 return StatusCode::FAILURE;
593 }
594 if ( m_twoPionAmplitudes != 4 )
595 {
596 log << MSG::ERROR
597 << " Neutral Pion Pairs only with full TwoPionAmplitudes == 4 ! \t You used "
598 << m_twoPionAmplitudes << endmsg;
599 return StatusCode::FAILURE;
600 }
601 if ( m_twoPionPhspAlg != 1 )
602 {
603 log << MSG::ERROR
604 << " Neutral Pion Pairs only with full TwoPionPhaseSpaceAlg == 1 ! \t You used "
605 << m_twoPionPhspAlg << endmsg;
606 return StatusCode::FAILURE;
607 }
608 log << MSG::INFO << "\tPhase space generation using algorithm: " << m_twoPionPhspAlg
609 << endmsg;
610 log << MSG::INFO << "\tMatrix element: \t|M_s + M_t + M_2g(full)|^2" << endmsg;
611
612 log << MSG::INFO << "The following conditions are applied:" << endmsg;
613 log << MSG::INFO << "\tPion Momentum: " << m_twoPionThetaMin << " < Theta [deg] < "
614 << m_twoPionThetaMax << endmsg;
615 log << MSG::INFO << "\tMissing Momentum: " << m_twoPionMissThetaMin << " < Theta [deg] < "
616 << m_twoPionMissThetaMax << endmsg;
617 }
618 else if ( m_finalState == 1 )
619 { //!- pi+pi-
620 if ( m_switchNLO )
621 { //!- Check switches for NLO
622 log << MSG::ERROR << "NLO not implemented for this Process!" << endmsg;
623 return StatusCode::FAILURE;
624 }
625 if ( m_twoPionPhspAlg < 0 || m_twoPionPhspAlg > 1 )
626 {
627 log << MSG::ERROR
628 << " Wrong choice of phase space generation algorithm: " << m_twoPionPhspAlg
629 << endmsg;
630 return StatusCode::FAILURE;
631 }
632
633 log << MSG::INFO << "\tPhase space generation using algorithm: " << m_twoPionPhspAlg
634 << endmsg;
635 log << MSG::INFO << "\tMatrix element: ";
636 switch ( m_twoPionAmplitudes )
637 { //!- 1 - s, 2 - s+t, 3 - s+t+g_Born, 4 - s+t+g_Full
638 case 1: log << MSG::INFO << "\t|M_s|^2" << endmsg; break;
639 case 2: log << MSG::INFO << "\t|M_s + M_t|^2" << endmsg; break;
640 case 3: log << MSG::INFO << "\t|M_s + M_t + M_2g(Born)|^2" << endmsg; break;
641 case 4: log << MSG::INFO << "\t|M_s + M_t + M_2g(full)|^2" << endmsg; break;
642 default: break;
643 }
644
645 switch ( m_twoPionFormFactor )
646 { //!- 0 - KS, 1 - GS, 2 - new GS
647 case 0: log << MSG::INFO << "\tUsing KS pion form factor" << endmsg; break;
648 case 1: log << MSG::INFO << "\tUsing GS pion form factor" << endmsg; break;
649 case 2: log << MSG::INFO << "\tUsing new GS pion form factor" << endmsg; break;
650 default:
651 log << MSG::WARNING
652 << "Undefined pion form factor switch! Using new GS pion form factor instead!"
653 << endmsg;
654 m_twoPionFormFactor = 2;
655 break;
656 }
657
658 log << MSG::INFO << "The following conditions are applied:" << endmsg;
659 log << MSG::INFO << "\tPion Momentum: " << m_twoPionThetaMin << " < Theta [deg] < "
660 << m_twoPionThetaMax << endmsg;
661 log << MSG::INFO << "\tMissing Momentum: " << m_twoPionMissThetaMin << " < Theta [deg] < "
662 << m_twoPionMissThetaMax << endmsg;
663 }
664 else if ( m_finalState >= 2 && m_finalState <= 4 )
665 { //!-Single Pseudoscalar Meson
666 if ( m_switchNLO )
667 { //!- Check switches for NLO and VacPol
668 log << MSG::INFO << "\tWith NLO corrections: mgamma = " << m_eps_ph * 0.51099906E-3
669 << endmsg;
670 if ( m_nloWithWeights )
671 { log << MSG::WARNING << "Output will contain weighted events!" << endmsg; }
672 else
673 { log << MSG::WARNING << "Attention! Events will be unweighted for output!" << endmsg; }
674 if ( m_switchVP )
675 {
676 log << MSG::INFO << "\tVacuum polarization included:" << endmsg;
677 log << MSG::INFO << "\thttp://www-com.physik.hu-berlin.de/~fjeger/software.html"
678 << endmsg;
679 }
680 else { log << MSG::INFO << "\tVacuum polarization NOT included:" << endmsg; }
681 }
682 else
683 {
684 log << MSG::INFO << "\tWithout NLO corrections!" << endmsg;
685 if ( m_switchVP )
686 {
687 log << MSG::ERROR << "Vacuum polarization can ONLY be included for NLO!" << endmsg;
688 return StatusCode::FAILURE;
689 }
690 }
691
692 log << MSG::INFO << "The cross section will be calculated" << endmsg;
693 log << MSG::INFO << "\tMatrix element: ";
694 switch ( m_mesonAmplitudes )
695 { //!- 1 - s, 2 - t, 3 - s+t
696 case 1:
697 log << MSG::INFO << "\t|M_s|^2" << endmsg;
698 if ( m_switchNLO == 1 )
699 {
700 log << MSG::ERROR << " NLO not implemented for s-channel!" << endmsg;
701 return StatusCode::FAILURE;
702 }
703 if ( m_TFF_ID == 7 || m_TFF_ID == 8 )
704 {
705 log << MSG::WARNING
706 << "TFF Models 7 and 8 were not compared to data in the only-s-channel "
707 "configuration!"
708 << endmsg;
709 }
710 break;
711 case 2: log << MSG::INFO << "\t|M_t|^2" << endmsg; break;
712 case 3: log << MSG::INFO << "\t|M_s + M_t|^2" << endmsg; break;
713 default: break;
714 }
715
716 log << MSG::INFO << "The following conditions are applied:" << endmsg;
717 log << MSG::INFO << "\te+ : " << m_posiThetaMin << " < Theta [deg] < " << m_posiThetaMax
718 << endmsg;
719 log << MSG::INFO << "\t " << m_posiEnergyMin << " < Energy [GeV] < "
720 << m_posiEnergyMin << endmsg;
721 log << MSG::INFO << "\te- : " << m_elecThetaMin << " < Theta [deg] < " << m_elecThetaMax
722 << endmsg;
723 log << MSG::INFO << "\t " << m_elecEnergyMin << " < Energy [GeV] < "
724 << m_elecEnergyMin << endmsg;
725
726 double cosTagAngleRad = fabs( cos( m_taggingAngle * TMath::DegToRad() ) );
727 switch ( m_taggingMode )
728 { //!- 1:Untagged , 2:Single, 3: double
729 case 1:
730 log << MSG::INFO << "\t Generating Untagged event configuration!" << endmsg;
731 log << MSG::INFO << "\t Accepting only events with: |cos(Theta Lepton)| > "
732 << cosTagAngleRad << "!" << endmsg;
733 break;
734 case 2:
735 log << MSG::INFO << "\t Generating Single Tagged event configuration!" << endmsg;
736 log << MSG::INFO << "\t Accepting only events with:" << endmsg;
737 log << MSG::INFO << "\t\t |cos(Theta e+/-)| > " << cosTagAngleRad
738 << " and |cos(Theta e-/+)| < " << cosTagAngleRad << endmsg;
739 break;
740 case 3:
741 log << MSG::INFO << "\t Generating Double Tagged event configuration!" << endmsg;
742 log << MSG::INFO << "\t Accepting only events with: |cos(Theta Lepton)| < "
743 << cosTagAngleRad << "!" << endmsg;
744 break;
745 case 11:
746 log << MSG::INFO << "\t Generating Untagged event configuration!" << endmsg;
747 log << MSG::INFO << "\t Accepting only events with: Q^2 (Lepton)| < " << m_taggingQsquare
748 << "!" << endmsg;
749 break;
750 case 12:
751 log << MSG::INFO << "\t Generating Single Tagged event configuration!" << endmsg;
752 log << MSG::INFO << "\t Accepting only events with:" << endmsg;
753 log << MSG::INFO << "\t\t Q^2 (e+/-) > " << m_taggingQsquare << " and Q^2 (e-/+) < "
754 << m_taggingQsquare << endmsg;
755 break;
756 case 13:
757 log << MSG::INFO << "\t Generating Double Tagged event configuration!" << endmsg;
758 log << MSG::INFO << "\t Accepting only events with: Q^2 (Lepton) > " << m_taggingQsquare
759 << "!" << endmsg;
760 break;
761 default:
762 if ( m_taggingMode > 0 )
763 {
764 log << MSG::INFO << "\t Unknown tagging mode selected!" << endmsg;
765 log << MSG::INFO << "\t Generating events without tagging configuration!" << endmsg;
766 }
767 m_taggingMode = 0;
768 break;
769 }
770
771 switch ( m_finalState )
772 {
773 case 2: log << MSG::INFO << "\tpi0 : "; break;
774 case 3: log << MSG::INFO << "\teta : "; break;
775 case 4: log << MSG::INFO << "\teta' : "; break;
776 default: break;
777 }
778 log << MSG::INFO << m_mesonThetaMin << " < Theta [deg] < " << m_mesonThetaMax << endmsg;
779 log << MSG::INFO << "\t " << m_mesonEnergyMin << " < Energy [GeV] < "
780 << m_mesonEnergyMin << endmsg;
781 }
782 else if ( m_finalState >= 5 && m_finalState <= 7 )
783 { //!- Chi_c0, Chi_c1, Chi_c2
784 if ( m_switchNLO )
785 { //!- Check switches for NLO and VacPol
786 log << MSG::ERROR << "NLO not implemented for this Process!" << endmsg;
787 return StatusCode::FAILURE;
788 }
789 log << MSG::INFO << "The following conditions are applied:" << endmsg;
790 log << MSG::INFO << "\te+ : " << m_posiThetaMin << " < Theta [deg] < " << m_posiThetaMax
791 << endmsg;
792 log << MSG::INFO << "\t " << m_posiEnergyMin << " < Energy [GeV] < "
793 << m_posiEnergyMin << endmsg;
794 log << MSG::INFO << "\te- : " << m_elecThetaMin << " < Theta [deg] < " << m_elecThetaMax
795 << endmsg;
796 log << MSG::INFO << "\t " << m_elecEnergyMin << " < Energy [GeV] < "
797 << m_elecEnergyMin << endmsg;
798 switch ( m_finalState )
799 {
800 case 5: log << MSG::INFO << "\tchi_c0 : "; break;
801 case 6: log << MSG::INFO << "\tchi_c1 : "; break;
802 case 7: log << MSG::INFO << "\tchi_c2 : "; break;
803 default: break;
804 }
805 log << MSG::INFO << m_chicjThetaMin << " < Theta [deg] < " << m_chicjThetaMax << endmsg;
806 log << MSG::INFO << "\t " << m_chicjEnergyMin << " < Energy [GeV] < "
807 << m_chicjEnergyMin << endmsg;
808 }
809 log << MSG::INFO
810 << "=============================================================" << endmsg;
811
812 return StatusCode::SUCCESS;
813}
ai
double ai[4]
Definition AbsCor.cxx:53
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
sprintf
sprintf(cut, "kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_" "pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
CHANNELSEL
#define CHANNELSEL
Definition EkharaDef.h:13
GET_FOURMOMENTA_CHICJ
#define GET_FOURMOMENTA_CHICJ(qcj)
GET_FINAL_CHICJ_INFO
#define GET_FINAL_CHICJ_INFO(chicjfinpar)
GET_FINAL_MESON_INFO
#define GET_FINAL_MESON_INFO(i, mfp)
GET_FOURMOMENTA_TWOPI
#define GET_FOURMOMENTA_TWOPI(pi1, pi2)
EKHARA
#define EKHARA(i)
Definition EkharaDef.h:57
GET_FOURMOMENTA_PHOTON
#define GET_FOURMOMENTA_PHOTON(kphp)
TAGGINGMODE
#define TAGGINGMODE
Definition EkharaDef.h:34
SWDIAG
#define SWDIAG
Definition EkharaDef.h:20
GET_FINAL_TWOPI_INFO
#define GET_FINAL_TWOPI_INFO(pipifinpar)
GET_FOURMOMENTA_PION
#define GET_FOURMOMENTA_PION(qpion)
NLOTYPE
#define NLOTYPE
Definition EkharaDef.h:48
PIONFFSW
#define PIONFFSW
Definition EkharaDef.h:27
GET_WEIGHT
#define GET_WEIGHT()
Definition EkharaDef.h:70
BOSS_INIT_EKHARA
#define BOSS_INIT_EKHARA(xpar)
Definition EkharaDef.h:63
DIAGNOSE
#define DIAGNOSE()
Definition EkharaDef.h:67
FFPARAMSET
#define FFPARAMSET
Definition EkharaDef.h:41
GET_FOURMOMENTA_LEPTONS
#define GET_FOURMOMENTA_LEPTONS(p1, p2, q1, q2)
GET_TWOPI_WEIGHTS
#define GET_TWOPI_WEIGHTS(weights)
Definition EkharaDef.h:73
pi2
character *LEPTONflag integer iresonances real pi2
Definition Ekhara/src/ekhara/common.h:3
McGenEventCol
ObjectVector< McGenEvent > McGenEventCol
Definition Event/GeneratorObject/include/GeneratorObject/McGenEvent.h:39
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition KarFin.h:34
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
EkharaRandom::setRandomEngine
static void setRandomEngine(CLHEP::HepRandomEngine *randomEngine)
Definition EkharaRandom.cxx:10
Ekhara::execute
StatusCode execute()
Definition Ekhara.cxx:179
Ekhara::Ekhara
Ekhara(const std::string &name, ISvcLocator *pSvcLocator)
Definition Ekhara.cxx:26
Ekhara::initialize
StatusCode initialize()
Definition Ekhara.cxx:103
Ekhara::checkAndPrintParameters
StatusCode checkAndPrintParameters()
Definition Ekhara.cxx:553
Ekhara::finalize
StatusCode finalize()
Definition Ekhara.cxx:372