BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara_pi0pi0pipi.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_pi0pi0pipi::nevtgen = 0;
57int EvtPhokhara_pi0pi0pipi::nphokharadecays = 0;
58EvtDecayBasePtr* EvtPhokhara_pi0pi0pipi::phokharadecays = 0;
59int EvtPhokhara_pi0pi0pipi::ntable = 0;
60
61int EvtPhokhara_pi0pi0pipi::ncommand = 0;
62int EvtPhokhara_pi0pi0pipi::lcommand = 0;
63std::string* EvtPhokhara_pi0pi0pipi::commands = 0;
64int EvtPhokhara_pi0pi0pipi::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_pi0pi0pipi::getName( std::string& model_name ) {
97
98 model_name = "PHOKHARA_pi0pi0pipi";
99}
100
102
104
106 m_pion = 2;
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(10)
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_pi0pi0pipi::commandName() { return std::string( "PhokharaPar" ); }
138
139void EvtPhokhara_pi0pi0pipi::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 : pi0pi0pi+pi- 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 std::cout << "FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate."
210 << std::endl;
211 //----
212storedEvents:
213 int more = 0;
214 int numstable = 0;
215 int numparton = 0;
216 EvtId evtnumstable[100]; //
217 EvtVector4R p4[20];
218
219 // except ISR photos, products depending on channel
220 if ( flags_.pion == 0 )
221 { // mu+ mu-
222 // mu+
223 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -13 );
224 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
225 ctes_.momenta[3][5] );
226 numstable++;
227 // mu -
228 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 13 );
229 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
230 ctes_.momenta[3][6] );
231 numstable++;
232 }
233 if ( flags_.pion == 1 )
234 { // pi+ pi-
235 // pi+
236 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
237 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
238 ctes_.momenta[3][5] );
239 numstable++;
240 // pi -
241 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
242 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
243 ctes_.momenta[3][6] );
244 numstable++;
245 }
246 if ( flags_.pion == 2 )
247 { // pi+ pi-2pi0
248 // pi0
249 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
250 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
251 ctes_.momenta[3][5] );
252 numstable++;
253 // pi0
254 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
255 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
256 ctes_.momenta[3][6] );
257 numstable++;
258 // pi-
259 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
260 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
261 ctes_.momenta[3][7] );
262 numstable++;
263 // pi +
264 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
265 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
266 ctes_.momenta[3][8] );
267 numstable++;
268 }
269 if ( flags_.pion == 3 )
270 { // 2(pi+ pi-)
271 // pi+
272 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
273 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
274 ctes_.momenta[3][5] );
275 numstable++;
276 // pi-
277 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
278 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
279 ctes_.momenta[3][6] );
280 numstable++;
281 // pi+
282 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
283 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
284 ctes_.momenta[3][7] );
285 numstable++;
286 // pi -
287 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
288 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
289 ctes_.momenta[3][8] );
290 numstable++;
291 }
292 if ( flags_.pion == 4 )
293 { // ppbar
294 // pbar
295 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2212 );
296 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
297 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],
302 ctes_.momenta[3][6] );
303 numstable++;
304 }
305 if ( flags_.pion == 5 )
306 { // nnbar
307 // pbar
308 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2112 );
309 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
310 ctes_.momenta[3][5] );
311 numstable++;
312 // p
313 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 2112 );
314 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
315 ctes_.momenta[3][6] );
316 numstable++;
317 }
318 if ( flags_.pion == 6 )
319 { // K+ K-
320 // K+
321 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 321 );
322 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
323 ctes_.momenta[3][5] );
324 numstable++;
325 // K -
326 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -321 );
327 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
328 ctes_.momenta[3][6] );
329 numstable++;
330 }
331 if ( flags_.pion == 7 )
332 { // K0K0bar
333 // Kbar
334 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 311 );
335 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
336 ctes_.momenta[3][5] );
337 numstable++;
338 // K0
339 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -311 );
340 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
341 ctes_.momenta[3][6] );
342 numstable++;
343 }
344 if ( flags_.pion == 8 )
345 { // pi+ pi-pi0
346 // pi+
347 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
348 p4[numstable].set( ctes_.momenta[0][5], ctes_.momenta[1][5], ctes_.momenta[2][5],
349 ctes_.momenta[3][5] );
350 numstable++;
351 // pi-
352 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
353 p4[numstable].set( ctes_.momenta[0][6], ctes_.momenta[1][6], ctes_.momenta[2][6],
354 ctes_.momenta[3][6] );
355 numstable++;
356 // pi0
357 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 111 );
358 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
359 ctes_.momenta[3][7] );
360 numstable++;
361 }
362 if ( flags_.pion == 9 )
363 { // Lambda Lambdabar-> pi+ pi- ppbar
364 // pi+
365 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 211 );
366 p4[numstable].set( ctes_.momenta[0][7], ctes_.momenta[1][7], ctes_.momenta[2][7],
367 ctes_.momenta[3][7] );
368 numstable++;
369 // pbar
370 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -2212 );
371 p4[numstable].set( ctes_.momenta[0][8], ctes_.momenta[1][8], ctes_.momenta[2][8],
372 ctes_.momenta[3][8] );
373 numstable++;
374 // pi-
375 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( -211 );
376 p4[numstable].set( ctes_.momenta[0][9], ctes_.momenta[1][9], ctes_.momenta[2][9],
377 ctes_.momenta[3][9] );
378 numstable++;
379 // p
380 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 2212 );
381 p4[numstable].set( ctes_.momenta[0][10], ctes_.momenta[1][10], ctes_.momenta[2][10],
382 ctes_.momenta[3][10] );
383 numstable++;
384 }
385
386 // ISR gamma
387 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 22 );
388 p4[numstable].set( ctes_.momenta[0][2], ctes_.momenta[1][2], ctes_.momenta[2][2],
389 ctes_.momenta[3][2] );
390 numstable++;
391 if ( ctes_.momenta[0][3] != 0 ) // second photon exists
392 {
393 evtnumstable[numstable] = EvtPDL::evtIdFromStdHep( 22 );
394 p4[numstable].set( ctes_.momenta[0][3], ctes_.momenta[1][3], ctes_.momenta[2][3],
395 ctes_.momenta[3][3] );
396 numstable++;
397 }
398
399 int channel = EvtDecayTable::inChannelList( p->getId(), numstable, evtnumstable );
400 more = ( channel != -1 );
401 if ( more )
402 {
403 std::cout << "Existence of mode " << channel
404 << " in exclusive decay list has the same final state as this one" << std::endl;
405 abort();
406 }
407
408 p->makeDaughters( numstable, evtnumstable );
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_pi0pi0pipi::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
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
457 m_pion = 2;
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}
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 set(int i, double d)