Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:59

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/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 // Constructor
0062 RPCRecHitReader::RPCRecHitReader(const edm::ParameterSet& pset) : _phi(0) {
0063   // Get the various input parameters
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 // 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;
0178   eventSetup.get<MuonGeometryRecord>().get(rpcGeom);
0179 
0180   // Get the RecHits collection :
0181   edm::Handle<RPCRecHitCollection> recHits;
0182   event.getByLabel(recHitLabel1, recHitLabel2, recHits);
0183 
0184   //----------------------------------------------------------------------
0185   //---------------  Loop over rechits -----------------------------------
0186 
0187   // Build iterator for rechits and loop :
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     // Find chamber with rechits in RPC
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   //-------- PRIMA STIMA EFFICIENZA ----------------------
0261 
0262   if (_mapTrig.size() > 4) {
0263     _trigger++;
0264     if (effX.size() > 0)
0265       _efficiencyBAD++;
0266 
0267     //-----------------------------------------------------
0268     //--------------- FIT SULLE TRACCE---------------------
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 // Destructor
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);