Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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 
0053 // Destructor
0054 DTRecHitReader::~DTRecHitReader() {
0055   if (debug)
0056     cout << "[DTRecHitReader] Destructor called" << endl;
0057 
0058   // Write the histos to file
0059   theFile->cd();
0060   hRHitPhi->Write();
0061   hRHitZ_W0->Write();
0062   hRHitZ_W1->Write();
0063   hRHitZ_W2->Write();
0064   hRHitZ_All->Write();
0065   theFile->Close();
0066 }
0067 
0068 // The real analysis
0069 void DTRecHitReader::analyze(const Event& event, const EventSetup& eventSetup) {
0070   cout << "--- [DTRecHitReader] Event analysed #Run: " << event.id().run() << " #Event: " << event.id().event() << endl;
0071   // Get the DT Geometry
0072   ESHandle<DTGeometry> dtGeom;
0073   eventSetup.get<MuonGeometryRecord>().get(dtGeom);
0074 
0075   // Get the rechit collection from the event
0076   Handle<DTRecHitCollection> dtRecHits;
0077   event.getByLabel(recHitLabel, dtRecHits);
0078 
0079   // Get the SimHit collection from the event
0080   Handle<PSimHitContainer> simHits;
0081 
0082   event.getByLabel(simHitLabel, "MuonDTHits", simHits);
0083 
0084   if (debug)
0085     cout << "   #SimHits: " << simHits->size() << endl;
0086 
0087   // Map simhits per wire
0088   map<DTWireId, vector<const PSimHit*> > simHitMap = mapSimHitsPerWire(simHits);
0089 
0090   // Iterate over all detunits
0091   DTRecHitCollection::id_iterator detUnitIt;
0092   for (detUnitIt = dtRecHits->id_begin(); detUnitIt != dtRecHits->id_end(); ++detUnitIt) {
0093     // Get the GeomDet from the setup
0094     const DTLayer* layer = dtGeom->layer(*detUnitIt);
0095 
0096     // Get the range for the corresponding LayerId
0097     DTRecHitCollection::range range = dtRecHits->get((*detUnitIt));
0098     // Loop over the rechits of this DetUnit
0099     for (DTRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
0100       // Get the wireId of the rechit
0101       DTWireId wireId = (*rechit).wireId();
0102 
0103       float xwire = layer->specificTopology().wirePosition(wireId.wire());
0104       if (fabs(xwire - (*rechit).localPosition().x()) > 0.00001) {
0105         cout << "  [DTRecHitReader]***Error in wire Position: xwire = " << xwire
0106              << " xRecHit = " << (*rechit).localPosition().x() << endl;
0107       }
0108 
0109       // Access to Right and left rechits
0110       pair<const DTRecHit1D*, const DTRecHit1D*> lrRecHits = (*rechit).componentRecHits();
0111 
0112       cout << "Left Hit x(cm): " << lrRecHits.first->localPosition().x() << endl;
0113       cout << "Right Hit x(cm): " << lrRecHits.second->localPosition().x() << endl;
0114 
0115       // Compute the rechit distance from wire
0116       float distFromWire =
0117           fabs((*rechit).localPosition(DTEnums::Left).x() - (*rechit).localPosition(DTEnums::Right).x()) / 2.;
0118 
0119       // Search the best mu simhit and compute its distance from the wire
0120       float simHitDistFromWire = 0;
0121       if (simHitMap.find(wireId) != simHitMap.end()) {
0122         const PSimHit* muSimHit = findBestMuSimHit(layer, wireId, simHitMap[wireId], distFromWire);
0123         // Check that a mu simhit is found
0124         if (muSimHit != 0) {
0125           // Compute the simhit distance from wire
0126           simHitDistFromWire = findSimHitDist(layer, wireId, muSimHit);
0127           // Fill the histos
0128           H1DRecHit* histo = 0;
0129           if (wireId.superlayer() == 1 || wireId.superlayer() == 3) {
0130             histo = hRHitPhi;
0131           } else if (wireId.superlayer() == 2) {
0132             hRHitZ_All->Fill(distFromWire, simHitDistFromWire);
0133             if (wireId.wheel() == 0) {
0134               histo = hRHitZ_W0;
0135             } else if (abs(wireId.wheel()) == 1) {
0136               histo = hRHitZ_W1;
0137             } else if (abs(wireId.wheel()) == 2) {
0138               histo = hRHitZ_W2;
0139             }
0140           }
0141           histo->Fill(distFromWire, simHitDistFromWire);
0142 
0143           if (fabs(distFromWire - simHitDistFromWire) > 2.1) {
0144             cout << "Warning: " << endl
0145                  << "  RecHit distance from wire is: " << distFromWire << endl
0146                  << "  SimHit distance from wire is: " << simHitDistFromWire << endl
0147                  << "  RecHit wire Id is: " << wireId << endl
0148                  << "  SimHit wire Id is: " << DTWireId((*muSimHit).detUnitId()) << endl;
0149             cout << "  Wire x: = " << xwire << endl
0150                  << "  RecHit x = " << (*rechit).localPosition().x() << endl
0151                  << "  SimHit x = " << (*muSimHit).localPosition().x() << endl;
0152           }
0153 
0154           // Some printout
0155           if (debug) {
0156             cout << "[DTRecHitReader]: " << endl
0157                  << "         WireId: " << wireId << endl
0158                  << "         1DRecHitPair local position (cm): " << (*rechit) << endl
0159                  << "         RecHit distance from wire (cm): " << distFromWire << endl
0160                  << "         Mu SimHit distance from wire (cm): " << simHitDistFromWire << endl;
0161           }
0162         }
0163       }
0164     }
0165   }
0166 }
0167 
0168 // Return a map between simhits of a layer and the wireId of their cell
0169 map<DTWireId, vector<const PSimHit*> > DTRecHitReader::mapSimHitsPerWire(const Handle<PSimHitContainer>& simhits) {
0170   map<DTWireId, vector<const PSimHit*> > hitWireMapResult;
0171 
0172   for (PSimHitContainer::const_iterator simhit = simhits->begin(); simhit != simhits->end(); simhit++) {
0173     hitWireMapResult[DTWireId((*simhit).detUnitId())].push_back(&(*simhit));
0174   }
0175 
0176   return hitWireMapResult;
0177 }
0178 
0179 const PSimHit* DTRecHitReader::findBestMuSimHit(const DTLayer* layer,
0180                                                 const DTWireId& wireId,
0181                                                 const vector<const PSimHit*>& simhits,
0182                                                 float recHitDistFromWire) {
0183   const PSimHit* retSimHit = 0;
0184   float tmp_distDiff = 999999;
0185   for (vector<const PSimHit*>::const_iterator simhit = simhits.begin(); simhit != simhits.end(); simhit++) {
0186     // Select muons
0187     if (abs((*simhit)->particleType()) == 13) {
0188       // Get the mu simhit closest to the rechit
0189       if (findSimHitDist(layer, wireId, *simhit) - recHitDistFromWire < tmp_distDiff) {
0190         tmp_distDiff = findSimHitDist(layer, wireId, *simhit) - recHitDistFromWire;
0191         retSimHit = (*simhit);
0192       }
0193     }
0194   }
0195   return retSimHit;
0196 }
0197 
0198 // Compute SimHit distance from wire
0199 double DTRecHitReader::findSimHitDist(const DTLayer* layer, const DTWireId& wireId, const PSimHit* hit) {
0200   float xwire = layer->specificTopology().wirePosition(wireId.wire());
0201   LocalPoint entryP = hit->entryPoint();
0202   LocalPoint exitP = hit->exitPoint();
0203   float xEntry = entryP.x() - xwire;
0204   float xExit = exitP.x() - xwire;
0205 
0206   return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));  //FIXME: check...
0207 }
0208 
0209 DEFINE_FWK_MODULE(DTRecHitReader);