00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenBase/EvtPatches.hh"
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtStringParticle.hh"
00026 #include "EvtGenBase/EvtDecayTable.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenModels/EvtLunda.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030 #include <string>
00031 #include "EvtGenBase/EvtId.hh"
00032 #include <iostream>
00033 #include <iomanip>
00034 #include <fstream>
00035 #include <string.h>
00036 #include <stdlib.h>
00037 #include <unistd.h>
00038 #include <stdio.h>
00039 using std::endl;
00040 using std::fstream;
00041 using std::ios;
00042 using std::ofstream;
00043 using std::resetiosflags;
00044 using std::setiosflags;
00045 using std::setw;
00046
00047
00048 int EvtLunda::nlundadecays=0;
00049 EvtDecayBasePtr* EvtLunda::lundadecays=0;
00050 int EvtLunda::ntable=0;
00051
00052 int EvtLunda::ncommand=0;
00053 int EvtLunda::lcommand=0;
00054 std::string* EvtLunda::commands=0;
00055 int EvtLunda::nevt=0;
00056
00057
00058 extern "C" {
00059 extern void lundar_(int *, int *, float *,int *,int *,int *,
00060 double *,double *,double *,double *);
00061 }
00062
00063 extern "C" {
00064 extern int lucomp_(int* kf);
00065 }
00066
00067 extern "C" {
00068 extern void lugive0_(const char *cnfgstr,int length);
00069 }
00070
00071
00072 extern "C" struct {
00073 double decaytag;
00074 } checktag_;
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 EvtLunda::EvtLunda(){}
00091 EvtLunda::~EvtLunda(){
00092
00093
00094 int i;
00095
00096
00097
00098
00099 if (nlundadecays==0) {
00100 delete [] commands;
00101 commands=0;
00102 return;
00103 }
00104
00105 for(i=0;i<nlundadecays;i++){
00106 if (lundadecays[i]==this){
00107 lundadecays[i]=lundadecays[nlundadecays-1];
00108 nlundadecays--;
00109 if (nlundadecays==0) {
00110 delete [] commands;
00111 commands=0;
00112 }
00113 return;
00114 }
00115 }
00116
00117 report(ERROR,"EvtGen") << "Error in destroying Lunda model!"<<endl;
00118
00119 }
00120
00121
00122 void EvtLunda::getName(std::string& model_name){
00123
00124 model_name="LUNDA";
00125
00126 }
00127
00128 EvtDecayBase* EvtLunda::clone(){
00129
00130 return new EvtLunda;
00131
00132 }
00133
00134
00135 void EvtLunda::initProbMax(){
00136
00137 noProbMax();
00138
00139 }
00140
00141
00142 void EvtLunda::init(){
00143
00144
00145 std::string strpdg=getenv("BESEVTGENROOT");
00146 strpdg +="/share/r_pdg.dat";
00147
00148 std::cout<<"mypdg= "<<strpdg<<std::endl;
00149
00150 if(getNArg()>2){std::cout<<"Parameter can be 1 or 2, You have "<<getNArg()<<std::endl; ::abort();}
00151
00152 if (getParentId().isAlias()){
00153
00154 report(ERROR,"EvtGen") << "EvtLunda finds that you are decaying the"<<endl
00155 << " aliased particle "
00156 << EvtPDL::name(getParentId()).c_str()
00157 << " with the Lunda model"<<endl
00158 << " this does not work, please modify decay table."
00159 << endl;
00160 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00161 ::abort();
00162
00163 }
00164
00165 store(this);
00166
00167 }
00168
00169
00170 std::string EvtLunda::commandName(){
00171
00172 return std::string("LundaPar");
00173
00174 }
00175
00176 void EvtLunda::command(std::string cmd){
00177
00178 if (ncommand==lcommand){
00179
00180 lcommand=10+2*lcommand;
00181
00182 std::string* newcommands=new std::string[lcommand];
00183
00184 int i;
00185
00186 for(i=0;i<ncommand;i++){
00187 newcommands[i]=commands[i];
00188 }
00189
00190 delete [] commands;
00191
00192 commands=newcommands;
00193
00194 }
00195
00196 commands[ncommand]=cmd;
00197
00198 ncommand++;
00199
00200 }
00201
00202
00203
00204 void EvtLunda::decay( EvtParticle *p){
00205
00206 static int iniflag=1;
00207
00208 static EvtId STRNG=EvtPDL::getId("string");
00209
00210 int istdheppar=EvtPDL::getStdHep(p->getId());
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 double mp=p->mass();
00229 float xmp=mp;
00230
00231
00232 EvtVector4R p4[20];
00233
00234 int i,more;
00235 int ip=EvtPDL::getStdHep(p->getId());
00236 int ndaugjs;
00237 static int kf[4000];
00238 EvtId evtnumstable[100],evtnumparton[100];
00239 int stableindex[100],partonindex[100];
00240 int numstable;
00241 int numparton;
00242 static int km[4000];
00243 EvtId type[MAX_DAUG];
00244
00245 int isr=1;
00246 if(getNArg()>0){ isr=getArg(0);}
00247
00248 static double px[4000],py[4000],pz[4000],e[4000];
00249 if (iniflag==1) lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
00250
00251 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
00252
00253 int count=0;
00254 bool mtg=0;
00255
00256 do{
00257
00258 iniflag=iniflag+1;
00259
00260 modeSelection:
00261 lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
00262
00263 mtg = checktag_.decaytag==1;
00264 if(mtg)std::cout<<"checktag_.decaytag="<<checktag_.decaytag<<std::endl;
00265
00266
00267
00268 LundaInit(0);
00269 numstable=0;
00270 numparton=0;
00271
00272 double esum=0;
00273
00274
00275
00276
00277
00278 for(i=0;i<ndaugjs;i++){
00279 if (abs(kf[i])==11 || kf[i]==92 || kf[i]==22) continue;
00280
00281 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
00282 report(ERROR,"EvtGen") << "Lunda returned particle:"<<kf[i]<<endl;
00283 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
00284 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
00285 report(ERROR,"EvtGen") << "The decay was of particle:"<<ip<<endl;
00286
00287 std::cout<<"xmp= "<<xmp<<std::endl;
00288 goto modeSelection;
00289 }
00290
00291
00292 if (abs(kf[i])<=6||kf[i]==21){
00293 partonindex[numparton]=i;
00294 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
00295 numparton++;
00296 }
00297 else{
00298 stableindex[numstable]=i;
00299 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
00300 numstable++;
00301 }
00302
00303 esum+=e[i];
00304
00305
00306
00307
00308 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
00309
00310 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
00311
00312 }
00313
00314 p4[i].set(e[i],px[i],py[i],pz[i]);
00315
00316
00317 }
00318
00319 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
00320
00321
00322
00323 more=(channel!=-1);
00324
00325 if(abs(esum-xmp)>0.001 ){
00326 std::cout<<"Unexpected final states from Lunda with total energy "<<esum<<" unequal to "<<xmp<<std::endl;
00327
00328 std::cout<<"xmp= "<<xmp<<std::endl;
00329 goto modeSelection;
00330 }
00331
00332 count++;
00333 }while( more && (count<10000) && mtg );
00334
00335
00336 if (count>9999) {
00337 report(INFO,"EvtGen") << "Too many loops in EvtLunda!!!"<<endl;
00338 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
00339 for(i=0;i<numstable;i++){
00340 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
00341 }
00342
00343 }
00344
00345
00346
00347 if (numparton==0){
00348
00349 p->makeDaughters(numstable,evtnumstable);
00350 int ndaugFound=0;
00351 for(i=0;i<numstable;i++){
00352 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
00353 ndaugFound++;
00354 }
00355 if ( ndaugFound == 0 ) {
00356 report(ERROR,"EvtGen") << "Lunda has failed to do a decay ";
00357 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
00358 assert(0);
00359 }
00360
00361 fixPolarizations(p);
00362
00363
00364
00365 return ;
00366
00367 }
00368 else{
00369
00370
00371
00372 EvtVector4R p4string(0.0,0.0,0.0,0.0);
00373
00374 for(i=0;i<numparton;i++){
00375 p4string+=p4[partonindex[i]];
00376 }
00377
00378 int nprimary=1;
00379 type[0]=STRNG;
00380 for(i=0;i<numstable;i++){
00381 if (km[stableindex[i]]==0){
00382 type[nprimary++]=evtnumstable[i];
00383 }
00384 }
00385
00386 p->makeDaughters(nprimary,type);
00387
00388 p->getDaug(0)->init(STRNG,p4string);
00389
00390 EvtVector4R p4partons[10];
00391
00392 for(i=0;i<numparton;i++){
00393 p4partons[i]=p4[partonindex[i]];
00394 }
00395
00396 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
00397
00398
00399
00400 nprimary=1;
00401
00402 for(i=0;i<numstable;i++){
00403
00404 if (km[stableindex[i]]==0){
00405 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
00406 }
00407 }
00408
00409
00410 int nsecond=0;
00411 for(i=0;i<numstable;i++){
00412 if (km[stableindex[i]]!=0){
00413 type[nsecond++]=evtnumstable[i];
00414 }
00415 }
00416
00417
00418 p->getDaug(0)->makeDaughters(nsecond,type);
00419
00420 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
00421 -p4string.get(2),-p4string.get(3));
00422
00423 nsecond=0;
00424 for(i=0;i<numstable;i++){
00425 if (km[stableindex[i]]!=0){
00426 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
00427 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
00428 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
00429 p->getDaug(0)->getDaug(nsecond)->decay();
00430 nsecond++;
00431 }
00432 }
00433
00434 if ( nsecond == 0 ) {
00435 report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
00436 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
00437 assert(0);
00438 }
00439
00440 fixPolarizations(p);
00441
00442 return ;
00443
00444 }
00445
00446 }
00447
00448 void EvtLunda::fixPolarizations(EvtParticle *p){
00449
00450
00451
00452 int ndaug=p->getNDaug();
00453
00454 int i;
00455
00456 static EvtId Jpsi=EvtPDL::getId("J/psi");
00457
00458 for(i=0;i<ndaug;i++){
00459 if(p->getDaug(i)->getId()==Jpsi){
00460
00461 EvtSpinDensity rho;
00462
00463 rho.SetDim(3);
00464 rho.Set(0,0,0.5);
00465 rho.Set(0,1,0.0);
00466 rho.Set(0,2,0.0);
00467
00468 rho.Set(1,0,0.0);
00469 rho.Set(1,1,1.0);
00470 rho.Set(1,2,0.0);
00471
00472 rho.Set(2,0,0.0);
00473 rho.Set(2,1,0.0);
00474 rho.Set(2,2,0.5);
00475
00476 EvtVector4R p4Psi=p->getDaug(i)->getP4();
00477
00478 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
00479 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
00480
00481
00482 p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
00483 setDaughterSpinDensity(i);
00484
00485 }
00486 }
00487
00488 }
00489
00490 void EvtLunda::store(EvtDecayBase* jsdecay){
00491
00492 if (nlundadecays==ntable){
00493
00494 EvtDecayBasePtr* newlundadecays=new EvtDecayBasePtr[2*ntable+10];
00495 int i;
00496 for(i=0;i<ntable;i++){
00497 newlundadecays[i]=lundadecays[i];
00498 }
00499 ntable=2*ntable+10;
00500 delete [] lundadecays;
00501 lundadecays=newlundadecays;
00502 }
00503
00504 lundadecays[nlundadecays++]=jsdecay;
00505
00506
00507
00508 }
00509
00510
00511 void EvtLunda::LundaInit(int dummy){
00512
00513 static int first=1;
00514 if (first){
00515
00516 first=0;
00517 for(int i=0;i<ncommand;i++)
00518 lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
00519 }
00520
00521 }
00522
00523 void EvtLunda::ExclusiveDecay(EvtParticle* p){
00524 EvtId daugs[8];
00525 int _ndaugs =4;
00526 daugs[0]=EvtPDL::getId(std::string("pi+"));
00527 daugs[1]=EvtPDL::getId(std::string("pi-"));
00528 daugs[2]=EvtPDL::getId(std::string("pi+"));
00529 daugs[3]=EvtPDL::getId(std::string("pi-"));
00530 checkA:
00531 p->makeDaughters(_ndaugs,daugs);
00532 double totMass=0;
00533 for(int i=0;i< _ndaugs;i++){
00534 EvtParticle* di=p->getDaug(i);
00535 double xmi=EvtPDL::getMass(di->getId());
00536 di->setMass(xmi);
00537 totMass+=xmi;
00538 }
00539 if(totMass>p->mass()) goto checkA;
00540 p->initializePhaseSpace( _ndaugs , daugs);
00541
00542
00543 }