DTRecHit/test/DTRecHitReader.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - INFN Torino
0005  */
0006
0007 #include "DTRecHitReader.h"
0008
0009 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0010
0011 #include "Geometry/DTGeometry/interface/DTLayer.h"
0012 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0013 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0014
0015 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0016
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021
0022 #include "TFile.h"
0023
0024 #include <iostream>
0025 #include <map>
0026
0027 using namespace std;
0028 using namespace edm;
0029
0030 // Constructor
0031 DTRecHitReader::DTRecHitReader(const ParameterSet& pset) {
0032   // Get the debug parameter for verbose output
0033   debug = pset.getUntrackedParameter<bool>("debug");
0034   rootFileName = pset.getUntrackedParameter<string>("rootFileName");
0035   simHitLabel = pset.getUntrackedParameter<string>("simHitLabel", "r");
0036   recHitLabel = pset.getUntrackedParameter<string>("recHitLabel", "rechitbuilder");
0037
0038   if (debug)
0039     cout << "[DTRecHitReader] Constructor called" << endl;
0040
0041   // Create the root file
0042   theFile = new TFile(rootFileName.c_str(), "RECREATE");
0043   theFile->cd();
0044
0045   // Book the histograms
0046   hRHitPhi = new H1DRecHit("RPhi");
0047   hRHitZ_W0 = new H1DRecHit("RZ_W0");
0048   hRHitZ_W1 = new H1DRecHit("RZ_W1");
0049   hRHitZ_W2 = new H1DRecHit("RZ_W2");
0050   hRHitZ_All = new H1DRecHit("RZ_All");
0051
0052   dtGeomToken_ = esConsumes();
0053 }
0054
0055 // Destructor
0056 DTRecHitReader::~DTRecHitReader() {
0057   if (debug)
0058     cout << "[DTRecHitReader] Destructor called" << endl;
0059
0060   // Write the histos to file
0061   theFile->cd();
0062   hRHitPhi->Write();
0063   hRHitZ_W0->Write();
0064   hRHitZ_W1->Write();
0065   hRHitZ_W2->Write();
0066   hRHitZ_All->Write();
0067   theFile->Close();
0068 }
0069
0070 // The real analysis
0071 void DTRecHitReader::analyze(const Event& event, const EventSetup& eventSetup) {
0072   cout << "--- [DTRecHitReader] Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;
0073   // Get the DT Geometry
0074   ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(dtGeomToken_);
0075   ;
0076
0077   // Get the rechit collection from the event
0078   Handle<DTRecHitCollection> dtRecHits;
0079   event.getByLabel(recHitLabel, dtRecHits);
0080
0081   // Get the SimHit collection from the event
0082   Handle<PSimHitContainer> simHits;
0083
0084   event.getByLabel(simHitLabel, "MuonDTHits", simHits);
0085
0086   if (debug)
0087     cout << "   #SimHits: " << simHits->size() << endl;
0088
0089   // Map simhits per wire
0090   map<DTWireId, vector<const PSimHit*> > simHitMap = mapSimHitsPerWire(simHits);
0091
0092   // Iterate over all detunits
0093   DTRecHitCollection::id_iterator detUnitIt;
0094   for (detUnitIt = dtRecHits->id_begin(); detUnitIt != dtRecHits->id_end(); ++detUnitIt) {
0095     // Get the GeomDet from the setup
0096     const DTLayer* layer = dtGeom->layer(*detUnitIt);
0097
0098     // Get the range for the corresponding LayerId
0099     DTRecHitCollection::range range = dtRecHits->get((*detUnitIt));
0100     // Loop over the rechits of this DetUnit
0101     for (DTRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
0102       // Get the wireId of the rechit
0103       DTWireId wireId = (*rechit).wireId();
0104
0105       float xwire = layer->specificTopology().wirePosition(wireId.wire());
0106       if (fabs(xwire - (*rechit).localPosition().x()) > 0.00001) {
0107         cout << "  [DTRecHitReader]***Error in wire Position: xwire = " << xwire
0108              << " xRecHit = " << (*rechit).localPosition().x() << endl;
0109       }
0110
0111       // Access to Right and left rechits
0112       pair<const DTRecHit1D*, const DTRecHit1D*> lrRecHits = (*rechit).componentRecHits();
0113
0114       cout << "Left Hit x(cm): " << lrRecHits.first->localPosition().x() << endl;
0115       cout << "Right Hit x(cm): " << lrRecHits.second->localPosition().x() << endl;
0116
0117       // Compute the rechit distance from wire
0118       float distFromWire =
0119           fabs((*rechit).localPosition(DTEnums::Left).x() - (*rechit).localPosition(DTEnums::Right).x()) / 2.;
0120
0121       // Search the best mu simhit and compute its distance from the wire
0122       float simHitDistFromWire = 0;
0123       if (simHitMap.find(wireId) != simHitMap.end()) {
0124         const PSimHit* muSimHit = findBestMuSimHit(layer, wireId, simHitMap[wireId], distFromWire);
0125         // Check that a mu simhit is found
0126         if (muSimHit != 0) {
0127           // Compute the simhit distance from wire
0128           simHitDistFromWire = findSimHitDist(layer, wireId, muSimHit);
0129           // Fill the histos
0130           H1DRecHit* histo = 0;
0131           if (wireId.superlayer() == 1 || wireId.superlayer() == 3) {
0132             histo = hRHitPhi;
0133           } else if (wireId.superlayer() == 2) {
0134             hRHitZ_All->Fill(distFromWire, simHitDistFromWire);
0135             if (wireId.wheel() == 0) {
0136               histo = hRHitZ_W0;
0137             } else if (abs(wireId.wheel()) == 1) {
0138               histo = hRHitZ_W1;
0139             } else if (abs(wireId.wheel()) == 2) {
0140               histo = hRHitZ_W2;
0141             }
0142           }
0143           histo->Fill(distFromWire, simHitDistFromWire);
0144
0145           if (fabs(distFromWire - simHitDistFromWire) > 2.1) {
0146             cout << "Warning: " << endl
0147                  << "  RecHit distance from wire is: " << distFromWire << endl
0148                  << "  SimHit distance from wire is: " << simHitDistFromWire << endl
0149                  << "  RecHit wire Id is: " << wireId << endl
0150                  << "  SimHit wire Id is: " << DTWireId((*muSimHit).detUnitId()) << endl;
0151             cout << "  Wire x: = " << xwire << endl
0152                  << "  RecHit x = " << (*rechit).localPosition().x() << endl
0153                  << "  SimHit x = " << (*muSimHit).localPosition().x() << endl;
0154           }
0155
0156           // Some printout
0157           if (debug) {
0158             cout << "[DTRecHitReader]: " << endl
0159                  << "         WireId: " << wireId << endl
0160                  << "         1DRecHitPair local position (cm): " << (*rechit) << endl
0161                  << "         RecHit distance from wire (cm): " << distFromWire << endl
0162                  << "         Mu SimHit distance from wire (cm): " << simHitDistFromWire << endl;
0163           }
0164         }
0165       }
0166     }
0167   }
0168 }
0169
0170 // Return a map between simhits of a layer and the wireId of their cell
0171 map<DTWireId, vector<const PSimHit*> > DTRecHitReader::mapSimHitsPerWire(const Handle<PSimHitContainer>& simhits) {
0172   map<DTWireId, vector<const PSimHit*> > hitWireMapResult;
0173
0174   for (PSimHitContainer::const_iterator simhit = simhits->begin(); simhit != simhits->end(); simhit++) {
0175     hitWireMapResult[DTWireId((*simhit).detUnitId())].push_back(&(*simhit));
0176   }
0177
0178   return hitWireMapResult;
0179 }
0180
0181 const PSimHit* DTRecHitReader::findBestMuSimHit(const DTLayer* layer,
0182                                                 const DTWireId& wireId,
0183                                                 const vector<const PSimHit*>& simhits,
0184                                                 float recHitDistFromWire) {
0185   const PSimHit* retSimHit = 0;
0186   float tmp_distDiff = 999999;
0187   for (vector<const PSimHit*>::const_iterator simhit = simhits.begin(); simhit != simhits.end(); simhit++) {
0188     // Select muons
0189     if (abs((*simhit)->particleType()) == 13) {
0190       // Get the mu simhit closest to the rechit
0191       if (findSimHitDist(layer, wireId, *simhit) - recHitDistFromWire < tmp_distDiff) {
0192         tmp_distDiff = findSimHitDist(layer, wireId, *simhit) - recHitDistFromWire;
0193         retSimHit = (*simhit);
0194       }
0195     }
0196   }
0197   return retSimHit;
0198 }
0199
0200 // Compute SimHit distance from wire
0201 double DTRecHitReader::findSimHitDist(const DTLayer* layer, const DTWireId& wireId, const PSimHit* hit) {
0202   float xwire = layer->specificTopology().wirePosition(wireId.wire());
0203   LocalPoint entryP = hit->entryPoint();
0204   LocalPoint exitP = hit->exitPoint();
0205   float xEntry = entryP.x() - xwire;
0206   float xExit = exitP.x() - xwire;
0207
0208   return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));  //FIXME: check...
0209 }
0210
0211 DEFINE_FWK_MODULE(DTRecHitReader);