Perform the real analysis.
123 map<DTChamberId, PSimHitContainer> simHitsPerCh;
124 for (
const auto &
simHit : *simHits) {
131 simHitsPerCh[chamberId].push_back(
simHit);
142 <<
" in this event, skipping!" << endl;
148 for (
auto &simHitsInChamber : simHitsPerCh) {
151 if (station == 4 && !(
local_)) {
163 int nMuSimHit = muSimHitPerWire.size();
168 cout <<
"=== Chamber " << chamberId <<
" has " << nMuSimHit <<
" SimHits" << endl;
175 if ((
DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() ==
176 (
DTWireId(inAndOutSimHit.second->detUnitId())).superLayer()) {
181 pair<LocalVector, LocalPoint> dirAndPosSimSegm =
184 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
185 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
194 float xSimSeg = simSegmLocalPos.
x();
195 float ySimSeg = simSegmLocalPos.
y();
197 float etaSimSeg = simSegmGlobalPos.
eta();
198 float phiSimSeg = simSegmGlobalPos.
phi();
200 double count_seg = 0;
203 cout <<
" Simulated segment: local direction " << simSegmLocalDir << endl
204 <<
" local position " << simSegmLocalPos << endl
205 <<
" alpha " << alphaSimSeg << endl
206 <<
" beta " << betaSimSeg << endl;
211 bool recHitFound =
false;
213 int nsegm =
distance(range.first, range.second);
215 cout <<
" Chamber: " << chamberId <<
" has " << nsegm <<
" 4D segments" << endl;
225 double deltaAlpha = 99999;
226 double deltaBeta = 99999;
231 if (station != 4 && (*segment4D).dimension() != 4) {
235 LocalVector recSegDirection = (*segment4D).localDirection();
236 LocalPoint recSegPosition = (*segment4D).localPosition();
239 float recSegAlpha = ab.first;
240 float recSegBeta = ab.second;
243 cout << &(*segment4D) <<
" RecSegment direction: " << recSegDirection << endl
244 <<
" position : " << (*segment4D).localPosition() << endl
245 <<
" alpha : " << recSegAlpha << endl
246 <<
" beta : " << recSegBeta << endl
247 <<
" nhits : " << (*segment4D).phiSegment()->recHits().size() <<
" " 248 << (((*segment4D).zSegment() !=
nullptr) ? (*segment4D).zSegment()->recHits().size() : 0) << endl;
251 float dAlphaRecSim = fabs(recSegAlpha - alphaSimSeg);
252 float dBetaRecSim = fabs(recSegBeta - betaSimSeg);
254 if ((fabs(recSegPosition.
x() - simSegmLocalPos.
x()) <
256 && ((fabs(recSegPosition.
y() - simSegmLocalPos.
y()) < 4) ||
257 (*segment4D).dimension() < 4)) {
260 if (fabs(dAlphaRecSim - deltaAlpha) <
epsilon) {
261 if (dBetaRecSim < deltaBeta) {
262 deltaAlpha = dAlphaRecSim;
263 deltaBeta = dBetaRecSim;
264 bestRecHit = &(*segment4D);
267 }
else if (dAlphaRecSim < deltaAlpha) {
268 deltaAlpha = dAlphaRecSim;
269 deltaBeta = dBetaRecSim;
270 bestRecHit = &(*segment4D);
278 cout << endl <<
"Chosen: " << bestRecHit << endl;
288 float alphaBestRHit = ab.first;
289 float betaBestRHit = ab.second;
303 float alphaBestRHitRZ = 0;
304 float alphaSimSegRZ = betaSimSeg;
320 simSegLocalPosRZTmp + simSegLocalDirRZ * (-simSegLocalPosRZTmp.
z() / (
cos(simSegLocalDirRZ.
theta())));
324 cout <<
"RZ SL: recPos " << bestRecHitLocalPosRZ <<
"recDir " << bestRecHitLocalDirRZ <<
"recAlpha " 325 << alphaBestRHitRZ << endl
326 <<
"RZ SL: simPos " << simSegLocalPosRZ <<
"simDir " << simSegLocalDirRZ <<
"simAlpha " 327 << alphaSimSegRZ << endl;
335 float t0theta = -999;
340 t0phi = phiSeg->
t0();
341 nHitPhi = phiSeg->
recHits().size();
345 t0theta = zedRecSeg->
t0();
346 nHitTheta = zedRecSeg->
recHits().size();
353 if (mirrorMinusWheels && wheel < 0) {
355 alphaBestRHit *= -1.;
366 }
else if (
abs(wheel) == 1) {
368 }
else if (
abs(wheel) == 2) {
373 float sigmaBetaBestRhit =
375 bestRecHitLocalDirErr.
yy()));
379 histo->
fill(alphaSimSeg,
384 bestRecHitLocalPos.
x(),
386 bestRecHitLocalPos.
y(),
389 bestRecHitLocalPosRZ.
x(),
390 simSegLocalPosRZ.
x(),
395 sqrt(bestRecHitLocalPosErr.
xx()),
396 sqrt(bestRecHitLocalPosErr.
yy()),
397 sigmaAlphaBestRhitRZ,
398 sqrt(bestRecHitLocalPosErrRZ.
xx()),
409 bestRecHitLocalPos.
x(),
411 bestRecHitLocalPos.
y(),
414 bestRecHitLocalPosRZ.
x(),
415 simSegLocalPosRZ.
x(),
420 sqrt(bestRecHitLocalPosErr.
xx()),
421 sqrt(bestRecHitLocalPosErr.
yy()),
422 sigmaAlphaBestRhitRZ,
423 sqrt(bestRecHitLocalPosErrRZ.
xx()),
430 histograms.h4DHitWS[
abs(wheel)][station - 1]->fill(alphaSimSeg,
435 bestRecHitLocalPos.
x(),
437 bestRecHitLocalPos.
y(),
440 bestRecHitLocalPosRZ.
x(),
441 simSegLocalPosRZ.
x(),
446 sqrt(bestRecHitLocalPosErr.
xx()),
447 sqrt(bestRecHitLocalPosErr.
yy()),
448 sigmaAlphaBestRhitRZ,
449 sqrt(bestRecHitLocalPosErrRZ.
xx()),
466 }
else if (
abs(wheel) == 1) {
468 }
else if (
abs(wheel) == 2) {
471 heff->
fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
473 etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound, count_seg);
476 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
std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
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.
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
C::const_iterator const_iterator
constant access iterator type
void fill(float etaSimSegm, float phiSimSegm, float xSimSegm, float ySimSegm, float alphaSimSegm, float betaSimSegm, bool fillRecHit, int nSeg)
std::map< DTWireId, const PSimHit * > mapMuSimHitsPerWire(const std::map< DTWireId, edm::PSimHitContainer > &simHitWireMap)
Create a map between the Mu SimHits and corresponding MuBarWireId ;.
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)
std::pair< const PSimHit *, const PSimHit * > findMuSimSegment(const std::map< DTWireId, const PSimHit * > &mapWireAndMuSimHit)
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)
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.