CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AlignmentMonitorMuonVsCurvature.cc
Go to the documentation of this file.
1 /*
2  * Package: CommonAlignmentProducer
3  * Class : AlignmentMonitorMuonVsCurvature
4  *
5  * Original Author: Jim Pivarski
6  * Created: Fri Feb 19 21:45:02 CET 2010
7  *
8  * $Id:$
9  */
10 
14 
21 
22 #include <sstream>
23 #include "TProfile.h"
24 #include "TH2F.h"
25 #include "TH1F.h"
26 
27 
29 {
30 public:
33 
34  void book() override;
35 
36  void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& iTrajTracks) override;
38 
39  void afterAlignment(const edm::EventSetup &iSetup) override {}
40 
41 private:
42 
44  double m_minTrackPt;
45  double m_minTrackP;
50  double m_maxDxy;
54  int m_layer;
56  bool m_doDT;
57  bool m_doCSC;
58 
59  enum {
60  kDeltaX = 0,
63  };
64 
65  // DT [wheel][station][sector]
68 
69  // CSC [endcap*station][ring][chamber]
72 
75 };
76 
77 
79 : AlignmentMonitorBase(cfg, "AlignmentMonitorMuonVsCurvature")
80 , m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag"))
81 , m_minTrackPt(cfg.getParameter<double>("minTrackPt"))
82 , m_minTrackP(cfg.getParameter<double>("minTrackP"))
83 , m_minTrackerHits(cfg.getParameter<int>("minTrackerHits"))
84 , m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2"))
85 , m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC"))
86 , m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers"))
87 , m_maxDxy(cfg.getParameter<double>("maxDxy"))
88 , m_minDT13Hits(cfg.getParameter<int>("minDT13Hits"))
89 , m_minDT2Hits(cfg.getParameter<int>("minDT2Hits"))
90 , m_minCSCHits(cfg.getParameter<int>("minCSCHits"))
91 , m_layer(cfg.getParameter<int>("layer"))
92 , m_propagator(cfg.getParameter<std::string>("propagator"))
93 , m_doDT(cfg.getParameter<bool>("doDT"))
94 , m_doCSC(cfg.getParameter<bool>("doCSC"))
95 {}
96 
97 
99 {
100  // DT
101  std::string wheelname[5] = {"wheelm2_", "wheelm1_", "wheelz_", "wheelp1_", "wheelp2_"};
102  if (m_doDT)
103  for (int wheel = -2; wheel <=2 ; wheel++)
104  for (int station = 1; station <= 4; station++)
105  for (int sector = 1; sector <= 14; sector++)
106  {
107  if (station != 4 && sector > 12) continue;
108 
109  char stationname[20];
110  sprintf(stationname,"st%d_", station);
111 
112  char sectorname[20];
113  sprintf(sectorname,"sector%02d_", sector);
114 
115  for (int component = 0; component < kNumComponents; component++)
116  {
117  std::stringstream th2f_name, tprofile_name;
118  th2f_name << "th2f_" << wheelname[wheel+2] <<stationname<< sectorname;
119  tprofile_name << "tprofile_" << wheelname[wheel+2] <<stationname<< sectorname;
120 
121  double yminmax = 50., xminmax = 0.05;
122  if (m_minTrackPt>0.) xminmax = 1./m_minTrackPt;
123  int ynbins = 50;
124  if (component == kDeltaX) {
125  th2f_name << "deltax";
126  tprofile_name << "deltax";
127  }
128  else if (component == kDeltaDxDz) {
129  th2f_name << "deltadxdz";
130  tprofile_name << "deltadxdz";
131  }
132 
133  th2f_wheel_st_sector[wheel+2][station-1][sector-1][component] =
134  book2D("/iterN/", th2f_name.str().c_str(), "", 30, -xminmax, xminmax, ynbins, -yminmax, yminmax);
135  tprofile_wheel_st_sector[wheel+2][station-1][sector-1][component] =
136  bookProfile("/iterN/", tprofile_name.str().c_str(), "", 30, -xminmax, xminmax);
137  }
138  }
139 
140  // CSC
141  std::string stname[8] = {"Ep_S1_", "Ep_S2_", "Ep_S3_", "Ep_S4_", "Em_S1_", "Em_S2_", "Em_S3_", "Em_S4_"};
142  if (m_doCSC)
143  for (int station = 0; station < 8; station++)
144  for (int ring = 1; ring <= 3; ring++)
145  for (int chamber = 1; chamber <= 36; chamber++)
146  {
147  int st = station%4+1;
148  if (st > 1 && ring > 2) continue; // only station 1 has more then 2 rings
149  if (st > 1 && ring == 1 && chamber > 18) continue; // ring 1 stations 1,2,3 have 18 chambers
150 
151  char ringname[20];
152  sprintf(ringname,"R%d_", ring);
153 
154  char chname[20];
155  sprintf(chname,"C%02d_", chamber);
156 
157  for (int component = 0; component < kNumComponents; component++)
158  {
159  std::stringstream componentname;
160  double yminmax = 50., xminmax = 0.05;
161  if (m_minTrackPt>0.) xminmax = 1./m_minTrackPt;
162  if (ring == 1) xminmax *= 0.5;
163  if (component == kDeltaX) {
164  componentname << "deltax";
165  }
166  else if (component == kDeltaDxDz) {
167  componentname << "deltadxdz";
168  }
169 
170  std::stringstream th2f_name, tprofile_name;
171  th2f_name << "th2f_" << stname[station] << ringname << chname << componentname.str();
172  tprofile_name << "tprofile_" << stname[station] << ringname << chname << componentname.str();
173 
174  th2f_st_ring_chamber[station][ring-1][chamber-1][component] =
175  book2D("/iterN/", th2f_name.str().c_str(), "", 30, -xminmax, xminmax, 100, -yminmax, yminmax);
176  tprofile_st_ring_chamber[station][ring-1][chamber-1][component] =
177  bookProfile("/iterN/", tprofile_name.str().c_str(), "", 30, -xminmax, xminmax);
178  }
179  }
180 
181  th1f_trackerRedChi2 = book1D("/iterN/", "trackerRedChi2", "Refit tracker reduced chi^2", 100, 0., 30.);
182  th1f_trackerRedChi2Diff = book1D("/iterN/", "trackerRedChi2Diff", "Fit-minus-refit tracker reduced chi^2", 100, -5., 5.);
183 }
184 
185 
187 {
189  iSetup.get<GlobalTrackingGeometryRecord>().get(globalGeometry);
190 
192  iEvent.getByLabel(m_beamSpotTag, beamSpot);
193 
194  if (m_muonCollectionTag.label().empty()) // use trajectories
195  {
196  for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end(); ++trajtrack)
197  {
198  const Trajectory* traj = (*trajtrack).first;
199  const reco::Track* track = (*trajtrack).second;
200 
201  if (track->pt() > m_minTrackPt && track->p() > m_minTrackP && fabs(track->dxy(beamSpot->position())) < m_maxDxy )
202  {
203  MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, traj, track, pNavigator(), 1000.);
204  processMuonResidualsFromTrack(muonResidualsFromTrack, traj );
205  } // end if track pT is within range
206  } // end loop over tracks
207  }
208  else
209  {
211  iEvent.getByLabel(m_muonCollectionTag, muons);
212 
213  for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon)
214  {
215  if ( !(muon->isTrackerMuon() && muon->innerTrack().isNonnull() ) ) continue;
216 
217  if (m_minTrackPt < muon->pt() && m_minTrackP < muon->p() && fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy)
218  {
219  MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
220  processMuonResidualsFromTrack(muonResidualsFromTrack);
221  }
222  }
223  }
224 }
225 
226 
228 {
229  if (mrft.trackerNumHits() < m_minTrackerHits) return;
230  if (!m_allowTIDTEC && mrft.contains_TIDTEC()) return;
231 
232  int nMuChambers = 0;
233  std::vector<DetId> chamberIds = mrft.chamberIds();
234  for (unsigned ch=0; ch < chamberIds.size(); ch++) if (chamberIds[ch].det() == DetId::Muon) nMuChambers++;
235  if (nMuChambers < m_minNCrossedChambers ) return;
236 
237  th1f_trackerRedChi2->Fill(mrft.trackerRedChi2());
239 
240  if (mrft.normalizedChi2() > m_maxTrackerRedChi2) return;
241 
242  double qoverpt = mrft.getTrack()->charge() / mrft.getTrack()->pt();
243  double qoverpz = 0.;
244  if (fabs(mrft.getTrack()->pz()) > 0.01) qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());
245 
246  for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId)
247  {
248  if (chamberId->det() != DetId::Muon ) continue;
249 
250  if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT)
251  {
252  DTChamberId dtid(chamberId->rawId());
254 
255  if (dt13 != NULL && dt13->numHits() >= m_minDT13Hits)
256  {
257  int wheel = dtid.wheel() + 2;
258  int station = dtid.station() -1;
259  int sector = dtid.sector() - 1;
260 
261  double resid_x = 10. * dt13->global_residual();
262  double resid_dxdz = 1000. * dt13->global_resslope();
263 
264  if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.)
265  {
266  th2f_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
267  tprofile_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
268  th2f_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
269  tprofile_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
270  }
271  } // if it's a good segment
272  } // if DT
273 
274  if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC)
275  {
276  CSCDetId cscid(chamberId->rawId());
278 
279  if (csc != NULL && csc->numHits() >= m_minCSCHits)
280  {
281  int station = 4*cscid.endcap() + cscid.station() - 5;
282  int ring = cscid.ring() - 1;
283  if (cscid.station()==1 && cscid.ring()==4) ring = 0; // join ME1/a to ME1/b
284  int chamber = cscid.chamber() - 1;
285 
286  double resid_x = 10. * csc->global_residual();
287  double resid_dxdz = 1000. * csc->global_resslope();
288 
289  if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.)
290  {
291  th2f_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
292  tprofile_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
293  th2f_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
294  tprofile_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
295  }
296  } // if it's a good segment
297  } // if CSC
298 
299  } // end loop over chamberIds
300 }
301 
double p() const
momentum vector magnitude
Definition: TrackBase.h:602
void afterAlignment(const edm::EventSetup &iSetup) override
AlignableNavigator * pNavigator()
tuple cfg
Definition: looper.py:293
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:548
TH2F * book2D(std::string dir, std::string name, std::string title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
TH2F * th2f_st_ring_chamber[8][3][36][kNumComponents]
#define NULL
Definition: scimark2.h:8
TProfile * tprofile_st_ring_chamber[8][3][36][kNumComponents]
TProfile * bookProfile(std::string dir, std::string name, std::string title, int nchX, double lowX, double highX, int nchY=1, double lowY=0., double highY=0., const char *option="s")
AlignmentMonitorMuonVsCurvature(const edm::ParameterSet &cfg)
const std::vector< DetId > chamberIds() const
void event(const edm::Event &iEvent, const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection &iTrajTracks) override
Called for each event (by &quot;run()&quot;): may be reimplemented.
int iEvent
Definition: GenABIO.cc:230
static const int CSC
Definition: MuonSubdetId.h:13
double pt() const
track transverse momentum
Definition: TrackBase.h:608
const reco::Track * getTrack()
void book() override
Book or retrieve histograms; MUST be reimplemented.
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection
const edm::InputTag m_beamSpotTag
TProfile * tprofile_wheel_st_sector[5][4][14][kNumComponents]
TH1F * book1D(std::string dir, std::string name, std::string title, int nchX, double lowX, double highX)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:420
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:626
const T & get() const
Definition: EventSetup.h:56
std::string const & label() const
Definition: InputTag.h:43
tuple muons
Definition: patZpeak.py:38
void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft, const Trajectory *traj=NULL)
int charge() const
track electric charge
Definition: TrackBase.h:554
#define DEFINE_EDM_PLUGIN(factory, type, name)
static const int DT
Definition: MuonSubdetId.h:12
TH2F * th2f_wheel_st_sector[5][4][14][kNumComponents]
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:578
MuonChamberResidual * chamberResidual(DetId chamberId, int type)