TEvtGen/Tauola/Tauola.cxx

   1 #include <fstream>
   2 #include <cstring>
   3 #include <vector>
   4 #include "TauolaLog.h"
   5 #include "Tauola.h"
   6 #include "TauolaEvent.h"
   7
   8 namespace Tauolapp
   9 {
  10
  11 int    Tauola::m_pdg_id          = 15;
  12 int    Tauola::m_firstDecayMode  = 0;
  13 int    Tauola::m_secondDecayMode = 0;
  14 bool   Tauola::m_rad             = true;
  15 double Tauola::m_rad_cut_off     = 0.001;
  16 double Tauola::m_iniphy          = 0.1;
  17 double Tauola::m_higgs_scalar_pseudoscalar_mix = M_PI/4;
  18 int    Tauola::m_higgs_scalar_pseudoscalar_pdg = 35;
  19 int    Tauola::m_helPlus  = 0;
  20 int    Tauola::m_helMinus = 0;
  21 double Tauola::m_wtEW     = 0.0;
  22 double Tauola::m_wtEW0    = 0.0;
  23 double Tauola::table11A[NS1][NCOS][4][4] = {{{{0.0}}}};
  24 double Tauola::table1A [NS1][NCOS][4][4] = {{{{0.0}}}};
  25 double Tauola::table2A [NS1][NCOS][4][4] = {{{{0.0}}}};
  26 double Tauola::wtable11A[NS1][NCOS]      = {{0.0}};
  27 double Tauola::wtable1A [NS1][NCOS]      = {{0.0}};
  28 double Tauola::wtable2A [NS1][NCOS]      = {{0.0}};
  29 double Tauola::w0table11A[NS1][NCOS]     = {{0.0}};
  30 double Tauola::w0table1A [NS1][NCOS]     = {{0.0}};
  31 double Tauola::w0table2A [NS1][NCOS]     = {{0.0}};
  32
  33 double Tauola::table11B[NS2][NCOS][4][4] = {{{{0.0}}}};
  34 double Tauola::table1B [NS2][NCOS][4][4] = {{{{0.0}}}};
  35 double Tauola::table2B [NS2][NCOS][4][4] = {{{{0.0}}}};
  36 double Tauola::wtable11B[NS2][NCOS]      = {{0.0}};
  37 double Tauola::wtable1B [NS2][NCOS]      = {{0.0}};
  38 double Tauola::wtable2B [NS2][NCOS]      = {{0.0}};
  39 double Tauola::w0table11B[NS2][NCOS]     = {{0.0}};
  40 double Tauola::w0table1B [NS2][NCOS]     = {{0.0}};
  41 double Tauola::w0table2B [NS2][NCOS]     = {{0.0}};
  42
  43 double Tauola::table11C[NS3][NCOS][4][4] = {{{{0.0}}}};
  44 double Tauola::table1C [NS3][NCOS][4][4] = {{{{0.0}}}};
  45 double Tauola::table2C [NS3][NCOS][4][4] = {{{{0.0}}}};
  46 double Tauola::wtable11C[NS3][NCOS]      = {{0.0}};
  47 double Tauola::wtable1C [NS3][NCOS]      = {{0.0}};
  48 double Tauola::wtable2C [NS3][NCOS]      = {{0.0}};
  49 double Tauola::w0table11C[NS3][NCOS]     = {{0.0}};
  50 double Tauola::w0table1C [NS3][NCOS]     = {{0.0}};
  51 double Tauola::w0table2C [NS3][NCOS]     = {{0.0}};
  52
  53 double Tauola::sminA = 0;
  54 double Tauola::smaxA = 0;
  55
  56 double Tauola::sminB = 0;
  57 double Tauola::smaxB = 0;
  58
  59 double Tauola::sminC = 0;
  60 double Tauola::smaxC = 0;
  61
  62 int    Tauola::ion[3] = {0};
  63 double Tauola::tau_lifetime = .08711;
  64 double Tauola::momentum_conservation_threshold   = 0.1;
  65
  66 Tauola::Particles Tauola::spin_correlation;
  67
  68 Tauola::MomentumUnits Tauola::momentumUnit = Tauola::DEFAULT_MOMENTUM;
  69 Tauola::LengthUnits   Tauola::lengthUnit   = Tauola::DEFAULT_LENGTH;
  70
  71 bool   Tauola::m_is_using_decay_one        = false;
  72 double Tauola::m_decay_one_polarization[3] = {0};
  73 void (*Tauola::m_decay_one_boost_routine)(TauolaParticle*,TauolaParticle*) = NULL;
  74
  75 int     Tauola::buf_incoming_pdg_id = 0;
  76 int     Tauola::buf_outgoing_pdg_id = 0;
  77 double  Tauola::buf_invariant_mass_squared = -1.;
  78 double  Tauola::buf_cosTheta               = 0.;
  79
  80 double  Tauola::buf_R[4][4] = {{0.0}}; //density matrix
  81
  82 double (*Tauola::randomDouble)()                               = Tauola::defaultRandomGenerator;
  83 void   (*Tauola::redefineTauPlusProperties)(TauolaParticle *)  = defaultRedPlus;
  84 void   (*Tauola::redefineTauMinusProperties)(TauolaParticle *) = defaultRedMinus;
  85
  86 /**************************************************************/
  87 void Tauola::setNewCurrents(int mode)
  88 {
  89   inirchl_(&mode);
  90 }
  91
  92 double Tauola::defaultRandomGenerator(){
  93   return rand()*1./RAND_MAX;
  94 }
  95
  96 void Tauola::setRandomGenerator(double (*gen)()){
  97   if(gen==NULL) randomDouble = defaultRandomGenerator;
  98   else          randomDouble = gen;
  99 }
 100
 101 void Tauola::defaultRedPlus(TauolaParticle *tau)  {}
 102 void Tauola::defaultRedMinus(TauolaParticle *tau) {}
 103
 104 void Tauola::setRedefineTauMinus( void (*fun)(TauolaParticle *) ){
 105   redefineTauMinusProperties=fun;
 106 }
 107
 108 void Tauola::setRedefineTauPlus ( void (*fun)(TauolaParticle *) ){
 109   redefineTauPlusProperties=fun;
 110 }
 111
 112 void Tauola::getBornKinematics(int *incoming_pdg_id, int *outgoing_pdg_id, double *invariant_mass_squared,double *cosTheta){
 113   *incoming_pdg_id = buf_incoming_pdg_id;
 114   *outgoing_pdg_id = buf_outgoing_pdg_id;
 115   *invariant_mass_squared = buf_invariant_mass_squared;
 116   *cosTheta               = buf_cosTheta;
 117   //  m_R[0][0] to be added in next step;
 118 }
 119
 120 void Tauola::setUnits(MomentumUnits m, LengthUnits l){
 121   Tauola::momentumUnit = m;
 122   Tauola::lengthUnit   = l;
 123 }
 124
 125 void Tauola::setTauLifetime(double t){
 126   tau_lifetime = t;
 127 }
 128
 129 void Tauola::initialize(){
 130   printf("\n");
 131   printf(" *************************************\n");
 132   printf(" *     TAUOLA C++ Interface v1.1.4   *\n");
 133   printf(" *-----------------------------------*\n");
 134   printf(" *                                   *\n");
 135   printf(" *   (c) Nadia    Davidson,   (1,2)  *\n");
 136   printf(" *       Gizo     Nanava,     (3)    *\n");
 137   printf(" *       Tomasz   Przedzinski,(4)    *\n");
 138   printf(" *       Elzbieta Richter-Was,(2,4)  *\n");
 139   printf(" *       Zbigniew Was         (2,5)  *\n");
 140   printf(" *                                   *\n");
 141   printf(" *  1) Unimelb, Melbourne, Australia *\n");
 142   printf(" *  2)     INP, Krakow, Poland       *\n");
 143   printf(" *  3) University Bonn, Germany      *\n");
 144   printf(" *  4)      UJ, Krakow, Poland       *\n");
 145   printf(" *  5)    CERN, Geneva, Switzerland  *\n");
 146   printf(" *************************************\n");
 147
 148   // Turn on all spin correlations
 149   spin_correlation.setAll(true);
 150
 151   // Ininitalize tauola-fortran
 152   f_interface_tauolaInitialize(m_pdg_id,m_firstDecayMode,
 153                                m_secondDecayMode,m_rad,
 154                                m_rad_cut_off, m_iniphy);
 155
 156   //---------------------------------------------------------------------------
 157   // Initialize SANC tables
 158   //---------------------------------------------------------------------------
 159   cout<<"Reading SANC input files."<<endl;
 160
 161   ifstream f("table1-1.txt");
 162
 163   if(!f.is_open()){
 164     cout<<"File 'table1-1.txt'  missing... skipped."<<endl;
 165   }
 166   else{
 167     string buf;
 168
 169     cout<<"Reading file 'table1-1.txt'..."<<endl;
 170
 171     int dbuf1,dbuf2,dbuf3,dbufcos;
 172     f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
 173
 174     // Check table sizes
 175     if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS)  {
 176       cout<<"mismatched NS1=   "<<dbuf1  <<" <--> "<<NS1<<endl;
 177       cout<<"           NS2=   "<<dbuf2  <<" <--> "<<NS2<<endl;
 178       cout<<"           NS3=   "<<dbuf3  <<" <--> "<<NS3<<endl;
 179       cout<<"           NCOS=  "<<dbufcos<<" <--> "<<NCOS<<endl;
 180       return;
 181     }
 182
 183     double buf1,buf2,buf3,buf4,buf5,buf6;
 184     f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
 185
 186     // Set ranges
 187     if(sminA==0.0){
 188       sminA=buf1;
 189       smaxA=buf2;
 190       sminB=buf3;
 191       smaxB=buf4;
 192       sminC=buf5;
 193       smaxC=buf6;
 194     }
 195
 196     // Check ranges
 197     if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
 198       cout<<"mismatched sminA=   "<<buf1<<" <--> "<<sminA<<endl;
 199       cout<<"           smaxA=   "<<buf2<<" <--> "<<smaxA<<endl;
 200       cout<<"           sminB=   "<<buf3<<" <--> "<<sminB<<endl;
 201       cout<<"           smaxB=   "<<buf4<<" <--> "<<smaxB<<endl;
 202       cout<<"           sminC=   "<<buf5<<" <--> "<<sminC<<endl;
 203       cout<<"           smaxC=   "<<buf6<<" <--> "<<smaxC<<endl;
 204       return;
 205     }
 206
 207     // Print out header
 208     while(!f.eof()){
 209       char head[255];
 210       f.getline(head,255);
 211       if(strcmp(head,"BeginRange1")==0) break;
 212       cout<<head<<endl;
 213     }
 214
 215     // Read table
 216     for (int i=0;i<NS1;i++){
 217       for (int j=0;j<NCOS;j++){
 218         for (int k=0;k<4;k++){
 219           for (int l=0;l<4;l++){
 220             f>>table1A[i][j][k][l];
 221           } // for(l)
 222         } // for(k)
 223         f>>wtable1A[i][j];
 224         f>>w0table1A[i][j];
 225       } // for(j)
 226     } // for(i)
 227
 228     // Find 2nd range
 229     while(!f.eof()){
 230       f>>buf;
 231       if(strcmp(buf.c_str(),"BeginRange2")==0) break;
 232     }
 233
 234     // Read table
 235     for (int i=0;i<NS2;i++){
 236       for (int j=0;j<NCOS;j++){
 237         for (int k=0;k<4;k++){
 238           for (int l=0;l<4;l++){
 239             f>>table1B[i][j][k][l];
 240           } // for(l)
 241         } // for(k)
 242         f>>wtable1B[i][j];
 243         f>>w0table1B[i][j];
 244       } // for(j)
 245     } // for(i)
 246
 247     // Find 3rd range
 248     while(!f.eof()){
 249       f>>buf;
 250       if(strcmp(buf.c_str(),"BeginRange3")==0) break;
 251     }
 252
 253     // Read table
 254     for (int i=0;i<NS3;i++){
 255       for (int j=0;j<NCOS;j++){
 256         for (int k=0;k<4;k++){
 257           for (int l=0;l<4;l++){
 258             f>>table1C[i][j][k][l];
 259           } // for(l)
 260         } // for(k)
 261         f>>wtable1C[i][j];
 262         f>>w0table1C[i][j];
 263       } // for(j)
 264     } // for(i)
 265
 266     // Check for proper file end
 267     f>>buf;
 268     if(buf.size() == 0 || strcmp(buf.c_str(),"End") != 0){
 269       cout<<"...incorrect file version or file incomplete/damaged!"<<endl;
 270
 271       // In case of the error - do not use tables
 272       table1A[0][0][0][0] = table1B[0][0][0][0] = table1C[0][0][0][0] = 0.0;
 273     }
 274   }  // if (file is open)
 275
 276   f.close();
 277   f.open("table2-2.txt");
 278
 279   if(!f.is_open()){
 280     cout<<"File 'table2-2.txt'  missing... skipped."<<endl;
 281   }
 282   else{
 283     string buf;
 284
 285     cout<<"Reading file 'table2-2.txt'..."<<endl;
 286
 287     int dbuf1,dbuf2,dbuf3,dbufcos;
 288     f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
 289
 290     // Check table sizes
 291     if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS)  {
 292       cout<<"mismatched NS1=   "<<dbuf1<<" <--> "<<NS1<<endl;
 293       cout<<"           NS2=   "<<dbuf2<<" <--> "<<NS2<<endl;
 294       cout<<"           NS3=   "<<dbuf3<<" <--> "<<NS3<<endl;
 295       cout<<"           NCOS=  "<<dbufcos<<" <--> "<<NCOS<<endl;
 296       return;
 297     }
 298
 299     double buf1,buf2,buf3,buf4,buf5,buf6;
 300     f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
 301
 302     // Set ranges
 303     if(sminA==0.0)
 304     {
 305       sminA=buf1;
 306       smaxA=buf2;
 307       sminB=buf3;
 308       smaxB=buf4;
 309       sminC=buf5;
 310       smaxC=buf6;
 311     }
 312
 313     // Check ranges
 314     if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
 315       cout<<"mismatched sminA=   "<<buf1<<" <--> "<<sminA<<endl;
 316       cout<<"           smaxA=   "<<buf2<<" <--> "<<smaxA<<endl;
 317       cout<<"           sminB=   "<<buf3<<" <--> "<<sminB<<endl;
 318       cout<<"           smaxB=   "<<buf4<<" <--> "<<smaxB<<endl;
 319       cout<<"           sminC=   "<<buf5<<" <--> "<<sminC<<endl;
 320       cout<<"           smaxC=   "<<buf6<<" <--> "<<smaxC<<endl;
 321       return;
 322     }
 323
 324     // Print out header
 325     while(!f.eof()){
 326       char head[255];
 327       f.getline(head,255);
 328       if(strcmp(head,"BeginRange1")==0) break;
 329       cout<<head<<endl;
 330     }
 331
 332     // Read table
 333     for (int i=0;i<NS1;i++){
 334       for (int j=0;j<NCOS;j++){
 335         for (int k=0;k<4;k++){
 336           for (int l=0;l<4;l++){
 337             f>>table2A[i][j][k][l];
 338           } // for(l)
 339         } // for(k)
 340         f>>wtable2A[i][j];
 341         f>>w0table2A[i][j];
 342       } // for(j)
 343     } // for(i)
 344
 345     // Find 2nd range
 346     while(!f.eof()){
 347       f>>buf;
 348       if(strcmp(buf.c_str(),"BeginRange2")==0) break;
 349     }
 350
 351     // Read table
 352     for (int i=0;i<NS2;i++){
 353       for (int j=0;j<NCOS;j++){
 354         for (int k=0;k<4;k++){
 355           for (int l=0;l<4;l++){
 356             f>>table2B[i][j][k][l];
 357           } // for(l)
 358         } // for(k)
 359         f>>wtable2B[i][j];
 360         f>>w0table2B[i][j];
 361       } // for(j)
 362     } // for(i)
 363
 364     // Find 3rd range
 365     while(!f.eof()){
 366       f>>buf;
 367       if(strcmp(buf.c_str(),"BeginRange3")==0) break;
 368     }
 369
 370     // Read table
 371     for (int i=0;i<NS3;i++){
 372       for (int j=0;j<NCOS;j++){
 373         for (int k=0;k<4;k++){
 374           for (int l=0;l<4;l++){
 375             f>>table2C[i][j][k][l];
 376           } // for(l)
 377         } // for(k)
 378         f>>wtable2C[i][j];
 379         f>>w0table2C[i][j];
 380       } // for(j)
 381     } // for(i)
 382
 383     // Check for proper file end
 384     f>>buf;
 385     if(buf.size()==0 || strcmp(buf.c_str(),"End")!=0){
 386       cout<<"...incorrect file version or file incomplete/damaged!"<<endl;
 387
 388       // In case of the error - do not use tables
 389       table2A[0][0][0][0] = table2B[0][0][0][0] = table2C[0][0][0][0] = 0.0;
 390     }
 391   } // if (file is open)
 392
 393   f.close();
 394   f.open("table11-11.txt");
 395
 396   if(!f.is_open()){
 397     cout<<"File 'table11-11.txt' missing... skipped."<<endl;
 398   }
 399   else{
 400     string buf;
 401
 402     cout<<"Reading file 'table11-11.txt'..."<<endl;
 403
 404     int dbuf1,dbuf2,dbuf3,dbufcos;
 405     f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
 406
 407     // Check table sizes
 408     if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS)  {
 409       cout<<"mismatched NS1=   "<<dbuf1<<" <--> "<<NS1<<endl;
 410       cout<<"           NS2=   "<<dbuf2<<" <--> "<<NS2<<endl;
 411       cout<<"           NS3=   "<<dbuf3<<" <--> "<<NS3<<endl;
 412       cout<<"           NCOS=  "<<dbufcos<<" <--> "<<NCOS<<endl;
 413       return;
 414     }
 415
 416     double buf1,buf2,buf3,buf4,buf5,buf6;
 417     f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
 418
 419     // Set ranges
 420     if(sminA==0.0)
 421     {
 422       sminA=buf1;
 423       smaxA=buf2;
 424       sminB=buf3;
 425       smaxB=buf4;
 426       sminC=buf5;
 427       smaxC=buf6;
 428     }
 429
 430     // Check ranges
 431     if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
 432       cout<<"mismatched sminA=   "<<buf1<<" <--> "<<sminA<<endl;
 433       cout<<"           smaxA=   "<<buf2<<" <--> "<<smaxA<<endl;
 434       cout<<"           sminB=   "<<buf3<<" <--> "<<sminB<<endl;
 435       cout<<"           smaxB=   "<<buf4<<" <--> "<<smaxB<<endl;
 436       cout<<"           sminC=   "<<buf5<<" <--> "<<sminC<<endl;
 437       cout<<"           smaxC=   "<<buf6<<" <--> "<<smaxC<<endl;
 438       return;
 439     }
 440
 441     // Print out header
 442     while(!f.eof()){
 443       char head[255];
 444       f.getline(head,255);
 445       if(strcmp(head,"BeginRange1")==0) break;
 446       cout<<head<<endl;
 447     }
 448
 449     // Read table
 450     for (int i=0;i<NS1;i++){
 451       for (int j=0;j<NCOS;j++){
 452         for (int k=0;k<4;k++){
 453           for (int l=0;l<4;l++){
 454             f>>table11A[i][j][k][l];
 455           } // for(l)
 456         } // for(k)
 457         f>>wtable11A[i][j];
 458         f>>w0table11A[i][j];
 459       } // for(j)
 460     } // for(i)
 461
 462     // Find 2nd range
 463     while(!f.eof()){
 464       f>>buf;
 465       if(strcmp(buf.c_str(),"BeginRange2")==0) break;
 466     }
 467
 468     // Read table
 469     for (int i=0;i<NS2;i++){
 470       for (int j=0;j<NCOS;j++){
 471         for (int k=0;k<4;k++){
 472           for (int l=0;l<4;l++){
 473             f>>table11B[i][j][k][l];
 474           } // for(l)
 475         } // for(k)
 476         f>>wtable11B[i][j];
 477         f>>w0table11B[i][j];
 478       } // for(j)
 479     } // for(i)
 480
 481     // Find 3rd range
 482     while(!f.eof()){
 483       f>>buf;
 484       if(strcmp(buf.c_str(),"BeginRange3")==0) break;
 485     }
 486
 487     // Read table
 488     for (int i=0;i<NS3;i++){
 489       for (int j=0;j<NCOS;j++){
 490         for (int k=0;k<4;k++){
 491           for (int l=0;l<4;l++){
 492             f>>table11C[i][j][k][l];
 493           } // for(l)
 494         } // for(k)
 495         f>>wtable11C[i][j];
 496         f>>w0table11C[i][j];
 497       } // for(j)
 498     } // for(i)
 499
 500     f>>buf;
 501     if(buf.size()==0 || strcmp(buf.c_str(),"End")!=0){
 502       cout<<"...incorrect file version or file incomplete/damaged!"<<endl;
 503
 504       // In case of the error - do not use tables
 505       table11A[0][0][0][0] = table11B[0][0][0][0] = table11C[0][0][0][0] = 0.0;
 506     }
 507   } // if (file is open)
 508
 509   f.close();
 510   cout<<endl;
 511 }
 512
 513 void Tauola::decayOne(TauolaParticle *tau, bool undecay, double polx, double poly, double polz)
 514 {
 515   if(!tau) return;
 516
 517   if(polx*polx+poly*poly+polz*polz>1)
 518   {
 519     Log::Warning()<<"decayOne(): ignoring wrong polarization vector: "<<polx<<" "<<poly<<" "<<polz<<endl;
 520     polx=poly=polz=0;
 521   }
 522
 523   // Let the interface know that we work in the decayOne mode
 524   m_is_using_decay_one = true;
 525
 526   m_decay_one_polarization[0] = polx;
 527   m_decay_one_polarization[1] = poly;
 528   m_decay_one_polarization[2] = polz;
 529
 530   // Undecay if needed
 531   if(tau->hasDaughters())
 532   {
 533     if(undecay) tau->undecay();
 534     else
 535     {
 536       m_is_using_decay_one = false;
 537       return;
 538     }
 539   }
 540
 541   std::vector<TauolaParticle *> list;
 542   list.push_back(tau);
 543
 544   // Decay single tau
 545   TauolaParticlePair t_pair(list);
 546   t_pair.decayTauPair();
 547   t_pair.checkMomentumConservation();
 548
 549   // Revert to normal mode
 550   m_is_using_decay_one = false;
 551 }
 552
 553 void Tauola::initialise(){
 554
 555   Log::Warning() <<"Deprecated routine 'Tauola::initialise'"<<endl;
 556   Log::Warning(0)<<"Use 'Tauola::initialize' instead."<<endl;
 557
 558   initialize();
 559
 560   // Deprecated routines:  initialise, setInitialisePhy,
 561   //                       f_interface_tauolaInitialise
 562 }
 563
 564 bool Tauola::isUsingDecayOne()
 565 {
 566   return m_is_using_decay_one;
 567 }
 568
 569 bool Tauola::isUsingDecayOneBoost()
 570 {
 571   return (bool) m_decay_one_boost_routine;
 572 }
 573
 574 void Tauola::setBoostRoutine(void (*boost)(TauolaParticle *, TauolaParticle *))
 575 {
 576   m_decay_one_boost_routine=boost;
 577 }
 578
 579 void Tauola::decayOneBoost(TauolaParticle *mother, TauolaParticle *target)
 580 {
 581   m_decay_one_boost_routine(mother,target);
 582 }
 583
 584 const double* Tauola::getDecayOnePolarization()
 585 {
 586   return m_decay_one_polarization;
 587 }
 588
 589 void Tauola::setDecayingParticle(int pdg_id){
 590   m_pdg_id=pdg_id;
 591 }
 592
 593 int Tauola::getDecayingParticle(){
 594   return abs(m_pdg_id);
 595 }
 596
 597 void Tauola::setSameParticleDecayMode(int firstDecayMode){
 598   m_firstDecayMode=firstDecayMode;
 599 }
 600
 601 void Tauola::setOppositeParticleDecayMode(int secondDecayMode){
 602   m_secondDecayMode=secondDecayMode;
 603 }
 604
 605 void Tauola::setRadiation(bool rad){
 606   m_rad=rad;
 607 }
 608
 609 void Tauola::setRadiationCutOff(double rad_cut_off){
 610   m_rad_cut_off=rad_cut_off;
 611 }
 612
 613
 614 void Tauola::setInitializePhy(double iniphy){
 615   m_iniphy=iniphy;
 616 }
 617
 618 void Tauola::setInitialisePhy(double iniphy){
 619
 620   Log::Warning() <<"Deprecated routine 'Tauola::setInitialisePhy'"<<endl;
 621   Log::Warning(0)<<"Use 'Tauola::setInitializePhy' instead."<<endl;
 622
 623   setInitializePhy(iniphy);
 624
 625   // Deprecated routines:  initialise, setInitialisePhy,
 626   //                       f_interface_tauolaInitialise
 627 }
 628
 629 void Tauola::setTauBr(int i, double value)
 630 {
 631   if(taubra_.nchan==0)
 632     Log::Warning()<<"setTauBr(): run Tauola::initialize() first."<<endl;
 633   else if(i<1 || i>taubra_.nchan || value<0.)
 634     Log::Warning()<<"setTauBr(): Invalid input. Value must be >= 0 and 0 < i <= "<<taubra_.nchan<<endl;
 635   else taubra_.gamprt[i-1]=(float)value;
 636 }
 637
 638 void Tauola::setTaukle(double bra1,double brk0, double brk0b, double brks)
 639 {
 640   if(bra1<0 || bra1>1 || brk0<0 ||brk0>1 || brk0b<0 || brk0b>1 || brks<0 ||brks>1)
 641   {
 642     Log::Warning()<<"setTaukle(): variables must be in range [0,1]. Ignored."<<endl;
 643     return;
 644   }
 645   taukle_.bra1 =(float)bra1;
 646   taukle_.brk0 =(float)brk0;
 647   taukle_.brk0b=(float)brk0b;
 648   taukle_.brks =(float)brks;
 649 }
 650
 651 double Tauola::getTauMass(){
 652   return f_getTauMass();
 653 }
 654
 655 double Tauola::getHiggsScalarPseudoscalarMixingAngle(){
 656   return m_higgs_scalar_pseudoscalar_mix;
 657 }
 658
 659 int Tauola::getHiggsScalarPseudoscalarPDG(){
 660   return m_higgs_scalar_pseudoscalar_pdg;
 661 }
 662
 663 /** set the mixing angle. coupling: tau~(cos(phi)+isin(phi)gamma5)tau */
 664 void Tauola::setHiggsScalarPseudoscalarMixingAngle(double angle){
 665   m_higgs_scalar_pseudoscalar_mix=angle;
 666 }
 667
 668 /** set the pdg code of the higgs particle which tauola should
 669     treat as a scalar-pseudoscalar mix  */
 670 void Tauola::setHiggsScalarPseudoscalarPDG(int pdg_code){
 671
 672   if (particleCharge(pdg_code)!=0.0){
 673     Log::Warning()<<"You want to use spin correlations of Higgs for particle of PDGID= "<<pdg_code<<endl
 674                   <<"This particle has charge="<<particleCharge(pdg_code)<<endl;
 675     Log::Fatal("This choice is not appropriate.",0);
 676   }
 677   m_higgs_scalar_pseudoscalar_pdg=pdg_code;
 678 }
 679
 680 int Tauola::getHelPlus(){
 681   return m_helPlus;
 682 }
 683 int Tauola::getHelMinus(){
 684   return m_helMinus;
 685 }
 686 double Tauola::getEWwt(){
 687   return m_wtEW;
 688 }
 689 double Tauola::getEWwt0(){
 690   return m_wtEW0;
 691 }
 692 void Tauola::setEWwt(double wt, double wt0)
 693 {
 694   m_wtEW  = wt;
 695   m_wtEW0 = wt0;
 696 }
 697 void Tauola::setHelicities(int minus, int plus)
 698 {
 699   m_helMinus = minus;
 700   m_helPlus  = plus;
 701 }
 702 void Tauola::setEtaK0sPi(int eta, int k, int pi)
 703 {
 704   ion[0] = pi;
 705   ion[1] = k;
 706   ion[2] = eta;
 707 }
 708
 709 void Tauola::summary()
 710 {
 711   int sign_type=100;
 712   double pol[4] = {0};
 713
 714   Log::Info()     <<"Tauola::summary(): We use old TAUOLA FORTRAN printout."<<endl;
 715   Log::Info(false)<<"As a consequence, there is a mismatch in printed TAUOLA version number."<<endl<<endl;
 716
 717   // Print summary taken from FORTRAN TAUOLA
 718   dekay_(&sign_type,pol);
 719 }
 720
 721 void Tauola::fill_val(int beg, int end, double* array, double value)
 722 {
 723   for (int i = beg; i < end; i++)
 724     array[i] = value;
 725 }
 726
 727 double Tauola::particleCharge(int idhep)
 728 {
 729   static double CHARGE[101] = { 0 };
 730   static int j=0;
 731   //--
 732   //--   Array 'SPIN' contains the spin  of  the first 100 particles accor-
 733   //--   ding to the PDG particle code...
 734
 735   if(j==0) // initialization
 736     {
 737       j=1;
 738       fill_val(0 ,  1, CHARGE, 0.0         );
 739       fill_val(1 ,  2, CHARGE,-0.3333333333);
 740       fill_val(2 ,  3, CHARGE, 0.6666666667);
 741       fill_val(3 ,  4, CHARGE,-0.3333333333);
 742       fill_val(4 ,  5, CHARGE, 0.6666666667);
 743       fill_val(5 ,  6, CHARGE,-0.3333333333);
 744       fill_val(6 ,  7, CHARGE, 0.6666666667);
 745       fill_val(7 ,  8, CHARGE,-0.3333333333);
 746       fill_val(8 ,  9, CHARGE, 0.6666666667);
 747       fill_val(9 , 11, CHARGE, 0.0         );
 748       fill_val(11 ,12, CHARGE,-1.0         );
 749       fill_val(12 ,13, CHARGE, 0.0         );
 750       fill_val(13 ,14, CHARGE,-1.0         );
 751       fill_val(14, 15, CHARGE, 0.0         );
 752       fill_val(15 ,16, CHARGE,-1.0         );
 753       fill_val(16, 17, CHARGE, 0.0         );
 754       fill_val(17 ,18, CHARGE,-1.0         );
 755       fill_val(18, 24, CHARGE, 0.0         );
 756       fill_val(24, 25, CHARGE, 1.0         );
 757       fill_val(25, 37, CHARGE, 0.0         );
 758       fill_val(37, 38, CHARGE, 1.0         );
 759       fill_val(38,101, CHARGE, 0.0         );
 760     }
 761
 762   int idabs=abs(idhep);
 763   double phoch=0.0;
 764
 765   //--
 766   //--   Charge of quark, lepton, boson etc....
 767   if (idabs<=100) phoch=CHARGE[idabs];
 768   else {
 769     int Q3= idabs/1000 % 10;
 770     int Q2= idabs/100  % 10;
 771     int Q1= idabs/10   % 10;
 772     if (Q3==0){
 773       //--
 774       //-- ...meson...
 775       if(Q2 % 2==0) phoch=CHARGE[Q2]-CHARGE[Q1];
 776       else          phoch=CHARGE[Q1]-CHARGE[Q2];
 777     }
 778     else{
 779       //--
 780       //--   ...diquarks or baryon.
 781       phoch=CHARGE[Q1]+CHARGE[Q2]+CHARGE[Q3];
 782     }
 783   }
 784   //--
 785   //--   Find the sign of the charge...
 786   if (idhep<0.0) phoch=-phoch;
 787   if (phoch*phoch<0.000001) phoch=0.0;
 788
 789   return phoch;
 790 }
 791
 792 } // namespace Tauolapp