CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonTrackResidualAnalyzer.cc
Go to the documentation of this file.
2 
3 // Collaborating Class Header
8 
10 
16 
21 
27 
28 
29 #include "TFile.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TDirectory.h"
33 
34 using namespace std;
35 using namespace edm;
36 
39 
40  // service parameters
41  ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
42  // the services
43  theService = new MuonServiceProxy(serviceParameters);
44 
45  theMuonTrackLabel = pset.getParameter<InputTag>("MuonTrack");
46  theSeedCollectionLabel = pset.getParameter<InputTag>("MuonSeed");
47 
48  cscSimHitLabel = pset.getParameter<InputTag>("CSCSimHit");
49  dtSimHitLabel = pset.getParameter<InputTag>("DTSimHit");
50  rpcSimHitLabel = pset.getParameter<InputTag>("RPCSimHit");
51 
53  out = pset.getUntrackedParameter<string>("rootFileName");
54  dirName_ = pset.getUntrackedParameter<std::string>("dirName");
55 
56  // Sim or Real
57  theDataType = pset.getParameter<InputTag>("DataType");
58  if(theDataType.label() != "RealData" && theDataType.label() != "SimData")
59  LogDebug("MuonTrackResidualAnalyzer")<<"Error in Data Type!!";
60 
61  theEtaRange = (EtaRange) pset.getParameter<int>("EtaRange");
62 
63  theUpdator = new KFUpdator();
64  theEstimator = new Chi2MeasurementEstimator(100000.);
65 
66  theMuonSimHitNumberPerEvent = 0;
67 }
68 
71  delete theUpdator;
72  delete theEstimator;
73  delete theService;
74 }
75 
76 // Operations
78  LogDebug("MuonTrackResidualAnalyzer")<<"Begin Job";
79 
81 
82  dbe_->cd();
83  InputTag algo = theMuonTrackLabel;
84  string dirName=dirName_;
85  if (algo.process()!="")
86  dirName+=algo.process()+"_";
87  if(algo.label()!="")
88  dirName+=algo.label()+"_";
89  if(algo.instance()!="")
90  dirName+=algo.instance()+"";
91  if (dirName.find("Tracks")<dirName.length()){
92  dirName.replace(dirName.find("Tracks"),6,"");
93  }
94  std::replace(dirName.begin(), dirName.end(), ':', '_');
95  dbe_->setCurrentFolder(dirName.c_str());
96 
97 
98  hDPtRef = dbe_->book1D("DeltaPtRef","P^{in}_{t}-P^{in ref}",10000,-20,20);
99 
100  // Resolution wrt the 1D Rec Hits
101  // h1DRecHitRes = new HResolution1DRecHit("TotalRec");
102 
103  // Resolution wrt the 1d Sim Hits
104  // h1DSimHitRes = new HResolution1DRecHit("TotalSim");
105 
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.);
110 }
111 
113  if ( out.size() != 0 && dbe_ ) dbe_->save(out);
114 }
115 
116 
118  LogDebug("MuonTrackResidualAnalyzer")<<"Analyze";
119 
120  // Update the services
121  theService->update(eventSetup);
123  theMuonSimHitNumberPerEvent = 0;
124 
125  // Get the SimHit collection from the event
126  Handle<PSimHitContainer> dtSimHits;
127  event.getByLabel(dtSimHitLabel.instance(),dtSimHitLabel.label(), dtSimHits);
128 
129  Handle<PSimHitContainer> cscSimHits;
130  event.getByLabel(cscSimHitLabel.instance(),cscSimHitLabel.label(), cscSimHits);
131 
132  Handle<PSimHitContainer> rpcSimHits;
133  event.getByLabel(rpcSimHitLabel.instance(),rpcSimHitLabel.label(), rpcSimHits);
134 
135  Handle<SimTrackContainer> simTracks;
136 
137  // FIXME Add the tracker one??
138 
139  // Map simhits per DetId
140  map<DetId, const PSimHit* > muonSimHitsPerId =
141  mapMuSimHitsPerId(dtSimHits,cscSimHits,rpcSimHits);
142 
143  hSimHitsPerTrack->Fill(theMuonSimHitNumberPerEvent);
144 
145  double etaSim=0;
146 
147  if(theDataType.label() == "SimData"){
148 
149  // Get the SimTrack collection from the event
150  event.getByLabel(theDataType.instance(),simTracks);
151 
152  // Loop over the Sim tracks
153  SimTrackContainer::const_iterator simTrack;
154 
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();
160 
161  }
162  }
163 
164 
165  // Get the RecTrack collection from the event
167  event.getByLabel(theMuonTrackLabel, muonTracks);
168 
169  reco::TrackCollection::const_iterator muonTrack;
170 
171  // Loop over the Rec tracks
172  for (muonTrack = muonTracks->begin(); muonTrack != muonTracks->end(); ++muonTrack) {
173 
174  reco::TransientTrack track(*muonTrack,&*theService->magneticField(),theService->trackingGeometry());
175 
176  TrajectoryStateOnSurface outerTSOS = track.outermostMeasurementState();
177  TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
178 
180 
181  trackingRecHit_iterator rhend = track.recHitsBegin()-1;
182  trackingRecHit_iterator rhbegin = track.recHitsEnd()-2;
183 
184  // SimHit Energy loss analysis
185  double momAtEntry = -150., momAtExit = -150.;
186 
187  if(theSimHitContainer.size()>1){
188 
189  const GeomDet *geomDetA = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.front()->detUnitId()));
190  double distA = geomDetA->toGlobal(theSimHitContainer.front()->localPosition()).mag();
191 
192  const GeomDet *geomDetB = theService->trackingGeometry()->idToDet(DetId(theSimHitContainer.back()->detUnitId()));
193  double distB = geomDetB->toGlobal(theSimHitContainer.back()->localPosition()).mag();
194 
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;
203 
204  momAtEntry = theSimHitContainer.front()->momentumAtEntry().perp();
205  momAtExit = theSimHitContainer.back()->momentumAtEntry().perp();
206  }
207 
208  trackingRecHit_iterator rhFirst = track.recHitsBegin();
209  trackingRecHit_iterator rhLast = track.recHitsEnd()-1;
210  map<DetId,const PSimHit*>::const_iterator itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
211  map<DetId,const PSimHit*>::const_iterator itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
212 
213  double momAtEntry2 = -150, momAtExit2 = -150.;
214  if (itFirst != muonSimHitsPerId.end() )
215  momAtEntry2 = itFirst->second->momentumAtEntry().perp();
216  else {
217  LogDebug("MuonTrackResidualAnalyzer")<<"No first sim hit found";
218  ++rhFirst;
219  itFirst = muonSimHitsPerId.find((*rhFirst)->geographicalId());
220  if (itFirst != muonSimHitsPerId.end() )
221  momAtEntry2 = itFirst->second->momentumAtEntry().perp();
222  else{
223  LogDebug("MuonTrackResidualAnalyzer")<<"No second sim hit found";
224  // continue;
225  }
226  }
227 
228  if (itLast != muonSimHitsPerId.end() )
229  momAtExit2 = itLast->second->momentumAtEntry().perp();
230  else {
231  LogDebug("MuonTrackResidualAnalyzer")<<"No last sim hit found";
232  --rhLast;
233  itLast = muonSimHitsPerId.find((*rhLast)->geographicalId());
234  if (itLast != muonSimHitsPerId.end() )
235  momAtExit2 = itLast->second->momentumAtEntry().perp();
236  else{
237  LogDebug("MuonTrackResidualAnalyzer")<<"No last but one sim hit found";
238  // continue;
239  }
240  }
241 
242  if(etaSim){
243  if(momAtEntry >=0 && momAtExit >= 0)
244  hDeltaPtVsEtaSim->Fill(etaSim,momAtEntry-momAtExit);
245  if(momAtEntry2 >=0 && momAtExit2 >= 0)
246  hDeltaPtVsEtaSim2->Fill(etaSim,momAtEntry2-momAtExit2);
247  }
248  else
249  LogDebug("MuonTrackResidualAnalyzer")<<"NO SimTrack'eta";
250  //
251 
252  // computeResolution(trajectoryBW,muonSimHitsPerId,h1DSimHitRes);
253  // computeResolution(smoothed,muonSimHitsPerId,h1DSimHitRes);
254 
255  }
256 }
257 
259  switch(theEtaRange){
260  case all:
261  return ( abs(eta) <= 2.4 ) ? true : false;
262  case barrel:
263  return ( abs(eta) < 1.1 ) ? true : false;
264  case endcap:
265  return ( abs(eta) >= 1.1 && abs(eta) <= 2.4 ) ? true : false;
266  default:
267  {LogDebug("MuonTrackResidualAnalyzer")<<"No correct Eta range selected!! "; return false;}
268  }
269 }
270 
271 // map the muon simhits by id
272 map<DetId,const PSimHit*>
274  Handle<PSimHitContainer> cscSimhits,
275  Handle<PSimHitContainer> rpcSimhits) {
276 
278 
279  map<DetId,const PSimHit*> hitIdMap;
280  theSimHitContainer.clear();
281 
282  mapMuSimHitsPerId(dtSimhits,hitIdMap);
283  mapMuSimHitsPerId(cscSimhits,hitIdMap);
284  mapMuSimHitsPerId(rpcSimhits,hitIdMap);
285 
286  if(theSimHitContainer.size() >1)
287  stable_sort(theSimHitContainer.begin(),theSimHitContainer.end(),RadiusComparatorInOut(theService->trackingGeometry()));
288 
289  LogDebug("MuonTrackResidualAnalyzer")<<"Sim Hit list";
290  int count=1;
291  for(vector<const PSimHit*>::const_iterator it = theSimHitContainer.begin();
292  it != theSimHitContainer.end(); ++it){
293  LogTrace("MuonTrackResidualAnalyzer")<<count
294  << " "
295  << " Process Type: " << (*it)->processType()
296  << " "
297  << debug.dumpMuonId(DetId( (*it)->detUnitId() ))<<endl;
298  }
299 
300  return hitIdMap;
301 }
302 
303 
305  map<DetId,const PSimHit*> &hitIdMap){
306 
307  for(PSimHitContainer::const_iterator simhit = simhits->begin();
308  simhit != simhits->end(); ++simhit) {
309 
310  if ( abs(simhit->particleType()) != 13 && theSimTkId != simhit->trackId()) continue;
311 
312  theSimHitContainer.push_back(&*simhit);
313  DetId id = DetId(simhit->detUnitId());
314 
315  if(id.subdetId() == MuonSubdetId::DT){
316  DTLayerId lId(id.rawId());
317  id = DetId(lId.rawId());
318  }
319 
320  map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(id);
321 
322  if (it == hitIdMap.end() )
323  hitIdMap[id] = &*simhit;
324  else
325  LogDebug("MuonTrackResidualAnalyzer")<<"TWO muons in the same sensible volume!!";
326 
327  ++theMuonSimHitNumberPerEvent;
328  }
329 }
330 
331 
333  map<DetId,const PSimHit*> &hitIdMap,
335 
337 
338  for(Trajectory::DataContainer::const_iterator datum = data.begin();
339  datum != data.end(); ++datum){
340 
341  GlobalPoint fitPoint = datum->updatedState().globalPosition();
342 
343  // FIXME!
344  // double errX = datum->updatedState().cartesianError().matrix()[0][0];
345  // double errY = datum->updatedState().cartesianError().matrix()[1][1];
346  // double errZ = datum->updatedState().cartesianError().matrix()[2][2];
347  //
348  double errX = datum->updatedState().localError().matrix_old()[3][3];
349  double errY = datum->updatedState().localError().matrix_old()[4][4];
350  double errZ = 1.;
351 
352  map<DetId,const PSimHit*>::const_iterator it = hitIdMap.find(datum->recHit()->geographicalId());
353 
354 
355  if (it == hitIdMap.end() ) continue; // FIXME! Put a counter
356 
357  const PSimHit* simhit = it->second;
358 
359  LocalPoint simHitPoint = simhit->localPosition();
360 
361  const GeomDet* geomDet = theService->trackingGeometry()->idToDet(DetId(simhit->detUnitId()));
362 
363  LocalPoint fitLocalPoint = geomDet->toLocal(fitPoint);
364 
365  LocalVector diff = fitLocalPoint-simHitPoint;
366 
367 
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;
375 
376  cout << "Eta direction: "<< simhit->momentumAtEntry().eta() <<" eta position: " << simHitPoint.eta() << endl;
377  cout << "Phi direction: "<< simhit->momentumAtEntry().phi() <<" phi position: " << simHitPoint.phi() << endl;
378 
379 
380  histos->Fill( simHitPoint.x(), simHitPoint.y(), simHitPoint.z(),
381  diff.x(), diff.y(), diff.z(),
382  errX, errY, errZ,
383  simhit->momentumAtEntry().eta(), simhit->momentumAtEntry().phi());
384  // simHitPoint.eta(), simHitPoint.phi() ); // FIXME!
385  }
386 
387 
388 }
#define LogDebug(id)
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.
Definition: DQMStore.cc:519
LocalVector momentumAtEntry() const
The momentum of the track that produced the hit, at entry point.
Definition: PSimHit.h:47
void cd(void)
go to top directory (ie. root)
Definition: DQMStore.cc:214
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
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")
Definition: DQMStore.cc:1883
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
T y() const
Definition: PV3DBase.h:57
#define abs(x)
Definition: mlp_lapack.h:159
def replace
Definition: linker.py:10
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:64
MuonTrackResidualAnalyzer(const edm::ParameterSet &pset)
Constructor.
T eta() const
std::string dumpMuonId(const DetId &id) const
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
DataContainer const & measurements() const
Definition: Trajectory.h:169
std::vector< TrajectoryMeasurement > DataContainer
Definition: Trajectory.h:42
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Local3DPoint localPosition() const
Definition: PSimHit.h:44
T z() const
Definition: PV3DBase.h:58
tuple result
Definition: query.py:137
tuple pset
Definition: CrabTask.py:85
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
#define LogTrace(id)
DQMStore * dbe_
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
tuple out
Definition: dbtoconf.py:99
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: DetId.h:20
dictionary histos
Definition: combine.py:3
std::string const & label() const
Definition: InputTag.h:25
T eta() const
Definition: PV3DBase.h:70
std::string const & process() const
Definition: InputTag.h:29
virtual ~MuonTrackResidualAnalyzer()
Destructor.
tuple cout
Definition: gather_cfg.py:41
static const int DT
Definition: MuonSubdetId.h:14
void Fill(double x, double y, double z, double dx, double dy, double dz, double errx, double erry, double errz, double eta, double phi)
Definition: Histograms.h:340
void showDirStructure(void) const
Definition: DQMStore.cc:2539
#define debug
Definition: MEtoEDMFormat.h:34
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:647
T x() const
Definition: PV3DBase.h:56
std::string const & instance() const
Definition: InputTag.h:26
void computeResolution(Trajectory &trajectory, std::map< DetId, const PSimHit * > &hitIdMap, HResolution1DRecHit *histos)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:237
unsigned int detUnitId() const
Definition: PSimHit.h:93