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