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