RPCRecHit/test/RPCRecHitReader.cc

0001 /*
0002  *  Program to calculate the global and local efficiency of
0003  *  RPC chambers using a software autotrigger and
0004  *  a linear reconstruction of tracks.
0005  *
0006  *  Author: R. Trentadue - Bari University
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 // Constructor
0061 RPCRecHitReader::RPCRecHitReader(const edm::ParameterSet& pset) : _phi(0) {
0062   rpcGeomToken_ = esConsumes();
0063   rpcGeomBRToken_ = esConsumes<edm::Transition::BeginRun>();
0064   // Get the various input parameters
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 // The Analysis  (the main)
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   // Get the RPC Geometry :
0177   edm::ESHandle<RPCGeometry> rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0178
0179   // Get the RecHits collection :
0180   edm::Handle<RPCRecHitCollection> recHits;
0181   event.getByLabel(recHitLabel1, recHitLabel2, recHits);
0182
0183   //----------------------------------------------------------------------
0184   //---------------  Loop over rechits -----------------------------------
0185
0186   // Build iterator for rechits and loop :
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     // Find chamber with rechits in RPC
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   //-------- PRIMA STIMA EFFICIENZA ----------------------
0260
0261   if (_mapTrig.size() > 4) {
0262     _trigger++;
0263     if (effX.size() > 0)
0264       _efficiencyBAD++;
0265
0266     //-----------------------------------------------------
0267     //--------------- FIT SULLE TRACCE---------------------
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 // Destructor
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);