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