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