114 int incoming_pdg_id=0;
115 int outgoing_pdg_id=0;
116 double invariant_mass_squared=-5.0;
120 for(
int x = 0; x < 4; x ++) {
121 for(
int y = 0; y < 4; y ++)
149 if(!Tauola::spin_correlation.HIGGS_H)
return;
151 double phi = Tauola::getHiggsScalarPseudoscalarMixingAngle();
154 double beta = sqrt(1-4*mass_ratio*mass_ratio);
155 double denominator = pow(cos(phi)*beta,2)+pow(sin(phi),2);
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;
171 if(!Tauola::spin_correlation.Z0)
break;
177 pz =
getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
183 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
187 if(!Tauola::spin_correlation.GAMMA)
break;
193 pz =
getZPolarization(&incoming_pdg_id, &outgoing_pdg_id, &invariant_mass_squared, &cosTheta);
199 recalculateRij(incoming_pdg_id, outgoing_pdg_id, invariant_mass_squared, cosTheta);
204 if(!Tauola::spin_correlation.HIGGS)
break;
213 if(!Tauola::spin_correlation.HIGGS_A)
break;
221 if(!Tauola::spin_correlation.HIGGS_PLUS)
break;
228 if(!Tauola::spin_correlation.HIGGS_MINUS)
break;
235 if(!Tauola::spin_correlation.W_PLUS)
break;
242 if(!Tauola::spin_correlation.W_MINUS)
break;
274 *invariant_mass_squared = pow(
m_mother->getMass(),2);
275 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
281 Log::Warning()<<
"Not enough mothers of Z to "
282 <<
"calculate cos(theta) between "
283 <<
"incoming about outgoing beam"
285 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
286 invariant_mass_squared, cosTheta);
290 Log::Error()<<
"tau+ or tau- not found in Z decay"<< endl;
303 vector<TauolaParticle*> extra_grandmothers;
313 vector<TauolaParticle*> extra_tau_siblings;
314 vector<TauolaParticle*> temp_production_particles;
322 vector<TauolaParticle*> extra_Z_siblings;
325 extra_Z_siblings.push_back(
m_grandmothers.at(0)->getDaughters().at(i));
330 std::vector<TauolaParticle*> effective_taus;
342 for(
int i=0; i<(int) extra_grandmothers.size(); i++)
346 for(
int i=0; i<(int) extra_Z_siblings.size(); i++)
350 for(
int i=0; i<(int) extra_tau_siblings.size() ; i++){
352 addToBeam(extra_tau_siblings.at(i),&effective_taus,0);
359 *invariant_mass_squared = pow(temp_mother->
getMass(),2);
364 double angle1,angle2,angle3;
384 double sintheta1 = sqrt(1-costheta1*costheta1);
385 double sintheta2 = sqrt(1-costheta2*costheta2);
386 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
392 if(abs(*incoming_pdg_id)==21 || abs(*incoming_pdg_id)==22)
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 )
403 Log::Error()<<
"Second class disaster: incoming_pdg_id = "<<*incoming_pdg_id<<endl;
410 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id,
411 invariant_mass_squared, cosTheta);
416 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
417 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
424 double angle1,angle2,angle3;
441 double sintheta1 = sqrt(1-costheta1*costheta1);
442 double sintheta2 = sqrt(1-costheta2*costheta2);
444 double avg = (costheta1*sintheta2+costheta2*sintheta1)/(sintheta1+sintheta2);
453 setBornKinematics(*incoming_pdg_id, *outgoing_pdg_id, *invariant_mass_squared, *cosTheta);
454 return 1-plzap0_(incoming_pdg_id,outgoing_pdg_id, invariant_mass_squared, cosTheta);
464 std::vector<TauolaParticle*> *candidates_same,
465 std::vector<TauolaParticle*> *candidates_opp){
473 double s0 =1.0/
getVirtuality(pcle,candidates_same->at(0),
false);
474 double s1 = 1.0/
getVirtuality(pcle,candidates_same->at(1),
false);
477 s_beam = candidates_same->at(0);
481 s_beam = candidates_same->at(1);
486 double o0 =1.0/
getVirtuality(pcle,candidates_opp->at(0),
true);
487 double o1 = 1.0/
getVirtuality(pcle,candidates_opp->at(1),
true);
490 o_beam = candidates_opp->at(0);
494 o_beam = candidates_opp->at(1);
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);
511 if((abs(pdg2)==21 || abs(pdg2)==22) && abs(pdg1)<17 && abs(pdg1)!=10 && abs(pdg1)!=9) o_beam->
setPdgID(-pdg1);
565 double h_tau_minus[4]={2,0,0,0};
566 double h_tau_plus[4]={2,0,0,0};
572 for(
double weight=0; weight < Tauola::randomDouble();){
575 Tauola::redefineTauMinusProperties(tau_minus);
583 Tauola::redefineTauPlusProperties(tau_plus);
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];
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]);
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);
618 double angle1,angle2,angle3;
653 double * rotation_angle2,
654 double * rotation_angle3,
656 vector<TauolaParticle *> grandmothers,
657 vector<TauolaParticle *> taus
663 for(
int i=0; i< (int) grandmothers.size(); i++)
664 grandmothers.at(i)->boostToRestFrame(mother);
667 for(
int i=0; i< (int) taus.size(); i++){
668 taus.at(i)->boostToRestFrame(mother);
671 taus.at(i)->boostDaughtersToRestFrame(mother);
689 if(grandmothers.size()<1){
690 *rotation_angle3 = 0;
708 double rotation_angle2,
709 double rotation_angle3,
711 vector<TauolaParticle *> grandmothers,
712 vector<TauolaParticle *> taus){
728 for(
int i=0; i< (int) grandmothers.size(); i++)
729 grandmothers.at(i)->boostFromRestFrame(mother);
732 for(
int i=0; i< (int) taus.size(); i++){
733 taus.at(i)->boostFromRestFrame(mother);
734 taus.at(i)->boostDaughtersFromRestFrame(mother);
955 if (abs(outgoing_pdg_id)!=15)
957 Log::Warning()<<
"interface was not used for taus pdg id="<<outgoing_pdg_id<<endl;
961 double (*T) [Tauola::NCOS][4][4] = NULL;
962 double (*Tw) [Tauola::NCOS] = NULL;
963 double (*Tw0)[Tauola::NCOS] = NULL;
969 switch(abs(incoming_pdg_id)){
971 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
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);
980 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
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);
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);
1004 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
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);
1013 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
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);
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);
1036 if(invariant_mass_squared<Tauola::smaxB && invariant_mass_squared>Tauola::sminB)
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);
1045 else if (invariant_mass_squared<Tauola::smaxC && invariant_mass_squared>Tauola::sminC)
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);
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);
1066 Log::Warning()<<
"interface was not used for proper beams pdg id="<<incoming_pdg_id<<endl;
1071 if (T[0][0][0][0]<0.5)
return;
1074 if (invariant_mass_squared <= exp(Tauola::sminA)){
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;
1083 for (
int i=0;i<4;i++)
1084 for (
int j=0;j<4;j++)
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;
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];
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];
1107 Tauola::setEWwt(w,w0);
1113 double remnant = 0.0;
1116 if(T==Tauola::table11A || T==Tauola::table1A || T== Tauola::table2A)
1118 while(log(invariant_mass_squared) > buf){
1122 remnant = (log(invariant_mass_squared)-(buf-step))/step;
1126 while(invariant_mass_squared > buf){
1130 remnant = (invariant_mass_squared-(buf-step))/step;
1135 double remnantc = 0.0;
1138 buf = cmin+1./Tauola::NCOS;
1139 while(cosTheta > buf){
1141 buf+=2./Tauola::NCOS;
1144 remnantc = (cosTheta-(buf-2./Tauola::NCOS))/(2./Tauola::NCOS);
1147 bool isUsingExtrapolation =
false;
1148 if(y >= Tauola::NCOS) { isUsingExtrapolation =
true; y = Tauola::NCOS-1; }
1149 if(y == 0) { isUsingExtrapolation =
true; y = 1; }
1152 double b11,b21,b12,b22;
1154 for (
int i=0; i<4; i++){
1155 for (
int j=0; j<4; j++){
1157 if(isUsingExtrapolation){
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];
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];
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];
1177 m_R[i][j] = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1178 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1185 if(isUsingExtrapolation){
1187 b11 = 2*Tw[x-1][0] - Tw[x-1][1];
1188 b21 = 2*Tw[x][0] - Tw[x][1];
1196 b12 = 2*Tw[x-1][y] - Tw[x-1][y-1];
1197 b22 = 2*Tw[x][y] - Tw[x][y-1];
1208 w = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1209 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1211 if(isUsingExtrapolation){
1213 b11 = 2*Tw0[x-1][0] - Tw0[x-1][1];
1214 b21 = 2*Tw0[x][0] - Tw0[x][1];
1221 b12 = 2*Tw0[x-1][y] - Tw0[x-1][y-1];
1222 b22 = 2*Tw0[x][y] - Tw0[x][y-1];
1226 b11 = Tw0[x-1][y-1];
1232 w0 = b11*(1-remnant)*(1-remnantc) + b21*remnant*(1-remnantc)
1233 + b12*(1-remnant)*remnantc + b22*remnant*remnantc;
1235 Tauola::setEWwt(w,w0);
TauolaParticle * getGrandmotherMinus(std::vector< TauolaParticle * > particles)
TauolaParticle * getGrandmotherPlus(std::vector< TauolaParticle * > particles)
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 * getTauMinus(std::vector< TauolaParticle * > particles)
TauolaParticle * getTauPlus(std::vector< TauolaParticle * > particles)
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)
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)
double getRotationAngle(int axis, int second_axis=Z_AXIS)