1 #include "TauSpinner/tau_reweight_lib.h" 2 #include "TauSpinner/nonSM.h" 4 #include "Tauola/TauolaParticlePair.h" 11 double CMSENE = 7000.0;
52 double f(
double x,
int ID,
double SS,
double cmsene)
69 return LHAPDF::xfx(x, sqrt(SS), ID)/x;
81 double sigborn(
int ID,
double SS,
double costhe)
95 SIGggHiggs=
disth_(&SS, &costhe, &IPOne, &IPOne)+
disth_(&SS, &costhe, &IPOne, &IMOne)+
96 disth_(&SS, &costhe, &IMOne, &IPOne)+
disth_(&SS, &costhe, &IMOne, &IMOne);
99 double PI=3.14159265358979324;
100 SIGggHiggs *= XMH * XMH * XMH *XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) + XGH*XGH*XMH*XMH);
103 if(IfHsimple==1) SIGggHiggs = Xnorm * XMH * XGH /PI/ ((SS-XMH*XMH)*(SS-XMH*XMH) +XGH*XGH*XMH*XMH);
111 if (ID==0)
return 0.0 ;
112 if (ID>0) initwk_( &ID, &tauID, &SS);
116 initwk_( &ID, &tauID, &SS);
123 return ( t_born_(&iZero, &SS, &costhe, &dOne , &dOne) + t_born_(&iZero, &SS, &costhe, &dOne , &dMOne)
124 + t_born_(&iZero, &SS, &costhe, &dMOne, &dOne) + t_born_(&iZero, &SS, &costhe, &dMOne, &dMOne))/SS/128.86674175/128.86674175/3.141559265358979324;
142 cout<<
" ------------------------------------------------------"<<endl;
143 cout<<
" TauSpinner v2.0.4"<<endl;
144 cout<<
" -----------------"<<endl;
145 cout<<
" 15.Aug.2017 "<<endl;
146 cout<<
" by Z. Czyczula (until 2015), T. Przedzinski, E. Richter-Was, Z. Was,"<<endl;
147 cout<<
" matrix elements implementations "<<endl;
148 cout<<
" also J. Kalinowski, W. Kotlarski and M. Bachmani"<<endl;
149 cout<<
" ------------------------------------------------------"<<endl;
150 cout<<
" Ipp - true for pp collision; otherwise polarization"<<endl;
151 cout<<
" of individual taus from Z/gamma* is set to 0.0"<<endl;
152 cout<<
" Ipp = "<<Ipp<<endl;
153 cout<<
" CMSENE - used in PDF calculations; only if Ipp = true"<<endl;
154 cout<<
" and only for Z/gamma*"<<endl;
155 cout<<
" CMSENE = "<<CMSENE<<endl;
156 cout<<
" Ipol - relevant for Z/gamma* decays "<<endl;
157 cout<<
" 0 - events generated without spin effects "<<endl;
158 cout<<
" 1 - events generated with all spin effects "<<endl;
159 cout<<
" 2 - events generated with spin correlations and <pol>=0 "<<endl;
160 cout<<
" 3 - events generated with spin correlations and"<<endl;
161 cout<<
" polarization but missing angular dependence of <pol>"<<endl;
162 cout<<
" Ipol = "<<Ipol<<endl;
163 cout<<
" Ipol - relevant for Z/gamma* decays "<<endl;
164 cout<<
" NOTE: For Ipol=0,1 algorithm is identical. "<<endl;
165 cout<<
" However in user program role of wt need change. "<<endl;
166 cout<<
" nonSM2 = "<<nonSM2<<endl;
167 cout<<
" 1/0 extra term in cross section, density matrix on/off "<<endl;
168 cout<<
" nonSMN = "<<nonSMN<<endl;
169 cout<<
" 1/0 extra term in cross section, for shapes only? on/off "<<endl;
170 cout<<
" note KEY - for options of matrix elements calculations "<<endl;
171 cout<<
" in cases of final states with two jets "<<endl;
172 cout<<
" ------------------------------------------------------ "<<endl;
183 relWTnonSM = _relWTnonSM;
195 Xnorm = normalization;
240 *normalization = Xnorm;
252 void setSpinOfSample(
int _Ipol)
306 Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
307 Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
308 Particle nu_tau(sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
310 vector<Particle> tau_daughters;
313 int tau_pdgid = sp_tau.pdgid();
316 for(
unsigned int i=0; i<sp_tau_daughters.size(); i++)
318 Particle pp(sp_tau_daughters[i].px(),
319 sp_tau_daughters[i].py(),
320 sp_tau_daughters[i].pz(),
321 sp_tau_daughters[i].e(),
322 sp_tau_daughters[i].pdgid() );
324 tau_daughters.push_back(pp);
327 double phi2 = 0.0, theta2 = 0.0;
337 double *HH =
calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
340 if ( abs(sp_X.pdgid()) == 24 ) { sign= 1.0; }
341 else if( abs(sp_X.pdgid()) == 37 ) { sign=-1.0; }
344 cout<<
"wrong sp_W/H.pdgid()="<<sp_X.pdgid()<<endl;
347 if (sp_X.pdgid() > 0 )
348 {WTamplitM = WTamplit;}
350 {WTamplitP = WTamplit;}
353 double WT = 1.0+sign*HH[2];
358 cout<<tau_pdgid<<
" -> ";
359 for(
unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<
" ";
360 cout<<
" (HH: "<<HH[0]<<
" "<<HH[1]<<
" "<<HH[2]<<
" "<<HH[3]<<
") WT: "<<WT<<endl;
365 printf(
"TauSpinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 0.0\n",WT);
370 printf(
"Tauspinner::calculateWeightFromParticlesWorHpn WT is: %13.10f. Setting WT = 2.0\n",WT);
402 vector<SimpleParticle> sp_tau_daughters;
407 if (sp_tau1.pdgid() == -15 )
411 sp_tau_daughters = sp_tau1_daughters;
417 sp_tau_daughters = sp_tau2_daughters;
426 Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
427 Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
428 Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
430 vector<Particle> tau_daughters;
433 int tau_pdgid = sp_tau.pdgid();
436 for(
unsigned int i=0; i<sp_tau_daughters.size(); i++)
438 Particle pp(sp_tau_daughters[i].px(),
439 sp_tau_daughters[i].py(),
440 sp_tau_daughters[i].pz(),
441 sp_tau_daughters[i].e(),
442 sp_tau_daughters[i].pdgid() );
444 tau_daughters.push_back(pp);
447 double phi2 = 0.0, theta2 = 0.0;
459 HHp =
calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
463 cout<<tau_pdgid<<
" -> ";
464 for(
unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<
" ";
465 cout<<
" (HHp: "<<HHp[0]<<
" "<<HHp[1]<<
" "<<HHp[2]<<
" "<<HHp[3]<<
") ";
469 WTamplitP = WTamplit;
476 if(sp_tau1.pdgid() == 15 )
480 sp_tau_daughters = sp_tau1_daughters;
486 sp_tau_daughters = sp_tau2_daughters;
493 Particle X ( sp_X.px(), sp_X.py(), sp_X.pz(), sp_X.e(), sp_X.pdgid() );
494 Particle tau ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
495 Particle nu_tau( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
497 vector<Particle> tau_daughters;
500 int tau_pdgid = sp_tau.pdgid();
503 for(
unsigned int i=0; i<sp_tau_daughters.size(); i++)
505 Particle pp(sp_tau_daughters[i].px(),
506 sp_tau_daughters[i].py(),
507 sp_tau_daughters[i].pz(),
508 sp_tau_daughters[i].e(),
509 sp_tau_daughters[i].pdgid() );
511 tau_daughters.push_back(pp);
514 double phi2 = 0.0, theta2 = 0.0;
525 HHm =
calculateHH(tau_pdgid, tau_daughters, phi2, theta2);
529 cout<<tau_pdgid<<
" -> ";
530 for(
unsigned int i=0;i<tau_daughters.size();i++) cout<<tau_daughters[i].pdgid()<<
" ";
531 cout<<
" (HHm: "<<HHm[0]<<
" "<<HHm[1]<<
" "<<HHm[2]<<
" "<<HHm[3]<<
") ";
535 WTamplitM = WTamplit;
544 if(sp_X.pdgid() == 25) { sign=-1.0; }
545 if(sp_X.pdgid() == 36) { sign=-1.0; }
546 if(sp_X.pdgid() ==553) { sign=-1.0; }
554 double S = sp_X.e()*sp_X.e() - sp_X.px()*sp_X.px() - sp_X.py()*sp_X.py() - sp_X.pz()*sp_X.pz();
570 WT = 1.0+corrX2*HHp[2]*HHm[2]+polX2*HHp[2]+polX2*HHm[2] + RzXX*HHp[0]*HHm[0] + RzYY*HHp[1]*HHm[1] + RzXY*HHp[0]*HHm[1] + RzYX*HHp[1]*HHm[0];
573 double RRR = Tauola::randomDouble();
575 if (RRR<(1.0+polX2)*(1.0+corrX2*HHp[2]*HHm[2]+HHp[2]+HHm[2])/(2.0+2.0*corrX2*HHp[2]*HHm[2]+2.0*polX2*HHp[2]+2.0*polX2*HHm[2])) Polari=-1.0;
579 WT = 1.0 + sign*HHp[2]*HHm[2] + RXX*HHp[0]*HHm[0] + RYY*HHp[1]*HHm[1] + RXY*HHp[0]*HHm[1] + RYX*HHp[1]*HHm[0];
582 double RRR = Tauola::randomDouble();
584 if (RRR<(1.0+sign*HHp[2]*HHm[2]+HHp[2]-HHm[2])/(2.0+2.0*sign*HHp[2]*HHm[2])) Polari=-1.0;
590 double S = sp_X.e()*sp_X.e() - sp_X.px()*sp_X.px() - sp_X.py()*sp_X.py() - sp_X.pz()*sp_X.pz();
599 WT = 1.0+sign*HHp[2]*HHm[2]+pol*HHp[2]+pol*HHm[2] + RzXX*R11*HHp[0]*HHm[0] - RzYY*R22*HHp[1]*HHm[1] + RzXY*R12*HHp[0]*HHm[1] + RzYX*R21*HHp[1]*HHm[0];
605 if(Ipol==2) WT = WT/(1.0+sign*HHp[2]*HHm[2]);
612 double polp1 =
plzap2(11,sp_tau.pdgid(),S,0.0);
613 double pol1 =(2*(1-polp1)-1) ;
614 WT = WT/(1.0+sign*HHp[2]*HHm[2]+pol1*HHp[2]+pol1*HHm[2]);
618 double RRR = Tauola::randomDouble();
620 if (RRR<(1.0+pol)*(1.0+sign*HHp[2]*HHm[2]+HHp[2]+HHm[2])/(2.0+2.0*sign*HHp[2]*HHm[2]+2.0*pol*HHp[2]+2.0*pol*HHm[2])) Polari=-1.0;
624 DEBUG( cout<<
" WT: "<<WT<<endl; )
627 printf(
"Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 0.0\n",WT);
632 printf(
"Tauspinner::calculateWeightFromParticlesH WT is: %13.10f. Setting WT = 4.0\n",WT);
636 if( WT>4.0 || WT<0.0)
638 cout<<
"Tauspinner::calculateWeightFromParticlesH ERROR: Z/gamma* or H, and WT not in range [0,4]."<<endl;
642 if (sign==-1.0 && nonSM2!=1) {
645 cout <<
"Tauspinner::calculateWeightFromParticlesH Setting WT to be 2.0" << endl;
649 if( sign==-1.0 && (WT>2.0 || WT<0.0) && nonSM2!=1)
651 cout<<
"Tauspinner::calculateWeightFromParticlesH ERROR: H and WT not in range [0,2]."<<endl;
669 Particle P_QQ( tau.px()+nu_tau.px(), tau.py()+nu_tau.py(), tau.pz()+nu_tau.pz(), tau.e()+nu_tau.e(), 0 );
676 tau.boostToRestFrame(P_QQ);
677 nu_tau.boostToRestFrame(P_QQ);
679 for(
unsigned int i=0; i<tau_daughters.size();i++)
680 tau_daughters[i].boostToRestFrame(P_QQ);
688 double phi = tau.getAnglePhi();
692 double theta = tau.getAngleTheta();
694 tau.rotateXZ(M_PI-theta);
696 nu_tau.rotateXY(-phi );
697 nu_tau.rotateXZ(M_PI-theta);
699 for(
unsigned int i=0; i<tau_daughters.size();i++)
701 tau_daughters[i].rotateXY(-phi );
702 tau_daughters[i].rotateXZ(M_PI-theta);
710 for(
unsigned int i=0; i<tau_daughters.size();i++)
711 tau_daughters[i].boostAlongZ(-tau.pz(),tau.e());
721 unsigned int i_stored = 0;
726 for(
unsigned int i=0; i<tau_daughters.size();i++)
728 if(abs(tau_daughters[i].pdgid())==16){
729 *phi2 = tau_daughters[i].getAnglePhi();
731 tau_daughters[i].rotateXY( -(*phi2) );
733 *theta2 = tau_daughters[i].getAngleTheta();
735 tau_daughters[i].rotateXZ( -(*theta2) );
742 for(
unsigned int i=0; i<tau_daughters.size();i++)
745 tau_daughters[i].rotateXY( -(*phi2) );
746 tau_daughters[i].rotateXZ( -(*theta2) );
764 double*
calculateHH(
int tau_pdgid, vector<Particle> &tau_daughters,
double phi,
double theta)
767 double *HH =
new double[4];
769 HH[0]=HH[1]=HH[2]=HH[3]=0.0;
774 for(
unsigned int i=0; i<tau_daughters.size(); i++)
775 pdgid.push_back( tau_daughters[i].pdgid() );
789 if( pdgid.size()==2 &&
791 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211) ) ||
792 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211) ) ||
793 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321) ) ||
794 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321) )
798 if(abs(pdgid[1])==321) channel = 6;
799 DEBUG( cout<<
"Channel "<<channel<<
" : "; )
807 const double AMTAU = 1.777;
809 double AMPI = sqrt(tau_daughters[1].e() *tau_daughters[1].e()
810 -tau_daughters[1].px()*tau_daughters[1].px()
811 -tau_daughters[1].py()*tau_daughters[1].py()
812 -tau_daughters[1].pz()*tau_daughters[1].pz());
815 double PXQ=AMTAU*tau_daughters[1].e();
816 double PXN=AMTAU*tau_daughters[0].e();
817 double QXN=tau_daughters[1].e()*tau_daughters[0].e()-tau_daughters[1].px()*tau_daughters[0].px()-tau_daughters[1].py()*tau_daughters[0].py()-tau_daughters[1].pz()*tau_daughters[0].pz();
818 double BRAK=(2*PXQ*QXN-AMPI*AMPI*PXN);
820 WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;
821 HH[0] = AMTAU*(2*tau_daughters[1].px()*QXN-tau_daughters[0].px()*AMPI*AMPI)/BRAK;
822 HH[1] = AMTAU*(2*tau_daughters[1].py()*QXN-tau_daughters[0].py()*AMPI*AMPI)/BRAK;
823 HH[2] = AMTAU*(2*tau_daughters[1].pz()*QXN-tau_daughters[0].pz()*AMPI*AMPI)/BRAK;
831 else if( pdgid.size()==3 &&
833 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 111) ) ||
834 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 111) ) ||
835 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321, 311) ) ||
836 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 311) ) ||
837 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321, 310) ) ||
838 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 310) ) ||
839 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321, 130) ) ||
840 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 130) )
846 DEBUG( cout<<
"Channel "<<channel<<
" : "; )
853 const double AMTAU = 1.777;
856 if(tau_daughters[2].pdgid() != 111) { MNUM=3; channel = 22;}
857 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
858 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
859 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
860 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
862 float HV[4] = { 0.0 };
864 dam2pi_( &MNUM, PT, PN, PIM1, PIM2, &LIT, HV );
879 else if( pdgid.size()==3 &&
881 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 130) ) ||
882 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 130) ) ||
883 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 310) ) ||
884 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 310) ) ||
885 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 311) ) ||
886 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 311) ) ||
887 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321, 111) ) ||
888 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 111) )
893 DEBUG( cout<<
"Channel "<<channel<<
" : "; )
900 const double AMTAU = 1.777;
903 QQ[0]=tau_daughters[1].e() -tau_daughters[2].e() ;
904 QQ[1]=tau_daughters[1].px()-tau_daughters[2].px();
905 QQ[2]=tau_daughters[1].py()-tau_daughters[2].py();
906 QQ[3]=tau_daughters[1].pz()-tau_daughters[2].pz();
909 PKS[0]=tau_daughters[1].e() +tau_daughters[2].e() ;
910 PKS[1]=tau_daughters[1].px()+tau_daughters[2].px();
911 PKS[2]=tau_daughters[1].py()+tau_daughters[2].py();
912 PKS[3]=tau_daughters[1].pz()+tau_daughters[2].pz();
915 double PKSD=PKS[0]*PKS[0]-PKS[1]*PKS[1]-PKS[2]*PKS[2]-PKS[3]*PKS[3];
916 double QQPKS=QQ[0]*PKS[0]-QQ[1]*PKS[1]-QQ[2]*PKS[2]-QQ[3]*PKS[3];
918 QQ[0]=QQ[0]-PKS[0]*QQPKS/PKSD;
919 QQ[1]=QQ[1]-PKS[1]*QQPKS/PKSD;
920 QQ[2]=QQ[2]-PKS[2]*QQPKS/PKSD;
921 QQ[3]=QQ[3]-PKS[3]*QQPKS/PKSD;
923 double PRODPQ=AMTAU*QQ[0];
924 double PRODNQ=tau_daughters[0].e() *QQ[0]
925 -tau_daughters[0].px()*QQ[1]
926 -tau_daughters[0].py()*QQ[2]
927 -tau_daughters[0].pz()*QQ[3];
928 double PRODPN=AMTAU*tau_daughters[0].e();
929 double QQ2 =QQ[0]*QQ[0]-QQ[1]*QQ[1]-QQ[2]*QQ[2]-QQ[3]*QQ[3];
932 double BRAK=(2*PRODPQ*PRODNQ-PRODPN*QQ2);
934 WTamplit = (1.16637E-5)*(1.16637E-5)*BRAK/2.;
935 HH[0]=AMTAU*(2*PRODNQ*QQ[1]-QQ2*tau_daughters[0].px())/BRAK;
936 HH[1]=AMTAU*(2*PRODNQ*QQ[2]-QQ2*tau_daughters[0].py())/BRAK;
937 HH[2]=AMTAU*(2*PRODNQ*QQ[3]-QQ2*tau_daughters[0].pz())/BRAK;
943 else if( pdgid.size()==3 &&
945 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 11,-12) ) ||
946 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16,-11, 12) )
949 DEBUG( cout<<
"Channel 1 : "; )
955 double XK0DEC = 0.01;
956 double XK[4] = { 0.0 };
957 double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
958 double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
959 double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
961 double HV[4] = { 0.0 };
965 QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
966 XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
968 dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &LIT, HV );
979 else if( pdgid.size()==4 &&
981 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 11,-12, 22) ) ||
982 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16,-11, 12, 22) )
985 DEBUG( cout<<
"Channel 1b : "; )
991 double XK0DEC = 0.01;
992 double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
993 double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
994 double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
995 double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
997 double HV[4] = { 0.0 };
1001 QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.511e-3*0.511e-3);
1002 XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1003 XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1005 if(XK0DEC > XK[3]/(XK[3]+XA[3]+QP[3]+XN[3])) XK0DEC=0.5*XK[3]/(XK[3]+XA[3]+QP[3]+XN[3]);
1007 dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &LIT, HV );
1018 else if( pdgid.size()==3 &&
1020 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 13,-14) ) ||
1021 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16,-13, 14) )
1025 DEBUG( cout<<
"Channel 2 : "; )
1031 double XK0DEC = 0.01;
1032 double XK[4] = { 0.0 };
1033 double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1034 double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1035 double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1036 double AMPLIT = 0.0;
1037 double HV[4] = { 0.0 };
1041 QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1042 XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1044 dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &LIT, HV );
1055 else if( pdgid.size()==4 &&
1057 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 13,-14, 22) ) ||
1058 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16,-13, 14, 22) )
1062 DEBUG( cout<<
"Channel 2b : "; )
1068 double XK0DEC = 0.01;
1069 double XK[4] = { tau_daughters[3].px(), tau_daughters[3].py(), tau_daughters[3].pz(), tau_daughters[3].e() };
1070 double XA[4] = { tau_daughters[2].px(), tau_daughters[2].py(), tau_daughters[2].pz(), tau_daughters[2].e() };
1071 double QP[4] = { tau_daughters[1].px(), tau_daughters[1].py(), tau_daughters[1].pz(), tau_daughters[1].e() };
1072 double XN[4] = { tau_daughters[0].px(), tau_daughters[0].py(), tau_daughters[0].pz(), tau_daughters[0].e() };
1073 double AMPLIT = 0.0;
1074 double HV[4] = { 0.0 };
1078 QP[3] = sqrt( QP[0]*QP[0] + QP[1]*QP[1] + QP[2]*QP[2] + 0.105659*0.105659);
1079 XA[3] = sqrt( XA[0]*XA[0] + XA[1]*XA[1] + XA[2]*XA[2] );
1080 XK[3] = sqrt( XK[0]*XK[0] + XK[1]*XK[1] + XK[2]*XK[2] );
1082 if(XK0DEC > XK[3]/(XK[3]+XA[3]+QP[3]+XN[3])) XK0DEC=0.5*XK[3]/(XK[3]+XA[3]+QP[3]+XN[3]);
1085 dampry_( &ITDKRC, &XK0DEC, XK, XA, QP, XN, &LIT, HV );
1096 else if( pdgid.size()==4 &&
1098 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 111, 111,-211) ) ||
1099 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 111, 111, 211) )
1102 DEBUG( cout<<
"Channel 5 : "; )
1107 const double AMTAU = 1.777;
1109 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1110 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1111 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1112 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1113 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1115 float HV[4] = { 0.0 };
1120 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1131 else if( pdgid.size()==4 &&
1133 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211,-211, 211) ) ||
1134 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 211,-211) )
1137 DEBUG( cout<<
"Channel 5 : "; )
1142 const double AMTAU = 1.777;
1144 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1145 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1146 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1147 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1148 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1150 float HV[4] = { 0.0 };
1155 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1190 else if( pdgid.size()==4 &&
1192 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, -321, -211, 321) ) ||
1193 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 211,-321) )
1196 DEBUG( cout<<
"Channel 5 : "; )
1201 const double AMTAU = 1.777;
1205 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1206 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1207 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1208 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1209 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1211 float HV[4] = { 0.0 };
1216 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1227 else if( pdgid.size()==4 &&
1229 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 311, -211, 311 ) ) ||
1230 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 311, 211, 311 ) ) ||
1231 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 311, -211, 310 ) ) ||
1232 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 311, 211, 310 ) ) ||
1233 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 311, -211, 130 ) ) ||
1234 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 311, 211, 130 ) ) ||
1235 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 310, -211, 311 ) ) ||
1236 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 310, 211, 311 ) ) ||
1237 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 310, -211, 310 ) ) ||
1238 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 310, 211, 310 ) ) ||
1239 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 310, -211, 130 ) ) ||
1240 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 310, 211, 130 ) ) ||
1241 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 130, -211, 311 ) ) ||
1242 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 130, 211, 311 ) ) ||
1243 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 130, -211, 310 ) ) ||
1244 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 130, 211, 310 ) ) ||
1245 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 130, -211, 130 ) ) ||
1246 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 130, 211, 130 ) )
1249 DEBUG( cout<<
"Channel 5 : "; )
1254 const double AMTAU = 1.777;
1257 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1258 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1259 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1260 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1261 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1263 float HV[4] = { 0.0 };
1268 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1279 else if( pdgid.size()==4 &&
1281 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, -321, 311, 111) ) ||
1282 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 311, 111) ) ||
1283 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, -321, 310, 111) ) ||
1284 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 310, 111) ) ||
1285 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, -321, 130, 111) ) ||
1286 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 130, 111) )
1289 DEBUG( cout<<
"Channel 5 : "; )
1294 const double AMTAU = 1.777;
1297 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1298 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1299 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1300 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1301 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1303 float HV[4] = { 0.0 };
1308 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1319 else if( pdgid.size()==4 &&
1321 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 111, 111,-321) ) ||
1322 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 111, 111, 321) )
1325 DEBUG( cout<<
"Channel 5 : "; )
1330 const double AMTAU = 1.777;
1333 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1334 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1335 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1336 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1337 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1339 float HV[4] = { 0.0 };
1344 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1354 else if( pdgid.size()==4 &&
1356 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-321,-211, 211) ) ||
1357 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 321, 211,-211) )
1360 DEBUG( cout<<
"Channel 5 : "; )
1365 const double AMTAU = 1.777;
1368 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1369 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1370 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1371 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1372 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1374 float HV[4] = { 0.0 };
1379 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1388 else if( pdgid.size()==4 &&
1390 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 311, 111) ) ||
1391 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 311, 111) ) ||
1392 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 310, 111) ) ||
1393 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 310, 111) ) ||
1394 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211, 130, 111) ) ||
1395 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 130, 111) )
1398 DEBUG( cout<<
"Channel 5 : "; )
1403 const double AMTAU = 1.777;
1406 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1407 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1408 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1409 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1410 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1412 float HV[4] = { 0.0 };
1417 damppk_( &MNUM, PT, PN, PIM1, PIM2, PIPL, &LIT, HV );
1427 else if( pdgid.size()==5 &&
1429 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16,-211,-211, 211, 111) ) ||
1430 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 211, 211,-211, 111) )
1433 DEBUG( cout<<
"Channel 8 : "; )
1436 const double AMTAU = 1.777;
1438 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1439 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1440 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1441 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1442 float PIZ [4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1443 float PIPL[4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1445 float HV[4] = { 0.0 };
1447 dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &LIT, HV );
1457 else if( pdgid.size()==5 &&
1459 ( tau_pdgid== 15 &&
channelMatch(tau_daughters, 16, 111, 111, 111,-211) ) ||
1460 ( tau_pdgid==-15 &&
channelMatch(tau_daughters,-16, 111, 111, 111, 211) )
1463 DEBUG( cout<<
"Channel 9 : "; )
1466 const double AMTAU = 1.777;
1468 float PT[4] = { 0.0, 0.0, 0.0, (float)AMTAU };
1469 float PN[4] = { (float)tau_daughters[0].px(), (float)tau_daughters[0].py(), (float)tau_daughters[0].pz(), (float)tau_daughters[0].e() };
1470 float PIM1[4] = { (float)tau_daughters[1].px(), (float)tau_daughters[1].py(), (float)tau_daughters[1].pz(), (float)tau_daughters[1].e() };
1471 float PIM2[4] = { (float)tau_daughters[2].px(), (float)tau_daughters[2].py(), (float)tau_daughters[2].pz(), (float)tau_daughters[2].e() };
1472 float PIZ [4] = { (float)tau_daughters[3].px(), (float)tau_daughters[3].py(), (float)tau_daughters[3].pz(), (float)tau_daughters[3].e() };
1473 float PIPL[4] = { (float)tau_daughters[4].px(), (float)tau_daughters[4].py(), (float)tau_daughters[4].pz(), (float)tau_daughters[4].e() };
1475 float HV[4] = { 0.0 };
1477 dam4pi_( &MNUM, PT, PN, PIM1, PIM2, PIZ, PIPL, &LIT, HV );
1487 DEBUG( cout<<tau_daughters.size()<<
"-part ???: "; )
1492 Particle HHbuf(HH[0], HH[1], HH[2], HH[3], 0);
1494 HHbuf.rotateXZ(theta);
1495 HHbuf.rotateXY(phi);
1517 Particle tau_plus ( sp_tau.px(), sp_tau.py(), sp_tau.pz(), sp_tau.e(), sp_tau.pdgid() );
1518 Particle tau_minus( sp_nu_tau.px(), sp_nu_tau.py(), sp_nu_tau.pz(), sp_nu_tau.e(), sp_nu_tau.pdgid() );
1520 Particle P_QQ( tau_plus.px()+tau_minus.px(), tau_plus.py()+tau_minus.py(), tau_plus.pz()+tau_minus.pz(), tau_plus.e()+tau_minus.e(), 0 );
1525 tau_plus. boostToRestFrame(P_QQ);
1526 tau_minus.boostToRestFrame(P_QQ);
1527 P_B1. boostToRestFrame(P_QQ);
1528 P_B2. boostToRestFrame(P_QQ);
1530 double costheta1 = (tau_plus.px()*P_B1.px() +tau_plus.py()*P_B1.py() +tau_plus.pz()*P_B1.pz() ) /
1531 sqrt(tau_plus.px()*tau_plus.px()+tau_plus.py()*tau_plus.py()+tau_plus.pz()*tau_plus.pz()) /
1532 sqrt(P_B1.px() *P_B1.px() +P_B1.py() *P_B1.py() +P_B1.pz() *P_B1.pz() );
1534 double costheta2 = (tau_minus.px()*P_B2.px() +tau_minus.py()*P_B2.py() +tau_minus.pz()*P_B2.pz() ) /
1535 sqrt(tau_minus.px()*tau_minus.px()+tau_minus.py()*tau_minus.py()+tau_minus.pz()*tau_minus.pz()) /
1536 sqrt(P_B2.px() *P_B2.px() +P_B2.py() *P_B2.py() +P_B2.pz() *P_B2.pz() );
1538 double sintheta1 = sqrt(1-costheta1*costheta1);
1539 double sintheta2 = sqrt(1-costheta2*costheta2);
1542 double costhe = (costheta1*sintheta2 + costheta2*sintheta1) / (sintheta1 + sintheta2);
1556 double x1x2 = SS/CMSENE/CMSENE;
1557 double x1Mx2 = P_QQ.pz()/CMSENE*2;
1559 double x1 = ( x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
1560 double x2 = ( -x1Mx2 + sqrt(x1Mx2*x1Mx2 + 4*x1x2) )/2;
1563 WID[0] = f(x1, 0,SS,CMSENE)*f(x2, 0,SS,CMSENE) * sigborn(0,SS, costhe);
1564 WID[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe);
1565 WID[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhe);
1566 WID[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe);
1567 WID[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhe);
1568 WID[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe);
1569 WID[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhe);
1570 WID[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe);
1571 WID[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhe);
1572 WID[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe);
1573 WID[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhe);
1576 for(
int i=0;i<=10;i++) sum+=WID[i];
1580 cout <<
"Tauspinner::calculateWeightFromParticlesH WARNING: sum of WID[0]-WID[10] is 0. Check LHAPDF configuration" << endl;
1584 if(relWTnonSM==0) WTnonSM=sum;
1588 WID2[0] = f(x1, 0,SS,CMSENE)*f(x2, 0,SS,CMSENE) * sigborn(0,SS, costhe) *
plweight(0,SS, costhe);
1589 WID2[1] = f(x1, 1,SS,CMSENE)*f(x2,-1,SS,CMSENE) * sigborn(1,SS, costhe) *
plweight(1,SS, costhe);
1590 WID2[2] = f(x1,-1,SS,CMSENE)*f(x2, 1,SS,CMSENE) * sigborn(1,SS,-costhe) *
plweight(1,SS,-costhe);
1591 WID2[3] = f(x1, 2,SS,CMSENE)*f(x2,-2,SS,CMSENE) * sigborn(2,SS, costhe) *
plweight(2,SS, costhe);
1592 WID2[4] = f(x1,-2,SS,CMSENE)*f(x2, 2,SS,CMSENE) * sigborn(2,SS,-costhe) *
plweight(2,SS,-costhe);
1593 WID2[5] = f(x1, 3,SS,CMSENE)*f(x2,-3,SS,CMSENE) * sigborn(3,SS, costhe) *
plweight(3,SS, costhe);
1594 WID2[6] = f(x1,-3,SS,CMSENE)*f(x2, 3,SS,CMSENE) * sigborn(3,SS,-costhe) *
plweight(3,SS,-costhe);
1595 WID2[7] = f(x1, 4,SS,CMSENE)*f(x2,-4,SS,CMSENE) * sigborn(4,SS, costhe) *
plweight(4,SS, costhe);
1596 WID2[8] = f(x1,-4,SS,CMSENE)*f(x2, 4,SS,CMSENE) * sigborn(4,SS,-costhe) *
plweight(4,SS,-costhe);
1597 WID2[9] = f(x1, 5,SS,CMSENE)*f(x2,-5,SS,CMSENE) * sigborn(5,SS, costhe) *
plweight(5,SS, costhe);
1598 WID2[10]= f(x1,-5,SS,CMSENE)*f(x2, 5,SS,CMSENE) * sigborn(5,SS,-costhe) *
plweight(5,SS,-costhe);
1601 for(
int i=0;i<=10;i++) sum2+=WID2[i];
1604 if(relWTnonSM==0) WTnonSM=sum2;
1612 if(IfHiggs && nonSM2==1) {
1613 double polp =
plzap2(0,15,S,costhe);
1614 pol += (2*(1-polp)-1);
1617 if(IfHiggs)
return NAN;
1620 for(
int i=0;i<=10;i++) WID[i]/=sum;
1628 for(
int i=1;i<=10;i++)
1632 double cost = costhe;
1635 if( ICC==2 || ICC==4 || ICC==6 || ICC==8 || ICC==10 ) cost = -cost;
1639 if( ICC==7 || ICC==8 ) ID = 4;
1640 if( ICC==1 || ICC==2 ) ID = 1;
1641 if( ICC==5 || ICC==6 ) ID = 3;
1642 if( ICC==9 || ICC==10 ) ID = 5;
1646 double polp =
plzap2(ID,tau_pdgid,S,cost);
1647 pol += (2*(1-polp)-1)*WID[i];
1654 pp.
m_R[1][1] = pp.
m_R[2][2] = 0.0;
1659 R11 += WID[i]*pp.
m_R[1][1];
1660 R22 += WID[i]*pp.
m_R[2][2];
1681 bool channelMatch(vector<Particle> &particles,
int p1,
int p2,
int p3,
int p4,
int p5,
int p6)
1686 for(
unsigned int i=0;i<particles.size();i++) list.push_back(particles[i].pdgid());
1689 int p[6] = { p1, p2, p3, p4, p5, p6 };
1693 for(
int i=0;i<6;i++)
1700 for(
unsigned int j=0;j<list.size(); j++)
1706 list.erase(list.begin()+j);
1711 if(!found)
return false;
1715 if(list.size()!=0)
return false;
1720 vector<Particle> newList;
1722 for(
int i=0;i<6;i++)
1727 for(
unsigned int j=0;j<particles.size(); j++)
1730 if(particles[j].pdgid()==p[i])
1732 newList.push_back(particles[j]);
1733 particles.erase(particles.begin()+j);
1739 particles = newList;
1754 double px=nu_tau.px()+tau.px();
1755 double py=nu_tau.py()+tau.py();
1756 double pz=nu_tau.pz()+tau.pz();
1757 double e =nu_tau.e ()+tau.e ();
1760 cout<<
"--------------------------------------------------------------------------------------------------------"<<endl;
1768 for(
unsigned int i=0; i<tau_daughters.size();i++) tau_daughters[i].
print();
1771 cout<<
"--------------------------------------------------------------------------------------------------------"<<endl;
1785 double px=0.0,py=0.0,pz=0.0,e=0.0;
1787 for(
unsigned int i=0; i<x.size();i++)
double calculateWeightFromParticlesH(SimpleParticle &sp_X, SimpleParticle &sp_tau1, SimpleParticle &sp_tau2, vector< SimpleParticle > &sp_tau1_daughters, vector< SimpleParticle > &sp_tau2_daughters)
double calculateWeightFromParticlesWorHpn(SimpleParticle &W, SimpleParticle &tau, SimpleParticle &nu_tau, vector< SimpleParticle > &tau_daughters)
void print(Particle &W, Particle &nu_tau, Particle &tau, vector< Particle > &tau_daughters)
double disth_(double *SVAR, double *COSTHE, int *TA, int *TB)
void setHiggsParameters(int jak, double mass, double width, double normalization)
void getZgamParametersTR(double &Rxx, double &Ryy, double &Rxy, double &Ryx)
double plzap2(int ide, int idf, double svar, double costhe)
void setRelWTnonSM(int _relWTnonSM)
void setZgamMultipliersTR(double Rxx, double Ryy, double Rxy, double Ryx)
void recalculateRij(int incoming_pdg_id, int outgoing_pdg_id, double invariant_mass_squared, double cosTheta)
void prepareKinematicForHH(Particle &tau, Particle &nu_tau, vector< Particle > &tau_daughters, double *phi2, double *theta2)
double plweight(int ide, double svar, double costhe)
double getLongitudinalPolarization(double, SimpleParticle &, SimpleParticle &)
double * calculateHH(int tau_pdgid, vector< Particle > &tau_daughters, double phi, double theta)
void setNonSMkey(int key)
void getHiggsParameters(double *mass, double *width, double *normalization)
bool channelMatch(vector< Particle > &particles, int p1, int p2=0, int p3=0, int p4=0, int p5=0, int p6=0)
Particle * vector_sum(vector< Particle > &x)
void setHiggsParametersTR(double Rxx, double Ryy, double Rxy, double Ryx)
void initialize_spinner(bool _Ipp, int _Ipol, int _nonSM2, int _nonSMN, double _CMSENE)