C++ Interface to Tauola
TauolaParticlePair.cxx
1#include "Tauola.h"
2#include "TauolaParticlePair.h"
3#include "Log.h"
4#include <stdlib.h>
5#include <math.h>
6#include <iostream>
7
8namespace Tauolapp
9{
10
11/** constructor. Get the mothers, grandmothers and siblings of the tau */
12TauolaParticlePair::TauolaParticlePair(std::vector<TauolaParticle*> &particle_list){
13 TauolaParticle *particle=particle_list.back();
14 particle_list.pop_back();
15 setBornKinematics(0,0,-1.,0.); //set the default born variables
16
17 // In case of decayOne() we need only the tau information
19 {
20 m_mother = 0;
21 m_mother_exists = false;
22 m_production_particles.push_back(particle);
23 m_final_particles.push_back(particle);
26 return;
27 }
28
29 //See if there are any mothers
30 std::vector<TauolaParticle *> temp_mothers = particle->findProductionMothers();
31
32 // Case that there are no mothers or grandmothers
33 if(temp_mothers.size()==0){
34 // NOTE: Not executed by release examples
35 // However, such cases were present if tests used older Pythia8.1 or so.
36 Log::Warning()<< "WARNING: Could not find taus mother or grandmothers. "
37 << "Ignoring spin effects" << std::endl;
38 m_mother = 0;
39 m_mother_exists = false;
40 m_production_particles.push_back(particle);
41 m_final_particles.push_back(particle->findLastSelf());
44 return;
45 }
46
47 //fill the sibling pointers
48 std::vector<TauolaParticle*> temp_daughters;
49 temp_daughters = temp_mothers.at(0)->getDaughters();
50 if(temp_daughters.size()==0)
51 Log::Fatal("WARNING: Something wrong with event structure or there is a bug in the TAUOLA interface.",6);
52
53 m_production_particles=temp_daughters;
54 m_final_particles.push_back(particle->findLastSelf());
55
56 //Second tau-like particle selection should match properties of the hard (exotic) process. At present, solution
57 //match one of the possible diagrams contributing to SM process. Rare case like tau+ tau+ nutau nutau will not be treted
58 //accordingly to any diagram of SM
59 for(signed int i=0; i < (int) m_production_particles.size(); i++)
60 {
61 //check if it has opposite PDGID sign
62 if(m_final_particles.at(0)->getPdgID()*m_production_particles.at(i)->getPdgID()>=0) continue;
63
66 {
67 //tau+ or tau- - check if it's on the particle_list
68 int j=-1;
69 for(j=0;j<(int)particle_list.size();j++)
70 if(m_production_particles.at(i)->getBarcode()==particle_list.at(j)->getBarcode()) break;
71 if(j>=(int)particle_list.size()) continue;
72
73 //exists on the list - add to m_final_particles and delete from particle_list
74 m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
75 particle_list.erase(particle_list.begin()+j);
76 }
77 else if(m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_NEUTRINO||
79 {
80 //neutrino - for now - just add to m_final_particles
81 m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
82 }
83 }
84
85
86 //fill the mother and grandmother pointers
87 if(temp_mothers.size()==1){ //one mother
88 m_mother_exists = true;
89 m_mother = temp_mothers.at(0);
92 }
93 else{ //no mother, but grandparents exist
94 m_mother_exists = false;
95 m_grandmothers = temp_mothers;
98 }
99
101
102 return;
103}
104
105
106/** The axis is defined by the boosting routine but our standard convention
107 is:
108 - Axis 3 is along the direction of the +ve tau,
109 - Axis 1 is perpendicular to the reaction plane (sign=??)
110 - Axis 2 is defined through the vector product so the system
111 is right handed (?? check). Axis 1,2 and 3 are parrellel for the 1st
112 and second tau. */
114 int incoming_pdg_id=0;
115 int outgoing_pdg_id=0;
116 double invariant_mass_squared=-5.0;
117 double cosTheta=3.0;
118
119 //initialize all elements of the density matrix to zero
120 for(int x = 0; x < 4; x ++) {
121 for(int y = 0; y < 4; y ++)
122 m_R[x][y] = 0;
123 }
124
125 m_R[0][0]=1;
126
128 {
129 const double *pol = Tauola::getDecayOnePolarization();
130
131 m_R[0][1]=pol[0];
132 m_R[0][2]=pol[1];
133 m_R[0][3]=pol[2];
134
135 m_R[1][0]=pol[0];
136 m_R[2][0]=pol[1];
137 m_R[3][0]=pol[2];
138 }
139
140 if(!m_mother) return;
141 // fill the matrix depending on mother
142
143
144 // do scalar-pseudoscalar mixed higgs case separately since
145 // switch doesn't allow non-constants in a case statement.
146 // very annoying!
147
149 if(!Tauola::spin_correlation.HIGGS_H) return;
150
151 double phi = Tauola::getHiggsScalarPseudoscalarMixingAngle();
152 double mass_ratio = Tauola::getTauMass()/m_mother->getMass();
153
154 double beta = sqrt(1-4*mass_ratio*mass_ratio);
155 double denominator = pow(cos(phi)*beta,2)+pow(sin(phi),2);
156
157 m_R[0][0]= 1;
158 m_R[1][1]= (pow(cos(phi)*beta,2)-pow(sin(phi),2))/denominator;
159 m_R[1][2]= 2*cos(phi)*sin(phi)*beta/denominator;
160 m_R[2][1]= -m_R[1][2];
161 m_R[2][2]= m_R[1][1];
162 m_R[3][3]= -1;
163 }
164 else {
165
166 double pz = 0.0;
167
168 switch(m_mother->getPdgID()){
169
171 if(!Tauola::spin_correlation.Z0) break;
172 // Here we calculate SVAR and COSTHE as well as IDE and IDF
173 // ANGULU(&IDE,&IDF,&SVAR,&COSTHE);
174 // this is ++ configuration of Fig 5 from HEPPH0101311:
175 //pol_tau_minus[2]=-1; pol_tau_plus[2]=1;
176
177 pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
178 m_R[0][0]=1;
179 m_R[0][3]=2*pz-1;
180 m_R[3][0]=2*pz-1;
181 m_R[3][3]=1;
182 // alternatively this may be overwritten if better solution exist
183 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
184 break;
185
187 if(!Tauola::spin_correlation.GAMMA) break;
188 // Here we calculate SVAR and COSTHE as well as IDE and IDF
189 // ANGULU(&IDE,&IDF,&SVAR,&COSTHE);
190 // this is ++ configuration of Fig 5 from HEPPH0101311:
191 //pol_tau_minus[2]=-1; pol_tau_plus[2]=1;
192
193 pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
194 m_R[0][0]=1;
195 m_R[0][3]=2*pz-1;
196 m_R[3][0]=2*pz-1;
197 m_R[3][3]=1;
198 // alternatively this may be overwritten if better solution exist
199 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
200 break;
201
202 //scalar higgs
204 if(!Tauola::spin_correlation.HIGGS) break;
205 m_R[0][0]=1;
206 m_R[1][1]=1;
207 m_R[2][2]=1;
208 m_R[3][3]=-1;
209 break;
210
211 //pseudoscalar higgs case
213 if(!Tauola::spin_correlation.HIGGS_A) break;
214 m_R[0][0]=1;
215 m_R[1][1]=-1;
216 m_R[2][2]=-1;
217 m_R[3][3]=-1;
218 break;
219
221 if(!Tauola::spin_correlation.HIGGS_PLUS) break;
222 m_R[0][0]=1;
223 m_R[0][3]=-1;
224 m_R[3][0]=-1;
225 break;
226
228 if(!Tauola::spin_correlation.HIGGS_MINUS) break;
229 m_R[0][0]=1;
230 m_R[0][3]=-1;
231 m_R[3][0]=-1;
232 break;
233
235 if(!Tauola::spin_correlation.W_PLUS) break;
236 m_R[0][0]=1;
237 m_R[0][3]=1; //tau minus (tau minus is on the -ve Z axis)
238 m_R[3][0]=1; //tau plus
239 break;
240
242 if(!Tauola::spin_correlation.W_MINUS) break;
243 m_R[0][0]=1;
244 m_R[0][3]=1; //tau minus (tau minus is on the -ve Z axis)
245 m_R[3][0]=1; //tau plus
246 break;
247
248 //ignore spin effects when mother is unknown
249 default:
250 m_R[0][0]=1;
251 break;
252 }
253 }
254
255}
256
257/**************************************************************/
258void TauolaParticlePair::setBornKinematics(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared,double cosTheta){
259 Tauola::buf_incoming_pdg_id=incoming_pdg_id;
260 Tauola::buf_outgoing_pdg_id=outgoing_pdg_id;
261 Tauola::buf_invariant_mass_squared=invariant_mass_squared;
262 Tauola::buf_cosTheta=cosTheta;
263 //cout<<"(TauolaParticlePair::Just to be sure:) "<<buf_incoming_pdg_id<<" "<<buf_outgoing_pdg_id<<" "<<buf_invariant_mass_squared<<" "<<buf_cosTheta<<endl;
264 // m_R[0][0] to be added in next step;
265}
266
267/**************************************************************/
268double TauolaParticlePair::getZPolarization(int *incoming_pdg_id, int *outgoing_pdg_id, double *invariant_mass_squared,double *cosTheta){
269
270 //defaults
271 *incoming_pdg_id = TauolaParticle::ELECTRON;
272 *cosTheta = 0 ;
273 *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
274 *invariant_mass_squared = pow(m_mother->getMass(),2);
275 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
276
277 //TRIVIAL CASE:
278 //if we don't know the incoming beams then
279 //return the average z polarisation
280 if(m_grandmothers.size()<2){
281 Log::Warning()<<"Not enough mothers of Z to "
282 <<"calculate cos(theta) between "
283 <<"incoming about outgoing beam"
284 << endl;
285 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
286 invariant_mass_squared, cosTheta);
287 }
288
290 Log::Error()<<"tau+ or tau- not found in Z decay"<< endl;
291 return 0;
292 }
293
294 //NOW CHECK FOR THE DIFFICULT EVENTS:
295 //case f1 + f2 + f3 -> Z -> tau tau
296 //case f1 + f2 -> Z + f3, Z-> tau tau or f1 + f2 -> tau tau f3
297 //case f1 + f2 -> Z -> tau tau gamma
298 if(m_grandmothers.size()>2 ||
299 (m_grandmothers.at(0)->getDaughters().size()>1 && m_mother_exists==true) ||
300 m_production_particles.size() > 2){
301
302 //make a vector of the extra grandmother particles
303 vector<TauolaParticle*> extra_grandmothers;
304 for(int i=0; i<(int) m_grandmothers.size(); i++){
305 // temp_grandmothers.push_back(m_grandmothers.at(i));
308 extra_grandmothers.push_back(m_grandmothers.at(i));
309 }
310
311 //make a vector of the tau siblings
312 //and copy all the production particle vector.
313 vector<TauolaParticle*> extra_tau_siblings;
314 vector<TauolaParticle*> temp_production_particles;
315 for(int i=0; i<(int) m_production_particles.size(); i++){
318 extra_tau_siblings.push_back(m_production_particles.at(i));
319 }
320
321 //make a vector of the Z's sibling
322 vector<TauolaParticle*> extra_Z_siblings;
323 for(int i=0; m_mother_exists && i<(int) m_grandmothers.at(0)->getDaughters().size(); i++){
324 if(m_grandmothers.at(0)->getDaughters().at(i)->getPdgID()!=TauolaParticle::Z0)
325 extra_Z_siblings.push_back(m_grandmothers.at(0)->getDaughters().at(i));
326 }
327
328 //make temporary particles for the effect beams
329 //and copy into a vector
330 std::vector<TauolaParticle*> effective_taus;
331 effective_taus.push_back(getTauPlus(m_production_particles)->clone());
332 effective_taus.push_back(getTauMinus(m_production_particles)->clone());
333
334 //copy grandmothers into the m_grandmothers vector since we want this to be perminante.
337 m_grandmothers.clear();
338 m_grandmothers.push_back(g1);
339 m_grandmothers.push_back(g2);
340
341 //loop over extra grandmothers
342 for(int i=0; i<(int) extra_grandmothers.size(); i++)
343 addToBeam(extra_grandmothers.at(i),&m_grandmothers,0);
344
345 //loop over siblings to the Z
346 for(int i=0; i<(int) extra_Z_siblings.size(); i++)
347 addToBeam(extra_Z_siblings.at(i),0,&m_grandmothers);
348
349 //loop over siblings of the taus
350 for(int i=0; i<(int) extra_tau_siblings.size() ; i++){
352 addToBeam(extra_tau_siblings.at(i),&effective_taus,0);
353 else
354 addToBeam(extra_tau_siblings.at(i),&effective_taus,&m_grandmothers);
355 }
356 //And we are finish making the effective income and outgoing beams
357
358 TauolaParticle * temp_mother = makeTemporaryMother(effective_taus);
359 *invariant_mass_squared = pow(temp_mother->getMass(),2);
360
361 //now we are ready to do the boosting, calculate the angle
362 //between incoming and outgoing, and get the polarisation of the Z
363
364 double angle1,angle2,angle3;
365 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3, temp_mother,
366 m_grandmothers, effective_taus);
367
368 /*double theta1 = -getGrandmotherPlus(m_grandmothers)->getRotationAngle(TauolaParticle::Y_AXIS);
369 double theta2 = -(M_PI+getGrandmotherMinus(m_grandmothers)->getRotationAngle(TauolaParticle::Y_AXIS));
370 *cosTheta = cos((theta1+theta2)/2.0); //just average the angles for now.*/
371
372 TauolaParticle *tM=getTauPlus(effective_taus);
374
375 double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
376 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
377 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
378 tM=getTauMinus(effective_taus);
380
381 double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
382 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
383 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
384 double sintheta1 = sqrt(1-costheta1*costheta1);
385 double sintheta2 = sqrt(1-costheta2*costheta2);
386 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
387
388 *cosTheta = avg;
389
390 *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
391
392 if(abs(*incoming_pdg_id)==21 || abs(*incoming_pdg_id)==22)
393 {
394 *incoming_pdg_id = -getGrandmotherMinus(m_grandmothers)->getPdgID();
395 //cout<<"INFO:\tgluon pdg id changed!"<<endl;
396 }
397
398 if( abs(*incoming_pdg_id)> 8 &&
399 abs(*incoming_pdg_id)!=11 &&
400 abs(*incoming_pdg_id)!=13 &&
401 abs(*incoming_pdg_id)!=15 )
402 {
403 Log::Error()<<"Second class disaster: incoming_pdg_id = "<<*incoming_pdg_id<<endl;
404
405 // Return avarage Z polarization
406 *incoming_pdg_id = TauolaParticle::ELECTRON;
407 *cosTheta = 0 ;
408 *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
409
410 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
411 invariant_mass_squared, cosTheta);
412 }
413
414 boostFromTauPairToLabFrame(angle1, angle2, angle3, temp_mother,
415 m_grandmothers, effective_taus);
416 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
417 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
418
419 }
420 //REGULAR CASE:
421 //f1 + f2 -> Z -> tau+ tau - or f1 f2 -> tau+ tau-
422 //This includes Z -> tau tau, tau -> gamma tau?
423
424 double angle1,angle2,angle3;
425 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
427
430 double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
431 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
432 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
433
436
437 double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
438 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
439 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
440
441 double sintheta1 = sqrt(1-costheta1*costheta1);
442 double sintheta2 = sqrt(1-costheta2*costheta2);
443
444 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
445
446 *cosTheta = avg;
447
448 *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
449
450 boostFromTauPairToLabFrame(angle1, angle2, angle3,
452
453 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta); // store for debugging
454 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
455 // return 0.5 - (-0.12 + 0.12*cosTheta)/2;
456}
457
458/** WHERE WE CALCULATE THE EFFECTIVE BEAMS **/
459/** This is where we decide which particle should be added into which beam,
460 add it and change the flavour if necessary. candidates_same
461 are on the same side of the vertex as the particle. This is needed
462 for negative the particle 4-momentum. **/
464 std::vector<TauolaParticle*> *candidates_same,
465 std::vector<TauolaParticle*> *candidates_opp){
466
467
468 double s=0.,o=-10.;
469 TauolaParticle * s_beam = NULL;
470 TauolaParticle * o_beam = NULL;
471
472 if(candidates_same){
473 double s0 =1.0/getVirtuality(pcle,candidates_same->at(0),false);
474 double s1 = 1.0/getVirtuality(pcle,candidates_same->at(1),false);
475 if(s0>s1){
476 s=s0;
477 s_beam = candidates_same->at(0);
478 }
479 else{
480 s=s1;
481 s_beam = candidates_same->at(1);
482 }
483 }
484 if(candidates_opp){
485
486 double o0 =1.0/getVirtuality(pcle,candidates_opp->at(0),true);
487 double o1 = 1.0/getVirtuality(pcle,candidates_opp->at(1),true);
488 if(o0>o1){
489 o=o0;
490 o_beam = candidates_opp->at(0);
491 }
492 else{
493 o=o1;
494 o_beam = candidates_opp->at(1);
495 }
496 }
497
498 if(s>o)
499 {
500 s_beam->add(pcle);
501 int pdg1 = pcle->getPdgID();
502 int pdg2 = s_beam->getPdgID();
503 if(abs(pdg2)==15) return;
504 if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) s_beam->setPdgID( pdg1);
505 }
506 else
507 {
508 o_beam->subtract(pcle);
509 int pdg1 = pcle->getPdgID();
510 int pdg2 = o_beam->getPdgID();
511 if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) o_beam->setPdgID(-pdg1);
512 }
513
514 //should we also do something with the flavours??
515
516}
517
519
520 //if one particle in a gluon and the other a lepton
521 if((p1->getPdgID()==TauolaParticle::GLUON&&abs(p2->getPdgID())>10&&abs(p2->getPdgID())<20) ||
522 (p2->getPdgID()==TauolaParticle::GLUON&&abs(p1->getPdgID())>10&&abs(p1->getPdgID())<20))
523 return -1;
524
525 //add some others...
526
527 //otherwise we calculate the virtuality:
528 double kp = p1->getE()*p2->getE() - p1->getPx()*p2->getPx()
529 - p1->getPy()*p2->getPy() - p1->getPz()*p2->getPz();
530 if(flip) // Note that we have 1/(p1+p2)^2 or 1/(p1-p2)^2 and p2^2=0
531 kp = kp - p1->getMass()*p1->getMass();
532
533 double q = 1;
535 if(abs(p2->getPdgID())==TauolaParticle::UP ||
536 abs(p2->getPdgID())==TauolaParticle::CHARM ||
537 abs(p2->getPdgID())==TauolaParticle::TOP)
538 q=2.0/3.0;
539 if(abs(p2->getPdgID())==TauolaParticle::DOWN ||
540 abs(p2->getPdgID())==TauolaParticle::STRANGE ||
541 abs(p2->getPdgID())==TauolaParticle::BOTTOM)
542 q=1.0/3.0;
543 }
545 if(abs(p1->getPdgID())==TauolaParticle::UP ||
546 abs(p1->getPdgID())==TauolaParticle::CHARM ||
547 abs(p1->getPdgID())==TauolaParticle::TOP)
548 q=2.0/3.0;
549 if(abs(p1->getPdgID())==TauolaParticle::DOWN ||
550 abs(p1->getPdgID())==TauolaParticle::STRANGE ||
551 abs(p1->getPdgID())==TauolaParticle::BOTTOM)
552 q=1.0/3.0;
553 }
554
555 return fabs(2*kp)/q;
556}
557
558/***********************************************************
559 ** TauolaParticlePair::decayTauPair().
560 ** Generate tau decay
561 ************************************************************/
563
564 //initalize h vectors in case the partner is not a tau
565 double h_tau_minus[4]={2,0,0,0}; //2 used to compensate for maximum weight
566 double h_tau_plus[4]={2,0,0,0}; //in the case when there is only 1 tau
567
570
571 //now calculate the spin weight
572 for(double weight=0; weight < Tauola::randomDouble();){
573 //call tauola decay and get polarimetric vectors
574 if(tau_minus){
575 Tauola::redefineTauMinusProperties(tau_minus);
576 tau_minus->decay();
577 h_tau_minus[0]=1;
578 h_tau_minus[1]=tau_minus->getPolarimetricX();
579 h_tau_minus[2]=tau_minus->getPolarimetricY();
580 h_tau_minus[3]=tau_minus->getPolarimetricZ();
581 }
582 if(tau_plus){
583 Tauola::redefineTauPlusProperties(tau_plus);
584 tau_plus->decay();
585 h_tau_plus[0]=1;
586 h_tau_plus[1]=tau_plus->getPolarimetricX();
587 h_tau_plus[2]=tau_plus->getPolarimetricY();
588 h_tau_plus[3]=tau_plus->getPolarimetricZ();
589 }
590 // cout<<" tau? tau? "<< tau_plus->getPdgID()<<" "<< tau_minus->getPdgID()<<endl;
591 // cout<<" przy uzyciu rrrRRRrrrRRR "<< m_R[0][0]<<" " <<m_R[3][3]<<" " << m_R[0][3]<<" " <<m_R[3][0] <<endl;
592 //calculate the weight
593 weight=0;
594 for(int i =0; i<4; i++){
595 for(int j=0; j<4; j++)
596 weight+=m_R[i][j]*h_tau_minus[i]*h_tau_plus[j];
597 }
598 weight = weight/4.0;
599 }
600 double wthel[4]={0};
601 wthel[0]=(h_tau_minus[0]+h_tau_minus[3])*(h_tau_plus[0]+h_tau_plus[3])*(m_R[0][0]+m_R[0][3]+m_R[3][0]+m_R[3][3]);
602 wthel[1]=(h_tau_minus[0]+h_tau_minus[3])*(h_tau_plus[0]-h_tau_plus[3])*(m_R[0][0]-m_R[0][3]+m_R[3][0]-m_R[3][3]);
603 wthel[2]=(h_tau_minus[0]-h_tau_minus[3])*(h_tau_plus[0]+h_tau_plus[3])*(m_R[0][0]+m_R[0][3]-m_R[3][0]-m_R[3][3]);
604 wthel[3]=(h_tau_minus[0]-h_tau_minus[3])*(h_tau_plus[0]-h_tau_plus[3])*(m_R[0][0]-m_R[0][3]-m_R[3][0]+m_R[3][3]);
605
606 double rn=Tauola::randomDouble();
607 if (rn>(wthel[0]+wthel[1]+wthel[2])/(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,-1);
608 else if (rn>(wthel[0]+wthel[1]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,+1);
609 else if (rn>(wthel[0]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities( 1,-1);
610 else Tauola::setHelicities( 1, 1);
611
612
613
614
615 //boost into the frame used to define the density matrices
616 //and add the tau's new daughters.
617
618 double angle1,angle2,angle3;
620
622 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
624
625 //add the accepted decay into the event record
626 if(tau_plus)
627 tau_plus->addDecayToEventRecord();
628 if(tau_minus)
629 tau_minus->addDecayToEventRecord();
630
632 boostFromTauPairToLabFrame(angle1,angle2,angle3,
634
635 // Apply final decay modification, once taus are in lab frame.
636 // At present 28.02.11, vertex position shift (in HepMC) is implemented.
637 // Thanks Sho Iwamoto for feedback
638 if(tau_plus)
639 tau_plus->decayEndgame();
640 if(tau_minus)
641 tau_minus->decayEndgame();
642
643}
644
645/***********************************************************
646 ** Below are the methods used for boosting and rotating
647 ** into another reference frame.
648 ************************************************************/
649
650/** Step 1. (Transformation A). Any modification to this method also requires a modification
651 to the inverse method boostFromTauPairFrameToLab (transformation A^-1).*/
653 double * rotation_angle2,
654 double * rotation_angle3,
655 TauolaParticle * mother,
656 vector<TauolaParticle *> grandmothers,
657 vector<TauolaParticle *> taus
658 ){ //in positive z axis (+1) //or negative z axis (-1)
659
660 /** boost all gradmothers and daughters (taus, neutrinos, etc,)
661 to the mothers rest frame */
662
663 for(int i=0; i< (int) grandmothers.size(); i++)
664 grandmothers.at(i)->boostToRestFrame(mother);
665 //If taus.size()==1 then taus[1] is the same as mum. Disaster to be avoided.
666 if(taus.size()!=1)
667 for(int i=0; i< (int) taus.size(); i++){
668 taus.at(i)->boostToRestFrame(mother);
669 // NOTE: Following line has no effect in release examples
670 // because taus don't have daughters at this step
671 taus.at(i)->boostDaughtersToRestFrame(mother);
672 }
673
674 /** rotate all particles so taus are on the z axis */
675 if(getTauPlus(taus)){ //if there's a tau+ use this to get the rotation angle
676 *rotation_angle1 = getTauPlus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
677 rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
678 *rotation_angle2 = getTauPlus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
679 rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
680 }
681 else{ //otherwise use the tau-
682 *rotation_angle1 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
683 rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
684 *rotation_angle2 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
685 rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
686 }
687
688 //now rotate incoming beams so they are aligned in the y-z plane
689 if(grandmothers.size()<1){ //if there are no grandmothers
690 *rotation_angle3 = 0;
691 return;
692 }
693
694 //the first grandmother will have a positive y component
695 if(getGrandmotherPlus(grandmothers))
697 else
699
700 rotateSystem(grandmothers,taus,*rotation_angle3,TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
701
702}
703
704/** Reverses boostFromLabtoMotherFrame. The three rotation angle must be provided.
705 Any modification to this would require a modification to boostFromLabToTauPairFrame
706 since this is the inverse transformation (Step 2: A^-1). */
708 double rotation_angle2,
709 double rotation_angle3,
710 TauolaParticle * mother,
711 vector<TauolaParticle *> grandmothers,
712 vector<TauolaParticle *> taus){
713
714
715 if(mother==0) //if there are no mothers
716 return;
717
718
719 //rotate grand mothers back away from the y-z plane
720 rotateSystem(grandmothers,taus,-rotation_angle3, TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
721
722 //rotate all so that taus are no longer on the z axis
723 rotateSystem(grandmothers,taus,-rotation_angle2,TauolaParticle::X_AXIS);
724 rotateSystem(grandmothers,taus,-rotation_angle1,TauolaParticle::Y_AXIS);
725
726
727 /*** boost grandmothers and daughters (taus) back into the lab frame */
728 for(int i=0; i< (int) grandmothers.size(); i++)
729 grandmothers.at(i)->boostFromRestFrame(mother);
730 //If taus.size()==1 then taus[1] is the same as mum. Disaster to be avoided.
731 if(taus.size()!=1)
732 for(int i=0; i< (int) taus.size(); i++){
733 taus.at(i)->boostFromRestFrame(mother);
734 taus.at(i)->boostDaughtersFromRestFrame(mother);
735 }
736}
737
738void TauolaParticlePair::rotateSystem(vector<TauolaParticle *> grandmothers,
739 vector<TauolaParticle *> taus,
740 double theta, int axis, int axis2){
741 for(int i=0; i< (int) grandmothers.size(); i++)
742 grandmothers.at(i)->rotate(axis,theta,axis2);
743 for(int i=0; i< (int) taus.size(); i++){
744 taus.at(i)->rotate(axis,theta,axis2);
745 taus.at(i)->rotateDaughters(axis,theta,axis2);
746 }
747}
748
749
750
751/***********************************************************
752 Some extra useful methods
753************************************************************/
755
756 //calculate mothers 4-momentum based on the daughters.
757 double e=0;
758 double px=0;
759 double py=0;
760 double pz=0;
761
762 bool tau_plus = false;
763 bool tau_minus = false;
764 bool tau_neutrino = false;
765 bool tau_antineutrino = false;
766
767 for(int d=0; d < (int) taus.size(); d++){
768 TauolaParticle * daughter = taus.at(d);
769 e+=daughter->getE();
770 px+=daughter->getPx();
771 py+=daughter->getPy();
772 pz+=daughter->getPz();
773 switch(daughter->getPdgID()){
775 tau_plus=true;
776 break;
778 tau_minus=true;
779 break;
781 tau_neutrino=true;
782 break;
784 tau_antineutrino=true;
785 break;
786 }
787 }
788 double m = e*e-px*px-py*py-pz*pz;
789 if(m<0)
790 m= -sqrt(-m);
791 else
792 m=sqrt(m);
793
794 //look for mothers type:
795 int pdg=0;
796
797 //Assume Z
798 if(tau_plus&&tau_minus)
800
801 //Assume W+
802 if(tau_plus&&tau_neutrino)
804
805 //Assume W-
806 if(tau_minus&&tau_antineutrino)
808
809 // now make the mother
810 return taus.at(0)->createNewParticle(pdg,2,m,px,py,pz,e);
811}
812
813// see if "particle" is one of the final taus in this tau pair
814// (based on particle barcode)
815// NOTE: Not executed by release examples
817
818 for(int i=0; i< (int) m_final_particles.size(); i++){
819 if(m_final_particles.at(i)->getBarcode()==particle->getBarcode())
820 return true;
821 }
822 return false;
823}
824
825TauolaParticle * TauolaParticlePair::getTauMinus(vector<TauolaParticle*> taus){
826 for(int i=0; i< (int) taus.size(); i++){
827 if(taus.at(i)->getPdgID()==TauolaParticle::TAU_MINUS)
828 return taus.at(i);
829 }
830 return 0;
831}
832
833TauolaParticle * TauolaParticlePair::getTauPlus(vector<TauolaParticle*> taus){
834 for(int i=0; i< (int) taus.size(); i++){
835 if(taus.at(i)->getPdgID()==TauolaParticle::TAU_PLUS)
836 return taus.at(i);
837 }
838 return 0;
839}
840
841TauolaParticle * TauolaParticlePair::getGrandmotherPlus(vector<TauolaParticle*> grandmothers){
842 //choose based on energy if there are more than one??
843 double e = -1e30;
844 int index = -1;
845 for(int i=0; i< (int) grandmothers.size(); i++){
846 if( (grandmothers.at(i)->getPdgID()<0 && grandmothers.at(i)->getPdgID()>-9) ||
847 grandmothers.at(i)->getPdgID()== TauolaParticle::POSITRON||
848 grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_PLUS){
849 if(e<grandmothers.at(i)->getE()){
850 e=grandmothers.at(i)->getE();
851 index=i;
852 }
853 }
854 }
855 if(index==-1)
856 {
857 for(int i=0; i< (int) grandmothers.size(); i++)
858 {
859 if(grandmothers.at(i)->getPdgID()==21 || grandmothers.at(i)->getPdgID()==22)
860 {
861 index=i;
862 e=grandmothers.at(i)->getE();
863 break;
864 }
865 }
866 }
867 if(index==-1) index=0;
868 if(e==0)
869 {
870 grandmothers.at(index)->print();
871 return 0;
872 }
873 else
874 return grandmothers.at(index);
875}
876
877TauolaParticle * TauolaParticlePair::getGrandmotherMinus(vector<TauolaParticle*> grandmothers){
878 //choose based on energy if there are more than one??
879 double e = -1e30;
880 int index = -1;
881 for(int i=0; i< (int) grandmothers.size(); i++){
882 if( (grandmothers.at(i)->getPdgID()>0 && grandmothers.at(i)->getPdgID()<9) ||
883 grandmothers.at(i)->getPdgID()== TauolaParticle::ELECTRON||
884 grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_MINUS){
885 if(e<grandmothers.at(i)->getE()){
886 e=grandmothers.at(i)->getE();
887 index=i;
888 }
889 }
890 }
891 if(index==-1)
892 {
893 for(int i=(int) grandmothers.size()-1; i>=0 ; i--)
894 {
895 if(grandmothers.at(i)->getPdgID()==21||grandmothers.at(i)->getPdgID()==22)
896 {
897 index=i;
898 e=grandmothers.at(i)->getE();
899 break;
900 }
901 }
902 }
903 if(index==-1) index=0;
904 if(e==0)
905 return 0;
906 else
907 return grandmothers.at(index);
908}
909
910
911
913
914 Log::RedirectOutput(Log::Info());
915 std::cout << "Daughters final:" << std::endl;
916 for(int i=0; i< (int) m_final_particles.size(); i++)
917 m_final_particles.at(i)->print();
918
919
920 std::cout << "Daughters at production:" << std::endl;
921 for(int i=0; i< (int) m_production_particles.size(); i++)
922 m_production_particles.at(i)->print();
923
924
925 std::cout << "Mother particle: " << std::endl;
926 if(m_mother)
927 m_mother->print();
928
929 std::cout << "Grandmother particles: " << std::endl;
930 for(int i=0; i< (int) m_grandmothers.size(); i++)
931 m_grandmothers.at(i)->print();
932
934}
935
936
938
939 for(int i=0; i< (int) m_final_particles.size(); i++)
940 m_final_particles.at(i)->checkMomentumConservation();
941
942 for(int i=0; i< (int) m_production_particles.size(); i++)
943 m_production_particles.at(i)->checkMomentumConservation();
944
945 if(m_mother)
947
948 for(int i=0; i< (int) m_grandmothers.size(); i++)
949 m_grandmothers.at(i)->checkMomentumConservation();
950
951}
952
953void TauolaParticlePair::recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta){
954
955 if (abs(outgoing_pdg_id)!=15)
956 {
957 Log::Warning()<<"interface was not used for taus pdg id="<<outgoing_pdg_id<<endl;
958 return;
959 }
960
961 double (*T) [Tauola::NCOS][4][4] = NULL;
962 double (*Tw) [Tauola::NCOS] = NULL;
963 double (*Tw0)[Tauola::NCOS] = NULL;
964 double smin = 0.0; // Tauola::sminA;
965 double smax = 0.0; // Tauola::smaxA;
966 double step = 0.0; // (smaxl-sminl)/(Tauola::NS1-1);
967
968 // Select table for appropriate incoming particles
969 switch(abs(incoming_pdg_id)){
970 case 11:
971 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
972 {
973 T = Tauola::table11B;
974 Tw = Tauola::wtable11B;
975 Tw0 = Tauola::w0table11B;
976 smin = Tauola::sminB;
977 smax = Tauola::smaxB;
978 step = (smax-smin)/(Tauola::NS2-1);
979 }
980 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
981 {
982 T = Tauola::table11C;
983 Tw = Tauola::wtable11C;
984 Tw0 = Tauola::w0table11C;
985 smin = Tauola::sminC;
986 smax = Tauola::smaxC;
987 step = (smax-smin)/(Tauola::NS3-1);
988 }
989 else
990 {
991 T = Tauola::table11A;
992 Tw = Tauola::wtable11A;
993 Tw0 = Tauola::w0table11A;
994 smin = Tauola::sminA;
995 smax = Tauola::smaxA;
996 step = (smax-smin)/(Tauola::NS1-1);
997 }
998 break;
999
1000 // NOTE: Use the same tables for 1, 3 and 5
1001 case 1:
1002 case 3:
1003 case 5:
1004 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
1005 {
1006 T = Tauola::table1B;
1007 Tw = Tauola::wtable1B;
1008 Tw0 = Tauola::w0table1B;
1009 smin = Tauola::sminB;
1010 smax = Tauola::smaxB;
1011 step = (smax-smin)/(Tauola::NS2-1);
1012 }
1013 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
1014 {
1015 T = Tauola::table1C;
1016 Tw = Tauola::wtable1C;
1017 Tw0 = Tauola::w0table1C;
1018 smin = Tauola::sminC;
1019 smax = Tauola::smaxC;
1020 step = (smax-smin)/(Tauola::NS3-1);
1021 }
1022 else
1023 {
1024 T = Tauola::table1A;
1025 Tw = Tauola::wtable1A;
1026 Tw0 = Tauola::w0table1A;
1027 smin = Tauola::sminA;
1028 smax = Tauola::smaxA;
1029 step = (smax-smin)/(Tauola::NS1-1);
1030 }
1031 break;
1032
1033 // NOTE: Use the same tables for 2 and 4
1034 case 2:
1035 case 4:
1036 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
1037 {
1038 T = Tauola::table2B;
1039 Tw = Tauola::wtable2B;
1040 Tw0 = Tauola::w0table2B;
1041 smin = Tauola::sminB;
1042 smax = Tauola::smaxB;
1043 step = (smax-smin)/(Tauola::NS2-1);
1044 }
1045 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
1046 {
1047 T = Tauola::table2C;
1048 Tw = Tauola::wtable2C;
1049 Tw0 = Tauola::w0table2C;
1050 smin = Tauola::sminC;
1051 smax = Tauola::smaxC;
1052 step = (smax-smin)/(Tauola::NS3-1);
1053 }
1054 else
1055 {
1056 T = Tauola::table2A;
1057 Tw = Tauola::wtable2A;
1058 Tw0 = Tauola::w0table2A;
1059 smin = Tauola::sminA;
1060 smax = Tauola::smaxA;
1061 step = (smax-smin)/(Tauola::NS1-1);
1062 }
1063 break;
1064
1065 default:
1066 Log::Warning()<<"interface was not used for proper beams pdg id="<<incoming_pdg_id<<endl;
1067 return;
1068 }
1069
1070 // Check if we have found a table for matrix recalculation
1071 if (T[0][0][0][0]<0.5) return;
1072
1073 // If mass is too small
1074 if (invariant_mass_squared <= exp(Tauola::sminA)){
1075 double mta = 1.777;
1076 double cosf = cosTheta;
1077 double s = invariant_mass_squared;
1078 double sinf = sqrt(1-cosf*cosf);
1079 double MM = sqrt(4e0*mta*mta/s);
1080 double xnorm = 1+cosf*cosf + MM*MM*sinf*sinf;
1081
1082 // Zero matrix
1083 for (int i=0;i<4;i++)
1084 for (int j=0;j<4;j++)
1085 m_R[i][j]=0.0;
1086
1087 m_R[0][0] = (1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
1088 m_R[1][1] = (-(1- MM*MM)*sinf*sinf)/xnorm;
1089 m_R[2][2] = ( (1+ MM*MM)*sinf*sinf)/xnorm;
1090 m_R[3][3] = (1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
1091 m_R[2][3] = (2*MM*sinf*cosf)/xnorm;
1092 m_R[3][2] = (2*MM*sinf*cosf)/xnorm;
1093
1094 // Get weights
1095 double w = 1.;
1096 double w0 = 1.;
1097
1098 if (Tauola::wtable2A[0][0]>0 ) w=Tauola::wtable2A[0][0];
1099 else if(Tauola::wtable1A[0][0]>0 ) w=Tauola::wtable1A[0][0];
1100 else if(Tauola::wtable11A[0][0]>0) w=Tauola::wtable11A[0][0];
1101
1102 if (Tauola::wtable2A[0][0]>0 ) w0=Tauola::w0table2A[0][0];
1103 else if(Tauola::wtable1A[0][0]>0 ) w0=Tauola::w0table1A[0][0];
1104 else if(Tauola::wtable11A[0][0]>0) w0=Tauola::w0table11A[0][0];
1105
1106 // Tauola::setEWwt(w/w0);
1107 Tauola::setEWwt(w,w0);
1108 return;
1109 } // if(mass too small)
1110
1111 int x = 0;
1112 double buf = smin;
1113 double remnant = 0.0;
1114
1115 // Interpolate s
1116 if(T==Tauola::table11A || T==Tauola::table1A || T== Tauola::table2A)
1117 {
1118 while(log(invariant_mass_squared) > buf){
1119 x++;
1120 buf += (step);
1121 }
1122 remnant = (log(invariant_mass_squared)-(buf-step))/step;
1123 }
1124 else
1125 {
1126 while(invariant_mass_squared > buf){
1127 x++;
1128 buf += step;
1129 }
1130 remnant = (invariant_mass_squared-(buf-step))/step;
1131 }
1132
1133 double cmin = -1.;
1134 int y = 0;
1135 double remnantc = 0.0;
1136
1137 // Interpolate cosTheta
1138 buf = cmin+1./Tauola::NCOS;
1139 while(cosTheta > buf){
1140 y++;
1141 buf+=2./Tauola::NCOS;
1142 }
1143
1144 remnantc = (cosTheta-(buf-2./Tauola::NCOS))/(2./Tauola::NCOS);
1145
1146 // Special cases at edges - extrapolation
1147 bool isUsingExtrapolation = false;
1148 if(y >= Tauola::NCOS) { isUsingExtrapolation = true; y = Tauola::NCOS-1; }
1149 if(y == 0) { isUsingExtrapolation = true; y = 1; }
1150
1151 // Bilinear interpolation
1152 double b11,b21,b12,b22;
1153
1154 for (int i=0; i<4; i++){
1155 for (int j=0; j<4; j++){
1156
1157 if(isUsingExtrapolation){
1158 if(y == 1){
1159 b11 = 2*T[x-1][0][i][j] - T[x-1][1][i][j];
1160 b21 = 2*T[x][0][i][j] - T[x][1][i][j];
1161 b12 = T[x-1][0][i][j];
1162 b22 = T[x][0][i][j];
1163 }
1164 else{
1165 b11 = T[x-1][y][i][j];
1166 b21 = T[x][y][i][j];
1167 b12 = 2*T[x-1][y][i][j] - T[x-1][y-1][i][j];
1168 b22 = 2*T[x][y][i][j] - T[x][y-1][i][j];
1169 }
1170 } // if(isUsingExtrapolation)
1171 else{
1172 b11 = T[x-1][y-1][i][j];
1173 b21 = T[x][y-1][i][j];
1174 b12 = T[x-1][y][i][j];
1175 b22 = T[x][y][i][j];
1176 }
1177 m_R[i][j] = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1178 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1179 } // for(j)
1180 } // for(i)
1181
1182 // Calculate electroweak weights
1183 double w,w0;
1184
1185 if(isUsingExtrapolation){
1186 if(y == 1){
1187 b11 = 2*Tw[x-1][0] - Tw[x-1][1];
1188 b21 = 2*Tw[x][0] - Tw[x][1];
1189 b12 = Tw[x-1][0];
1190 b22 = Tw[x][0];
1191 }
1192 else
1193 {
1194 b11 = Tw[x-1][y];
1195 b21 = Tw[x][y];
1196 b12 = 2*Tw[x-1][y] - Tw[x-1][y-1];
1197 b22 = 2*Tw[x][y] - Tw[x][y-1];
1198 }
1199 } // if(isUsingExtrapolation)
1200 else
1201 {
1202 b11 = Tw[x-1][y-1];
1203 b21 = Tw[x][y-1];
1204 b12 = Tw[x-1][y];
1205 b22 = Tw[x][y];
1206 }
1207
1208 w = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1209 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1210
1211 if(isUsingExtrapolation){
1212 if(y == 1){
1213 b11 = 2*Tw0[x-1][0] - Tw0[x-1][1];
1214 b21 = 2*Tw0[x][0] - Tw0[x][1];
1215 b12 = Tw0[x-1][0];
1216 b22 = Tw0[x][0];
1217 }
1218 else{
1219 b11 = Tw0[x-1][y];
1220 b21 = Tw0[x][y];
1221 b12 = 2*Tw0[x-1][y] - Tw0[x-1][y-1];
1222 b22 = 2*Tw0[x][y] - Tw0[x][y-1];
1223 }
1224 } // if (isUsingExtrapolation)
1225 else{
1226 b11 = Tw0[x-1][y-1];
1227 b21 = Tw0[x][y-1];
1228 b12 = Tw0[x-1][y];
1229 b22 = Tw0[x][y];
1230 }
1231
1232 w0 = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1233 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1234
1235 Tauola::setEWwt(w,w0);
1236}
1237
1238} // namespace Tauolapp
static void AddDecay(int type)
Definition Log.cxx:25
static void Fatal(string text, unsigned short int code=0)
static void RevertOutput()
static void RedirectOutput(void(*func)(), ostream &where=*out)
Definition Log.cxx:93
TauolaParticle * getGrandmotherMinus(std::vector< TauolaParticle * > particles)
TauolaParticle * getGrandmotherPlus(std::vector< TauolaParticle * > particles)
bool contains(TauolaParticle *particle)
void addToBeam(TauolaParticle *pcle, std::vector< TauolaParticle * > *candidates_same, std::vector< TauolaParticle * > *candidates_opp)
void rotateSystem(vector< TauolaParticle * > grandmothers, vector< TauolaParticle * > taus, double theta, int axis, int axis2=TauolaParticle::Z_AXIS)
void recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta)
TauolaParticle * makeTemporaryMother(vector< TauolaParticle * > taus)
TauolaParticle * getTauMinus(std::vector< TauolaParticle * > particles)
TauolaParticle * getTauPlus(std::vector< TauolaParticle * > particles)
std::vector< TauolaParticle * > m_grandmothers
static void setBornKinematics(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta)
double getVirtuality(TauolaParticle *p1, TauolaParticle *p2, bool flip)
void boostFromLabToTauPairFrame(double *rotation_angle1, double *rotation_angle2, double *rotation_angle3, TauolaParticle *mother, vector< TauolaParticle * > grandmothers, vector< TauolaParticle * > taus)
std::vector< TauolaParticle * > m_production_particles
std::vector< TauolaParticle * > m_final_particles
void boostFromTauPairToLabFrame(double rotation_angle1, double rotation_angle2, double rotation_angle3, TauolaParticle *mother, vector< TauolaParticle * > grandmothers, vector< TauolaParticle * > taus)
double getZPolarization(int *incoming_pdg_id, int *outgoing_pdg_id, double *invMass, double *cosTheta)
virtual void print()=0
virtual int getBarcode()=0
virtual int getPdgID()=0
virtual double getPz()=0
TauolaParticle * clone()
void subtract(TauolaParticle *)
std::vector< TauolaParticle * > findProductionMothers()
virtual double getE()=0
virtual double getPy()=0
virtual double getPx()=0
virtual void setPdgID(int pdg_id)=0
double getRotationAngle(int axis, int second_axis=Z_AXIS)
TauolaParticle * findLastSelf()
void add(TauolaParticle *)
static double getTauMass()
Definition Tauola.cxx:668
static int getHiggsScalarPseudoscalarPDG()
Definition Tauola.cxx:676
static bool isUsingDecayOneBoost()
Definition Tauola.cxx:586
static bool isUsingDecayOne()
Definition Tauola.cxx:581
static const double * getDecayOnePolarization()
Definition Tauola.cxx:601