BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Phokhara.cxx
Go to the documentation of this file.
1//*****************************************************************************
2//
3// Phokhara.cxx
4//
5// Algorithm runs event generator Phokhara (hep-ph/0710.4227v1)
6// and stores output to transient store
7//
8// Nov. 2007 A. Zhemchugov, initial version written for BES3
9//*****************************************************************************
10
11// Header for this module:-
12#include "Phokhara.h"
13#include "PhokharaDef.h"
14
15// Framework Related Headers:-
16#include "CLHEP/Vector/LorentzVector.h"
17#include "HepMC/GenEvent.h"
18#include "HepMC/GenParticle.h"
19#include "HepMC/GenVertex.h"
20
21#include "GaudiKernel/ISvcLocator.h"
22#include "GaudiKernel/MsgStream.h"
23#include "GaudiKernel/SmartDataPtr.h"
24
25#include "BesRndmGenSvc/IBesRndmGenSvc.h"
26#include "CLHEP/Random/RandomEngine.h"
27#include "GeneratorObject/McGenEvent.h"
28
29#include "cfortran/cfortran.h"
30
31#include <GaudiKernel/StatusCode.h>
32#include <fstream>
33#include <iostream>
34#include <stdlib.h>
35
36using namespace std;
37using namespace CLHEP;
38DECLARE_COMPONENT( Phokhara )
39Phokhara::Phokhara( const string& name, ISvcLocator* pSvcLocator )
40 : Algorithm( name, pSvcLocator ) {
41 m_initSeed = ( 0 );
42
43 m_tagged = 0;
44
45 declareProperty( "InitialSeed", m_initSeed = 0 );
46
47 declareProperty( "InitialEvents", m_nm = 50000 ); // # of events to determine the maximum
48 declareProperty( "NLO", m_nlo = 1 ); // Born(0), NLO(1)
49 declareProperty( "SoftPhotonCutoff", m_w = 0.0001 ); // soft photon cutoff
50 declareProperty( "Channel", m_pion = 0 ); // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
51 // 2pi+2pi-(3),ppbar(4),nnbar(5),
52 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
53 // Lamb Lambbar->pi-pi+ppbar(9)
54 declareProperty( "FSR", m_fsr = 1 ); // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
55 declareProperty( "FSRNLO", m_fsrnlo = 1 ); // yes(1), no(0)
56 declareProperty( "NarrowRes", m_NarrowRes = 1 ); // none(0) jpsi(1) psip(2)
57 declareProperty( "KaonFormfactor", m_FF_Kaon = 1 ); // constrained (0), unconstrained (1),
58 // Kuhn-Khodjamirian-Bruch (2)
59 declareProperty( "VacuumPolarization", m_ivac = 0 ); // yes(1), no(0)
60 // declareProperty( "TaggedPhotons", m_tagged = 0 ); // tagged photons(0), untagged
61 // photons(1)
62 declareProperty( "PionFormfactor", m_FF_Pion = 0 ); // KS Pionformfactor(0), GS
63 // Pionformfactor old(1) and new(2)
64 declareProperty( "F0model", m_f0_model = 0 ); // f0+f0(600): KK model(0), no structure(1),
65 // no f0+f0(600)(2), f0 KLOE(3)
66 declareProperty( "Ecm", m_E = 3.097 ); // CMS-energy (GeV)
67 declareProperty( "CutTotalInvMass",
68 m_q2min = 0.0 ); // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
69 declareProperty( "CutHadInvMass", m_q2_min_c = 0.0447 ); // minimal inv. mass squared of the
70 // hadrons(muons)(GeV^2)
71 declareProperty( "MaxHadInvMass", m_q2_max_c = 1.06 ); // maximal inv. mass squared of the
72 // hadrons(muons)(GeV^2)
73 declareProperty( "MinPhotonEnergy",
74 m_gmin = 0.05 ); // minimal photon energy/missing energy (GeV)
75 declareProperty( "MinPhotonAngle",
76 m_phot1cut = 0.0 ); // minimal photon angle/missing momentum angle (grad)
77 declareProperty( "MaxPhotonAngle",
78 m_phot2cut = 180.0 ); // maximal photon angle/missing momentum angle (grad)
79 declareProperty( "MinHadronsAngle", m_pi1cut = 0.0 ); // minimal hadrons(muons) angle (grad)
80 declareProperty( "MaxHadronsAngle",
81 m_pi2cut = 180.0 ); // maximal hadrons(muons) angle (grad)
82
83 ievent = 0;
84}
85
86StatusCode Phokhara::initialize() {
87
88 MsgStream log( msgSvc(), name() );
89
90 log << MSG::INFO << "Phokhara initialize" << endmsg;
91
92 int i, s_seed[105];
93 long int k, j;
94
95 if ( m_initSeed == 0 ) // if seed is not set in the jobptions file, set it using
96 // BesRndmGenSvc
97 {
98 IBesRndmGenSvc* p_BesRndmGenSvc;
99 static const bool CREATEIFNOTTHERE( true );
100 StatusCode RndmStatus = service( "BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE );
101 if ( !RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc )
102 {
103 log << MSG::ERROR << " Could not initialize Random Number Service" << endmsg;
104 return RndmStatus;
105 }
106
107 // Save the random number seeds in the event
108 HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine( "PHOKHARA" );
109 const long* s = engine->getSeeds();
110 m_initSeed = s[0];
111 }
112
113 log << MSG::INFO << "Set initial seed " << m_initSeed << endmsg;
114
115 // --- initialize the seed ------
116 /* ifstream seeds("seed.dat");
117 if( seeds.is_open() )
118 {
119 int ii=0;
120 while( ! seeds.eof() )
121 seeds >> s_seed[ii++];
122 }
123 else
124 cerr << "Cannot open file seed.dat for reading" << endl;
125 */
126
127 RLXDINIT( 1, m_initSeed );
128
129 // RLXDRESETF(s_seed);
130 // rlxd_reset(s_seed);
131
132 // --- input parameters ----------------------------
133 MAXIMA.iprint = 0;
134
135 FLAGS.nlo = m_nlo;
136 FLAGS.pion = m_pion;
137 FLAGS.fsr = m_fsr;
138 FLAGS.fsrnlo = m_fsrnlo;
139 FLAGS.ivac = m_ivac;
140 // FLAGS.tagged = m_tagged;
141 FLAGS.FF_pion = m_FF_Pion;
142 FLAGS.f0_model = m_f0_model;
143 FLAGS.FF_kaon = m_FF_Kaon;
144 FLAGS.narr_res = m_NarrowRes;
145
146 CTES.Sp = m_E * m_E;
147 ;
148
149 CUTS.w = m_w;
150 CUTS.q2min = m_q2min;
151 CUTS.q2_min_c = m_q2_min_c;
152 CUTS.q2_max_c = m_q2_max_c;
153 CUTS.gmin = m_gmin;
154 CUTS.phot1cut = m_phot1cut;
155 CUTS.phot2cut = m_phot2cut;
156 CUTS.pi1cut = m_pi1cut;
157 CUTS.pi2cut = m_pi2cut;
158
159 INPUT();
160 // --- print run data ------------------------------
161 cout << "-------------------------------------------------------------" << endl;
162 if ( FLAGS.pion == 0 ) cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
163 else if ( FLAGS.pion == 1 ) cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
164 else if ( FLAGS.pion == 2 )
165 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
166 else if ( FLAGS.pion == 3 ) cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
167 else if ( FLAGS.pion == 4 ) cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
168 else if ( FLAGS.pion == 5 ) cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
169 else if ( FLAGS.pion == 6 ) cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
170 else if ( FLAGS.pion == 7 ) cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
171 else if ( FLAGS.pion == 8 )
172 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
173 else if ( FLAGS.pion == 9 )
174 {
175 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
176 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
177 }
178 else cout << " PHOKHARA 7.0: not yet implemented" << endl;
179
180 // --------------------------------
181 cout << "--------------------------------------------------------------" << endl;
182 printf( " %s %f %s\n", "cms total energy = ", sqrt( CTES.Sp ),
183 " GeV " );
184 // if (FLAGS.tagged == 0) {
185 if ( ( CUTS.gmin / 2.0 / CTES.ebeam ) < 0.0098 )
186 {
187 cerr << " minimal missing energy set to small" << endl;
188 return StatusCode::FAILURE;
189 }
190 printf( " %s %f %s\n", "minimal tagged photon energy = ", CUTS.gmin, " GeV " );
191 printf( " %s %f,%f\n", "angular cuts on tagged photon = ", CUTS.phot1cut,
192 CUTS.phot2cut );
193 /* }
194 else {
195 if((CUTS.gmin/2.0/CTES.ebeam) < 0.0098){
196 cerr << " minimal missing energy set to small" << endl;
197 return StatusCode::FAILURE;
198 }
199 printf(" %s %f %s\n","minimal missing energy = ",CUTS.gmin," GeV ");
200 printf(" %s %f,%f\n","angular cuts on missing momentum = ",CUTS.phot1cut,
201 CUTS.phot2cut);
202 }
203 */
204
205 // --------------------------------
206 if ( FLAGS.pion == 0 )
207 printf( " %s %f,%f\n", "angular cuts on muons = ", CUTS.pi1cut,
208 CUTS.pi2cut );
209 else if ( FLAGS.pion == 4 )
210 printf( " %s %f,%f\n", "angular cuts on protons = ", CUTS.pi1cut,
211 CUTS.pi2cut );
212 else if ( FLAGS.pion == 5 )
213 printf( " %s %f,%f\n", "angular cuts on neutrons = ", CUTS.pi1cut,
214 CUTS.pi2cut );
215 else if ( ( FLAGS.pion == 6 ) || ( FLAGS.pion == 7 ) )
216 printf( " %s %f,%f\n", "angular cuts on kaons = ", CUTS.pi1cut,
217 CUTS.pi2cut );
218 else if ( FLAGS.pion == 9 )
219 printf( " %s %f,%f\n", "angular cuts on pions and protons = ", CUTS.pi1cut,
220 CUTS.pi2cut );
221 else
222 printf( " %s %f,%f\n", "angular cuts on pions = ", CUTS.pi1cut,
223 CUTS.pi2cut );
224
225 // if (FLAGS.tagged == 0) {
226 if ( FLAGS.pion == 0 )
227 printf( " %s %f %s\n", "min. muons-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
228 else if ( FLAGS.pion == 4 )
229 printf( " %s %f %s\n", "min. protons-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
230 else if ( FLAGS.pion == 5 )
231 printf( " %s %f %s\n", "min. neutrons-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
232 else if ( ( FLAGS.pion == 6 ) || ( FLAGS.pion == 7 ) )
233 printf( " %s %f %s\n", "min. kaons-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
234 else if ( FLAGS.pion == 9 )
235 printf( " %s %f %s\n", " min. lambdas-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
236 else
237 printf( " %s %f %s\n", "min. pions-tagged photon inv.mass^2 = ", CUTS.q2min, " GeV^2" );
238 // }
239 // --- set cuts ------------------------------------
240 // if (FLAGS.tagged == 0) {
241 cos1min = cos( CUTS.phot2cut * CTES.pi / 180.0 ); // photon1 angle cuts in the
242 cos1max = cos( CUTS.phot1cut * CTES.pi / 180.0 ); // LAB rest frame
243 // } else {
244 // cos1min = -1.0;
245 // cos1max = 1.0;
246 // CUTS.gmin = CUTS.gmin/2.0;
247 // }
248 cos2min = -1.0; // photon2 angle limits
249 cos2max = 1.0; //
250 cos3min = -1.0; // hadrons/muons angle limits
251 cos3max = 1.0; // in their rest frame
252 if ( FLAGS.pion == 0 ) // virtual photon energy cut
253 qqmin = 4.0 * CTES.mmu * CTES.mmu;
254 else if ( FLAGS.pion == 1 ) qqmin = 4.0 * CTES.mpi * CTES.mpi;
255 else if ( FLAGS.pion == 2 )
256 qqmin = 4.0 * ( CTES.mpi + CTES.mpi0 ) * ( CTES.mpi + CTES.mpi0 );
257 else if ( FLAGS.pion == 3 ) qqmin = 16.0 * CTES.mpi * CTES.mpi;
258 else if ( FLAGS.pion == 4 ) qqmin = 4.0 * CTES.mp * CTES.mp;
259 else if ( FLAGS.pion == 5 ) qqmin = 4.0 * CTES.mnt * CTES.mnt;
260 else if ( FLAGS.pion == 6 ) qqmin = 4.0 * CTES.mKp * CTES.mKp;
261 else if ( FLAGS.pion == 7 ) qqmin = 4.0 * CTES.mKn * CTES.mKn;
262 else if ( FLAGS.pion == 8 )
263 qqmin = ( 2.0 * CTES.mpi + CTES.mpi0 ) * ( 2.0 * CTES.mpi + CTES.mpi0 );
264 else if ( FLAGS.pion == 9 ) qqmin = 4.0 * CTES.mlamb * CTES.mlamb;
265 // else
266 // continue;
267
268 qqmax = CTES.Sp - 2.0 * sqrt( CTES.Sp ) * CUTS.gmin; // if only one photon
269 if ( CUTS.q2_max_c < qqmax ) qqmax = CUTS.q2_max_c; // external cuts
270
271 // -------------------
272 if ( ( CUTS.q2_min_c > qqmin ) &&
273 ( CUTS.q2_min_c <
274 ( CTES.Sp * ( 1.0 - 2.0 * ( CUTS.gmin / sqrt( CTES.Sp ) + CUTS.w ) ) ) ) )
275 qqmin = CUTS.q2_min_c;
276 else
277 {
278 cerr << "------------------------------" << endl;
279 cerr << " Q^2_min TOO SMALL" << endl;
280 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
281 cerr << "------------------------------" << endl;
282 }
283 // -------------------
284 if ( qqmax <= qqmin )
285 {
286 cerr << " Q^2_max to small " << endl;
287 cerr << " Q^2_max = " << qqmax << endl;
288 cerr << " Q^2_min = " << qqmin << endl;
289 return StatusCode::FAILURE;
290 }
291
292 // -------------------
293 if ( FLAGS.pion == 0 )
294 {
295 printf( " %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin, " GeV^2" );
296 printf( " %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax, " GeV^2" );
297 }
298 else if ( FLAGS.pion == 1 )
299 {
300 printf( " %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin, " GeV^2" );
301 printf( " %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax, " GeV^2" );
302 }
303 else if ( FLAGS.pion == 4 )
304 {
305 printf( " %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin, " GeV^2" );
306 printf( " %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax, " GeV^2" );
307 }
308 else if ( FLAGS.pion == 5 )
309 {
310 printf( " %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin, " GeV^2" );
311 printf( " %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax, " GeV^2" );
312 }
313 else if ( ( FLAGS.pion == 6 ) || ( FLAGS.pion == 7 ) )
314 {
315 printf( " %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin, " GeV^2" );
316 printf( " %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax, " GeV^2" );
317 }
318 else if ( FLAGS.pion == 8 )
319 {
320 printf( " %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin, " GeV^2" );
321 printf( " %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax, " GeV^2" );
322 }
323 else if ( FLAGS.pion == 9 )
324 {
325 printf( " %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin, " GeV^2" );
326 printf( " %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax, " GeV^2" );
327 }
328 else
329 {
330 printf( " %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin, " GeV^2" );
331 printf( " %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax, " GeV^2" );
332 }
333 // -------------------
334 if ( FLAGS.nlo == 0 )
335 {
336 cout << "Born" << endl;
337 if ( FLAGS.fsrnlo != 0 )
338 {
339 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
340 return StatusCode::FAILURE;
341 }
342 }
343 // -------------------
344 if ( ( FLAGS.pion == 9 ) && ( FLAGS.nlo != 0 ) )
345 {
346 cerr << "WRONG NLO flag - only Born allowed for Lambdas" << endl;
347 cerr << "If you feel that you need NLO" << endl;
348 cerr << "please contact the authors" << endl;
349 return StatusCode::FAILURE;
350 }
351 // -------------------
352 if ( FLAGS.nlo == 1 )
353 printf( " %s %f\n", "NLO: soft photon cutoff w = ", CUTS.w );
354 if ( ( FLAGS.pion <= 1 ) || ( FLAGS.pion == 6 ) )
355 {
356
357 if ( !( ( FLAGS.fsr == 1 ) || ( FLAGS.fsr == 2 ) || ( FLAGS.fsrnlo == 0 ) ||
358 ( FLAGS.fsr == 1 ) || ( FLAGS.fsrnlo == 1 ) || ( FLAGS.fsr == 0 ) ||
359 ( FLAGS.fsrnlo == 0 ) ) )
360 {
361 cerr << "WRONG combination of FSR, FSRNLO flags" << endl;
362 return StatusCode::FAILURE;
363 }
364
365 // ------------------
366 if ( FLAGS.fsr == 0 ) cout << "ISR only" << endl;
367 else if ( FLAGS.fsr == 1 ) cout << "ISR+FSR" << endl;
368 else if ( FLAGS.fsr == 2 )
369 {
370 if ( FLAGS.nlo == 0 ) cout << "ISR+INT+FSR" << endl;
371 else
372 {
373 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
374 return StatusCode::FAILURE;
375 }
376 }
377 else
378 {
379 cerr << "WRONG FSR flag" << FLAGS.fsr << endl;
380 return StatusCode::FAILURE;
381 }
382
383 if ( FLAGS.fsrnlo == 1 ) cout << "IFSNLO included" << endl;
384 }
385 else
386 {
387 if ( ( FLAGS.fsr == 0 ) && ( FLAGS.fsrnlo == 0 ) ) cout << "ISR only" << endl;
388 else
389 {
390 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
391 return StatusCode::FAILURE;
392 }
393 }
394
395 // ------------------
396 if ( FLAGS.ivac == 0 ) { cout << "Vacuum polarization is NOT included" << endl; }
397 else if ( FLAGS.ivac == 1 )
398 {
399 cout << "Vacuum polarization by Fred Jegerlehner "
400 "(http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)"
401 << endl;
402 }
403 else if ( FLAGS.ivac == 2 )
404 {
405 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner"
406 << endl;
407 }
408 else
409 {
410 cout << "WRONG vacuum polarization switch" << endl;
411 return StatusCode::FAILURE;
412 }
413
414 // -----------------
415 if ( FLAGS.pion == 1 )
416 {
417 if ( FLAGS.FF_pion == 0 ) cout << "Kuhn-Santamaria PionFormFactor" << endl;
418 else if ( FLAGS.FF_pion == 1 ) cout << "Gounaris-Sakurai PionFormFactor old" << endl;
419 else if ( FLAGS.FF_pion == 2 ) cout << "Gounaris-Sakurai PionFormFactor new" << endl;
420 else
421 {
422 cout << "WRONG PionFormFactor switch" << endl;
423 return StatusCode::FAILURE;
424 }
425
426 if ( FLAGS.fsr != 0 )
427 {
428 if ( FLAGS.f0_model == 0 ) cout << "f0+f0(600): K+K- model" << endl;
429 else if ( FLAGS.f0_model == 1 ) cout << "f0+f0(600): \"no structure\" model" << endl;
430 else if ( FLAGS.f0_model == 2 ) cout << "NO f0+f0(600)" << endl;
431 else if ( FLAGS.f0_model == 3 )
432 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
433 else
434 {
435 cout << "WRONG f0+f0(600) switch" << endl;
436 return StatusCode::FAILURE;
437 }
438 }
439 }
440
441 //-------
442 if ( ( FLAGS.pion == 6 ) || ( FLAGS.pion == 7 ) )
443 {
444 if ( FLAGS.FF_kaon == 0 ) cout << "constrained KaonFormFactor" << endl;
445 else if ( FLAGS.FF_kaon == 1 ) { cout << "unconstrained KaonFormFactor" << endl; }
446 else if ( FLAGS.FF_kaon == 2 )
447 { cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl; }
448 else
449 {
450 cout << "WRONG KaonFormFactor switch" << endl;
451 return StatusCode::FAILURE;
452 }
453 }
454 // --------------------
455 if ( ( FLAGS.pion == 0 ) || ( FLAGS.pion == 1 ) || ( FLAGS.pion == 6 ) ||
456 ( FLAGS.pion == 7 ) )
457 {
458 if ( FLAGS.narr_res == 1 ) { cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl; }
459 else if ( FLAGS.narr_res == 2 )
460 { cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl; }
461 else if ( FLAGS.narr_res != 0 )
462 {
463 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
464 return StatusCode::FAILURE;
465 }
466 }
467 // ------
468
469 cout << "--------------------------------------------------------------" << endl;
470 //
471 // =================================================
472 // --- finding the maximum -------------------------
473 for ( i = 0; i < 2; i++ )
474 {
475 MAXIMA.Mmax[i] = 1.0;
476 MAXIMA.gross[i] = 0.0;
477 MAXIMA.klein[i] = 0.0;
478 }
479
480 if ( FLAGS.nlo == 0 ) MAXIMA.Mmax[1] = 0.0; // only 1 photon events generated
481
482 MAXIMA.tr[0] = 0.0;
483 MAXIMA.tr[1] = 0.0;
484 MAXIMA.count[0] = 0.0;
485 MAXIMA.count[1] = 0.0;
486
487 // =================================================
488 // --- beginning the MC loop event generation ------
489 for ( j = 1; j <= m_nm; j++ )
490 {
491 RANLXDF( Ar_r, 1 );
492 Ar[1] = Ar_r[0];
493 if ( Ar[1] <= ( MAXIMA.Mmax[0] / ( MAXIMA.Mmax[0] + MAXIMA.Mmax[1] ) ) )
494 {
495 MAXIMA.count[0] = MAXIMA.count[0] + 1.0;
496 GEN_1PH( 1, qqmin, qqmax, cos1min, cos1max, cos3min, cos3max );
497 }
498 else
499 {
500 MAXIMA.count[1] = MAXIMA.count[1] + 1.0;
501 GEN_2PH( 1, qqmin, cos1min, cos1max, cos2min, cos2max, cos3min, cos3max );
502 }
503 }
504 // --- end of the MC loop --------------------------
505 // =================================================
506 // --- for the second run ---
507 MAXIMA.Mmax[0] = MAXIMA.gross[0] + .05 * sqrt( MAXIMA.gross[0] * MAXIMA.gross[0] );
508 MAXIMA.Mmax[1] =
509 MAXIMA.gross[1] + ( .03 + .02 * CTES.Sp ) * sqrt( MAXIMA.gross[1] * MAXIMA.gross[1] );
510 if ( ( FLAGS.pion == 1 ) && ( FLAGS.fsrnlo == 1 ) ) MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * 1.5;
511 if ( ( FLAGS.pion == 0 ) && ( FLAGS.fsrnlo == 1 ) ) MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * 1.5;
512
513 if ( ( FLAGS.pion == 0 ) && ( FLAGS.fsr == 1 ) && ( FLAGS.fsrnlo == 0 ) )
514 MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * 1.2;
515
516 if ( ( FLAGS.pion == 2 ) || ( FLAGS.pion == 3 ) )
517 {
518 MAXIMA.Mmax[0] = MAXIMA.Mmax[0] * 1.1;
519 MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * 1.1;
520 }
521
522 if ( FLAGS.pion == 8 )
523 {
524 MAXIMA.Mmax[0] = MAXIMA.Mmax[0] * 1.08;
525 MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * 1.1;
526 }
527 // --- end of the second run -----------------------
528
529 MAXIMA.tr[0] = 0.0;
530 MAXIMA.tr[1] = 0.0;
531 MAXIMA.count[0] = 0.0;
532 MAXIMA.count[1] = 0.0;
533
534 return StatusCode::SUCCESS;
535}
536
537StatusCode Phokhara::execute() {
538 MsgStream log( msgSvc(), name() );
539
540 int ntrials = 0;
541 int tr_old[2];
542 tr_old[0] = (int)MAXIMA.tr[0];
543 tr_old[1] = (int)MAXIMA.tr[1];
544
545 while ( ntrials < 10000 )
546 {
547 ievent++;
548 RANLXDF( Ar_r, 1 );
549 Ar[1] = Ar_r[0];
550
551 if ( Ar[1] <= ( MAXIMA.Mmax[0] / ( MAXIMA.Mmax[0] + MAXIMA.Mmax[1] ) ) )
552 {
553 MAXIMA.count[0] = MAXIMA.count[0] + 1.0;
554 GEN_1PH( 2, qqmin, qqmax, cos1min, cos1max, cos3min, cos3max );
555 }
556 else
557 {
558 MAXIMA.count[1] = MAXIMA.count[1] + 1.0;
559 GEN_2PH( 2, qqmin, cos1min, cos1max, cos2min, cos2max, cos3min, cos3max );
560 }
561
562 if ( ( (int)MAXIMA.tr[0] + (int)MAXIMA.tr[1] ) >
563 ( tr_old[0] + tr_old[1] ) ) // event accepted after cuts
564 {
565 storeParticles().ignore();
566 return StatusCode::SUCCESS;
567 }
568 ntrials++;
569 }
570
571 log << MSG::FATAL << "Could not satisfy cuts after " << ntrials << "trials. Terminate."
572 << endmsg;
573
574 return StatusCode::FAILURE;
575}
576
577StatusCode Phokhara::storeParticles() {
578 MsgStream log( msgSvc(), name() );
579 int npart = 0;
580
581 // Fill event information
582 GenEvent* evt = new GenEvent( 1, 1 );
583
584 GenVertex* prod_vtx = new GenVertex();
585 // prod_vtx->add_particle_out( p );
586 evt->add_vertex( prod_vtx );
587
588 // incoming beam e+
589 GenParticle* p =
590 new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][0], CTES.momenta[2][0],
591 CTES.momenta[3][0], CTES.momenta[0][0] ),
592 -11, 3 );
593 p->suggest_barcode( ++npart );
594 prod_vtx->add_particle_in( p );
595
596 // incoming beam e-
597 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][1], CTES.momenta[2][1],
598 CTES.momenta[3][1], CTES.momenta[0][1] ),
599 11, 3 );
600 p->suggest_barcode( ++npart );
601 prod_vtx->add_particle_in( p );
602
603 // gamma
604 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][2], CTES.momenta[2][2],
605 CTES.momenta[3][2], CTES.momenta[0][2] ),
606 22, 1 );
607 p->suggest_barcode( ++npart );
608 prod_vtx->add_particle_out( p );
609
610 if ( CTES.momenta[0][3] != 0 ) // second photon exists
611 {
612 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][3], CTES.momenta[2][3],
613 CTES.momenta[3][3], CTES.momenta[0][3] ),
614 22, 1 );
615 p->suggest_barcode( ++npart );
616 prod_vtx->add_particle_out( p );
617 }
618
619 // other products depending on channel
620 if ( FLAGS.pion == 0 )
621 { // mu+ mu-
622 // mu+
623 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
624 CTES.momenta[3][5], CTES.momenta[0][5] ),
625 -13, 1 );
626 p->suggest_barcode( ++npart );
627 prod_vtx->add_particle_out( p );
628 // mu -
629 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
630 CTES.momenta[3][6], CTES.momenta[0][6] ),
631 13, 1 );
632 p->suggest_barcode( ++npart );
633 prod_vtx->add_particle_out( p );
634 }
635
636 if ( FLAGS.pion == 1 )
637 { // pi+ pi-
638 // pi+
639 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
640 CTES.momenta[3][5], CTES.momenta[0][5] ),
641 211, 1 );
642 p->suggest_barcode( ++npart );
643 prod_vtx->add_particle_out( p );
644 // pi-
645 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
646 CTES.momenta[3][6], CTES.momenta[0][6] ),
647 -211, 1 );
648 p->suggest_barcode( ++npart );
649 prod_vtx->add_particle_out( p );
650 }
651
652 if ( FLAGS.pion == 2 )
653 { // pi+ pi- pi0 pi0
654 // pi0
655 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
656 CTES.momenta[3][5], CTES.momenta[0][5] ),
657 111, 1 );
658 p->suggest_barcode( ++npart );
659 prod_vtx->add_particle_out( p );
660 // pi0
661 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
662 CTES.momenta[3][6], CTES.momenta[0][6] ),
663 111, 1 );
664 p->suggest_barcode( ++npart );
665 prod_vtx->add_particle_out( p );
666 // pi-
667 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7],
668 CTES.momenta[3][7], CTES.momenta[0][7] ),
669 -211, 1 );
670 p->suggest_barcode( ++npart );
671 prod_vtx->add_particle_out( p );
672 // pi+
673 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8],
674 CTES.momenta[3][8], CTES.momenta[0][8] ),
675 211, 1 );
676 p->suggest_barcode( ++npart );
677 prod_vtx->add_particle_out( p );
678 }
679
680 if ( FLAGS.pion == 3 )
681 { // pi+ pi- pi+ pi-
682 // pi+
683 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
684 CTES.momenta[3][5], CTES.momenta[0][5] ),
685 211, 1 );
686 p->suggest_barcode( ++npart );
687 prod_vtx->add_particle_out( p );
688 // pi-
689 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
690 CTES.momenta[3][6], CTES.momenta[0][6] ),
691 -211, 1 );
692 p->suggest_barcode( ++npart );
693 prod_vtx->add_particle_out( p );
694 // pi-
695 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7],
696 CTES.momenta[3][7], CTES.momenta[0][7] ),
697 -211, 1 );
698 p->suggest_barcode( ++npart );
699 prod_vtx->add_particle_out( p );
700 // pi+
701 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8],
702 CTES.momenta[3][8], CTES.momenta[0][8] ),
703 211, 1 );
704 p->suggest_barcode( ++npart );
705 prod_vtx->add_particle_out( p );
706 }
707
708 if ( FLAGS.pion == 4 )
709 { // p pbar
710 // pbar
711 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
712 CTES.momenta[3][5], CTES.momenta[0][5] ),
713 -2212, 1 );
714 p->suggest_barcode( ++npart );
715 prod_vtx->add_particle_out( p );
716 // p
717 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
718 CTES.momenta[3][6], CTES.momenta[0][6] ),
719 2212, 1 );
720 p->suggest_barcode( ++npart );
721 prod_vtx->add_particle_out( p );
722 }
723
724 if ( FLAGS.pion == 5 )
725 { // n nbar
726 // nbar
727 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
728 CTES.momenta[3][5], CTES.momenta[0][5] ),
729 -2112, 1 );
730 p->suggest_barcode( ++npart );
731 prod_vtx->add_particle_out( p );
732 // n
733 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
734 CTES.momenta[3][6], CTES.momenta[0][6] ),
735 2112, 1 );
736 p->suggest_barcode( ++npart );
737 prod_vtx->add_particle_out( p );
738 }
739
740 if ( FLAGS.pion == 6 )
741 { // K+ K-
742 // pbar
743 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
744 CTES.momenta[3][5], CTES.momenta[0][5] ),
745 321, 1 );
746 p->suggest_barcode( ++npart );
747 prod_vtx->add_particle_out( p );
748 // p
749 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
750 CTES.momenta[3][6], CTES.momenta[0][6] ),
751 -321, 1 );
752 p->suggest_barcode( ++npart );
753 prod_vtx->add_particle_out( p );
754 }
755
756 if ( FLAGS.pion == 7 )
757 { // K0 K0bar
758 // pbar
759 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
760 CTES.momenta[3][5], CTES.momenta[0][5] ),
761 311, 1 );
762 p->suggest_barcode( ++npart );
763 prod_vtx->add_particle_out( p );
764 // p
765 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
766 CTES.momenta[3][6], CTES.momenta[0][6] ),
767 -311, 1 );
768 p->suggest_barcode( ++npart );
769 prod_vtx->add_particle_out( p );
770 }
771
772 if ( FLAGS.pion == 8 )
773 { // pi+ pi- pi0
774 // pi+
775 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
776 CTES.momenta[3][5], CTES.momenta[0][5] ),
777 211, 1 );
778 p->suggest_barcode( ++npart );
779 prod_vtx->add_particle_out( p );
780 // pi-
781 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
782 CTES.momenta[3][6], CTES.momenta[0][6] ),
783 -211, 1 );
784 p->suggest_barcode( ++npart );
785 prod_vtx->add_particle_out( p );
786 // pi0
787 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7],
788 CTES.momenta[3][7], CTES.momenta[0][7] ),
789 111, 1 );
790 p->suggest_barcode( ++npart );
791 prod_vtx->add_particle_out( p );
792 }
793
794 if ( FLAGS.pion == 9 )
795 { // Lambda Lambdabar Pi+ Pi- P Pbar
796 // Lambdabar
797 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5],
798 CTES.momenta[3][5], CTES.momenta[0][5] ),
799 -3122, 2 );
800 p->suggest_barcode( ++npart );
801 prod_vtx->add_particle_out( p );
802 // Lambda
803 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6],
804 CTES.momenta[3][6], CTES.momenta[0][6] ),
805 3122, 2 );
806 p->suggest_barcode( ++npart );
807 prod_vtx->add_particle_out( p );
808 // pi+
809 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7],
810 CTES.momenta[3][7], CTES.momenta[0][7] ),
811 211, 1 );
812 p->suggest_barcode( ++npart );
813 prod_vtx->add_particle_out( p );
814 // pbar
815 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8],
816 CTES.momenta[3][8], CTES.momenta[0][8] ),
817 -2212, 1 );
818 p->suggest_barcode( ++npart );
819 prod_vtx->add_particle_out( p );
820 // pi-
821 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][9], CTES.momenta[2][9],
822 CTES.momenta[3][9], CTES.momenta[0][9] ),
823 -211, 1 );
824 p->suggest_barcode( ++npart );
825 prod_vtx->add_particle_out( p );
826 // p
827 p = new GenParticle( CLHEP::HepLorentzVector( CTES.momenta[1][10], CTES.momenta[2][10],
828 CTES.momenta[3][10], CTES.momenta[0][10] ),
829 2212, 1 );
830 p->suggest_barcode( ++npart );
831 prod_vtx->add_particle_out( p );
832 }
833
834 if ( log.level() < MSG::INFO ) { evt->print(); }
835
836 // Check if the McCollection already exists
837 SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
838 if ( anMcCol != 0 )
839 {
840 // Add event to existing collection
841 MsgStream log( msgSvc(), name() );
842 log << MSG::INFO << "Add McGenEvent to existing collection" << endmsg;
843 McGenEvent* mcEvent = new McGenEvent( evt );
844 anMcCol->push_back( mcEvent );
845 }
846 else
847 {
848 // Create Collection and add to the transient store
849 McGenEventCol* mcColl = new McGenEventCol;
850 McGenEvent* mcEvent = new McGenEvent( evt );
851 mcColl->push_back( mcEvent );
852 StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
853 if ( sc != StatusCode::SUCCESS )
854 {
855 log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
856 delete mcColl;
857 delete evt;
858 delete mcEvent;
859 return StatusCode::FAILURE;
860 }
861 else
862 {
863 log << MSG::INFO << "McGenEventCol created and " << npart
864 << " particles stored in McGenEvent" << endmsg;
865 }
866 }
867
868 return StatusCode::SUCCESS;
869}
870
871StatusCode Phokhara::finalize() {
872 MsgStream log( msgSvc(), name() );
873 // =================================================
874 if ( FLAGS.pion == 9 )
875 MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * ( 1.0 + LAMBDA_PAR.alpha_lamb ) *
876 ( 1.0 + LAMBDA_PAR.alpha_lamb ) * LAMBDA_PAR.ratio_lamb *
877 LAMBDA_PAR.ratio_lamb;
878
879 // --- value of the cross section ------------------
880 if ( FLAGS.nlo == 0 )
881 {
882 sigma = MAXIMA.Mmax[0] / MAXIMA.count[0] * MAXIMA.tr[0];
883 dsigm = MAXIMA.Mmax[0] *
884 sqrt( ( MAXIMA.tr[0] / MAXIMA.count[0] -
885 ( MAXIMA.tr[0] / MAXIMA.count[0] ) * ( MAXIMA.tr[0] / MAXIMA.count[0] ) ) /
886 MAXIMA.count[0] );
887 }
888 else
889 {
890 sigma1 = MAXIMA.Mmax[0] / MAXIMA.count[0] * MAXIMA.tr[0];
891 sigma2 = MAXIMA.Mmax[1] / MAXIMA.count[1] * MAXIMA.tr[1];
892 dsigm1 =
893 MAXIMA.Mmax[0] *
894 sqrt( ( MAXIMA.tr[0] / MAXIMA.count[0] -
895 ( MAXIMA.tr[0] / MAXIMA.count[0] ) * ( MAXIMA.tr[0] / MAXIMA.count[0] ) ) /
896 MAXIMA.count[0] );
897 dsigm2 =
898 MAXIMA.Mmax[1] *
899 sqrt( ( MAXIMA.tr[1] / MAXIMA.count[1] -
900 ( MAXIMA.tr[1] / MAXIMA.count[1] ) * ( MAXIMA.tr[1] / MAXIMA.count[1] ) ) /
901 MAXIMA.count[1] );
902
903 sigma = sigma1 + sigma2;
904 dsigm = sqrt( dsigm1 * dsigm1 + dsigm2 * dsigm2 );
905 }
906
907 // --- output --------------------------------------
908 cout << "-------------------------------------------------------------" << endl;
909 cout << " PHOKHARA 7.0 Final Statistics " << endl;
910 cout << "-------------------------------------------------------------" << endl;
911 cout << int( MAXIMA.tr[0] + MAXIMA.tr[1] ) << " total events accepted of " << endl;
912 cout << int( ievent ) << " total events generated" << endl;
913 cout << int( MAXIMA.tr[0] ) << " one photon events accepted of " << endl;
914 cout << int( MAXIMA.count[0] ) << " events generated" << endl;
915 cout << int( MAXIMA.tr[1] ) << " two photon events accepted of " << endl;
916 cout << int( MAXIMA.count[1] ) << " events generated" << endl;
917 cout << endl;
918 if ( FLAGS.nlo != 0 ) cout << "sigma1(nbarn) = " << sigma1 << " +- " << dsigm1 << endl;
919 if ( FLAGS.nlo != 0 ) cout << "sigma2(nbarn) = " << sigma2 << " +- " << dsigm2 << endl;
920 cout << "sigma (nbarn) = " << sigma << " +- " << dsigm << endl;
921 cout << endl;
922 cout << "maximum1 = " << MAXIMA.gross[0] << " minimum1 = " << MAXIMA.klein[0] << endl;
923 cout << "Mmax1 = " << MAXIMA.Mmax[0] << endl;
924 cout << "maximum2 = " << MAXIMA.gross[1] << " minimum2 = " << MAXIMA.klein[1] << endl;
925 cout << "Mmax2 = " << MAXIMA.Mmax[1] << endl;
926 cout << "-------------------------------------------------------------" << endl;
927
928 log << MSG::INFO << "Phokhara finalized" << endmsg;
929
930 return StatusCode::SUCCESS;
931}
INPUT
#define INPUT
Definition BesBdkRc.cxx:31
FLAGS
#define FLAGS
Definition BesBdkRc.cxx:78
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
McGenEventCol
ObjectVector< McGenEvent > McGenEventCol
Definition Event/GeneratorObject/include/GeneratorObject/McGenEvent.h:39
RLXDINIT
#define RLXDINIT(LUXURY, SEED)
Definition EvtPhokharaDef.hh:88
GEN_1PH
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
Definition EvtPhokharaDef.hh:112
RANLXDF
#define RANLXDF(AR, VAL)
Definition EvtPhokharaDef.hh:109
GEN_2PH
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
s
XmlRpcServer s
Definition HelloServer.cpp:11
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
CTES
#define CTES
Definition PhokharaDef.h:22
LAMBDA_PAR
#define LAMBDA_PAR
Definition PhokharaDef.h:59
MAXIMA
#define MAXIMA
Definition PhokharaDef.h:82
CUTS
#define CUTS
Definition PhokharaDef.h:31
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
IBesRndmGenSvc
manage multiple CLHEP random engines as named streams
Definition Control/BesRndmGenSvc/include/BesRndmGenSvc/IBesRndmGenSvc.h:23
IBesRndmGenSvc::GetEngine
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.
Phokhara::Phokhara
Phokhara(const string &name, ISvcLocator *pSvcLocator)
Definition Phokhara.cxx:39
Phokhara::initialize
StatusCode initialize()
Definition Phokhara.cxx:86
Phokhara::storeParticles
StatusCode storeParticles()
Definition Phokhara.cxx:577
Phokhara::finalize
StatusCode finalize()
Definition Phokhara.cxx:871
Phokhara::execute
StatusCode execute()
Definition Phokhara.cxx:537