DTSegment/test/DTClusAnalyzer.cc

0001 /******* \class DTClusAnalyzer *******
0002  *
0003  * Description:
0004  *
0005  *  detailed description
0006  *
0007  * \author : Stefano Lacaprara - INFN LNL <stefano.lacaprara@pd.infn.it>
0008  *
0009  * Modification:
0010  *
0011  *********************************/
0012
0013 /* This Class Header */
0014 #include "RecoLocalMuon/DTSegment/test/DTClusAnalyzer.h"
0015
0016 /* Collaborating Class Header */
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023
0024 using namespace edm;
0025
0026 #include "TFile.h"
0027 #include "TH1F.h"
0028 #include "TH2F.h"
0029
0030 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0031 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0032
0033 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0034 #include "DataFormats/DTRecHit/interface/DTRecClusterCollection.h"
0035 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0036 /* C++ Headers */
0037 #include <iostream>
0038 #include <cmath>
0039 using namespace std;
0040
0041 /* ====================================================================== */
0042
0043 /* Constructor */
0044 DTClusAnalyzer::DTClusAnalyzer(const ParameterSet& pset) : _ev(0) {
0045   theDtGeomToken = esConsumes();
0046   // Get the debug parameter for verbose output
0047   debug = pset.getUntrackedParameter<bool>("debug");
0048   theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0049
0050   // the name of the clus rec hits collection
0051   theRecClusLabel = pset.getParameter<string>("recClusLabel");
0052   theRecClusToken = consumes(edm::InputTag(theRecClusLabel));
0053
0054   // the name of the 1D rec hits collection
0055   theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0056   theRecHits1DToken = consumes(edm::InputTag(theRecHits1DLabel));
0057
0058   // the name of the 2D rec hits collection
0059   theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0060   theRecHits2DToken = consumes(edm::InputTag(theRecHits2DLabel));
0061
0062   // Create the root file
0063   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0064   bool dirStat = TH1::AddDirectoryStatus();
0065   TH1::AddDirectory(kTRUE);
0066
0067   /// DT histos
0068   new TH1F("hnClus", "Num 1d clus ", 50, 0., 50.);
0069   new TH1F("hClusPosX", "Local Pos Cluster :x(cm)", 100, -300., 300.);
0070   new TH1F("hClusRms", "Width clus ", 100, 0., 20.);
0071   new TH1F("hClusNHits", "# Hits clus ", 30, 0., 30.);
0072
0073   new TH1F("hnHit", "Num 1d Hits ", 100, 0., 100.);
0074
0075   new TH1F("hnSeg", "Num 1d seg ", 50, 0., 50.);
0076   new TH1F("hSegPosX", "Local Pos seg :x(cm)", 100, -300., 300.);
0077   new TH1F("hSegRms", "Width seg ", 100, 0., 1.);
0078   new TH1F("hSegNHits", "# Hits seg ", 10, 0., 10.);
0079
0080   new TH2F("hnClusVsSegs", "# clus vs # segs", 30, 0, 30, 30, 0, 30);
0081   new TH2F("hnClusVsHits", "# clus vs # hits", 100, 0, 100, 30, 0, 30);
0082
0083   // per SL
0084   new TH1F("hnHitSL", "Num 1d hits per SL", 100, 0., 100.);
0085   new TH1F("hnClusSL", "Num 1d clus per SL", 10, 0., 10.);
0086   new TH1F("hnSegSL", "Num 2d seg per SL", 10, 0., 10.);
0087
0088   new TH2F("hnClusVsSegsSL", "# clus vs # segs per SL", 30, 0, 30, 30, 0, 30);
0089   new TH2F("hnClusVsHitsSL", "# clus vs # hits per SL", 100, 0, 100, 30, 0, 30);
0090
0091   new TH1F("hClusSegDistSL", "#Delta x (clus segs) per SL", 100, -20, 20);
0092   new TH2F("hClusVsSegPosSL", "X (clus vs segs) per SL", 100, -300, 300, 100, -300, 300);
0093
0094   new TH2F("hClusVsSegHitSL", "#hits (clus vs segs) per SL", 20, 0, 20, 20, 0, 20);
0095   TH1::AddDirectory(dirStat);
0096 }
0097
0098 /* Destructor */
0099 DTClusAnalyzer::~DTClusAnalyzer() {
0100   theFile->cd();
0101   theFile->Write();
0102   theFile->Close();
0103 }
0104
0105 /* Operations */
0106 void DTClusAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0107   _ev++;
0108
0109   static int j = 1;
0110   if ((_ev % j) == 0) {
0111     if ((_ev / j) == 9)
0112       j *= 10;
0113     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev << endl;
0114   }
0115
0116   // Get the DT Geometry
0117   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDtGeomToken);
0118
0119   // Get the 1D clusters from the event --------------
0120   Handle<DTRecClusterCollection> dtClusters = event.getHandle(theRecClusToken);
0121
0122   // Get the 1D rechits from the event --------------
0123   Handle<DTRecHitCollection> dtRecHits = event.getHandle(theRecHits1DToken);
0124
0125   // Get the 2D rechit collection from the event -------------------
0126   edm::Handle<DTRecSegment2DCollection> segs2d = event.getHandle(theRecHits2DToken);
0127
0128   // only clusters
0129   int nClus = dtClusters->size();
0130   histo("hnClus")->Fill(nClus);
0131   for (DTRecClusterCollection::const_iterator clus = dtClusters->begin(); clus != dtClusters->end(); ++clus) {
0132     histo("hClusPosX")->Fill((*clus).localPosition().x());
0133     histo("hClusRms")->Fill(sqrt((*clus).localPositionError().xx()));
0134     histo("hClusNHits")->Fill((*clus).nHits());
0135   }
0136
0137   // only segments
0138   int nSeg = segs2d->size();
0139   histo("hnSeg")->Fill(nSeg);
0140   for (DTRecSegment2DCollection::const_iterator seg = segs2d->begin(); seg != segs2d->end(); ++seg) {
0141     histo("hSegPosX")->Fill((*seg).localPosition().x());
0142     histo("hSegRms")->Fill(sqrt((*seg).localPositionError().xx()));
0143     histo("hSegNHits")->Fill((*seg).recHits().size());
0144   }
0145
0146   int nHit = dtRecHits->size();
0147   histo("hnHit")->Fill(nHit);
0148
0149   // clus vs segs, hits and all this
0150   histo2d("hnClusVsSegs")->Fill(nSeg, nClus);
0151   histo2d("hnClusVsHits")->Fill(nHit, nClus);
0152
0153   // loop over SL and get hits, clus and segs2d for each
0154   const std::vector<const DTSuperLayer*>& sls = dtGeom->superLayers();
0155   for (auto sl = sls.begin(); sl != sls.end(); ++sl) {
0156     DTSuperLayerId slid((*sl)->id());
0157
0158     // Hits 1d in this sl
0159     // SL: apparently I have to loop over the 4 layers!! What the Fuck!
0160     int nHitSL = 0;
0161     for (int il = 1; il <= 4; ++il) {
0162       DTRecHitCollection::range hitSL = dtRecHits->get(DTLayerId(slid, il));
0163       nHitSL += hitSL.second - hitSL.first;
0164     }
0165
0166     // Cluster 1d in this sl
0167     DTRecClusterCollection::range clusSL = dtClusters->get(slid);
0168     int nClusSL = clusSL.second - clusSL.first;
0169
0170     // Segment 2d in this sl
0171     DTRecSegment2DCollection::range segSL = segs2d->get(slid);
0172     int nSegSL = segSL.second - segSL.first;
0173
0174     histo("hnHitSL")->Fill(nHitSL);
0175     histo("hnClusSL")->Fill(nClusSL);
0176     histo("hnSegSL")->Fill(nSegSL);
0177     // clus vs segs, hits and all this
0178     histo2d("hnClusVsSegsSL")->Fill(nSegSL, nClusSL);
0179     histo2d("hnClusVsHitsSL")->Fill(nHitSL, nClusSL);
0180
0181     for (DTRecClusterCollection::const_iterator clus = clusSL.first; clus != clusSL.second; ++clus) {
0182       float minDist = 99999.;
0183       LocalPoint clusPos = (*clus).localPosition();
0184       DTRecSegment2DCollection::const_iterator closestSeg = segSL.second;
0185       for (DTRecSegment2DCollection::const_iterator seg = segSL.first; seg != segSL.second; ++seg) {
0186         LocalPoint segPos = (*seg).localPosition();
0187         float dist = (clusPos - segPos).mag2();
0188         if (dist < minDist) {
0189           minDist = dist;
0190           closestSeg = seg;
0191         }
0192       }
0193       if (closestSeg != segSL.second) {
0194         histo("hClusSegDistSL")->Fill((clusPos - (*closestSeg).localPosition()).x());
0195         histo2d("hClusVsSegPosSL")->Fill((*closestSeg).localPosition().x(), clusPos.x());
0196         histo2d("hClusVsSegHitSL")->Fill((*closestSeg).recHits().size(), (*clus).nHits());
0197       }
0198     }
0199   }
0200 }
0201
0202 TH1F* DTClusAnalyzer::histo(const string& name) const {
0203   if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0204     return h;
0205   else
0206     throw cms::Exception("DTSegAnalyzer") << " Not a TH1F " << name;
0207 }
0208
0209 TH2F* DTClusAnalyzer::histo2d(const string& name) const {
0210   if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0211     return h;
0212   else
0213     throw cms::Exception("DTSegAnalyzer") << " Not a TH2F " << name;
0214 }
0215
0216 DEFINE_FWK_MODULE(DTClusAnalyzer);