File indexing completed on 2021-02-14 14:24:59
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "RPCRecHitReader.h"
0010
0011 #include <memory>
0012 #include <fstream>
0013 #include <iostream>
0014 #include <string>
0015 #include <cmath>
0016 #include <cmath>
0017 #include <vector>
0018 #include <iomanip>
0019 #include <set>
0020 #include <stdio.h>
0021
0022 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0026 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0027 #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
0028 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0029
0030 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0031 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0032 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0033 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0034 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0035 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0036 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0037
0038 #include "FWCore/Framework/interface/MakerMacros.h"
0039 #include "FWCore/Framework/interface/Frameworkfwd.h"
0040 #include "FWCore/Framework/interface/Event.h"
0041 #include "FWCore/Framework/interface/EventSetup.h"
0042 #include "FWCore/Framework/interface/ESHandle.h"
0043
0044 #include "TFile.h"
0045 #include "TVector3.h"
0046 #include "TH1F.h"
0047 #include "TH2F.h"
0048 #include "TF1.h"
0049 #include "TF2.h"
0050 #include "TVectorT.h"
0051 #include "TGraph.h"
0052
0053 #include <map>
0054
0055 Double_t linearF(Double_t* x, Double_t* par) {
0056 Double_t y = 0.;
0057 y = par[0] * (*x) + par[1];
0058 return y;
0059 }
0060
0061
0062 RPCRecHitReader::RPCRecHitReader(const edm::ParameterSet& pset) : _phi(0) {
0063
0064 fOutputFileName = pset.getUntrackedParameter<std::string>("HistOutFile");
0065 recHitLabel1 = pset.getUntrackedParameter<std::string>("recHitLabel1");
0066 recHitLabel2 = pset.getUntrackedParameter<std::string>("recHitLabel2");
0067
0068 region = pset.getUntrackedParameter<int>("region", 0);
0069 wheel = pset.getUntrackedParameter<int>("wheel", 1);
0070 sector = pset.getUntrackedParameter<int>("sector", 10);
0071 station = pset.getUntrackedParameter<int>("station");
0072 layer = pset.getUntrackedParameter<int>("layer");
0073 subsector = pset.getUntrackedParameter<int>("subsector");
0074
0075 _trigConfig.clear();
0076
0077 _trigRPC1 = pset.getUntrackedParameter<bool>("trigRPC1");
0078 _trigRPC2 = pset.getUntrackedParameter<bool>("trigRPC2");
0079 _trigRPC3 = pset.getUntrackedParameter<bool>("trigRPC3");
0080 _trigRPC4 = pset.getUntrackedParameter<bool>("trigRPC4");
0081 _trigRPC5 = pset.getUntrackedParameter<bool>("trigRPC5");
0082 _trigRPC6 = pset.getUntrackedParameter<bool>("trigRPC6");
0083
0084 _trigConfig.push_back(_trigRPC1);
0085 _trigConfig.push_back(_trigRPC2);
0086 _trigConfig.push_back(_trigRPC3);
0087 _trigConfig.push_back(_trigRPC4);
0088 _trigConfig.push_back(_trigRPC5);
0089 _trigConfig.push_back(_trigRPC6);
0090 }
0091
0092 void RPCRecHitReader::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0093 edm::ESHandle<RPCGeometry> rpcGeo;
0094 iSetup.get<MuonGeometryRecord>().get(rpcGeo);
0095
0096 RPCDetId id(region, wheel, station, sector, layer, subsector, 3);
0097 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeo->roll(id));
0098 _rollEff = roll;
0099
0100 fOutputFile = new TFile(fOutputFileName.c_str(), "RECREATE");
0101 fout = new std::fstream("RecHitOut.dat", std::ios::out);
0102
0103 _mapLayer[0] = -413.675;
0104 _mapLayer[1] = -448.675;
0105 _mapLayer[2] = -494.975;
0106 _mapLayer[3] = -529.975;
0107 _mapLayer[4] = -602.150;
0108 _mapLayer[5] = -704.550;
0109
0110 unsigned int layer = 0;
0111
0112 for (unsigned int i = 0; i < _trigConfig.size(); ++i) {
0113 if (_trigConfig[i] == false)
0114 layer = i;
0115 }
0116
0117 yLayer = _mapLayer[layer];
0118
0119 if (sector != 10) {
0120 GlobalPoint cntr10, cntr11;
0121 for (RPCGeometry::DetContainer::const_iterator it = rpcGeo->dets().begin(); it < rpcGeo->dets().end(); it++) {
0122 RPCRoll const* ir = dynamic_cast<const RPCRoll*>(*it);
0123 RPCDetId id = ir->id();
0124
0125 const Surface& bSurface = ir->surface();
0126
0127 if (id.region() == region && id.ring() == wheel && id.sector() == 10 && id.station() == 1 && id.layer() == 1) {
0128 LocalPoint orgn(0, 0, 0);
0129 cntr10 = bSurface.toGlobal(orgn);
0130 }
0131 if (id.region() == region && id.ring() == wheel && id.sector() == sector && id.station() == 1 &&
0132 id.layer() == 1) {
0133 LocalPoint orng(0, 0, 0);
0134 cntr11 = bSurface.toGlobal(orng);
0135 }
0136 }
0137
0138 float radius = 413.675;
0139 float crd = sqrt(std::pow((cntr10.x() - cntr11.x()), 2) + std::pow((cntr10.y() - cntr11.y()), 2));
0140 _phi = 2 * asin(crd / (2 * radius));
0141 }
0142
0143 *fout << "Angolo di rotazione = " << _phi << std::endl;
0144
0145 _trigger = 0;
0146 _triggerGOOD = 0;
0147 _spurious = 0;
0148 _spuriousPeak = 0;
0149 _efficiencyBAD = 0;
0150 _efficiencyGOOD = 0;
0151 _efficiencyBEST = 0;
0152
0153 histoXY = new TH2F("HistoXY", "Histoxy", 300, -300, 300, 300, -900, -300);
0154 histoSlope = new TH1F("Slope", "Slope", 100, -100, 100);
0155 histoChi2 = new TH1F("Chi2", "Chi2", 100, 0, 10000);
0156 histoRes = new TH1F("Residui", "Residui", 500, -100, 100);
0157 histoRes1 = new TH1F("Single Cluster Residuals", "Single Cluster Residuals", 500, -100, 100);
0158 histoRes2 = new TH1F("Double Cluster Residuals", "Double Cluster Residuals", 500, -100, 100);
0159 histoRes3 = new TH1F("Triple Cluster Residuals", "Triple Cluster Residuals", 500, -100, 100);
0160 histoPool1 = new TH1F("Single Cluster Size Pools", "Single Cluster Size Pools", 500, -100, 100);
0161 histoPool2 = new TH1F("Double Cluster Size Pools", "Double Cluster Size Pools", 500, -100, 100);
0162 histoPool3 = new TH1F("Triple Cluster Size Pools", "Triple Cluster Size Pools", 500, -100, 100);
0163
0164 histoExpectedOcc = new TH1F("ExpectedOcc", "ExpectedOcc", 100, -0.5, 99.5);
0165 histoRealOcc = new TH1F("RealOcc", "RealOcc", 100, -0.5, 99.5);
0166 histoLocalEff = new TH1F("LocalEff", "LocalEff", 100, -0.5, 99.5);
0167
0168 return;
0169 }
0170
0171
0172 void RPCRecHitReader::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0173 if (event.id().event() % 100 == 0)
0174 std::cout << " Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event() << std::endl;
0175
0176
0177 edm::ESHandle<RPCGeometry> rpcGeom;
0178 eventSetup.get<MuonGeometryRecord>().get(rpcGeom);
0179
0180
0181 edm::Handle<RPCRecHitCollection> recHits;
0182 event.getByLabel(recHitLabel1, recHitLabel2, recHits);
0183
0184
0185
0186
0187
0188 RPCRecHitCollection::const_iterator recIt;
0189
0190 std::vector<float> globalX, globalY, globalZ;
0191 std::vector<float> effX, effY, effZ;
0192 std::vector<float> errX, clus, rescut;
0193
0194 std::map<int, bool> _mapTrig;
0195
0196 TF1* func = new TF1("linearF", linearF, -300, 300, 2);
0197 func->SetParameters(0., 0.);
0198 func->SetParNames("angCoef", "interc");
0199
0200 for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
0201
0202 RPCDetId id = (RPCDetId)(*recIt).rpcId();
0203 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(id));
0204 const Surface& bSurface = roll->surface();
0205
0206 if ((roll->isForward()))
0207 return;
0208
0209 LocalPoint rhitlocal = (*recIt).localPosition();
0210 LocalError locerr = (*recIt).localPositionError();
0211 GlobalPoint rhitglob = bSurface.toGlobal(rhitlocal);
0212
0213 float x = 0, y = 0, z = 0;
0214
0215 if (id.sector() > 10 || (1 <= id.sector() && id.sector() <= 4)) {
0216 x = rhitglob.x() * cos(-_phi) - rhitglob.y() * sin(-_phi);
0217 y = rhitglob.y() * cos(-_phi) + rhitglob.x() * sin(-_phi);
0218 z = rhitglob.z();
0219 } else if (5 <= id.sector() && id.sector() <= 9) {
0220 x = rhitglob.x() * cos(_phi) - rhitglob.y() * sin(_phi);
0221 y = rhitglob.y() * cos(_phi) + rhitglob.x() * sin(_phi);
0222 z = rhitglob.z();
0223 } else if (id.sector() == 10) {
0224 x = rhitglob.x();
0225 y = rhitglob.y();
0226 z = rhitglob.z();
0227 }
0228
0229 if (_trigConfig[layerRecHit(*recIt)] == true && id.sector() == sector) {
0230 _mapTrig[layerRecHit(*recIt)] = true;
0231
0232 globalX.push_back(x);
0233 globalY.push_back(y);
0234 globalZ.push_back(z);
0235
0236 } else if (_trigConfig[layerRecHit(*recIt)] == false && id.sector() == sector) {
0237 effX.push_back(rhitglob.x());
0238 effY.push_back(rhitglob.y());
0239 effZ.push_back(rhitglob.z());
0240 errX.push_back(locerr.xx());
0241 clus.push_back((*recIt).clusterSize());
0242
0243 } else {
0244 continue;
0245 }
0246 }
0247
0248 if (_mapTrig.size() == 0)
0249 return;
0250
0251 char folder[128];
0252 sprintf(folder, "HistoXYFit_%llu", event.id().event());
0253 TH1F* histoXYFit = new TH1F(folder, folder, 300, -300, 300);
0254
0255 for (unsigned int i = 0; i < globalX.size(); ++i) {
0256 histoXY->Fill(globalX[i], globalY[i]);
0257 histoXYFit->Fill(globalX[i], globalY[i]);
0258 }
0259
0260
0261
0262 if (_mapTrig.size() > 4) {
0263 _trigger++;
0264 if (effX.size() > 0)
0265 _efficiencyBAD++;
0266
0267
0268
0269
0270 histoXYFit->Fit("linearF", "r");
0271 histoSlope->Fill(func->GetParameter(0));
0272 histoChi2->Fill(func->GetChisquare());
0273
0274 if (func->GetChisquare() < 0.5) {
0275 _triggerGOOD++;
0276
0277 float prepoint = ((yLayer)-func->GetParameter(1)) / func->GetParameter(0);
0278
0279 LocalPoint expPoint(prepoint, 0, 0);
0280 float expStrip = _rollEff->strip(expPoint);
0281
0282 histoExpectedOcc->Fill(expStrip);
0283
0284 if (effX.size() > 0) {
0285 _efficiencyGOOD++;
0286
0287 unsigned int k = 0;
0288 for (unsigned int i = 0; i < effX.size(); ++i) {
0289 float res = (effX[i] - prepoint);
0290 float errcl = errX[i] * clus[i];
0291 float pools = res / errX[i];
0292
0293 histoRes->Fill(res);
0294 if (clus[i] == 1) {
0295 histoRes1->Fill(res);
0296 histoPool1->Fill(pools);
0297 } else if (clus[i] == 2) {
0298 histoRes2->Fill(res);
0299 histoPool2->Fill(pools);
0300 } else if (clus[i] == 3) {
0301 histoRes3->Fill(res);
0302 histoPool3->Fill(pools);
0303 }
0304
0305 if (fabs(res) > errcl) {
0306 _spurious++;
0307 *fout << "Il residuo è maggiore di " << errcl << " ";
0308 *fout << " #Event: " << event.id().event() << std::endl;
0309 } else if (fabs(res) <= errcl) {
0310 k++;
0311 rescut.push_back(res);
0312 if (k == 1)
0313 histoRealOcc->Fill(expStrip);
0314 }
0315 }
0316
0317 if (rescut.size() > 1) {
0318 _spuriousPeak += rescut.size() - 1;
0319 *fout << "Ci sono più RecHit = " << rescut.size() << " ";
0320 *fout << " #Event: " << event.id().event() << std::endl;
0321 }
0322 } else {
0323 *fout << "Camera non efficiente!"
0324 << " "
0325 << " #Event: " << event.id().event() << std::endl;
0326 }
0327 }
0328 } else {
0329 *fout << "Non esiste RecHit appartenente a piani di interesse!"
0330 << " ";
0331 *fout << " #Event: " << event.id().event() << std::endl;
0332 }
0333 }
0334
0335 unsigned int RPCRecHitReader::layerRecHit(RPCRecHit rechit) {
0336 unsigned int layer = 0;
0337 RPCDetId id = (RPCDetId)(rechit).rpcId();
0338
0339 if (id.station() == 1 && id.layer() == 1)
0340 layer = 0;
0341 if (id.station() == 1 && id.layer() == 2)
0342 layer = 1;
0343 if (id.station() == 2 && id.layer() == 1)
0344 layer = 2;
0345 if (id.station() == 2 && id.layer() == 2)
0346 layer = 3;
0347 if (id.station() == 3 && id.layer() == 1)
0348 layer = 4;
0349 if (id.station() == 4 && id.layer() == 1)
0350 layer = 5;
0351
0352 return layer;
0353 }
0354
0355 void RPCRecHitReader::endJob() {
0356 TF1* bell = new TF1("Gaussiana", "gaus", -50, 50);
0357 histoRes->Fit("Gaussiana", "r");
0358 float sgmp = bell->GetParameter(1) + 3 * (bell->GetParameter(2));
0359 float sgmm = bell->GetParameter(1) - 3 * (bell->GetParameter(2));
0360
0361 int binmax = histoRes->GetXaxis()->FindBin(sgmp);
0362 int binmin = histoRes->GetXaxis()->FindBin(sgmm);
0363 _efficiencyBEST = histoRes->Integral(binmin, binmax);
0364
0365 _efficiencyBEST -= _spuriousPeak;
0366
0367 *fout << "Media = " << bell->GetParameter(1) << " Deviazione standard = " << bell->GetParameter(2) << std::endl;
0368 *fout << "Taglio a 3 sigma" << std::endl;
0369
0370 for (unsigned int i = 1; i <= 100; ++i) {
0371 if (histoExpectedOcc->GetBinContent(i) != 0) {
0372 float eff = histoRealOcc->GetBinContent(i) / histoExpectedOcc->GetBinContent(i);
0373 float erreff = sqrt(eff * (1 - eff) / histoExpectedOcc->GetBinContent(i));
0374 histoLocalEff->SetBinContent(i, eff);
0375 histoLocalEff->SetBinError(i, erreff);
0376 }
0377 }
0378
0379 histoRes1->Fit("Gaussiana", "r");
0380 histoRes2->Fit("Gaussiana", "r");
0381 histoRes3->Fit("Gaussiana", "r");
0382
0383 histoPool1->Fit("Gaussiana", "r");
0384 histoPool2->Fit("Gaussiana", "r");
0385 histoPool3->Fit("Gaussiana", "r");
0386
0387 fOutputFile->Write();
0388 fOutputFile->Close();
0389 return;
0390 }
0391
0392
0393 RPCRecHitReader::~RPCRecHitReader() {
0394 *fout << "Trigger Complessivi = " << _trigger << std::endl;
0395 *fout << "Trigger GOOD = " << _triggerGOOD << std::endl;
0396 *fout << "Efficiency BAD Counter = " << _efficiencyBAD << std::endl;
0397 *fout << "Efficiency GOOD Counter = " << _efficiencyGOOD << std::endl;
0398 *fout << "Efficiency BEST Counter = " << _efficiencyBEST << std::endl;
0399 *fout << "Spurious counter = " << _spurious << std::endl;
0400
0401 *fout << "Efficienza BAD = " << _efficiencyBAD / _trigger << std::endl;
0402 *fout << "Efficienza GOOD = " << _efficiencyGOOD / _triggerGOOD << std::endl;
0403 *fout << "Efficienza BEST = " << _efficiencyBEST / _triggerGOOD << std::endl;
0404 }
0405
0406 DEFINE_FWK_MODULE(RPCRecHitReader);
