Perform the real analysis.
119 map<DTChamberId, PSimHitContainer > simHitsPerCh;
120 for (
const auto &
simHit : *simHits) {
127 simHitsPerCh[chamberId].push_back(
simHit);
137 cout <<
"[DTSegment4DQuality]**Warning: no 4D Segments with label: " <<
segment4DLabel_ <<
" in this event, skipping!" << endl;
143 for (
auto & simHitsInChamber : simHitsPerCh) {
147 if (station == 4 && !(
local_) ) {
continue;
158 int nMuSimHit = muSimHitPerWire.size();
163 cout <<
"=== Chamber " << chamberId <<
" has " << nMuSimHit <<
" SimHits" << endl;
170 if ((
DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() == (
DTWireId(inAndOutSimHit.second->detUnitId())).superLayer()) {
continue;
175 chamberId,&(*dtGeom));
177 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
178 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
187 float xSimSeg = simSegmLocalPos.
x();
188 float ySimSeg = simSegmLocalPos.
y();
190 float etaSimSeg = simSegmGlobalPos.
eta();
191 float phiSimSeg = simSegmGlobalPos.
phi();
193 double count_seg = 0;
196 cout <<
" Simulated segment: local direction " << simSegmLocalDir << endl
197 <<
" local position " << simSegmLocalPos << endl
198 <<
" alpha " << alphaSimSeg << endl
199 <<
" beta " << betaSimSeg << endl;
204 bool recHitFound =
false;
206 int nsegm =
distance(range.first, range.second);
208 cout <<
" Chamber: " << chamberId <<
" has " << nsegm
209 <<
" 4D segments" << endl;
219 double deltaAlpha = 99999;
220 double deltaBeta = 99999;
224 segment4D!= range.second;
228 if (station!= 4 && (*segment4D).dimension() != 4) {
232 LocalVector recSegDirection = (*segment4D).localDirection();
233 LocalPoint recSegPosition = (*segment4D).localPosition();
236 float recSegAlpha = ab.first;
237 float recSegBeta = ab.second;
240 cout << &(*segment4D)
241 <<
" RecSegment direction: " << recSegDirection << endl
242 <<
" position : " << (*segment4D).localPosition() << endl
243 <<
" alpha : " << recSegAlpha << endl
244 <<
" beta : " << recSegBeta << endl
245 <<
" nhits : " << (*segment4D).phiSegment()->recHits().size() <<
" " << (((*segment4D).zSegment()!=
nullptr)?(*segment4D).zSegment()->recHits().size():0)
250 float dAlphaRecSim = fabs(recSegAlpha - alphaSimSeg);
251 float dBetaRecSim = fabs(recSegBeta - betaSimSeg);
254 if ((fabs(recSegPosition.
x()-simSegmLocalPos.
x())<4)
255 && ((fabs(recSegPosition.
y()-simSegmLocalPos.
y())<4)||(*segment4D).dimension()<4)) {
258 if (fabs(dAlphaRecSim-deltaAlpha) <
epsilon) {
259 if (dBetaRecSim < deltaBeta) {
260 deltaAlpha = dAlphaRecSim;
261 deltaBeta = dBetaRecSim;
262 bestRecHit = &(*segment4D);
265 }
else if (dAlphaRecSim < deltaAlpha) {
266 deltaAlpha = dAlphaRecSim;
267 deltaBeta = dBetaRecSim;
268 bestRecHit = &(*segment4D);
276 cout << endl <<
"Chosen: " << bestRecHit << endl;
286 float alphaBestRHit = ab.first;
287 float betaBestRHit = ab.second;
300 float alphaBestRHitRZ = 0;
301 float alphaSimSegRZ = betaSimSeg;
317 simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.
z()/(
cos(simSegLocalDirRZ.
theta())));
322 "RZ SL: recPos " << bestRecHitLocalPosRZ <<
323 "recDir " << bestRecHitLocalDirRZ <<
324 "recAlpha " << alphaBestRHitRZ << endl <<
325 "RZ SL: simPos " << simSegLocalPosRZ <<
326 "simDir " << simSegLocalDirRZ <<
327 "simAlpha " << alphaSimSegRZ << endl ;
335 float t0theta = -999;
340 t0phi = phiSeg->
t0();
341 nHitPhi = phiSeg->
recHits().size();
345 t0theta = zedRecSeg->
t0();
346 nHitTheta = zedRecSeg->
recHits().size();
352 if (mirrorMinusWheels && wheel<0) {
364 }
else if (
abs(wheel) == 1) {
366 }
else if (
abs(wheel) == 2) {
374 histo->
fill(alphaSimSeg,
379 bestRecHitLocalPos.
x(),
381 bestRecHitLocalPos.
y(),
384 bestRecHitLocalPosRZ.
x(),
385 simSegLocalPosRZ.
x(),
390 sqrt(bestRecHitLocalPosErr.
xx()),
391 sqrt(bestRecHitLocalPosErr.
yy()),
392 sigmaAlphaBestRhitRZ,
393 sqrt(bestRecHitLocalPosErrRZ.
xx()),
394 nHitPhi, nHitTheta, t0phi, t0theta);
401 bestRecHitLocalPos.
x(),
403 bestRecHitLocalPos.
y(),
406 bestRecHitLocalPosRZ.
x(),
407 simSegLocalPosRZ.
x(),
412 sqrt(bestRecHitLocalPosErr.
xx()),
413 sqrt(bestRecHitLocalPosErr.
yy()),
414 sigmaAlphaBestRhitRZ,
415 sqrt(bestRecHitLocalPosErrRZ.
xx()),
416 nHitPhi, nHitTheta, t0phi, t0theta);
423 bestRecHitLocalPos.
x(),
425 bestRecHitLocalPos.
y(),
428 bestRecHitLocalPosRZ.
x(),
429 simSegLocalPosRZ.
x(),
434 sqrt(bestRecHitLocalPosErr.
xx()),
435 sqrt(bestRecHitLocalPosErr.
yy()),
436 sigmaAlphaBestRhitRZ,
437 sqrt(bestRecHitLocalPosErrRZ.
xx()),
438 nHitPhi, nHitTheta, t0phi, t0theta);
451 }
else if (
abs(wheel) == 1) {
453 }
else if (
abs(wheel) == 2) {
456 heff->
fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
457 histograms.hEff_All->fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
459 histograms.hEffWS[
abs(wheel)][station-1]->fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
std::pair< double, double > findSegmentAlphaAndBeta(const LocalVector &direction)
LocalError localPositionError() const override
local position error in SL frame
LocalPoint localPosition() const override
Local position in Chamber frame.
std::pair< const_iterator, const_iterator > range
iterator range
const DTChamber * chamber(const DTChamberId &id) const
Return a DTChamber given its id.
std::map< DTWireId, const PSimHit * > mapMuSimHitsPerWire(const std::map< DTWireId, edm::PSimHitContainer > &simHitWireMap)
Create a map between the Mu SimHits and corresponding MuBarWireId ;.
LocalVector localDirection() const override
Local direction in Chamber frame.
const DTChamberRecSegment2D * phiSegment() const
The superPhi segment: 0 if no phi projection available.
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Geom::Phi< T > phi() const
def setup(process, global_tag, zero_tesla=False)
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
edm::InputTag segment4DLabel_
edm::EDGetTokenT< DTRecSegment4DCollection > segment4DToken_
Geom::Theta< T > theta() const
std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
void fill(float etaSimSegm, float phiSimSegm, float xSimSegm, float ySimSegm, float alphaSimSegm, float betaSimSegm, bool fillRecHit, int nSeg)
Cos< T >::type cos(const T &t)
LocalError localPositionError() const override
Local position error in Chamber frame.
A set of histograms for efficiency 4D RecHits (producer)
Abs< T >::type abs(const T &t)
LocalError localDirectionError() const override
Local direction error in the Chamber frame.
DTSuperLayerId superLayerId() const
The id of the superlayer on which reside the segment.
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
LocalPoint localPosition() const override
local position in SL frame
void fill(float simDirectionAlpha, float recDirectionAlpha, float simDirectionBeta, float recDirectionBeta, float simX, float recX, float simY, float recY, float simEta, float simPhi, float recYRZ, float simYRZ, float recBetaRZ, float simBetaRZ, float sigmaAlpha, float sigmaBeta, float sigmaX, float sigmaY, float sigmaBetaRZ, float sigmaYRZ, int nHitsPhi, int nHitsTheta, float t0Phi, float t0Theta)
std::pair< LocalVector, LocalPoint > findMuSimSegmentDirAndPos(const std::pair< const PSimHit *, const PSimHit * > &inAndOutSimHit, const DetId detId, const DTGeometry *muonGeom)
Find direction and position of a segment (in local RF) from outer and inner mu SimHit in the RF of ob...
std::pair< const PSimHit *, const PSimHit * > findMuSimSegment(const std::map< DTWireId, const PSimHit * > &mapWireAndMuSimHit)
Find Innermost and outermost SimHit from Mu in a SL (they identify a simulated segment) ...
LocalVector localDirection() const override
the local direction in SL frame
std::vector< const TrackingRecHit * > recHits() const override
Access to component RecHits (if any)
std::vector< PSimHit > PSimHitContainer
double sigmaAngle(double Angle, double sigma2TanAngle)
double t0() const
Get the segment t0 (if recomputed, 0 is returned otherwise)
int station() const
Return the station number.
LocalError localDirectionError() const override
the local direction error (xx,xy,yy) in SL frame: only xx is not 0.
edm::EDGetTokenT< edm::PSimHitContainer > simHitToken_
int wheel() const
Return the wheel number.
const DTSuperLayer * superLayer(const DTSuperLayerId &id) const
Return a DTSuperLayer given its id.