BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of models developed at BES collaboration
5// based on the EvtGen framework. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/BesCopyright
9// Copyright (A) 2006 Ping Rong-Gang @IHEP
10//
11// Module: EvtPhokhara.cc
12// the necessary file: phokharar.F
13// fist.inc,gen.inc mix.inc stdhep.inc
14// Description: Modified Lund model at tau-charm energy level, see
15// PHYSICAL REVIEW D, VOLUME 62, 034003
16// Modification history:
17//
18// Ping R.-G. Jan.25, 2010 Module created
19// The random engine RLU0 is unified with RLU for BesEvtGen
20//------------------------------------------------------------------------
21//
22#include "EvtPhokhara.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtId.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include "../EvtGenBase/EvtStringParticle.hh"
30#include "EvtPhokharaDef.hh"
31#include <fstream>
32#include <iomanip>
33#include <iostream>
34#include <stdio.h>
35#include <stdlib.h>
36#include <string.h>
37#include <string>
38#include <unistd.h>
39
40#include "BesRndmGenSvc/IBesRndmGenSvc.h"
41#include "CLHEP/Random/RandomEngine.h"
42#include "GeneratorObject/McGenEvent.h"
43#include "cfortran/cfortran.h"
44
45using namespace std;
46using namespace CLHEP;
47
48using std::endl;
49using std::fstream;
50using std::ios;
51using std::ofstream;
52using std::resetiosflags;
53using std::setiosflags;
54using std::setw;
55
56int EvtPhokhara::nphokharadecays = 0;
57EvtDecayBasePtr* EvtPhokhara::phokharadecays = 0;
58int EvtPhokhara::ntable = 0;
59
60int EvtPhokhara::ncommand = 0;
61int EvtPhokhara::lcommand = 0;
62std::string* EvtPhokhara::commands = 0;
63int EvtPhokhara::nevt = 0;
64
65EvtPhokhara::EvtPhokhara() {}
66EvtPhokhara::~EvtPhokhara() {
67
68 // print out message
69 // =================================================
70 if ( flags_.pion == 9 )
71 maxima_.Mmax[1] = maxima_.Mmax[1] * ( 1.0 + lambda_par_.alpha_lamb ) *
72 ( 1.0 + lambda_par_.alpha_lamb ) * lambda_par_.ratio_lamb *
73 lambda_par_.ratio_lamb;
74
75 // --- value of the cross section ------------------
76 if ( flags_.nlo == 0 )
77 {
78 sigma = maxima_.Mmax[0] / maxima_.count[0] * maxima_.tr[0];
79 dsigm =
80 maxima_.Mmax[0] *
81 sqrt( ( maxima_.tr[0] / maxima_.count[0] -
82 ( maxima_.tr[0] / maxima_.count[0] ) * ( maxima_.tr[0] / maxima_.count[0] ) ) /
83 maxima_.count[0] );
84 }
85 else
86 {
87 sigma1 = maxima_.Mmax[0] / maxima_.count[0] * maxima_.tr[0];
88 sigma2 = maxima_.Mmax[1] / maxima_.count[1] * maxima_.tr[1];
89 dsigm1 =
90 maxima_.Mmax[0] *
91 sqrt( ( maxima_.tr[0] / maxima_.count[0] -
92 ( maxima_.tr[0] / maxima_.count[0] ) * ( maxima_.tr[0] / maxima_.count[0] ) ) /
93 maxima_.count[0] );
94 dsigm2 =
95 maxima_.Mmax[1] *
96 sqrt( ( maxima_.tr[1] / maxima_.count[1] -
97 ( maxima_.tr[1] / maxima_.count[1] ) * ( maxima_.tr[1] / maxima_.count[1] ) ) /
98 maxima_.count[1] );
99
100 sigma = sigma1 + sigma2;
101 dsigm = sqrt( dsigm1 * dsigm1 + dsigm2 * dsigm2 );
102 }
103 // --- output --------------------------------------
104 cout << "-------------------------------------------------------------" << endl;
105 cout << " PHOKHARA 7.0 Final Statistics " << endl;
106 cout << "-------------------------------------------------------------" << endl;
107 cout << int( maxima_.tr[0] + maxima_.tr[1] ) << " total events accepted of " << endl;
108 cout << int( ievent ) << " total events generated" << endl;
109 cout << int( maxima_.tr[0] ) << " one photon events accepted of " << endl;
110 cout << int( maxima_.count[0] ) << " events generated" << endl;
111 cout << int( maxima_.tr[1] ) << " two photon events accepted of " << endl;
112 cout << int( maxima_.count[1] ) << " events generated" << endl;
113 cout << endl;
114 if ( flags_.nlo != 0 ) cout << "sigma1(nbarn) = " << sigma1 << " +- " << dsigm1 << endl;
115 if ( flags_.nlo != 0 ) cout << "sigma2(nbarn) = " << sigma2 << " +- " << dsigm2 << endl;
116 cout << "sigma (nbarn) = " << sigma << " +- " << dsigm << endl;
117 cout << endl;
118 cout << "maximum1 = " << maxima_.gross[0] << " minimum1 = " << maxima_.klein[0] << endl;
119 cout << "Mmax1 = " << maxima_.Mmax[0] << endl;
120 cout << "maximum2 = " << maxima_.gross[1] << " minimum2 = " << maxima_.klein[1] << endl;
121 cout << "Mmax2 = " << maxima_.Mmax[1] << endl;
122 cout << "-------------------------------------------------------------" << endl;
123
124 int i;
125 // the deletion of commands is really uggly!
126
127 if ( nphokharadecays == 0 )
128 {
129 delete[] commands;
130 commands = 0;
131 return;
132 }
133
134 for ( i = 0; i < nphokharadecays; i++ )
135 {
136 if ( phokharadecays[i] == this )
137 {
138 phokharadecays[i] = phokharadecays[nphokharadecays - 1];
139 nphokharadecays--;
140 if ( nphokharadecays == 0 )
141 {
142 delete[] commands;
143 commands = 0;
144 }
145 return;
146 }
147 }
148
149 report( ERROR, "EvtGen" ) << "Error in destroying Phokhara model!" << endl;
150}
151
152void EvtPhokhara::getName( std::string& model_name ) { model_name = "PHOKHARA"; }
153
154EvtDecayBase* EvtPhokhara::clone() { return new EvtPhokhara; }
155
156void EvtPhokhara::initProbMax() { noProbMax(); }
157
158void EvtPhokhara::init_mode( EvtParticle* p ) {
159 m_pion = getArg( 0 );
160 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
161 // 2pi+2pi-(3),ppbar(4),nnbar(5),
162 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
163 // Lamb Lambbar->pi-pi+ppbar(9)
164 if ( m_pion < 0 || m_pion > 9 )
165 {
166 std::cout << "mode index for phokhar 0~9, but you give " << m_pion << std::endl;
167 abort();
168 }
169 EvtId myvpho = EvtPDL::getId( "vpho" );
170 m_E = p->mass(); // EvtPDL::getMeanMass(myvpho);
171 ///======== list parameters to be initialized
172 m_tagged = 0;
173 m_nm = 50000; // # of events to determine the maximum
174 m_nlo = 1; // Born(0), NLO(1)
175 m_w = 0.0001; // soft photon cutoff
176 m_fsr = 1; // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
177 m_fsrnlo = 1; // yes(1), no(0)
178 m_NarrowRes = 1; // none(0) jpsi(1) psip(2)
179 m_FF_Kaon = 1; // constrained (0), unconstrained (1), Kuhn-Khodjamirian-Bruch (2)
180 m_ivac = 0; // yes(1), no(0)
181 m_FF_Pion = 0; // KS Pionformfactor(0), GS Pionformfactor old(1) and new(2)
182 m_f0_model = 0; // f0+f0(600): KK model(0), no structure(1), no f0+f0(600)(2), f0 KLOE(3)
183 m_q2min = 0.0; // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
184 m_q2_min_c = 0.0447; // minimal inv. mass squared of the hadrons(muons)(GeV^2)
185 m_q2_max_c = m_E * m_E; // maximal inv. mass squared of the hadrons(muons)(GeV^2)
186 m_gmin = 0.001; // minimal photon energy/missing energy (GeV)
187 m_phot1cut = 0.0; // maximal photon angle/missing momentum angle (grad)
188 m_phot2cut = 180.0; // maximal photon angle/missing momentum angle (grad)
189 m_pi1cut = 0.0; // minimal hadrons(muons) angle (grad)
190 m_pi2cut = 180.0; // maximal hadrons(muons) angle (grad)
191
192 if ( !( m_pion == 0 || m_pion == 1 || m_pion == 6 ) )
193 {
194 m_fsr = 0;
195 m_fsrnlo = 0;
196 }
197 if ( m_pion == 9 ) { m_nlo = 0; }
198 // --- input parameter initialization -----------
199 m_initSeed = 123456789;
200 RLXDINIT( 1, m_initSeed );
201 maxima_.iprint = 0;
202 flags_.nlo = m_nlo;
203 flags_.pion = m_pion;
204 flags_.fsr = m_fsr;
205 flags_.fsrnlo = m_fsrnlo;
206 flags_.ivac = m_ivac;
207 flags_.FF_pion = m_FF_Pion;
208 flags_.f0_model = m_f0_model;
209 flags_.FF_kaon = m_FF_Kaon;
210 flags_.narr_res = m_NarrowRes;
211
212 ctes_.Sp = m_E * m_E;
213 ;
214
215 cuts_.w = m_w;
216 cuts_.q2min = m_q2min;
217 cuts_.q2_min_c = m_q2_min_c;
218 cuts_.q2_max_c = m_q2_max_c;
219 cuts_.gmin = m_gmin;
220 cuts_.phot1cut = m_phot1cut;
221 cuts_.phot2cut = m_phot2cut;
222 cuts_.pi1cut = m_pi1cut;
223 cuts_.pi2cut = m_pi2cut;
224
225 INPUT();
226
227 // --- print run data ------------------------------
228 cout << "-------------------------------------------------------------" << endl;
229 if ( flags_.pion == 0 ) cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
230 else if ( flags_.pion == 1 ) cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
231 else if ( flags_.pion == 2 )
232 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
233 else if ( flags_.pion == 3 ) cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
234 else if ( flags_.pion == 4 ) cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
235 else if ( flags_.pion == 5 ) cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
236 else if ( flags_.pion == 6 ) cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
237 else if ( flags_.pion == 7 ) cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
238 else if ( flags_.pion == 8 )
239 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
240 else if ( flags_.pion == 9 )
241 {
242 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
243 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
244 }
245 else cout << " PHOKHARA 7.0: not yet implemented" << endl;
246
247 // --------------------------------
248 cout << "--------------------------------------------------------------" << endl;
249 printf( " %s %f %s\n", "cms total energy = ", sqrt( ctes_.Sp ),
250 " GeV " );
251
252 // if((cuts_.gmin/2.0/ctes_.ebeam) < 0.000098){
253 if ( cuts_.gmin < 0.001 )
254 {
255 cerr << " minimal missing energy set too small" << endl;
256 abort();
257 }
258 printf( " %s %f %s\n", "minimal tagged photon energy = ", cuts_.gmin, " GeV " );
259 printf( " %s %f,%f\n", "angular cuts on tagged photon = ", cuts_.phot1cut,
260 cuts_.phot2cut );
261
262 // --------------------------------
263 if ( flags_.pion == 0 )
264 printf( " %s %f,%f\n", "angular cuts on muons = ", cuts_.pi1cut,
265 cuts_.pi2cut );
266 else if ( flags_.pion == 4 )
267 printf( " %s %f,%f\n", "angular cuts on protons = ", cuts_.pi1cut,
268 cuts_.pi2cut );
269 else if ( flags_.pion == 5 )
270 printf( " %s %f,%f\n", "angular cuts on neutrons = ", cuts_.pi1cut,
271 cuts_.pi2cut );
272 else if ( ( flags_.pion == 6 ) || ( flags_.pion == 7 ) )
273 printf( " %s %f,%f\n", "angular cuts on kaons = ", cuts_.pi1cut,
274 cuts_.pi2cut );
275 else if ( flags_.pion == 9 )
276 printf( " %s %f,%f\n", "angular cuts on pions and protons = ", cuts_.pi1cut,
277 cuts_.pi2cut );
278 else
279 printf( " %s %f,%f\n", "angular cuts on pions = ", cuts_.pi1cut,
280 cuts_.pi2cut );
281 if ( flags_.pion == 0 )
282 printf( " %s %f %s\n", "min. muons-tagged photon inv.mass^2 = ", cuts_.q2min,
283 " GeV^2" );
284 else if ( flags_.pion == 4 )
285 printf( " %s %f %s\n", "min. protons-tagged photon inv.mass^2 = ", cuts_.q2min,
286 " GeV^2" );
287 else if ( flags_.pion == 5 )
288 printf( " %s %f %s\n", "min. neutrons-tagged photon inv.mass^2 = ", cuts_.q2min,
289 " GeV^2" );
290 else if ( ( flags_.pion == 6 ) || ( flags_.pion == 7 ) )
291 printf( " %s %f %s\n", "min. kaons-tagged photon inv.mass^2 = ", cuts_.q2min,
292 " GeV^2" );
293 else if ( flags_.pion == 9 )
294 printf( " %s %f %s\n", " min. lambdas-tagged photon inv.mass^2 = ", cuts_.q2min,
295 " GeV^2" );
296 else
297 printf( " %s %f %s\n", "min. pions-tagged photon inv.mass^2 = ", cuts_.q2min,
298 " GeV^2" );
299 cos1min = cos( cuts_.phot2cut * ctes_.pi / 180.0 ); // photon1 angle cuts in the
300 cos1max = cos( cuts_.phot1cut * ctes_.pi / 180.0 ); // LAB rest frame
301 cos2min = -1.0; // photon2 angle limits
302 cos2max = 1.0; //
303 cos3min = -1.0; // hadrons/muons angle limits
304 cos3max = 1.0; // in their rest frame
305 if ( flags_.pion == 0 ) // virtual photon energy cut
306 qqmin = 4.0 * ctes_.mmu * ctes_.mmu;
307 else if ( flags_.pion == 1 ) qqmin = 4.0 * ctes_.mpi * ctes_.mpi;
308 else if ( flags_.pion == 2 )
309 qqmin = 4.0 * ( ctes_.mpi + ctes_.mpi0 ) * ( ctes_.mpi + ctes_.mpi0 );
310 else if ( flags_.pion == 3 ) qqmin = 16.0 * ctes_.mpi * ctes_.mpi;
311 else if ( flags_.pion == 4 ) qqmin = 4.0 * ctes_.mp * ctes_.mp;
312 else if ( flags_.pion == 5 ) qqmin = 4.0 * ctes_.mnt * ctes_.mnt;
313 else if ( flags_.pion == 6 ) qqmin = 4.0 * ctes_.mKp * ctes_.mKp;
314 else if ( flags_.pion == 7 ) qqmin = 4.0 * ctes_.mKn * ctes_.mKn;
315 else if ( flags_.pion == 8 )
316 qqmin = ( 2.0 * ctes_.mpi + ctes_.mpi0 ) * ( 2.0 * ctes_.mpi + ctes_.mpi0 );
317 else if ( flags_.pion == 9 ) qqmin = 4.0 * ctes_.mlamb * ctes_.mlamb;
318 qqmax = ctes_.Sp - 2.0 * sqrt( ctes_.Sp ) * cuts_.gmin; // if only one photon
319 if ( cuts_.q2_max_c < qqmax ) qqmax = cuts_.q2_max_c; // external cuts
320
321 // -------------------
322 if ( ( cuts_.q2_min_c > qqmin ) &&
323 ( cuts_.q2_min_c <
324 ( ctes_.Sp * ( 1.0 - 2.0 * ( cuts_.gmin / sqrt( ctes_.Sp ) + cuts_.w ) ) ) ) )
325 qqmin = cuts_.q2_min_c;
326 else
327 {
328 cerr << "------------------------------" << endl;
329 cerr << " Q^2_min TOO SMALL" << endl;
330 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
331 cerr << "------------------------------" << endl;
332 }
333 // -------------------
334 if ( qqmax <= qqmin )
335 {
336 cerr << " Q^2_max to small " << endl;
337 cerr << " Q^2_max = " << qqmax << endl;
338 cerr << " Q^2_min = " << qqmin << endl;
339 abort();
340 }
341
342 // -------------------
343 if ( flags_.pion == 0 )
344 {
345 printf( " %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin, " GeV^2" );
346 printf( " %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax, " GeV^2" );
347 }
348 else if ( flags_.pion == 1 )
349 {
350 printf( " %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin, " GeV^2" );
351 printf( " %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax, " GeV^2" );
352 }
353 else if ( flags_.pion == 4 )
354 {
355 printf( " %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin, " GeV^2" );
356 printf( " %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax, " GeV^2" );
357 }
358 else if ( flags_.pion == 5 )
359 {
360 printf( " %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin, " GeV^2" );
361 printf( " %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax, " GeV^2" );
362 }
363 else if ( ( flags_.pion == 6 ) || ( flags_.pion == 7 ) )
364 {
365 printf( " %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin, " GeV^2" );
366 printf( " %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax, " GeV^2" );
367 }
368 else if ( flags_.pion == 8 )
369 {
370 printf( " %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin, " GeV^2" );
371 printf( " %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax, " GeV^2" );
372 }
373 else if ( flags_.pion == 9 )
374 {
375 printf( " %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin, " GeV^2" );
376 printf( " %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax, " GeV^2" );
377 }
378 else
379 {
380 printf( " %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin, " GeV^2" );
381 printf( " %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax, " GeV^2" );
382 }
383 // -------------------
384 if ( flags_.nlo == 0 )
385 {
386 cout << "Born" << endl;
387 if ( flags_.fsrnlo != 0 )
388 {
389 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
390 abort();
391 }
392 }
393 // -------------------
394 if ( ( flags_.pion == 9 ) && ( flags_.nlo != 0 ) )
395 {
396 cerr << "WRONG NLO flag - only Born allowed for Lambdas" << endl;
397 cerr << "If you feel that you need NLO" << endl;
398 cerr << "please contact the authors" << endl;
399 abort();
400 }
401 // -------------------
402 if ( flags_.nlo == 1 )
403 printf( " %s %f\n", "NLO: soft photon cutoff w = ", cuts_.w );
404 if ( ( flags_.pion <= 1 ) || ( flags_.pion == 6 ) )
405 {
406
407 if ( !( ( flags_.fsr == 1 ) || ( flags_.fsr == 2 ) || ( flags_.fsrnlo == 0 ) ||
408 ( flags_.fsr == 1 ) || ( flags_.fsrnlo == 1 ) || ( flags_.fsr == 0 ) ||
409 ( flags_.fsrnlo == 0 ) ) )
410 {
411 cerr << "WRONG combination of FSR, FSRNLO flags" << endl;
412 abort();
413 }
414 // ------------------
415 if ( flags_.fsr == 0 ) cout << "ISR only" << endl;
416 else if ( flags_.fsr == 1 ) cout << "ISR+FSR" << endl;
417 else if ( flags_.fsr == 2 )
418 {
419 if ( flags_.nlo == 0 ) cout << "ISR+INT+FSR" << endl;
420 else
421 {
422 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
423 abort();
424 }
425 }
426 else
427 {
428 cerr << "WRONG FSR flag" << flags_.fsr << endl;
429 abort();
430 }
431 if ( flags_.fsrnlo == 1 ) cout << "IFSNLO included" << endl;
432 }
433 else
434 {
435 if ( ( flags_.fsr == 0 ) && ( flags_.fsrnlo == 0 ) ) cout << "ISR only" << endl;
436 else
437 {
438 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
439 abort();
440 }
441 }
442 // ------------------
443 if ( flags_.ivac == 0 ) { cout << "Vacuum polarization is NOT included" << endl; }
444 else if ( flags_.ivac == 1 )
445 {
446 cout << "Vacuum polarization by Fred Jegerlehner "
447 "(http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)"
448 << endl;
449 }
450 else if ( flags_.ivac == 2 )
451 {
452 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner"
453 << endl;
454 }
455 else
456 {
457 cout << "WRONG vacuum polarization switch" << endl;
458 abort();
459 }
460
461 // -----------------
462 if ( flags_.pion == 1 )
463 {
464 if ( flags_.FF_pion == 0 ) cout << "Kuhn-Santamaria PionFormFactor" << endl;
465 else if ( flags_.FF_pion == 1 ) cout << "Gounaris-Sakurai PionFormFactor old" << endl;
466 else if ( flags_.FF_pion == 2 ) cout << "Gounaris-Sakurai PionFormFactor new" << endl;
467 else
468 {
469 cout << "WRONG PionFormFactor switch" << endl;
470 abort();
471 }
472 if ( flags_.fsr != 0 )
473 {
474 if ( flags_.f0_model == 0 ) cout << "f0+f0(600): K+K- model" << endl;
475 else if ( flags_.f0_model == 1 ) cout << "f0+f0(600): \"no structure\" model" << endl;
476 else if ( flags_.f0_model == 2 ) cout << "NO f0+f0(600)" << endl;
477 else if ( flags_.f0_model == 3 )
478 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
479 else
480 {
481 cout << "WRONG f0+f0(600) switch" << endl;
482 abort();
483 }
484 }
485 }
486 //-------
487 if ( ( flags_.pion == 6 ) || ( flags_.pion == 7 ) )
488 {
489 if ( flags_.FF_kaon == 0 ) cout << "constrained KaonFormFactor" << endl;
490 else if ( flags_.FF_kaon == 1 ) { cout << "unconstrained KaonFormFactor" << endl; }
491 else if ( flags_.FF_kaon == 2 )
492 { cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl; }
493 else
494 {
495 cout << "WRONG KaonFormFactor switch" << endl;
496 abort();
497 }
498 }
499 // --------------------
500 if ( ( flags_.pion == 0 ) || ( flags_.pion == 1 ) || ( flags_.pion == 6 ) ||
501 ( flags_.pion == 7 ) )
502 {
503 if ( flags_.narr_res == 1 ) { cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl; }
504 else if ( flags_.narr_res == 2 )
505 { cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl; }
506 else if ( flags_.narr_res != 0 )
507 {
508 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
509 abort();
510 }
511 }
512 // ------
513 cout << "--------------------------------------------------------------" << endl;
514 //
515 // =================================================
516 // --- finding the maximum -------------------------
517 for ( int i = 0; i < 2; i++ )
518 {
519 maxima_.Mmax[i] = 1.0;
520 maxima_.gross[i] = 0.0;
521 maxima_.klein[i] = 0.0;
522 }
523
524 if ( flags_.nlo == 0 ) maxima_.Mmax[1] = 0.0; // only 1 photon events generated
525
526 maxima_.tr[0] = 0.0;
527 maxima_.tr[1] = 0.0;
528 maxima_.count[0] = 0.0;
529 maxima_.count[1] = 0.0;
530
531 // =================================================
532 // --- beginning the MC loop event generation ------
533 for ( int j = 1; j <= m_nm; j++ )
534 {
535 RANLXDF( Ar_r, 1 );
536 Ar[1] = Ar_r[0];
537 if ( Ar[1] <= ( maxima_.Mmax[0] / ( maxima_.Mmax[0] + maxima_.Mmax[1] ) ) )
538 {
539 maxima_.count[0] = maxima_.count[0] + 1.0;
540 GEN_1PH( 1, qqmin, qqmax, cos1min, cos1max, cos3min, cos3max );
541 }
542 else
543 {
544 maxima_.count[1] = maxima_.count[1] + 1.0;
545 GEN_2PH( 1, qqmin, cos1min, cos1max, cos2min, cos2max, cos3min, cos3max );
546 }
547 }
548 // --- end of the MC loop --------------------------
549 // =================================================
550 // --- for the second run ---
551 maxima_.Mmax[0] = maxima_.gross[0] + .05 * sqrt( maxima_.gross[0] * maxima_.gross[0] );
552 maxima_.Mmax[1] = maxima_.gross[1] +
553 ( .03 + .02 * ctes_.Sp ) * sqrt( maxima_.gross[1] * maxima_.gross[1] );
554 theMmax[m_pion][0] = maxima_.Mmax[0];
555 theMmax[m_pion][1] = maxima_.Mmax[1];
556
557 if ( ( flags_.pion == 1 ) && ( flags_.fsrnlo == 1 ) )
558 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.5;
559 if ( ( flags_.pion == 0 ) && ( flags_.fsrnlo == 1 ) )
560 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.5;
561
562 if ( ( flags_.pion == 0 ) && ( flags_.fsr == 1 ) && ( flags_.fsrnlo == 0 ) )
563 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.2;
564
565 if ( ( flags_.pion == 2 ) || ( flags_.pion == 3 ) )
566 {
567 maxima_.Mmax[0] = maxima_.Mmax[0] * 1.1;
568 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.1;
569 }
570
571 if ( flags_.pion == 8 )
572 {
573 maxima_.Mmax[0] = maxima_.Mmax[0] * 1.08;
574 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.1;
575 }
576 // --- end of the second run -----------------------
577
578 maxima_.tr[0] = 0.0;
579 maxima_.tr[1] = 0.0;
580 maxima_.count[0] = 0.0;
581 maxima_.count[1] = 0.0;
582}
583
584void EvtPhokhara::init() {
585 checkNDaug( 0 );
586 checkNArg( 1 );
587 nevtgen.resize( 10 );
588 theMmax.resize( 10 );
589 for ( int i = 0; i <= 10; i++ )
590 {
591 theMmax[i].resize( 2 );
592 nevtgen[i] = 0;
593 }
594
595 Vfs.push_back( " mu+mu-" );
596 Vfs.push_back( " pi+pi-" );
597 Vfs.push_back( " 2pi0pi+pi-" );
598 Vfs.push_back( " 2pi+2pi-" );
599 Vfs.push_back( " ppbar" );
600 Vfs.push_back( " nnbar" );
601 Vfs.push_back( " K+K-" );
602 Vfs.push_back( " K0K0bar" );
603 Vfs.push_back( " pi+pi-pi0" );
604 Vfs.push_back( " Lambda Lambdabar" );
605
606 std::string locvp = getenv( "BESEVTGENROOT" );
607 system( "cat $BESEVTGENROOT/share/phokhara.param>phokhara.param" );
608
609 if ( getParentId().isAlias() )
610 {
611
612 report( ERROR, "EvtGen" ) << "EvtPhokhara finds that you are decaying the" << endl
613 << " aliased particle " << EvtPDL::name( getParentId() ).c_str()
614 << " with the Phokhara model" << endl
615 << " this does not work, please modify decay table." << endl;
616 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
617 ::abort();
618 }
619
620 // store(this);
621}
622
623std::string EvtPhokhara::commandName() { return std::string( "PhokharaPar" ); }
624
625void EvtPhokhara::command( std::string cmd ) {
626
627 if ( ncommand == lcommand )
628 {
629
630 lcommand = 10 + 2 * lcommand;
631
632 std::string* newcommands = new std::string[lcommand];
633
634 int i;
635
636 for ( i = 0; i < ncommand; i++ ) { newcommands[i] = commands[i]; }
637
638 delete[] commands;
639
640 commands = newcommands;
641 }
642
643 commands[ncommand] = cmd;
644
645 ncommand++;
646}
647
648void EvtPhokhara::decay( EvtParticle* p ) {
649 EvtId myvpho = EvtPDL::getId( "vpho" );
650 if ( p->getId() != myvpho )
651 {
652 std::cout << "Parent particle is required to be vpho for Phokhara model" << std::endl;
653 abort();
654 }
655 m_pion = getArg( 0 );
656 if ( nevtgen[m_pion] == 0 ) { init_mode( p ); }
657 else { init_evt( p ); }
658 std::cout << "PHOKHARA : " << Vfs[m_pion] << " mode " << std::endl;
659
660 int istdheppar = EvtPDL::getStdHep( p->getId() );
661 int ntrials = 0;
662 int tr_old[2];
663 tr_old[0] = (int)maxima_.tr[0];
664 tr_old[1] = (int)maxima_.tr[1];
665
666 while ( ntrials < 10000 )
667 {
668 ievent++;
669 RANLXDF( Ar_r, 1 );
670 Ar[1] = Ar_r[0];
671
672 if ( Ar[1] <= ( maxima_.Mmax[0] / ( maxima_.Mmax[0] + maxima_.Mmax[1] ) ) )
673 {
674 maxima_.count[0] = maxima_.count[0] + 1.0;
675 GEN_1PH( 2, qqmin, qqmax, cos1min, cos1max, cos3min, cos3max );
676 }
677 else
678 {
679 maxima_.count[1] = maxima_.count[1] + 1.0;
680 GEN_2PH( 2, qqmin, cos1min, cos1max, cos2min, cos2max, cos3min, cos3max );
681 }
682
683 if ( ( (int)maxima_.tr[0] + (int)maxima_.tr[1] ) >
684 ( tr_old[0] + tr_old[1] ) ) // event accepted after cuts
685 { goto storedEvents; }
686 ntrials++;
687 }
688
689 std::cout << "FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate."
690 << std::endl;
691 //----
692storedEvents:
693 int more = 0;
694 int numstable = 0;
695 int numparton = 0;
696 EvtId evtnumstable[100]; //
697 EvtVector4R p4[20];
698
699 // except ISR photos, products depending on channel
700 if ( flags_.pion == 0 )
701 { // mu+ mu-
702 // mu+
703 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -13 );
704 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
705 ctes_.momenta[3][5] );
706 numstable++;
707 // mu -
708 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 13 );
709 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
710 ctes_.momenta[3][6] );
711 numstable++;
712 }
713 if ( flags_.pion == 1 )
714 { // pi+ pi-
715 // pi+
716 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
717 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
718 ctes_.momenta[3][5] );
719 numstable++;
720 // pi -
721 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
722 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
723 ctes_.momenta[3][6] );
724 numstable++;
725 }
726 if ( flags_.pion == 2 )
727 { // pi+ pi-2pi0
728 // pi0
729 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
730 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
731 ctes_.momenta[3][5] );
732 numstable++;
733 // pi0
734 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
735 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
736 ctes_.momenta[3][6] );
737 numstable++;
738 // pi-
739 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
740 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
741 ctes_.momenta[3][7] );
742 numstable++;
743 // pi +
744 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
745 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
746 ctes_.momenta[3][8] );
747 numstable++;
748 }
749 if ( flags_.pion == 3 )
750 { // 2(pi+ pi-)
751 // pi+
752 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
753 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
754 ctes_.momenta[3][5] );
755 numstable++;
756 // pi-
757 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
758 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
759 ctes_.momenta[3][6] );
760 numstable++;
761 // pi+
762 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
763 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
764 ctes_.momenta[3][7] );
765 numstable++;
766 // pi -
767 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
768 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
769 ctes_.momenta[3][8] );
770 numstable++;
771 }
772 if ( flags_.pion == 4 )
773 { // ppbar
774 // pbar
775 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2212 );
776 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
777 ctes_.momenta[3][5] );
778 numstable++;
779 // p
780 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 2212 );
781 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
782 ctes_.momenta[3][6] );
783 numstable++;
784 }
785 if ( flags_.pion == 5 )
786 { // nnbar
787 // pbar
788 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2112 );
789 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
790 ctes_.momenta[3][5] );
791 numstable++;
792 // p
793 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 2112 );
794 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
795 ctes_.momenta[3][6] );
796 numstable++;
797 }
798 if ( flags_.pion == 6 )
799 { // K+ K-
800 // K+
801 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 321 );
802 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
803 ctes_.momenta[3][5] );
804 numstable++;
805 // K -
806 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -321 );
807 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
808 ctes_.momenta[3][6] );
809 numstable++;
810 }
811 if ( flags_.pion == 7 )
812 { // K0K0bar
813 // Kbar
814 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 310 );
815 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
816 ctes_.momenta[3][5] );
817 numstable++;
818 // K0
819 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 130 );
820 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
821 ctes_.momenta[3][6] );
822 numstable++;
823 }
824 if ( flags_.pion == 8 )
825 { // pi+ pi-pi0
826 // pi+
827 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
828 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
829 ctes_.momenta[3][5] );
830 numstable++;
831 // pi-
832 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
833 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
834 ctes_.momenta[3][6] );
835 numstable++;
836 // pi0
837 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
838 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
839 ctes_.momenta[3][7] );
840 numstable++;
841 }
842 if ( flags_.pion == 9 )
843 { // Lambda Lambdabar-> pi+ pi- ppbar
844 // pi+
845 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
846 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
847 ctes_.momenta[3][7] );
848 numstable++;
849 // pbar
850 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2212 );
851 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
852 ctes_.momenta[3][8] );
853 numstable++;
854 // pi-
855 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
856 p4[numstable].set( ctes_.momenta[0][9], ctes_.momenta[1][9], ctes_.momenta[2][9],
857 ctes_.momenta[3][9] );
858 numstable++;
859 // p
860 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 2212 );
861 p4[numstable].set( ctes_.momenta[0][10], ctes_.momenta[1][10], ctes_.momenta[2][10],
862 ctes_.momenta[3][10] );
863 numstable++;
864 }
865
866 // ISR gamma
867 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 22 );
868 p4[numstable].set( ctes_.momenta[0][2], ctes_.momenta[1][2], ctes_.momenta[2][2],
869 ctes_.momenta[3][2] );
870 numstable++;
871 if ( ctes_.momenta[0][3] != 0 ) // second photon exists
872 {
873 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 22 );
874 p4[numstable].set( ctes_.momenta[0][3], ctes_.momenta[1][3], ctes_.momenta[2][3],
875 ctes_.momenta[3][3] );
876 numstable++;
877 }
878
879 int channel = EvtDecayTable::inChannelList( p->getId(), numstable, evtnumstable );
880 more = ( channel != -1 );
881 if ( more )
882 {
883 std::cout << "Existence of mode " << channel
884 << " in exclusive decay list has the same final state as this one" << std::endl;
885 abort();
886 }
887
888 p->makeDaughters( numstable, evtnumstable );
889
890 int ndaugFound = 0;
891 EvtVector4R SUMP4( 0, 0, 0, 0 );
892 for ( int i = 0; i < numstable; i++ )
893 {
894 p->getDaug( i )->init( evtnumstable[i], p4[i] );
895 ndaugFound++;
896 }
897 if ( ndaugFound == 0 )
898 {
899 report( ERROR, "EvtGen" ) << "Phokhara has failed to do a decay ";
900 report( ERROR, "EvtGen" ) << EvtPDL::name( p->getId() ).c_str() << " " << p->mass()
901 << endl;
902 assert( 0 );
903 }
904
905 nevtgen[m_pion]++;
906 return;
907}
908
909void EvtPhokhara::store( EvtDecayBase* jsdecay ) {
910
911 if ( nphokharadecays == ntable )
912 {
913
914 EvtDecayBasePtr* newphokharadecays = new EvtDecayBasePtr[2 * ntable + 10];
915 int i;
916 for ( i = 0; i < ntable; i++ ) { newphokharadecays[i] = phokharadecays[i]; }
917 ntable = 2 * ntable + 10;
918 delete[] phokharadecays;
919 phokharadecays = newphokharadecays;
920 }
921
922 phokharadecays[nphokharadecays++] = jsdecay;
923}
924
925void EvtPhokhara::PhokharaInit( int dummy ) {
926 static int first = 1;
927 if ( first )
928 {
929
930 first = 0;
931 // for(int i=0;i<ncommand;i++)
932 // lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
933 }
934}
935
936void EvtPhokhara::init_evt( EvtParticle* p ) {
937 m_pion = getArg( 0 );
938 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
939 // 2pi+2pi-(3),ppbar(4),nnbar(5),
940 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
941 // Lamb Lambbar->pi-pi+ppbar(9)
942 if ( m_pion < 0 || m_pion > 9 )
943 {
944 std::cout << "mode index for phokhar 0~9, but you give " << m_pion << std::endl;
945 abort();
946 }
947 EvtId myvpho = EvtPDL::getId( "vpho" );
948 m_E = p->mass(); // EvtPDL::getMeanMass(myvpho);
949 ///======== list parameters to be initialized
950 m_tagged = 0;
951 m_nm = 50000; // # of events to determine the maximum
952 m_nlo = 1; // Born(0), NLO(1)
953 m_w = 0.0001; // soft photon cutoff
954 m_fsr = 1; // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
955 m_fsrnlo = 1; // yes(1), no(0)
956 m_NarrowRes = 1; // none(0) jpsi(1) psip(2)
957 m_FF_Kaon = 1; // constrained (0), unconstrained (1), Kuhn-Khodjamirian-Bruch (2)
958 m_ivac = 0; // yes(1), no(0)
959 m_FF_Pion = 0; // KS Pionformfactor(0), GS Pionformfactor old(1) and new(2)
960 m_f0_model = 0; // f0+f0(600): KK model(0), no structure(1), no f0+f0(600)(2), f0 KLOE(3)
961 m_q2min = 0.0; // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
962 m_q2_min_c = 0.0447; // minimal inv. mass squared of the hadrons(muons)(GeV^2)
963 m_q2_max_c = m_E * m_E; // maximal inv. mass squared of the hadrons(muons)(GeV^2)
964 m_gmin = 0.001; // minimal photon energy/missing energy (GeV)
965 m_phot1cut = 0.0; // maximal photon angle/missing momentum angle (grad)
966 m_phot2cut = 180.0; // maximal photon angle/missing momentum angle (grad)
967 m_pi1cut = 0.0; // minimal hadrons(muons) angle (grad)
968 m_pi2cut = 180.0; // maximal hadrons(muons) angle (grad)
969
970 if ( !( m_pion == 0 || m_pion == 1 || m_pion == 6 ) )
971 {
972 m_fsr = 0;
973 m_fsrnlo = 0;
974 }
975 if ( m_pion == 9 ) { m_nlo = 0; }
976 // --- input parameter initialization -----------
977 maxima_.iprint = 0;
978 flags_.nlo = m_nlo;
979 flags_.pion = m_pion;
980 flags_.fsr = m_fsr;
981 flags_.fsrnlo = m_fsrnlo;
982 flags_.ivac = m_ivac;
983 flags_.FF_pion = m_FF_Pion;
984 flags_.f0_model = m_f0_model;
985 flags_.FF_kaon = m_FF_Kaon;
986 flags_.narr_res = m_NarrowRes;
987
988 ctes_.Sp = m_E * m_E;
989 ;
990
991 cuts_.w = m_w;
992 cuts_.q2min = m_q2min;
993 cuts_.q2_min_c = m_q2_min_c;
994 cuts_.q2_max_c = m_q2_max_c;
995 cuts_.gmin = m_gmin;
996 cuts_.phot1cut = m_phot1cut;
997 cuts_.phot2cut = m_phot2cut;
998 cuts_.pi1cut = m_pi1cut;
999 cuts_.pi2cut = m_pi2cut;
1000
1001 INPUT();
1002
1003 cos1min = cos( cuts_.phot2cut * ctes_.pi / 180.0 ); // photon1 angle cuts in the
1004 cos1max = cos( cuts_.phot1cut * ctes_.pi / 180.0 ); // LAB rest frame
1005 cos2min = -1.0; // photon2 angle limits
1006 cos2max = 1.0; //
1007 cos3min = -1.0; // hadrons/muons angle limits
1008 cos3max = 1.0; // in their rest frame
1009 if ( flags_.pion == 0 ) // virtual photon energy cut
1010 qqmin = 4.0 * ctes_.mmu * ctes_.mmu;
1011 else if ( flags_.pion == 1 ) qqmin = 4.0 * ctes_.mpi * ctes_.mpi;
1012 else if ( flags_.pion == 2 )
1013 qqmin = 4.0 * ( ctes_.mpi + ctes_.mpi0 ) * ( ctes_.mpi + ctes_.mpi0 );
1014 else if ( flags_.pion == 3 ) qqmin = 16.0 * ctes_.mpi * ctes_.mpi;
1015 else if ( flags_.pion == 4 ) qqmin = 4.0 * ctes_.mp * ctes_.mp;
1016 else if ( flags_.pion == 5 ) qqmin = 4.0 * ctes_.mnt * ctes_.mnt;
1017 else if ( flags_.pion == 6 ) qqmin = 4.0 * ctes_.mKp * ctes_.mKp;
1018 else if ( flags_.pion == 7 ) qqmin = 4.0 * ctes_.mKn * ctes_.mKn;
1019 else if ( flags_.pion == 8 )
1020 qqmin = ( 2.0 * ctes_.mpi + ctes_.mpi0 ) * ( 2.0 * ctes_.mpi + ctes_.mpi0 );
1021 else if ( flags_.pion == 9 ) qqmin = 4.0 * ctes_.mlamb * ctes_.mlamb;
1022 qqmax = ctes_.Sp - 2.0 * sqrt( ctes_.Sp ) * cuts_.gmin; // if only one photon
1023 if ( cuts_.q2_max_c < qqmax ) qqmax = cuts_.q2_max_c; // external cuts
1024
1025 // -------------------
1026 if ( ( cuts_.q2_min_c > qqmin ) &&
1027 ( cuts_.q2_min_c <
1028 ( ctes_.Sp * ( 1.0 - 2.0 * ( cuts_.gmin / sqrt( ctes_.Sp ) + cuts_.w ) ) ) ) )
1029 qqmin = cuts_.q2_min_c;
1030 else {}
1031
1032 // =================================================
1033 // --- finding the maximum -------------------------
1034 for ( int i = 0; i < 2; i++ )
1035 {
1036 maxima_.Mmax[i] = 1.0;
1037 maxima_.gross[i] = 0.0;
1038 maxima_.klein[i] = 0.0;
1039 }
1040
1041 if ( flags_.nlo == 0 ) maxima_.Mmax[1] = 0.0; // only 1 photon events generated
1042
1043 maxima_.tr[0] = 0.0;
1044 maxima_.tr[1] = 0.0;
1045 maxima_.count[0] = 0.0;
1046 maxima_.count[1] = 0.0;
1047
1048 // =================================================
1049 // --- for the second run ---
1050 maxima_.Mmax[0] = theMmax[m_pion][0];
1051 maxima_.Mmax[1] = theMmax[m_pion][1];
1052 if ( ( flags_.pion == 1 ) && ( flags_.fsrnlo == 1 ) )
1053 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.5;
1054 if ( ( flags_.pion == 0 ) && ( flags_.fsrnlo == 1 ) )
1055 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.5;
1056
1057 if ( ( flags_.pion == 0 ) && ( flags_.fsr == 1 ) && ( flags_.fsrnlo == 0 ) )
1058 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.2;
1059
1060 if ( ( flags_.pion == 2 ) || ( flags_.pion == 3 ) )
1061 {
1062 maxima_.Mmax[0] = maxima_.Mmax[0] * 1.1;
1063 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.1;
1064 }
1065
1066 if ( flags_.pion == 8 )
1067 {
1068 maxima_.Mmax[0] = maxima_.Mmax[0] * 1.08;
1069 maxima_.Mmax[1] = maxima_.Mmax[1] * 1.1;
1070 }
1071 // --- end of the second run -----------------------
1072
1073 maxima_.tr[0] = 0.0;
1074 maxima_.tr[1] = 0.0;
1075 maxima_.count[0] = 0.0;
1076 maxima_.count[1] = 0.0;
1077}
INPUT
#define INPUT
Definition BesBdkRc.cxx:31
EvtDecayBasePtr
EvtDecayBase * EvtDecayBasePtr
Definition EvtDecayTable.hh:32
ctes_
struct @053254170326070136226344307237142165176240334330 ctes_
RLXDINIT
#define RLXDINIT(LUXURY, SEED)
Definition EvtPhokharaDef.hh:88
maxima_
struct @027003056066344010031101102046265032310161340072 maxima_
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
lambda_par_
struct @152360376016154317326265314214112361240371075021 lambda_par_
GEN_2PH
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
flags_
struct @366025114004024134344043225357106161032107165361 flags_
cuts_
struct @276363222274272223263152166363125067340201005217 cuts_
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
ERROR
@ ERROR
Definition EvtReport.hh:49
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
EvtDecayBase::getArg
double getArg(int j)
Definition EvtDecayBase.cc:520
EvtDecayBase::getParentId
EvtId getParentId()
Definition EvtDecayBase.hh:63
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayTable::inChannelList
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition EvtDecayTable.cc:831
EvtId
Definition EvtId.hh:27
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:61
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition EvtPDL.cc:232
EvtPDL::name
static std::string name(EvtId i)
Definition EvtPDL.hh:70
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition EvtParticle.cc:1150
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::mass
double mass() const
Definition EvtParticle.cc:114
EvtPhokhara::clone
EvtDecayBase * clone()
Definition EvtPhokhara.cc:154
EvtPhokhara::command
void command(std::string cmd)
Definition EvtPhokhara.cc:625
EvtPhokhara::~EvtPhokhara
virtual ~EvtPhokhara()
Definition EvtPhokhara.cc:66
EvtPhokhara::init_mode
void init_mode(EvtParticle *p)
Definition EvtPhokhara.cc:158
EvtPhokhara::decay
void decay(EvtParticle *p)
Definition EvtPhokhara.cc:648
EvtPhokhara::commandName
std::string commandName()
Definition EvtPhokhara.cc:623
EvtPhokhara::init_evt
void init_evt(EvtParticle *p)
Definition EvtPhokhara.cc:936
EvtPhokhara::PhokharaInit
void PhokharaInit(int dummy)
Definition EvtPhokhara.cc:925
EvtPhokhara::initProbMax
void initProbMax()
Definition EvtPhokhara.cc:156
EvtPhokhara::getName
void getName(std::string &name)
Definition EvtPhokhara.cc:152
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:173