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