12 const float minPFRecHitFrac = 1E-4;
24 std::cout <<
"entering test link by rechit function" << std::endl;
30 auto clusterZ = cluster.
position().Z();
35 double horesolscale = 1.0;
42 double tracketa = 999999.9;
43 double trackphi = 999999.9;
44 double track_X = 999999.9;
45 double track_Y = 999999.9;
46 double track_Z = 999999.9;
54 switch (cluster.
layer()) {
61 std::cout <<
"Fetching Ecal Resolution Maps" << std::endl;
87 std::cout <<
"Fetching Hcal Resolution Maps" << std::endl;
103 dHPhi = dHPhi - 2. *
M_PI;
104 else if (dHPhi < -
M_PI)
105 dHPhi = dHPhi + 2. *
M_PI;
106 tracketa = atHCAL.
positionREP().Eta() + 0.1 * dHEta;
107 trackphi = atHCAL.
positionREP().Phi() + 0.1 * dHPhi;
124 std::cout <<
"Fetching HO Resolution Maps" << std::endl;
144 if (fabs(track_Z) > 700.25)
172 std::cout <<
"No link by rechit possible for pre-shower yet!" << std::endl;
186 if (!hcal && fabs(track_Z) < 300.)
188 if (track_Z * clusterZ < 0.)
194 if (!hcal && fabs(track_Z) > 300.)
202 std::cout <<
"test link by rechit " << dist <<
" " << std::endl;
204 std::cout <<
" clustereta " << clustereta <<
" clusterphi " << clusterphi <<
" tracketa " << tracketa
205 <<
" trackphi " << trackphi << std::endl;
214 const std::vector<reco::PFRecHitFraction>& fracs = cluster.
recHitFractions();
216 bool linkedbyrechit =
false;
218 for (
unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
220 double fraction = fracs[rhit].fraction();
221 if (fraction < minPFRecHitFrac)
227 const auto& rechit_cluster = *rh;
228 const auto& posxyz = rechit_cluster.position();
229 const auto& posrep = rechit_cluster.positionREP();
232 const auto& cornersxyz = rechit_cluster.getCornersXYZ();
233 const auto& corners = rechit_cluster.getCornersREP();
235 if (barrel || hcal) {
243 if (rhsizePhi >
M_PI)
244 rhsizePhi = 2. *
M_PI - rhsizePhi;
246 const double mult = horesolscale * (1.50 + 0.5 / fracs.size());
247 rhsizeEta = rhsizeEta * mult + 0.2 *
std::abs(dHEta);
248 rhsizePhi = rhsizePhi * mult + 0.2 * fabs(dHPhi);
251 const double mult = 2.00 + 1.0 / (fracs.size() *
std::min(1., 0.5 * trackPt));
258 std::cout << rhit <<
" Hcal RecHit=" << posrep.Eta() <<
" " << posrep.Phi() <<
" " << rechit_cluster.energy()
260 for (
unsigned jc = 0; jc < 4; ++jc)
261 std::cout <<
"corners " << jc <<
" " << corners[jc].
eta() <<
" " << corners[jc].phi() << std::endl;
263 std::cout <<
"RecHit SizeEta=" << rhsizeEta <<
" SizePhi=" << rhsizePhi << std::endl;
270 double deta = fabs(posrep.eta() - tracketa);
271 double dphi = fabs(posrep.phi() - trackphi);
273 dphi = 2. *
M_PI - dphi;
277 std::cout <<
"distance=" << deta <<
" " << dphi <<
" ";
278 if (deta < (0.5 * rhsizeEta) && dphi < (0.5 * rhsizePhi))
279 std::cout <<
" link here !" << std::endl;
285 if (deta < (0.5 * rhsizeEta) && dphi < (0.5 * rhsizePhi)) {
286 linkedbyrechit =
true;
293 const auto& posxyz = rechit_cluster.position();
295 std::cout <<
"RH " << posxyz.x() <<
" " << posxyz.y() << std::endl;
297 std::cout <<
"TRACK " << track_X <<
" " << track_Y << std::endl;
304 for (
unsigned jc = 0; jc < 4; ++jc) {
305 const auto& cornerposxyz = cornersxyz[jc];
306 const double mult = (1.00 + 0.50 / (fracs.size() *
std::min(1., 0.5 * trackPt)));
307 x[3 - jc] = cornerposxyz.x() + (cornerposxyz.x() - posxyz.x()) * mult;
308 y[3 - jc] = cornerposxyz.y() + (cornerposxyz.y() - posxyz.y()) * mult;
312 std::cout <<
"corners " << jc <<
" " << cornerposxyz.x() <<
" " << cornerposxyz.y() << std::endl;
324 bool isinside = TMath::IsInside(track_X, track_Y, 5, x, y);
327 linkedbyrechit =
true;
334 if (linkedbyrechit) {
337 std::cout <<
"Track and Cluster LINKED BY RECHIT" << std::endl;
368 static const double resPSpitch = 0.19 /
sqrt(12.);
369 static const double resPSlength = 6.1 /
sqrt(12.);
378 std::cout <<
"entering test link by rechit function for ECAL and PS" << std::endl;
382 double zECAL = clusterECAL.
position().Z();
383 double xECAL = clusterECAL.
position().X();
384 double yECAL = clusterECAL.
position().Y();
387 double zPS = clusterPS.
position().Z();
388 double xPS = clusterPS.
position().X();
389 double yPS = clusterPS.
position().Y();
391 if (zECAL * zPS < 0.)
396 switch (clusterPS.
layer()) {
399 deltaX = resPSpitch * sqr12;
400 deltaY = resPSlength * sqr12;
404 deltaY = resPSpitch * sqr12;
405 deltaX = resPSlength * sqr12;
411 auto deltaXY =
BVector2D(deltaX, deltaY).
v * 0.5;
413 auto zCorr = zPS / zECAL;
414 const std::vector<reco::PFRecHitFraction>& fracs = clusterECAL.
recHitFractions();
415 bool linkedbyrechit =
false;
417 for (
unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
418 const auto& rh = fracs[rhit].recHitRef();
419 double fraction = fracs[rhit].fraction();
420 if (fraction < minPFRecHitFrac)
434 std::cout <<
"Ecal rechit " << posxy.x() <<
" " << posxy.y() << std::endl;
435 std::cout <<
"PS cluster " << xPS <<
" " << yPS << std::endl;
441 for (
unsigned jc = 0; jc < 4; ++jc) {
444 auto dist = (cornerpos - posxy);
448 auto xy = cornerpos + (dist * (fivepercent2D + one2D / adist) * deltaXY);
460 std::cout <<
"corners " << jc <<
" " << cornerpos.x() <<
" " << x[3 - jc] <<
" " << cornerpos.y() <<
" "
461 << y[3 - jc] << std::endl;
473 bool isinside = TMath::IsInside(xPS, yPS, 5, x, y);
476 linkedbyrechit =
true;
482 if (linkedbyrechit) {
485 std::cout <<
"Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
487 constexpr
double scale = 1. / 1000.;
488 double dist =
computeDist(xECAL * scale, yECAL * scale, xPS * scale, yPS * scale,
false);
498 const auto& posxyzEM = clusterHFEM.
position();
499 const auto& posxyzHAD = clusterHFHAD.
position();
501 double sameZ = posxyzEM.Z() * posxyzHAD.Z();
505 double dX = posxyzEM.X() - posxyzHAD.X();
506 double dY = posxyzEM.Y() - posxyzHAD.Y();
508 double dist2 = dX * dX + dY * dY;
512 double dist =
sqrt(dist2);
520 auto phicor = etaPhi ?
normalizedPhi(phi1 - phi2) : phi1 - phi2;
521 auto etadiff = eta1 -
eta2;
527 return std::sqrt(etadiff * etadiff + phicor * phicor);
PFLayer::Layer layer() const
cluster layer, see PFLayer.h in this directory
const REPPoint & positionREP() const
trajectory position in (rho, eta, phi) base
const math::XYZPoint & position() const
cluster centroid position
reconstructed track used as an input to particle flow
Basic2DVector< T > xy() const
static double computeDist(double eta1, double phi1, double eta2, double phi2, bool etaPhi=true)
computes a chisquare
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
static double testECALAndPSByRecHit(const reco::PFCluster &clusterECAL, const reco::PFCluster &clusterPS, bool debug=false)
const math::XYZPoint & position() const
cartesian position (x, y, z)
constexpr T normalizedPhi(T phi)
PositionType const & position() const
rechit cell centre x, y, z
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Basic2DVector< double >::MathVector Vector2D
Point of closest approach from beam axis (initial point in the case of PFSimParticle) ...
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Abs< T >::type abs(const T &t)
const reco::PFTrajectoryPoint & extrapolatedPoint(unsigned layerid) const
bool isNull() const
Checks for null.
Basic2DVector< T > xy() const
const math::XYZTLorentzVector & momentum() const
4-momenta quadrivector
bool isValid() const
is this point valid ?
static double testHFEMAndHFHADByRecHit(const reco::PFCluster &clusterHFEM, const reco::PFCluster &clusterHFHAD, bool debug=false)
test association between HFEM and HFHAD, by rechit
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Basic2DVector< double > BVector2D
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
static double testTrackAndClusterByRecHit(const reco::PFRecTrack &track, const reco::PFCluster &cluster, bool isBrem=false, bool debug=false)
CornersVec const & getCornersXYZ() const
rechit corners