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