6 #include "TauolaEvent.h"
10 bool Tauola::m_is_initialized =
false;
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}};
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}};
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}};
53 double Tauola::sminA = 0;
54 double Tauola::smaxA = 0;
56 double Tauola::sminB = 0;
57 double Tauola::smaxB = 0;
59 double Tauola::sminC = 0;
60 double Tauola::smaxC = 0;
62 int Tauola::ion[3] = {0};
63 double Tauola::tau_lifetime = .08711;
64 double Tauola::momentum_conservation_threshold = 0.1;
66 Tauola::Particles Tauola::spin_correlation;
69 Tauola::LengthUnits Tauola::lengthUnit = Tauola::DEFAULT_LENGTH;
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.;
80 double Tauola::buf_R[4][4] = {{0.0}};
83 void (*Tauola::redefineTauPlusProperties)(
TauolaParticle *) = defaultRedPlus;
84 void (*Tauola::redefineTauMinusProperties)(
TauolaParticle *) = defaultRedMinus;
93 return rand()*1./RAND_MAX;
98 else randomDouble = gen;
104 void Tauola::setRedefineTauMinus(
void (*fun)(
TauolaParticle *) ){
105 redefineTauMinusProperties=fun;
108 void Tauola::setRedefineTauPlus (
void (*fun)(
TauolaParticle *) ){
109 redefineTauPlusProperties=fun;
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;
121 Tauola::momentumUnit = m;
122 Tauola::lengthUnit = l;
131 printf(
" *************************************\n");
132 printf(
" * TAUOLA C++ Interface v1.1.8 *\n");
133 printf(
" *-----------------------------------*\n");
135 printf(
" * (c) Nadia Davidson, (1,2) *\n");
136 printf(
" * (c) 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");
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");
148 if(m_is_initialized) {
151 printf(
" *WARNING: Tauola already initialized*\n");
152 printf(
" * reinitialization skipped *\n");
155 printf(
" *************************************\n");
159 m_is_initialized =
true;
162 spin_correlation.setAll(
true);
165 f_interface_tauolaInitialize(m_pdg_id,m_firstDecayMode,
166 m_secondDecayMode,m_rad,
167 m_rad_cut_off, m_iniphy);
172 cout<<
"Reading SANC input files."<<endl;
174 ifstream f(
"table1-1.txt");
177 cout<<
"File 'table1-1.txt' missing... skipped."<<endl;
182 cout<<
"Reading file 'table1-1.txt'..."<<endl;
184 int dbuf1,dbuf2,dbuf3,dbufcos;
185 f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
188 if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS) {
189 cout<<
"mismatched NS1= "<<dbuf1 <<
" <--> "<<NS1<<endl;
190 cout<<
" NS2= "<<dbuf2 <<
" <--> "<<NS2<<endl;
191 cout<<
" NS3= "<<dbuf3 <<
" <--> "<<NS3<<endl;
192 cout<<
" NCOS= "<<dbufcos<<
" <--> "<<NCOS<<endl;
196 double buf1,buf2,buf3,buf4,buf5,buf6;
197 f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
210 if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
211 cout<<
"mismatched sminA= "<<buf1<<
" <--> "<<sminA<<endl;
212 cout<<
" smaxA= "<<buf2<<
" <--> "<<smaxA<<endl;
213 cout<<
" sminB= "<<buf3<<
" <--> "<<sminB<<endl;
214 cout<<
" smaxB= "<<buf4<<
" <--> "<<smaxB<<endl;
215 cout<<
" sminC= "<<buf5<<
" <--> "<<sminC<<endl;
216 cout<<
" smaxC= "<<buf6<<
" <--> "<<smaxC<<endl;
224 if(strcmp(head,
"BeginRange1")==0)
break;
229 for (
int i=0;i<NS1;i++){
230 for (
int j=0;j<NCOS;j++){
231 for (
int k=0;k<4;k++){
232 for (
int l=0;l<4;l++){
233 f>>table1A[i][j][k][l];
244 if(strcmp(buf.c_str(),
"BeginRange2")==0)
break;
248 for (
int i=0;i<NS2;i++){
249 for (
int j=0;j<NCOS;j++){
250 for (
int k=0;k<4;k++){
251 for (
int l=0;l<4;l++){
252 f>>table1B[i][j][k][l];
263 if(strcmp(buf.c_str(),
"BeginRange3")==0)
break;
267 for (
int i=0;i<NS3;i++){
268 for (
int j=0;j<NCOS;j++){
269 for (
int k=0;k<4;k++){
270 for (
int l=0;l<4;l++){
271 f>>table1C[i][j][k][l];
281 if(buf.size() == 0 || strcmp(buf.c_str(),
"End") != 0){
282 cout<<
"...incorrect file version or file incomplete/damaged!"<<endl;
285 table1A[0][0][0][0] = table1B[0][0][0][0] = table1C[0][0][0][0] = 0.0;
290 f.open(
"table2-2.txt");
293 cout<<
"File 'table2-2.txt' missing... skipped."<<endl;
298 cout<<
"Reading file 'table2-2.txt'..."<<endl;
300 int dbuf1,dbuf2,dbuf3,dbufcos;
301 f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
304 if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS) {
305 cout<<
"mismatched NS1= "<<dbuf1<<
" <--> "<<NS1<<endl;
306 cout<<
" NS2= "<<dbuf2<<
" <--> "<<NS2<<endl;
307 cout<<
" NS3= "<<dbuf3<<
" <--> "<<NS3<<endl;
308 cout<<
" NCOS= "<<dbufcos<<
" <--> "<<NCOS<<endl;
312 double buf1,buf2,buf3,buf4,buf5,buf6;
313 f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
327 if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
328 cout<<
"mismatched sminA= "<<buf1<<
" <--> "<<sminA<<endl;
329 cout<<
" smaxA= "<<buf2<<
" <--> "<<smaxA<<endl;
330 cout<<
" sminB= "<<buf3<<
" <--> "<<sminB<<endl;
331 cout<<
" smaxB= "<<buf4<<
" <--> "<<smaxB<<endl;
332 cout<<
" sminC= "<<buf5<<
" <--> "<<sminC<<endl;
333 cout<<
" smaxC= "<<buf6<<
" <--> "<<smaxC<<endl;
341 if(strcmp(head,
"BeginRange1")==0)
break;
346 for (
int i=0;i<NS1;i++){
347 for (
int j=0;j<NCOS;j++){
348 for (
int k=0;k<4;k++){
349 for (
int l=0;l<4;l++){
350 f>>table2A[i][j][k][l];
361 if(strcmp(buf.c_str(),
"BeginRange2")==0)
break;
365 for (
int i=0;i<NS2;i++){
366 for (
int j=0;j<NCOS;j++){
367 for (
int k=0;k<4;k++){
368 for (
int l=0;l<4;l++){
369 f>>table2B[i][j][k][l];
380 if(strcmp(buf.c_str(),
"BeginRange3")==0)
break;
384 for (
int i=0;i<NS3;i++){
385 for (
int j=0;j<NCOS;j++){
386 for (
int k=0;k<4;k++){
387 for (
int l=0;l<4;l++){
388 f>>table2C[i][j][k][l];
398 if(buf.size()==0 || strcmp(buf.c_str(),
"End")!=0){
399 cout<<
"...incorrect file version or file incomplete/damaged!"<<endl;
402 table2A[0][0][0][0] = table2B[0][0][0][0] = table2C[0][0][0][0] = 0.0;
407 f.open(
"table11-11.txt");
410 cout<<
"File 'table11-11.txt' missing... skipped."<<endl;
415 cout<<
"Reading file 'table11-11.txt'..."<<endl;
417 int dbuf1,dbuf2,dbuf3,dbufcos;
418 f>>buf>>dbuf1>>dbuf2>>dbuf3>>dbufcos;
421 if(dbuf1!=NS1 || dbuf2!=NS2 || dbuf3!=NS3 || dbufcos!=NCOS) {
422 cout<<
"mismatched NS1= "<<dbuf1<<
" <--> "<<NS1<<endl;
423 cout<<
" NS2= "<<dbuf2<<
" <--> "<<NS2<<endl;
424 cout<<
" NS3= "<<dbuf3<<
" <--> "<<NS3<<endl;
425 cout<<
" NCOS= "<<dbufcos<<
" <--> "<<NCOS<<endl;
429 double buf1,buf2,buf3,buf4,buf5,buf6;
430 f>>buf>>buf1>>buf2>>buf3>>buf4>>buf5>>buf6;
444 if(buf1!=sminA || buf2!=smaxA || buf3!=sminB || buf4!=smaxB || buf5!=sminC || buf6!=smaxC) {
445 cout<<
"mismatched sminA= "<<buf1<<
" <--> "<<sminA<<endl;
446 cout<<
" smaxA= "<<buf2<<
" <--> "<<smaxA<<endl;
447 cout<<
" sminB= "<<buf3<<
" <--> "<<sminB<<endl;
448 cout<<
" smaxB= "<<buf4<<
" <--> "<<smaxB<<endl;
449 cout<<
" sminC= "<<buf5<<
" <--> "<<sminC<<endl;
450 cout<<
" smaxC= "<<buf6<<
" <--> "<<smaxC<<endl;
458 if(strcmp(head,
"BeginRange1")==0)
break;
463 for (
int i=0;i<NS1;i++){
464 for (
int j=0;j<NCOS;j++){
465 for (
int k=0;k<4;k++){
466 for (
int l=0;l<4;l++){
467 f>>table11A[i][j][k][l];
478 if(strcmp(buf.c_str(),
"BeginRange2")==0)
break;
482 for (
int i=0;i<NS2;i++){
483 for (
int j=0;j<NCOS;j++){
484 for (
int k=0;k<4;k++){
485 for (
int l=0;l<4;l++){
486 f>>table11B[i][j][k][l];
497 if(strcmp(buf.c_str(),
"BeginRange3")==0)
break;
501 for (
int i=0;i<NS3;i++){
502 for (
int j=0;j<NCOS;j++){
503 for (
int k=0;k<4;k++){
504 for (
int l=0;l<4;l++){
505 f>>table11C[i][j][k][l];
514 if(buf.size()==0 || strcmp(buf.c_str(),
"End")!=0){
515 cout<<
"...incorrect file version or file incomplete/damaged!"<<endl;
518 table11A[0][0][0][0] = table11B[0][0][0][0] = table11C[0][0][0][0] = 0.0;
530 if(polx*polx+poly*poly+polz*polz>1)
532 Log::Warning()<<
"decayOne(): ignoring wrong polarization vector: "<<polx<<
" "<<poly<<
" "<<polz<<endl;
554 std::vector<TauolaParticle *> list;
568 Log::Warning() <<
"Deprecated routine 'Tauola::initialise'"<<endl;
569 Log::Warning(0)<<
"Use 'Tauola::initialize' instead."<<endl;
578 return m_is_initialized;
611 return abs(m_pdg_id);
615 m_firstDecayMode=firstDecayMode;
619 m_secondDecayMode=secondDecayMode;
627 m_rad_cut_off=rad_cut_off;
637 Log::Warning() <<
"Deprecated routine 'Tauola::setInitialisePhy'"<<endl;
638 Log::Warning(0)<<
"Use 'Tauola::setInitializePhy' instead."<<endl;
649 Log::Warning()<<
"setTauBr(): run Tauola::initialize() first."<<endl;
650 else if(i<1 || i>taubra_.nchan || value<0.)
651 Log::Warning()<<
"setTauBr(): Invalid input. Value must be >= 0 and 0 < i <= "<<taubra_.nchan<<endl;
652 else taubra_.gamprt[i-1]=(float)value;
655 void Tauola::setTaukle(
double bra1,
double brk0,
double brk0b,
double brks)
657 if(bra1<0 || bra1>1 || brk0<0 ||brk0>1 || brk0b<0 || brk0b>1 || brks<0 ||brks>1)
659 Log::Warning()<<
"setTaukle(): variables must be in range [0,1]. Ignored."<<endl;
662 taukle_.bra1 =(float)bra1;
663 taukle_.brk0 =(float)brk0;
664 taukle_.brk0b=(float)brk0b;
665 taukle_.brks =(float)brks;
669 return f_getTauMass();
672 double Tauola::getHiggsScalarPseudoscalarMixingAngle(){
673 return m_higgs_scalar_pseudoscalar_mix;
677 return m_higgs_scalar_pseudoscalar_pdg;
682 m_higgs_scalar_pseudoscalar_mix=angle;
690 Log::Warning()<<
"You want to use spin correlations of Higgs for particle of PDGID= "<<pdg_code<<endl
692 Log::Fatal(
"This choice is not appropriate.",0);
694 m_higgs_scalar_pseudoscalar_pdg=pdg_code;
697 int Tauola::getHelPlus(){
700 int Tauola::getHelMinus(){
703 double Tauola::getEWwt(){
706 double Tauola::getEWwt0(){
709 void Tauola::setEWwt(
double wt,
double wt0)
714 void Tauola::setHelicities(
int minus,
int plus)
719 void Tauola::setEtaK0sPi(
int eta,
int k,
int pi)
726 void Tauola::summary()
731 Log::Info() <<
"Tauola::summary(): We use old TAUOLA FORTRAN printout."<<endl;
732 Log::Info(
false)<<
"As a consequence, there is a mismatch in printed TAUOLA version number."<<endl<<endl;
735 dekay_(&sign_type,pol);
740 for (
int i = beg; i < end; i++)
746 static double CHARGE[101] = { 0 };
756 fill_val(1 , 2, CHARGE,-0.3333333333);
757 fill_val(2 , 3, CHARGE, 0.6666666667);
758 fill_val(3 , 4, CHARGE,-0.3333333333);
759 fill_val(4 , 5, CHARGE, 0.6666666667);
760 fill_val(5 , 6, CHARGE,-0.3333333333);
761 fill_val(6 , 7, CHARGE, 0.6666666667);
762 fill_val(7 , 8, CHARGE,-0.3333333333);
763 fill_val(8 , 9, CHARGE, 0.6666666667);
779 int idabs=abs(idhep);
784 if (idabs<=100) phoch=CHARGE[idabs];
786 int Q3= idabs/1000 % 10;
787 int Q2= idabs/100 % 10;
788 int Q1= idabs/10 % 10;
792 if(Q2 % 2==0) phoch=CHARGE[Q2]-CHARGE[Q1];
793 else phoch=CHARGE[Q1]-CHARGE[Q2];
798 phoch=CHARGE[Q1]+CHARGE[Q2]+CHARGE[Q3];
803 if (idhep<0.0) phoch=-phoch;
804 if (phoch*phoch<0.000001) phoch=0.0;
static void(* m_decay_one_boost_routine)(TauolaParticle *, TauolaParticle *)
static void setBoostRoutine(void(*boost)(TauolaParticle *, TauolaParticle *))
static double defaultRandomGenerator()
static void setRadiationCutOff(double rad_cut_off)
void checkMomentumConservation()
static bool getIsTauolaIni()
static void setHiggsScalarPseudoscalarPDG(int pdg_id)
static void setNewCurrents(int mode)
static const double * getDecayOnePolarization()
static double getTauMass()
static int getDecayingParticle()
Abstract base class for particle in the event. This class also handles boosting.
static void setHiggsScalarPseudoscalarMixingAngle(double angle)
static void setInitialisePhy(double iniphy)
static void decayOne(TauolaParticle *tau, bool undecay=false, double polx=0, double poly=0, double polz=0)
static double m_decay_one_polarization[3]
static int getHiggsScalarPseudoscalarPDG()
static void setSameParticleDecayMode(int firstDecayMode)
static void setTauBr(int i, double value)
static void setTauLifetime(double t)
static void fill_val(int beg, int end, double *array, double value)
static void decayOneBoost(TauolaParticle *mother, TauolaParticle *target)
static void setDecayingParticle(int pdg_id)
static bool isUsingDecayOneBoost()
static bool m_is_using_decay_one
static void setOppositeParticleDecayMode(int secondDecayMode)
static void setRandomGenerator(double(*gen)())
static bool isUsingDecayOne()
static void setInitializePhy(double iniphy)
static double particleCharge(int idhep)
static void setRadiation(bool rad)
static void setUnits(MomentumUnits m, LengthUnits l)
static void Fatal(string text, unsigned short int code=0)