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