]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/Tauola/Tauola.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes for object target dependencies
[u/mrichter/AliRoot.git] / TEvtGen / Tauola / Tauola.cxx
CommitLineData
0ca57c2f 1#include <fstream>
2#include <cstring>
3#include <vector>
e1938fed 4#include "TauolaLog.h"
0ca57c2f 5#include "Tauola.h"
6#include "TauolaEvent.h"
7
8namespace Tauolapp
9{
10
11int Tauola::m_pdg_id = 15;
12int Tauola::m_firstDecayMode = 0;
13int Tauola::m_secondDecayMode = 0;
14bool Tauola::m_rad = true;
15double Tauola::m_rad_cut_off = 0.001;
16double Tauola::m_iniphy = 0.1;
17double Tauola::m_higgs_scalar_pseudoscalar_mix = M_PI/4;
18int Tauola::m_higgs_scalar_pseudoscalar_pdg = 35;
19int Tauola::m_helPlus = 0;
20int Tauola::m_helMinus = 0;
21double Tauola::m_wtEW = 0.0;
22double Tauola::m_wtEW0 = 0.0;
23double Tauola::table11A[NS1][NCOS][4][4] = {{{{0.0}}}};
24double Tauola::table1A [NS1][NCOS][4][4] = {{{{0.0}}}};
25double Tauola::table2A [NS1][NCOS][4][4] = {{{{0.0}}}};
26double Tauola::wtable11A[NS1][NCOS] = {{0.0}};
27double Tauola::wtable1A [NS1][NCOS] = {{0.0}};
28double Tauola::wtable2A [NS1][NCOS] = {{0.0}};
29double Tauola::w0table11A[NS1][NCOS] = {{0.0}};
30double Tauola::w0table1A [NS1][NCOS] = {{0.0}};
31double Tauola::w0table2A [NS1][NCOS] = {{0.0}};
32
33double Tauola::table11B[NS2][NCOS][4][4] = {{{{0.0}}}};
34double Tauola::table1B [NS2][NCOS][4][4] = {{{{0.0}}}};
35double Tauola::table2B [NS2][NCOS][4][4] = {{{{0.0}}}};
36double Tauola::wtable11B[NS2][NCOS] = {{0.0}};
37double Tauola::wtable1B [NS2][NCOS] = {{0.0}};
38double Tauola::wtable2B [NS2][NCOS] = {{0.0}};
39double Tauola::w0table11B[NS2][NCOS] = {{0.0}};
40double Tauola::w0table1B [NS2][NCOS] = {{0.0}};
41double Tauola::w0table2B [NS2][NCOS] = {{0.0}};
42
43double Tauola::table11C[NS3][NCOS][4][4] = {{{{0.0}}}};
44double Tauola::table1C [NS3][NCOS][4][4] = {{{{0.0}}}};
45double Tauola::table2C [NS3][NCOS][4][4] = {{{{0.0}}}};
46double Tauola::wtable11C[NS3][NCOS] = {{0.0}};
47double Tauola::wtable1C [NS3][NCOS] = {{0.0}};
48double Tauola::wtable2C [NS3][NCOS] = {{0.0}};
49double Tauola::w0table11C[NS3][NCOS] = {{0.0}};
50double Tauola::w0table1C [NS3][NCOS] = {{0.0}};
51double Tauola::w0table2C [NS3][NCOS] = {{0.0}};
52
53double Tauola::sminA = 0;
54double Tauola::smaxA = 0;
55
56double Tauola::sminB = 0;
57double Tauola::smaxB = 0;
58
59double Tauola::sminC = 0;
60double Tauola::smaxC = 0;
61
62int Tauola::ion[3] = {0};
63double Tauola::tau_lifetime = .08711;
64double Tauola::momentum_conservation_threshold = 0.1;
65
66Tauola::Particles Tauola::spin_correlation;
67
68Tauola::MomentumUnits Tauola::momentumUnit = Tauola::DEFAULT_MOMENTUM;
69Tauola::LengthUnits Tauola::lengthUnit = Tauola::DEFAULT_LENGTH;
70
71bool Tauola::m_is_using_decay_one = false;
72double Tauola::m_decay_one_polarization[3] = {0};
73void (*Tauola::m_decay_one_boost_routine)(TauolaParticle*,TauolaParticle*) = NULL;
74
75int Tauola::buf_incoming_pdg_id = 0;
76int Tauola::buf_outgoing_pdg_id = 0;
77double Tauola::buf_invariant_mass_squared = -1.;
78double Tauola::buf_cosTheta = 0.;
79
80double Tauola::buf_R[4][4] = {{0.0}}; //density matrix
81
82double (*Tauola::randomDouble)() = Tauola::defaultRandomGenerator;
83void (*Tauola::redefineTauPlusProperties)(TauolaParticle *) = defaultRedPlus;
84void (*Tauola::redefineTauMinusProperties)(TauolaParticle *) = defaultRedMinus;
85
86/**************************************************************/
87void Tauola::setNewCurrents(int mode)
88{
89 inirchl_(&mode);
90}
91
92double Tauola::defaultRandomGenerator(){
93 return rand()*1./RAND_MAX;
94}
95
96void Tauola::setRandomGenerator(double (*gen)()){
97 if(gen==NULL) randomDouble = defaultRandomGenerator;
98 else randomDouble = gen;
99}
100
101void Tauola::defaultRedPlus(TauolaParticle *tau) {}
102void Tauola::defaultRedMinus(TauolaParticle *tau) {}
103
104void Tauola::setRedefineTauMinus( void (*fun)(TauolaParticle *) ){
105 redefineTauMinusProperties=fun;
106}
107
108void Tauola::setRedefineTauPlus ( void (*fun)(TauolaParticle *) ){
109 redefineTauPlusProperties=fun;
110}
111
112void 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
120void Tauola::setUnits(MomentumUnits m, LengthUnits l){
121 Tauola::momentumUnit = m;
122 Tauola::lengthUnit = l;
123}
124
125void Tauola::setTauLifetime(double t){
126 tau_lifetime = t;
127}
128
129void 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
513void 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
553void 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
564bool Tauola::isUsingDecayOne()
565{
566 return m_is_using_decay_one;
567}
568
569bool Tauola::isUsingDecayOneBoost()
570{
571 return (bool) m_decay_one_boost_routine;
572}
573
574void Tauola::setBoostRoutine(void (*boost)(TauolaParticle *, TauolaParticle *))
575{
576 m_decay_one_boost_routine=boost;
577}
578
579void Tauola::decayOneBoost(TauolaParticle *mother, TauolaParticle *target)
580{
581 m_decay_one_boost_routine(mother,target);
582}
583
584const double* Tauola::getDecayOnePolarization()
585{
586 return m_decay_one_polarization;
587}
588
589void Tauola::setDecayingParticle(int pdg_id){
590 m_pdg_id=pdg_id;
591}
592
593int Tauola::getDecayingParticle(){
594 return abs(m_pdg_id);
595}
596
597void Tauola::setSameParticleDecayMode(int firstDecayMode){
598 m_firstDecayMode=firstDecayMode;
599}
600
601void Tauola::setOppositeParticleDecayMode(int secondDecayMode){
602 m_secondDecayMode=secondDecayMode;
603}
604
605void Tauola::setRadiation(bool rad){
606 m_rad=rad;
607}
608
609void Tauola::setRadiationCutOff(double rad_cut_off){
610 m_rad_cut_off=rad_cut_off;
611}
612
613
614void Tauola::setInitializePhy(double iniphy){
615 m_iniphy=iniphy;
616}
617
618void 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
629void 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
638void 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
651double Tauola::getTauMass(){
652 return f_getTauMass();
653}
654
655double Tauola::getHiggsScalarPseudoscalarMixingAngle(){
656 return m_higgs_scalar_pseudoscalar_mix;
657}
658
659int Tauola::getHiggsScalarPseudoscalarPDG(){
660 return m_higgs_scalar_pseudoscalar_pdg;
661}
662
663/** set the mixing angle. coupling: tau~(cos(phi)+isin(phi)gamma5)tau */
664void 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 */
670void 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
680int Tauola::getHelPlus(){
681 return m_helPlus;
682}
683int Tauola::getHelMinus(){
684 return m_helMinus;
685}
686double Tauola::getEWwt(){
687 return m_wtEW;
688}
689double Tauola::getEWwt0(){
690 return m_wtEW0;
691}
692void Tauola::setEWwt(double wt, double wt0)
693{
694 m_wtEW = wt;
695 m_wtEW0 = wt0;
696}
697void Tauola::setHelicities(int minus, int plus)
698{
699 m_helMinus = minus;
700 m_helPlus = plus;
701}
702void Tauola::setEtaK0sPi(int eta, int k, int pi)
703{
704 ion[0] = pi;
705 ion[1] = k;
706 ion[2] = eta;
707}
708
709void 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
721void 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
727double 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