BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara_etapgam.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
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_etapgam::nevtgen=0;
58int EvtPhokhara_etapgam::nphokharadecays=0;
59EvtDecayBasePtr* EvtPhokhara_etapgam::phokharadecays=0;
60int EvtPhokhara_etapgam::ntable=0;
61
62int EvtPhokhara_etapgam::ncommand=0;
63int EvtPhokhara_etapgam::lcommand=0;
64std::string* EvtPhokhara_etapgam::commands=0;
65int EvtPhokhara_etapgam::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_etapgam::getName(std::string& model_name){
96
97 model_name="PHOKHARA_etapgam";
98
99}
100
106
107
109
110 noProbMax();
111
112}
113
114
116 m_pion=15;
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 // pi+pi-eta
122#include "Phokhara_init_mode.txt"
123}
124
125
126
128 checkNArg(0);
129
130 std::string locvp=getenv("BESEVTGENROOT");
131 system("cat $BESEVTGENROOT/share/phokhara_10.0.param>phokhara_10.0.param");
132 system("cat $BESEVTGENROOT/share/phokhara_10.0.fferr>phokhara_10.0.fferr");
133 system("cat $BESEVTGENROOT/share/phokhara_10.0.ffwarn>phokhara_10.0.ffwarn");
134
135
136 if (getParentId().isAlias()){
137
138 report(ERROR,"EvtGen") << "EvtPhokhara finds that you are decaying the"<<endl
139 << " aliased particle "
140 << EvtPDL::name(getParentId()).c_str()
141 << " with the Phokhara model"<<endl
142 << " this does not work, please modify decay table."
143 << endl;
144 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
145 ::abort();
146
147 }
148
149 store(this);
150
151}
152
153
155
156 return std::string("PhokharaPar");
157
158}
159
160void EvtPhokhara_etapgam::command(std::string cmd){
161
162 if (ncommand==lcommand){
163
164 lcommand=10+2*lcommand;
165
166 std::string* newcommands=new std::string[lcommand];
167
168 int i;
169
170 for(i=0;i<ncommand;i++){
171 newcommands[i]=commands[i];
172 }
173
174 delete [] commands;
175
176 commands=newcommands;
177
178 }
179
180 commands[ncommand]=cmd;
181
182 ncommand++;
183
184}
185
186
187
189 EvtId myvpho=EvtPDL::getId("vpho");
190 if(p->getId()!=myvpho) {std::cout<<"Parent particle is required to be vpho for Phokhara model"<<std::endl;abort();}
191 if(nevtgen==0) {init_mode(p);}
192 else{init_evt(p);}
193
194 std::cout<<"PHOKHARA : eta' gamma mode "<<std::endl;
195 int istdheppar=EvtPDL::getStdHep(p->getId());
196 int ntrials = 0;
197 int tr_old[3];
198 tr_old[0] = (int)maxima_.tr[0];
199 tr_old[1] = (int)maxima_.tr[1];
200 tr_old[2] = (int)maxima_.tr[2];
201
202 while( ntrials < 1000000)
203 {
204 ievent++;
205 RANLXDF(Ar_r,1);
206 Ar[1] = Ar_r[0];
207
208 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]+maxima_.Mmax[2]))) {
209 maxima_.count[0] = maxima_.count[0]+1.0;
210 GEN_0PH(2,qqmin,ctes_.Sp,cos3min,cos3max);
211 }else
212 if (Ar[1] <= ( (maxima_.Mmax[0]+maxima_.Mmax[1])/(maxima_.Mmax[0]+maxima_.Mmax[1]+maxima_.Mmax[2]))) {
213 maxima_.count[1] = maxima_.count[1]+1.0;
214 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
215 }
216 else {
217 maxima_.count[2] = maxima_.count[2]+1.0;
218 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
219 }
220
221 if( ((int)maxima_.tr[0]+(int)maxima_.tr[1]+(int)maxima_.tr[2]) > (tr_old[0]+tr_old[1]+tr_old[2]) ) // event accepted after cuts
222 {
223 goto storedEvents;
224 }
225 ntrials ++;
226 }
227
228 std::cout <<"FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate." <<std::endl;
229 //----
230 storedEvents:
231 int more=0;
232 int numstable=0;
233 int numparton=0;
234 EvtId evtnumstable[100];//
235 EvtVector4R p4[20];
236
237 // except ISR photos, products depending on channel
238 if (flags_.pion == 0) { // mu+ mu-
239 // mu+
240 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-13);
241 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
242 numstable++;
243 // mu -
244 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(13);
245 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
246 numstable++;
247 }
248 if (flags_.pion == 1) { // pi+ pi-
249 // pi+
250 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
251 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
252 numstable++;
253 // pi -
254 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
255 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
256 numstable++;
257 }
258 if (flags_.pion == 2) { // pi+ pi-2pi0
259 // pi0
260 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
261 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
262 numstable++;
263 // pi0
264 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
265 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
266 numstable++;
267 // pi-
268 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
269 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
270 numstable++;
271 // pi +
272 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
273 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
274 numstable++;
275 }
276 if (flags_.pion == 3) { // 2(pi+ pi-)
277 // pi+
278 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
279 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
280 numstable++;
281 // pi-
282 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
283 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
284 numstable++;
285 // pi+
286 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
287 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
288 numstable++;
289 // pi -
290 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
291 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
292 numstable++;
293 }
294 if (flags_.pion == 4) { // ppbar
295 // pbar
296 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
297 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
298 numstable++;
299 // p
300 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
301 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
302 numstable++;
303 }
304 if (flags_.pion == 5) { // nnbar
305 // pbar
306 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2112);
307 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
308 numstable++;
309 // p
310 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2112);
311 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
312 numstable++;
313 }
314 if (flags_.pion == 6) { // K+ K-
315 // K+
316 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(321);
317 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
318 numstable++;
319 // K -
320 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-321);
321 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
322 numstable++;
323 }
324 if (flags_.pion == 7) { // K0K0bar
325 // Kbar
326 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(311);
327 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
328 numstable++;
329 // K0
330 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-311);
331 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
332 numstable++;
333 }
334 if (flags_.pion == 8) { // pi+ pi-pi0
335 // pi+
336 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
337 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
338 numstable++;
339 // pi-
340 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
341 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
342 numstable++;
343 // pi0
344 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
345 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
346 numstable++;
347 }
348 if (flags_.pion == 9) { //Lambda Lambdabar-> pi+ pi- ppbar
349 // pi+
350 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
351 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
352 numstable++;
353 // pbar
354 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
355 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
356 numstable++;
357 // pi-
358 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
359 p4[numstable].set(ctes_.momenta[0][9],ctes_.momenta[1][9], ctes_.momenta[2][9], ctes_.momenta[3][9]);
360 numstable++;
361 // p
362 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
363 p4[numstable].set(ctes_.momenta[0][10],ctes_.momenta[1][10], ctes_.momenta[2][10], ctes_.momenta[3][10]);
364 numstable++;
365 }
366 if (flags_.pion == 13) { // pi0 gamma
367 // pi0
368 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
369 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
370 numstable++;
371 // gamma
372 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
373 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
374 numstable++;
375 }
376 if (flags_.pion == 14) { // eta gamma
377 // eta
378 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(221);
379 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
380 numstable++;
381 // gamma
382 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
383 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
384 numstable++;
385 }
386 if (flags_.pion == 15) { // eta' gamma
387 // eta'
388 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(331);
389 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
390 numstable++;
391 // gamma
392 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
393 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
394 numstable++;
395 }
396
397 // ISR gamma
398 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
399 p4[numstable].set(ctes_.momenta[0][2],ctes_.momenta[1][2], ctes_.momenta[2][2], ctes_.momenta[3][2]);
400 numstable++;
401 if( ctes_.momenta[0][3] != 0 ) // second photon exists
402 {
403 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
404 p4[numstable].set(ctes_.momenta[0][3],ctes_.momenta[1][3], ctes_.momenta[2][3], ctes_.momenta[3][3]);
405 numstable++;
406 }
407
408 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
409 more=(channel!=-1);
410 if(more) {std::cout<<"Existence of mode "<<channel<<" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
411
412 p->makeDaughters(numstable,evtnumstable);
413 //double weight = p->initializePhaseSpace(getNDaug(),getDaugs());
414
415 int ndaugFound=0;
416 EvtVector4R SUMP4(0,0,0,0);
417 for(int i=0;i<numstable;i++){
418 p->getDaug(i)->init(evtnumstable[i],p4[i]);
419 ndaugFound++;
420 }
421 if ( ndaugFound == 0 ) {
422 report(ERROR,"EvtGen") << "Phokhara has failed to do a decay ";
423 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
424 assert(0);
425 }
426
427 nevtgen++;
428 return ;
429
430}
431
432
433
434void EvtPhokhara_etapgam::store(EvtDecayBase* jsdecay){
435
436 if (nphokharadecays==ntable){
437
438 EvtDecayBasePtr* newphokharadecays=new EvtDecayBasePtr[2*ntable+10];
439 int i;
440 for(i=0;i<ntable;i++){
441 newphokharadecays[i]=phokharadecays[i];
442 }
443 ntable=2*ntable+10;
444 delete [] phokharadecays;
445 phokharadecays=newphokharadecays;
446 }
447
448 phokharadecays[nphokharadecays++]=jsdecay;
449
450
451
452}
453
454
456 static int first=1;
457 if (first){
458
459 first=0;
460 //for(int i=0;i<ncommand;i++)
461 // lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
462 }
463
464}
465
466
467
469 m_pion=15;
470 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
471 // 2pi+2pi-(3),ppbar(4),nnbar(5),
472 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
473 // Lamb Lambbar->pi-pi+ppbar(9)
474 // pi+pi-eta
475#include "Phokhara_init_evt.txt"
476}
477
478
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
#define RANLXDF(AR, VAL)
#define GEN_0PH(I, QQMIN, SP, COS3MIN, COS3MAX)
struct @15 ctes_
struct @21 maxima_
struct @22 flags_
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:36
@ ERROR
Definition EvtReport.hh:49
EvtId getParentId()
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
EvtParticle * getDaug(int i)
double mass() const
void command(std::string cmd)
void getName(std::string &name)
void init_evt(EvtParticle *p)
void decay(EvtParticle *p)
void init_mode(EvtParticle *p)
double double double * p4
Definition qcdloop1.h:77