File indexing completed on 2024-04-06 12:26:15
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/MessageLogger/interface/MessageLogger.h"
0024 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0025 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0026 #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h"
0027 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0028
0029 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0030 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0031 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0032 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0033 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0034 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0035 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0036
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/Frameworkfwd.h"
0039 #include "FWCore/Framework/interface/Event.h"
0040 #include "FWCore/Framework/interface/EventSetup.h"
0041 #include "FWCore/Framework/interface/ESHandle.h"
0042
0043 #include "TFile.h"
0044 #include "TVector3.h"
0045 #include "TH1F.h"
0046 #include "TH2F.h"
0047 #include "TF1.h"
0048 #include "TF2.h"
0049 #include "TVectorT.h"
0050 #include "TGraph.h"
0051
0052 #include <map>
0053
0054 Double_t linearF(Double_t* x, Double_t* par) {
0055 Double_t y = 0.;
0056 y = par[0] * (*x) + par[1];
0057 return y;
0058 }
0059
0060
0061 RPCRecHitReader::RPCRecHitReader(const edm::ParameterSet& pset) : _phi(0) {
0062 rpcGeomToken_ = esConsumes();
0063 rpcGeomBRToken_ = esConsumes<edm::Transition::BeginRun>();
0064
0065 fOutputFileName = pset.getUntrackedParameter<std::string>("HistOutFile");
0066 recHitLabel1 = pset.getUntrackedParameter<std::string>("recHitLabel1");
0067 recHitLabel2 = pset.getUntrackedParameter<std::string>("recHitLabel2");
0068
0069 region = pset.getUntrackedParameter<int>("region", 0);
0070 wheel = pset.getUntrackedParameter<int>("wheel", 1);
0071 sector = pset.getUntrackedParameter<int>("sector", 10);
0072 station = pset.getUntrackedParameter<int>("station");
0073 layer = pset.getUntrackedParameter<int>("layer");
0074 subsector = pset.getUntrackedParameter<int>("subsector");
0075
0076 _trigConfig.clear();
0077
0078 _trigRPC1 = pset.getUntrackedParameter<bool>("trigRPC1");
0079 _trigRPC2 = pset.getUntrackedParameter<bool>("trigRPC2");
0080 _trigRPC3 = pset.getUntrackedParameter<bool>("trigRPC3");
0081 _trigRPC4 = pset.getUntrackedParameter<bool>("trigRPC4");
0082 _trigRPC5 = pset.getUntrackedParameter<bool>("trigRPC5");
0083 _trigRPC6 = pset.getUntrackedParameter<bool>("trigRPC6");
0084
0085 _trigConfig.push_back(_trigRPC1);
0086 _trigConfig.push_back(_trigRPC2);
0087 _trigConfig.push_back(_trigRPC3);
0088 _trigConfig.push_back(_trigRPC4);
0089 _trigConfig.push_back(_trigRPC5);
0090 _trigConfig.push_back(_trigRPC6);
0091 }
0092
0093 void RPCRecHitReader::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0094 edm::ESHandle<RPCGeometry> rpcGeo = iSetup.getHandle(rpcGeomBRToken_);
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 = eventSetup.getHandle(rpcGeomToken_);
0178
0179
0180 edm::Handle<RPCRecHitCollection> recHits;
0181 event.getByLabel(recHitLabel1, recHitLabel2, recHits);
0182
0183
0184
0185
0186
0187 RPCRecHitCollection::const_iterator recIt;
0188
0189 std::vector<float> globalX, globalY, globalZ;
0190 std::vector<float> effX, effY, effZ;
0191 std::vector<float> errX, clus, rescut;
0192
0193 std::map<int, bool> _mapTrig;
0194
0195 TF1* func = new TF1("linearF", linearF, -300, 300, 2);
0196 func->SetParameters(0., 0.);
0197 func->SetParNames("angCoef", "interc");
0198
0199 for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
0200
0201 RPCDetId id = (RPCDetId)(*recIt).rpcId();
0202 const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(id));
0203 const Surface& bSurface = roll->surface();
0204
0205 if ((roll->isForward()))
0206 return;
0207
0208 LocalPoint rhitlocal = (*recIt).localPosition();
0209 LocalError locerr = (*recIt).localPositionError();
0210 GlobalPoint rhitglob = bSurface.toGlobal(rhitlocal);
0211
0212 float x = 0, y = 0, z = 0;
0213
0214 if (id.sector() > 10 || (1 <= id.sector() && id.sector() <= 4)) {
0215 x = rhitglob.x() * cos(-_phi) - rhitglob.y() * sin(-_phi);
0216 y = rhitglob.y() * cos(-_phi) + rhitglob.x() * sin(-_phi);
0217 z = rhitglob.z();
0218 } else if (5 <= id.sector() && id.sector() <= 9) {
0219 x = rhitglob.x() * cos(_phi) - rhitglob.y() * sin(_phi);
0220 y = rhitglob.y() * cos(_phi) + rhitglob.x() * sin(_phi);
0221 z = rhitglob.z();
0222 } else if (id.sector() == 10) {
0223 x = rhitglob.x();
0224 y = rhitglob.y();
0225 z = rhitglob.z();
0226 }
0227
0228 if (_trigConfig[layerRecHit(*recIt)] == true && id.sector() == sector) {
0229 _mapTrig[layerRecHit(*recIt)] = true;
0230
0231 globalX.push_back(x);
0232 globalY.push_back(y);
0233 globalZ.push_back(z);
0234
0235 } else if (_trigConfig[layerRecHit(*recIt)] == false && id.sector() == sector) {
0236 effX.push_back(rhitglob.x());
0237 effY.push_back(rhitglob.y());
0238 effZ.push_back(rhitglob.z());
0239 errX.push_back(locerr.xx());
0240 clus.push_back((*recIt).clusterSize());
0241
0242 } else {
0243 continue;
0244 }
0245 }
0246
0247 if (_mapTrig.size() == 0)
0248 return;
0249
0250 char folder[128];
0251 sprintf(folder, "HistoXYFit_%llu", event.id().event());
0252 TH1F* histoXYFit = new TH1F(folder, folder, 300, -300, 300);
0253
0254 for (unsigned int i = 0; i < globalX.size(); ++i) {
0255 histoXY->Fill(globalX[i], globalY[i]);
0256 histoXYFit->Fill(globalX[i], globalY[i]);
0257 }
0258
0259
0260
0261 if (_mapTrig.size() > 4) {
0262 _trigger++;
0263 if (effX.size() > 0)
0264 _efficiencyBAD++;
0265
0266
0267
0268
0269 histoXYFit->Fit("linearF", "r");
0270 histoSlope->Fill(func->GetParameter(0));
0271 histoChi2->Fill(func->GetChisquare());
0272
0273 if (func->GetChisquare() < 0.5) {
0274 _triggerGOOD++;
0275
0276 float prepoint = ((yLayer)-func->GetParameter(1)) / func->GetParameter(0);
0277
0278 LocalPoint expPoint(prepoint, 0, 0);
0279 float expStrip = _rollEff->strip(expPoint);
0280
0281 histoExpectedOcc->Fill(expStrip);
0282
0283 if (effX.size() > 0) {
0284 _efficiencyGOOD++;
0285
0286 unsigned int k = 0;
0287 for (unsigned int i = 0; i < effX.size(); ++i) {
0288 float res = (effX[i] - prepoint);
0289 float errcl = errX[i] * clus[i];
0290 float pools = res / errX[i];
0291
0292 histoRes->Fill(res);
0293 if (clus[i] == 1) {
0294 histoRes1->Fill(res);
0295 histoPool1->Fill(pools);
0296 } else if (clus[i] == 2) {
0297 histoRes2->Fill(res);
0298 histoPool2->Fill(pools);
0299 } else if (clus[i] == 3) {
0300 histoRes3->Fill(res);
0301 histoPool3->Fill(pools);
0302 }
0303
0304 if (fabs(res) > errcl) {
0305 _spurious++;
0306 *fout << "Il residuo è maggiore di " << errcl << " ";
0307 *fout << " #Event: " << event.id().event() << std::endl;
0308 } else if (fabs(res) <= errcl) {
0309 k++;
0310 rescut.push_back(res);
0311 if (k == 1)
0312 histoRealOcc->Fill(expStrip);
0313 }
0314 }
0315
0316 if (rescut.size() > 1) {
0317 _spuriousPeak += rescut.size() - 1;
0318 *fout << "Ci sono più RecHit = " << rescut.size() << " ";
0319 *fout << " #Event: " << event.id().event() << std::endl;
0320 }
0321 } else {
0322 *fout << "Camera non efficiente!"
0323 << " "
0324 << " #Event: " << event.id().event() << std::endl;
0325 }
0326 }
0327 } else {
0328 *fout << "Non esiste RecHit appartenente a piani di interesse!"
0329 << " ";
0330 *fout << " #Event: " << event.id().event() << std::endl;
0331 }
0332 }
0333
0334 unsigned int RPCRecHitReader::layerRecHit(RPCRecHit rechit) {
0335 unsigned int layer = 0;
0336 RPCDetId id = (RPCDetId)(rechit).rpcId();
0337
0338 if (id.station() == 1 && id.layer() == 1)
0339 layer = 0;
0340 if (id.station() == 1 && id.layer() == 2)
0341 layer = 1;
0342 if (id.station() == 2 && id.layer() == 1)
0343 layer = 2;
0344 if (id.station() == 2 && id.layer() == 2)
0345 layer = 3;
0346 if (id.station() == 3 && id.layer() == 1)
0347 layer = 4;
0348 if (id.station() == 4 && id.layer() == 1)
0349 layer = 5;
0350
0351 return layer;
0352 }
0353
0354 void RPCRecHitReader::endJob() {
0355 TF1* bell = new TF1("Gaussiana", "gaus", -50, 50);
0356 histoRes->Fit("Gaussiana", "r");
0357 float sgmp = bell->GetParameter(1) + 3 * (bell->GetParameter(2));
0358 float sgmm = bell->GetParameter(1) - 3 * (bell->GetParameter(2));
0359
0360 int binmax = histoRes->GetXaxis()->FindBin(sgmp);
0361 int binmin = histoRes->GetXaxis()->FindBin(sgmm);
0362 _efficiencyBEST = histoRes->Integral(binmin, binmax);
0363
0364 _efficiencyBEST -= _spuriousPeak;
0365
0366 *fout << "Media = " << bell->GetParameter(1) << " Deviazione standard = " << bell->GetParameter(2) << std::endl;
0367 *fout << "Taglio a 3 sigma" << std::endl;
0368
0369 for (unsigned int i = 1; i <= 100; ++i) {
0370 if (histoExpectedOcc->GetBinContent(i) != 0) {
0371 float eff = histoRealOcc->GetBinContent(i) / histoExpectedOcc->GetBinContent(i);
0372 float erreff = sqrt(eff * (1 - eff) / histoExpectedOcc->GetBinContent(i));
0373 histoLocalEff->SetBinContent(i, eff);
0374 histoLocalEff->SetBinError(i, erreff);
0375 }
0376 }
0377
0378 histoRes1->Fit("Gaussiana", "r");
0379 histoRes2->Fit("Gaussiana", "r");
0380 histoRes3->Fit("Gaussiana", "r");
0381
0382 histoPool1->Fit("Gaussiana", "r");
0383 histoPool2->Fit("Gaussiana", "r");
0384 histoPool3->Fit("Gaussiana", "r");
0385
0386 fOutputFile->Write();
0387 fOutputFile->Close();
0388 return;
0389 }
0390
0391
0392 RPCRecHitReader::~RPCRecHitReader() {
0393 *fout << "Trigger Complessivi = " << _trigger << std::endl;
0394 *fout << "Trigger GOOD = " << _triggerGOOD << std::endl;
0395 *fout << "Efficiency BAD Counter = " << _efficiencyBAD << std::endl;
0396 *fout << "Efficiency GOOD Counter = " << _efficiencyGOOD << std::endl;
0397 *fout << "Efficiency BEST Counter = " << _efficiencyBEST << std::endl;
0398 *fout << "Spurious counter = " << _spurious << std::endl;
0399
0400 *fout << "Efficienza BAD = " << _efficiencyBAD / _trigger << std::endl;
0401 *fout << "Efficienza GOOD = " << _efficiencyGOOD / _triggerGOOD << std::endl;
0402 *fout << "Efficienza BEST = " << _efficiencyBEST / _triggerGOOD << std::endl;
0403 }
0404
0405 DEFINE_FWK_MODULE(RPCRecHitReader);
