17 if(
jet.m() == 0.0 )
return PF;
18 if(
jet.constituents().empty() )
return PF;
20 vector<fastjet::PseudoJet> constit_pseudojets =
jet.constituents();
22 TMatrixDSym MomentumTensor(2);
25 double eta0=
jet.eta();
33 double RotationMatrix[3][3];
35 for(
int i=0;
i<3;
i++) {
36 for(
int j=0; j<3; j++) {
37 RotationMatrix[
i][j] = 0.;
41 const double mag3 = sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1]+ nvec[2]*nvec[2]);
42 const double mag2 = sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1]);
49 const double inv_mag3 = 1. / mag3;
51 double ctheta0 = nvec[2]*inv_mag3;
52 double stheta0 =
mag2*inv_mag3;
53 double cphi0 = (
mag2>0.) ? nvec[0]/
mag2:0.;
54 double sphi0 = (
mag2>0.) ? nvec[1]/
mag2:0.;
56 RotationMatrix[0][0] =- ctheta0*cphi0;
57 RotationMatrix[0][1] =- ctheta0*sphi0;
58 RotationMatrix[0][2] = stheta0;
59 RotationMatrix[1][0] = sphi0;
60 RotationMatrix[1][1] =- 1.*cphi0;
61 RotationMatrix[1][2] = 0.;
62 RotationMatrix[2][0] = stheta0*cphi0;
63 RotationMatrix[2][1] = stheta0*sphi0;
64 RotationMatrix[2][2] = ctheta0;
69 const fastjet::PseudoJet &
cp = *cit;
70 TLorentzVector
p = TLorentzVector(
cp.px(),
cp.py(),
cp.pz(),
cp.e());
71 double n=1./(
cp.e()*
jet.m());
72 double px_rot = RotationMatrix[0][0] * (
p.Px())+RotationMatrix[0][1]
73 * (
p.Py())+RotationMatrix[0][2]*(
p.Pz());
74 double py_rot = RotationMatrix[1][0] * (
p.Px())+RotationMatrix[1][1]
75 * (
p.Py())+RotationMatrix[1][2]*(
p.Pz());
76 double pz_rot = RotationMatrix[2][0] * (
p.Px())+RotationMatrix[2][1]
77 * (
p.Py())+RotationMatrix[2][2]*(
p.Pz());
80 prot.SetPxPyPzE(px_rot, py_rot, pz_rot,
p.E() );
82 MomentumTensor(0,0) +=
n * prot.Px() * prot.Px();
83 MomentumTensor(0,1) +=
n * prot.Py() * prot.Px();
84 MomentumTensor(1,0) +=
n * prot.Px() * prot.Py();
85 MomentumTensor(1,1) +=
n * prot.Py() * prot.Py();
88 TMatrixDSymEigen eigen(MomentumTensor);
89 TVectorD
Lambda = eigen.GetEigenValues();
92 if ( fabs(den) < 1.e-20 )
return PF;