BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtLunda.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: EvtLunda.cc
12// the necessary file: lundar.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//
27#include "EvtGenBase/EvtPDL.hh"
30#include <string>
31#include "EvtGenBase/EvtId.hh"
32#include <iostream>
33#include <iomanip>
34#include <fstream>
35#include <string.h>
36#include <stdlib.h>
37#include <unistd.h>
38#include <stdio.h>
39using std::endl;
40using std::fstream;
41using std::ios;
42using std::ofstream;
43using std::resetiosflags;
44using std::setiosflags;
45using std::setw;
46
47
48int EvtLunda::nlundadecays=0;
49 EvtDecayBasePtr* EvtLunda::lundadecays=0;
50int EvtLunda::ntable=0;
51
52int EvtLunda::ncommand=0;
53int EvtLunda::lcommand=0;
54std::string* EvtLunda::commands=0;
55int EvtLunda::nevt=0;
56
57
58extern "C" {
59 extern void lundar_(int *, int *, float *,int *,int *,int *,
60 double *,double *,double *,double *);
61}
62
63extern "C" {
64 extern int lucomp_(int* kf);
65}
66
67extern "C" {
68 extern void lugive0_(const char *cnfgstr,int length);
69}
70
71//COMMON/CHECKTAG/DECAYTAG
72extern "C" struct {
73double decaytag;
75
76/*
77extern struct
78{
79 char mypdg[200];
80}mypdgfile_;
81*/
82
83/*
84extern struct
85{
86 int xpdg; // pdg code for e+e- -> X particle, see subroutine LUBECN( , ) in luarlw.F,
87}nores_;
88*/
89
92
93
94 int i;
95
96
97 //the deletion of commands is really uggly!
98
99 if (nlundadecays==0) {
100 delete [] commands;
101 commands=0;
102 return;
103 }
104
105 for(i=0;i<nlundadecays;i++){
106 if (lundadecays[i]==this){
107 lundadecays[i]=lundadecays[nlundadecays-1];
108 nlundadecays--;
109 if (nlundadecays==0) {
110 delete [] commands;
111 commands=0;
112 }
113 return;
114 }
115 }
116
117 report(ERROR,"EvtGen") << "Error in destroying Lunda model!"<<endl;
118
119}
120
121
122void EvtLunda::getName(std::string& model_name){
123
124 model_name="LUNDA";
125
126}
127
129
130 return new EvtLunda;
131
132}
133
134
136
137 noProbMax();
138
139}
140
141
143
144// checkNArg(1);
145 std::string strpdg=getenv("BESEVTGENROOT");
146 strpdg +="/share/r_pdg.dat";
147 //strcpy(mypdgfile_.mypdg, strpdg.c_str());
148 std::cout<<"mypdg= "<<strpdg<<std::endl;
149
150 if(getNArg()>2){std::cout<<"Parameter can be 1 or 2, You have "<<getNArg()<<std::endl; ::abort();}
151
152 if (getParentId().isAlias()){
153
154 report(ERROR,"EvtGen") << "EvtLunda finds that you are decaying the"<<endl
155 << " aliased particle "
156 << EvtPDL::name(getParentId()).c_str()
157 << " with the Lunda model"<<endl
158 << " this does not work, please modify decay table."
159 << endl;
160 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
161 ::abort();
162
163 }
164
165 store(this);
166
167}
168
169
171
172 return std::string("LundaPar");
173
174}
175
176void EvtLunda::command(std::string cmd){
177 if (ncommand==lcommand){
178 lcommand=10+2*lcommand;
179 std::string* newcommands=new std::string[lcommand];
180 int i;
181 for(i=0;i<ncommand;i++){
182 newcommands[i]=commands[i];
183 }
184 delete [] commands;
185 commands=newcommands;
186 }
187 commands[ncommand]=cmd;
188 ncommand++;
189}
190
191
192
194
195 static int iniflag=1;
196
197 static EvtId STRNG=EvtPDL::getId("string");
198
199 int istdheppar=EvtPDL::getStdHep(p->getId());
200
201 /*
202 int Xpdg=0;
203 if(getNArg()==1){
204 Xpdg = getArg(0);
205 if(Xpdg == 100443) Xpdg = 30443; //chage the curent pdg code to jetset pdg code for psiprime
206 }
207 nores_.xpdg = Xpdg;
208 */
209
210/*
211 if (pycomp_(&istdheppar)==0){
212 report(ERROR,"EvtGen") << "Lunda can not decay:"
213 <<EvtPDL::name(p->getId()).c_str()<<endl;
214 return;
215 }
216*/
217 double mp=p->mass();
218 float xmp=mp;
219// std::cout<<"float xmp="<<xmp<<std::endl;
220
221 EvtVector4R p4[50];
222
223 int i,more;
224 int ip=EvtPDL::getStdHep(p->getId());
225 int ndaugjs;
226 static int kf[4000];
227 EvtId evtnumstable[100],evtnumparton[100];
228 int stableindex[100],partonindex[100];
229 int numstable;
230 int numparton;
231 static int km[4000];
232 EvtId type[MAX_DAUG];
233
234 int isr=1; //open ISR (default)
235 if(getNArg()>0){ isr=getArg(0);}
236
237 static double px[4000],py[4000],pz[4000],e[4000];
238 if (iniflag==1) lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
239
240 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
241
242 int count=0;
243 bool mtg=0;
244
245 do{
246 //report(INFO,"EvtGen") << "calling lunda " << ip<< " " << mp <<endl;
247 iniflag=iniflag+1; //to count the event number
248 //std::cout<<"Event number by Lunda: "<<iniflag<<std::endl;
249modeSelection:
250 lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
251
252 mtg = checktag_.decaytag==1;
253 if(mtg)std::cout<<"checktag_.decaytag="<<checktag_.decaytag<<std::endl;
254 //if the Ecms is too low, Lunda fails to decay, so change to 3pi decay exclusively
255 //if(mtg) {ExclusiveDecay(p); return;} //see SUBROUTINE LUBEGN(KFL,ECM) in luarlw.F
256
257 LundaInit(0); // allow user to set LundaPar in the decay list
258 numstable=0;
259 numparton=0;
260 //report(INFO,"EvtGen") << "found some daughters " << ndaugjs << endl;
261 double esum=0;
262 //for debugging
263 //for(int i=0;i<ndaugjs;i++){
264 //std::cout<<"ndaugjs,kf,km,px,py,pz,e: "<<i<<", "<<km[i]<<", "<<EvtPDL::name(EvtPDL::evtIdFromStdHep(kf[i]))<<", "<<px[i]<<" ,"<<py[i]<<", "<<pz[i]<<", "<<e[i]<<std::endl;
265 //}
266
267 for(i=0;i<ndaugjs;i++){
268 if (abs(kf[i])==11 || kf[i]==92 || kf[i]==22) continue; //fill out the unstatble particle
269 //std::cout<<i<<", "<<km[i]<<", "<<kf[i]<<", "<<EvtPDL::name(EvtPDL::evtIdFromStdHep(kf[i]))<<" "<<px[i]<<" ,"<<py[i]<<", "<<pz[i]<<", "<<e[i]<<std::endl; //for debugging
270 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
271 report(ERROR,"EvtGen") << "Lunda returned particle:"<<kf[i]<<endl;
272 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
273 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
274 report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
275 //xmp=(1+0.01)*xmp;
276 std::cout<<"xmp= "<<xmp<<std::endl;
277 goto modeSelection;
278 }
279
280 //sort out the partons
281 if (abs(kf[i])<=6||kf[i]==21){
282 partonindex[numparton]=i;
283 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
284 numparton++;
285 }
286 else{
287 stableindex[numstable]=i;
288 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
289 numstable++;
290 }
291
292 esum+=e[i];
293 // have to protect against negative mass^2 for massless particles
294 // i.e. neutrinos and photons.
295 // this is uggly but I need to fix it right now....
296
297 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
298
299 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
300
301 }
302
303 p4[i].set(e[i],px[i],py[i],pz[i]);
304
305
306 }
307
308 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
309
310 // Test the branching fraction of lunda
311 // the specified decay mode is put as the 0-th channel with specifing mother particle
312 more=(channel!=-1);
313 //debugging
314 if(abs(esum-xmp)>0.001 ){
315 std::cout<<"Unexpected final states from Lunda with total energy "<<esum<<" unequal to "<<xmp<<std::endl;
316 //xmp=(1+0.01)*xmp;
317 std::cout<<"xmp= "<<xmp<<std::endl;
318 goto modeSelection;
319 }
320
321 count++;
322 }while( more && (count<10000) && mtg );
323 // }while( more && (count<10000) );
324
325 if (count>9999) {
326 report(INFO,"EvtGen") << "Too many loops in EvtLunda!!!"<<endl;
327 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
328 for(i=0;i<numstable;i++){
329 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
330 }
331
332 }
333
334
335
336 if (numparton==0){
337
338 p->makeDaughters(numstable,evtnumstable);
339 int ndaugFound=0;
340 for(i=0;i<numstable;i++){
341 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
342 ndaugFound++;
343 }
344 if ( ndaugFound == 0 ) {
345 report(ERROR,"EvtGen") << "Lunda has failed to do a decay ";
346 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
347 assert(0);
348 }
349
350 fixPolarizations(p);
351 //debugging
352 // p->printTree();
353
354 return ;
355
356 }
357 else{
358
359 //have partons in LUNDA
360
361 EvtVector4R p4string(0.0,0.0,0.0,0.0);
362
363 for(i=0;i<numparton;i++){
364 p4string+=p4[partonindex[i]];
365 }
366
367 int nprimary=1;
368 type[0]=STRNG;
369 for(i=0;i<numstable;i++){
370 if (km[stableindex[i]]==0){
371 type[nprimary++]=evtnumstable[i];
372 }
373 }
374
375 p->makeDaughters(nprimary,type);
376
377 p->getDaug(0)->init(STRNG,p4string);
378
379 EvtVector4R p4partons[10];
380
381 for(i=0;i<numparton;i++){
382 p4partons[i]=p4[partonindex[i]];
383 }
384
385 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
386
387
388
389 nprimary=1;
390
391 for(i=0;i<numstable;i++){
392
393 if (km[stableindex[i]]==0){
394 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
395 }
396 }
397
398
399 int nsecond=0;
400 for(i=0;i<numstable;i++){
401 if (km[stableindex[i]]!=0){
402 type[nsecond++]=evtnumstable[i];
403 }
404 }
405
406
407 p->getDaug(0)->makeDaughters(nsecond,type);
408
409 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
410 -p4string.get(2),-p4string.get(3));
411
412 nsecond=0;
413 for(i=0;i<numstable;i++){
414 if (km[stableindex[i]]!=0){
415 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
416 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
417 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
418 p->getDaug(0)->getDaug(nsecond)->decay();
419 nsecond++;
420 }
421 }
422
423 if ( nsecond == 0 ) {
424 report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
425 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
426 assert(0);
427 }
428
429 fixPolarizations(p);
430
431 return ;
432
433 }
434
435}
436
437void EvtLunda::fixPolarizations(EvtParticle *p){
438
439 //special case for now to handle the J/psi polarization
440
441 int ndaug=p->getNDaug();
442
443 int i;
444
445 static EvtId Jpsi=EvtPDL::getId("J/psi");
446
447 for(i=0;i<ndaug;i++){
448 if(p->getDaug(i)->getId()==Jpsi){
449
450 EvtSpinDensity rho;
451
452 rho.SetDim(3);
453 rho.Set(0,0,0.5);
454 rho.Set(0,1,0.0);
455 rho.Set(0,2,0.0);
456
457 rho.Set(1,0,0.0);
458 rho.Set(1,1,1.0);
459 rho.Set(1,2,0.0);
460
461 rho.Set(2,0,0.0);
462 rho.Set(2,1,0.0);
463 rho.Set(2,2,0.5);
464
465 EvtVector4R p4Psi=p->getDaug(i)->getP4();
466
467 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
468 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
469
470
473
474 }
475 }
476
477}
478
479void EvtLunda::store(EvtDecayBase* jsdecay){
480
481 if (nlundadecays==ntable){
482
483 EvtDecayBasePtr* newlundadecays=new EvtDecayBasePtr[2*ntable+10];
484 int i;
485 for(i=0;i<ntable;i++){
486 newlundadecays[i]=lundadecays[i];
487 }
488 ntable=2*ntable+10;
489 delete [] lundadecays;
490 lundadecays=newlundadecays;
491 }
492
493 lundadecays[nlundadecays++]=jsdecay;
494
495
496
497}
498
499
500void EvtLunda::LundaInit(int dummy){
501
502 static int first=1;
503 if (first){
504
505 first=0;
506 for(int i=0;i<ncommand;i++)
507 lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
508 }
509
510}
511
513 EvtId daugs[8];
514 int _ndaugs =4;
515 daugs[0]=EvtPDL::getId(std::string("pi+"));
516 daugs[1]=EvtPDL::getId(std::string("pi-"));
517 daugs[2]=EvtPDL::getId(std::string("pi+"));
518 daugs[3]=EvtPDL::getId(std::string("pi-"));
519 checkA:
520 p->makeDaughters(_ndaugs,daugs);
521 double totMass=0;
522 for(int i=0;i< _ndaugs;i++){
523 EvtParticle* di=p->getDaug(i);
524 double xmi=EvtPDL::getMass(di->getId());
525 di->setMass(xmi);
526 totMass+=xmi;
527 }
528 if(totMass>p->mass()) goto checkA;
529 p->initializePhaseSpace( _ndaugs , daugs);
530
531
532}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
int lucomp_(int *kf)
void lundar_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *)
struct @9 checktag_
void lugive0_(const char *cnfgstr, int length)
double decaytag
Definition: EvtLunda.cc:73
const int MAX_DAUG
Definition: EvtParticle.hh:38
DOUBLE_PRECISION count[3]
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
@ ERROR
Definition: EvtReport.hh:49
@ INFO
Definition: EvtReport.hh:52
const double alpha
double getArg(int j)
void noProbMax()
EvtId getParentId()
Definition: EvtDecayBase.hh:60
void setDaughterSpinDensity(int daughter)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition: EvtId.hh:27
void LundaInit(int dummy)
Definition: EvtLunda.cc:500
void initProbMax()
Definition: EvtLunda.cc:135
std::string commandName()
Definition: EvtLunda.cc:170
void decay(EvtParticle *p)
Definition: EvtLunda.cc:193
void ExclusiveDecay(EvtParticle *p)
Definition: EvtLunda.cc:512
virtual ~EvtLunda()
Definition: EvtLunda.cc:91
void command(std::string cmd)
Definition: EvtLunda.cc:176
EvtDecayBase * clone()
Definition: EvtLunda.cc:128
void init()
Definition: EvtLunda.cc:142
void getName(std::string &name)
Definition: EvtLunda.cc:122
EvtLunda()
Definition: EvtLunda.cc:90
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
void setMass(double m)
Definition: EvtParticle.hh:372
void setSpinDensityForwardHelicityBasis(const EvtSpinDensity &rho)
Definition: EvtParticle.cc:178
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void decay()
Definition: EvtParticle.cc:404
EvtId getId() const
Definition: EvtParticle.cc:113
void setDiagonalSpinDensity()
Definition: EvtParticle.cc:133
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
void deleteDaughters(bool keepChannel=false)
Definition: EvtParticle.cc:540
double mass() const
Definition: EvtParticle.cc:127
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
void Set(int i, int j, const EvtComplex &rhoij)
void SetDim(int n)
double get(int i) const
Definition: EvtVector4R.hh:179
double d3mag() const
Definition: EvtVector4R.cc:186
void set(int i, double d)
Definition: EvtVector4R.hh:183
const double mp
Definition: incllambda.cxx:45