BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara_pipi.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_pipi.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_pipi::nevtgen = 0;
57int EvtPhokhara_pipi::nphokharadecays = 0;
58EvtDecayBasePtr* EvtPhokhara_pipi::phokharadecays = 0;
59int EvtPhokhara_pipi::ntable = 0;
60
61int EvtPhokhara_pipi::ncommand = 0;
62int EvtPhokhara_pipi::lcommand = 0;
63std::string* EvtPhokhara_pipi::commands = 0;
64int EvtPhokhara_pipi::nevt = 0;
65
66EvtPhokhara_pipi::EvtPhokhara_pipi() {}
67EvtPhokhara_pipi::~EvtPhokhara_pipi() {
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_pipi::getName( std::string& model_name ) { model_name = "PHOKHARA_pipi"; }
97
98EvtDecayBase* EvtPhokhara_pipi::clone() { return new EvtPhokhara_pipi; }
99
100void EvtPhokhara_pipi::initProbMax() { noProbMax(); }
101
102void EvtPhokhara_pipi::init_mode( EvtParticle* p ) {
103 m_pion = 1;
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 // pi+pi-eta
109#include "Phokhara_init_mode.txt"
110}
111
112void EvtPhokhara_pipi::init() {
113 checkNArg( 0 );
114
115 std::string locvp = getenv( "BESEVTGENROOT" );
116 system( "cat $BESEVTGENROOT/share/phokhara_10.0.param>phokhara_10.0.param" );
117 system( "cat $BESEVTGENROOT/share/phokhara_10.0.fferr>phokhara_10.0.fferr" );
118 system( "cat $BESEVTGENROOT/share/phokhara_10.0.ffwarn>phokhara_10.0.ffwarn" );
119
120 if ( getParentId().isAlias() )
121 {
122
123 report( ERROR, "EvtGen" ) << "EvtPhokhara finds that you are decaying the" << endl
124 << " aliased particle " << EvtPDL::name( getParentId() ).c_str()
125 << " with the Phokhara model" << endl
126 << " this does not work, please modify decay table." << endl;
127 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
128 ::abort();
129 }
130
131 store( this );
132}
133
134std::string EvtPhokhara_pipi::commandName() { return std::string( "PhokharaPar" ); }
135
136void EvtPhokhara_pipi::command( std::string cmd ) {
137
138 if ( ncommand == lcommand )
139 {
140
141 lcommand = 10 + 2 * lcommand;
142
143 std::string* newcommands = new std::string[lcommand];
144
145 int i;
146
147 for ( i = 0; i < ncommand; i++ ) { newcommands[i] = commands[i]; }
148
149 delete[] commands;
150
151 commands = newcommands;
152 }
153
154 commands[ncommand] = cmd;
155
156 ncommand++;
157}
158
159void EvtPhokhara_pipi::decay( EvtParticle* p ) {
160 EvtId myvpho = EvtPDL::getId( "vpho" );
161 if ( p->getId() != myvpho )
162 {
163 std::cout << "Parent particle is required to be vpho for Phokhara model" << std::endl;
164 abort();
165 }
166 if ( nevtgen == 0 ) { init_mode( p ); }
167 else { init_evt( p ); }
168
169 std::cout << "PHOKHARA : pi+pi- mode " << std::endl;
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
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_pipi::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_pipi::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_pipi::init_evt( EvtParticle* p ) {
456 m_pion = 1;
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 // pi+pi-eta
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_pipi.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_pipi::PhokharaInit
void PhokharaInit(int dummy)
Definition EvtPhokhara_pipi.cc:444
EvtPhokhara_pipi::init_evt
void init_evt(EvtParticle *p)
Definition EvtPhokhara_pipi.cc:455
EvtPhokhara_pipi::commandName
std::string commandName()
Definition EvtPhokhara_pipi.cc:134
EvtPhokhara_pipi::initProbMax
void initProbMax()
Definition EvtPhokhara_pipi.cc:100
EvtPhokhara_pipi::getName
void getName(std::string &name)
Definition EvtPhokhara_pipi.cc:96
EvtPhokhara_pipi::init
void init()
Definition EvtPhokhara_pipi.cc:112
EvtPhokhara_pipi::init_mode
void init_mode(EvtParticle *p)
Definition EvtPhokhara_pipi.cc:102
EvtPhokhara_pipi::command
void command(std::string cmd)
Definition EvtPhokhara_pipi.cc:136
EvtPhokhara_pipi::~EvtPhokhara_pipi
virtual ~EvtPhokhara_pipi()
Definition EvtPhokhara_pipi.cc:67
EvtPhokhara_pipi::clone
EvtDecayBase * clone()
Definition EvtPhokhara_pipi.cc:98
EvtPhokhara_pipi::EvtPhokhara_pipi
EvtPhokhara_pipi()
Definition EvtPhokhara_pipi.cc:66
EvtPhokhara_pipi::decay
void decay(EvtParticle *p)
Definition EvtPhokhara_pipi.cc:159
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:173