32 #include "TDirectory.h"
58 if(theDataType.label() !=
"RealData" && theDataType.label() !=
"SimData")
59 LogDebug(
"MuonTrackResidualAnalyzer")<<
"Error in Data Type!!";
66 theMuonSimHitNumberPerEvent = 0;
78 LogDebug(
"MuonTrackResidualAnalyzer")<<
"Begin Job";
88 dirName+=algo.
label()+
"_";
91 if (dirName.find(
"Tracks")<dirName.length()){
92 dirName.replace(dirName.find(
"Tracks"),6,
"");
98 hDPtRef =
dbe_->
book1D(
"DeltaPtRef",
"P^{in}_{t}-P^{in ref}",10000,-20,20);
106 hSimHitsPerTrack =
dbe_->
book1D(
"SimHitsPerTrack",
"Number of sim hits per track",55,0,55);
107 hSimHitsPerTrackVsEta =
dbe_->
book2D(
"SimHitsPerTrackVsEta",
"Number of sim hits per track VS #eta",120,-3.,3.,55,0,55);
108 hDeltaPtVsEtaSim =
dbe_->
book2D(
"DeltaPtVsEtaSim",
"#Delta P_{t} vs #eta gen, sim quantity",120,-3.,3.,500,-250.,250.);
109 hDeltaPtVsEtaSim2 =
dbe_->
book2D(
"DeltaPtVsEtaSim2",
"#Delta P_{t} vs #eta gen, sim quantity",120,-3.,3.,500,-250.,250.);
118 LogDebug(
"MuonTrackResidualAnalyzer")<<
"Analyze";
121 theService->update(eventSetup);
123 theMuonSimHitNumberPerEvent = 0;
127 event.getByLabel(dtSimHitLabel.instance(),dtSimHitLabel.label(), dtSimHits);
130 event.getByLabel(cscSimHitLabel.instance(),cscSimHitLabel.label(), cscSimHits);
133 event.getByLabel(rpcSimHitLabel.instance(),rpcSimHitLabel.label(), rpcSimHits);
140 map<DetId, const PSimHit* > muonSimHitsPerId =
141 mapMuSimHitsPerId(dtSimHits,cscSimHits,rpcSimHits);
143 hSimHitsPerTrack->Fill(theMuonSimHitNumberPerEvent);
147 if(theDataType.label() ==
"SimData"){
150 event.getByLabel(theDataType.instance(),simTracks);
153 SimTrackContainer::const_iterator simTrack;
155 for (simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack)
156 if (
abs((*simTrack).type()) == 13){
157 hSimHitsPerTrackVsEta->Fill((*simTrack).momentum().eta(),theMuonSimHitNumberPerEvent);
158 etaSim = (*simTrack).momentum().eta();
159 theSimTkId = (*simTrack).trackId();
167 event.getByLabel(theMuonTrackLabel, muonTracks);
169 reco::TrackCollection::const_iterator muonTrack;
172 for (muonTrack = muonTracks->begin(); muonTrack != muonTracks->end(); ++muonTrack) {
185 double momAtEntry = -150., momAtExit = -150.;
187 if(theSimHitContainer.size()>1){
189 const GeomDet *geomDetA = theService->trackingGeometry()->idToDet(
DetId(theSimHitContainer.front()->detUnitId()));
190 double distA = geomDetA->
toGlobal(theSimHitContainer.front()->localPosition()).
mag();
192 const GeomDet *geomDetB = theService->trackingGeometry()->idToDet(
DetId(theSimHitContainer.back()->detUnitId()));
193 double distB = geomDetB->
toGlobal(theSimHitContainer.back()->localPosition()).
mag();
195 LogTrace(
"MuonTrackResidualAnalyzer")<<
"Inner SimHit: "<<theSimHitContainer.front()->particleType()
196 <<
" Pt: "<<theSimHitContainer.front()->momentumAtEntry().perp()
197 <<
" E: "<<theSimHitContainer.front()->momentumAtEntry().perp()
198 <<
" R: "<<distA<<endl;
199 LogTrace(
"MuonTrackResidualAnalyzer")<<
"Outer SimHit: "<<theSimHitContainer.back()->particleType()
200 <<
" Pt: "<<theSimHitContainer.back()->momentumAtEntry().perp()
201 <<
" E: "<<theSimHitContainer.front()->momentumAtEntry().perp()
202 <<
" R: "<<distB<<endl;
204 momAtEntry = theSimHitContainer.front()->momentumAtEntry().perp();
205 momAtExit = theSimHitContainer.back()->momentumAtEntry().perp();
210 map<DetId,const PSimHit*>::const_iterator itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
211 map<DetId,const PSimHit*>::const_iterator itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
213 double momAtEntry2 = -150, momAtExit2 = -150.;
214 if (itFirst != muonSimHitsPerId.end() )
215 momAtEntry2 = itFirst->second->momentumAtEntry().perp();
217 LogDebug(
"MuonTrackResidualAnalyzer")<<
"No first sim hit found";
219 itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
220 if (itFirst != muonSimHitsPerId.end() )
221 momAtEntry2 = itFirst->second->momentumAtEntry().perp();
223 LogDebug(
"MuonTrackResidualAnalyzer")<<
"No second sim hit found";
228 if (itLast != muonSimHitsPerId.end() )
229 momAtExit2 = itLast->second->momentumAtEntry().perp();
231 LogDebug(
"MuonTrackResidualAnalyzer")<<
"No last sim hit found";
233 itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
234 if (itLast != muonSimHitsPerId.end() )
235 momAtExit2 = itLast->second->momentumAtEntry().perp();
237 LogDebug(
"MuonTrackResidualAnalyzer")<<
"No last but one sim hit found";
243 if(momAtEntry >=0 && momAtExit >= 0)
244 hDeltaPtVsEtaSim->Fill(etaSim,momAtEntry-momAtExit);
245 if(momAtEntry2 >=0 && momAtExit2 >= 0)
246 hDeltaPtVsEtaSim2->Fill(etaSim,momAtEntry2-momAtExit2);
249 LogDebug(
"MuonTrackResidualAnalyzer")<<
"NO SimTrack'eta";
261 return (
abs(eta) <= 2.4 ) ?
true :
false;
263 return (
abs(eta) < 1.1 ) ?
true :
false;
265 return (
abs(eta) >= 1.1 &&
abs(eta) <= 2.4 ) ?
true :
false;
267 {
LogDebug(
"MuonTrackResidualAnalyzer")<<
"No correct Eta range selected!! ";
return false;}
272 map<DetId,const PSimHit*>
279 map<DetId,const PSimHit*> hitIdMap;
280 theSimHitContainer.clear();
282 mapMuSimHitsPerId(dtSimhits,hitIdMap);
283 mapMuSimHitsPerId(cscSimhits,hitIdMap);
284 mapMuSimHitsPerId(rpcSimhits,hitIdMap);
286 if(theSimHitContainer.size() >1)
287 stable_sort(theSimHitContainer.begin(),theSimHitContainer.end(),
RadiusComparatorInOut(theService->trackingGeometry()));
289 LogDebug(
"MuonTrackResidualAnalyzer")<<
"Sim Hit list";
291 for(vector<const PSimHit*>::const_iterator it = theSimHitContainer.begin();
292 it != theSimHitContainer.end(); ++it){
293 LogTrace(
"MuonTrackResidualAnalyzer")<<count
295 <<
" Process Type: " << (*it)->processType()
305 map<DetId,const PSimHit*> &hitIdMap){
307 for(PSimHitContainer::const_iterator simhit = simhits->begin();
308 simhit != simhits->end(); ++simhit) {
310 if (
abs(simhit->particleType()) != 13 && theSimTkId != simhit->trackId())
continue;
312 theSimHitContainer.push_back(&*simhit);
320 map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(
id);
322 if (it == hitIdMap.end() )
323 hitIdMap[
id] = &*simhit;
325 LogDebug(
"MuonTrackResidualAnalyzer")<<
"TWO muons in the same sensible volume!!";
327 ++theMuonSimHitNumberPerEvent;
333 map<DetId,const PSimHit*> &hitIdMap,
338 for(Trajectory::DataContainer::const_iterator datum = data.begin();
339 datum != data.end(); ++datum){
341 GlobalPoint fitPoint = datum->updatedState().globalPosition();
348 double errX = datum->updatedState().localError().matrix_old()[3][3];
349 double errY = datum->updatedState().localError().matrix_old()[4][4];
352 map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(datum->recHit()->geographicalId());
355 if (it == hitIdMap.end() )
continue;
357 const PSimHit* simhit = it->second;
368 cout <<
"SimHit position "<< simHitPoint << endl;
369 cout <<
"Fit position "<< fitLocalPoint << endl;
370 cout <<
"Fit position2 "<< datum->updatedState().localPosition() << endl;
371 cout <<
"Errors on the fit position: (" << errX <<
"," << errY <<
"," << errZ <<
")"<<endl;
372 cout <<
"Resolution on x: " << diff.
x()/
abs(simHitPoint.x()) << endl;
373 cout <<
"Resolution on y: " << diff.
y()/
abs(simHitPoint.y()) << endl;
374 cout <<
"Resolution on z: " << diff.
z()/
abs(simHitPoint.z()) << endl;
376 cout <<
"Eta direction: "<< simhit->
momentumAtEntry().
eta() <<
" eta position: " << simHitPoint.eta() << endl;
377 cout <<
"Phi direction: "<< simhit->
momentumAtEntry().
phi() <<
" phi position: " << simHitPoint.phi() << endl;
380 histos->
Fill( simHitPoint.x(), simHitPoint.y(), simHitPoint.z(),
381 diff.
x(), diff.
y(), diff.
z(),
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
std::map< DetId, const PSimHit * > mapMuSimHitsPerId(edm::Handle< edm::PSimHitContainer > dtSimhits, edm::Handle< edm::PSimHitContainer > cscSimhits, edm::Handle< edm::PSimHitContainer > rpcSimhits)
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
bool isInTheAcceptance(double eta)
void cd(void)
go to top directory (ie. root)
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Geom::Phi< T > phi() const
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
MuonTrackResidualAnalyzer(const edm::ParameterSet &pset)
Constructor.
std::string dumpMuonId(const DetId &id) const
uint32_t rawId() const
get the raw id
DataContainer const & measurements() const
std::vector< TrajectoryMeasurement > DataContainer
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Local3DPoint localPosition() const
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
tuple Chi2MeasurementEstimator
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
std::vector< ConstRecHitPointer > ConstRecHitContainer
virtual ~MuonTrackResidualAnalyzer()
Destructor.
void Fill(double x, double y, double z, double dx, double dy, double dz, double errx, double erry, double errz, double eta, double phi)
void showDirStructure(void) const
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
void computeResolution(Trajectory &trajectory, std::map< DetId, const PSimHit * > &hitIdMap, HResolution1DRecHit *histos)
void setCurrentFolder(const std::string &fullpath)
unsigned int detUnitId() const