1
- /* From SimpleFits Package
2
- * Designed an written by
3
- * author: Ian M. Nugent
4
- * Humboldt Foundations
5
- */
6
1
#include " RecoTauTag/ImpactParameter/interface/MultiProngTauSolver.h"
7
- #include < iostream>
8
2
#include " TMatrixTSym.h"
3
+ #include " TVectorT.h"
4
+ #include " FWCore/MessageLogger/interface/MessageLogger.h"
9
5
10
6
using namespace tauImpactParameter ;
11
7
12
- void MultiProngTauSolver::quadratic (double & x_plus,double & x_minus,double a, double b, double c){
8
+ void MultiProngTauSolver::quadratic (double & x_plus,double & x_minus,double a, double b, double c, bool &isReal ){
13
9
double R=b*b-4 *a*c;
14
- if (R<0 ){R=0 ;}
15
- x_minus=(-b+sqrt (R))/(2.0 *a); // opposite sign is smaller
16
- x_plus=(-b-sqrt (R))/(2.0 *a);
10
+ isReal=true ;
11
+ if (R<0 ){isReal=false ;}// flag cases when R<0 but compute quadratic equation with |R|
12
+ x_minus=(-b+sqrt (fabs (R)))/(2.0 *a); // opposite sign is smaller
13
+ x_plus=(-b-sqrt (fabs (R)))/(2.0 *a);
17
14
}
18
15
19
- void MultiProngTauSolver::analyticESolver (TLorentzVector& nu_plus,TLorentzVector& nu_minus, const TLorentzVector& A1 ){
16
+ void MultiProngTauSolver::analyticESolver (TLorentzVector& nu_plus,TLorentzVector& nu_minus,const TLorentzVector &A1, bool &isReal ){
20
17
double a=(A1.Pz ()*A1.Pz ())/(A1.E ()*A1.E ())-1.0 ;
21
18
double K=(PDGInfo::tau_mass ()*PDGInfo::tau_mass ()-A1.M2 ()-2.0 *A1.Pt ()*A1.Pt ())/(2.0 *A1.E ());
22
19
double b=2.0 *K*A1.Pz ()/A1.E ();
23
20
double c=K*K-A1.Pt ()*A1.Pt ();
24
21
double z_plus (0 ),z_minus (0 );
25
- quadratic (z_plus,z_minus,a,b,c);
26
- nu_plus. SetPxPyPzE (-A1.Px (),-A1.Py (),z_plus,sqrt (z_plus*z_plus+A1.Pt ()*A1.Pt ()));
27
- nu_minus. SetPxPyPzE (-A1.Px (),-A1.Py (),z_minus,sqrt (z_minus*z_minus+A1.Pt ()*A1.Pt ()));
22
+ quadratic (z_plus,z_minus,a,b,c,isReal );
23
+ nu_plus= TLorentzVector (-A1.Px (),-A1.Py (),z_plus,sqrt (z_plus*z_plus+A1.Pt ()*A1.Pt ()));
24
+ nu_minus= TLorentzVector (-A1.Px (),-A1.Py (),z_minus,sqrt (z_minus*z_minus+A1.Pt ()*A1.Pt ()));
28
25
}
29
26
30
- void MultiProngTauSolver::numericalESolver (TLorentzVector& nu_plus,TLorentzVector& nu_minus, const TLorentzVector& A1){
27
+ void MultiProngTauSolver::numericalESolver (TLorentzVector& nu_plus,TLorentzVector& nu_minus,const TLorentzVector& A1, bool &isReal ){
31
28
double rmin (-100 ), rmax (100 ), step (0.01 ), mtau2 (PDGInfo::tau_mass ()*PDGInfo::tau_mass ()), z1 (-999 ), z2 (-999 ), zmin (-999 ), min (9999 ), prev (9999 );
32
29
double z=rmin;
33
30
TLorentzVector nu,tau;
@@ -42,25 +39,24 @@ void MultiProngTauSolver::numericalESolver(TLorentzVector& nu_plus,TLorentzVecto
42
39
z+=step;
43
40
}
44
41
if (z1!=-999 && z2!=-999 ){
45
- nu_plus. SetPxPyPzE (-A1.Px (),-A1.Py (),z1,sqrt (z1*z1+A1.Pt ()*A1.Pt ()));
46
- nu_minus. SetPxPyPzE (-A1.Px (),-A1.Py (),z2,sqrt (z2*z2+A1.Pt ()*A1.Pt ()));
42
+ nu_plus= TLorentzVector (-A1.Px (),-A1.Py (),z1,sqrt (z1*z1+A1.Pt ()*A1.Pt ()));
43
+ nu_minus= TLorentzVector (-A1.Px (),-A1.Py (),z2,sqrt (z2*z2+A1.Pt ()*A1.Pt ()));
47
44
}
48
45
else {
49
- nu_plus. SetPxPyPzE (-A1.Px (),-A1.Py (),zmin,sqrt (zmin*zmin+A1.Pt ()*A1.Pt ()));
50
- nu_minus. SetPxPyPzE (-A1.Px (),-A1.Py (),zmin,sqrt (zmin*zmin+A1.Pt ()*A1.Pt ()));
46
+ nu_plus= TLorentzVector (-A1.Px (),-A1.Py (),zmin,sqrt (zmin*zmin+A1.Pt ()*A1.Pt ()));
47
+ nu_minus= TLorentzVector (-A1.Px (),-A1.Py (),zmin,sqrt (zmin*zmin+A1.Pt ()*A1.Pt ()));
51
48
}
52
49
}
53
50
54
51
void MultiProngTauSolver::solveByRotation (const TVector3& TauDir,const TLorentzVector& A1, TLorentzVector& Tau_plus,TLorentzVector& Tau_minus,
55
- TLorentzVector& nu_plus,TLorentzVector& nu_minus, bool rotateback){
52
+ TLorentzVector & nu_plus,TLorentzVector & nu_minus, bool &isReal, bool rotateback){
56
53
TLorentzVector A1rot=A1;
57
54
double phi (TauDir.Phi ()),theta (TauDir.Theta ());
58
55
A1rot.RotateZ (-phi);
59
56
A1rot.RotateY (-theta);
60
57
// ///////////////////////////////////////////////////
61
- // numericalESolver(nu_plus,nu_minus,A1rot); // for debugging AnalyticESolver (slow)
62
- // analyticESolver(nu_plus,nu_minus,A1rot);
63
- analyticESolver (nu_plus,nu_minus,A1rot);
58
+ // NumericalESolver(nu_plus,nu_minus,A1rot); // for debugging analyticESolver (slow)
59
+ analyticESolver (nu_plus,nu_minus,A1rot,isReal);
64
60
// ///////////////////////////////////////////////////
65
61
if (rotateback){
66
62
nu_plus.RotateY (theta);
@@ -78,13 +74,21 @@ void MultiProngTauSolver::solveByRotation(const TVector3& TauDir,const TLorentzV
78
74
}
79
75
80
76
bool MultiProngTauSolver::setTauDirectionatThetaGJMax (const TLorentzVector& a1, double & theta,double & phi,double scale){
81
- double thetaGJMax_a1 =thetaGJMax (a1);
82
- double dtheta=(theta-a1.Theta ());
83
- double dphi=fmod (fabs (phi-a1.Phi ()),2 *TMath::Pi ());if (phi<a1.Phi ())dphi*=-1.0 ;
84
- double dphitheta=sqrt (dtheta*dtheta+dphi*dphi);
85
- if (thetaGJMax_a1<dphitheta || scale<0 ){
86
- theta=a1.Theta ()+dtheta*(thetaGJMax_a1/dphitheta)*fabs (scale);
87
- phi=a1.Phi ()+dphi*(thetaGJMax_a1/dphitheta)*fabs (scale);
77
+ double thetaGJMaxvar =thetaGJMax (a1);
78
+ TVector3 a1v (a1.Vect ()); if (a1v.Mag ()!=0 ) a1v*=1 /a1v.Mag ();
79
+ TVector3 tau (cos (phi)*sin (theta),sin (phi)*sin (theta),cos (theta));
80
+ double dphitheta=acos (a1v.Dot (tau)/(a1v.Mag ()*tau.Mag ()));
81
+ if (thetaGJMaxvar<dphitheta || scale<0 ){
82
+ if (scale<0 ) scale=1.0 ;
83
+ double a=(thetaGJMaxvar/dphitheta)-(1 -scale);
84
+ double b=1 -(thetaGJMaxvar/dphitheta)+(1 -scale);
85
+ edm::LogInfo (" RecoTauTag/ImpactParameter" ) << " SetTauDirectionatThetaGJMax before GF " << thetaGJMaxvar << " dot " << acos (a1v.Dot (tau)/(a1v.Mag ()*tau.Mag ())) << " a1 phi " << a1v.Phi () << " tau phi " << tau.Phi () << " a1 theta " <<a1v.Theta () << " tau theta " << tau.Theta () ;
86
+ tau*=a;
87
+ a1v*=b;
88
+ tau+=a1v;
89
+ theta=tau.Theta ();
90
+ phi=tau.Phi ();
91
+ edm::LogInfo (" RecoTauTag/ImpactParameter" ) << " SetTauDirectionatThetaGJMax GF " << thetaGJMaxvar << " dot " << acos (a1v.Dot (tau)/(a1v.Mag ()*tau.Mag ())) << " phi " << phi << " theta " << theta ;
88
92
return true ;
89
93
}
90
94
return false ;
@@ -105,14 +109,15 @@ LorentzVectorParticle MultiProngTauSolver::estimateNu(const LorentzVectorParticl
105
109
double theta=tauFlghtDir.Theta ();
106
110
double phi=tauFlghtDir.Phi ();
107
111
setTauDirectionatThetaGJMax (lorentzA1,theta,phi);
108
- startingtauFlghtDir. SetMagThetaPhi ( 1.0 , theta, phi);
112
+ startingtauFlghtDir= TVector3 ( sin ( theta)* cos (phi), sin (theta)* sin ( phi), cos (theta) );
109
113
}
110
114
TLorentzVector tau1,tau2,nu1,nu2;
111
- solveByRotation (startingtauFlghtDir,lorentzA1,tau1,tau2,nu1,nu2);
115
+ bool isReal;
116
+ solveByRotation (startingtauFlghtDir,lorentzA1,tau1,tau2,nu1,nu2,isReal);
112
117
if (ambiguity==plus){ nuGuess=nu1; tau=tau1; }
113
118
if (ambiguity==minus){ nuGuess=nu2; tau=tau1; }
114
119
if (ambiguity==zero){ nuGuess=nu1; tau=tau1; }
115
- TVectorT<double > par (LorentzVectorParticle::NLorentzandVertexPar, 10 );
120
+ TVectorT<double > par (LorentzVectorParticle::NLorentzandVertexPar);
116
121
par (LorentzVectorParticle::vx)=a1.parameter (LorentzVectorParticle::vx);
117
122
par (LorentzVectorParticle::vy)=a1.parameter (LorentzVectorParticle::vy);
118
123
par (LorentzVectorParticle::vz)=a1.parameter (LorentzVectorParticle::vz);
@@ -135,11 +140,19 @@ LorentzVectorParticle MultiProngTauSolver::estimateNu(const LorentzVectorParticl
135
140
return LorentzVectorParticle (par,Cov,PDGInfo::nu_tau,0 ,a1.bField ());
136
141
}
137
142
138
- TVectorT<double > MultiProngTauSolver::rotateToTauFrame (const TVectorT<double >& inpar){
143
+ TVectorT<double > MultiProngTauSolver::rotateToTauFrame (const TVectorT<double > & inpar){
139
144
TVectorT<double > outpar (3 );
140
145
TVector3 res (inpar (0 ),inpar (1 ),inpar (2 ));
141
- TVector3 Uz;Uz. SetMagThetaPhi ( 1 , inpar (4 ), inpar (3 ));
146
+ TVector3 Uz ( sin ( inpar (4 ))* cos ( inpar ( 3 )), sin ( inpar (4 ))* sin ( inpar ( 3 )), cos ( inpar ( 4 ) ));
142
147
res.RotateUz (Uz);
148
+ /* double phi=inpar(3,0);
149
+ double theta=inpar(4,0);
150
+ res.RotateZ(-phi);
151
+ TVector3 Y(0,1,0);
152
+ TVector3 thetadir=res.Cross(Y);
153
+ thetadir.RotateY(-theta);
154
+ res.RotateY(-theta);
155
+ res.RotateZ(thetadir.Phi());*/
143
156
outpar (0 )=res.X ();
144
157
outpar (1 )=res.Y ();
145
158
outpar (2 )=res.Z ();
0 commit comments