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