00001 #include "Tauola.h"
00002 #include "TauolaParticlePair.h"
00003 #include "Log.h"
00004 #include <stdlib.h>
00005 #include <math.h>
00006 #include <iostream>
00007
00008 namespace Tauolapp
00009 {
00010
00011
00012 TauolaParticlePair::TauolaParticlePair(std::vector<TauolaParticle*> &particle_list){
00013 TauolaParticle *particle=particle_list.back();
00014 particle_list.pop_back();
00015 setBornKinematics(0,0,-1.,0.);
00016
00017
00018 if(Tauola::isUsingDecayOne())
00019 {
00020 m_mother = 0;
00021 m_mother_exists = false;
00022 m_production_particles.push_back(particle);
00023 m_final_particles.push_back(particle);
00024 initializeDensityMatrix();
00025 Log::AddDecay(0);
00026 return;
00027 }
00028
00029
00030 std::vector<TauolaParticle *> temp_mothers = particle->findProductionMothers();
00031
00032
00033 if(temp_mothers.size()==0){
00034
00035
00036 Log::Warning()<< "WARNING: Could not find taus mother or grandmothers. "
00037 << "Ignoring spin effects" << std::endl;
00038 m_mother = 0;
00039 m_mother_exists = false;
00040 m_production_particles.push_back(particle);
00041 m_final_particles.push_back(particle->findLastSelf());
00042 initializeDensityMatrix();
00043 Log::AddDecay(1);
00044 return;
00045 }
00046
00047
00048 std::vector<TauolaParticle*> temp_daughters;
00049 temp_daughters = temp_mothers.at(0)->getDaughters();
00050 if(temp_daughters.size()==0)
00051 Log::Fatal("WARNING: Something wrong with event structure or there is a bug in the TAUOLA interface.",6);
00052
00053 m_production_particles=temp_daughters;
00054 m_final_particles.push_back(particle->findLastSelf());
00055
00056
00057
00058
00059 for(signed int i=0; i < (int) m_production_particles.size(); i++)
00060 {
00061
00062 if(m_final_particles.at(0)->getPdgID()*m_production_particles.at(i)->getPdgID()>=0) continue;
00063
00064 if(m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_PLUS||
00065 m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_MINUS)
00066 {
00067
00068 int j=-1;
00069 for(j=0;j<(int)particle_list.size();j++)
00070 if(m_production_particles.at(i)->getBarcode()==particle_list.at(j)->getBarcode()) break;
00071 if(j>=(int)particle_list.size()) continue;
00072
00073
00074 m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
00075 particle_list.erase(particle_list.begin()+j);
00076 }
00077 else if(m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_NEUTRINO||
00078 m_production_particles.at(i)->getPdgID()==TauolaParticle::TAU_ANTINEUTRINO)
00079 {
00080
00081 m_final_particles.push_back(m_production_particles.at(i)->findLastSelf());
00082 }
00083 }
00084
00085
00086
00087 if(temp_mothers.size()==1){
00088 m_mother_exists = true;
00089 m_mother = temp_mothers.at(0);
00090 m_grandmothers = m_mother->findProductionMothers();
00091 Log::AddDecay(3);
00092 }
00093 else{
00094 m_mother_exists = false;
00095 m_grandmothers = temp_mothers;
00096 m_mother = makeTemporaryMother(m_production_particles);
00097 Log::AddDecay(2);
00098 }
00099
00100 initializeDensityMatrix();
00101
00102 return;
00103 }
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 void TauolaParticlePair::initializeDensityMatrix(){
00114 int incoming_pdg_id=0;
00115 int outgoing_pdg_id=0;
00116 double invariant_mass_squared=-5.0;
00117 double cosTheta=3.0;
00118
00119
00120 for(int x = 0; x < 4; x ++) {
00121 for(int y = 0; y < 4; y ++)
00122 m_R[x][y] = 0;
00123 }
00124
00125 m_R[0][0]=1;
00126
00127 if(Tauola::isUsingDecayOne())
00128 {
00129 const double *pol = Tauola::getDecayOnePolarization();
00130
00131 m_R[0][1]=pol[0];
00132 m_R[0][2]=pol[1];
00133 m_R[0][3]=pol[2];
00134
00135 m_R[1][0]=pol[0];
00136 m_R[2][0]=pol[1];
00137 m_R[3][0]=pol[2];
00138 }
00139
00140 if(!m_mother) return;
00141
00142
00143
00144
00145
00146
00147
00148 if(m_mother->getPdgID()==Tauola::getHiggsScalarPseudoscalarPDG()){
00149 if(!Tauola::spin_correlation.HIGGS_H) return;
00150
00151 double phi = Tauola::getHiggsScalarPseudoscalarMixingAngle();
00152 double mass_ratio = Tauola::getTauMass()/m_mother->getMass();
00153
00154 double beta = sqrt(1-4*mass_ratio*mass_ratio);
00155 double denominator = pow(cos(phi)*beta,2)+pow(sin(phi),2);
00156
00157 m_R[0][0]= 1;
00158 m_R[1][1]= (pow(cos(phi)*beta,2)-pow(sin(phi),2))/denominator;
00159 m_R[1][2]= 2*cos(phi)*sin(phi)*beta/denominator;
00160 m_R[2][1]= -m_R[1][2];
00161 m_R[2][2]= m_R[1][1];
00162 m_R[3][3]= -1;
00163 }
00164 else {
00165
00166 double pz = 0.0;
00167
00168 switch(m_mother->getPdgID()){
00169
00170 case TauolaParticle::Z0:
00171 if(!Tauola::spin_correlation.Z0) break;
00172
00173
00174
00175
00176
00177 pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
00178 m_R[0][0]=1;
00179 m_R[0][3]=2*pz-1;
00180 m_R[3][0]=2*pz-1;
00181 m_R[3][3]=1;
00182
00183 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
00184 break;
00185
00186 case TauolaParticle::GAMMA:
00187 if(!Tauola::spin_correlation.GAMMA) break;
00188
00189
00190
00191
00192
00193 pz = getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
00194 m_R[0][0]=1;
00195 m_R[0][3]=2*pz-1;
00196 m_R[3][0]=2*pz-1;
00197 m_R[3][3]=1;
00198
00199 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
00200 break;
00201
00202
00203 case TauolaParticle::HIGGS:
00204 if(!Tauola::spin_correlation.HIGGS) break;
00205 m_R[0][0]=1;
00206 m_R[1][1]=1;
00207 m_R[2][2]=1;
00208 m_R[3][3]=-1;
00209 break;
00210
00211
00212 case TauolaParticle::HIGGS_A:
00213 if(!Tauola::spin_correlation.HIGGS_A) break;
00214 m_R[0][0]=1;
00215 m_R[1][1]=-1;
00216 m_R[2][2]=-1;
00217 m_R[3][3]=-1;
00218 break;
00219
00220 case TauolaParticle::HIGGS_PLUS:
00221 if(!Tauola::spin_correlation.HIGGS_PLUS) break;
00222 m_R[0][0]=1;
00223 m_R[0][3]=-1;
00224 m_R[3][0]=-1;
00225 break;
00226
00227 case TauolaParticle::HIGGS_MINUS:
00228 if(!Tauola::spin_correlation.HIGGS_MINUS) break;
00229 m_R[0][0]=1;
00230 m_R[0][3]=-1;
00231 m_R[3][0]=-1;
00232 break;
00233
00234 case TauolaParticle::W_PLUS:
00235 if(!Tauola::spin_correlation.W_PLUS) break;
00236 m_R[0][0]=1;
00237 m_R[0][3]=1;
00238 m_R[3][0]=1;
00239 break;
00240
00241 case TauolaParticle::W_MINUS:
00242 if(!Tauola::spin_correlation.W_MINUS) break;
00243 m_R[0][0]=1;
00244 m_R[0][3]=1;
00245 m_R[3][0]=1;
00246 break;
00247
00248
00249 default:
00250 m_R[0][0]=1;
00251 break;
00252 }
00253 }
00254
00255 }
00256
00257
00258 void TauolaParticlePair::setBornKinematics(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared,double cosTheta){
00259 Tauola::buf_incoming_pdg_id=incoming_pdg_id;
00260 Tauola::buf_outgoing_pdg_id=outgoing_pdg_id;
00261 Tauola::buf_invariant_mass_squared=invariant_mass_squared;
00262 Tauola::buf_cosTheta=cosTheta;
00263
00264
00265 }
00266
00267
00268 double TauolaParticlePair::getZPolarization(int *incoming_pdg_id, int *outgoing_pdg_id, double *invariant_mass_squared,double *cosTheta){
00269
00270
00271 *incoming_pdg_id = TauolaParticle::ELECTRON;
00272 *cosTheta = 0 ;
00273 *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
00274 *invariant_mass_squared = pow(m_mother->getMass(),2);
00275 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
00276
00277
00278
00279
00280 if(m_grandmothers.size()<2){
00281 Log::Warning()<<"Not enough mothers of Z to "
00282 <<"calculate cos(theta) between "
00283 <<"incoming about outgoing beam"
00284 << endl;
00285 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
00286 invariant_mass_squared, cosTheta);
00287 }
00288
00289 if(!getTauPlus(m_production_particles)||!getTauMinus(m_production_particles)){
00290 Log::Error()<<"tau+ or tau- not found in Z decay"<< endl;
00291 return 0;
00292 }
00293
00294
00295
00296
00297
00298 if(m_grandmothers.size()>2 ||
00299 (m_grandmothers.at(0)->getDaughters().size()>1 && m_mother_exists==true) ||
00300 m_production_particles.size() > 2){
00301
00302
00303 vector<TauolaParticle*> extra_grandmothers;
00304 for(int i=0; i<(int) m_grandmothers.size(); i++){
00305
00306 if(m_grandmothers.at(i)!=getGrandmotherPlus(m_grandmothers)&&
00307 m_grandmothers.at(i)!=getGrandmotherMinus(m_grandmothers))
00308 extra_grandmothers.push_back(m_grandmothers.at(i));
00309 }
00310
00311
00312
00313 vector<TauolaParticle*> extra_tau_siblings;
00314 vector<TauolaParticle*> temp_production_particles;
00315 for(int i=0; i<(int) m_production_particles.size(); i++){
00316 if(m_production_particles.at(i)!=getTauPlus(m_production_particles)&&
00317 m_production_particles.at(i)!=getTauMinus(m_production_particles))
00318 extra_tau_siblings.push_back(m_production_particles.at(i));
00319 }
00320
00321
00322 vector<TauolaParticle*> extra_Z_siblings;
00323 for(int i=0; m_mother_exists && i<(int) m_grandmothers.at(0)->getDaughters().size(); i++){
00324 if(m_grandmothers.at(0)->getDaughters().at(i)->getPdgID()!=TauolaParticle::Z0)
00325 extra_Z_siblings.push_back(m_grandmothers.at(0)->getDaughters().at(i));
00326 }
00327
00328
00329
00330 std::vector<TauolaParticle*> effective_taus;
00331 effective_taus.push_back(getTauPlus(m_production_particles)->clone());
00332 effective_taus.push_back(getTauMinus(m_production_particles)->clone());
00333
00334
00335 TauolaParticle * g1 = getGrandmotherPlus(m_grandmothers)->clone();
00336 TauolaParticle * g2 = getGrandmotherMinus(m_grandmothers)->clone();
00337 m_grandmothers.clear();
00338 m_grandmothers.push_back(g1);
00339 m_grandmothers.push_back(g2);
00340
00341
00342 for(int i=0; i<(int) extra_grandmothers.size(); i++)
00343 addToBeam(extra_grandmothers.at(i),&m_grandmothers,0);
00344
00345
00346 for(int i=0; i<(int) extra_Z_siblings.size(); i++)
00347 addToBeam(extra_Z_siblings.at(i),0,&m_grandmothers);
00348
00349
00350 for(int i=0; i<(int) extra_tau_siblings.size() ; i++){
00351 if(m_mother_exists)
00352 addToBeam(extra_tau_siblings.at(i),&effective_taus,0);
00353 else
00354 addToBeam(extra_tau_siblings.at(i),&effective_taus,&m_grandmothers);
00355 }
00356
00357
00358 TauolaParticle * temp_mother = makeTemporaryMother(effective_taus);
00359 *invariant_mass_squared = pow(temp_mother->getMass(),2);
00360
00361
00362
00363
00364 double angle1,angle2,angle3;
00365 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3, temp_mother,
00366 m_grandmothers, effective_taus);
00367
00368
00369
00370
00371
00372 TauolaParticle *tM=getTauPlus(effective_taus);
00373 TauolaParticle *gM=getGrandmotherPlus(m_grandmothers);
00374
00375 double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
00376 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
00377 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
00378 tM=getTauMinus(effective_taus);
00379 gM=getGrandmotherMinus(m_grandmothers);
00380
00381 double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
00382 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
00383 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
00384 double sintheta1 = sqrt(1-costheta1*costheta1);
00385 double sintheta2 = sqrt(1-costheta2*costheta2);
00386 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
00387
00388 *cosTheta = avg;
00389
00390 *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
00391
00392 if(abs(*incoming_pdg_id)==21 || abs(*incoming_pdg_id)==22)
00393 {
00394 *incoming_pdg_id = -getGrandmotherMinus(m_grandmothers)->getPdgID();
00395
00396 }
00397
00398 if( abs(*incoming_pdg_id)> 8 &&
00399 abs(*incoming_pdg_id)!=11 &&
00400 abs(*incoming_pdg_id)!=13 &&
00401 abs(*incoming_pdg_id)!=15 )
00402 {
00403 Log::Error()<<"Second class disaster: incoming_pdg_id = "<<*incoming_pdg_id<<endl;
00404
00405
00406 *incoming_pdg_id = TauolaParticle::ELECTRON;
00407 *cosTheta = 0 ;
00408 *outgoing_pdg_id = TauolaParticle::TAU_PLUS;
00409
00410 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
00411 invariant_mass_squared, cosTheta);
00412 }
00413
00414 boostFromTauPairToLabFrame(angle1, angle2, angle3, temp_mother,
00415 m_grandmothers, effective_taus);
00416 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
00417 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
00418
00419 }
00420
00421
00422
00423
00424 double angle1,angle2,angle3;
00425 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
00426 m_mother,m_grandmothers,m_production_particles);
00427
00428 TauolaParticle *tM=getTauPlus(m_production_particles);
00429 TauolaParticle *gM=getGrandmotherPlus(m_grandmothers);
00430 double costheta1=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
00431 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
00432 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
00433
00434 tM=getTauMinus(m_production_particles);
00435 gM=getGrandmotherMinus(m_grandmothers);
00436
00437 double costheta2=(tM->getPx()*gM->getPx()+tM->getPy()*gM->getPy()+tM->getPz()*gM->getPz())/
00438 sqrt(tM->getPx()*tM->getPx()+tM->getPy()*tM->getPy()+tM->getPz()*tM->getPz())/
00439 sqrt(gM->getPx()*gM->getPx()+gM->getPy()*gM->getPy()+gM->getPz()*gM->getPz());
00440
00441 double sintheta1 = sqrt(1-costheta1*costheta1);
00442 double sintheta2 = sqrt(1-costheta2*costheta2);
00443
00444 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
00445
00446 *cosTheta = avg;
00447
00448 *incoming_pdg_id = getGrandmotherPlus(m_grandmothers)->getPdgID();
00449
00450 boostFromTauPairToLabFrame(angle1, angle2, angle3,
00451 m_mother,m_grandmothers,m_production_particles);
00452
00453 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
00454 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
00455
00456 }
00457
00458
00459
00460
00461
00462
00463 void TauolaParticlePair::addToBeam(TauolaParticle * pcle,
00464 std::vector<TauolaParticle*> *candidates_same,
00465 std::vector<TauolaParticle*> *candidates_opp){
00466
00467
00468 double s=0.,o=-10.;
00469 TauolaParticle * s_beam = NULL;
00470 TauolaParticle * o_beam = NULL;
00471
00472 if(candidates_same){
00473 double s0 =1.0/getVirtuality(pcle,candidates_same->at(0),false);
00474 double s1 = 1.0/getVirtuality(pcle,candidates_same->at(1),false);
00475 if(s0>s1){
00476 s=s0;
00477 s_beam = candidates_same->at(0);
00478 }
00479 else{
00480 s=s1;
00481 s_beam = candidates_same->at(1);
00482 }
00483 }
00484 if(candidates_opp){
00485
00486 double o0 =1.0/getVirtuality(pcle,candidates_opp->at(0),true);
00487 double o1 = 1.0/getVirtuality(pcle,candidates_opp->at(1),true);
00488 if(o0>o1){
00489 o=o0;
00490 o_beam = candidates_opp->at(0);
00491 }
00492 else{
00493 o=o1;
00494 o_beam = candidates_opp->at(1);
00495 }
00496 }
00497
00498 if(s>o)
00499 {
00500 s_beam->add(pcle);
00501 int pdg1 = pcle->getPdgID();
00502 int pdg2 = s_beam->getPdgID();
00503 if(abs(pdg2)==15) return;
00504 if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) s_beam->setPdgID( pdg1);
00505 }
00506 else
00507 {
00508 o_beam->subtract(pcle);
00509 int pdg1 = pcle->getPdgID();
00510 int pdg2 = o_beam->getPdgID();
00511 if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) o_beam->setPdgID(-pdg1);
00512 }
00513
00514
00515
00516 }
00517
00518 double TauolaParticlePair::getVirtuality(TauolaParticle * p1, TauolaParticle*p2, bool flip){
00519
00520
00521 if((p1->getPdgID()==TauolaParticle::GLUON&&abs(p2->getPdgID())>10&&abs(p2->getPdgID())<20) ||
00522 (p2->getPdgID()==TauolaParticle::GLUON&&abs(p1->getPdgID())>10&&abs(p1->getPdgID())<20))
00523 return -1;
00524
00525
00526
00527
00528 double kp = p1->getE()*p2->getE() - p1->getPx()*p2->getPx()
00529 - p1->getPy()*p2->getPy() - p1->getPz()*p2->getPz();
00530 if(flip)
00531 kp = kp - p1->getMass()*p1->getMass();
00532
00533 double q = 1;
00534 if(p1->getPdgID()==TauolaParticle::GAMMA){
00535 if(abs(p2->getPdgID())==TauolaParticle::UP ||
00536 abs(p2->getPdgID())==TauolaParticle::CHARM ||
00537 abs(p2->getPdgID())==TauolaParticle::TOP)
00538 q=2.0/3.0;
00539 if(abs(p2->getPdgID())==TauolaParticle::DOWN ||
00540 abs(p2->getPdgID())==TauolaParticle::STRANGE ||
00541 abs(p2->getPdgID())==TauolaParticle::BOTTOM)
00542 q=1.0/3.0;
00543 }
00544 if(p2->getPdgID()==TauolaParticle::GAMMA){
00545 if(abs(p1->getPdgID())==TauolaParticle::UP ||
00546 abs(p1->getPdgID())==TauolaParticle::CHARM ||
00547 abs(p1->getPdgID())==TauolaParticle::TOP)
00548 q=2.0/3.0;
00549 if(abs(p1->getPdgID())==TauolaParticle::DOWN ||
00550 abs(p1->getPdgID())==TauolaParticle::STRANGE ||
00551 abs(p1->getPdgID())==TauolaParticle::BOTTOM)
00552 q=1.0/3.0;
00553 }
00554
00555 return fabs(2*kp)/q;
00556 }
00557
00558
00559
00560
00561
00562 void TauolaParticlePair::decayTauPair(){
00563
00564
00565 double h_tau_minus[4]={2,0,0,0};
00566 double h_tau_plus[4]={2,0,0,0};
00567
00568 TauolaParticle * tau_minus = getTauMinus(m_final_particles);
00569 TauolaParticle * tau_plus = getTauPlus(m_final_particles);
00570
00571
00572 for(double weight=0; weight < Tauola::randomDouble();){
00573
00574 if(tau_minus){
00575 Tauola::redefineTauMinusProperties(tau_minus);
00576 tau_minus->decay();
00577 h_tau_minus[0]=1;
00578 h_tau_minus[1]=tau_minus->getPolarimetricX();
00579 h_tau_minus[2]=tau_minus->getPolarimetricY();
00580 h_tau_minus[3]=tau_minus->getPolarimetricZ();
00581 }
00582 if(tau_plus){
00583 Tauola::redefineTauPlusProperties(tau_plus);
00584 tau_plus->decay();
00585 h_tau_plus[0]=1;
00586 h_tau_plus[1]=tau_plus->getPolarimetricX();
00587 h_tau_plus[2]=tau_plus->getPolarimetricY();
00588 h_tau_plus[3]=tau_plus->getPolarimetricZ();
00589 }
00590
00591
00592
00593 weight=0;
00594 for(int i =0; i<4; i++){
00595 for(int j=0; j<4; j++)
00596 weight+=m_R[i][j]*h_tau_minus[i]*h_tau_plus[j];
00597 }
00598 weight = weight/4.0;
00599 }
00600 double wthel[4]={0};
00601 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]);
00602 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]);
00603 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]);
00604 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]);
00605
00606 double rn=Tauola::randomDouble();
00607 if (rn>(wthel[0]+wthel[1]+wthel[2])/(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,-1);
00608 else if (rn>(wthel[0]+wthel[1]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities(-1,+1);
00609 else if (rn>(wthel[0]) /(wthel[0]+wthel[1]+wthel[2]+wthel[3])) Tauola::setHelicities( 1,-1);
00610 else Tauola::setHelicities( 1, 1);
00611
00612
00613
00614
00615
00616
00617
00618 double angle1,angle2,angle3;
00619 TauolaParticle * mum = makeTemporaryMother(m_final_particles);
00620
00621 if(!Tauola::isUsingDecayOneBoost())
00622 boostFromLabToTauPairFrame(&angle1, &angle2, &angle3,
00623 mum,m_grandmothers,m_final_particles);
00624
00625
00626 if(tau_plus)
00627 tau_plus->addDecayToEventRecord();
00628 if(tau_minus)
00629 tau_minus->addDecayToEventRecord();
00630
00631 if(!Tauola::isUsingDecayOneBoost())
00632 boostFromTauPairToLabFrame(angle1,angle2,angle3,
00633 mum,m_grandmothers,m_final_particles);
00634
00635
00636
00637
00638 if(tau_plus)
00639 tau_plus->decayEndgame();
00640 if(tau_minus)
00641 tau_minus->decayEndgame();
00642
00643 }
00644
00645
00646
00647
00648
00649
00650
00651
00652 void TauolaParticlePair::boostFromLabToTauPairFrame(double * rotation_angle1,
00653 double * rotation_angle2,
00654 double * rotation_angle3,
00655 TauolaParticle * mother,
00656 vector<TauolaParticle *> grandmothers,
00657 vector<TauolaParticle *> taus
00658 ){
00659
00660
00661
00662
00663 for(int i=0; i< (int) grandmothers.size(); i++)
00664 grandmothers.at(i)->boostToRestFrame(mother);
00665
00666 if(taus.size()!=1)
00667 for(int i=0; i< (int) taus.size(); i++){
00668 taus.at(i)->boostToRestFrame(mother);
00669
00670
00671 taus.at(i)->boostDaughtersToRestFrame(mother);
00672 }
00673
00674
00675 if(getTauPlus(taus)){
00676 *rotation_angle1 = getTauPlus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
00677 rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
00678 *rotation_angle2 = getTauPlus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
00679 rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
00680 }
00681 else{
00682 *rotation_angle1 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::Y_AXIS);
00683 rotateSystem(grandmothers,taus,*rotation_angle1,TauolaParticle::Y_AXIS);
00684 *rotation_angle2 = M_PI+getTauMinus(taus)->getRotationAngle(TauolaParticle::X_AXIS);
00685 rotateSystem(grandmothers,taus,*rotation_angle2,TauolaParticle::X_AXIS);
00686 }
00687
00688
00689 if(grandmothers.size()<1){
00690 *rotation_angle3 = 0;
00691 return;
00692 }
00693
00694
00695 if(getGrandmotherPlus(grandmothers))
00696 *rotation_angle3 = getGrandmotherPlus(grandmothers)->getRotationAngle(TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
00697 else
00698 *rotation_angle3 = M_PI+getGrandmotherMinus(grandmothers)->getRotationAngle(TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
00699
00700 rotateSystem(grandmothers,taus,*rotation_angle3,TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
00701
00702 }
00703
00704
00705
00706
00707 void TauolaParticlePair::boostFromTauPairToLabFrame(double rotation_angle1,
00708 double rotation_angle2,
00709 double rotation_angle3,
00710 TauolaParticle * mother,
00711 vector<TauolaParticle *> grandmothers,
00712 vector<TauolaParticle *> taus){
00713
00714
00715 if(mother==0)
00716 return;
00717
00718
00719
00720 rotateSystem(grandmothers,taus,-rotation_angle3, TauolaParticle::X_AXIS,TauolaParticle::Y_AXIS);
00721
00722
00723 rotateSystem(grandmothers,taus,-rotation_angle2,TauolaParticle::X_AXIS);
00724 rotateSystem(grandmothers,taus,-rotation_angle1,TauolaParticle::Y_AXIS);
00725
00726
00727
00728 for(int i=0; i< (int) grandmothers.size(); i++)
00729 grandmothers.at(i)->boostFromRestFrame(mother);
00730
00731 if(taus.size()!=1)
00732 for(int i=0; i< (int) taus.size(); i++){
00733 taus.at(i)->boostFromRestFrame(mother);
00734 taus.at(i)->boostDaughtersFromRestFrame(mother);
00735 }
00736 }
00737
00738 void TauolaParticlePair::rotateSystem(vector<TauolaParticle *> grandmothers,
00739 vector<TauolaParticle *> taus,
00740 double theta, int axis, int axis2){
00741 for(int i=0; i< (int) grandmothers.size(); i++)
00742 grandmothers.at(i)->rotate(axis,theta,axis2);
00743 for(int i=0; i< (int) taus.size(); i++){
00744 taus.at(i)->rotate(axis,theta,axis2);
00745 taus.at(i)->rotateDaughters(axis,theta,axis2);
00746 }
00747 }
00748
00749
00750
00751
00752
00753
00754 TauolaParticle * TauolaParticlePair::makeTemporaryMother(vector<TauolaParticle *> taus){
00755
00756
00757 double e=0;
00758 double px=0;
00759 double py=0;
00760 double pz=0;
00761
00762 bool tau_plus = false;
00763 bool tau_minus = false;
00764 bool tau_neutrino = false;
00765 bool tau_antineutrino = false;
00766
00767 for(int d=0; d < (int) taus.size(); d++){
00768 TauolaParticle * daughter = taus.at(d);
00769 e+=daughter->getE();
00770 px+=daughter->getPx();
00771 py+=daughter->getPy();
00772 pz+=daughter->getPz();
00773 switch(daughter->getPdgID()){
00774 case TauolaParticle::TAU_PLUS:
00775 tau_plus=true;
00776 break;
00777 case TauolaParticle::TAU_MINUS:
00778 tau_minus=true;
00779 break;
00780 case TauolaParticle::TAU_NEUTRINO:
00781 tau_neutrino=true;
00782 break;
00783 case TauolaParticle::TAU_ANTINEUTRINO:
00784 tau_antineutrino=true;
00785 break;
00786 }
00787 }
00788 double m = e*e-px*px-py*py-pz*pz;
00789 if(m<0)
00790 m= -sqrt(-m);
00791 else
00792 m=sqrt(m);
00793
00794
00795 int pdg=0;
00796
00797
00798 if(tau_plus&&tau_minus)
00799 pdg=TauolaParticle::Z0 ;
00800
00801
00802 if(tau_plus&&tau_neutrino)
00803 pdg=TauolaParticle::W_PLUS;
00804
00805
00806 if(tau_minus&&tau_antineutrino)
00807 pdg=TauolaParticle::W_MINUS;
00808
00809
00810 return taus.at(0)->createNewParticle(pdg,2,m,px,py,pz,e);
00811 }
00812
00813
00814
00815
00816 bool TauolaParticlePair::contains(TauolaParticle * particle){
00817
00818 for(int i=0; i< (int) m_final_particles.size(); i++){
00819 if(m_final_particles.at(i)->getBarcode()==particle->getBarcode())
00820 return true;
00821 }
00822 return false;
00823 }
00824
00825 TauolaParticle * TauolaParticlePair::getTauMinus(vector<TauolaParticle*> taus){
00826 for(int i=0; i< (int) taus.size(); i++){
00827 if(taus.at(i)->getPdgID()==TauolaParticle::TAU_MINUS)
00828 return taus.at(i);
00829 }
00830 return 0;
00831 }
00832
00833 TauolaParticle * TauolaParticlePair::getTauPlus(vector<TauolaParticle*> taus){
00834 for(int i=0; i< (int) taus.size(); i++){
00835 if(taus.at(i)->getPdgID()==TauolaParticle::TAU_PLUS)
00836 return taus.at(i);
00837 }
00838 return 0;
00839 }
00840
00841 TauolaParticle * TauolaParticlePair::getGrandmotherPlus(vector<TauolaParticle*> grandmothers){
00842
00843 double e = -1e30;
00844 int index = -1;
00845 for(int i=0; i< (int) grandmothers.size(); i++){
00846 if( (grandmothers.at(i)->getPdgID()<0 && grandmothers.at(i)->getPdgID()>-9) ||
00847 grandmothers.at(i)->getPdgID()== TauolaParticle::POSITRON||
00848 grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_PLUS){
00849 if(e<grandmothers.at(i)->getE()){
00850 e=grandmothers.at(i)->getE();
00851 index=i;
00852 }
00853 }
00854 }
00855 if(index==-1)
00856 {
00857 for(int i=0; i< (int) grandmothers.size(); i++)
00858 {
00859 if(grandmothers.at(i)->getPdgID()==21 || grandmothers.at(i)->getPdgID()==22)
00860 {
00861 index=i;
00862 e=grandmothers.at(i)->getE();
00863 break;
00864 }
00865 }
00866 }
00867 if(index==-1) index=0;
00868 if(e==0)
00869 {
00870 grandmothers.at(index)->print();
00871 return 0;
00872 }
00873 else
00874 return grandmothers.at(index);
00875 }
00876
00877 TauolaParticle * TauolaParticlePair::getGrandmotherMinus(vector<TauolaParticle*> grandmothers){
00878
00879 double e = -1e30;
00880 int index = -1;
00881 for(int i=0; i< (int) grandmothers.size(); i++){
00882 if( (grandmothers.at(i)->getPdgID()>0 && grandmothers.at(i)->getPdgID()<9) ||
00883 grandmothers.at(i)->getPdgID()== TauolaParticle::ELECTRON||
00884 grandmothers.at(i)->getPdgID()== TauolaParticle::MUON_MINUS){
00885 if(e<grandmothers.at(i)->getE()){
00886 e=grandmothers.at(i)->getE();
00887 index=i;
00888 }
00889 }
00890 }
00891 if(index==-1)
00892 {
00893 for(int i=(int) grandmothers.size()-1; i>=0 ; i--)
00894 {
00895 if(grandmothers.at(i)->getPdgID()==21||grandmothers.at(i)->getPdgID()==22)
00896 {
00897 index=i;
00898 e=grandmothers.at(i)->getE();
00899 break;
00900 }
00901 }
00902 }
00903 if(index==-1) index=0;
00904 if(e==0)
00905 return 0;
00906 else
00907 return grandmothers.at(index);
00908 }
00909
00910
00911
00912 void TauolaParticlePair::print(){
00913
00914 Log::RedirectOutput(Log::Info());
00915 std::cout << "Daughters final:" << std::endl;
00916 for(int i=0; i< (int) m_final_particles.size(); i++)
00917 m_final_particles.at(i)->print();
00918
00919
00920 std::cout << "Daughters at production:" << std::endl;
00921 for(int i=0; i< (int) m_production_particles.size(); i++)
00922 m_production_particles.at(i)->print();
00923
00924
00925 std::cout << "Mother particle: " << std::endl;
00926 if(m_mother)
00927 m_mother->print();
00928
00929 std::cout << "Grandmother particles: " << std::endl;
00930 for(int i=0; i< (int) m_grandmothers.size(); i++)
00931 m_grandmothers.at(i)->print();
00932
00933 Log::RevertOutput();
00934 }
00935
00936
00937 void TauolaParticlePair::checkMomentumConservation(){
00938
00939 for(int i=0; i< (int) m_final_particles.size(); i++)
00940 m_final_particles.at(i)->checkMomentumConservation();
00941
00942 for(int i=0; i< (int) m_production_particles.size(); i++)
00943 m_production_particles.at(i)->checkMomentumConservation();
00944
00945 if(m_mother)
00946 m_mother->checkMomentumConservation();
00947
00948 for(int i=0; i< (int) m_grandmothers.size(); i++)
00949 m_grandmothers.at(i)->checkMomentumConservation();
00950
00951 }
00952
00953 void TauolaParticlePair::recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta){
00954
00955 if (abs(outgoing_pdg_id)!=15)
00956 {
00957 Log::Warning()<<"interface was not used for taus pdg id="<<outgoing_pdg_id<<endl;
00958 return;
00959 }
00960
00961 double (*T) [Tauola::NCOS][4][4] = NULL;
00962 double (*Tw) [Tauola::NCOS] = NULL;
00963 double (*Tw0)[Tauola::NCOS] = NULL;
00964 double smin = 0.0;
00965 double smax = 0.0;
00966 double step = 0.0;
00967
00968
00969 switch(abs(incoming_pdg_id)){
00970 case 11:
00971 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
00972 {
00973 T = Tauola::table11B;
00974 Tw = Tauola::wtable11B;
00975 Tw0 = Tauola::w0table11B;
00976 smin = Tauola::sminB;
00977 smax = Tauola::smaxB;
00978 step = (smax-smin)/(Tauola::NS2-1);
00979 }
00980 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
00981 {
00982 T = Tauola::table11C;
00983 Tw = Tauola::wtable11C;
00984 Tw0 = Tauola::w0table11C;
00985 smin = Tauola::sminC;
00986 smax = Tauola::smaxC;
00987 step = (smax-smin)/(Tauola::NS3-1);
00988 }
00989 else
00990 {
00991 T = Tauola::table11A;
00992 Tw = Tauola::wtable11A;
00993 Tw0 = Tauola::w0table11A;
00994 smin = Tauola::sminA;
00995 smax = Tauola::smaxA;
00996 step = (smax-smin)/(Tauola::NS1-1);
00997 }
00998 break;
00999
01000
01001 case 1:
01002 case 3:
01003 case 5:
01004 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
01005 {
01006 T = Tauola::table1B;
01007 Tw = Tauola::wtable1B;
01008 Tw0 = Tauola::w0table1B;
01009 smin = Tauola::sminB;
01010 smax = Tauola::smaxB;
01011 step = (smax-smin)/(Tauola::NS2-1);
01012 }
01013 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
01014 {
01015 T = Tauola::table1C;
01016 Tw = Tauola::wtable1C;
01017 Tw0 = Tauola::w0table1C;
01018 smin = Tauola::sminC;
01019 smax = Tauola::smaxC;
01020 step = (smax-smin)/(Tauola::NS3-1);
01021 }
01022 else
01023 {
01024 T = Tauola::table1A;
01025 Tw = Tauola::wtable1A;
01026 Tw0 = Tauola::w0table1A;
01027 smin = Tauola::sminA;
01028 smax = Tauola::smaxA;
01029 step = (smax-smin)/(Tauola::NS1-1);
01030 }
01031 break;
01032
01033
01034 case 2:
01035 case 4:
01036 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
01037 {
01038 T = Tauola::table2B;
01039 Tw = Tauola::wtable2B;
01040 Tw0 = Tauola::w0table2B;
01041 smin = Tauola::sminB;
01042 smax = Tauola::smaxB;
01043 step = (smax-smin)/(Tauola::NS2-1);
01044 }
01045 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
01046 {
01047 T = Tauola::table2C;
01048 Tw = Tauola::wtable2C;
01049 Tw0 = Tauola::w0table2C;
01050 smin = Tauola::sminC;
01051 smax = Tauola::smaxC;
01052 step = (smax-smin)/(Tauola::NS3-1);
01053 }
01054 else
01055 {
01056 T = Tauola::table2A;
01057 Tw = Tauola::wtable2A;
01058 Tw0 = Tauola::w0table2A;
01059 smin = Tauola::sminA;
01060 smax = Tauola::smaxA;
01061 step = (smax-smin)/(Tauola::NS1-1);
01062 }
01063 break;
01064
01065 default:
01066 Log::Warning()<<"interface was not used for proper beams pdg id="<<incoming_pdg_id<<endl;
01067 return;
01068 }
01069
01070
01071 if (T[0][0][0][0]<0.5) return;
01072
01073
01074 if (invariant_mass_squared <= exp(Tauola::sminA)){
01075 double mta = 1.777;
01076 double cosf = cosTheta;
01077 double s = invariant_mass_squared;
01078 double sinf = sqrt(1-cosf*cosf);
01079 double MM = sqrt(4e0*mta*mta/s);
01080 double xnorm = 1+cosf*cosf + MM*MM*sinf*sinf;
01081
01082
01083 for (int i=0;i<4;i++)
01084 for (int j=0;j<4;j++)
01085 m_R[i][j]=0.0;
01086
01087 m_R[0][0] = (1+cosf*cosf + MM*MM*sinf*sinf)/xnorm;
01088 m_R[1][1] = (-(1- MM*MM)*sinf*sinf)/xnorm;
01089 m_R[2][2] = ( (1+ MM*MM)*sinf*sinf)/xnorm;
01090 m_R[3][3] = (1+cosf*cosf - MM*MM*sinf*sinf)/xnorm;
01091 m_R[2][3] = (2*MM*sinf*cosf)/xnorm;
01092 m_R[3][2] = (2*MM*sinf*cosf)/xnorm;
01093
01094
01095 double w = 1.;
01096 double w0 = 1.;
01097
01098 if (Tauola::wtable2A[0][0]>0 ) w=Tauola::wtable2A[0][0];
01099 else if(Tauola::wtable1A[0][0]>0 ) w=Tauola::wtable1A[0][0];
01100 else if(Tauola::wtable11A[0][0]>0) w=Tauola::wtable11A[0][0];
01101
01102 if (Tauola::wtable2A[0][0]>0 ) w0=Tauola::w0table2A[0][0];
01103 else if(Tauola::wtable1A[0][0]>0 ) w0=Tauola::w0table1A[0][0];
01104 else if(Tauola::wtable11A[0][0]>0) w0=Tauola::w0table11A[0][0];
01105
01106
01107 Tauola::setEWwt(w,w0);
01108 return;
01109 }
01110
01111 int x = 0;
01112 double buf = smin;
01113 double remnant = 0.0;
01114
01115
01116 if(T==Tauola::table11A || T==Tauola::table1A || T== Tauola::table2A)
01117 {
01118 while(log(invariant_mass_squared) > buf){
01119 x++;
01120 buf += (step);
01121 }
01122 remnant = (log(invariant_mass_squared)-(buf-step))/step;
01123 }
01124 else
01125 {
01126 while(invariant_mass_squared > buf){
01127 x++;
01128 buf += step;
01129 }
01130 remnant = (invariant_mass_squared-(buf-step))/step;
01131 }
01132
01133 double cmin = -1.;
01134 int y = 0;
01135 double remnantc = 0.0;
01136
01137
01138 buf = cmin+1./Tauola::NCOS;
01139 while(cosTheta > buf){
01140 y++;
01141 buf+=2./Tauola::NCOS;
01142 }
01143
01144 remnantc = (cosTheta-(buf-2./Tauola::NCOS))/(2./Tauola::NCOS);
01145
01146
01147 bool isUsingExtrapolation = false;
01148 if(y >= Tauola::NCOS) { isUsingExtrapolation = true; y = Tauola::NCOS-1; }
01149 if(y == 0) { isUsingExtrapolation = true; y = 1; }
01150
01151
01152 double b11,b21,b12,b22;
01153
01154 for (int i=0; i<4; i++){
01155 for (int j=0; j<4; j++){
01156
01157 if(isUsingExtrapolation){
01158 if(y == 1){
01159 b11 = 2*T[x-1][0][i][j] - T[x-1][1][i][j];
01160 b21 = 2*T[x][0][i][j] - T[x][1][i][j];
01161 b12 = T[x-1][0][i][j];
01162 b22 = T[x][0][i][j];
01163 }
01164 else{
01165 b11 = T[x-1][y][i][j];
01166 b21 = T[x][y][i][j];
01167 b12 = 2*T[x-1][y][i][j] - T[x-1][y-1][i][j];
01168 b22 = 2*T[x][y][i][j] - T[x][y-1][i][j];
01169 }
01170 }
01171 else{
01172 b11 = T[x-1][y-1][i][j];
01173 b21 = T[x][y-1][i][j];
01174 b12 = T[x-1][y][i][j];
01175 b22 = T[x][y][i][j];
01176 }
01177 m_R[i][j] = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
01178 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
01179 }
01180 }
01181
01182
01183 double w,w0;
01184
01185 if(isUsingExtrapolation){
01186 if(y == 1){
01187 b11 = 2*Tw[x-1][0] - Tw[x-1][1];
01188 b21 = 2*Tw[x][0] - Tw[x][1];
01189 b12 = Tw[x-1][0];
01190 b22 = Tw[x][0];
01191 }
01192 else
01193 {
01194 b11 = Tw[x-1][y];
01195 b21 = Tw[x][y];
01196 b12 = 2*Tw[x-1][y] - Tw[x-1][y-1];
01197 b22 = 2*Tw[x][y] - Tw[x][y-1];
01198 }
01199 }
01200 else
01201 {
01202 b11 = Tw[x-1][y-1];
01203 b21 = Tw[x][y-1];
01204 b12 = Tw[x-1][y];
01205 b22 = Tw[x][y];
01206 }
01207
01208 w = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
01209 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
01210
01211 if(isUsingExtrapolation){
01212 if(y == 1){
01213 b11 = 2*Tw0[x-1][0] - Tw0[x-1][1];
01214 b21 = 2*Tw0[x][0] - Tw0[x][1];
01215 b12 = Tw0[x-1][0];
01216 b22 = Tw0[x][0];
01217 }
01218 else{
01219 b11 = Tw0[x-1][y];
01220 b21 = Tw0[x][y];
01221 b12 = 2*Tw0[x-1][y] - Tw0[x-1][y-1];
01222 b22 = 2*Tw0[x][y] - Tw0[x][y-1];
01223 }
01224 }
01225 else{
01226 b11 = Tw0[x-1][y-1];
01227 b21 = Tw0[x][y-1];
01228 b12 = Tw0[x-1][y];
01229 b22 = Tw0[x][y];
01230 }
01231
01232 w0 = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
01233 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
01234
01235 Tauola::setEWwt(w,w0);
01236 }
01237
01238 }