17 std::cout<<
"entering test link by rechit function"<<std::endl;
23 double clusterX = cluster.
position().X();
24 double clusterY = cluster.
position().Y();
25 double clusterZ = cluster.
position().Z();
29 double distance = 999999.9;
39 double tracketa = 999999.9;
40 double trackphi = 999999.9;
41 double track_X = 999999.9;
42 double track_Y = 999999.9;
43 double track_Z = 999999.9;
48 double trackPt = isBrem ? 999. :
sqrt(atVertex.
momentum().Vect().Perp2());
52 switch (cluster.
layer()) {
57 std::cout <<
"Fetching Ecal Resolution Maps"
61 if( ! atECAL.
isValid() )
return -1.;
70 =
std::sqrt( (track_X-clusterX)*(track_X-clusterX)
71 +(track_Y-clusterY)*(track_Y-clusterY)
72 +(track_Z-clusterZ)*(track_Z-clusterZ)
81 std::cout <<
"Fetching Hcal Resolution Maps"
93 if( ! atHCAL.
isValid() )
return -1.;
98 if ( dHPhi >
M_PI ) dHPhi = dHPhi - 2.*
M_PI;
99 else if ( dHPhi < -
M_PI ) dHPhi = dHPhi + 2.*
M_PI;
106 = -
std::sqrt( (track_X-clusterX)*(track_X-clusterX)
107 +(track_Y-clusterY)*(track_Y-clusterY)
108 +(track_Z-clusterZ)*(track_Z-clusterZ)
119 std::cout <<
"No link by rechit possible for pre-shower yet!"
135 if ( !hcal && fabs(track_Z) < 300. )
return -1.;
136 if ( track_Z * clusterZ < 0. )
return -1.;
141 if ( !hcal && fabs(track_Z) > 300. )
return -1.;
151 if(debug)
std::cout<<
"test link by rechit "<< dist <<
" "<<std::endl;
154 <<
" clusterphi " << clusterphi
155 <<
" tracketa " << tracketa
156 <<
" trackphi " << trackphi << std::endl;
165 const std::vector< reco::PFRecHitFraction >&
168 bool linkedbyrechit =
false;
170 for(
unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
173 double fraction = fracs[rhit].fraction();
174 if(fraction < 1E-4)
continue;
185 const std::vector< math::XYZPoint >&
187 const std::vector<reco::PFRecHit::REPPoint>& corners =
189 assert(corners.size() == 4);
191 if( barrel || hcal ){
198 = fabs(corners[0].
Eta() - corners[2].
Eta());
200 = fabs(corners[0].
Phi() - corners[2].
Phi());
201 if ( rhsizePhi >
M_PI ) rhsizePhi = 2.*
M_PI - rhsizePhi;
203 rhsizeEta = rhsizeEta * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHEta);
204 rhsizePhi = rhsizePhi * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHPhi);
207 rhsizeEta *= 2.00 + 1.0/fracs.size()/
std::min(1.,trackPt/2.);
208 rhsizePhi *= 2.00 + 1.0/fracs.size()/
std::min(1.,trackPt/2.);
214 << posrep.Eta() <<
" "
215 << posrep.Phi() <<
" "
216 << rechit_cluster.
energy()
218 for (
unsigned jc=0; jc<4; ++jc )
220 <<
" "<<corners[jc].Phi()<<std::endl;
222 std::cout <<
"RecHit SizeEta=" << rhsizeEta
223 <<
" SizePhi=" << rhsizePhi << std::endl;
230 double deta = fabs(posrep.Eta() - tracketa);
231 double dphi = fabs(posrep.Phi() - trackphi);
232 if ( dphi >
M_PI ) dphi = 2.*
M_PI - dphi;
239 if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.))
240 std::cout <<
" link here !" << std::endl;
245 if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.)){
246 linkedbyrechit =
true;
270 for (
unsigned jc=0; jc<4; ++jc ) {
272 x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X())
273 * (1.00+0.50/fracs.size()/
std::min(1.,trackPt/2.));
274 y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y())
275 * (1.00+0.50/fracs.size()/
std::min(1.,trackPt/2.));
280 <<
" " << cornerposxyz.X()
281 <<
" " << cornerposxyz.Y()
294 bool isinside = TMath::IsInside(track_X,
299 linkedbyrechit =
true;
306 if( linkedbyrechit ) {
309 std::cout <<
"Track and Cluster LINKED BY RECHIT" << std::endl;
345 static double resPSpitch = 0.19/
sqrt(12.);
346 static double resPSlength = 6.1/
sqrt(12.);
355 std::cout<<
"entering test link by rechit function for ECAL and PS"<<std::endl;
359 double zECAL = clusterECAL.
position().Z();
360 double xECAL = clusterECAL.
position().X();
361 double yECAL = clusterECAL.
position().Y();
364 double zPS = clusterPS.
position().Z();
365 double xPS = clusterPS.
position().X();
366 double yPS = clusterPS.
position().Y();
368 if (zECAL*zPS <0.)
return -1.;
372 switch (clusterPS.
layer()) {
375 deltaX = resPSpitch * sqr12;
376 deltaY = resPSlength * sqr12;
380 deltaY = resPSpitch * sqr12;
381 deltaX = resPSlength * sqr12;
388 const std::vector< reco::PFRecHitFraction >& fracs = clusterECAL.
recHitFractions();
389 bool linkedbyrechit =
false;
391 for(
unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
394 double fraction = fracs[rhit].fraction();
395 if(fraction < 1E-4)
continue;
402 const std::vector< math::XYZPoint >& corners = rechit_cluster.
getCornersXYZ();
403 assert(corners.size() == 4);
408 std::cout <<
"Ecal rechit " << posxyz.X() <<
" " << posxyz.Y() << std::endl;
409 std::cout <<
"PS cluster " << xPS <<
" " << yPS << std::endl;
415 for (
unsigned jc=0; jc<4; ++jc ) {
419 x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
420 y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
425 <<
" " << cornerpos.X() <<
" " << x[jc]
426 <<
" " << cornerpos.Y() <<
" " << y[jc]
439 bool isinside = TMath::IsInside(xPS,yPS,5,x,y);
442 linkedbyrechit =
true;
448 if( linkedbyrechit ) {
449 if( debug )
std::cout <<
"Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
450 double dist =
computeDist( xECAL/1000.,yECAL/1000.,
451 xPS/1000. ,yPS/1000);
469 double dX = posxyzEM.X()-posxyzHAD.X();
470 double dY = posxyzEM.Y()-posxyzHAD.Y();
471 double sameZ = posxyzEM.Z()*posxyzHAD.Z();
473 if(sameZ<0)
return -1.;
475 double dist2 = dX*dX + dY*dY;
479 double dist =
sqrt( dist2 );
489 double eta2,
double phi2 ) {
497 double dist =
std::sqrt( (eta1 - eta2)*(eta1 - eta2)
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
ROOT::Math::PositionVector3D< ROOT::Math::CylindricalEta3D< Double32_t > > REPPoint
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)
static double computeDist(double eta1, double phi1, double eta2, double phi2)
computes a chisquare
const std::vector< math::XYZPoint > & getCornersXYZ() const
rechit corners
bool isNull() const
Checks for null.
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Point of closest approach from beam axis (initial point in the case of PFSimParticle) ...
const REPPoint & positionREP() const
cluster position: rho, eta, phi
const reco::PFTrajectoryPoint & extrapolatedPoint(unsigned layerid) const
const math::XYZPoint & position() const
is seed ?
const REPPointVector & getCornersREP() const
rechit corners
const math::XYZTLorentzVector & momentum() const
4-momenta quadrivector
bool isValid() const
is this point valid ?
XYZPointD XYZPoint
point in space with cartesian internal representation
static double testHFEMAndHFHADByRecHit(const reco::PFCluster &clusterHFEM, const reco::PFCluster &clusterHFHAD, bool debug=false)
test association between HFEM and HFHAD, by rechit
double energy() const
rechit energy
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
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)
const REPPoint & positionREP() const
rechit cell centre rho, eta, phi. call calculatePositionREP before !