CMS 3D CMS Logo

MuonSeedsAnalyzer.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * $Date: 2008/03/25
5  18:37:05 $
6  * \author G. Mila - INFN Torino
7  */
8 
10 
19 
21 
24 
29 
30 #include <string>
31 
32 
33 using namespace std;
34 using namespace edm;
35 
36 
38  parameters = pSet;
39 
40  theService = new MuonServiceProxy(parameters.getParameter<ParameterSet>("ServiceParameters"));
41 
42  theSeedsCollectionLabel_ = consumes<TrajectorySeedCollection>(parameters.getParameter<InputTag>("SeedCollection"));
43 
44  seedHitBin = parameters.getParameter<int>("RecHitBin");
45  seedHitMin = parameters.getParameter<double>("RecHitMin");
46  seedHitMax = parameters.getParameter<double>("RecHitMax");
47  PhiBin = parameters.getParameter<int>("PhiBin");
48  PhiMin = parameters.getParameter<double>("PhiMin");
49  PhiMax = parameters.getParameter<double>("PhiMax");
50  EtaBin = parameters.getParameter<int>("EtaBin");
51  EtaMin = parameters.getParameter<double>("EtaMin");
52  EtaMax = parameters.getParameter<double>("EtaMax");
53  ThetaBin = parameters.getParameter<int>("ThetaBin");
54  ThetaMin = parameters.getParameter<double>("ThetaMin");
55  ThetaMax = parameters.getParameter<double>("ThetaMax");
56  seedPtBin = parameters.getParameter<int>("seedPtBin");
57  seedPtMin = parameters.getParameter<double>("seedPtMin");
58  seedPtMax = parameters.getParameter<double>("seedPtMax");
59  seedPxyzBin = parameters.getParameter<int>("seedPxyzBin");
60  seedPxyzMin = parameters.getParameter<double>("seedPxyzMin");
61  seedPxyzMax = parameters.getParameter<double>("seedPxyzMax");
62  pErrBin = parameters.getParameter<int>("pErrBin");
63  pErrMin = parameters.getParameter<double>("pErrMin");
64  pErrMax = parameters.getParameter<double>("pErrMax");
65  pxyzErrBin = parameters.getParameter<int>("pxyzErrBin");
66  pxyzErrMin = parameters.getParameter<double>("pxyzErrMin");
67  pxyzErrMax = parameters.getParameter<double>("pxyzErrMax");
68  phiErrBin = parameters.getParameter<int>("phiErrBin");
69  phiErrMin = parameters.getParameter<double>("phiErrMin");
70  phiErrMax = parameters.getParameter<double>("phiErrMax");
71  etaErrBin = parameters.getParameter<int>("etaErrBin");
72  etaErrMin = parameters.getParameter<double>("etaErrMin");
73  etaErrMax = parameters.getParameter<double>("etaErrMax");
74 }
76  delete theService;
77 }
79 
80  ibooker.cd();
81  ibooker.setCurrentFolder("Muons/MuonSeedsAnalyzer");
82 
83  string histname = "NumberOfRecHitsPerSeed_";
84  NumberOfRecHitsPerSeed = ibooker.book1D(histname, "Number of seed recHits", seedHitBin, seedHitMin, seedHitMax);
85 
86  histname = "seedPhi_";
87  seedPhi = ibooker.book1D(histname, "Seed #phi", PhiBin, PhiMin, PhiMax);
88  seedPhi->setAxisTitle("rad");
89 
90  histname = "seedEta_";
91  seedEta = ibooker.book1D(histname, "Seed #eta", EtaBin, EtaMin, EtaMax);
92 
93  histname = "seedTheta_";
94  seedTheta = ibooker.book1D(histname, "Seed #theta", ThetaBin, ThetaMin, ThetaMax);
95  seedTheta->setAxisTitle("rad");
96 
97  histname = "seedPt_";
98  seedPt = ibooker.book1D(histname, "Seed p_{t}", seedPtBin, seedPtMin, seedPtMax);
99  seedPt->setAxisTitle("GeV");
100 
101  histname = "seedPx_";
102  seedPx = ibooker.book1D(histname, "Seed p_{x}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
103  seedPx->setAxisTitle("GeV");
104  histname = "seedPy_";
105  seedPy = ibooker.book1D(histname, "Seed p_{y}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
106  seedPy->setAxisTitle("GeV");
107  histname = "seedPz_";
108  seedPz = ibooker.book1D(histname, "Seed p_{z}", seedPxyzBin, seedPxyzMin, seedPxyzMax);
109  seedPz->setAxisTitle("GeV");
110 
111  histname = "seedPtErrOverPt_";
112  seedPtErr = ibooker.book1D(histname, "Seed p_{t}Err/p_{t}", pErrBin, pErrMin, pErrMax);
113  histname = "seedPtErrOverPtVsPhi_";
114  seedPtErrVsPhi = ibooker.book2D(histname, "Seed p_{t}Err/p_{t} vs #phi", PhiBin, PhiMin, PhiMax, pErrBin, pErrMin, pErrMax);
115  seedPtErrVsPhi->setAxisTitle("rad",2);
116  histname = "seedPtErrOverPtVsEta_";
117  seedPtErrVsEta = ibooker.book2D(histname, "Seed p_{t}Err/p_{t} vs #eta", EtaBin, EtaMin, EtaMax, pErrBin, pErrMin, pErrMax);
118  histname = "seedPtErrOverPtVsPt_";
119  seedPtErrVsPt = ibooker.book2D(histname, "Seed p_{t}Err/p_{t} vs p_{t}", seedPtBin/5, seedPtMin, seedPtMax, pErrBin, pErrMin, pErrMax);
120  seedPtErrVsPt->setAxisTitle("GeV",2);
121  histname = "seedPErrOverP_";
122  seedPErr = ibooker.book1D(histname, "Seed pErr/p", pErrBin, pErrMin, pErrMax);
123  histname = "seedPErrOverPVsPhi_";
124  seedPErrVsPhi = ibooker.book2D(histname, "Seed pErr/p vs #phi", PhiBin, PhiMin, PhiMax, pErrBin, pErrMin, pErrMax);
125  seedPErrVsPhi->setAxisTitle("rad",2);
126  histname = "seedPErrOverPVsEta_";
127  seedPErrVsEta = ibooker.book2D(histname, "Seed pErr/p vs #eta", EtaBin, EtaMin, EtaMax, pErrBin, pErrMin, pErrMax);
128  histname = "seedPErrOverPVsPt_";
129  seedPErrVsPt = ibooker.book2D(histname, "Seed pErr/p vs p_{t}", seedPtBin/5, seedPtMin, seedPtMax, pErrBin, pErrMin, pErrMax);
130  seedPErrVsPt->setAxisTitle("GeV",2);
131 
132  histname = "seedPxErrOverPx_";
133  seedPxErr = ibooker.book1D(histname, "Seed p_{x}Err/p_{x}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
134  histname = "seedPyErrOverPy_";
135  seedPyErr = ibooker.book1D(histname, "Seed p_{y}Err/p_{y}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
136  histname = "seedPzErrOverPz_";
137  seedPzErr = ibooker.book1D(histname, "Seed p_{z}Err/p_{z}", pxyzErrBin, pxyzErrMin, pxyzErrMax);
138 
139 
140  histname = "seedPhiErr_";
141  seedPhiErr = ibooker.book1D(histname, "Seed #phi error", phiErrBin, phiErrMin, phiErrMax);
142 
143  histname = "seedEtaErr_";
144  seedEtaErr = ibooker.book1D(histname, "Seed #eta error", etaErrBin, etaErrMin, etaErrMax);
145 
146 }
147 
149  theService->update(iSetup);
150 
151  // Take the seeds container
153  iEvent.getByToken(theSeedsCollectionLabel_, seeds);
154 
155  // if not valid, skip
156  if (!seeds.isValid()) return;
157 
158  for(TrajectorySeedCollection::const_iterator seed = seeds->begin(); seed != seeds->end(); ++seed){
159  // const TrajectorySeed sd = *seed;
160 
161  // Get the Trajectory State on Det (persistent version of a TSOS) from the seed
162  PTrajectoryStateOnDet pTSOD = seed->startingState();
163 
164  // Transform it in a TrajectoryStateOnSurface
165  DetId seedDetId(pTSOD.detId());
166  const GeomDet* gdet = theService->trackingGeometry()->idToDet( seedDetId );
167  TrajectoryStateOnSurface seedTSOS = trajectoryStateTransform::transientState(pTSOD, &(gdet->surface()), &*(theService)->magneticField());
169  double partialPterror = errors(3,3)*pow(seedTSOS.globalMomentum().x(),2) + errors(4,4)*pow(seedTSOS.globalMomentum().y(),2);
170 
171  LogTrace(metname)<<"[MuonSeedAnalyzer] Filling the histos";
172 
173  // nhits
174  LogTrace(metname)<<"Number od recHits per seed: "<<seed->nHits();
175  NumberOfRecHitsPerSeed->Fill(seed->nHits());
176 
177  // pt
178  LogTrace(metname)<<"seed momentum: "<<seedTSOS.globalMomentum().perp();
179  seedPt->Fill(seedTSOS.globalMomentum().perp());
180 
181  // px
182  LogTrace(metname)<<"seed px: "<<seedTSOS.globalMomentum().x();
183  seedPx->Fill(seedTSOS.globalMomentum().x());
184 
185  // py
186  LogTrace(metname)<<"seed py: "<<seedTSOS.globalMomentum().y();
187  seedPy->Fill(seedTSOS.globalMomentum().y());
188 
189  // pz
190  LogTrace(metname)<<"seed pz: "<<seedTSOS.globalMomentum().z();
191  seedPz->Fill(seedTSOS.globalMomentum().z());
192 
193  // phi
194  LogTrace(metname)<<"seed phi: "<<seedTSOS.globalMomentum().phi();
195  seedPhi->Fill(seedTSOS.globalMomentum().phi());
196 
197  // theta
198  LogTrace(metname)<<"seed theta: "<<seedTSOS.globalMomentum().theta();
199  seedTheta->Fill(seedTSOS.globalMomentum().theta());
200 
201  // eta
202  LogTrace(metname)<<"seed eta: "<<seedTSOS.globalMomentum().eta();
203  seedEta->Fill(seedTSOS.globalMomentum().eta());
204 
205  // pt err
206  LogTrace(metname)<<"seed pt error: "<<sqrt(partialPterror)/seedTSOS.globalMomentum().perp();
207  seedPtErr->Fill(sqrt(partialPterror)/seedTSOS.globalMomentum().perp());
208 
209  // ptErr/pt Vs phi
210  seedPtErrVsPhi->Fill(seedTSOS.globalMomentum().phi(),
211  sqrt(partialPterror)/seedTSOS.globalMomentum().perp());
212  // ptErr/pt Vs eta
213  seedPtErrVsEta->Fill(seedTSOS.globalMomentum().eta(),
214  sqrt(partialPterror)/seedTSOS.globalMomentum().perp());
215  // ptErr/pt Vs pt
216  seedPtErrVsPt->Fill(seedTSOS.globalMomentum().perp(),
217  sqrt(partialPterror)/seedTSOS.globalMomentum().perp());
218 
219  // px err
220  LogTrace(metname)<<"seed px error: "<<sqrt(errors(3,3))/seedTSOS.globalMomentum().x();
221  seedPxErr->Fill(sqrt(errors(3,3))/seedTSOS.globalMomentum().x());
222 
223  // py err
224  LogTrace(metname)<<"seed py error: "<<sqrt(errors(4,4))/seedTSOS.globalMomentum().y();
225  seedPyErr->Fill(sqrt(errors(4,4))/seedTSOS.globalMomentum().y());
226 
227  // pz err
228  LogTrace(metname)<<"seed pz error: "<<sqrt(errors(5,5))/seedTSOS.globalMomentum().z();
229  seedPzErr->Fill(sqrt(errors(5,5))/seedTSOS.globalMomentum().z());
230 
231  // p err
232  LogTrace(metname)<<"seed p error: "<<sqrt(partialPterror+errors(5,5)*pow(seedTSOS.globalMomentum().z(),2))/seedTSOS.globalMomentum().mag();
233  seedPErr->Fill(sqrt(partialPterror+errors(5,5)*pow(seedTSOS.globalMomentum().z(),2))/seedTSOS.globalMomentum().mag());
234 
235  // pErr/p Vs phi
236  seedPErrVsPhi->Fill(seedTSOS.globalMomentum().phi(),
237  sqrt(partialPterror+errors(5,5)*pow(seedTSOS.globalMomentum().z(),2))/seedTSOS.globalMomentum().mag());
238  // pErr/p Vs eta
239  seedPErrVsEta->Fill(seedTSOS.globalMomentum().eta(),
240  sqrt(partialPterror+errors(5,5)*pow(seedTSOS.globalMomentum().z(),2))/seedTSOS.globalMomentum().mag());
241  // pErr/p Vs pt
242  seedPErrVsPt->Fill(seedTSOS.globalMomentum().perp(),
243  sqrt(partialPterror+errors(5,5)*pow(seedTSOS.globalMomentum().z(),2))/seedTSOS.globalMomentum().mag());
244 
245  // phi err
246  LogTrace(metname)<<"seed phi error: "<<sqrt(seedTSOS.curvilinearError().matrix()(2,2));
247  seedPhiErr->Fill(sqrt(seedTSOS.curvilinearError().matrix()(2,2)));
248 
249  // eta err
250  LogTrace(metname)<<"seed eta error: "<<sqrt(seedTSOS.curvilinearError().matrix()(1,1))*abs(sin(seedTSOS.globalMomentum().theta()));
251  seedEtaErr->Fill(sqrt(seedTSOS.curvilinearError().matrix()(1,1))*abs(sin(seedTSOS.globalMomentum().theta())));
252  }
253 }
254 
255 
T perp() const
Definition: PV3DBase.h:72
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
const std::string metname
void cd(void)
Definition: DQMStore.cc:269
const CurvilinearTrajectoryError & curvilinearError() const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
const CartesianTrajectoryError cartesianError() const
T y() const
Definition: PV3DBase.h:63
const double EtaMax[kNumberCalorimeter]
Geom::Theta< T > theta() const
Definition: PV3DBase.h:75
int iEvent
Definition: GenABIO.cc:230
T mag() const
Definition: PV3DBase.h:67
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
unsigned int detId() const
const double EtaMin[kNumberCalorimeter]
bool isValid() const
Definition: HandleBase.h:74
#define LogTrace(id)
const AlgebraicSymMatrix66 & matrix() const
Definition: DetId.h:18
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
TrajectoryStateOnSurface transientState(const PTrajectoryStateOnDet &ts, const Surface *surface, const MagneticField *field)
virtual ~MuonSeedsAnalyzer()
Destructor.
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
T eta() const
Definition: PV3DBase.h:76
void analyze(const edm::Event &, const edm::EventSetup &) override
HLT enums.
const AlgebraicSymMatrix55 & matrix() const
GlobalVector globalMomentum() const
Definition: errors.py:1
T x() const
Definition: PV3DBase.h:62
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: Run.h:42
MuonSeedsAnalyzer(const edm::ParameterSet &)
Constructor.
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override