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