File indexing completed on 2024-04-06 12:26:14
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <memory>
0018 #include <fstream>
0019 #include <sys/time.h>
0020 #include <string>
0021 #include <sstream>
0022 #include <iostream>
0023 #include <iomanip>
0024
0025
0026 #include "TFile.h"
0027 #include "TH1F.h"
0028
0029
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/ESHandle.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036
0037 #include <DataFormats/GEMRecHit/interface/ME0SegmentCollection.h>
0038 #include <DataFormats/GEMRecHit/interface/ME0RecHitCollection.h>
0039 #include "DataFormats/GEMDigi/interface/ME0DigiPreRecoCollection.h"
0040 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0041 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0042
0043 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0044 #include <Geometry/GEMGeometry/interface/ME0EtaPartition.h>
0045 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0046 #include <DataFormats/MuonDetId/interface/ME0DetId.h>
0047 #include "DataFormats/Math/interface/deltaPhi.h"
0048
0049
0050
0051
0052
0053 class TestME0SegmentAnalyzer : public edm::one::EDAnalyzer<> {
0054 public:
0055 explicit TestME0SegmentAnalyzer(const edm::ParameterSet&);
0056 ~TestME0SegmentAnalyzer();
0057
0058 private:
0059 virtual void analyze(const edm::Event&, const edm::EventSetup&);
0060
0061
0062 edm::ESGetToken<ME0Geometry, MuonGeometryRecord> me0Geom_Token;
0063
0064 edm::EDGetTokenT<ME0SegmentCollection> ME0Segment_Token;
0065 edm::EDGetTokenT<ME0RecHitCollection> ME0RecHit_Token;
0066 edm::EDGetTokenT<ME0DigiPreRecoCollection> ME0Digi_Token;
0067
0068 std::string rootFileName;
0069
0070 std::unique_ptr<TFile> outputfile;
0071
0072 std::unique_ptr<TH1F> ME0_recdR;
0073 std::unique_ptr<TH1F> ME0_recdPhi;
0074 std::unique_ptr<TH1F> ME0_segdR;
0075 std::unique_ptr<TH1F> ME0_segdPhi;
0076
0077 std::unique_ptr<TH1F> ME0_fitchi2;
0078 std::unique_ptr<TH1F> ME0_rhmult;
0079 std::unique_ptr<TH1F> ME0_rhmultb;
0080 std::unique_ptr<TH1F> ME0_sgmult;
0081 std::unique_ptr<TH1F> ME0_shtime;
0082 std::unique_ptr<TH1F> ME0_rhtime;
0083 std::unique_ptr<TH1F> ME0_sgtime;
0084 std::unique_ptr<TH1F> ME0_sgterr;
0085
0086 std::unique_ptr<TH1F> ME0_Residuals_x;
0087 std::unique_ptr<TH1F> ME0_Residuals_l1_x;
0088 std::unique_ptr<TH1F> ME0_Residuals_l2_x;
0089 std::unique_ptr<TH1F> ME0_Residuals_l3_x;
0090 std::unique_ptr<TH1F> ME0_Residuals_l4_x;
0091 std::unique_ptr<TH1F> ME0_Residuals_l5_x;
0092 std::unique_ptr<TH1F> ME0_Residuals_l6_x;
0093 std::unique_ptr<TH1F> ME0_Pull_x;
0094 std::unique_ptr<TH1F> ME0_Pull_l1_x;
0095 std::unique_ptr<TH1F> ME0_Pull_l2_x;
0096 std::unique_ptr<TH1F> ME0_Pull_l3_x;
0097 std::unique_ptr<TH1F> ME0_Pull_l4_x;
0098 std::unique_ptr<TH1F> ME0_Pull_l5_x;
0099 std::unique_ptr<TH1F> ME0_Pull_l6_x;
0100 std::unique_ptr<TH1F> ME0_Residuals_y;
0101 std::unique_ptr<TH1F> ME0_Residuals_l1_y;
0102 std::unique_ptr<TH1F> ME0_Residuals_l2_y;
0103 std::unique_ptr<TH1F> ME0_Residuals_l3_y;
0104 std::unique_ptr<TH1F> ME0_Residuals_l4_y;
0105 std::unique_ptr<TH1F> ME0_Residuals_l5_y;
0106 std::unique_ptr<TH1F> ME0_Residuals_l6_y;
0107 std::unique_ptr<TH1F> ME0_Pull_y;
0108 std::unique_ptr<TH1F> ME0_Pull_l1_y;
0109 std::unique_ptr<TH1F> ME0_Pull_l2_y;
0110 std::unique_ptr<TH1F> ME0_Pull_l3_y;
0111 std::unique_ptr<TH1F> ME0_Pull_l4_y;
0112 std::unique_ptr<TH1F> ME0_Pull_l5_y;
0113 std::unique_ptr<TH1F> ME0_Pull_l6_y;
0114 };
0115
0116
0117
0118
0119
0120
0121 TestME0SegmentAnalyzer::TestME0SegmentAnalyzer(const edm::ParameterSet& iConfig)
0122
0123 {
0124
0125 me0Geom_Token = esConsumes();
0126 ME0Segment_Token = consumes<ME0SegmentCollection>(edm::InputTag("me0Segments", "", "ME0RERECO"));
0127 ME0RecHit_Token = consumes<ME0RecHitCollection>(edm::InputTag("me0RecHits"));
0128 ME0Digi_Token = consumes<ME0DigiPreRecoCollection>(edm::InputTag("simMuonME0PseudoReDigis"));
0129
0130 rootFileName = iConfig.getUntrackedParameter<std::string>("RootFileName");
0131 outputfile.reset(TFile::Open(rootFileName.c_str(), "RECREATE"));
0132
0133 ME0_recdR = std::unique_ptr<TH1F>(new TH1F("rechitdR", "rechidR", 50, -10., 10.));
0134 ME0_recdPhi = std::unique_ptr<TH1F>(new TH1F("rechitdphi", "rechidphi", 50, -0.005, 0.005));
0135 ME0_segdR = std::unique_ptr<TH1F>(new TH1F("segmentdR", "segmentdR", 50, -10., 10.));
0136 ME0_segdPhi = std::unique_ptr<TH1F>(new TH1F("segmentdphi", "segmentdphi", 50, -0.1, 0.1));
0137 ME0_fitchi2 = std::unique_ptr<TH1F>(new TH1F("chi2Vsndf", "chi2Vsndf", 50, 0., 100.));
0138 ME0_rhmult = std::unique_ptr<TH1F>(new TH1F("rhmulti", "rhmulti", 11, -0.5, 10.5));
0139 ME0_rhmultb = std::unique_ptr<TH1F>(new TH1F("rhmultib", "rhmultib", 11, -0.5, 10.5));
0140 ME0_sgmult = std::unique_ptr<TH1F>(new TH1F("sgmult", "sgmult", 11, -0.5, 10.5));
0141 ME0_rhtime = std::unique_ptr<TH1F>(new TH1F("rhtime", "rhtime", 100, -125., 125.));
0142 ME0_sgtime = std::unique_ptr<TH1F>(new TH1F("sgtime", "sgtime", 100, -125., 125.));
0143 ME0_sgterr = std::unique_ptr<TH1F>(new TH1F("sgterr", "sgterr", 100, 0., 10.));
0144 ME0_Residuals_x = std::unique_ptr<TH1F>(new TH1F("xME0Res", "xME0Res", 100, -0.5, 0.5));
0145 ME0_Residuals_l1_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l1", "xME0Res_l1", 100, -0.5, 0.5));
0146 ME0_Residuals_l2_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l2", "xME0Res_l2", 100, -0.5, 0.5));
0147 ME0_Residuals_l3_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l3", "xME0Res_l3", 100, -0.5, 0.5));
0148 ME0_Residuals_l4_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l4", "xME0Res_l4", 100, -0.5, 0.5));
0149 ME0_Residuals_l5_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l5", "xME0Res_l5", 100, -0.5, 0.5));
0150 ME0_Residuals_l6_x = std::unique_ptr<TH1F>(new TH1F("xME0Res_l6", "xME0Res_l6", 100, -0.5, 0.5));
0151 ME0_Pull_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull", "xME0Pull", 100, -5., 5.));
0152 ME0_Pull_l1_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l1", "xME0Pull_l1", 100, -5., 5.));
0153 ME0_Pull_l2_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l2", "xME0Pull_l2", 100, -5., 5.));
0154 ME0_Pull_l3_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l3", "xME0Pull_l3", 100, -5., 5.));
0155 ME0_Pull_l4_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l4", "xME0Pull_l4", 100, -5., 5.));
0156 ME0_Pull_l5_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l5", "xME0Pull_l5", 100, -5., 5.));
0157 ME0_Pull_l6_x = std::unique_ptr<TH1F>(new TH1F("xME0Pull_l6", "xME0Pull_l6", 100, -5., 5.));
0158 ME0_Residuals_y = std::unique_ptr<TH1F>(new TH1F("yME0Res", "yME0Res", 100, -5., 5.));
0159 ME0_Residuals_l1_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l1", "yME0Res_l1", 100, -5., 5.));
0160 ME0_Residuals_l2_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l2", "yME0Res_l2", 100, -5., 5.));
0161 ME0_Residuals_l3_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l3", "yME0Res_l3", 100, -5., 5.));
0162 ME0_Residuals_l4_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l4", "yME0Res_l4", 100, -5., 5.));
0163 ME0_Residuals_l5_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l5", "yME0Res_l5", 100, -5., 5.));
0164 ME0_Residuals_l6_y = std::unique_ptr<TH1F>(new TH1F("yME0Res_l6", "yME0Res_l6", 100, -5., 5.));
0165 ME0_Pull_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull", "yME0Pull", 100, -5., 5.));
0166 ME0_Pull_l1_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l1", "yME0Pull_l1", 100, -5., 5.));
0167 ME0_Pull_l2_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l2", "yME0Pull_l2", 100, -5., 5.));
0168 ME0_Pull_l3_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l3", "yME0Pull_l3", 100, -5., 5.));
0169 ME0_Pull_l4_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l4", "yME0Pull_l4", 100, -5., 5.));
0170 ME0_Pull_l5_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l5", "yME0Pull_l5", 100, -5., 5.));
0171 ME0_Pull_l6_y = std::unique_ptr<TH1F>(new TH1F("yME0Pull_l6", "yME0Pull_l6", 100, -5., 5.));
0172 }
0173
0174 TestME0SegmentAnalyzer::~TestME0SegmentAnalyzer() {
0175 ME0_recdR->Write();
0176 ME0_recdPhi->Write();
0177 ME0_segdR->Write();
0178 ME0_segdPhi->Write();
0179 ME0_fitchi2->Write();
0180 ME0_rhmult->Write();
0181 ME0_rhmultb->Write();
0182 ME0_sgmult->Write();
0183 ME0_rhtime->Write();
0184 ME0_sgtime->Write();
0185 ME0_sgterr->Write();
0186 ME0_Residuals_x->Write();
0187 ME0_Residuals_l1_x->Write();
0188 ME0_Residuals_l2_x->Write();
0189 ME0_Residuals_l3_x->Write();
0190 ME0_Residuals_l4_x->Write();
0191 ME0_Residuals_l5_x->Write();
0192 ME0_Residuals_l6_x->Write();
0193 ME0_Pull_x->Write();
0194 ME0_Pull_l1_x->Write();
0195 ME0_Pull_l2_x->Write();
0196 ME0_Pull_l3_x->Write();
0197 ME0_Pull_l4_x->Write();
0198 ME0_Pull_l5_x->Write();
0199 ME0_Pull_l6_x->Write();
0200 ME0_Residuals_y->Write();
0201 ME0_Residuals_l1_y->Write();
0202 ME0_Residuals_l2_y->Write();
0203 ME0_Residuals_l3_y->Write();
0204 ME0_Residuals_l4_y->Write();
0205 ME0_Residuals_l5_y->Write();
0206 ME0_Residuals_l6_y->Write();
0207 ME0_Pull_y->Write();
0208 ME0_Pull_l1_y->Write();
0209 ME0_Pull_l2_y->Write();
0210 ME0_Pull_l3_y->Write();
0211 ME0_Pull_l4_y->Write();
0212 ME0_Pull_l5_y->Write();
0213 ME0_Pull_l6_y->Write();
0214 }
0215
0216
0217
0218
0219
0220
0221 void TestME0SegmentAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0222 auto me0Geom = iSetup.getHandle(me0Geom_Token);
0223
0224
0225
0226
0227 edm::Handle<ME0SegmentCollection> me0Segment;
0228 iEvent.getByToken(ME0Segment_Token, me0Segment);
0229 edm::Handle<ME0RecHitCollection> me0RecHit;
0230 iEvent.getByToken(ME0RecHit_Token, me0RecHit);
0231 edm::Handle<ME0DigiPreRecoCollection> me0Digi;
0232 iEvent.getByToken(ME0Digi_Token, me0Digi);
0233
0234
0235 ME0DigiPreRecoCollection::DigiRangeIterator me0dgIt;
0236 for (me0dgIt = me0Digi->begin(); me0dgIt != me0Digi->end(); ++me0dgIt) {
0237 const ME0DetId me0Id = (*me0dgIt).first;
0238 const ME0DigiPreRecoCollection::Range& digiRange = (*me0dgIt).second;
0239 std::cout << " Original DIGI DET ID " << me0Id << " # digis = " << (digiRange.second - digiRange.first + 1)
0240 << std::endl;
0241
0242
0243 for (ME0DigiPreRecoCollection::const_iterator digi = digiRange.first; digi != digiRange.second; digi++) {
0244 std::cout << "x= " << digi->x() << " , y= " << digi->y() << " timing " << digi->tof() << std::endl;
0245 }
0246 }
0247
0248 std::cout << "Number of rec hit " << me0RecHit->size() << std::endl;
0249 float rlmax = 0.;
0250 float rlmin = 999999.;
0251 float plmin = 999999.;
0252 float plmax = 999999.;
0253 int lmin = 100;
0254 int lmax = 0;
0255 for (auto recHit = me0RecHit->begin(); recHit != me0RecHit->end(); recHit++) {
0256 ME0DetId id = recHit->me0Id();
0257 auto roll = me0Geom->etaPartition(id);
0258 std::cout << " Original ME0DetID " << id << " Timing " << recHit->tof() << std::endl;
0259 int layer = id.layer();
0260 if (layer < lmin) {
0261 lmin = layer;
0262 rlmin = (roll->toGlobal(recHit->localPosition())).perp();
0263 plmin = (roll->toGlobal(recHit->localPosition())).barePhi();
0264 }
0265 if (layer > lmax) {
0266 lmax = layer;
0267 rlmax = (roll->toGlobal(recHit->localPosition())).perp();
0268 plmax = (roll->toGlobal(recHit->localPosition())).barePhi();
0269 }
0270 }
0271 std::cout << " Radius max min and delta " << rlmax << " " << rlmin << " " << rlmax - rlmin << std::endl;
0272 std::cout << " Phi max min and delta " << plmax << " " << plmin << " " << plmax - plmin << std::endl;
0273 ME0_recdR->Fill(rlmax - rlmin);
0274 ME0_recdPhi->Fill(fabs(deltaPhi(plmax, plmin)));
0275
0276 std::cout << "Number of Segments " << me0Segment->size() << std::endl;
0277 ME0_sgmult->Fill(me0Segment->size());
0278 float hmax = 0;
0279 for (auto me0s = me0Segment->begin(); me0s != me0Segment->end(); me0s++) {
0280
0281 ME0DetId id = me0s->me0DetId();
0282
0283 ME0_sgtime->Fill(me0s->time());
0284 ME0_sgterr->Fill(me0s->timeErr());
0285 std::cout << " Original ME0DetID " << id << std::endl;
0286 auto roll = me0Geom->etaPartition(id);
0287 std::cout << " Global Segment Position " << roll->toGlobal(me0s->localPosition()) << std::endl;
0288 auto segLP = me0s->localPosition();
0289 auto segLD = me0s->localDirection();
0290 std::cout << " Local Direction theta = " << segLD.theta() << " phi=" << segLD.phi() << std::endl;
0291 ME0_fitchi2->Fill(me0s->chi2() * 1.0 / me0s->degreesOfFreedom());
0292 std::cout << " Chi2 = " << me0s->chi2() << " ndof = " << me0s->degreesOfFreedom()
0293 << " ==> chi2/ndof = " << me0s->chi2() * 1.0 / me0s->degreesOfFreedom() << " Timing " << me0s->time()
0294 << " +- " << me0s->timeErr() << std::endl;
0295
0296 auto me0rhs = me0s->specificRecHits();
0297 std::cout << " ME0 Ensemble Det Id " << id << " Number of RecHits " << me0rhs.size() << std::endl;
0298 ME0_rhmult->Fill(me0rhs.size());
0299 if (me0rhs.size() > hmax)
0300 hmax = me0rhs.size();
0301
0302 for (auto rh = me0rhs.begin(); rh != me0rhs.end(); rh++) {
0303 auto me0id = rh->me0Id();
0304 ME0_rhtime->Fill(rh->tof());
0305 auto rhr = me0Geom->etaPartition(me0id);
0306 auto rhLP = rh->localPosition();
0307 auto erhLEP = rh->localPositionError();
0308 auto rhGP = rhr->toGlobal(rhLP);
0309 auto rhLPSegm = roll->toLocal(rhGP);
0310 float xe = segLP.x() + segLD.x() * (rhLPSegm.z() - segLP.z()) / segLD.z();
0311 float ye = segLP.y() + segLD.y() * (rhLPSegm.z() - segLP.z()) / segLD.z();
0312 float ze = rhLPSegm.z();
0313 LocalPoint extrPoint(xe, ye, ze);
0314 auto extSegm = rhr->toLocal(roll->toGlobal(extrPoint));
0315 std::cout << " ME0 Layer Id " << rh->me0Id() << " error on the local point " << erhLEP
0316 << "\n-> Ensemble Rest Frame RH local position " << rhLPSegm << " Segment extrapolation "
0317 << extrPoint << "\n-> Layer Rest Frame RH local position " << rhLP << " Segment extrapolation "
0318 << extSegm << "\n-> Global Position rechit " << rhGP << " Segm Extrapolation "
0319 << roll->toGlobal(extrPoint) << "\n-> Timing " << rh->tof() << std::endl;
0320 ME0_Residuals_x->Fill(rhLP.x() - extSegm.x());
0321 ME0_Residuals_y->Fill(rhLP.y() - extSegm.y());
0322 ME0_Pull_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0323 ME0_Pull_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0324 switch (me0id.layer()) {
0325 case 1:
0326 ME0_Residuals_l1_x->Fill(rhLP.x() - extSegm.x());
0327 ME0_Residuals_l1_y->Fill(rhLP.y() - extSegm.y());
0328 ME0_Pull_l1_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0329 ME0_Pull_l1_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0330 break;
0331 case 2:
0332 ME0_Residuals_l2_x->Fill(rhLP.x() - extSegm.x());
0333 ME0_Residuals_l2_y->Fill(rhLP.y() - extSegm.y());
0334 ME0_Pull_l2_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0335 ME0_Pull_l2_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0336 break;
0337 case 3:
0338 ME0_Residuals_l3_x->Fill(rhLP.x() - extSegm.x());
0339 ME0_Residuals_l3_y->Fill(rhLP.y() - extSegm.y());
0340 ME0_Pull_l3_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0341 ME0_Pull_l3_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0342 break;
0343 case 4:
0344 ME0_Residuals_l4_x->Fill(rhLP.x() - extSegm.x());
0345 ME0_Residuals_l4_y->Fill(rhLP.y() - extSegm.y());
0346 ME0_Pull_l4_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0347 ME0_Pull_l4_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0348 break;
0349 case 5:
0350 ME0_Residuals_l5_x->Fill(rhLP.x() - extSegm.x());
0351 ME0_Residuals_l5_y->Fill(rhLP.y() - extSegm.y());
0352 ME0_Pull_l5_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0353 ME0_Pull_l5_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0354 break;
0355 case 6:
0356 ME0_Residuals_l6_x->Fill(rhLP.x() - extSegm.x());
0357 ME0_Residuals_l6_y->Fill(rhLP.y() - extSegm.y());
0358 ME0_Pull_l6_x->Fill((rhLP.x() - extSegm.x()) / sqrt(erhLEP.xx()));
0359 ME0_Pull_l6_y->Fill((rhLP.y() - extSegm.y()) / sqrt(erhLEP.yy()));
0360 break;
0361 default:
0362 std::cout << " Unphysical ME0 layer " << me0id << std::endl;
0363 }
0364 }
0365 std::cout << "\n" << std::endl;
0366 }
0367 std::cout << "\n" << std::endl;
0368 ME0_rhmultb->Fill(hmax);
0369 }
0370
0371
0372 DEFINE_FWK_MODULE(TestME0SegmentAnalyzer);