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/EvtTauola.hh"
00029 #include "EvtGenBase/EvtReport.hh"
00030 #include "EvtGenBase/EvtDecayBase.hh"
00031 #include <string>
00032 #include "EvtGenBase/EvtId.hh"
00033 #include <iostream>
00034 #include <iomanip>
00035 #include <fstream>
00036 #include <string.h>
00037 #include <stdlib.h>
00038 #include <unistd.h>
00039 #include <stdio.h>
00040 using std::endl;
00041 using std::fstream;
00042 using std::ios;
00043 using std::ofstream;
00044 using std::resetiosflags;
00045 using std::setiosflags;
00046 using std::setw;
00047
00048
00049 int EvtTauola::ntauoladecays=0;
00050 EvtDecayBasePtr* EvtTauola::tauoladecays=0;
00051 int EvtTauola::ntable=0;
00052
00053 int EvtTauola::ncommand=0;
00054 int EvtTauola::lcommand=0;
00055 std::string* EvtTauola::commands=0;
00056
00057
00058 extern "C" {
00059 extern void dectes_(int *, int *,int *,int *,int *,
00060 double *,double *,double *,double *);
00061 }
00062
00063
00064 EvtTauola::EvtTauola(){}
00065
00066 EvtTauola::~EvtTauola(){
00067
00068
00069 int i;
00070
00071
00072
00073
00074 if (ntauoladecays==0) {
00075 delete [] commands;
00076 commands=0;
00077 return;
00078 }
00079
00080 for(i=0;i<ntauoladecays;i++){
00081 if (tauoladecays[i]==this){
00082 tauoladecays[i]=tauoladecays[ntauoladecays-1];
00083 ntauoladecays--;
00084 if (ntauoladecays==0) {
00085 delete [] commands;
00086 commands=0;
00087 }
00088 return;
00089 }
00090 }
00091
00092 report(ERROR,"EvtGen") << "Error in destroying Tauola model!"<<endl;
00093
00094 }
00095
00096
00097 void EvtTauola::getName(std::string& model_name){
00098
00099 model_name="TAUOLA";
00100
00101 }
00102
00103 EvtDecayBase* EvtTauola::clone(){
00104
00105 return new EvtTauola;
00106
00107 }
00108
00109
00110 void EvtTauola::initProbMax(){
00111
00112 noProbMax();
00113
00114 }
00115
00116
00117 void EvtTauola::init(){
00118
00119
00120
00121
00122 if (getParentId().isAlias()){
00123
00124 report(ERROR,"EvtGen") << "EvtTauola finds that you are decaying the"<<endl
00125 << " aliased particle "
00126 << EvtPDL::name(getParentId()).c_str()
00127 << " with the Tauola model"<<endl
00128 << " this does not work, please modify decay table."
00129 << endl;
00130 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00131 ::abort();
00132
00133 }
00134
00135 store(this);
00136
00137
00138 }
00139
00140
00141 std::string EvtTauola::commandName(){
00142
00143 return std::string("TauolaPar");
00144
00145 }
00146
00147 void EvtTauola::command(std::string cmd){
00148
00149 if (ncommand==lcommand){
00150
00151 lcommand=10+2*lcommand;
00152
00153 std::string* newcommands=new std::string[lcommand];
00154
00155 int i;
00156
00157 for(i=0;i<ncommand;i++){
00158 newcommands[i]=commands[i];
00159 }
00160
00161 delete [] commands;
00162
00163 commands=newcommands;
00164
00165 }
00166
00167 commands[ncommand]=cmd;
00168
00169 ncommand++;
00170
00171 }
00172
00173
00174
00175 void EvtTauola::decay( EvtParticle *p){
00176
00177 static int iniflag=0;
00178
00179 static EvtId STRNG=EvtPDL::getId("string");
00180
00181 int istdheppar=EvtPDL::getStdHep(p->getId());
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 EvtVector4R p4[20];
00192
00193
00194 int i,more;
00195 int idtau=EvtPDL::getStdHep(p->getId());
00196 int ndaugjs;
00197 static int kf[20];
00198 EvtId evtnumstable[20],evtnumparton[20];
00199 int stableindex[20],partonindex[20];
00200 int numstable;
00201 int numparton;
00202 static int km[20];
00203 EvtId type[MAX_DAUG];
00204
00205 static double px[20],py[20],pz[20],e[20];
00206
00207
00208 if (iniflag==0) dectes_(&iniflag,&idtau,&ndaugjs,kf,km,px,py,pz,e);
00209
00210 if ( p->getNDaug() != 0 ) { p->deleteDaughters(true);}
00211
00212 int count=0;
00213
00214 do{
00215
00216 iniflag=iniflag+1;
00217 dectes_(&iniflag,&idtau,&ndaugjs,kf,km,px,py,pz,e);
00218
00219 numstable=0;
00220 numparton=0;
00221
00222 for(i=0;i<ndaugjs;i++){
00223
00224
00225 if (EvtPDL::evtIdFromStdHep(kf[i])==EvtId(-1,-1)) {
00226 report(ERROR,"EvtGen") << "Tauola returned particle:"<<kf[i]<<endl;
00227 report(ERROR,"EvtGen") << "This can not be translated to evt number"<<endl;
00228 report(ERROR,"EvtGen") << "and the decay will be rejected!"<<endl;
00229 report(ERROR,"EvtGen") << "The decay was of particle:"<<idtau<<endl;
00230
00231 }
00232
00233
00234 if (abs(kf[i])<=6||kf[i]==21){
00235 partonindex[numparton]=i;
00236 evtnumparton[numparton]=EvtPDL::evtIdFromStdHep(kf[i]);
00237 numparton++;
00238 }
00239 else{
00240 stableindex[numstable]=i;
00241 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(kf[i]);
00242 numstable++;
00243 }
00244
00245
00246
00247
00248
00249
00250 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
00251
00252 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
00253
00254 }
00255
00256 p4[i].set(e[i],px[i],py[i],pz[i]);
00257
00258
00259 }
00260
00261 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
00262
00263
00264 more=(channel!=-1);
00265
00266
00267 count++;
00268
00269 }while( more && (count<10000) );
00270
00271 if (count>9999) {
00272 report(INFO,"EvtGen") << "Too many loops in EvtTauola!!!"<<endl;
00273 report(INFO,"EvtGen") << "Parent:"<<EvtPDL::name(getParentId()).c_str()<<endl;
00274 for(i=0;i<numstable;i++){
00275 report(INFO,"EvtGen") << "Daug("<<i<<")"<<EvtPDL::name(evtnumstable[i]).c_str()<<endl;
00276 }
00277
00278 }
00279
00280
00281
00282 if (numparton==0){
00283
00284 p->makeDaughters(numstable,evtnumstable);
00285 int ndaugFound=0;
00286 for(i=0;i<numstable;i++){
00287 p->getDaug(i)->init(evtnumstable[i],p4[stableindex[i]]);
00288 ndaugFound++;
00289 }
00290 if ( ndaugFound == 0 ) {
00291 report(ERROR,"EvtGen") << "Tauola has failed to do a decay ";
00292 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
00293 assert(0);
00294 }
00295
00296 fixPolarizations(p);
00297
00298 return ;
00299
00300 }
00301 else{
00302
00303
00304
00305 EvtVector4R p4string(0.0,0.0,0.0,0.0);
00306
00307 for(i=0;i<numparton;i++){
00308 p4string+=p4[partonindex[i]];
00309 }
00310
00311 int nprimary=1;
00312 type[0]=STRNG;
00313 for(i=0;i<numstable;i++){
00314 if (km[stableindex[i]]==0){
00315 type[nprimary++]=evtnumstable[i];
00316 }
00317 }
00318
00319 p->makeDaughters(nprimary,type);
00320
00321 p->getDaug(0)->init(STRNG,p4string);
00322
00323 EvtVector4R p4partons[10];
00324
00325 for(i=0;i<numparton;i++){
00326 p4partons[i]=p4[partonindex[i]];
00327 }
00328
00329 ((EvtStringParticle*)p->getDaug(0))->initPartons(numparton,p4partons,evtnumparton);
00330
00331
00332
00333 nprimary=1;
00334
00335 for(i=0;i<numstable;i++){
00336
00337 if (km[stableindex[i]]==0){
00338 p->getDaug(nprimary++)->init(evtnumstable[i],p4[stableindex[i]]);
00339 }
00340 }
00341
00342
00343 int nsecond=0;
00344 for(i=0;i<numstable;i++){
00345 if (km[stableindex[i]]!=0){
00346 type[nsecond++]=evtnumstable[i];
00347 }
00348 }
00349
00350
00351 p->getDaug(0)->makeDaughters(nsecond,type);
00352
00353 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
00354 -p4string.get(2),-p4string.get(3));
00355
00356 nsecond=0;
00357 for(i=0;i<numstable;i++){
00358 if (km[stableindex[i]]!=0){
00359 p4[stableindex[i]]=boostTo(p4[stableindex[i]],p4stringboost);
00360 p->getDaug(0)->getDaug(nsecond)->init(evtnumstable[i],p4[stableindex[i]]);
00361 p->getDaug(0)->getDaug(nsecond)->setDiagonalSpinDensity();
00362 p->getDaug(0)->getDaug(nsecond)->decay();
00363 nsecond++;
00364 }
00365 }
00366
00367 if ( nsecond == 0 ) {
00368 report(ERROR,"EvtGen") << "Jetset has failed to do a decay ";
00369 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass() <<endl;
00370 assert(0);
00371 }
00372
00373 fixPolarizations(p);
00374
00375 return ;
00376
00377 }
00378
00379 }
00380
00381 void EvtTauola::fixPolarizations(EvtParticle *p){
00382
00383
00384
00385 int ndaug=p->getNDaug();
00386
00387 int i;
00388
00389
00390 static EvtId tau=getParentId();
00391
00392 for(i=0;i<ndaug;i++){
00393 if(p->getDaug(i)->getId()==tau){
00394
00395 EvtSpinDensity rho;
00396
00397 rho.SetDim(2);
00398 rho.Set(0,0,1.0);
00399 rho.Set(0,1,0.0);
00400
00401 rho.Set(1,0,0.0);
00402 rho.Set(1,1,1.0);
00403
00404
00405 EvtVector4R p4Psi=p->getDaug(i)->getP4();
00406
00407 double alpha=atan2(p4Psi.get(2),p4Psi.get(1));
00408 double beta=acos(p4Psi.get(3)/p4Psi.d3mag());
00409
00410
00411 p->getDaug(i)->setSpinDensityForwardHelicityBasis(rho,alpha,beta,0.0);
00412 setDaughterSpinDensity(i);
00413
00414 }
00415 }
00416
00417 }
00418
00419 void EvtTauola::store(EvtDecayBase* jsdecay){
00420
00421 if (ntauoladecays==ntable){
00422
00423 EvtDecayBasePtr* newtauoladecays=new EvtDecayBasePtr[2*ntable+10];
00424 int i;
00425 for(i=0;i<ntable;i++){
00426 newtauoladecays[i]=tauoladecays[i];
00427 }
00428 ntable=2*ntable+10;
00429 delete [] tauoladecays;
00430 tauoladecays=newtauoladecays;
00431 }
00432
00433 tauoladecays[ntauoladecays++]=jsdecay;
00434
00435
00436
00437 }
00438
00439