Perform the real analysis.
161 map<DTChamberId, PSimHitContainer > simHitsPerCh;
162 for(PSimHitContainer::const_iterator simHit = simHits->begin();
163 simHit != simHits->end(); simHit++){
166 if (
abs((*simHit).particleType())==13) {
168 DTChamberId chamberId = (((
DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId();
170 simHitsPerCh[chamberId].push_back(*simHit);
180 <<
" in this event, skipping!" << endl;
185 for (map<DTChamberId, PSimHitContainer>::const_iterator simHitsInChamber = simHitsPerCh.begin(); simHitsInChamber != simHitsPerCh.end(); ++simHitsInChamber) {
189 if (station == 4 && !(
local) )
continue;
190 int wheel = chamberId.
wheel();
199 int nMuSimHit = muSimHitPerWire.size();
204 cout <<
"=== Chamber " << chamberId <<
" has " << nMuSimHit <<
" SimHits" << endl;
210 if ((
DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() == (
DTWireId(inAndOutSimHit.second->detUnitId())).superLayer())
continue;
214 chamberId,&(*dtGeom));
216 LocalVector simSegmLocalDir = dirAndPosSimSegm.first;
217 LocalPoint simSegmLocalPos = dirAndPosSimSegm.second;
218 const DTChamber* chamber = dtGeom->chamber(chamberId);
226 float xSimSeg = simSegmLocalPos.
x();
227 float ySimSeg = simSegmLocalPos.
y();
229 float etaSimSeg = simSegmGlobalPos.
eta();
230 float phiSimSeg = simSegmGlobalPos.
phi();
232 double count_seg = 0;
235 cout<<
" Simulated segment: local direction "<<simSegmLocalDir<<endl
236 <<
" local position "<<simSegmLocalPos<<endl
237 <<
" alpha "<<alphaSimSeg<<endl
238 <<
" beta "<<betaSimSeg<<endl;
242 bool recHitFound =
false;
244 int nsegm =
distance(range.first, range.second);
246 cout <<
" Chamber: " << chamberId <<
" has " << nsegm
247 <<
" 4D segments" << endl;
256 double deltaAlpha = 99999;
257 double deltaBeta = 99999;
261 segment4D!=range.second;
265 if(station!=4 && (*segment4D).dimension() != 4) {
269 LocalVector recSegDirection = (*segment4D).localDirection();
270 LocalPoint recSegPosition = (*segment4D).localPosition();
273 float recSegAlpha = ab.first;
274 float recSegBeta = ab.second;
277 cout << &(*segment4D)
278 <<
" RecSegment direction: " << recSegDirection << endl
279 <<
" position : " << (*segment4D).localPosition() << endl
280 <<
" alpha : " << recSegAlpha << endl
281 <<
" beta : " << recSegBeta << endl
282 <<
" nhits : " << (*segment4D).phiSegment()->recHits().size() <<
" " << (((*segment4D).zSegment()!=0)?(*segment4D).zSegment()->recHits().size():0)
286 float dAlphaRecSim=fabs(recSegAlpha - alphaSimSeg);
287 float dBetaRecSim =fabs(recSegBeta - betaSimSeg);
290 if ((fabs(recSegPosition.
x()-simSegmLocalPos.
x())<4)
291 && ((fabs(recSegPosition.
y()-simSegmLocalPos.
y())<4)||(*segment4D).dimension()<4)) {
294 if (fabs(dAlphaRecSim-deltaAlpha) <
epsilon) {
295 if (dBetaRecSim < deltaBeta) {
296 deltaAlpha = dAlphaRecSim;
297 deltaBeta = dBetaRecSim;
298 bestRecHit = &(*segment4D);
301 }
else if (dAlphaRecSim < deltaAlpha) {
302 deltaAlpha = dAlphaRecSim;
303 deltaBeta = dBetaRecSim;
304 bestRecHit = &(*segment4D);
311 if (
debug)
cout << endl <<
"Chosen: " << bestRecHit << endl;
320 float alphaBestRHit = ab.first;
321 float betaBestRHit = ab.second;
334 float alphaBestRHitRZ = 0;
335 float alphaSimSegRZ = betaSimSeg;
351 simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.
z()/(
cos(simSegLocalDirRZ.
theta())));
355 "RZ SL: recPos " << bestRecHitLocalPosRZ <<
356 "recDir " << bestRecHitLocalDirRZ <<
357 "recAlpha " << alphaBestRHitRZ << endl <<
358 "RZ SL: simPos " << simSegLocalPosRZ <<
359 "simDir " << simSegLocalDirRZ <<
360 "simAlpha " << alphaSimSegRZ << endl ;
367 float t0theta = -999;
372 t0phi = phiSeg->
t0();
373 nHitPhi = phiSeg->
recHits().size();
377 t0theta = zedRecSeg->
t0();
378 nHitTheta = zedRecSeg->
recHits().size();
391 if (mirrorMinusWheels && wheel<0){
403 else if(
abs(wheel) == 1)
405 else if(
abs(wheel) == 2)
412 histo->
Fill(alphaSimSeg,
417 bestRecHitLocalPos.
x(),
419 bestRecHitLocalPos.
y(),
422 bestRecHitLocalPosRZ.
x(),
423 simSegLocalPosRZ.
x(),
428 sqrt(bestRecHitLocalPosErr.
xx()),
429 sqrt(bestRecHitLocalPosErr.
yy()),
430 sigmaAlphaBestRhitRZ,
431 sqrt(bestRecHitLocalPosErrRZ.
xx()),
432 nHitPhi,nHitTheta,t0phi,t0theta
440 bestRecHitLocalPos.
x(),
442 bestRecHitLocalPos.
y(),
445 bestRecHitLocalPosRZ.
x(),
446 simSegLocalPosRZ.
x(),
451 sqrt(bestRecHitLocalPosErr.
xx()),
452 sqrt(bestRecHitLocalPosErr.
yy()),
453 sigmaAlphaBestRhitRZ,
454 sqrt(bestRecHitLocalPosErrRZ.
xx()),
455 nHitPhi,nHitTheta,t0phi,t0theta
463 bestRecHitLocalPos.
x(),
465 bestRecHitLocalPos.
y(),
468 bestRecHitLocalPosRZ.
x(),
469 simSegLocalPosRZ.
x(),
474 sqrt(bestRecHitLocalPosErr.
xx()),
475 sqrt(bestRecHitLocalPosErr.
yy()),
476 sigmaAlphaBestRhitRZ,
477 sqrt(bestRecHitLocalPosErrRZ.
xx()),
478 nHitPhi,nHitTheta,t0phi,t0theta
491 else if(
abs(wheel) == 1)
493 else if(
abs(wheel) == 2)
495 heff->
Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg);
496 hEff_All->
Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg);
497 if (
local)
hEffWS[
abs(wheel)][station-1]->
Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg);
virtual LocalError localPositionError() const
local position error in SL frame
void Fill(float etaSimSegm, float phiSimSegm, float xSimSegm, float ySimSegm, float alphaSimSegm, float betaSimSegm, bool fillRecHit, int nSeg)
std::pair< const_iterator, const_iterator > range
iterator range
static double sigmaAngle(double Angle, double sigma2TanAngle)
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
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
virtual LocalError localDirectionError() const
the local direction error (xx,xy,yy) in SL frame: only xx is not 0.
edm::EDGetTokenT< DTRecSegment4DCollection > segment4DToken_
static std::map< DTWireId, const PSimHit * > mapMuSimHitsPerWire(const std::map< DTWireId, edm::PSimHitContainer > &simHitWireMap)
Create a map between the Mu SimHits and corresponding MuBarWireId ;.
Geom::Theta< T > theta() const
virtual LocalError localPositionError() const
Local position error in Chamber frame.
static 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) ...
virtual LocalVector localDirection() const
Local direction in Chamber frame.
edm::InputTag segment4DLabel
HRes4DHit * h4DHitWS[3][4]
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)
virtual std::vector< const TrackingRecHit * > recHits() const
Access to component RecHits (if any)
Cos< T >::type cos(const T &t)
A set of histograms for efficiency 4D RecHits.
Abs< T >::type abs(const T &t)
virtual LocalPoint localPosition() const
Local position in Chamber frame.
static std::pair< double, double > findSegmentAlphaAndBeta(const LocalVector &direction)
DTSuperLayerId superLayerId() const
The id of the superlayer on which reside the segment.
virtual LocalPoint localPosition() const
local position in SL frame
const DTSLRecSegment2D * zSegment() const
The Z segment: 0 if not zed projection available.
static 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...
virtual LocalVector localDirection() const
the local direction in SL frame
std::vector< PSimHit > PSimHitContainer
static std::map< DTWireId, edm::PSimHitContainer > mapSimHitsPerWire(const edm::PSimHitContainer &simhits)
double t0() const
Get the segment t0 (if recomputed, 0 is returned otherwise)
int station() const
Return the station number.
edm::EDGetTokenT< edm::PSimHitContainer > simHitToken_
int wheel() const
Return the wheel number.
virtual LocalError localDirectionError() const
Local direction error in the Chamber frame.