BOSS 7.0.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara_ppbar.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_ppbar::nevtgen=0;
58int EvtPhokhara_ppbar::nphokharadecays=0;
59EvtDecayBasePtr* EvtPhokhara_ppbar::phokharadecays=0;
60int EvtPhokhara_ppbar::ntable=0;
61
62int EvtPhokhara_ppbar::ncommand=0;
63int EvtPhokhara_ppbar::lcommand=0;
64std::string* EvtPhokhara_ppbar::commands=0;
65int EvtPhokhara_ppbar::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_ppbar::getName(std::string& model_name){
96
97 model_name="PHOKHARA_ppbar";
98
99}
100
102
103 return new EvtPhokhara_ppbar;
104
105}
106
107
109
110 noProbMax();
111
112}
113
114
116 m_pion=4;
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 = 123456789;
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
529void EvtPhokhara_ppbar::command(std::string cmd){
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 std::cout<<"PHOKHARA : ppbar mode "<<std::endl;
563
564
565 int istdheppar=EvtPDL::getStdHep(p->getId());
566 int ntrials = 0;
567 int tr_old[2];
568 tr_old[0] = (int)maxima_.tr[0];
569 tr_old[1] = (int)maxima_.tr[1];
570
571 while( ntrials < 10000)
572 {
573 ievent++;
574 RANLXDF(Ar_r,1);
575 Ar[1] = Ar_r[0];
576
577 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
578 maxima_.count[0] = maxima_.count[0]+1.0;
579 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
580 }
581 else {
582 maxima_.count[1] = maxima_.count[1]+1.0;
583 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
584 }
585
586 if( ((int)maxima_.tr[0]+(int)maxima_.tr[1]) > (tr_old[0]+tr_old[1]) ) // event accepted after cuts
587 {
588 goto storedEvents;
589 }
590 ntrials ++;
591 }
592
593 std::cout <<"FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate." <<std::endl;
594 //----
595 storedEvents:
596 int more=0;
597 int numstable=0;
598 int numparton=0;
599 EvtId evtnumstable[100];//
600 EvtVector4R p4[20];
601
602 // except ISR photos, products depending on channel
603 if (flags_.pion == 0) { // mu+ mu-
604 // mu+
605 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-13);
606 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
607 numstable++;
608 // mu -
609 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(13);
610 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
611 numstable++;
612 }
613 if (flags_.pion == 1) { // pi+ pi-
614 // pi+
615 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
616 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
617 numstable++;
618 // pi -
619 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
620 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
621 numstable++;
622 }
623 if (flags_.pion == 2) { // pi+ pi-2pi0
624 // pi0
625 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
626 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
627 numstable++;
628 // pi0
629 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
630 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
631 numstable++;
632 // pi-
633 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
634 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
635 numstable++;
636 // pi +
637 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
638 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
639 numstable++;
640 }
641 if (flags_.pion == 3) { // 2(pi+ pi-)
642 // pi+
643 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
644 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
645 numstable++;
646 // pi-
647 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
648 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
649 numstable++;
650 // pi+
651 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
652 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
653 numstable++;
654 // pi -
655 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
656 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
657 numstable++;
658 }
659 if (flags_.pion == 4) { // ppbar
660 // pbar
661 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
662 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
663 numstable++;
664 // p
665 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
666 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
667 numstable++;
668 }
669 if (flags_.pion == 5) { // nnbar
670 // pbar
671 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2112);
672 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
673 numstable++;
674 // p
675 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2112);
676 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
677 numstable++;
678 }
679 if (flags_.pion == 6) { // K+ K-
680 // K+
681 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(321);
682 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
683 numstable++;
684 // K -
685 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-321);
686 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
687 numstable++;
688 }
689 if (flags_.pion == 7) { // K0K0bar
690 // Kbar
691 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(311);
692 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
693 numstable++;
694 // K0
695 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-311);
696 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
697 numstable++;
698 }
699 if (flags_.pion == 8) { // pi+ pi-pi0
700 // pi+
701 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
702 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
703 numstable++;
704 // pi-
705 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
706 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
707 numstable++;
708 // pi0
709 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
710 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
711 numstable++;
712 }
713 if (flags_.pion == 9) { //Lambda Lambdabar-> pi+ pi- ppbar
714 // pi+
715 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
716 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
717 numstable++;
718 // pbar
719 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
720 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
721 numstable++;
722 // pi-
723 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
724 p4[numstable].set(ctes_.momenta[0][9],ctes_.momenta[1][9], ctes_.momenta[2][9], ctes_.momenta[3][9]);
725 numstable++;
726 // p
727 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
728 p4[numstable].set(ctes_.momenta[0][10],ctes_.momenta[1][10], ctes_.momenta[2][10], ctes_.momenta[3][10]);
729 numstable++;
730 }
731
732 // ISR gamma
733 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
734 p4[numstable].set(ctes_.momenta[0][2],ctes_.momenta[1][2], ctes_.momenta[2][2], ctes_.momenta[3][2]);
735 numstable++;
736 if( ctes_.momenta[0][3] != 0 ) // second photon exists
737 {
738 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
739 p4[numstable].set(ctes_.momenta[0][3],ctes_.momenta[1][3], ctes_.momenta[2][3], ctes_.momenta[3][3]);
740 numstable++;
741 }
742
743 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
744 more=(channel!=-1);
745 if(more) {std::cout<<"Existence of mode "<<channel<<" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
746
747 p->makeDaughters(numstable,evtnumstable);
748 //double weight = p->initializePhaseSpace(getNDaug(),getDaugs());
749
750 int ndaugFound=0;
751 EvtVector4R SUMP4(0,0,0,0);
752 for(int i=0;i<numstable;i++){
753 p->getDaug(i)->init(evtnumstable[i],p4[i]);
754 ndaugFound++;
755 }
756 if ( ndaugFound == 0 ) {
757 report(ERROR,"EvtGen") << "Phokhara has failed to do a decay ";
758 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
759 assert(0);
760 }
761
762 nevtgen++;
763 return ;
764
765}
766
767
768
769void EvtPhokhara_ppbar::store(EvtDecayBase* jsdecay){
770
771 if (nphokharadecays==ntable){
772
773 EvtDecayBasePtr* newphokharadecays=new EvtDecayBasePtr[2*ntable+10];
774 int i;
775 for(i=0;i<ntable;i++){
776 newphokharadecays[i]=phokharadecays[i];
777 }
778 ntable=2*ntable+10;
779 delete [] phokharadecays;
780 phokharadecays=newphokharadecays;
781 }
782
783 phokharadecays[nphokharadecays++]=jsdecay;
784
785
786
787}
788
789
791 static int first=1;
792 if (first){
793
794 first=0;
795 //for(int i=0;i<ncommand;i++)
796 // lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
797 }
798
799}
800
801
802
804 m_pion=4;
805 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
806 // 2pi+2pi-(3),ppbar(4),nnbar(5),
807 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
808 // Lamb Lambbar->pi-pi+ppbar(9)
809 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
810 EvtId myvpho=EvtPDL::getId("vpho");
811 m_E = p->mass();//EvtPDL::getMeanMass(myvpho);
812///======== list parameters to be initialized
813 m_tagged = 0;
814 m_nm = 50000 ; // # of events to determine the maximum
815 m_nlo = 1; // Born(0), NLO(1)
816 m_w = 0.0001; // soft photon cutoff
817 m_fsr = 1; // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
818 m_fsrnlo = 1 ; // yes(1), no(0)
819 m_NarrowRes = 0 ;// none(0) jpsi(1) psip(2)
820 m_FF_Kaon = 1 ; // constrained (0), unconstrained (1), Kuhn-Khodjamirian-Bruch (2)
821 m_ivac = 0; // yes(1), no(0)
822 m_FF_Pion = 0 ; // KS Pionformfactor(0), GS Pionformfactor old(1) and new(2)
823 m_f0_model = 0 ; // f0+f0(600): KK model(0), no structure(1), no f0+f0(600)(2), f0 KLOE(3)
824 m_q2min = 0.0; // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
825 m_q2_min_c = 0.0447 ; // minimal inv. mass squared of the hadrons(muons)(GeV^2)
826 m_q2_max_c = m_E*m_E; // maximal inv. mass squared of the hadrons(muons)(GeV^2)
827 m_gmin = 0.001; // minimal photon energy/missing energy (GeV)
828 m_phot1cut = 0.0; // maximal photon angle/missing momentum angle (grad)
829 m_phot2cut = 180.0; // maximal photon angle/missing momentum angle (grad)
830 m_pi1cut = 0.0 ; // minimal hadrons(muons) angle (grad)
831 m_pi2cut = 180.0; // maximal hadrons(muons) angle (grad)
832
833 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
834 if( m_pion==9 ){m_nlo = 0 ;}
835 // --- input parameter initialization -----------
836 maxima_.iprint = 0;
837 flags_.nlo = m_nlo;
838 flags_.pion = m_pion;
839 flags_.fsr = m_fsr;
840 flags_.fsrnlo = m_fsrnlo;
841 flags_.ivac = m_ivac;
842 flags_.FF_pion = m_FF_Pion;
843 flags_.f0_model = m_f0_model;
844 flags_.FF_kaon = m_FF_Kaon;
845 flags_.narr_res = m_NarrowRes;
846
847 ctes_.Sp = m_E*m_E; ;
848
849 cuts_.w = m_w;
850 cuts_.q2min = m_q2min;
851 cuts_.q2_min_c = m_q2_min_c;
852 cuts_.q2_max_c = m_q2_max_c;
853 cuts_.gmin = m_gmin;
854 cuts_.phot1cut = m_phot1cut;
855 cuts_.phot2cut = m_phot2cut;
856 cuts_.pi1cut = m_pi1cut;
857 cuts_.pi2cut = m_pi2cut;
858
859 INPUT();
860
861 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0); // photon1 angle cuts in the
862 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0); // LAB rest frame
863 cos2min = -1.0; // photon2 angle limits
864 cos2max = 1.0; //
865 cos3min = -1.0; // hadrons/muons angle limits
866 cos3max = 1.0; // in their rest frame
867 if (flags_.pion == 0) // virtual photon energy cut
868 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
869 else if (flags_.pion == 1)
870 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
871 else if (flags_.pion == 2)
872 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
873 else if (flags_.pion == 3)
874 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
875 else if (flags_.pion == 4)
876 qqmin = 4.0*ctes_.mp*ctes_.mp;
877 else if (flags_.pion == 5)
878 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
879 else if (flags_.pion == 6)
880 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
881 else if (flags_.pion == 7)
882 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
883 else if (flags_.pion == 8)
884 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
885 else if (flags_.pion == 9)
886 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
887 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin; // if only one photon
888 if (cuts_.q2_max_c < qqmax)
889 qqmax=cuts_.q2_max_c; // external cuts
890
891 // -------------------
892 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
893 qqmin = cuts_.q2_min_c;
894 else {
895 }
896
897
898 // =================================================
899// --- finding the maximum -------------------------
900 for(int i = 0; i<2; i++)
901 {
902 maxima_.Mmax[i] = 1.0;
903 maxima_.gross[i] = 0.0;
904 maxima_.klein[i] = 0.0;
905 }
906
907 if (flags_.nlo == 0)
908 maxima_.Mmax[1]=0.0; // only 1 photon events generated
909
910 maxima_.tr[0] = 0.0;
911 maxima_.tr[1] = 0.0;
912 maxima_.count[0] = 0.0;
913 maxima_.count[1] = 0.0;
914
915 // =================================================
916 // --- for the second run ---
917 maxima_.Mmax[0] = theMmax0;
918 maxima_.Mmax[1] = theMmax1;
919 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
920 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
921 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
922 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
923
924 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
925 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
926
927 if((flags_.pion == 2) || (flags_.pion == 3)){
928 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
929 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
930 }
931
932 if(flags_.pion == 8){
933 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
934 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
935 }
936// --- end of the second run -----------------------
937
938 maxima_.tr[0] = 0.0;
939 maxima_.tr[1] = 0.0;
940 maxima_.count[0] = 0.0;
941 maxima_.count[1] = 0.0;
942}
943
944
#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 command(std::string cmd)
EvtDecayBase * clone()
void decay(EvtParticle *p)
void init_evt(EvtParticle *p)
void getName(std::string &name)
void PhokharaInit(int dummy)
virtual ~EvtPhokhara_ppbar()
std::string commandName()
void init_mode(EvtParticle *p)
void set(int i, double d)
Definition: EvtVector4R.hh:183