764double*
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;
1656 pp.recalculateRij(ID,tau_pdgid,S,cost);
1659 R11 += WID[i]*pp.m_R[1][1];
1660 R22 += WID[i]*pp.m_R[2][2];