Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-10 22:24:54 

0001 /******* \class STAnalyzer *******
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/STAnalyzer.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 "DataFormats/LTCDigi/interface/LTCDigi.h"
0031 
0032 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0033 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0034 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0035 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0036 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0037 
0038 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0039 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0040 #include "MagneticField/Engine/interface/MagneticField.h"
0041 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0044 
0045 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0046 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0047 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0048 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0049 
0050 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 
0053 /* C++ Headers */
0054 #include <iostream>
0055 #include <cmath>
0056 using namespace std;
0057 
0058 /* ====================================================================== */
0059 
0060 /* Constructor */
0061 STAnalyzer::STAnalyzer(const ParameterSet& pset) : _ev(0) {
0062   if (debug)
0063     cout << "STAnalyzer::STAnalyzer" << endl;
0064 
0065   // Get the debug parameter for verbose output
0066   debug = pset.getUntrackedParameter<bool>("debug");
0067   theRootFileName = pset.getUntrackedParameter<string>("rootFileName");
0068 
0069   // the name of the 1D rec hits collection
0070   theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");
0071 
0072   // the name of the 2D rec hits collection
0073   theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");
0074 
0075   // the name of the 4D rec hits collection
0076   theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");
0077 
0078   // the name of the STA rec hits collection
0079   theSTAMuonLabel = pset.getParameter<string>("SALabel");
0080 
0081   thePropagatorName = pset.getParameter<std::string>("PropagatorName");
0082   thePropagator = 0;
0083   thePropagatorToken = esConsumes(edm::ESInputTag("", thePropagatorName));
0084 
0085   doSA = pset.getParameter<bool>("doSA");
0086 
0087   // Create the root file
0088   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0089   bool dirStat = TH1::AddDirectoryStatus();
0090   TH1::AddDirectory(kTRUE);
0091 
0092   // CosmicMuon
0093   new TH1F("hNSA", "Num SA tracks in events", 6, -0.5, 5.5);
0094   new TH1F("hNSADifferentSectors", "Num SA tracks with hits in different sectors", 6, -0.5, 5.5);
0095   new TH1F("hNhitsSA", "Num hits in SA tracks", 50, .5, 50.5);
0096   new TH1F("hChi2SA", "#chi^{2}/NDoF for SA tracks", 100, 0, 100.);
0097   new TH1F("hPIPSA", "p for SA tracks @ IP", 100, 0., 100);
0098   new TH1F("hPtIPSA", "pt for SA tracks @ IP", 100, 0., 100);
0099   new TH1F("hPhiIPSA", "#phi for SA tracks @ IP", 100, -M_PI_2, M_PI_2);
0100   new TH1F("hEtaIPSA", "#eta for SA tracks @ IP", 100, -2.5, 2.5);
0101   new TH1F("hPInnerTSOSSA", "p for SA tracks @ innermost TSOS", 100, 0., 100);
0102   new TH1F("hPtInnerTSOSSA", "pt for SA tracks @ innermost TSOS", 100, 0., 100);
0103   new TH1F("hPhiInnerTSOSSA", "#phi for SA tracks @ innermost TSOS", 100, -M_PI_2, M_PI_2);
0104   new TH1F("hEtaInnerTSOSSA", "#eta for SA tracks @ innermost TSOS", 100, -2.5, 2.5);
0105   new TH1F("hInnerRSA", "Radius of innermost TSOS for SA tracks", 100, 400, 1000.);
0106   new TH1F("hOuterRSA", "Radius of outermost TSOS for SA tracks", 100, 400, 1000.);
0107   new TH2F("hInnerOuterRSA", "Radius of outermost TSOS vs innermost one for SA tracks", 40, 300, 1000., 40, 300, 1000.);
0108 
0109   new TH2F("hNSAVsNHits", "Num 1d Hits vs N SA", 100, -0.5, 39.5, 5, -0.5, 4.5);
0110   new TH2F("hNSAVsNSegs2D", "Num 2d Segs vs N SA", 20, -0.5, 19.5, 5, -0.5, 4.5);
0111   new TH2F("hNSAVsNSegs4D", "Num 4d Segs vs N SA", 10, -0.5, 9.5, 5, -0.5, 4.5);
0112 
0113   new TH2F("hNHitsSAVsNHits", "Num 1d Hits vs N Hits SA", 40, -0.5, 39.5, 100, 0, 100);
0114   new TH2F("hNHitsSAVsNSegs2D", "Num 2d Segs vs N Hits SA", 20, -0.5, 19.5, 100, 0, 100);
0115   new TH2F("hNHitsSAVsNSegs4D", "Num 4d Segs vs N Hits SA", 10, -0.5, 9.5, 100, 0, 100);
0116 
0117   new TH1F("hSANLayers", "Num Layers used SA", 50, -0.5, 49.5);
0118   new TH1F("hSANSLs", "Num SuperLayers used SA", 15, -0.5, 14.5);
0119   new TH1F("hSANLayersPerSL", "Num Layers per SL used SA", 15, -0.5, 14.5);
0120   new TH1F("hSANChambers", "Num Chambers used SA", 6, -0.5, 5.5);
0121   new TH1F("hSANLayersPerChamber", "Num Layers per SL used SA", 15, -0.5, 14.5);
0122   new TH1F("hSANSLsPerChamber", "Num SLs per Chamber used SA", 5, -0.5, 4.5);
0123 
0124   new TH2F("hSANLayersVsNChambers", "Num Layers Vs Num Chambers used SA", 6, -0.5, 5.5, 50, -0.5, 49.5);
0125   new TH2F("hSANSLsVsNChambers", "Num SuperLayers Vs Num Chambers used SA", 6, -0.5, 5.5, 16, -0.5, 15.5);
0126 
0127   new TH2F("hSANHitsVsNLayers", "Num Hits vs Num Layers used SA", 50, -0.5, 49.5, 50, 0., 50.);
0128   new TH2F("hSANHitsVsNSLs", "Num Hits vs Num SuperLayers used SA", 16, -0.5, 15.5, 50, 0., 50.);
0129   new TH2F("hSANHitsVsNChambers", "Num Hits vs Num Chambers used SA", 6, -0.5, 5.5, 50, 0., 50.);
0130 
0131   new TH2F("hHitsPosXYSA", "Hits position (x,y) SA", 100, -800, 800, 100, -800, 800);
0132   new TH2F("hHitsPosXZSA", "Hits position (x,z) SA", 100, -800, 800, 100, -670, 670);
0133   new TH2F("hHitsPosYZSA", "Hits position (y,z) SA", 100, -670, 670, 100, -800, 800);
0134 
0135   new TH1F("nHitsSAInChamber", "Num STA recHits in each chamber", 4, 0.5, 4.5);
0136   new TH1F("nHitsSAInSL", "Num STA recHits in each sl", 50, 0.5, 50.5);
0137   new TH1F("nHitsSAInLayer", "Num STA recHits in each layer", 75, 25.5, 100.5);
0138 
0139   new TH1F("hHitsLostChamber", "Num hits in det w/o associated hits chamber", 30, 0., 30.);
0140   new TH1F("hHitsLostSL", "Num hits in det w/o associated hits SL", 30, 0., 30.);
0141   new TH1F("hHitsLostLayer", "Num hits in det w/o associated hits Layer", 30, 0., 30.);
0142 
0143   new TH1F("hMinDistChamberFound",
0144            "Min distance between extrap post and closest hit in det w associated hits chamber",
0145            100,
0146            0.,
0147            300.);
0148   new TH1F("hMinDistChamberNotFound",
0149            "Min distance between extrap post and closest hit in det w/o associated hits chamber",
0150            100,
0151            0.,
0152            300.);
0153   new TH1F("hMinPhiChamberFound",
0154            "Min phi between extrap post and closest hit in det w associated hits chamber",
0155            100,
0156            -1.,
0157            1.);
0158   new TH1F("hMinPhiChamberNotFound",
0159            "Min phi between extrap post and closest hit in det w/o associated hits chamber",
0160            100,
0161            -1.,
0162            1.);
0163   new TH1F("hMinThetaChamberFound",
0164            "Min theta between extrap post and closest hit in det w associated hits chamber",
0165            100,
0166            -1.,
0167            1.);
0168   new TH1F("hMinThetaChamberNotFound",
0169            "Min theta between extrap post and closest hit in det w/o associated hits chamber",
0170            100,
0171            -1.,
0172            1.);
0173 
0174   new TH1F(
0175       "hMinDistSLFound", "Min distance between extrap post and closest hit in det w associated hits SL", 100, 0., 300.);
0176   new TH1F("hMinDistSLNotFound",
0177            "Min distance between extrap post and closest hit in det w/o associated hits SL",
0178            100,
0179            0.,
0180            300.);
0181   new TH1F(
0182       "hMinThetaSLFound", "Min theta between extrap post and closest hit in det w associated hits SL", 100, -1., 1.);
0183   new TH1F("hMinThetaSLNotFound",
0184            "Min theta between extrap post and closest hit in det w/o associated hits SL",
0185            100,
0186            -1.,
0187            1.);
0188 
0189   new TH1F("hMinDistLayerFound",
0190            "Min distance between extrap post and closest hit in det w associated hits Layer",
0191            100,
0192            0.,
0193            300.);
0194   new TH1F("hMinDistLayerNotFound",
0195            "Min distance between extrap post and closest hit in det w/o associated hits Layer",
0196            100,
0197            0.,
0198            300.);
0199 
0200   TH1::AddDirectory(dirStat);
0201 
0202   theDtGeomTokenBR = esConsumes<edm::Transition::BeginRun>();
0203   theDtGeomToken = esConsumes();
0204   theMGFieldToken = esConsumes();
0205   theTrackingGeometryToken = esConsumes();
0206 }
0207 
0208 /* Destructor */
0209 STAnalyzer::~STAnalyzer() {
0210   theFile->cd();
0211   theFile->Write();
0212   theFile->Close();
0213 }
0214 
0215 void STAnalyzer::beginRun(const edm::Run& run, const EventSetup& setup) {
0216   // Get the DT Geometry
0217   ESHandle<DTGeometry> dtGeom = setup.getHandle(theDtGeomTokenBR);
0218 
0219   if (firstPass) {
0220     const std::vector<const DTChamber*>& chs = dtGeom->chambers();
0221     for (auto ch = chs.begin(); ch != chs.end(); ++ch)
0222       hitsPerChamber[(*ch)->id()] = 0;
0223   }
0224 
0225   firstPass = false;
0226 }
0227 
0228 /* Operations */
0229 void STAnalyzer::beginJob() {
0230   if (debug)
0231     cout << "STAnalyzer::beginJob" << endl;
0232 }
0233 
0234 void STAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0235   if (debug)
0236     cout << "STAnalyzer::analyze" << endl;
0237 
0238   if (debug)
0239     cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() << " Num " << _ev++ << endl;
0240 
0241   if (doSA)
0242     analyzeSATrack(event, eventSetup);
0243 }
0244 
0245 void STAnalyzer::analyzeSATrack(const Event& event, const EventSetup& eventSetup) {
0246   if (debug)
0247     cout << "STAnalyzer::analyzeSATrack" << endl;
0248   if (!thePropagator) {
0249     ESHandle<Propagator> prop = eventSetup.getHandle(thePropagatorToken);
0250     thePropagator = prop->clone();
0251     thePropagator->setPropagationDirection(anyDirection);
0252   }
0253 
0254   MuonPatternRecoDumper muonDumper;
0255 
0256   // Get the 4D rechit collection from the event -------------------
0257   edm::Handle<DTRecSegment4DCollection> segs;
0258   event.getByLabel(theRecHits4DLabel, segs);
0259 
0260   // Get the 2D rechit collection from the event -------------------
0261   edm::Handle<DTRecSegment2DCollection> segs2d;
0262   event.getByLabel(theRecHits2DLabel, segs2d);
0263 
0264   // Get the 1D rechits from the event --------------
0265   Handle<DTRecHitCollection> dtRecHits;
0266   event.getByLabel(theRecHits1DLabel, dtRecHits);
0267 
0268   // Get the RecTrack collection from the event
0269   Handle<reco::TrackCollection> staTracks;
0270   event.getByLabel(theSTAMuonLabel, staTracks);
0271   if (debug)
0272     cout << "SA " << staTracks->size() << endl;
0273   // Get the RecTrack collection from the event
0274   std::vector<Handle<reco::TrackCollection> > tracks;
0275   event.getManyByType(tracks);
0276   std::vector<Handle<reco::TrackCollection> >::iterator i(tracks.begin()), end(tracks.end());
0277   for (; i != end; ++i) {
0278     if (debug)
0279       cout << "ALL Tracks id " << (*i).id() << endl;
0280     if (debug)
0281       cout << "ALL Tracks size " << (*i).product()->size() << endl;
0282   }
0283 
0284   ESHandle<MagneticField> theMGField = eventSetup.getHandle(theMGFieldToken);
0285 
0286   ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(theTrackingGeometryToken);
0287 
0288   reco::TrackCollection::const_iterator staTrack;
0289 
0290   double recPt = 0.;
0291   histo("hNSA")->Fill(staTracks->size());
0292   histo2d("hNSAVsNHits")->Fill(dtRecHits->size(), staTracks->size());
0293   histo2d("hNSAVsNSegs2D")->Fill(segs2d->size(), staTracks->size());
0294   histo2d("hNSAVsNSegs4D")->Fill(segs->size(), staTracks->size());
0295 
0296   if (debug && staTracks->size())
0297     cout << endl << "R:E " << event.id().run() << ":" << event.id().event() << " SA " << staTracks->size() << endl;
0298 
0299   for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack) {
0300     reco::TransientTrack track(*staTrack, &*theMGField, theTrackingGeometry);
0301 
0302     if (debug) {
0303       cout << muonDumper.dumpFTS(track.impactPointTSCP().theState());
0304 
0305       recPt = track.impactPointTSCP().momentum().perp();
0306       cout << " p: " << track.impactPointTSCP().momentum().mag() << " pT: " << recPt << endl;
0307       cout << " normalized chi2: " << track.normalizedChi2() << endl;
0308     }
0309 
0310     histo("hChi2SA")->Fill(track.normalizedChi2());
0311 
0312     histo("hPIPSA")->Fill(track.impactPointTSCP().momentum().mag());
0313     histo("hPtIPSA")->Fill(recPt);
0314     histo("hPhiIPSA")->Fill(track.impactPointTSCP().momentum().phi());
0315     histo("hEtaIPSA")->Fill(track.impactPointTSCP().momentum().eta());
0316 
0317     TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
0318     if (debug) {
0319       cout << "Inner TSOS:" << endl;
0320       cout << muonDumper.dumpTSOS(innerTSOS);
0321       cout << " p: " << innerTSOS.globalMomentum().mag() << " pT: " << innerTSOS.globalMomentum().perp() << endl;
0322     }
0323 
0324     // try to extrapolate using thePropagator
0325 
0326     // Get a surface (here a cylinder of radius 1290mm) ECAL
0327     float radiusECAL = 129.0;  // radius in centimeter
0328     Cylinder::PositionType pos0;
0329     Cylinder::RotationType rot0;
0330     const Cylinder::CylinderPointer ecal = Cylinder::build(radiusECAL, pos0, rot0);
0331     //cout << "Cyl " << ecal->radius() << endl;
0332 
0333     TrajectoryStateOnSurface tsosAtEcal = thePropagator->propagate(*innerTSOS.freeState(), *ecal);
0334     // if (tsosAtEcal.isValid())
0335     //   cout << "extrap to ECAL " << tsosAtEcal.globalPosition() << " r " <<
0336     //     tsosAtEcal.globalPosition().perp() << endl;
0337     // else
0338     //   cout << "Extrapolation to ECAL failed" << endl;
0339 
0340     // Get a surface (here a cylinder of radius 1811 mm) HCAL
0341     float radiusHCAL = 181.1;  // radius in centimeter
0342     const Cylinder::CylinderPointer hcal = Cylinder::build(radiusHCAL, pos0, rot0);
0343     //cout << "Cyl " << hcal->radius() << endl;
0344 
0345     TrajectoryStateOnSurface tsosAtHcal = thePropagator->propagate(*innerTSOS.freeState(), *hcal);
0346     // if (tsosAtHcal.isValid())
0347     //   cout << "extrap to HCAL " << tsosAtHcal.globalPosition() << " r " <<
0348     //     tsosAtHcal.globalPosition().perp() << endl;
0349     // else
0350     //   cout << "Extrapolation to HCAL failed" << endl;
0351 
0352     histo("hPInnerTSOSSA")->Fill(innerTSOS.globalMomentum().mag());
0353     histo("hPtInnerTSOSSA")->Fill(innerTSOS.globalMomentum().perp());
0354     histo("hPhiInnerTSOSSA")->Fill(innerTSOS.globalMomentum().phi());
0355     histo("hEtaInnerTSOSSA")->Fill(innerTSOS.globalMomentum().eta());
0356 
0357     histo("hNhitsSA")->Fill(staTrack->recHitsSize());
0358 
0359     histo2d("hNHitsSAVsNHits")->Fill(dtRecHits->size(), staTrack->recHitsSize());
0360     histo2d("hNHitsSAVsNSegs2D")->Fill(segs2d->size(), staTrack->recHitsSize());
0361     histo2d("hNHitsSAVsNSegs4D")->Fill(segs->size(), staTrack->recHitsSize());
0362 
0363     histo("hInnerRSA")->Fill(sqrt(staTrack->innerPosition().perp2()));
0364     histo("hOuterRSA")->Fill(sqrt(staTrack->outerPosition().perp2()));
0365     histo2d("hInnerOuterRSA")->Fill(sqrt(staTrack->innerPosition().perp2()), sqrt(staTrack->outerPosition().perp2()));
0366 
0367     if (debug)
0368       cout << "RecHits:" << endl;
0369 
0370     trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
0371     trackingRecHit_iterator rhend = staTrack->recHitsEnd();
0372 
0373     // zero's the maps
0374     for (std::map<DTChamberId, int>::iterator h = hitsPerChamber.begin(); h != hitsPerChamber.end(); ++h)
0375       (*h).second = 0;
0376     for (std::map<DTSuperLayerId, int>::iterator h = hitsPerSL.begin(); h != hitsPerSL.end(); ++h)
0377       (*h).second = 0;
0378     for (std::map<DTLayerId, int>::iterator h = hitsPerLayer.begin(); h != hitsPerLayer.end(); ++h)
0379       (*h).second = 0;
0380 
0381     int firstHitWheel = 0;   // the Wheel of first hit
0382     int lastHitWheel = 0;    // the Wheel of last hit
0383     int firstHitSector = 0;  // the sector of first hit
0384     int lastHitSector = 0;   // the sector of last hit
0385     for (trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit) {
0386       const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0387       GlobalPoint gpos = geomDet->toGlobal((*recHit)->localPosition());
0388       if (debug)
0389         cout << "r: " << gpos.perp() << " z: " << gpos.z() << endl;
0390       histo2d("hHitsPosXYSA")->Fill(gpos.x(), gpos.y());
0391       histo2d("hHitsPosXZSA")->Fill(gpos.z(), gpos.x());
0392       histo2d("hHitsPosYZSA")->Fill(gpos.z(), gpos.y());
0393 
0394       if (const DTLayer* lay = dynamic_cast<const DTLayer*>(geomDet)) {
0395         if (firstHitSector == 0)
0396           firstHitSector = lay->id().sector();
0397         lastHitSector = lay->id().sector();
0398         if (firstHitWheel == 0)
0399           firstHitWheel = lay->id().wheel();
0400         lastHitWheel = lay->id().wheel();
0401 
0402         if (debug)
0403           cout << "Layer " << lay->id() << endl;
0404         histo("nHitsSAInChamber")->Fill(lay->id().station());
0405         hitsPerChamber[lay->id()]++;
0406         histo("nHitsSAInSL")->Fill(lay->id().station() * 10 + lay->id().superlayer());
0407         hitsPerSL[lay->id()]++;
0408         histo("nHitsSAInLayer")->Fill(lay->id().station() * 20 + lay->id().superlayer() * 5 + lay->id().layer());
0409         hitsPerLayer[lay->id()]++;
0410       }
0411       if (const DTSuperLayer* lay = dynamic_cast<const DTSuperLayer*>(geomDet))
0412         cout << "SuperLayer " << lay->id() << endl;
0413       if (const DTChamber* lay = dynamic_cast<const DTChamber*>(geomDet))
0414         cout << "Chamber " << lay->id() << endl;
0415 
0416       TrajectoryStateOnSurface extraptsos = thePropagator->propagate(innerTSOS, geomDet->specificSurface());
0417     }
0418     if (debug) {
0419       cout << "PerChamber " << muonDumper.dumpTSOS(innerTSOS) << endl;
0420       for (std::map<DTChamberId, int>::iterator h = hitsPerChamber.begin(); h != hitsPerChamber.end(); ++h)
0421         if ((*h).second)
0422           cout << (*h).first << ":" << (*h).second << endl;
0423       cout << "=====" << endl;
0424 
0425       cout << "PerSL " << endl;
0426       for (std::map<DTSuperLayerId, int>::iterator h = hitsPerSL.begin(); h != hitsPerSL.end(); ++h)
0427         if ((*h).second)
0428           cout << (*h).first << ":" << (*h).second << endl;
0429       cout << "=====" << endl;
0430       cout << "PerLayer " << endl;
0431       for (std::map<DTLayerId, int>::iterator h = hitsPerLayer.begin(); h != hitsPerLayer.end(); ++h)
0432         if ((*h).second)
0433           cout << (*h).first << ":" << (*h).second << endl;
0434       cout << "=====" << endl;
0435     }
0436 
0437     // Get the DT Geometry
0438     ESHandle<DTGeometry> dtGeom = eventSetup.getHandle(theDtGeomToken);
0439 
0440     // Get the 1D rechits from the event --------------
0441     Handle<DTRecHitCollection> hits1d;
0442     event.getByLabel(theRecHits1DLabel, hits1d);
0443 
0444     // Get the 2D rechit collection from the event -------------------
0445     edm::Handle<DTRecSegment2DCollection> segs2d;
0446     event.getByLabel(theRecHits2DLabel, segs2d);
0447 
0448     // Get the 4D rechit collection from the event -------------------
0449     edm::Handle<DTRecSegment4DCollection> segs;
0450     event.getByLabel(theRecHits4DLabel, segs);
0451 
0452     // Now loop on chamber of relevant sector and look if there are hits in all
0453     // of them
0454     if (firstHitWheel == lastHitWheel &&
0455         ((firstHitSector == lastHitSector) || (firstHitSector == 10 && lastHitSector == 14) ||
0456          (firstHitSector == 14 && lastHitSector == 10) || (firstHitSector == 13 && lastHitSector == 4) ||
0457          (firstHitSector == 4 && lastHitSector == 13))) {
0458       int nLay = 0, nSL = 0, nCh = 0;
0459       for (int c = 1; c <= 5; ++c) {
0460         // cout << "c, lh, fh " << c <<" " << firstHitSector << " " <<
0461         //   lastHitSector << endl;
0462         if (c == 5 && firstHitSector == 14)
0463           firstHitSector = 10;
0464         else if (c == 5 && lastHitSector == 14)
0465           lastHitSector = 10;
0466         else if (c == 5 && firstHitSector == 13)
0467           firstHitSector = 4;
0468         else if (c == 5 && lastHitSector == 13)
0469           lastHitSector = 4;
0470         else if (c == 5)
0471           continue;
0472         int cc = c;
0473         if (c == 5)
0474           cc = 4;
0475 
0476         int sect = min(firstHitSector, lastHitSector);
0477         DTChamberId chid(firstHitWheel, cc, sect);
0478         const DTChamber* ch = dtGeom->chamber(chid);
0479         if (hitsPerChamber[chid] == 0) {
0480           missingHit(dtGeom, segs, ch, innerTSOS, false);
0481         } else {
0482           missingHit(dtGeom, segs, ch, innerTSOS, true);
0483           nCh++;
0484           //cout << "Hits in " << chid << " = " << hitsPerChamber[chid]<< endl;
0485         }
0486         int nLayPerCh = 0, nSLPerCh = 0;
0487         for (int sl = 1; sl <= 3; ++sl) {
0488           if (sl == 2 && cc == 4)
0489             continue;
0490           DTSuperLayerId slid(chid, sl);
0491           const DTSuperLayer* dtsl = dtGeom->superLayer(slid);
0492           if (hitsPerSL[slid] == 0) {
0493             //cout << "Hits missing in " << slid << endl;
0494             missingHit(dtGeom, segs2d, dtsl, innerTSOS, false);
0495           } else {
0496             missingHit(dtGeom, segs2d, dtsl, innerTSOS, true);
0497             nSL++;
0498             nSLPerCh++;
0499             //cout << "Hits in " << slid << " = " << hitsPerSL[slid]<< endl;
0500           }
0501           int nLayPerSL = 0;
0502           for (int l = 1; l <= 4; ++l) {
0503             DTLayerId lid(slid, l);
0504             const DTLayer* lay = dtGeom->layer(lid);
0505             if (hitsPerLayer[lid] == 0) {
0506               //cout << "Hits missing in " << lid << endl;
0507               missingHit(dtGeom, hits1d, lay, innerTSOS, false);
0508             } else {
0509               missingHit(dtGeom, hits1d, lay, innerTSOS, true);
0510               nLay++;
0511               nLayPerCh++;
0512               nLayPerSL++;
0513               //cout << "Hits in " << lid << " = " << hitsPerLayer[lid]<< endl;
0514             }
0515           }
0516           histo("hSANLayersPerSL")->Fill(nLayPerSL);
0517         }
0518         histo("hSANLayersPerChamber")->Fill(nLayPerCh);
0519         histo("hSANSLsPerChamber")->Fill(nSLPerCh);
0520       }
0521 
0522       // how many different Chambers, SL, Layers are used in a SA?
0523       histo("hSANLayers")->Fill(nLay);
0524       histo("hSANSLs")->Fill(nSL);
0525       histo("hSANChambers")->Fill(nCh);
0526       histo2d("hSANLayersVsNChambers")->Fill(nCh, nLay);
0527       histo2d("hSANSLsVsNChambers")->Fill(nCh, nSL);
0528       // n det vs num hits used
0529       histo2d("hSANHitsVsNLayers")->Fill(nLay, staTrack->recHitsSize());
0530       histo2d("hSANHitsVsNSLs")->Fill(nSL, staTrack->recHitsSize());
0531       histo2d("hSANHitsVsNChambers")->Fill(nCh, staTrack->recHitsSize());
0532       if (nCh > 4)
0533         cout << endl << "R:E " << event.id().run() << ":" << event.id().event() << " nCh " << nCh << endl;
0534     } else {  // different wheel/sector
0535       // cout << "First/Last hit in different wheel/sector " << firstHitSector
0536       //    << " " << lastHitSector<< endl;
0537       histo("hNSADifferentSectors")->Fill(1);
0538     }
0539   }
0540   if (debug)
0541     cout << "---" << endl;
0542 }
0543 
0544 template <typename T, typename C>
0545 void STAnalyzer::missingHit(const ESHandle<DTGeometry>& dtGeom,
0546                             const Handle<C>& segs,
0547                             const T* ch,
0548                             const TrajectoryStateOnSurface& startTsos,
0549                             bool found) {
0550   if (!ch)
0551     return;
0552   TrajectoryStateOnSurface extraptsos = thePropagator->propagate(startTsos, ch->specificSurface());
0553   //cout << "extrap " << extraptsos.isValid() << endl;
0554   if (extraptsos.isValid()) {
0555     //cout << "extraptsos " << extraptsos.localPosition() << endl;
0556     bool inside = ch->specificSurface().bounds().inside(extraptsos.localPosition());
0557     //cout << "Is extrap inside? " << inside << endl;
0558     if (inside) {
0559       // Here I should loop on all hits of this chamber (if any) and get
0560       // the closest
0561       typename C::range segsch = segs->get(ch->id());
0562       //int nSegsCh=segsch.second-segsch.first;
0563       //cout << "Hits in chamber " << segsch.second-segsch.first << endl;
0564       fillMinDist(segsch, ch, extraptsos, found);
0565     }
0566   }
0567 }
0568 
0569 void STAnalyzer::fillMinDist(const DTRecSegment4DCollection::range segsch,
0570                              const DTChamber* ch,
0571                              const TrajectoryStateOnSurface& extraptsos,
0572                              bool found) {
0573   int nSegsCh = segsch.second - segsch.first;
0574   histo("hHitsLostChamber")->Fill(nSegsCh);
0575   if (nSegsCh) {
0576     LocalPoint extrapPos = extraptsos.localPosition();
0577     LocalVector extrapDir = extraptsos.localDirection();
0578     //cout << "Extrap pos " << extrapPos << endl;
0579     double minDist = 99999.;
0580     double minPhi = 0., minTheta = 0.;
0581     for (DTRecSegment4DCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0582       //cout << "Hit pos " << hit->localPosition() << endl;
0583       LocalVector dist = hit->localPosition() - extrapPos;
0584       //cout << "dist " << dist << " " << dist.perp() << endl;
0585       if (dist.perp() < minDist) {
0586         minDist = dist.perp();
0587         minPhi = hit->localDirection().phi() - extrapDir.phi();
0588         minTheta = hit->localDirection().theta() - extrapDir.theta();
0589       }
0590     }
0591     string inOut = (found) ? "Found" : "NotFound";
0592     histo(string("hMinDistChamber") + inOut)->Fill(minDist);
0593     histo(string("hMinPhiChamber") + inOut)->Fill(minPhi);
0594     histo(string("hMinThetaChamber") + inOut)->Fill(minTheta);
0595   }
0596 }
0597 
0598 void STAnalyzer::fillMinDist(const DTRecSegment2DCollection::range segsch,
0599                              const DTSuperLayer* ch,
0600                              const TrajectoryStateOnSurface& extraptsos,
0601                              bool found) {
0602   int nSegsCh = segsch.second - segsch.first;
0603   histo("hHitsLostSL")->Fill(nSegsCh);
0604   if (nSegsCh) {
0605     LocalPoint extrapPos = extraptsos.localPosition();
0606     LocalVector extrapDir = extraptsos.localDirection();
0607     //cout << "Extrap pos " << extrapPos << endl;
0608     double minDist = 99999.;
0609     double minTheta = 0.;
0610     for (DTRecSegment2DCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0611       //cout << "Hit pos " << hit->localPosition() << endl;
0612       LocalVector dist = hit->localPosition() - extrapPos;
0613       //cout << "dist " << dist << " " << dist.perp() << endl;
0614       if (dist.perp() < minDist) {
0615         minDist = dist.x();
0616         minTheta = hit->localDirection().theta() - extrapDir.theta();
0617       }
0618     }
0619     string inOut = (found) ? "Found" : "NotFound";
0620     histo(string("hMinDistSL") + inOut)->Fill(minDist);
0621     histo(string("hMinThetaSL") + inOut)->Fill(minTheta);
0622   }
0623 }
0624 
0625 void STAnalyzer::fillMinDist(const DTRecHitCollection::range segsch,
0626                              const DTLayer* ch,
0627                              const TrajectoryStateOnSurface& extraptsos,
0628                              bool found) {
0629   int nSegsCh = segsch.second - segsch.first;
0630   histo("hHitsLostLayer")->Fill(nSegsCh);
0631   if (nSegsCh) {
0632     LocalPoint extrapPos = extraptsos.localPosition();
0633     //cout << "Extrap pos " << extrapPos << endl;
0634     double minDist = 99999.;
0635     for (DTRecHitCollection::const_iterator hit = segsch.first; hit != segsch.second; ++hit) {
0636       //cout << "Hit pos " << hit->localPosition() << endl;
0637       LocalVector dist = hit->localPosition() - extrapPos;
0638       //cout << "dist " << dist << " " << dist.perp() << endl;
0639       if (dist.perp() < minDist)
0640         minDist = dist.x();
0641     }
0642     string inOut = (found) ? "Found" : "NotFound";
0643     histo(string("hMinDistLayer") + inOut)->Fill(minDist);
0644   }
0645 }
0646 
0647 TH1F* STAnalyzer::histo(const string& name) const {
0648   if (TH1F* h = dynamic_cast<TH1F*>(theFile->Get(name.c_str())))
0649     return h;
0650   else
0651     throw cms::Exception("STAnalyzer") << " Not a TH1F " << name;
0652 }
0653 
0654 TH2F* STAnalyzer::histo2d(const string& name) const {
0655   if (TH2F* h = dynamic_cast<TH2F*>(theFile->Get(name.c_str())))
0656     return h;
0657   else
0658     throw cms::Exception("STAnalyzer") << " Not a TH2F " << name;
0659 }
0660 
0661 string STAnalyzer::toString(const DTLayerId& id) const {
0662   stringstream result;
0663   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer() << "_Lay"
0664          << id.layer();
0665   return result.str();
0666 }
0667 
0668 string STAnalyzer::toString(const DTSuperLayerId& id) const {
0669   stringstream result;
0670   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station() << "_Sl" << id.superLayer();
0671   return result.str();
0672 }
0673 
0674 string STAnalyzer::toString(const DTChamberId& id) const {
0675   stringstream result;
0676   result << "_Wh" << id.wheel() << "_Sec" << id.sector() << "_St" << id.station();
0677   return result.str();
0678 }
0679 
0680 template <class T>
0681 string STAnalyzer::hName(const string& s, const T& id) const {
0682   string name(toString(id));
0683   stringstream hName;
0684   hName << s << name;
0685   return hName.str();
0686 }
0687 
0688 void STAnalyzer::createTH1F(const std::string& name,
0689                             const std::string& title,
0690                             const std::string& suffix,
0691                             int nbin,
0692                             const double& binMin,
0693                             const double& binMax) const {
0694   stringstream hName;
0695   stringstream hTitle;
0696   hName << name << suffix;
0697   hTitle << title << suffix;
0698   new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin, binMin, binMax);
0699 }
0700 
0701 void STAnalyzer::createTH2F(const std::string& name,
0702                             const std::string& title,
0703                             const std::string& suffix,
0704                             int nBinX,
0705                             const double& binXMin,
0706                             const double& binXMax,
0707                             int nBinY,
0708                             const double& binYMin,
0709                             const double& binYMax) const {
0710   stringstream hName;
0711   stringstream hTitle;
0712   hName << name << suffix;
0713   hTitle << title << suffix;
0714   new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX, binXMin, binXMax, nBinY, binYMin, binYMax);
0715 }
0716 
0717 DEFINE_FWK_MODULE(STAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!