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