diff --git a/CalibTracker/SiStripHitResolution/BuildFile.xml b/CalibTracker/SiStripHitResolution/BuildFile.xml index 355e76fa0de11..7c6c794054e2c 100644 --- a/CalibTracker/SiStripHitResolution/BuildFile.xml +++ b/CalibTracker/SiStripHitResolution/BuildFile.xml @@ -1,32 +1,34 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibTracker/SiStripHitResolution/interface/HitResol.h b/CalibTracker/SiStripHitResolution/interface/HitResol.h index 6d146c875328a..61bdfe67f21f3 100644 --- a/CalibTracker/SiStripHitResolution/interface/HitResol.h +++ b/CalibTracker/SiStripHitResolution/interface/HitResol.h @@ -1,144 +1,182 @@ +#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" +#include "CalibTracker/Records/interface/SiStripQualityRcd.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/DetId/interface/DetIdCollection.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/Scalers/interface/LumiScalers.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" #include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" #include "FWCore/Framework/interface/Event.h" -#include "DataFormats/Common/interface/Handle.h" #include "FWCore/Framework/interface/EventSetup.h" -#include "DataFormats/GeometryVector/interface/GlobalPoint.h" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" - +#include "FWCore/Framework/interface/one/EDAnalyzer.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h" -#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "MagneticField/Engine/interface/MagneticField.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" +#include "RecoTracker/SingleTrackPattern/interface/CosmicTrajectoryBuilder.h" #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h" -#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" -#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h" +#include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h" +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" #include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h" #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h" -#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" -#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" -#include "MagneticField/Engine/interface/MagneticField.h" #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "RecoTracker/SingleTrackPattern/interface/CosmicTrajectoryBuilder.h" -#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h" -#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h" -#include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" -#include "DataFormats/DetId/interface/DetIdCollection.h" - -#include "DataFormats/Scalers/interface/LumiScalers.h" -#include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h" - -#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" +#include "TRandom2.h" +#include "TTree.h" #include "TROOT.h" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" + #include -#include "TTree.h" #include #include #include -#include "Riostream.h" -#include "TRandom2.h" class TrackerTopology; -class HitResol : public edm::EDAnalyzer { - public: +class HitResol : public edm::one::EDAnalyzer { +public: explicit HitResol(const edm::ParameterSet& conf); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); double checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters, double xx, double xerr); bool isDoubleSided(unsigned int iidd, const TrackerTopology* tTopo) const; bool check2DPartner(unsigned int iidd, const std::vector& traj); - ~HitResol() override; + ~HitResol() override = default; unsigned int checkLayer(unsigned int iidd, const TrackerTopology* tTopo); -// float getSimHitRes(const GeomDetUnit * det, const LocalVector& trackdirection, const TrackingRecHit& recHit, float & trackWidth, float* pitch, LocalVector& drift); - void getSimHitRes(const GeomDetUnit * det, const LocalVector& trackdirection, const TrackingRecHit& recHit, float & trackWidth, float* pitch, LocalVector& drift); + void getSimHitRes(const GeomDetUnit* det, + const LocalVector& trackdirection, + const TrackingRecHit& recHit, + float& trackWidth, + float* pitch, + LocalVector& drift); double getSimpleRes(const TrajectoryMeasurement* traj1); - bool getPairParameters(const MagneticField* magField_, AnalyticalPropagator& propagator,const TrajectoryMeasurement* traj1, const TrajectoryMeasurement* traj2, float & pairPath, float & hitDX, float & trackDX, float & trackDXE, float & trackParamX, float &trackParamY , float & trackParamDXDZ, float &trackParamDYDZ , float & trackParamXE, float &trackParamYE, float & trackParamDXDZE, float &trackParamDYDZE); + bool getPairParameters(const MagneticField* magField_, + AnalyticalPropagator& propagator, + const TrajectoryMeasurement* traj1, + const TrajectoryMeasurement* traj2, + float& pairPath, + float& hitDX, + float& trackDX, + float& trackDXE, + float& trackParamX, + float& trackParamY, + float& trackParamDXDZ, + float& trackParamDYDZ, + float& trackParamXE, + float& trackParamYE, + float& trackParamDXDZE, + float& trackParamDYDZE); typedef std::vector TrajectoryCollection; - private: +private: void beginJob() override; - void endJob() override; + void endJob() override; void analyze(const edm::Event& e, const edm::EventSetup& c) override; - // ----------member data --------------------------- + // ----------member data --------------------------- + // ED tokens const edm::EDGetTokenT scalerToken_; - const edm::EDGetTokenT > commonModeToken_; - - bool addLumi_; - bool addCommonMode_; - bool cutOnTracks_; - unsigned int trackMultiplicityCut_; - bool useFirstMeas_; - bool useLastMeas_; - bool useAllHitsFromTracksWithMissingHits_; - double MomentumCut_; - unsigned int UsePairsOnly_; - - const edm::EDGetTokenT< reco::TrackCollection > combinatorialTracks_token_; - const edm::EDGetTokenT< std::vector > trajectories_token_; - const edm::EDGetTokenT< TrajTrackAssociationCollection > trajTrackAsso_token_; - const edm::EDGetTokenT< std::vector > tjToken_; + const edm::EDGetTokenT combinatorialTracks_token_; + const edm::EDGetTokenT > tjToken_; const edm::EDGetTokenT tkToken_; - const edm::EDGetTokenT< edmNew::DetSetVector > clusters_token_; - const edm::EDGetTokenT digis_token_; - const edm::EDGetTokenT trackerEvent_token_; - edm::ParameterSet conf_; - + // ES tokens + const edm::ESGetToken topoToken_; + const edm::ESGetToken geomToken_; + const edm::ESGetToken cpeToken_; + const edm::ESGetToken siStripQualityToken_; + const edm::ESGetToken magFieldToken_; + + // configuration parameters + const bool addLumi_; + const bool DEBUG_; + const bool cutOnTracks_; + const double momentumCut_; + const int compSettings_; + const unsigned int usePairsOnly_; + const unsigned int layers_; + const unsigned int trackMultiplicityCut_; + + // output file TTree* reso; TTree* treso; + std::map histos2d_; - int events,EventTrackCKF; - - int compSettings; - unsigned int layers; - bool DEBUG; - unsigned int whatlayer; - - double MomentumCut; - - -// Tree declarations -//Hit Resolution Ntuple Content - float mymom ; - int numHits ; - float ProbTrackChi2 ; - unsigned int iidd1 ; - float mypitch1 ; - unsigned int clusterWidth ; - float expWidth ; - float atEdge ; - float simpleRes ; - unsigned int iidd2 ; - unsigned int clusterWidth_2 ; - float expWidth_2 ; - float atEdge_2 ; - - float pairPath ; - float hitDX ; - float trackDX ; - float trackDXE ; - float trackParamX ; - float trackParamY ; - float trackParamDXDZ ; - float trackParamDYDZ ; - float trackParamXE ; - float trackParamYE ; - float trackParamDXDZE ; - float trackParamDYDZE ; - unsigned int pairsOnly ; - float track_momentum ; - float track_eta ; - float track_trackChi2 ; + // conversion + static constexpr float cmToUm = 10000.f; -}; + // counters + int events, EventTrackCKF; + // Tree declarations + // Hit Resolution Ntuple Content + float mymom; + int numHits; + int NumberOf_tracks; + float ProbTrackChi2; + float StripCPE1_smp_pos_error; + float StripCPE2_smp_pos_error; + float StripErrorSquared1; + float StripErrorSquared2; + float uerr2; + float uerr2_2; + unsigned int clusterWidth; + unsigned int clusterWidth_2; + unsigned int clusterCharge; + unsigned int clusterCharge_2; + unsigned int iidd1; + unsigned int iidd2; + unsigned int pairsOnly; + float pairPath; + float mypitch1; + float mypitch2; + float expWidth; + float expWidth_2; + float driftAlpha; + float driftAlpha_2; + float thickness; + float thickness_2; + float trackWidth; + float trackWidth_2; + float atEdge; + float atEdge_2; + float simpleRes; + float hitDX; + float trackDX; + float trackDXE; + float trackParamX; + float trackParamY; + float trackParamXE; + float trackParamYE; + float trackParamDXDZ; + float trackParamDYDZ; + float trackParamDXDZE; + float trackParamDYDZE; + float track_momentum; + float track_pt; + float track_eta; + float track_width; + float track_phi; + float track_trackChi2; + float N1; + float N2; + float N1uProj; + float N2uProj; + int Nstrips; + int Nstrips_2; +}; //#endif diff --git a/CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h b/CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h new file mode 100644 index 0000000000000..9dc5e3b338d9f --- /dev/null +++ b/CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h @@ -0,0 +1,61 @@ +// -*- C++ -*- +// +// Package: CalibTracker/SiStripHitResolution +// Class: SiStripOverlapHit +// +/**\class SiStripOverlapHit SiStripOverlapHit.h CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.cc + + Description: A pair of hits on overlaping modules + + Implementation: + Designed for CPE studies. Includes methods to compute residuals, etc. +*/ +// +// Original Author: Christophe Delaere +// Created: Fri, 20 Sep 2019 14:45:00 GMT +// +// +#ifndef CalibTracker_SiStripHitResolution_SiStripOverlapHit_H +#define CalibTracker_SiStripHitResolution_SiStripOverlapHit_H + +#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" +#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" + +class SiStripOverlapHit { +public: + using RecHitPointer = TrackingRecHit::RecHitPointer; + using ConstRecHitPointer = TrackingRecHit::ConstRecHitPointer; + + // constructes an overlap from 2 hits and a track. Hits are internally sorted inside-out + explicit SiStripOverlapHit(TrajectoryMeasurement const& measA, TrajectoryMeasurement const& measB); + // destructor + virtual ~SiStripOverlapHit(){}; + + // access to indivitual hits and to the trajectory state + inline ConstRecHitPointer const& hitA() const { return measA_.recHit(); } + inline ConstRecHitPointer const& hitB() const { return measB_.recHit(); } + inline ConstRecHitPointer const& hit(unsigned int hit) const { return hit ? hitB() : hitA(); } + TrajectoryStateOnSurface const& trajectoryStateOnSurface(unsigned int hit = 0, bool updated = true) const; + + // utilities + + // track local angle + double getTrackLocalAngle(unsigned int hit) const; // 0: average, 1: hit A, 2: hit B + // raw local distance between hit and strajectory state + double offset(unsigned int hit) const; + // distance between the two hits in the "trajectory frame" + double shift() const; + // absolute global distance between the hits (useful to debug pair finding) + double distance(bool fromTrajectory = false) const; + // global position (averaged over the two hits) + GlobalPoint position() const; + +private: + TrajectoryMeasurement measA_; + TrajectoryMeasurement measB_; +}; + +#endif diff --git a/CalibTracker/SiStripHitResolution/macros/ResolutionPlots.cc b/CalibTracker/SiStripHitResolution/macros/ResolutionPlots.cc index 027861260e76b..c353176b49ed8 100644 --- a/CalibTracker/SiStripHitResolution/macros/ResolutionPlots.cc +++ b/CalibTracker/SiStripHitResolution/macros/ResolutionPlots.cc @@ -1,11 +1,10 @@ #include - -void ResolutionPlots_HistoMaker(const std::string& unit){ - +void ResolutionPlots_HistoMaker(const std::string& unit) { int n = 20; - float numbers[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0}; + float numbers[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, + 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0}; std::string YAxisTitle; std::string HistoTitle; @@ -15,52 +14,51 @@ void ResolutionPlots_HistoMaker(const std::string& unit){ float UL[20] = {0}; float NonUL[20] = {0}; - if(unit == "centimetres"){ - - float UL_Array[] = {0.00137706, 0.00142804, 0.00223064, 0.00313546, 0.00528442, - 0.00328381, 0.00365521, 0.00362465, 0.00404834, 0.00304523, - 0.00154502, 0.00103244, 0.0108947, 0, 0.0030021, - 0.00229157, 0.00244635, 0.00576464, 0.00461518, 0.00464276}; - - float NonUL_Array[] = {0.00211815, 0.00243601, 0.002224, 0.00223777, 0.00284721, - 0.00306343, 0.00310523, 0.00334504, 0.00333095, 0.00326225, - 0.00331065, 0.00237821, 0.00233033, 0, 0.00285122, - 0.00287095, 0.00226949, 0.0030951 , 0.00284769, 0.00285031}; - - for(int i = 0; i < 20; i++){UL[i] = UL_Array[i]; NonUL[i] = NonUL_Array[i];} - - YAxisTitle = "Resolution [cm]"; - HistoTitle = "Resolution values for UL and non-UL samples, in centimetres"; - OutputFile = "ResolutionComparison_ULAndNonUL_Centimetres.pdf"; - RootFile = "ResolutionComparison_ULAndNonUL_Centimetres.root"; - - } - else if(unit == "pitch units"){ - - float UL_Array[] = {0.172131, 0.143091, 0.161787, 0.126722, 0.17909, - 0.197285, 0.180396, 0.170818, 0.182258, 0.166405, - 0.0844271, 0.0846207, 0.0400775, 0, 0.120119, - 0.171899, 0.160656, 0.18299, 0.177929, 0.178037}; - - - float NonUL_Array[] = {0.153758, 0.151801, 0.14859, 0.148245, 0.147986, - 0.146962, 0.147919, 0.147431, 0.146219, 0.145619, - 0.14549, 0.147042, 0.147267, 0, 0.146873, - 0.153169, 0.151639, 0.146694, 0.148681, 0.148683}; - - for(int i = 0; i < 20; i++){UL[i] = UL_Array[i]; NonUL[i] = NonUL_Array[i];} - - YAxisTitle = "Resolution [pitch units]"; - HistoTitle = "Resolution values for the UL and non-UL samples in pitch units"; - OutputFile = "ResolutionComparison_ULAndNonUL_PitchUnits.pdf"; - RootFile = "ResolutionComparison_ULAndNonUL_PitchUnits.root"; - + if (unit == "centimetres") { + float UL_Array[] = {0.00137706, 0.00142804, 0.00223064, 0.00313546, 0.00528442, 0.00328381, 0.00365521, + 0.00362465, 0.00404834, 0.00304523, 0.00154502, 0.00103244, 0.0108947, 0, + 0.0030021, 0.00229157, 0.00244635, 0.00576464, 0.00461518, 0.00464276}; + + float NonUL_Array[] = {0.00211815, 0.00243601, 0.002224, 0.00223777, 0.00284721, 0.00306343, 0.00310523, + 0.00334504, 0.00333095, 0.00326225, 0.00331065, 0.00237821, 0.00233033, 0, + 0.00285122, 0.00287095, 0.00226949, 0.0030951, 0.00284769, 0.00285031}; + + for (int i = 0; i < 20; i++) { + UL[i] = UL_Array[i]; + NonUL[i] = NonUL_Array[i]; + } + + YAxisTitle = "Resolution [cm]"; + HistoTitle = "Resolution values for UL and non-UL samples, in centimetres"; + OutputFile = "ResolutionComparison_ULAndNonUL_Centimetres.pdf"; + RootFile = "ResolutionComparison_ULAndNonUL_Centimetres.root"; + + } else if (unit == "pitch units") { + float UL_Array[] = {0.172131, 0.143091, 0.161787, 0.126722, 0.17909, 0.197285, 0.180396, + 0.170818, 0.182258, 0.166405, 0.0844271, 0.0846207, 0.0400775, 0, + 0.120119, 0.171899, 0.160656, 0.18299, 0.177929, 0.178037}; + + float NonUL_Array[] = {0.153758, 0.151801, 0.14859, 0.148245, 0.147986, 0.146962, 0.147919, + 0.147431, 0.146219, 0.145619, 0.14549, 0.147042, 0.147267, 0, + 0.146873, 0.153169, 0.151639, 0.146694, 0.148681, 0.148683}; + + for (int i = 0; i < 20; i++) { + UL[i] = UL_Array[i]; + NonUL[i] = NonUL_Array[i]; + } + + YAxisTitle = "Resolution [pitch units]"; + HistoTitle = "Resolution values for the UL and non-UL samples in pitch units"; + OutputFile = "ResolutionComparison_ULAndNonUL_PitchUnits.pdf"; + RootFile = "ResolutionComparison_ULAndNonUL_PitchUnits.root"; + + } else { + std::cout << "ERROR: Unit must be centimetres or pitch units" << std::endl; } - else{std::cout << "ERROR: Unit must be centimetres or pitch units" << std::endl;} - auto c1 = new TCanvas("c1","c1",800, 600); + auto c1 = new TCanvas("c1", "c1", 800, 600); - TGraph * gr1 = new TGraph(n, numbers, UL); + TGraph* gr1 = new TGraph(n, numbers, UL); gr1->SetName("UL samples"); gr1->SetTitle("UL samples"); gr1->SetMarkerStyle(21); @@ -69,7 +67,7 @@ void ResolutionPlots_HistoMaker(const std::string& unit){ gr1->SetLineWidth(4); gr1->SetFillStyle(0); - TGraph * gr2 = new TGraph(n, numbers, NonUL); + TGraph* gr2 = new TGraph(n, numbers, NonUL); gr2->SetName("Non-UL samples"); gr2->SetTitle("Non-UL samples"); gr2->SetMarkerStyle(22); @@ -79,54 +77,48 @@ void ResolutionPlots_HistoMaker(const std::string& unit){ gr2->SetLineWidth(0); gr2->SetFillStyle(0); - - TMultiGraph * mg = new TMultiGraph(); + TMultiGraph* mg = new TMultiGraph(); mg->Add(gr1); mg->Add(gr2); mg->GetHistogram()->SetTitle(HistoTitle.c_str()); mg->GetHistogram()->SetTitleOffset(0.05); mg->GetHistogram()->GetXaxis()->SetTitle("Region"); mg->GetHistogram()->GetXaxis()->SetTitleOffset(2.5); - c1->SetBottomMargin(0.2); - - - mg->GetHistogram()->GetXaxis()->SetBinLabel( 1.0, "TIB L1"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(2.0), "TIB L2"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(3.0), "TIB L3"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(4.0), "TIB L4"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(5.0), "Side TID"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(6.0), "Wheel TID"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(7.0), "Ring TID"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(8.0), "TOB L1"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(9.0), "TOB L2"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(10.0), "TOB L3"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(11.0), "TOB L4"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(12.0), "TOB L5"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(13.0), "TOB L6"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(14.0), "Side TEC"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(15.0), "Wheel TEC"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(16.0), "Ring TEC"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(17.0), "TIB (All)"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(18.0), "TOB (All)"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( mg->GetHistogram()->GetXaxis()->FindBin(19.0), "TID (All)"); - mg->GetHistogram()->GetXaxis()->SetBinLabel( 100, "TEC (All)"); - + c1->SetBottomMargin(0.2); + + mg->GetHistogram()->GetXaxis()->SetBinLabel(1.0, "TIB L1"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(2.0), "TIB L2"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(3.0), "TIB L3"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(4.0), "TIB L4"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(5.0), "Side TID"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(6.0), "Wheel TID"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(7.0), "Ring TID"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(8.0), "TOB L1"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(9.0), "TOB L2"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(10.0), "TOB L3"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(11.0), "TOB L4"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(12.0), "TOB L5"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(13.0), "TOB L6"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(14.0), "Side TEC"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(15.0), "Wheel TEC"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(16.0), "Ring TEC"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(17.0), "TIB (All)"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(18.0), "TOB (All)"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(mg->GetHistogram()->GetXaxis()->FindBin(19.0), "TID (All)"); + mg->GetHistogram()->GetXaxis()->SetBinLabel(100, "TEC (All)"); mg->GetHistogram()->GetYaxis()->SetTitle(YAxisTitle.c_str()); mg->Draw("ALP"); c1->BuildLegend(); mg->SaveAs(OutputFile.c_str()); - TFile * output = new TFile(RootFile.c_str(), "RECREATE"); + TFile* output = new TFile(RootFile.c_str(), "RECREATE"); output->cd(); mg->Write(); output->Close(); - } -void ResolutionPlots(){ - +void ResolutionPlots() { ResolutionPlots_HistoMaker("centimetres"); ResolutionPlots_HistoMaker("pitch units"); - } diff --git a/CalibTracker/SiStripHitResolution/macros/Resolutions.cc b/CalibTracker/SiStripHitResolution/macros/Resolutions.cc index eaa0d93b8abde..b7b720ea7cf05 100644 --- a/CalibTracker/SiStripHitResolution/macros/Resolutions.cc +++ b/CalibTracker/SiStripHitResolution/macros/Resolutions.cc @@ -1,4 +1,4 @@ -using namespace ROOT; +using namespace ROOT; using ROOT::RDF::RNode; using floats = ROOT::VecOps::RVec; using ints = ROOT::VecOps::RVec; @@ -16,10 +16,7 @@ std::string InputFileString; std::string HitResoFileName; std::string GaussianFitsFileName; - - -void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& UL){ - +void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& UL) { std::string CutFlowReportString; std::string DoubleDiffString; std::string HitDXString; @@ -28,57 +25,69 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& std::string ClusterW1String = "clusterW1"; std::string ClusterW2String = "clusterW2"; - switch(UL){ - case 0: switch(Unit_Int){ - case 0: GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO.root"; - HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO.txt"; - CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO.txt"; - DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch"; - HitDXString = "hitDX_OverPitch"; - TrackDXString = "trackDX_OverPitch"; - TrackDXEString = "trackDXE_OverPitch"; - break; - - case 1: GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO.root"; - HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO.txt"; - CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO.txt"; - DoubleDiffString = "hitDX-trackDX"; - HitDXString = "hitDX"; - TrackDXString = "trackDX"; - TrackDXEString = "trackDXE"; - break; - - default: std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl; break; - } - - InputFileString = "hitresol_ALCARECO.root"; - break; - - case 1: switch(Unit_Int){ - case 0: GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO_UL.root"; - HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO_UL.txt"; - CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO_UL.txt"; - DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch"; - HitDXString = "hitDX_OverPitch"; - TrackDXString = "trackDX_OverPitch"; - TrackDXEString = "trackDXE_OverPitch"; - break; - - case 1: GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO_UL.root"; - HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO_UL.txt"; - CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO_UL.txt"; - DoubleDiffString = "hitDX-trackDX"; - HitDXString = "hitDX"; - TrackDXString = "trackDX"; - TrackDXEString = "trackDXE"; - break; - - default: std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl; break; - } - - InputFileString = "hitresol_ALCARECO_UL.root"; - break; - default: std::cout << "The UL input parameter must be set to 0 (for ALCARECO) or 1 (for UL ALCARECO)." << std::endl; break; + switch (UL) { + case 0: + switch (Unit_Int) { + case 0: + GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO.root"; + HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO.txt"; + CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO.txt"; + DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch"; + HitDXString = "hitDX_OverPitch"; + TrackDXString = "trackDX_OverPitch"; + TrackDXEString = "trackDXE_OverPitch"; + break; + + case 1: + GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO.root"; + HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO.txt"; + CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO.txt"; + DoubleDiffString = "hitDX-trackDX"; + HitDXString = "hitDX"; + TrackDXString = "trackDX"; + TrackDXEString = "trackDXE"; + break; + + default: + std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl; + break; + } + + InputFileString = "hitresol_ALCARECO.root"; + break; + + case 1: + switch (Unit_Int) { + case 0: + GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO_UL.root"; + HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO_UL.txt"; + CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO_UL.txt"; + DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch"; + HitDXString = "hitDX_OverPitch"; + TrackDXString = "trackDX_OverPitch"; + TrackDXEString = "trackDXE_OverPitch"; + break; + + case 1: + GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO_UL.root"; + HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO_UL.txt"; + CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO_UL.txt"; + DoubleDiffString = "hitDX-trackDX"; + HitDXString = "hitDX"; + TrackDXString = "trackDX"; + TrackDXEString = "trackDXE"; + break; + + default: + std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl; + break; + } + + InputFileString = "hitresol_ALCARECO_UL.root"; + break; + default: + std::cout << "The UL input parameter must be set to 0 (for ALCARECO) or 1 (for UL ALCARECO)." << std::endl; + break; } //opening the root file @@ -86,214 +95,287 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& int RegionInt = 0; - if(region == "TIB_L1"){RegionInt = 1;} - else if(region == "TIB_L2"){RegionInt = 2;} - else if(region == "TIB_L3"){RegionInt = 3;} - else if(region == "TIB_L4"){RegionInt = 4;} - else if(region == "Side_TID"){RegionInt = 5;} - else if(region == "Wheel_TID"){RegionInt = 6;} - else if(region == "Ring_TID"){RegionInt = 7;} - else if(region == "TOB_L1"){RegionInt = 8;} - else if(region == "TOB_L2"){RegionInt = 9;} - else if(region == "TOB_L3"){RegionInt = 10;} - else if(region == "TOB_L4"){RegionInt = 11;} - else if(region == "TOB_L5"){RegionInt = 12;} - else if(region == "TOB_L6"){RegionInt = 13;} - else if(region == "Side_TEC"){RegionInt = 14;} - else if(region == "Wheel_TEC"){RegionInt = 15;} - else if(region == "Ring_TEC"){RegionInt = 16;} - else if(region == "TIB_All"){RegionInt = 17;} - else if(region == "TOB_All"){RegionInt = 18;} - else if(region == "TID_All"){RegionInt = 19;} - else if(region == "TEC_All"){RegionInt = 20;} - else if(region == "Pixel_Barrel"){RegionInt = 21;} - else if(region == "Pixel_EndcapDisk"){RegionInt = 22;} - else{std::cout << "Error: The tracker region " << region << " was chosen. Please choose a region out of: TIB L1, TIB L2, TIB L3, TIB L4, Side TID, Wheel TID, Ring TID, TOB L1, TOB L2, TOB L3, TOB L4, TOB L5, TOB L6, Side TEC, Wheel TEC or Ring TEC." << std::endl; return 0;} - - + if (region == "TIB_L1") { + RegionInt = 1; + } else if (region == "TIB_L2") { + RegionInt = 2; + } else if (region == "TIB_L3") { + RegionInt = 3; + } else if (region == "TIB_L4") { + RegionInt = 4; + } else if (region == "Side_TID") { + RegionInt = 5; + } else if (region == "Wheel_TID") { + RegionInt = 6; + } else if (region == "Ring_TID") { + RegionInt = 7; + } else if (region == "TOB_L1") { + RegionInt = 8; + } else if (region == "TOB_L2") { + RegionInt = 9; + } else if (region == "TOB_L3") { + RegionInt = 10; + } else if (region == "TOB_L4") { + RegionInt = 11; + } else if (region == "TOB_L5") { + RegionInt = 12; + } else if (region == "TOB_L6") { + RegionInt = 13; + } else if (region == "Side_TEC") { + RegionInt = 14; + } else if (region == "Wheel_TEC") { + RegionInt = 15; + } else if (region == "Ring_TEC") { + RegionInt = 16; + } else if (region == "TIB_All") { + RegionInt = 17; + } else if (region == "TOB_All") { + RegionInt = 18; + } else if (region == "TID_All") { + RegionInt = 19; + } else if (region == "TEC_All") { + RegionInt = 20; + } else if (region == "Pixel_Barrel") { + RegionInt = 21; + } else if (region == "Pixel_EndcapDisk") { + RegionInt = 22; + } else { + std::cout << "Error: The tracker region " << region + << " was chosen. Please choose a region out of: TIB L1, TIB L2, TIB L3, TIB L4, Side TID, Wheel TID, " + "Ring TID, TOB L1, TOB L2, TOB L3, TOB L4, TOB L5, TOB L6, Side TEC, Wheel TEC or Ring TEC." + << std::endl; + return 0; + } //Lambda function to filter the detID for different layers - auto SubDet_Function{[&RegionInt](const int& detID1_input, const int& detID2_input){ - - bool OutputBool = 0; - - switch(RegionInt){ - - case 1: {OutputBool = (((detID1_input>>25)&0x7) == 3) && ((detID1_input>>14)&0x7) == 1 && - (((detID2_input>>25)&0x7) == 3) && ((detID2_input>>14)&0x7) == 1; //TIB L1 - break;} - - case 2: {OutputBool = (((detID1_input>>25)&0x7) == 3) && (((detID1_input>>14)&0x7) == 2) && - (((detID2_input>>25)&0x7) == 3) && (((detID2_input>>14)&0x7) == 2); //TIB L2 - break;} - - case 3: {OutputBool = (((detID1_input>>25)&0x7) == 3) && (((detID1_input>>14)&0x7) == 3) && - (((detID2_input>>25)&0x7) == 3) && (((detID2_input>>14)&0x7) == 3); //TIB L3 - break;} - - case 4: {OutputBool = (((detID1_input>>25)&0x7) == 3) && (((detID1_input>>14)&0x7) == 4) && - (((detID2_input>>25)&0x7) == 3) && (((detID2_input>>14)&0x7) == 4); //TIB L4 - break;} - - - case 5: {OutputBool = ( (((detID1_input>>13)&0x3) == 1) && (((detID2_input>>13)&0x3) == 1) ) || - ( (((detID1_input>>13)&0x3) == 2) && (((detID2_input>>13)&0x3) == 2) ); //TID Side (1 -> TID-, 2 -> TID+) - - break;} - - case 6: {OutputBool = (((detID1_input>>11)&0x3) == 2) && (((detID2_input>>11)&0x3) == 2); //TID Wheel - - break;} - - case 7: {OutputBool = ( (((detID1_input>>9)&0x3) == 2) && (((detID2_input>>9)&0x3) == 2) ); //TID Ring - - break;} - - - case 8: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 1) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 1); //TOB L1 - break;} - - case 9: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 2) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 2); //TOB L2 - break;} - - case 10: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 3) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 3); //TOB L3 - break;} - - case 11: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 4) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 4); //TOB L4 - break;} - - case 12: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 5) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 5); //TOB L5 - break;} - - case 13: {OutputBool = (((detID1_input>>25)&0x7) == 5) && (((detID1_input>>14)&0x7) == 6) && - (((detID2_input>>25)&0x7) == 5) && (((detID2_input>>14)&0x7) == 6); //TOB L6 - break;} - - - case 14: {OutputBool = ( (((detID1_input>>18)&0x3) == 1) && (((detID2_input>>18)&0x3) == 1) ) || - ( (((detID1_input>>18)&0x3) == 2) && (((detID2_input>>18)&0x3) == 2) ); //Side TEC (1 -> back, 2 -> front) - break;} - - case 15: {OutputBool = (((detID1_input>>14)&0xF) == 4) && (((detID2_input>>14)&0xF) == 4); //Wheel TEC - break;} - - case 16: {OutputBool = (((detID1_input>>5)&0x7) == 3) && (((detID2_input>>5)&0x7) == 3); //Ring TEC - - break;} - - case 17: {OutputBool = ( (((detID1_input>>25)&0x7) == 3) && (((detID2_input>>25)&0x7) == 3) ); //All TIB - - break;} - - case 18: {OutputBool = ( (((detID1_input>>25)&0x7) == 5) && (((detID2_input>>25)&0x7) == 5) ); //All TOB - - break;} - - case 19: {OutputBool = ( (((detID1_input>>13)&0x3) == 1) && (((detID2_input>>13)&0x7) == 1) ) || - ( (((detID1_input>>13)&0x3) == 2) && (((detID2_input>>13)&0x7) == 2) ) || - ( (((detID1_input>>11)&0x3) == 2) && (((detID2_input>>11)&0x3) == 2) ) || - ( (((detID1_input>>9)&0x3) == 2) && (((detID2_input>>9)&0x3) == 2) ) || - ( (((detID1_input>>7)&0x3) == 1) && (((detID2_input>>7)&0x3) == 1) ) || - ( (((detID1_input>>7)&0x3) == 2) && (((detID2_input>>7)&0x3) == 2) ) || - ( (((detID1_input>>2)&0x1F) == 5) && (((detID2_input>>2)&0x1F) == 5) ) || - ( (((detID1_input>>0)&0x3) == 0) && (((detID2_input>>0)&0x3) == 0) ) || - ( (((detID1_input>>0)&0x3) == 1) && (((detID2_input>>0)&0x3) == 1) ) || - ( (((detID1_input>>0)&0x3) == 2) && (((detID2_input>>0)&0x3) == 2) ); //All TID - - break;} - - case 20: { - - OutputBool = ( (((detID1_input>>18)&0x3) == 1) && (((detID2_input>>18)&0x3) == 1) ) || - ( (((detID1_input>>18)&0x3) == 2) && (((detID2_input>>18)&0x3) == 2) ) || - ( (((detID1_input>>14)&0xF) == 4) && (((detID2_input>>14)&0xF) == 4) ) || - ( (((detID1_input>>12)&0x3) == 1) && (((detID2_input>>12)&0x3) == 1) ) || - ( (((detID1_input>>12)&0x3) == 2) && (((detID2_input>>12)&0x3) == 2) ) || - ( (((detID1_input>>8)&0xF) == 4) && (((detID2_input>>8)&0xF) == 4) ) || - ( (((detID1_input>>5)&0x7) == 3) && (((detID2_input>>5)&0x7) == 3) ) || - ( (((detID1_input>>2)&0x7) == 3) && (((detID2_input>>2)&0x7) == 3) ) || - ( (((detID1_input>>0)&0x3) == 1) && (((detID2_input>>0)&0x3) == 1) ) || - ( (((detID1_input>>0)&0x3) == 2) && (((detID2_input>>0)&0x3) == 2) ) || - ( (((detID1_input>>0)&0x3) == 3) && (((detID2_input>>0)&0x3) == 3) ); //All TEC - - break;} - - case 21: {OutputBool = (((detID1_input>>20)&0xF) == 4) && (((detID2_input>>20)&0xF) == 4); //pixel barrel (phase 1) - break;} - - case 22: {OutputBool = (((detID1_input>>18)&0xF) == 4) && (((detID2_input>>18)&0xF) == 4); //pixel endcap disk (phase 1) - break;} - - } - - return OutputBool; - + auto SubDet_Function{[&RegionInt](const int& detID1_input, const int& detID2_input) { + bool OutputBool = 0; + + switch (RegionInt) { + case 1: { + OutputBool = (((detID1_input >> 25) & 0x7) == 3) && ((detID1_input >> 14) & 0x7) == 1 && + (((detID2_input >> 25) & 0x7) == 3) && ((detID2_input >> 14) & 0x7) == 1; //TIB L1 + break; + } + + case 2: { + OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 2) && + (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 2); //TIB L2 + break; + } + + case 3: { + OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 3) && + (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 3); //TIB L3 + break; + } + + case 4: { + OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 4) && + (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 4); //TIB L4 + break; + } + + case 5: { + OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x3) == 1)) || + ((((detID1_input >> 13) & 0x3) == 2) && + (((detID2_input >> 13) & 0x3) == 2)); //TID Side (1 -> TID-, 2 -> TID+) + + break; + } + + case 6: { + OutputBool = (((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2); //TID Wheel + + break; + } + + case 7: { + OutputBool = ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2)); //TID Ring + + break; + } + + case 8: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 1) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 1); //TOB L1 + break; + } + + case 9: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 2) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 2); //TOB L2 + break; + } + + case 10: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 3) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 3); //TOB L3 + break; + } + + case 11: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 4) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 4); //TOB L4 + break; + } + + case 12: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 5) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 5); //TOB L5 + break; + } + + case 13: { + OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 6) && + (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 6); //TOB L6 + break; + } + + case 14: { + OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) || + ((((detID1_input >> 18) & 0x3) == 2) && + (((detID2_input >> 18) & 0x3) == 2)); //Side TEC (1 -> back, 2 -> front) + break; + } + + case 15: { + OutputBool = (((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4); //Wheel TEC + break; + } + + case 16: { + OutputBool = (((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3); //Ring TEC + + break; + } + + case 17: { + OutputBool = ((((detID1_input >> 25) & 0x7) == 3) && (((detID2_input >> 25) & 0x7) == 3)); //All TIB + + break; + } + + case 18: { + OutputBool = ((((detID1_input >> 25) & 0x7) == 5) && (((detID2_input >> 25) & 0x7) == 5)); //All TOB + + break; + } + + case 19: { + OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x7) == 1)) || + ((((detID1_input >> 13) & 0x3) == 2) && (((detID2_input >> 13) & 0x7) == 2)) || + ((((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2)) || + ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2)) || + ((((detID1_input >> 7) & 0x3) == 1) && (((detID2_input >> 7) & 0x3) == 1)) || + ((((detID1_input >> 7) & 0x3) == 2) && (((detID2_input >> 7) & 0x3) == 2)) || + ((((detID1_input >> 2) & 0x1F) == 5) && (((detID2_input >> 2) & 0x1F) == 5)) || + ((((detID1_input >> 0) & 0x3) == 0) && (((detID2_input >> 0) & 0x3) == 0)) || + ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) || + ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2)); //All TID + + break; + } + + case 20: { + OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) || + ((((detID1_input >> 18) & 0x3) == 2) && (((detID2_input >> 18) & 0x3) == 2)) || + ((((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4)) || + ((((detID1_input >> 12) & 0x3) == 1) && (((detID2_input >> 12) & 0x3) == 1)) || + ((((detID1_input >> 12) & 0x3) == 2) && (((detID2_input >> 12) & 0x3) == 2)) || + ((((detID1_input >> 8) & 0xF) == 4) && (((detID2_input >> 8) & 0xF) == 4)) || + ((((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3)) || + ((((detID1_input >> 2) & 0x7) == 3) && (((detID2_input >> 2) & 0x7) == 3)) || + ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) || + ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2)) || + ((((detID1_input >> 0) & 0x3) == 3) && (((detID2_input >> 0) & 0x3) == 3)); //All TEC + + break; + } + + case 21: { + OutputBool = + (((detID1_input >> 20) & 0xF) == 4) && (((detID2_input >> 20) & 0xF) == 4); //pixel barrel (phase 1) + break; + } + + case 22: { + OutputBool = + (((detID1_input >> 18) & 0xF) == 4) && (((detID2_input >> 18) & 0xF) == 4); //pixel endcap disk (phase 1) + break; + } + } + + return OutputBool; }}; //Function for expressing the hit resolution in either micrometres or pitch units. - auto Pitch_Function{[&Unit_Int](const float& pitch, const float& input){ - - float InputOverPitch = input/pitch; - return InputOverPitch; - + auto Pitch_Function{[&Unit_Int](const float& pitch, const float& input) { + float InputOverPitch = input / pitch; + return InputOverPitch; }}; - //Defining columns needed for the unit conversion into pitch units, and applying the filter for the subdetector auto dataframe = d.Define("hitDX_OverPitch", Pitch_Function, {"pitch1", "hitDX"}) - .Define("trackDX_OverPitch", Pitch_Function, {"pitch1", "trackDX"}) - .Define("trackDXE_OverPitch", Pitch_Function, {"pitch1", "trackDXE"}) - .Filter(SubDet_Function, {"detID1", "detID2"}, "Subdetector filter"); + .Define("trackDX_OverPitch", Pitch_Function, {"pitch1", "trackDX"}) + .Define("trackDXE_OverPitch", Pitch_Function, {"pitch1", "trackDXE"}) + .Filter(SubDet_Function, {"detID1", "detID2"}, "Subdetector filter"); //Implementing selection criteria that were not implemented in HitResol.cc - auto PairPathCriteriaFunction{[&RegionInt](const float& pairPath_input){ - - if((RegionInt > 0 && RegionInt < 5) || (RegionInt > 7 || RegionInt < 13) || (RegionInt == 17) || (RegionInt == 18)){return abs(pairPath_input) < 7;} //for TIB and TOB - else if(RegionInt == 21 || RegionInt == 22){return abs(pairPath_input) < 2;} //for pixels - else{return abs(pairPath_input) < 20;}//for everything else (max value is 15cm so this will return all events anyway) + auto PairPathCriteriaFunction{[&RegionInt](const float& pairPath_input) { + if ((RegionInt > 0 && RegionInt < 5) || (RegionInt > 7 || RegionInt < 13) || (RegionInt == 17) || + (RegionInt == 18)) { + return abs(pairPath_input) < 7; + } //for TIB and TOB + else if (RegionInt == 21 || RegionInt == 22) { + return abs(pairPath_input) < 2; + } //for pixels + else { + return abs(pairPath_input) < 20; + } //for everything else (max value is 15cm so this will return all events anyway) }}; - auto MomentaFunction{[&RegionInt](const float& momentum_input){ - - if(RegionInt == 21 || RegionInt == 22){return momentum_input > 5;} //pixels - else{return momentum_input > 15;} //strips + auto MomentaFunction{[&RegionInt](const float& momentum_input) { + if (RegionInt == 21 || RegionInt == 22) { + return momentum_input > 5; + } //pixels + else { + return momentum_input > 15; + } //strips }}; - auto dataframe_filtered = dataframe.Filter(PairPathCriteriaFunction, {"pairPath"}, "Pair path criterion filter") - .Filter(MomentaFunction, {"momentum"}, "Momentum criterion filter") - .Filter("trackChi2 > 0.001", "chi2 criterion filter") - .Filter("numHits > 6", "numHits filter") - .Filter("trackDXE < 0.0025", "trackDXE filter") - .Filter("(clusterW1 == clusterW2) && clusterW1 <= 4 && clusterW2 <= 4", "cluster filter"); + auto dataframe_filtered = + dataframe.Filter(PairPathCriteriaFunction, {"pairPath"}, "Pair path criterion filter") + .Filter(MomentaFunction, {"momentum"}, "Momentum criterion filter") + .Filter("trackChi2 > 0.001", "chi2 criterion filter") + .Filter("numHits > 6", "numHits filter") + .Filter("trackDXE < 0.0025", "trackDXE filter") + .Filter("(clusterW1 == clusterW2) && clusterW1 <= 4 && clusterW2 <= 4", "cluster filter"); //Creating histograms for the difference between the two hit positions, the difference between the two predicted positions and for the double difference //hitDX = the difference in the hit positions for the pair - //trackDX = the difference in the track positions for the pair + //trackDX = the difference in the track positions for the pair auto HistoName_DoubleDiff = "DoubleDifference_" + region; auto HistoName_HitDX = "HitDX_" + region; - auto HistoName_TrackDX = "TrackDX_" + region; + auto HistoName_TrackDX = "TrackDX_" + region; auto HistoName_TrackDXE = "TrackDXE_" + region; auto HistoName_ClusterW1 = "ClusterW1_" + region; auto HistoName_ClusterW2 = "ClusterW2_" + region; - - auto h_DoubleDifference = dataframe_filtered.Define(HistoName_DoubleDiff, DoubleDiffString).Histo1D({HistoName_DoubleDiff.c_str(), HistoName_DoubleDiff.c_str(), 40, -0.5, 0.5}, HistoName_DoubleDiff); + auto h_DoubleDifference = + dataframe_filtered.Define(HistoName_DoubleDiff, DoubleDiffString) + .Histo1D({HistoName_DoubleDiff.c_str(), HistoName_DoubleDiff.c_str(), 40, -0.5, 0.5}, HistoName_DoubleDiff); auto h_hitDX = dataframe_filtered.Define(HistoName_HitDX, HitDXString).Histo1D(HistoName_HitDX); auto h_trackDX = dataframe_filtered.Define(HistoName_TrackDX, TrackDXString).Histo1D(HistoName_TrackDX); auto h_trackDXE = dataframe_filtered.Define(HistoName_TrackDXE, TrackDXEString).Histo1D(HistoName_TrackDXE); - + auto h_clusterW1 = dataframe_filtered.Define(HistoName_ClusterW1, ClusterW1String).Histo1D(HistoName_ClusterW1); auto h_clusterW2 = dataframe_filtered.Define(HistoName_ClusterW2, ClusterW2String).Histo1D(HistoName_ClusterW2); //Applying gaussian fits, taking the resolutions and squaring them h_DoubleDifference->Fit("gaus"); - + auto double_diff_StdDev = h_DoubleDifference->GetStdDev(); auto hitDX_StdDev = h_hitDX->GetStdDev(); auto trackDX_StdDev = h_trackDX->GetStdDev(); @@ -301,7 +383,7 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& auto sigma2_MeasMinusPred = pow(double_diff_StdDev, 2); auto sigma2_Meas = pow(hitDX_StdDev, 2); - auto sigma2_Pred = pow(trackDX_StdDev, 2); + auto sigma2_Pred = pow(trackDX_StdDev, 2); auto sigma2_PredError = pow(trackDXE_Mean, 2); DoubleDifferenceVector.push_back(sigma2_MeasMinusPred); @@ -309,9 +391,8 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& TrackDXVector.push_back(sigma2_Pred); TrackDXEVector.push_back(sigma2_PredError); - //Saving the histograms with gaussian fits applied to an output root file - TFile * output = new TFile(GaussianFitsFileName.c_str(), "UPDATE"); + TFile* output = new TFile(GaussianFitsFileName.c_str(), "UPDATE"); h_DoubleDifference->Write(); h_hitDX->Write(); @@ -319,16 +400,16 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& h_trackDXE->Write(); h_clusterW1->Write(); h_clusterW2->Write(); - + output->Close(); - + //Calculating the hit resolution; auto numerator = sigma2_MeasMinusPred - sigma2_PredError; - auto HitResolution = sqrt( numerator/2 ); + auto HitResolution = sqrt(numerator / 2); HitResolutionVector.push_back(HitResolution); - //Printing the resolution + //Printing the resolution std::cout << '\n' << std::endl; std::cout << "The hit resolution for tracker region " << region << " is: " << HitResolution << std::endl; std::cout << '\n' << std::endl; @@ -339,45 +420,42 @@ void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& CutFlowReport.open(CutFlowReportString.c_str()); - for(auto&& cutInfo: allCutsReport){ - CutFlowReport << cutInfo.GetName() << '\t' << cutInfo.GetAll() << '\t' << cutInfo.GetPass() << '\t' << cutInfo.GetEff() << " %" << std::endl; + for (auto&& cutInfo : allCutsReport) { + CutFlowReport << cutInfo.GetName() << '\t' << cutInfo.GetAll() << '\t' << cutInfo.GetPass() << '\t' + << cutInfo.GetEff() << " %" << std::endl; } - } - - - -void Resolutions(){ - +void Resolutions() { int UnitInteger = 0; int ULInteger = 0; + vector LayerNames = {"TIB_L1", "TIB_L2", "TIB_L3", "TIB_L4", "Side_TID", "Wheel_TID", + "Ring_TID", "TOB_L1", "TOB_L2", "TOB_L3", "TOB_L4", "TOB_L5", + "TOB_L6", "Side_TEC", "Wheel_TEC", "Ring_TEC", "TIB_All", "TOB_All", + "TID_All", "TEC_All", "Pixel_Barrel", "Pixel_EndcapDisk"}; - vector LayerNames = {"TIB_L1", "TIB_L2", "TIB_L3", "TIB_L4", - "Side_TID", "Wheel_TID", "Ring_TID", "TOB_L1", - "TOB_L2", "TOB_L3", "TOB_L4", "TOB_L5", - "TOB_L6", "Side_TEC", "Wheel_TEC", "Ring_TEC", - "TIB_All", "TOB_All", "TID_All", "TEC_All", - "Pixel_Barrel", "Pixel_EndcapDisk"}; - - - for(int i = 0; i < LayerNames.size(); i++){ResolutionsCalculator(LayerNames.at(i), UnitInteger, ULInteger);} + for (int i = 0; i < LayerNames.size(); i++) { + ResolutionsCalculator(LayerNames.at(i), UnitInteger, ULInteger); + } std::ofstream HitResoTextFile; HitResoTextFile.open(HitResoFileName); auto Width = 28; - HitResoTextFile << std::right << "Layer " << std::setw(Width) << " Resolution " << std::setw(Width) << " sigma2_HitDX " << std::setw(Width) << " sigma2_trackDX " << std::setw(Width) << " sigma2_trackDXE " << std::setw(Width) << " sigma2_DoubleDifference " << std::endl; - - for(int i = 0; i < HitResolutionVector.size(); i++){ - HitResoTextFile << std::right << LayerNames.at(i) << std::setw(Width) << HitResolutionVector.at(i) << std::setw(Width) << HitDXVector.at(i) << std::setw(Width) << TrackDXVector.at(i) << std::setw(Width) << TrackDXEVector.at(i) << std::setw(Width) << DoubleDifferenceVector.at(i) << std::endl; + HitResoTextFile << std::right << "Layer " << std::setw(Width) << " Resolution " << std::setw(Width) + << " sigma2_HitDX " << std::setw(Width) << " sigma2_trackDX " << std::setw(Width) + << " sigma2_trackDXE " << std::setw(Width) << " sigma2_DoubleDifference " << std::endl; + for (int i = 0; i < HitResolutionVector.size(); i++) { + HitResoTextFile << std::right << LayerNames.at(i) << std::setw(Width) << HitResolutionVector.at(i) + << std::setw(Width) << HitDXVector.at(i) << std::setw(Width) << TrackDXVector.at(i) + << std::setw(Width) << TrackDXEVector.at(i) << std::setw(Width) << DoubleDifferenceVector.at(i) + << std::endl; } - - system("mkdir HitResolutionValues; mkdir GaussianFits; mkdir CutFlowReports; mv CutFlowReport_* CutFlowReports/; mv HitResolutionValues_* HitResolutionValues/; mv GaussianFits_* GaussianFits/;"); - - + system( + "mkdir HitResolutionValues; mkdir GaussianFits; mkdir CutFlowReports; mv CutFlowReport_* CutFlowReports/; mv " + "HitResolutionValues_* HitResolutionValues/; mv GaussianFits_* GaussianFits/;"); } diff --git a/CalibTracker/SiStripHitResolution/plugins/BuildFile.xml b/CalibTracker/SiStripHitResolution/plugins/BuildFile.xml index 329f10240f767..90ac2222d56a9 100644 --- a/CalibTracker/SiStripHitResolution/plugins/BuildFile.xml +++ b/CalibTracker/SiStripHitResolution/plugins/BuildFile.xml @@ -1,33 +1,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/CalibTracker/SiStripHitResolution/plugins/SiStripCPEAnalyzer.cc b/CalibTracker/SiStripHitResolution/plugins/SiStripCPEAnalyzer.cc new file mode 100644 index 0000000000000..229e0738a050a --- /dev/null +++ b/CalibTracker/SiStripHitResolution/plugins/SiStripCPEAnalyzer.cc @@ -0,0 +1,346 @@ +// -*- C++ -*- +// +// Package: CalibTracker/SiStripHitResolution +// Class: SiStripCPEAnalyzer +// +/**\class SiStripCPEAnalyzer SiStripCPEAnalyzer.cc CalibTracker/SiStripHitResolution/plugins/SiStripCPEAnalyzer.cc + + Description: CPE analysis (resolution param, etc.) + + Implementation: + Used for Strip Tracker DPG studies +*/ +// +// Original Author: Christophe Delaere +// Created: Thu, 12 Sep 2019 13:51:00 GMT +// +// + +// system include files +#include +#include +#include + +// user include files +#include "CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" +#include "RecoLocalTracker/Records/interface/TkStripCPERecord.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h" +#include "TrackingTools/PatternTools/interface/Trajectory.h" + +#include "TH1.h" +#include "TTree.h" + +// +// class declaration +// + +using reco::TrackCollection; +typedef std::vector TrajectoryCollection; + +struct CPEBranch { + float x, y, z; + float distance, mdistance, shift, offsetA, offsetB, angle; + float x1r, x2r, x1m, x2m, y1m, y2m; +}; + +struct TrackBranch { + float px, py, pz, pt, eta, phi, charge; +}; + +struct GeometryBranch { + int subdet, moduleGeometry, stereo, layer, side, ring; + float pitch; + int detid; +}; + +struct ClusterBranch { + int strips[11]; + int size, firstStrip; + float barycenter; +}; + +class SiStripCPEAnalyzer : public edm::one::EDAnalyzer { +public: + explicit SiStripCPEAnalyzer(const edm::ParameterSet&); + ~SiStripCPEAnalyzer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void analyze(const edm::Event&, const edm::EventSetup&) override; + static bool goodMeasurement(TrajectoryMeasurement const& m); + + // ----------member data --------------------------- + + // tokens for event data + const edm::EDGetTokenT tracksToken_; //used to select what tracks to read from configuration file + const edm::EDGetTokenT + trajsToken_; //used to select what trajectories to read from configuration file + const edm::EDGetTokenT tjTagToken_; //association map between tracks and trajectories + const edm::EDGetTokenT > clusToken_; //clusters + + // tokens for ES data + const edm::ESInputTag cpeTag_; + const edm::ESGetToken cpeToken_; + const edm::ESGetToken topoToken_; + const edm::ESGetToken geomToken_; + + const StripClusterParameterEstimator* parameterestimator_; //CPE + + // to fill the tree + TTree* tree_; + CPEBranch treeBranches_; + TrackBranch trackBranches_; + GeometryBranch geom1Branches_; + GeometryBranch geom2Branches_; + ClusterBranch cluster1Branches_; + ClusterBranch cluster2Branches_; +}; + +// +// constructors and destructor +// +SiStripCPEAnalyzer::SiStripCPEAnalyzer(const edm::ParameterSet& iConfig) + : tracksToken_(consumes(iConfig.getUntrackedParameter("tracks"))), + trajsToken_(consumes(iConfig.getUntrackedParameter("trajectories"))), + tjTagToken_( + consumes(iConfig.getUntrackedParameter("association"))), + clusToken_( + consumes >(iConfig.getUntrackedParameter("clusters"))), + cpeTag_(iConfig.getParameter("StripCPE")), + cpeToken_(esConsumes(cpeTag_)), + topoToken_(esConsumes()), + geomToken_(esConsumes()) { + //now do what ever initialization is needed + usesResource("TFileService"); + edm::Service fs; + //histo = fs->make("charge" , "Charges" , 3 , -1 , 2 ); + tree_ = fs->make("CPEanalysis", "CPE analysis tree"); + tree_->Branch( + "Overlaps", &treeBranches_, "x:y:z:distance:mdistance:shift:offsetA:offsetB:angle:x1r:x2r:x1m:x2m:y1m:y2m"); + tree_->Branch("Tracks", &trackBranches_, "px:py:pz:pt:eta:phi:charge"); + tree_->Branch("Cluster1", &cluster1Branches_, "strips[11]/I:size/I:firstStrip/I:barycenter/F"); + tree_->Branch("Cluster2", &cluster2Branches_, "strips[11]/I:size/I:firstStrip/I:barycenter/F"); + tree_->Branch("Geom1", &geom1Branches_, "subdet/I:moduleGeometry/I:stereo/I:layer/I:side/I:ring/I:pitch/F:detid/I"); + tree_->Branch("Geom2", &geom2Branches_, "subdet/I:moduleGeometry/I:stereo/I:layer/I:side/I:ring/I:pitch/F:detid/I"); +} + +// destructor +SiStripCPEAnalyzer::~SiStripCPEAnalyzer() {} + +// +// member functions +// +// ------------ method called for each event ------------ +void SiStripCPEAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + using namespace edm; + + // load the CPE, geometry and topology + parameterestimator_ = &iSetup.getData(cpeToken_); //CPE + const TrackerGeometry* tracker = &iSetup.getData(geomToken_); + const TrackerTopology* const tTopo = &iSetup.getData(topoToken_); + + // prepare the output + std::vector hitpairs; + + // loop on trajectories and tracks + //for(const auto& tt : iEvent.get(tjTagToken_) ) { + edm::Handle trajTrackAssociations; + iEvent.getByToken(tjTagToken_, trajTrackAssociations); + for (const auto& tt : *trajTrackAssociations) { + auto& traj = *(tt.key); + auto& track = *(tt.val); + // track quantities + trackBranches_.px = track.px(); + trackBranches_.py = track.py(); + trackBranches_.pz = track.pz(); + trackBranches_.pt = track.pt(); + trackBranches_.eta = track.eta(); + trackBranches_.phi = track.phi(); + trackBranches_.charge = track.charge(); + // loop on measurements + for (auto it = traj.measurements().begin(); it != traj.measurements().end(); ++it) { + auto& meas = *it; + + // only OT measurements on valid hits (not glued det) + if (!goodMeasurement(meas)) + continue; + + // restrict the search for compatible hits in the same layer (measurements are sorted by layer) + auto layerRangeEnd = it + 1; + for (; layerRangeEnd < traj.measurements().end(); ++layerRangeEnd) { + if (layerRangeEnd->layer()->seqNum() != meas.layer()->seqNum()) + break; + } + + // now look for a matching hit in that range. + auto meas2it = std::find_if(it + 1, layerRangeEnd, [meas](const TrajectoryMeasurement& m) -> bool { + return goodMeasurement(m) && (m.recHit()->rawId() & 0x3) == (meas.recHit()->rawId() & 0x3); + }); + + // check if we found something, build a pair object and add to vector + if (meas2it != layerRangeEnd) { + auto& meas2 = *meas2it; + hitpairs.push_back(SiStripOverlapHit(meas, meas2)); + } + } + } + + // load clusters + Handle > clusters; + iEvent.getByToken(clusToken_, clusters); + + // At this stage we will have a vector of pairs of measurement. Fill a ntuple and some histograms. + for (const auto& pair : hitpairs) { + //TODO basically everything below is done twice. -> can be factorized. + + // store generic information about the pair + treeBranches_.x = pair.position().x(); + treeBranches_.y = pair.position().y(); + treeBranches_.z = pair.position().z(); + treeBranches_.distance = pair.distance(false); + treeBranches_.mdistance = pair.distance(true); + treeBranches_.shift = pair.shift(); + treeBranches_.offsetA = pair.offset(0); + treeBranches_.offsetB = pair.offset(1); + treeBranches_.angle = pair.getTrackLocalAngle(0); + treeBranches_.x1r = pair.hitA()->localPosition().x(); + treeBranches_.x1m = pair.trajectoryStateOnSurface(0, false).localPosition().x(); + treeBranches_.y1m = pair.trajectoryStateOnSurface(0, false).localPosition().y(); + treeBranches_.x2r = pair.hitB()->localPosition().x(); + treeBranches_.x2m = pair.trajectoryStateOnSurface(1, false).localPosition().x(); + treeBranches_.y2m = pair.trajectoryStateOnSurface(1, false).localPosition().y(); + + // load the clusters for the detectors + auto detset1 = (*clusters)[pair.hitA()->rawId()]; + auto detset2 = (*clusters)[pair.hitB()->rawId()]; + + // look for the proper cluster + const GeomDetUnit* du; + du = tracker->idToDetUnit(pair.hitA()->rawId()); + auto cluster1 = std::min_element( + detset1.begin(), detset1.end(), [du, this](SiStripCluster const& cl1, SiStripCluster const& cl2) { + return (fabs(parameterestimator_->localParameters(cl1, *du).first.x() - treeBranches_.x1r) < + fabs(parameterestimator_->localParameters(cl2, *du).first.x() - treeBranches_.x1r)); + }); + du = tracker->idToDetUnit(pair.hitB()->rawId()); + auto cluster2 = std::min_element( + detset2.begin(), detset2.end(), [du, this](SiStripCluster const& cl1, SiStripCluster const& cl2) { + return (fabs(parameterestimator_->localParameters(cl1, *du).first.x() - treeBranches_.x2r) < + fabs(parameterestimator_->localParameters(cl2, *du).first.x() - treeBranches_.x2r)); + }); + + // store the amplitudes centered around the maximum + auto amplitudes1 = cluster1->amplitudes(); + auto amplitudes2 = cluster2->amplitudes(); + auto max1 = std::max_element(amplitudes1.begin(), amplitudes1.end()); + auto max2 = std::max_element(amplitudes2.begin(), amplitudes2.end()); + for (unsigned int i = 0; i < 11; ++i) { + cluster1Branches_.strips[i] = 0.; + cluster2Branches_.strips[i] = 0.; + } + cluster1Branches_.size = amplitudes1.size(); + cluster2Branches_.size = amplitudes2.size(); + cluster1Branches_.firstStrip = cluster1->firstStrip(); + cluster1Branches_.barycenter = cluster1->barycenter(); + cluster2Branches_.firstStrip = cluster2->firstStrip(); + cluster2Branches_.barycenter = cluster2->barycenter(); + + int cnt = 0; + int offset = 5 - (max1 - amplitudes1.begin()); + for (auto& s : amplitudes1) { + if ((offset + cnt) >= 0 && (offset + cnt) < 11) { + cluster1Branches_.strips[offset + cnt] = s; + } + cnt++; + } + cnt = 0; + offset = 5 - (max2 - amplitudes2.begin()); + for (auto& s : amplitudes2) { + if ((offset + cnt) >= 0 && (offset + cnt) < 11) { + cluster2Branches_.strips[offset + cnt] = s; + } + cnt++; + } + + // store information about the geometry (for both sensors) + DetId detid1 = pair.hitA()->geographicalId(); + DetId detid2 = pair.hitB()->geographicalId(); + geom1Branches_.detid = detid1.rawId(); + geom2Branches_.detid = detid2.rawId(); + geom1Branches_.subdet = detid1.subdetId(); + geom2Branches_.subdet = detid2.subdetId(); + //geom1Branches_.moduleGeometry = tTopo->moduleGeometry(detid1); + //geom2Branches_.moduleGeometry = tTopo->moduleGeometry(detid2); + geom1Branches_.stereo = tTopo->isStereo(detid1); + geom2Branches_.stereo = tTopo->isStereo(detid2); + geom1Branches_.layer = tTopo->layer(detid1); + geom2Branches_.layer = tTopo->layer(detid2); + geom1Branches_.side = tTopo->side(detid1); + geom2Branches_.side = tTopo->side(detid2); + if (geom1Branches_.subdet == StripSubdetector::TID) { + geom1Branches_.ring = tTopo->tidRing(detid1); + } else if (geom1Branches_.subdet == StripSubdetector::TEC) { + geom1Branches_.ring = tTopo->tecRing(detid1); + } else { + geom1Branches_.ring = 0; + } + if (geom2Branches_.subdet == StripSubdetector::TID) { + geom2Branches_.ring = tTopo->tidRing(detid2); + } else if (geom2Branches_.subdet == StripSubdetector::TEC) { + geom2Branches_.ring = tTopo->tecRing(detid2); + } else { + geom2Branches_.ring = 0; + } + + const StripGeomDetUnit* theStripDet; + theStripDet = dynamic_cast(tracker->idToDetUnit(pair.hitA()->rawId())); + geom1Branches_.pitch = + theStripDet->specificTopology().localPitch(pair.trajectoryStateOnSurface(0, false).localPosition()); + theStripDet = dynamic_cast(tracker->idToDetUnit(pair.hitB()->rawId())); + geom2Branches_.pitch = + theStripDet->specificTopology().localPitch(pair.trajectoryStateOnSurface(1, false).localPosition()); + + //fill the tree. + tree_->Fill(); // we fill one entry per overlap, to track info is multiplicated + } +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void SiStripCPEAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.addUntracked("tracks", edm::InputTag("generalTracks")); + desc.addUntracked("trajectories", edm::InputTag("generalTracks")); + desc.addUntracked("association", edm::InputTag("generalTracks")); + desc.addUntracked("clusters", edm::InputTag("siStripClusters")); + desc.add("StripCPE", edm::ESInputTag("StripCPEfromTrackAngleESProducer", "StripCPEfromTrackAngle")); + descriptions.addWithDefaultLabel(desc); +} + +bool SiStripCPEAnalyzer::goodMeasurement(TrajectoryMeasurement const& m) { + return m.recHit()->isValid() && // valid + m.recHit()->geographicalId().subdetId() > 2 && // not IT + (m.recHit()->rawId() & 0x3) != 0 && // not glued DetLayer + m.recHit()->getType() == 0; // only valid hits (redundant?) +} + +//define this as a plug-in +DEFINE_FWK_MODULE(SiStripCPEAnalyzer); diff --git a/CalibTracker/SiStripHitResolution/plugins/SiStripHitResolFromCalibTree.cc b/CalibTracker/SiStripHitResolution/plugins/SiStripHitResolFromCalibTree.cc old mode 100755 new mode 100644 index abdf495397683..1ab682d8e6f97 --- a/CalibTracker/SiStripHitResolution/plugins/SiStripHitResolFromCalibTree.cc +++ b/CalibTracker/SiStripHitResolution/plugins/SiStripHitResolFromCalibTree.cc @@ -6,58 +6,53 @@ #include #include -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/ParameterSet/interface/FileInPath.h" - -#include "UserCode/SiStripHitResolution/interface/HitResol.h" -#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h" +#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" +#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h" +#include "CalibTracker/Records/interface/SiStripQualityRcd.h" +#include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h" +#include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h" +#include "CalibTracker/SiStripHitResolution/interface/HitResol.h" +#include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h" +#include "CommonTools/TrackerMap/interface/TrackerMap.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" #include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/DetId/interface/DetIdCollection.h" +#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "DataFormats/GeometryVector/interface/GlobalVector.h" #include "DataFormats/GeometryVector/interface/LocalVector.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "DataFormats/TrackReco/interface/DeDxData.h" #include "DataFormats/TrackReco/interface/Track.h" -#include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/TrackReco/interface/TrackExtra.h" -#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDetType.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" -#include "TrackingTools/Records/interface/TransientRecHitRecord.h" -#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" -#include "UserCode/SiStripHitResolution/interface/TrajectoryAtInvalidHit.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" -#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h" -#include "DataFormats/TrackReco/interface/DeDxData.h" -#include "DataFormats/DetId/interface/DetIdCollection.h" -#include "TrackingTools/DetLayers/interface/DetLayer.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" - #include "RecoTracker/Record/interface/CkfComponentsRecord.h" -#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h" -#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h" -#include "CalibTracker/Records/interface/SiStripQualityRcd.h" -#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" -#include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h" -#include "DataFormats/Common/interface/DetSetVector.h" -#include "DataFormats/Common/interface/DetSetVectorNew.h" -#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" - -#include "DataFormats/MuonReco/interface/Muon.h" -#include "DataFormats/MuonReco/interface/MuonFwd.h" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" -#include "CommonTools/TrackerMap/interface/TrackerMap.h" -#include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" #include "TFile.h" #include "TCanvas.h" @@ -81,7 +76,7 @@ using namespace edm; using namespace reco; using namespace std; -struct hit{ +struct hit { double x; double y; double z; @@ -89,125 +84,133 @@ struct hit{ }; class SiStripHitResolFromCalibTree : public ConditionDBWriter { - public: - explicit SiStripHitResolFromCalibTree(const edm::ParameterSet&); - ~SiStripHitResolFromCalibTree() override; - - private: - void algoBeginJob(const edm::EventSetup&) override; - void algoEndJob() override; - void algoAnalyze(const edm::Event& e, const edm::EventSetup& c) override; - void SetBadComponents(int i, int component,SiStripQuality::BadComponent& BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4]); - void makeTKMap(bool autoTagging); - void makeHotColdMaps(); - void makeSQLite(); - void totalStatistics(); - void makeSummary(); - void makeSummaryVsBx(); - void ComputeEff(vector< TH1F* > &vhfound, vector< TH1F* > &vhtotal, string name); - void makeSummaryVsLumi(); - void makeSummaryVsCM(); - TString GetLayerName(Long_t k); - TString GetLayerSideName(Long_t k); - float calcPhi(float x, float y); - - edm::Service fs; - SiStripDetInfoFileReader* reader; - edm::FileInPath FileInPath_; - SiStripQuality* quality_; - std::unique_ptr getNewObject() override; - - TTree* CalibTree; - vector CalibTreeFilenames; - float threshold; - unsigned int nModsMin; - unsigned int doSummary; - string _badModulesFile; - bool _autoIneffModTagging; - unsigned int _clusterMatchingMethod; - float _ResXSig; - float _clusterTrajDist; - float _stripsApvEdge; - bool _useOnlyHighPurityTracks; - unsigned int _bunchx; - unsigned int _spaceBetweenTrains; - bool _useCM; - bool _showEndcapSides; - bool _showRings; - bool _showTOB6TEC9; - bool _showOnlyGoodModules; - float _tkMapMin; - float _effPlotMin; - TString _title; - - unsigned int nTEClayers; - - TH1F* bxHisto; - TH1F* instLumiHisto; - TH1F* PUHisto; - - // for association of informations of the hitEff tree and the event infos tree - map< pair< unsigned int, unsigned int>, array > eventInfos; - - vector hits[23]; - vector HotColdMaps; - map< unsigned int, pair< unsigned int, unsigned int> > modCounter[23]; - TrackerMap *tkmap; - TrackerMap *tkmapbad; - TrackerMap *tkmapeff; - TrackerMap *tkmapnum; - TrackerMap *tkmapden; - long layerfound[23]; - long layertotal[23]; - map< unsigned int, vector > layerfound_perBx; - map< unsigned int, vector > layertotal_perBx; - vector< TH1F* > layerfound_vsLumi; - vector< TH1F* > layertotal_vsLumi; - vector< TH1F* > layerfound_vsPU; - vector< TH1F* > layertotal_vsPU; - vector< TH1F* > layerfound_vsCM; - vector< TH1F* > layertotal_vsCM; - int goodlayertotal[35]; - int goodlayerfound[35]; - int alllayertotal[35]; - int alllayerfound[35]; - map< unsigned int, double > BadModules; +public: + explicit SiStripHitResolFromCalibTree(const edm::ParameterSet&); + ~SiStripHitResolFromCalibTree() override; + +private: + void algoBeginJob(const edm::EventSetup&) override; + void algoEndJob() override; + void algoAnalyze(const edm::Event& e, const edm::EventSetup& c) override; + void SetBadComponents(int i, + int component, + SiStripQuality::BadComponent& BC, + std::stringstream ssV[4][19], + int NBadComponent[4][19][4]); + void makeTKMap(bool autoTagging); + void makeHotColdMaps(); + void makeSQLite(); + void totalStatistics(); + void makeSummary(); + void makeSummaryVsBx(); + void ComputeEff(vector& vhfound, vector& vhtotal, string name); + void makeSummaryVsLumi(); + void makeSummaryVsCM(); + TString GetLayerName(Long_t k); + TString GetLayerSideName(Long_t k); + float calcPhi(float x, float y); + + edm::Service fs; + SiStripDetInfo _detInfo; + edm::FileInPath FileInPath_; + SiStripQuality* quality_; + std::unique_ptr getNewObject() override; + + TTree* CalibTree; + vector CalibTreeFilenames; + float threshold; + unsigned int nModsMin; + unsigned int doSummary; + string _badModulesFile; + bool _autoIneffModTagging; + unsigned int _clusterMatchingMethod; + float _ResXSig; + float _clusterTrajDist; + float _stripsApvEdge; + bool _useOnlyHighPurityTracks; + unsigned int _bunchx; + unsigned int _spaceBetweenTrains; + bool _useCM; + bool _showEndcapSides; + bool _showRings; + bool _showTOB6TEC9; + bool _showOnlyGoodModules; + float _tkMapMin; + float _effPlotMin; + TString _title; + + edm::ESGetToken _tkGeomToken; + edm::ESGetToken _tTopoToken; + + unsigned int nTEClayers; + + TH1F* bxHisto; + TH1F* instLumiHisto; + TH1F* PUHisto; + + // for association of informations of the hitEff tree and the event infos tree + map, array > eventInfos; + + vector hits[23]; + vector HotColdMaps; + map > modCounter[23]; + TrackerMap* tkmap; + TrackerMap* tkmapbad; + TrackerMap* tkmapeff; + TrackerMap* tkmapnum; + TrackerMap* tkmapden; + long layerfound[23]; + long layertotal[23]; + map > layerfound_perBx; + map > layertotal_perBx; + vector layerfound_vsLumi; + vector layertotal_vsLumi; + vector layerfound_vsPU; + vector layertotal_vsPU; + vector layerfound_vsCM; + vector layertotal_vsCM; + int goodlayertotal[35]; + int goodlayerfound[35]; + int alllayertotal[35]; + int alllayerfound[35]; + map BadModules; }; -SiStripHitResolFromCalibTree::SiStripHitResolFromCalibTree(const edm::ParameterSet& conf) : - ConditionDBWriter(conf), - FileInPath_("CalibTracker/SiStripCommon/data/SiStripDetInfo.dat") -{ +SiStripHitResolFromCalibTree::SiStripHitResolFromCalibTree(const edm::ParameterSet& conf) + : ConditionDBWriter(conf), FileInPath_("CalibTracker/SiStripCommon/data/SiStripDetInfo.dat") { CalibTreeFilenames = conf.getUntrackedParameter >("CalibTreeFilenames"); threshold = conf.getParameter("Threshold"); nModsMin = conf.getParameter("nModsMin"); doSummary = conf.getParameter("doSummary"); - _badModulesFile = conf.getUntrackedParameter("BadModulesFile", ""); + _badModulesFile = conf.getUntrackedParameter("BadModulesFile", ""); _autoIneffModTagging = conf.getUntrackedParameter("AutoIneffModTagging", false); - _clusterMatchingMethod = conf.getUntrackedParameter("ClusterMatchingMethod",0); - _ResXSig = conf.getUntrackedParameter("ResXSig",-1); - _clusterTrajDist = conf.getUntrackedParameter("ClusterTrajDist",64.0); - _stripsApvEdge = conf.getUntrackedParameter("StripsApvEdge",10.0); + _clusterMatchingMethod = conf.getUntrackedParameter("ClusterMatchingMethod", 0); + _ResXSig = conf.getUntrackedParameter("ResXSig", -1); + _clusterTrajDist = conf.getUntrackedParameter("ClusterTrajDist", 64.0); + _stripsApvEdge = conf.getUntrackedParameter("StripsApvEdge", 10.0); _useOnlyHighPurityTracks = conf.getUntrackedParameter("UseOnlyHighPurityTracks", true); - _bunchx = conf.getUntrackedParameter("BunchCrossing",0); - _spaceBetweenTrains = conf.getUntrackedParameter("SpaceBetweenTrains",25); - _useCM = conf.getUntrackedParameter("UseCommonMode",false); - _showEndcapSides = conf.getUntrackedParameter("ShowEndcapSides",true); - _showRings = conf.getUntrackedParameter("ShowRings",false); - _showTOB6TEC9 = conf.getUntrackedParameter("ShowTOB6TEC9",false); - _showOnlyGoodModules = conf.getUntrackedParameter("ShowOnlyGoodModules",false); - _tkMapMin = conf.getUntrackedParameter("TkMapMin",0.9); - _effPlotMin = conf.getUntrackedParameter("EffPlotMin",0.9); - _title = conf.getParameter("Title"); - reader = new SiStripDetInfoFileReader(FileInPath_.fullPath()); - - nTEClayers = 9; // number of wheels - if(_showRings) nTEClayers = 7; // number of rings - - quality_ = new SiStripQuality; + _bunchx = conf.getUntrackedParameter("BunchCrossing", 0); + _spaceBetweenTrains = conf.getUntrackedParameter("SpaceBetweenTrains", 25); + _useCM = conf.getUntrackedParameter("UseCommonMode", false); + _showEndcapSides = conf.getUntrackedParameter("ShowEndcapSides", true); + _showRings = conf.getUntrackedParameter("ShowRings", false); + _showTOB6TEC9 = conf.getUntrackedParameter("ShowTOB6TEC9", false); + _showOnlyGoodModules = conf.getUntrackedParameter("ShowOnlyGoodModules", false); + _tkMapMin = conf.getUntrackedParameter("TkMapMin", 0.9); + _effPlotMin = conf.getUntrackedParameter("EffPlotMin", 0.9); + _title = conf.getParameter("Title"); + _tkGeomToken = esConsumes(); + _tTopoToken = esConsumes(); + _detInfo = SiStripDetInfoFileReader::read(FileInPath_.fullPath()); + + nTEClayers = 9; // number of wheels + if (_showRings) + nTEClayers = 7; // number of rings + + quality_ = new SiStripQuality(_detInfo); } -SiStripHitResolFromCalibTree::~SiStripHitResolFromCalibTree() { } +SiStripHitResolFromCalibTree::~SiStripHitResolFromCalibTree() {} void SiStripHitResolFromCalibTree::algoBeginJob(const edm::EventSetup&) { //I have no idea what goes here @@ -216,57 +219,49 @@ void SiStripHitResolFromCalibTree::algoBeginJob(const edm::EventSetup&) { void SiStripHitResolFromCalibTree::algoEndJob() { //Still have no idea what goes here - } void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::EventSetup& c) { - - edm::ESHandle tracker; - c.get().get( tracker ); - const TrackerGeometry * tkgeom=&(* tracker); - - //Retrieve tracker topology from geometry - edm::ESHandle tTopoHandle; - c.get().get(tTopoHandle); - const TrackerTopology* const tTopo = tTopoHandle.product(); + const auto& tkgeom = c.getData(_tkGeomToken); + const auto& tTopo = c.getData(_tTopoToken); // read bad modules to mask ifstream badModules_file; set badModules_list; - if(!_badModulesFile.empty()) { - badModules_file.open(_badModulesFile.c_str()); - uint32_t badmodule_detid; - int mods, fiber1, fiber2, fiber3; - if(badModules_file.is_open()) { + if (!_badModulesFile.empty()) { + badModules_file.open(_badModulesFile.c_str()); + uint32_t badmodule_detid; + int mods, fiber1, fiber2, fiber3; + if (badModules_file.is_open()) { string line; - while ( getline (badModules_file,line) ) { - if(badModules_file.eof()) continue; - stringstream ss(line); - ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3; - if(badmodule_detid!=0 && mods==1 && (fiber1==1 || fiber2==1 || fiber3==1) ) - badModules_list.insert(badmodule_detid); - } + while (getline(badModules_file, line)) { + if (badModules_file.eof()) + continue; + stringstream ss(line); + ss >> badmodule_detid >> mods >> fiber1 >> fiber2 >> fiber3; + if (badmodule_detid != 0 && mods == 1 && (fiber1 == 1 || fiber2 == 1 || fiber3 == 1)) + badModules_list.insert(badmodule_detid); + } badModules_file.close(); - } + } } - if(!badModules_list.empty()) cout<<"Remove additionnal bad modules from the analysis: "<::iterator itBadMod; - for (itBadMod=badModules_list.begin(); itBadMod!=badModules_list.end(); ++itBadMod) - cout<<" "<<*itBadMod<make("bx","bx",3600,0,3600); - instLumiHisto = fs->make("instLumi","inst. lumi.",250,0,25000); - PUHisto = fs->make("PU","PU",200,0,200); - - - for(int l=0; l < 35; l++) { + + bxHisto = fs->make("bx", "bx", 3600, 0, 3600); + instLumiHisto = fs->make("instLumi", "inst. lumi.", 250, 0, 25000); + PUHisto = fs->make("PU", "PU", 200, 0, 200); + + for (int l = 0; l < 35; l++) { goodlayertotal[l] = 0; goodlayerfound[l] = 0; alllayertotal[l] = 0; - alllayerfound[l] = 0; + alllayerfound[l] = 0; } TH1F* PredPlots_m[23]; @@ -282,61 +277,78 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E std::string PlotTitleTrajX; std::string FileNameEnding; - if(UnitString == "cm"){PlotTitleClusX = "clusX [cm]"; PlotTitleTrajX = "trajX [cm]"; FileNameEnding = "CM";} - else if(UnitString == "strip unit"){PlotTitleClusX = "clusX [strip unit]"; PlotTitleTrajX = "trajX [strip unit]"; FileNameEnding = "StripUnit";} - else{std::cout << "ERROR: Unit must either be cm or strip unit." << std::endl; } - - for(Long_t ilayer = 0; ilayer <23; ilayer++) { - - MeasPlots_m[ilayer] = fs->make(Form("MeasPlots_m_layer_%i",(int)(ilayer)),GetLayerName(ilayer),1000,-5., 5.); - MeasPlots_m[ilayer]->GetXaxis()->SetTitle("clusX [cm]"); - MeasPlots_p[ilayer] = fs->make(Form("MeasPlots_p_layer_%i",(int)(ilayer)),GetLayerName(ilayer),1000,-400,400); - MeasPlots_p[ilayer]->GetXaxis()->SetTitle("clusX [strip]"); - - PredPlots_m[ilayer] = fs->make(Form("PredPlots_m_layer_%i",(int)(ilayer)),GetLayerName(ilayer),1000,-5., 5.); - PredPlots_m[ilayer]->GetXaxis()->SetTitle("trajX [cm]"); - PredPlots_p[ilayer] = fs->make(Form("PredPlots_p_layer_%i",(int)(ilayer)),GetLayerName(ilayer),1000,-400,400); - PredPlots_p[ilayer]->GetXaxis()->SetTitle("trajX [strip]"); - - ResidPlots_m[ilayer] = fs->make(Form("ResidPlots_m_layer_%i",(int)(ilayer)),GetLayerName(ilayer),250,-.125, .125); - ResidPlots_m[ilayer]->GetXaxis()->SetTitle("trajX [cm]"); - ResidPlots_p[ilayer] = fs->make(Form("ResidPlots_p_layer_%i",(int)(ilayer)),GetLayerName(ilayer),200,-10.,10); - ResidPlots_p[ilayer]->GetXaxis()->SetTitle("trajX [strip]"); - - - layerfound_vsLumi.push_back( fs->make(Form("layerfound_vsLumi_layer_%i",(int)(ilayer)),GetLayerName(ilayer),100,0,25000)); - layertotal_vsLumi.push_back( fs->make(Form("layertotal_vsLumi_layer_%i",(int)(ilayer)),GetLayerName(ilayer),100,0,25000)); - layerfound_vsPU.push_back( fs->make(Form("layerfound_vsPU_layer_%i",(int)(ilayer)),GetLayerName(ilayer),45,0,90)); - layertotal_vsPU.push_back( fs->make(Form("layertotal_vsPU_layer_%i",(int)(ilayer)),GetLayerName(ilayer),45,0,90)); - - if(_useCM) { - layerfound_vsCM.push_back( fs->make(Form("layerfound_vsCM_layer_%i",(int)(ilayer)),GetLayerName(ilayer),20,0,400)); - layertotal_vsCM.push_back( fs->make(Form("layertotal_vsCM_layer_%i",(int)(ilayer)),GetLayerName(ilayer),20,0,400)); - } - layertotal[ilayer] = 0; - layerfound[ilayer] = 0; + if (UnitString == "cm") { + PlotTitleClusX = "clusX [cm]"; + PlotTitleTrajX = "trajX [cm]"; + FileNameEnding = "CM"; + } else if (UnitString == "strip unit") { + PlotTitleClusX = "clusX [strip unit]"; + PlotTitleTrajX = "trajX [strip unit]"; + FileNameEnding = "StripUnit"; + } else { + std::cout << "ERROR: Unit must either be cm or strip unit." << std::endl; + } + for (Long_t ilayer = 0; ilayer < 23; ilayer++) { + MeasPlots_m[ilayer] = + fs->make(Form("MeasPlots_m_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 1000, -5., 5.); + MeasPlots_m[ilayer]->GetXaxis()->SetTitle("clusX [cm]"); + MeasPlots_p[ilayer] = + fs->make(Form("MeasPlots_p_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 1000, -400, 400); + MeasPlots_p[ilayer]->GetXaxis()->SetTitle("clusX [strip]"); + + PredPlots_m[ilayer] = + fs->make(Form("PredPlots_m_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 1000, -5., 5.); + PredPlots_m[ilayer]->GetXaxis()->SetTitle("trajX [cm]"); + PredPlots_p[ilayer] = + fs->make(Form("PredPlots_p_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 1000, -400, 400); + PredPlots_p[ilayer]->GetXaxis()->SetTitle("trajX [strip]"); + + ResidPlots_m[ilayer] = + fs->make(Form("ResidPlots_m_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 250, -.125, .125); + ResidPlots_m[ilayer]->GetXaxis()->SetTitle("trajX [cm]"); + ResidPlots_p[ilayer] = + fs->make(Form("ResidPlots_p_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 200, -10., 10); + ResidPlots_p[ilayer]->GetXaxis()->SetTitle("trajX [strip]"); + + layerfound_vsLumi.push_back( + fs->make(Form("layerfound_vsLumi_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 100, 0, 25000)); + layertotal_vsLumi.push_back( + fs->make(Form("layertotal_vsLumi_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 100, 0, 25000)); + layerfound_vsPU.push_back( + fs->make(Form("layerfound_vsPU_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 45, 0, 90)); + layertotal_vsPU.push_back( + fs->make(Form("layertotal_vsPU_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 45, 0, 90)); + + if (_useCM) { + layerfound_vsCM.push_back( + fs->make(Form("layerfound_vsCM_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 20, 0, 400)); + layertotal_vsCM.push_back( + fs->make(Form("layertotal_vsCM_layer_%i", (int)(ilayer)), GetLayerName(ilayer), 20, 0, 400)); + } + layertotal[ilayer] = 0; + layerfound[ilayer] = 0; } - if(!_autoIneffModTagging) cout << "A module is bad if efficiency < " << threshold << " and has at least " << nModsMin << " nModsMin." << endl; - else cout << "A module is bad if efficiency < the avg in layer - " << threshold << " and has at least " << nModsMin << " nModsMin." << endl; - - + if (!_autoIneffModTagging) + cout << "A module is bad if efficiency < " << threshold << " and has at least " << nModsMin << " nModsMin." << endl; + else + cout << "A module is bad if efficiency < the avg in layer - " << threshold << " and has at least " << nModsMin + << " nModsMin." << endl; + unsigned int run, evt, bx; double instLumi, PU; //Open the ROOT Calib Tree - for( unsigned int ifile=0; ifile < CalibTreeFilenames.size(); ifile++) { + for (unsigned int ifile = 0; ifile < CalibTreeFilenames.size(); ifile++) { + cout << "Loading file: " << CalibTreeFilenames[ifile] << endl; + TFile* CalibTreeFile = TFile::Open(CalibTreeFilenames[ifile].c_str(), "READ"); - cout<<"Loading file: "<cd("eventInfo"); - } - catch(exception& e) { + } catch (exception& e) { cout << "No event infos tree" << endl; } TTree* EventTree = (TTree*)(gDirectory->Get("tree")); @@ -346,8 +358,7 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E TLeaf* BunchLf; TLeaf* InstLumiLf; TLeaf* PULf; - if(EventTree){ - + if (EventTree) { cout << "Found event infos tree" << endl; runLf = EventTree->GetLeaf("run"); @@ -358,9 +369,10 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E PULf = EventTree->GetLeaf("PU"); int nevt = EventTree->GetEntries(); - if(nevt) foundEventInfos=true; + if (nevt) + foundEventInfos = true; - for(int j =0; j < nevt; j++) { + for (int j = 0; j < nevt; j++) { EventTree->GetEntry(j); run = runLf->GetValue(); evt = evtLf->GetValue(); @@ -372,62 +384,61 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E instLumiHisto->Fill(instLumi); PUHisto->Fill(PU); - eventInfos[ make_pair(run, evt) ] = array {{ (double) bx, instLumi, PU}}; + eventInfos[make_pair(run, evt)] = array{{(double)bx, instLumi, PU}}; } - } - - - // Get hit infos - CalibTreeFile->cd("anEff"); - CalibTree = (TTree*)(gDirectory->Get("traj")); - - runLf = CalibTree->GetLeaf("run"); - evtLf = CalibTree->GetLeaf("event"); - TLeaf* BadLf = CalibTree->GetLeaf("ModIsBad"); - TLeaf* sistripLf = CalibTree->GetLeaf("SiStripQualBad"); - TLeaf* idLf = CalibTree->GetLeaf("Id"); - TLeaf* acceptLf = CalibTree->GetLeaf("withinAcceptance"); - TLeaf* layerLf = CalibTree->GetLeaf("layer"); - //TLeaf* nHitsLf = CalibTree->GetLeaf("nHits"); - TLeaf* highPurityLf = CalibTree->GetLeaf("highPurity"); - TLeaf* xLf = CalibTree->GetLeaf("TrajGlbX"); - TLeaf* yLf = CalibTree->GetLeaf("TrajGlbY"); - TLeaf* zLf = CalibTree->GetLeaf("TrajGlbZ"); - TLeaf* ResXSigLf = CalibTree->GetLeaf("ResXSig"); - TLeaf* TrajLocXLf = CalibTree->GetLeaf("TrajLocX"); - TLeaf* TrajLocYLf = CalibTree->GetLeaf("TrajLocY"); - TLeaf* ClusterLocXLf = CalibTree->GetLeaf("ClusterLocX"); - BunchLf = CalibTree->GetLeaf("bunchx"); - InstLumiLf = CalibTree->GetLeaf("instLumi"); - PULf = CalibTree->GetLeaf("PU"); - TLeaf* CMLf = nullptr; - if(_useCM) CMLf = CalibTree->GetLeaf("commonMode"); - - int nevents = CalibTree->GetEntries(); - cout << "Successfully loaded analyze function with " << nevents << " events!\n"; - - - map< pair< unsigned int, unsigned int>, array >::iterator itEventInfos; - - - //Loop through all of the events - for(int j =0; j < nevents; j++) { + + // Get hit infos + CalibTreeFile->cd("anEff"); + CalibTree = (TTree*)(gDirectory->Get("traj")); + + runLf = CalibTree->GetLeaf("run"); + evtLf = CalibTree->GetLeaf("event"); + TLeaf* BadLf = CalibTree->GetLeaf("ModIsBad"); + TLeaf* sistripLf = CalibTree->GetLeaf("SiStripQualBad"); + TLeaf* idLf = CalibTree->GetLeaf("Id"); + TLeaf* acceptLf = CalibTree->GetLeaf("withinAcceptance"); + TLeaf* layerLf = CalibTree->GetLeaf("layer"); + //TLeaf* nHitsLf = CalibTree->GetLeaf("nHits"); + TLeaf* highPurityLf = CalibTree->GetLeaf("highPurity"); + TLeaf* xLf = CalibTree->GetLeaf("TrajGlbX"); + TLeaf* yLf = CalibTree->GetLeaf("TrajGlbY"); + TLeaf* zLf = CalibTree->GetLeaf("TrajGlbZ"); + TLeaf* ResXSigLf = CalibTree->GetLeaf("ResXSig"); + TLeaf* TrajLocXLf = CalibTree->GetLeaf("TrajLocX"); + TLeaf* TrajLocYLf = CalibTree->GetLeaf("TrajLocY"); + TLeaf* ClusterLocXLf = CalibTree->GetLeaf("ClusterLocX"); + BunchLf = CalibTree->GetLeaf("bunchx"); + InstLumiLf = CalibTree->GetLeaf("instLumi"); + PULf = CalibTree->GetLeaf("PU"); + TLeaf* CMLf = nullptr; + if (_useCM) + CMLf = CalibTree->GetLeaf("commonMode"); + + int nevents = CalibTree->GetEntries(); + cout << "Successfully loaded analyze function with " << nevents << " events!\n"; + + map, array >::iterator itEventInfos; + + //Loop through all of the events + for (int j = 0; j < nevents; j++) { CalibTree->GetEntry(j); - run = (unsigned int)runLf->GetValue(); - evt = (unsigned int)evtLf->GetValue(); + run = (unsigned int)runLf->GetValue(); + evt = (unsigned int)evtLf->GetValue(); unsigned int isBad = (unsigned int)BadLf->GetValue(); unsigned int quality = (unsigned int)sistripLf->GetValue(); unsigned int id = (unsigned int)idLf->GetValue(); unsigned int accept = (unsigned int)acceptLf->GetValue(); unsigned int layer_wheel = (unsigned int)layerLf->GetValue(); unsigned int layer = layer_wheel; - if(_showRings && layer >10) { // use rings instead of wheels - if(layer<14) layer = 10 + ((id>>9)&0x3); //TID 3 disks and also 3 rings -> use the same container - else layer = 13 + ((id>>5)&0x7); //TEC + if (_showRings && layer > 10) { // use rings instead of wheels + if (layer < 14) + layer = 10 + ((id >> 9) & 0x3); //TID 3 disks and also 3 rings -> use the same container + else + layer = 13 + ((id >> 5) & 0x7); //TEC } //unsigned int nHits = (unsigned int)nHitsLf->GetValue(); - bool highPurity = (bool) highPurityLf->GetValue(); + bool highPurity = (bool)highPurityLf->GetValue(); double x = xLf->GetValue(); double y = yLf->GetValue(); double z = zLf->GetValue(); @@ -439,70 +450,75 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E double stripTrajMid; double stripCluster; bool badquality = false; - - instLumi = 0; - PU = 0; - - // if no special tree with event infos, they may be stored in the hit eff tree - if(!foundEventInfos){ - bx = (unsigned int)BunchLf->GetValue(); - if(InstLumiLf!=nullptr) instLumi = InstLumiLf->GetValue(); // branch not filled by default - if(PULf!=nullptr) PU = PULf->GetValue(); // branch not filled by default - } - int CM = -100; - if(_useCM) CM = CMLf->GetValue(); - - - // Get infos from eventInfos if they exist - if(foundEventInfos){ - - itEventInfos = eventInfos.find( make_pair(run, evt) ); - if(itEventInfos!=eventInfos.end()){ - bx = itEventInfos->second[0]; - instLumi = itEventInfos->second[1]; - PU = itEventInfos->second[2]; - } - } + instLumi = 0; + PU = 0; + // if no special tree with event infos, they may be stored in the hit eff tree + if (!foundEventInfos) { + bx = (unsigned int)BunchLf->GetValue(); + if (InstLumiLf != nullptr) + instLumi = InstLumiLf->GetValue(); // branch not filled by default + if (PULf != nullptr) + PU = PULf->GetValue(); // branch not filled by default + } + int CM = -100; + if (_useCM) + CM = CMLf->GetValue(); + + // Get infos from eventInfos if they exist + if (foundEventInfos) { + itEventInfos = eventInfos.find(make_pair(run, evt)); + if (itEventInfos != eventInfos.end()) { + bx = itEventInfos->second[0]; + instLumi = itEventInfos->second[1]; + PU = itEventInfos->second[2]; + } + } //We have two things we want to do, both an XY color plot, and the efficiency measurement //First, ignore anything that isn't in acceptance and isn't good quality - - if(_bunchx > 0 && _bunchx != bx) continue; + + if (_bunchx > 0 && _bunchx != bx) + continue; //if(quality == 1 || accept != 1 || nHits < 8) continue; - if(accept != 1) continue; - if(_useOnlyHighPurityTracks && !highPurity) continue; - if(quality == 1) badquality = true; + if (accept != 1) + continue; + if (_useOnlyHighPurityTracks && !highPurity) + continue; + if (quality == 1) + badquality = true; // don't compute efficiencies in modules from TOB6 and TEC9 - if(!_showTOB6TEC9 && (layer_wheel==10 || layer_wheel==22)) continue; - - // don't use bad modules given in the bad module list - itBadMod = badModules_list.find(id); - if(itBadMod!=badModules_list.end()) continue; + if (!_showTOB6TEC9 && (layer_wheel == 10 || layer_wheel == 22)) + continue; + // don't use bad modules given in the bad module list + itBadMod = badModules_list.find(id); + if (itBadMod != badModules_list.end()) + continue; //Now that we have a good event, we need to look at if we expected it or not, and the location //if we didn't //Fill the missing hit information first bool badflag = false; - // By default uses the old matching method - if(_ResXSig < 0) { - if(isBad == 1) badflag = true; // isBad set to false in the tree when resxsig<999.0 - } - else { - if(isBad == 1 || resxsig > _ResXSig) badflag = true; + // By default uses the old matching method + if (_ResXSig < 0) { + if (isBad == 1) + badflag = true; // isBad set to false in the tree when resxsig<999.0 + } else { + if (isBad == 1 || resxsig > _ResXSig) + badflag = true; } - // Conversion of positions in strip unit - int nstrips = -9; - float Pitch = -9.0; + // Conversion of positions in strip unit + int nstrips = -9; + float Pitch = -9.0; //For converting from pitch units into micrometres - //const StripGeomDetUnit * conversion = dynamic_cast(tkgeom->idToDetUnit(DetId)); + //const StripGeomDetUnit * conversion = dynamic_cast(tkgeom->idToDetUnit(DetId)); //auto SpecTopo = conversion->specificTopology(); //const StripGeomDetUnit * conversion=(const StripGeomDetUnit*)tkgeom->idToDetUnit(DetId); @@ -512,41 +528,40 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E //std::cout << "SpecTopo = " << SpecTopo << std::endl; //std::cout << '\n' << std::endl; - - if (resxsig==1000.0) { // special treatment, no GeomDetUnit associated in some cases when no cluster found - Pitch = 0.0205; // maximum - nstrips = 768; // maximum - stripTrajMid = TrajLocX/Pitch + nstrips/2.0 ; - stripCluster = ClusterLocX/Pitch + nstrips/2.0 ; + if (resxsig == 1000.0) { // special treatment, no GeomDetUnit associated in some cases when no cluster found + Pitch = 0.0205; // maximum + nstrips = 768; // maximum + stripTrajMid = TrajLocX / Pitch + nstrips / 2.0; + stripCluster = ClusterLocX / Pitch + nstrips / 2.0; + } else { + DetId ClusterDetId(id); + const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)tkgeom.idToDetUnit(ClusterDetId); + const StripTopology& Topo = stripdet->specificTopology(); + nstrips = Topo.nstrips(); + Pitch = stripdet->surface().bounds().width() / Topo.nstrips(); + stripTrajMid = TrajLocX / Pitch + nstrips / 2.0; //layer01->10 + stripCluster = ClusterLocX / Pitch + nstrips / 2.0; + + // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module + // for correct comparison with cluster position + float hbedge = 0; + float htedge = 0; + float hapoth = 0; + if (layer >= 11) { + const BoundPlane& plane = stripdet->surface(); + const TrapezoidalPlaneBounds* trapezoidalBounds( + dynamic_cast(&(plane.bounds()))); + std::array const& parameters = (*trapezoidalBounds).parameters(); + hbedge = parameters[0]; + htedge = parameters[1]; + hapoth = parameters[3]; + TrajLocXMid = TrajLocX / (1 + (htedge - hbedge) * TrajLocY / (htedge + hbedge) / + hapoth); // radialy extrapolated x loc position at middle + stripTrajMid = TrajLocXMid / Pitch + nstrips / 2.0; + } } - else { - DetId ClusterDetId(id); - const StripGeomDetUnit * stripdet=(const StripGeomDetUnit*)tkgeom->idToDetUnit(ClusterDetId); - const StripTopology& Topo = stripdet->specificTopology(); - nstrips = Topo.nstrips(); - Pitch = stripdet->surface().bounds().width() / Topo.nstrips(); - stripTrajMid = TrajLocX/Pitch + nstrips/2.0 ; //layer01->10 - stripCluster = ClusterLocX/Pitch + nstrips/2.0 ; - - // For trapezoidal modules: extrapolation of x trajectory position to the y middle of the module - // for correct comparison with cluster position - float hbedge = 0; - float htedge = 0; - float hapoth = 0; - if(layer>=11) { - const BoundPlane& plane = stripdet->surface(); - const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast(&(plane.bounds()))); - std::array const & parameters = (*trapezoidalBounds).parameters(); - hbedge = parameters[0]; - htedge = parameters[1]; - hapoth = parameters[3]; - TrajLocXMid = TrajLocX / (1 + (htedge-hbedge)*TrajLocY/(htedge+hbedge)/hapoth) ; // radialy extrapolated x loc position at middle - stripTrajMid = TrajLocXMid/Pitch + nstrips/2.0 ; - } - } - - - /* + + /* double MeasPlotsVariable; double PredPlotsVariable; @@ -555,142 +570,153 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E else{std::cout << "ERROR: Unit must be cm or strip unit" << std::endl;} */ - if(!badquality && layer<23) { - if(resxsig!=1000.0){ - MeasPlots_m[layer]->Fill(ClusterLocX); - MeasPlots_p[layer]->Fill(ClusterLocX/Pitch); - PredPlots_m[layer]->Fill(TrajLocX); - PredPlots_p[layer]->Fill(TrajLocX/Pitch); - ResidPlots_m[layer]->Fill(ClusterLocX-TrajLocX); - ResidPlots_p[layer]->Fill((ClusterLocX-TrajLocX)/Pitch); - } else { - MeasPlots_m[layer]->Fill(1000); - MeasPlots_p[layer]->Fill(1000); - PredPlots_m[layer]->Fill(1000); - PredPlots_p[layer]->Fill(1000); - } - } - - // New matching methods - int tapv = -9; - int capv = -9; - float stripInAPV = 64.; - - if ( _clusterMatchingMethod >=1 ) { - badflag = false; // reset - if(resxsig == 1000.0) { // default value when no cluster found in the module - badflag = true; // consider the module inefficient in this case - } - else{ - if (_clusterMatchingMethod==2 || _clusterMatchingMethod==4) { // check the distance between cluster and trajectory position - if ( abs(stripCluster - stripTrajMid) > _clusterTrajDist ) badflag = true; - } - if (_clusterMatchingMethod==3 || _clusterMatchingMethod==4) { // cluster and traj have to be in the same APV (don't take edges into accounts) - tapv = (int) stripTrajMid/128; - capv = (int) stripCluster/128; - stripInAPV = stripTrajMid-tapv*128; - - if(stripInAPV<_stripsApvEdge || stripInAPV>128-_stripsApvEdge) continue; - if(tapv != capv) badflag = true; - } - } - } - - - - if(badflag && !badquality) { - hit temphit; - temphit.x = x; - temphit.y = y; - temphit.z = z; - temphit.id = id; - hits[layer].push_back(temphit); - } - pair newgoodpair (1,1); - pair newbadpair (1,0); + if (!badquality && layer < 23) { + if (resxsig != 1000.0) { + MeasPlots_m[layer]->Fill(ClusterLocX); + MeasPlots_p[layer]->Fill(ClusterLocX / Pitch); + PredPlots_m[layer]->Fill(TrajLocX); + PredPlots_p[layer]->Fill(TrajLocX / Pitch); + ResidPlots_m[layer]->Fill(ClusterLocX - TrajLocX); + ResidPlots_p[layer]->Fill((ClusterLocX - TrajLocX) / Pitch); + } else { + MeasPlots_m[layer]->Fill(1000); + MeasPlots_p[layer]->Fill(1000); + PredPlots_m[layer]->Fill(1000); + PredPlots_p[layer]->Fill(1000); + } + } + + // New matching methods + int tapv = -9; + int capv = -9; + float stripInAPV = 64.; + + if (_clusterMatchingMethod >= 1) { + badflag = false; // reset + if (resxsig == 1000.0) { // default value when no cluster found in the module + badflag = true; // consider the module inefficient in this case + } else { + if (_clusterMatchingMethod == 2 || + _clusterMatchingMethod == 4) { // check the distance between cluster and trajectory position + if (abs(stripCluster - stripTrajMid) > _clusterTrajDist) + badflag = true; + } + if (_clusterMatchingMethod == 3 || + _clusterMatchingMethod == + 4) { // cluster and traj have to be in the same APV (don't take edges into accounts) + tapv = (int)stripTrajMid / 128; + capv = (int)stripCluster / 128; + stripInAPV = stripTrajMid - tapv * 128; + + if (stripInAPV < _stripsApvEdge || stripInAPV > 128 - _stripsApvEdge) + continue; + if (tapv != capv) + badflag = true; + } + } + } + + if (badflag && !badquality) { + hit temphit; + temphit.x = x; + temphit.y = y; + temphit.z = z; + temphit.id = id; + hits[layer].push_back(temphit); + } + pair newgoodpair(1, 1); + pair newbadpair(1, 0); //First, figure out if the module already exists in the map of maps - map< unsigned int, pair< unsigned int, unsigned int> >::iterator it = modCounter[layer].find(id); - if(!badquality) { - if(it == modCounter[layer].end()) { - if(badflag) modCounter[layer][id] = newbadpair; - else modCounter[layer][id] = newgoodpair; - } - else { + map >::iterator it = modCounter[layer].find(id); + if (!badquality) { + if (it == modCounter[layer].end()) { + if (badflag) + modCounter[layer][id] = newbadpair; + else + modCounter[layer][id] = newgoodpair; + } else { ((*it).second.first)++; - if(!badflag) ((*it).second.second)++; - } - - if(layerfound_perBx.find(bx)==layerfound_perBx.end()) { - layerfound_perBx[bx] = vector(23, 0); - layertotal_perBx[bx] = vector(23, 0); - } - if(!badflag) layerfound_perBx[bx][layer]++; - layertotal_perBx[bx][layer]++; - - if(!badflag) layerfound_vsLumi[layer]->Fill(instLumi); - layertotal_vsLumi[layer]->Fill(instLumi); - if(!badflag) layerfound_vsPU[layer]->Fill(PU); - layertotal_vsPU[layer]->Fill(PU); - - if(_useCM){ - if(!badflag) layerfound_vsCM[layer]->Fill(CM); - layertotal_vsCM[layer]->Fill(CM); - } - - //Have to do the decoding for which side to go on (ugh) - if(layer <= 10) { - if(!badflag) goodlayerfound[layer]++; + if (!badflag) + ((*it).second.second)++; + } + + if (layerfound_perBx.find(bx) == layerfound_perBx.end()) { + layerfound_perBx[bx] = vector(23, 0); + layertotal_perBx[bx] = vector(23, 0); + } + if (!badflag) + layerfound_perBx[bx][layer]++; + layertotal_perBx[bx][layer]++; + + if (!badflag) + layerfound_vsLumi[layer]->Fill(instLumi); + layertotal_vsLumi[layer]->Fill(instLumi); + if (!badflag) + layerfound_vsPU[layer]->Fill(PU); + layertotal_vsPU[layer]->Fill(PU); + + if (_useCM) { + if (!badflag) + layerfound_vsCM[layer]->Fill(CM); + layertotal_vsCM[layer]->Fill(CM); + } + + //Have to do the decoding for which side to go on (ugh) + if (layer <= 10) { + if (!badflag) + goodlayerfound[layer]++; goodlayertotal[layer]++; - } - else if(layer > 10 && layer < 14) { - if( ((id>>13)&0x3) == 1) { - if(!badflag) goodlayerfound[layer]++; - goodlayertotal[layer]++; - } - else if( ((id>>13)&0x3) == 2) { - if(!badflag) goodlayerfound[layer+3]++; - goodlayertotal[layer+3]++; - } - } - else if(layer > 13 && layer <= 22) { - if( ((id>>18)&0x3) == 1) { - if(!badflag) goodlayerfound[layer+3]++; - goodlayertotal[layer+3]++; - } - else if( ((id>>18)&0x3) == 2) { - if(!badflag) goodlayerfound[layer+3+nTEClayers]++; - goodlayertotal[layer+3+nTEClayers]++; - } - } + } else if (layer > 10 && layer < 14) { + if (((id >> 13) & 0x3) == 1) { + if (!badflag) + goodlayerfound[layer]++; + goodlayertotal[layer]++; + } else if (((id >> 13) & 0x3) == 2) { + if (!badflag) + goodlayerfound[layer + 3]++; + goodlayertotal[layer + 3]++; + } + } else if (layer > 13 && layer <= 22) { + if (((id >> 18) & 0x3) == 1) { + if (!badflag) + goodlayerfound[layer + 3]++; + goodlayertotal[layer + 3]++; + } else if (((id >> 18) & 0x3) == 2) { + if (!badflag) + goodlayerfound[layer + 3 + nTEClayers]++; + goodlayertotal[layer + 3 + nTEClayers]++; + } + } } //Do the one where we don't exclude bad modules! - if(layer <= 10) { - if(!badflag) alllayerfound[layer]++; - alllayertotal[layer]++; - } - else if(layer > 10 && layer < 14) { - if( ((id>>13)&0x3) == 1) { - if(!badflag) alllayerfound[layer]++; + if (layer <= 10) { + if (!badflag) + alllayerfound[layer]++; + alllayertotal[layer]++; + } else if (layer > 10 && layer < 14) { + if (((id >> 13) & 0x3) == 1) { + if (!badflag) + alllayerfound[layer]++; alllayertotal[layer]++; - } - else if( ((id>>13)&0x3) == 2) { - if(!badflag) alllayerfound[layer+3]++; - alllayertotal[layer+3]++; - } + } else if (((id >> 13) & 0x3) == 2) { + if (!badflag) + alllayerfound[layer + 3]++; + alllayertotal[layer + 3]++; + } + } else if (layer > 13 && layer <= 22) { + if (((id >> 18) & 0x3) == 1) { + if (!badflag) + alllayerfound[layer + 3]++; + alllayertotal[layer + 3]++; + } else if (((id >> 18) & 0x3) == 2) { + if (!badflag) + alllayerfound[layer + 3 + nTEClayers]++; + alllayertotal[layer + 3 + nTEClayers]++; + } } - else if(layer > 13 && layer <= 22) { - if( ((id>>18)&0x3) == 1) { - if(!badflag) alllayerfound[layer+3]++; - alllayertotal[layer+3]++; - } - else if( ((id>>18)&0x3) == 2) { - if(!badflag) alllayerfound[layer+3+nTEClayers]++; - alllayertotal[layer+3+nTEClayers]++; - } - } //At this point, both of our maps are loaded with the correct information } - }// go to next CalibTreeFile + } // go to next CalibTreeFile makeHotColdMaps(); makeTKMap(_autoIneffModTagging); @@ -699,350 +725,365 @@ void SiStripHitResolFromCalibTree::algoAnalyze(const edm::Event& e, const edm::E makeSummary(); makeSummaryVsBx(); makeSummaryVsLumi(); - if(_useCM) makeSummaryVsCM(); - + if (_useCM) + makeSummaryVsCM(); + //////////////////////////////////////////////////////////////////////// //try to write out what's in the quality record ///////////////////////////////////////////////////////////////////////////// - int NTkBadComponent[4]; //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips - int NBadComponent[4][19][4]; + int NTkBadComponent[4]; //k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips + int NBadComponent[4][19][4]; //legend: NBadComponent[i][j][k]= SubSystem i, layer/disk/wheel j, BadModule/Fiber/Apv k // i: 0=TIB, 1=TID, 2=TOB, 3=TEC // k: 0=BadModule, 1=BadFiber, 2=BadApv, 3=BadStrips std::stringstream ssV[4][19]; - - for(int i=0;i<4;++i){ - NTkBadComponent[i]=0; - for(int j=0;j<19;++j){ + + for (int i = 0; i < 4; ++i) { + NTkBadComponent[i] = 0; + for (int j = 0; j < 19; ++j) { ssV[i][j].str(""); - for(int k=0;k<4;++k) - NBadComponent[i][j][k]=0; + for (int k = 0; k < 4; ++k) + NBadComponent[i][j][k] = 0; } } - - + std::vector BC = quality_->getBadComponentList(); - - for (size_t i=0;i>2)&0x1 )+ ( (BC[i].BadFibers>>1)&0x1 ) + ( (BC[i].BadFibers)&0x1 ); + if (BC[i].BadFibers) + NTkBadComponent[1] += ((BC[i].BadFibers >> 2) & 0x1) + ((BC[i].BadFibers >> 1) & 0x1) + ((BC[i].BadFibers) & 0x1); if (BC[i].BadApvs) - NTkBadComponent[2]+= ( (BC[i].BadApvs>>5)&0x1 )+ ( (BC[i].BadApvs>>4)&0x1 ) + ( (BC[i].BadApvs>>3)&0x1 ) + - ( (BC[i].BadApvs>>2)&0x1 )+ ( (BC[i].BadApvs>>1)&0x1 ) + ( (BC[i].BadApvs)&0x1 ); - + NTkBadComponent[2] += ((BC[i].BadApvs >> 5) & 0x1) + ((BC[i].BadApvs >> 4) & 0x1) + ((BC[i].BadApvs >> 3) & 0x1) + + ((BC[i].BadApvs >> 2) & 0x1) + ((BC[i].BadApvs >> 1) & 0x1) + ((BC[i].BadApvs) & 0x1); + //&&&&&&&&&&&&&&&&& //Single SubSystem //&&&&&&&&&&&&&&&&& - + int component; DetId a(BC[i].detid); - if ( a.subdetId() == StripSubdetector::TIB ){ + if (a.subdetId() == StripSubdetector::TIB) { //&&&&&&&&&&&&&&&&& //TIB //&&&&&&&&&&&&&&&&& - - component=tTopo->tibLayer(BC[i].detid); - SetBadComponents(0, component, BC[i], ssV, NBadComponent); - - } else if ( a.subdetId() == StripSubdetector::TID ) { + + component = tTopo.tibLayer(BC[i].detid); + SetBadComponents(0, component, BC[i], ssV, NBadComponent); + + } else if (a.subdetId() == StripSubdetector::TID) { //&&&&&&&&&&&&&&&&& //TID //&&&&&&&&&&&&&&&&& - - component=tTopo->tidSide(BC[i].detid)==2?tTopo->tidWheel(BC[i].detid):tTopo->tidWheel(BC[i].detid)+3; - SetBadComponents(1, component, BC[i], ssV, NBadComponent); - - } else if ( a.subdetId() == StripSubdetector::TOB ) { + + component = tTopo.tidSide(BC[i].detid) == 2 ? tTopo.tidWheel(BC[i].detid) : tTopo.tidWheel(BC[i].detid) + 3; + SetBadComponents(1, component, BC[i], ssV, NBadComponent); + + } else if (a.subdetId() == StripSubdetector::TOB) { //&&&&&&&&&&&&&&&&& //TOB //&&&&&&&&&&&&&&&&& - - component=tTopo->tobLayer(BC[i].detid); - SetBadComponents(2, component, BC[i], ssV, NBadComponent); - - } else if ( a.subdetId() == StripSubdetector::TEC ) { + + component = tTopo.tobLayer(BC[i].detid); + SetBadComponents(2, component, BC[i], ssV, NBadComponent); + + } else if (a.subdetId() == StripSubdetector::TEC) { //&&&&&&&&&&&&&&&&& //TEC //&&&&&&&&&&&&&&&&& - - component=tTopo->tecSide(BC[i].detid)==2?tTopo->tecWheel(BC[i].detid):tTopo->tecWheel(BC[i].detid)+9; - SetBadComponents(3, component, BC[i], ssV, NBadComponent); - - } + + component = tTopo.tecSide(BC[i].detid) == 2 ? tTopo.tecWheel(BC[i].detid) : tTopo.tecWheel(BC[i].detid) + 9; + SetBadComponents(3, component, BC[i], ssV, NBadComponent); + } } - + //&&&&&&&&&&&&&&&&&& // Single Strip Info //&&&&&&&&&&&&&&&&&& - float percentage=0; - + float percentage = 0; + SiStripQuality::RegistryIterator rbegin = quality_->getRegistryVectorBegin(); - SiStripQuality::RegistryIterator rend = quality_->getRegistryVectorEnd(); - - for (SiStripBadStrip::RegistryIterator rp=rbegin; rp != rend; ++rp) { - unsigned int detid=rp->detid; - - int subdet=-999; int component=-999; + SiStripQuality::RegistryIterator rend = quality_->getRegistryVectorEnd(); + + for (SiStripBadStrip::RegistryIterator rp = rbegin; rp != rend; ++rp) { + unsigned int detid = rp->detid; + + int subdet = -999; + int component = -999; DetId a(detid); - if ( a.subdetId() == 3 ){ - subdet=0; - component=tTopo->tibLayer(detid); - } else if ( a.subdetId() == 4 ) { - subdet=1; - component=tTopo->tidSide(detid)==2?tTopo->tidWheel(detid):tTopo->tidWheel(detid)+3; - } else if ( a.subdetId() == 5 ) { - subdet=2; - component=tTopo->tobLayer(detid); - } else if ( a.subdetId() == 6 ) { - subdet=3; - component=tTopo->tecSide(detid)==2?tTopo->tecWheel(detid):tTopo->tecWheel(detid)+9; - } - - SiStripQuality::Range sqrange = SiStripQuality::Range( quality_->getDataVectorBegin()+rp->ibegin , quality_->getDataVectorBegin()+rp->iend ); - - percentage=0; - for(int it=0;itdecode( *(sqrange.first+it) ).range; - NTkBadComponent[3]+=range; - NBadComponent[subdet][0][3]+=range; - NBadComponent[subdet][component][3]+=range; - percentage+=range; + if (a.subdetId() == 3) { + subdet = 0; + component = tTopo.tibLayer(detid); + } else if (a.subdetId() == 4) { + subdet = 1; + component = tTopo.tidSide(detid) == 2 ? tTopo.tidWheel(detid) : tTopo.tidWheel(detid) + 3; + } else if (a.subdetId() == 5) { + subdet = 2; + component = tTopo.tobLayer(detid); + } else if (a.subdetId() == 6) { + subdet = 3; + component = tTopo.tecSide(detid) == 2 ? tTopo.tecWheel(detid) : tTopo.tecWheel(detid) + 9; + } + + SiStripQuality::Range sqrange = + SiStripQuality::Range(quality_->getDataVectorBegin() + rp->ibegin, quality_->getDataVectorBegin() + rp->iend); + + percentage = 0; + for (int it = 0; it < sqrange.second - sqrange.first; it++) { + unsigned int range = quality_->decode(*(sqrange.first + it)).range; + NTkBadComponent[3] += range; + NBadComponent[subdet][0][3] += range; + NBadComponent[subdet][component][3] += range; + percentage += range; } - if(percentage!=0) - percentage/=128.*reader->getNumberOfApvsAndStripLength(detid).first; - if(percentage>1) - edm::LogError("SiStripQualityStatistics") << "PROBLEM detid " << detid << " value " << percentage<< std::endl; + if (percentage != 0) + percentage /= 128. * _detInfo.getNumberOfApvsAndStripLength(detid).first; + if (percentage > 1) + edm::LogError("SiStripQualityStatistics") << "PROBLEM detid " << detid << " value " << percentage << std::endl; } //&&&&&&&&&&&&&&&&&& // printout //&&&&&&&&&&&&&&&&&& - std::ofstream ResolutionValues; - int RunNumInt = e.id().run(); - std::string RunNumString = std::to_string(RunNumInt); - std::cout << "RunNumString" << RunNumString << std::endl; - std::string ResolutionTextFileString = "ResolutionValues_" + RunNumString + "_" + FileNameEnding + ".txt"; - - - ResolutionValues.open(ResolutionTextFileString.c_str()); - - for(Long_t ilayer = 0; ilayer <23; ilayer++) { + std::ofstream ResolutionValues; + int RunNumInt = e.id().run(); + std::string RunNumString = std::to_string(RunNumInt); + std::cout << "RunNumString" << RunNumString << std::endl; + std::string ResolutionTextFileString = "ResolutionValues_" + RunNumString + "_" + FileNameEnding + ".txt"; - //Calculating and printing out the resolution values + ResolutionValues.open(ResolutionTextFileString.c_str()); - float Meas = MeasPlots_p[ilayer]->GetStdDev(); - float Pred = PredPlots_p[ilayer]->GetStdDev(); + for (Long_t ilayer = 0; ilayer < 23; ilayer++) { + //Calculating and printing out the resolution values - float PredMinusMeas = pow(Meas, 2) + pow(Pred, 2); //width^2= sigma(deltaX_pred)^2 + sigma(deltaX_hit)^2 + float Meas = MeasPlots_p[ilayer]->GetStdDev(); + float Pred = PredPlots_p[ilayer]->GetStdDev(); - float Resolution = sqrt( Pred / 2 ); + float PredMinusMeas = pow(Meas, 2) + pow(Pred, 2); //width^2= sigma(deltaX_pred)^2 + sigma(deltaX_hit)^2 + float Resolution = sqrt(Pred / 2); - //Saving the resolution values to a text file - ResolutionValues << '\n' << "Resolution for layer number " << ilayer << " ("<< GetLayerName(ilayer) << ")" << " is: " << Resolution << '\n' - << "Double difference of the measured and predicted position between the two sensors under consideration for layer number " << ilayer << " ("<< GetLayerName(ilayer) << ")" << " is: " << PredMinusMeas << '\n' - << "The difference between the two positions of the hit measurements for layer number " << ilayer << " ("<< GetLayerName(ilayer) << ")" << " is: " << Meas << '\n' << std::endl; + //Saving the resolution values to a text file + ResolutionValues << '\n' + << "Resolution for layer number " << ilayer << " (" << GetLayerName(ilayer) << ")" + << " is: " << Resolution << '\n' + << "Double difference of the measured and predicted position between the two sensors under " + "consideration for layer number " + << ilayer << " (" << GetLayerName(ilayer) << ")" + << " is: " << PredMinusMeas << '\n' + << "The difference between the two positions of the hit measurements for layer number " << ilayer + << " (" << GetLayerName(ilayer) << ")" + << " is: " << Meas << '\n' + << std::endl; + } - } - - cout << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event() << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n"; + cout << "\n-----------------\nNew IOV starting from run " << e.id().run() << " event " << e.id().event() + << " lumiBlock " << e.luminosityBlock() << " time " << e.time().value() << "\n-----------------\n"; cout << "\n-----------------\nGlobal Info\n-----------------"; - cout << "\nBadComponent \t Modules \tFibers \tApvs\tStrips\n----------------------------------------------------------------"; - cout << "\nTracker:\t\t"<cd(); gStyle->SetPalette(1); gStyle->SetCanvasColor(kWhite); gStyle->SetOptStat(0); //Here we make the hot/cold color maps that we love so very much //Already have access to the data as a private variable - //Create all of the histograms in the TFileService - TH2F *temph2; - for(Long_t maplayer = 1; maplayer <=22; maplayer++) { - //Initialize all of the histograms - if(maplayer > 0 && maplayer <= 4) { - //We are in the TIB - temph2 = fs->make(Form("%s%i","TIB",(int)(maplayer)),"TIB",100,-1,361,100,-100,100); - temph2->GetXaxis()->SetTitle("Phi"); - temph2->GetXaxis()->SetBinLabel(1,TString("360")); - temph2->GetXaxis()->SetBinLabel(50,TString("180")); - temph2->GetXaxis()->SetBinLabel(100,TString("0")); - temph2->GetYaxis()->SetTitle("Global Z"); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); - } - else if(maplayer > 4 && maplayer <= 10) { - //We are in the TOB - temph2 = fs->make(Form("%s%i","TOB",(int)(maplayer-4)),"TOB",100,-1,361,100,-120,120); - temph2->GetXaxis()->SetTitle("Phi"); - temph2->GetXaxis()->SetBinLabel(1,TString("360")); - temph2->GetXaxis()->SetBinLabel(50,TString("180")); - temph2->GetXaxis()->SetBinLabel(100,TString("0")); - temph2->GetYaxis()->SetTitle("Global Z"); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); - } - else if(maplayer > 10 && maplayer <= 13) { - //We are in the TID - //Split by +/- - temph2 = fs->make(Form("%s%i","TID-",(int)(maplayer-10)),"TID-",100,-100,100,100,-100,100); - temph2->GetXaxis()->SetTitle("Global Y"); - temph2->GetXaxis()->SetBinLabel(1,TString("+Y")); - temph2->GetXaxis()->SetBinLabel(50,TString("0")); - temph2->GetXaxis()->SetBinLabel(100,TString("-Y")); - temph2->GetYaxis()->SetTitle("Global X"); - temph2->GetYaxis()->SetBinLabel(1,TString("-X")); - temph2->GetYaxis()->SetBinLabel(50,TString("0")); - temph2->GetYaxis()->SetBinLabel(100,TString("+X")); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); - temph2 = fs->make(Form("%s%i","TID+",(int)(maplayer-10)),"TID+",100,-100,100,100,-100,100); - temph2->GetXaxis()->SetTitle("Global Y"); - temph2->GetXaxis()->SetBinLabel(1,TString("+Y")); - temph2->GetXaxis()->SetBinLabel(50,TString("0")); - temph2->GetXaxis()->SetBinLabel(100,TString("-Y")); - temph2->GetYaxis()->SetTitle("Global X"); - temph2->GetYaxis()->SetBinLabel(1,TString("-X")); - temph2->GetYaxis()->SetBinLabel(50,TString("0")); - temph2->GetYaxis()->SetBinLabel(100,TString("+X")); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); - } - else if(maplayer > 13) { - //We are in the TEC - //Split by +/- - temph2 = fs->make(Form("%s%i","TEC-",(int)(maplayer-13)),"TEC-",100,-120,120,100,-120,120); - temph2->GetXaxis()->SetTitle("Global Y"); - temph2->GetXaxis()->SetBinLabel(1,TString("+Y")); - temph2->GetXaxis()->SetBinLabel(50,TString("0")); - temph2->GetXaxis()->SetBinLabel(100,TString("-Y")); - temph2->GetYaxis()->SetTitle("Global X"); - temph2->GetYaxis()->SetBinLabel(1,TString("-X")); - temph2->GetYaxis()->SetBinLabel(50,TString("0")); - temph2->GetYaxis()->SetBinLabel(100,TString("+X")); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); - temph2 = fs->make(Form("%s%i","TEC+",(int)(maplayer-13)),"TEC+",100,-120,120,100,-120,120); - temph2->GetXaxis()->SetTitle("Global Y"); - temph2->GetXaxis()->SetBinLabel(1,TString("+Y")); - temph2->GetXaxis()->SetBinLabel(50,TString("0")); - temph2->GetXaxis()->SetBinLabel(100,TString("-Y")); - temph2->GetYaxis()->SetTitle("Global X"); - temph2->GetYaxis()->SetBinLabel(1,TString("-X")); - temph2->GetYaxis()->SetBinLabel(50,TString("0")); - temph2->GetYaxis()->SetBinLabel(100,TString("+X")); - temph2->SetOption("colz"); - HotColdMaps.push_back(temph2); + //Create all of the histograms in the TFileService + TH2F* temph2; + for (Long_t maplayer = 1; maplayer <= 22; maplayer++) { + //Initialize all of the histograms + if (maplayer > 0 && maplayer <= 4) { + //We are in the TIB + temph2 = fs->make(Form("%s%i", "TIB", (int)(maplayer)), "TIB", 100, -1, 361, 100, -100, 100); + temph2->GetXaxis()->SetTitle("Phi"); + temph2->GetXaxis()->SetBinLabel(1, TString("360")); + temph2->GetXaxis()->SetBinLabel(50, TString("180")); + temph2->GetXaxis()->SetBinLabel(100, TString("0")); + temph2->GetYaxis()->SetTitle("Global Z"); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); + } else if (maplayer > 4 && maplayer <= 10) { + //We are in the TOB + temph2 = fs->make(Form("%s%i", "TOB", (int)(maplayer - 4)), "TOB", 100, -1, 361, 100, -120, 120); + temph2->GetXaxis()->SetTitle("Phi"); + temph2->GetXaxis()->SetBinLabel(1, TString("360")); + temph2->GetXaxis()->SetBinLabel(50, TString("180")); + temph2->GetXaxis()->SetBinLabel(100, TString("0")); + temph2->GetYaxis()->SetTitle("Global Z"); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); + } else if (maplayer > 10 && maplayer <= 13) { + //We are in the TID + //Split by +/- + temph2 = fs->make(Form("%s%i", "TID-", (int)(maplayer - 10)), "TID-", 100, -100, 100, 100, -100, 100); + temph2->GetXaxis()->SetTitle("Global Y"); + temph2->GetXaxis()->SetBinLabel(1, TString("+Y")); + temph2->GetXaxis()->SetBinLabel(50, TString("0")); + temph2->GetXaxis()->SetBinLabel(100, TString("-Y")); + temph2->GetYaxis()->SetTitle("Global X"); + temph2->GetYaxis()->SetBinLabel(1, TString("-X")); + temph2->GetYaxis()->SetBinLabel(50, TString("0")); + temph2->GetYaxis()->SetBinLabel(100, TString("+X")); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); + temph2 = fs->make(Form("%s%i", "TID+", (int)(maplayer - 10)), "TID+", 100, -100, 100, 100, -100, 100); + temph2->GetXaxis()->SetTitle("Global Y"); + temph2->GetXaxis()->SetBinLabel(1, TString("+Y")); + temph2->GetXaxis()->SetBinLabel(50, TString("0")); + temph2->GetXaxis()->SetBinLabel(100, TString("-Y")); + temph2->GetYaxis()->SetTitle("Global X"); + temph2->GetYaxis()->SetBinLabel(1, TString("-X")); + temph2->GetYaxis()->SetBinLabel(50, TString("0")); + temph2->GetYaxis()->SetBinLabel(100, TString("+X")); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); + } else if (maplayer > 13) { + //We are in the TEC + //Split by +/- + temph2 = fs->make(Form("%s%i", "TEC-", (int)(maplayer - 13)), "TEC-", 100, -120, 120, 100, -120, 120); + temph2->GetXaxis()->SetTitle("Global Y"); + temph2->GetXaxis()->SetBinLabel(1, TString("+Y")); + temph2->GetXaxis()->SetBinLabel(50, TString("0")); + temph2->GetXaxis()->SetBinLabel(100, TString("-Y")); + temph2->GetYaxis()->SetTitle("Global X"); + temph2->GetYaxis()->SetBinLabel(1, TString("-X")); + temph2->GetYaxis()->SetBinLabel(50, TString("0")); + temph2->GetYaxis()->SetBinLabel(100, TString("+X")); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); + temph2 = fs->make(Form("%s%i", "TEC+", (int)(maplayer - 13)), "TEC+", 100, -120, 120, 100, -120, 120); + temph2->GetXaxis()->SetTitle("Global Y"); + temph2->GetXaxis()->SetBinLabel(1, TString("+Y")); + temph2->GetXaxis()->SetBinLabel(50, TString("0")); + temph2->GetXaxis()->SetBinLabel(100, TString("-Y")); + temph2->GetYaxis()->SetTitle("Global X"); + temph2->GetYaxis()->SetBinLabel(1, TString("-X")); + temph2->GetYaxis()->SetBinLabel(50, TString("0")); + temph2->GetYaxis()->SetBinLabel(100, TString("+X")); + temph2->SetOption("colz"); + HotColdMaps.push_back(temph2); } } - for(Long_t mylayer = 1; mylayer <= 22; mylayer++) { + for (Long_t mylayer = 1; mylayer <= 22; mylayer++) { //Determine what kind of plot we want to write out //Loop through the entirety of each layer //Create an array of the histograms vector::const_iterator iter; - for(iter = hits[mylayer].begin(); iter != hits[mylayer].end(); iter++) { + for (iter = hits[mylayer].begin(); iter != hits[mylayer].end(); iter++) { //Looping over the particular layer //Fill by 360-x to get the proper location to compare with TKMaps of phi //Also global xy is messed up - if(mylayer > 0 && mylayer <= 4) { + if (mylayer > 0 && mylayer <= 4) { //We are in the TIB float phi = calcPhi(iter->x, iter->y); - HotColdMaps[mylayer - 1]->Fill(360.-phi,iter->z,1.); - } - else if(mylayer > 4 && mylayer <= 10) { + HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.); + } else if (mylayer > 4 && mylayer <= 10) { //We are in the TOB - float phi = calcPhi(iter->x,iter->y); - HotColdMaps[mylayer - 1]->Fill(360.-phi,iter->z,1.); - } - else if(mylayer > 10 && mylayer <= 13) { + float phi = calcPhi(iter->x, iter->y); + HotColdMaps[mylayer - 1]->Fill(360. - phi, iter->z, 1.); + } else if (mylayer > 10 && mylayer <= 13) { //We are in the TID //There are 2 different maps here - int side = (((iter->id)>>13) & 0x3); - if(side == 1) HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(-iter->y,iter->x,1.); - else if(side == 2) HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(-iter->y,iter->x,1.); - //if(side == 1) HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(iter->x,iter->y,1.); + int side = (((iter->id) >> 13) & 0x3); + if (side == 1) + HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(-iter->y, iter->x, 1.); + else if (side == 2) + HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(-iter->y, iter->x, 1.); + //if(side == 1) HotColdMaps[(mylayer - 1) + (mylayer - 11)]->Fill(iter->x,iter->y,1.); //else if(side == 2) HotColdMaps[(mylayer - 1) + (mylayer - 10)]->Fill(iter->x,iter->y,1.); - } - else if(mylayer > 13) { + } else if (mylayer > 13) { //We are in the TEC //There are 2 different maps here - int side = (((iter->id)>>18) & 0x3); - if(side == 1) HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(-iter->y,iter->x,1.); - else if(side == 2) HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(-iter->y,iter->x,1.); + int side = (((iter->id) >> 18) & 0x3); + if (side == 1) + HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(-iter->y, iter->x, 1.); + else if (side == 2) + HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(-iter->y, iter->x, 1.); //if(side == 1) HotColdMaps[(mylayer + 2) + (mylayer - 14)]->Fill(iter->x,iter->y,1.); //else if(side == 2) HotColdMaps[(mylayer + 2) + (mylayer - 13)]->Fill(iter->x,iter->y,1.); } @@ -1051,94 +1092,96 @@ void SiStripHitResolFromCalibTree::makeHotColdMaps() { cout << "Finished HotCold Map Generation\n"; } -void SiStripHitResolFromCalibTree::makeTKMap(bool autoTagging=false) { +void SiStripHitResolFromCalibTree::makeTKMap(bool autoTagging = false) { cout << "Entering TKMap generation!\n"; tkmap = new TrackerMap(" Detector Inefficiency "); tkmapbad = new TrackerMap(" Inefficient Modules "); tkmapeff = new TrackerMap(_title.Data()); tkmapnum = new TrackerMap(" Detector numerator "); tkmapden = new TrackerMap(" Detector denominator "); - + double myeff, mynum, myden; - double eff_limit=0; - - for(Long_t i = 1; i <= 22; i++) { - //Loop over every layer, extracting the information from - //the map of the efficiencies - layertotal[i] = 0; - layerfound[i] = 0; - TH1F* hEffInLayer = fs->make(Form("eff_layer%i", int(i)),Form("Module efficiency in layer %i", int(i)), 201,0,1.005); - + double eff_limit = 0; + + for (Long_t i = 1; i <= 22; i++) { + //Loop over every layer, extracting the information from + //the map of the efficiencies + layertotal[i] = 0; + layerfound[i] = 0; + TH1F* hEffInLayer = + fs->make(Form("eff_layer%i", int(i)), Form("Module efficiency in layer %i", int(i)), 201, 0, 1.005); + map >::const_iterator ih; - for( ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) { + for (ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) { //We should be in the layer in question, and looping over all of the modules in said layer //Generate the list for the TKmap, and the bad module list - mynum = (double)(((*ih).second).second); - myden = (double)(((*ih).second).first); - if(myden>0) myeff = mynum/myden; - else myeff=0; - hEffInLayer->Fill(myeff); - - if(!autoTagging) { - if ( (myden >= nModsMin) && (myeff < threshold) ) { + mynum = (double)(((*ih).second).second); + myden = (double)(((*ih).second).first); + if (myden > 0) + myeff = mynum / myden; + else + myeff = 0; + hEffInLayer->Fill(myeff); + + if (!autoTagging) { + if ((myden >= nModsMin) && (myeff < threshold)) { //We have a bad module, put it in the list! - BadModules[(*ih).first] = myeff; - tkmapbad->fillc((*ih).first,255,0,0); - cout << "Layer " << i <<" ("<< GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff << " , " << mynum << "/" << myden << endl; - } - else { + BadModules[(*ih).first] = myeff; + tkmapbad->fillc((*ih).first, 255, 0, 0); + cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff + << " , " << mynum << "/" << myden << endl; + } else { //Fill the bad list with empty results for every module - tkmapbad->fillc((*ih).first,255,255,255); - } - if(myden < nModsMin ) { - cout << "Layer " << i <<" ("<< GetLayerName(i) << ") module " << (*ih).first << " is under occupancy at " << myden << endl; - } + tkmapbad->fillc((*ih).first, 255, 255, 255); + } + if (myden < nModsMin) { + cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " is under occupancy at " + << myden << endl; + } } - + //Put any module into the TKMap - tkmap->fill((*ih).first,1.-myeff); - tkmapeff->fill((*ih).first,myeff); - tkmapnum->fill((*ih).first,mynum); - tkmapden->fill((*ih).first,myden); - + tkmap->fill((*ih).first, 1. - myeff); + tkmapeff->fill((*ih).first, myeff); + tkmapnum->fill((*ih).first, mynum); + tkmapden->fill((*ih).first, myden); + //Add the number of hits in the layer layertotal[i] += long(myden); layerfound[i] += long(mynum); } - - if(autoTagging) { - - //Compute threshold to use for each layer - hEffInLayer->GetXaxis()->SetRange(3, hEffInLayer->GetNbinsX()+1); // Remove from the avg modules below 1% - eff_limit = hEffInLayer->GetMean()-threshold; - cout << "Layer " << i << " threshold for bad modules: " << eff_limit << endl; - hEffInLayer->GetXaxis()->SetRange(1, hEffInLayer->GetNbinsX()+1); - - - for( ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) { - // Second loop over modules to tag inefficient ones - mynum = (double)(((*ih).second).second); - myden = (double)(((*ih).second).first); - if(myden>0) myeff = mynum/myden; - else myeff=0; - if ( (myden >= nModsMin) && (myeff < eff_limit) ) { + + if (autoTagging) { + //Compute threshold to use for each layer + hEffInLayer->GetXaxis()->SetRange(3, hEffInLayer->GetNbinsX() + 1); // Remove from the avg modules below 1% + eff_limit = hEffInLayer->GetMean() - threshold; + cout << "Layer " << i << " threshold for bad modules: " << eff_limit << endl; + hEffInLayer->GetXaxis()->SetRange(1, hEffInLayer->GetNbinsX() + 1); + + for (ih = modCounter[i].begin(); ih != modCounter[i].end(); ih++) { + // Second loop over modules to tag inefficient ones + mynum = (double)(((*ih).second).second); + myden = (double)(((*ih).second).first); + if (myden > 0) + myeff = mynum / myden; + else + myeff = 0; + if ((myden >= nModsMin) && (myeff < eff_limit)) { //We have a bad module, put it in the list! - BadModules[(*ih).first] = myeff; - tkmapbad->fillc((*ih).first,255,0,0); - cout << "Layer " << i <<" ("<< GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff << " , " << mynum << "/" << myden << endl; - } - else { + BadModules[(*ih).first] = myeff; + tkmapbad->fillc((*ih).first, 255, 0, 0); + cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " efficiency: " << myeff + << " , " << mynum << "/" << myden << endl; + } else { //Fill the bad list with empty results for every module - tkmapbad->fillc((*ih).first,255,255,255); - } - if(myden < nModsMin ) { - cout << "Layer " << i <<" ("<< GetLayerName(i) << ") module " << (*ih).first << " layer " << i << " is under occupancy at " << myden << endl; - } + tkmapbad->fillc((*ih).first, 255, 255, 255); + } + if (myden < nModsMin) { + cout << "Layer " << i << " (" << GetLayerName(i) << ") module " << (*ih).first << " layer " << i + << " is under occupancy at " << myden << endl; + } } - - } - - + } } tkmap->save(true, 0, 0, "SiStripHitResolTKMap.png"); tkmapbad->save(true, 0, 0, "SiStripHitResolTKMapBad.png"); @@ -1148,30 +1191,29 @@ void SiStripHitResolFromCalibTree::makeTKMap(bool autoTagging=false) { cout << "Finished TKMap Generation\n"; } - void SiStripHitResolFromCalibTree::makeSQLite() { //Generate the SQLite file for use in the Database of the bad modules! cout << "Entering SQLite file generation!\n"; std::vector BadStripList; unsigned short NStrips; unsigned int id1; - SiStripQuality* pQuality = new SiStripQuality; + std::unique_ptr pQuality = std::make_unique(_detInfo); //This is the list of the bad strips, use to mask out entire APVs //Now simply go through the bad hit list and mask out things that //are bad! - map< unsigned int, double >::const_iterator it; - for(it = BadModules.begin(); it != BadModules.end(); it++) { + map::const_iterator it; + for (it = BadModules.begin(); it != BadModules.end(); it++) { //We need to figure out how many strips are in this particular module //To Mask correctly! - NStrips=reader->getNumberOfApvsAndStripLength((*it).first).first*128; + NStrips = _detInfo.getNumberOfApvsAndStripLength((*it).first).first * 128; cout << "Number of strips module " << (*it).first << " is " << NStrips << endl; - BadStripList.push_back(pQuality->encode(0,NStrips,0)); + BadStripList.push_back(pQuality->encode(0, NStrips, 0)); //Now compact into a single bad module - id1=(unsigned int)(*it).first; + id1 = (unsigned int)(*it).first; cout << "ID1 shoudl match list of modules above " << id1 << endl; - quality_->compact(id1,BadStripList); - SiStripQuality::Range range(BadStripList.begin(),BadStripList.end()); - quality_->put(id1,range); + quality_->compact(id1, BadStripList); + SiStripQuality::Range range(BadStripList.begin(), BadStripList.end()); + quality_->put(id1, range); BadStripList.clear(); } //Fill all the bad components now @@ -1185,101 +1227,130 @@ void SiStripHitResolFromCalibTree::totalStatistics() { double layereff; int subdetfound[5]; int subdettotal[5]; - - for(Long_t i=1; i<5; i++) {subdetfound[i]=0; subdettotal[i]=0;} - - for(Long_t i=1; i<=22; i++) { - layereff = double(layerfound[i])/double(layertotal[i]); - cout << "Layer " << i << " (" << GetLayerName(i) << ") has total efficiency " << layereff << " " << layerfound[i] << "/" << layertotal[i] << endl; + + for (Long_t i = 1; i < 5; i++) { + subdetfound[i] = 0; + subdettotal[i] = 0; + } + + for (Long_t i = 1; i <= 22; i++) { + layereff = double(layerfound[i]) / double(layertotal[i]); + cout << "Layer " << i << " (" << GetLayerName(i) << ") has total efficiency " << layereff << " " << layerfound[i] + << "/" << layertotal[i] << endl; totalfound += layerfound[i]; totaltotal += layertotal[i]; - if(i<5) {subdetfound[1]+=layerfound[i]; subdettotal[1]+=layertotal[i];} - if(i>=5 && i<11) {subdetfound[2]+=layerfound[i]; subdettotal[2]+=layertotal[i];} - if(i>=11 && i<14) {subdetfound[3]+=layerfound[i]; subdettotal[3]+=layertotal[i];} - if(i>=14) {subdetfound[4]+=layerfound[i]; subdettotal[4]+=layertotal[i];} - + if (i < 5) { + subdetfound[1] += layerfound[i]; + subdettotal[1] += layertotal[i]; + } + if (i >= 5 && i < 11) { + subdetfound[2] += layerfound[i]; + subdettotal[2] += layertotal[i]; + } + if (i >= 11 && i < 14) { + subdetfound[3] += layerfound[i]; + subdettotal[3] += layertotal[i]; + } + if (i >= 14) { + subdetfound[4] += layerfound[i]; + subdettotal[4] += layertotal[i]; + } } - - cout << "The total efficiency is " << double(totalfound)/double(totaltotal) << endl; - cout << " TIB: " << double(subdetfound[1])/subdettotal[1] <<" "<< subdetfound[1]<<"/"<make("found","found",nLayers+1,0,nLayers+1); - TH1F *all = fs->make("all","all",nLayers+1,0,nLayers+1); - TH1F *found2 = fs->make("found2","found2",nLayers+1,0,nLayers+1); - TH1F *all2 = fs->make("all2","all2",nLayers+1,0,nLayers+1); + + TH1F* found = fs->make("found", "found", nLayers + 1, 0, nLayers + 1); + TH1F* all = fs->make("all", "all", nLayers + 1, 0, nLayers + 1); + TH1F* found2 = fs->make("found2", "found2", nLayers + 1, 0, nLayers + 1); + TH1F* all2 = fs->make("all2", "all2", nLayers + 1, 0, nLayers + 1); // first bin only to keep real data off the y axis so set to -1 - found->SetBinContent(0,-1); - all->SetBinContent(0,1); - + found->SetBinContent(0, -1); + all->SetBinContent(0, 1); + // new ROOT version: TGraph::Divide don't handle null or negative values - for (Long_t i=1; i< nLayers+2; ++i) { - found->SetBinContent(i,1e-6); - all->SetBinContent(i,1); - found2->SetBinContent(i,1e-6); - all2->SetBinContent(i,1); + for (Long_t i = 1; i < nLayers + 2; ++i) { + found->SetBinContent(i, 1e-6); + all->SetBinContent(i, 1); + found2->SetBinContent(i, 1e-6); + all2->SetBinContent(i, 1); } - - TCanvas *c7 =new TCanvas("c7"," test ",10,10,800,600); + + TCanvas* c7 = new TCanvas("c7", " test ", 10, 10, 800, 600); c7->SetFillColor(0); c7->SetGrid(); - int nLayers_max=nLayers+1; // barrel+endcap - if(!_showEndcapSides) nLayers_max=11; // barrel - for (Long_t i=1; i< nLayers_max; ++i) { - cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] << " B = " << goodlayertotal[i] << endl; + int nLayers_max = nLayers + 1; // barrel+endcap + if (!_showEndcapSides) + nLayers_max = 11; // barrel + for (Long_t i = 1; i < nLayers_max; ++i) { + cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] << " B = " << goodlayertotal[i] + << endl; if (goodlayertotal[i] > 5) { - found->SetBinContent(i,goodlayerfound[i]); - all->SetBinContent(i,goodlayertotal[i]); + found->SetBinContent(i, goodlayerfound[i]); + all->SetBinContent(i, goodlayertotal[i]); } - - cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i] << " B = " << alllayertotal[i] << endl; + + cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i] << " B = " << alllayertotal[i] + << endl; if (alllayertotal[i] > 5) { - found2->SetBinContent(i,alllayerfound[i]); - all2->SetBinContent(i,alllayertotal[i]); + found2->SetBinContent(i, alllayerfound[i]); + all2->SetBinContent(i, alllayertotal[i]); } - } // endcap - merging sides - if(!_showEndcapSides) { - for (Long_t i=11; i< 14; ++i) { // TID disks - cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i]+goodlayerfound[i+3] << " B = " << goodlayertotal[i]+goodlayertotal[i+3] << endl; - if (goodlayertotal[i]+goodlayertotal[i+3] > 5) { - found->SetBinContent(i,goodlayerfound[i]+goodlayerfound[i+3]); - all->SetBinContent(i,goodlayertotal[i]+goodlayertotal[i+3]); + if (!_showEndcapSides) { + for (Long_t i = 11; i < 14; ++i) { // TID disks + cout << "Fill only good modules layer " << i << ": S = " << goodlayerfound[i] + goodlayerfound[i + 3] + << " B = " << goodlayertotal[i] + goodlayertotal[i + 3] << endl; + if (goodlayertotal[i] + goodlayertotal[i + 3] > 5) { + found->SetBinContent(i, goodlayerfound[i] + goodlayerfound[i + 3]); + all->SetBinContent(i, goodlayertotal[i] + goodlayertotal[i + 3]); } - cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i]+alllayerfound[i+3] << " B = " << alllayertotal[i]+alllayertotal[i+3] << endl; - if (alllayertotal[i]+alllayertotal[i+3] > 5) { - found2->SetBinContent(i,alllayerfound[i]+alllayerfound[i+3]); - all2->SetBinContent(i,alllayertotal[i]+alllayertotal[i+3]); + cout << "Filling all modules layer " << i << ": S = " << alllayerfound[i] + alllayerfound[i + 3] + << " B = " << alllayertotal[i] + alllayertotal[i + 3] << endl; + if (alllayertotal[i] + alllayertotal[i + 3] > 5) { + found2->SetBinContent(i, alllayerfound[i] + alllayerfound[i + 3]); + all2->SetBinContent(i, alllayertotal[i] + alllayertotal[i + 3]); } - } - for (Long_t i=17; i< 17+nTEClayers; ++i) { // TEC disks - cout << "Fill only good modules layer " << i-3 << ": S = " << goodlayerfound[i]+goodlayerfound[i+nTEClayers] << " B = " << goodlayertotal[i]+goodlayertotal[i+nTEClayers] << endl; - if (goodlayertotal[i]+goodlayertotal[i+nTEClayers] > 5) { - found->SetBinContent(i-3,goodlayerfound[i]+goodlayerfound[i+nTEClayers]); - all->SetBinContent(i-3,goodlayertotal[i]+goodlayertotal[i+nTEClayers]); + } + for (Long_t i = 17; i < 17 + nTEClayers; ++i) { // TEC disks + cout << "Fill only good modules layer " << i - 3 + << ": S = " << goodlayerfound[i] + goodlayerfound[i + nTEClayers] + << " B = " << goodlayertotal[i] + goodlayertotal[i + nTEClayers] << endl; + if (goodlayertotal[i] + goodlayertotal[i + nTEClayers] > 5) { + found->SetBinContent(i - 3, goodlayerfound[i] + goodlayerfound[i + nTEClayers]); + all->SetBinContent(i - 3, goodlayertotal[i] + goodlayertotal[i + nTEClayers]); } - cout << "Filling all modules layer " << i-3 << ": S = " << alllayerfound[i]+alllayerfound[i+nTEClayers] << " B = " << alllayertotal[i]+alllayertotal[i+nTEClayers] << endl; - if (alllayertotal[i]+alllayertotal[i+nTEClayers] > 5) { - found2->SetBinContent(i-3,alllayerfound[i]+alllayerfound[i+nTEClayers]); - all2->SetBinContent(i-3,alllayertotal[i]+alllayertotal[i+nTEClayers]); + cout << "Filling all modules layer " << i - 3 << ": S = " << alllayerfound[i] + alllayerfound[i + nTEClayers] + << " B = " << alllayertotal[i] + alllayertotal[i + nTEClayers] << endl; + if (alllayertotal[i] + alllayertotal[i + nTEClayers] > 5) { + found2->SetBinContent(i - 3, alllayerfound[i] + alllayerfound[i + nTEClayers]); + all2->SetBinContent(i - 3, alllayertotal[i] + alllayertotal[i + nTEClayers]); } - } + } } found->Sumw2(); @@ -1288,20 +1359,20 @@ void SiStripHitResolFromCalibTree::makeSummary() { found2->Sumw2(); all2->Sumw2(); - TGraphAsymmErrors *gr = fs->make(nLayers+1); + TGraphAsymmErrors* gr = fs->make(nLayers + 1); gr->SetName("eff_good"); - gr->BayesDivide(found,all); + gr->BayesDivide(found, all); - TGraphAsymmErrors *gr2 = fs->make(nLayers+1); + TGraphAsymmErrors* gr2 = fs->make(nLayers + 1); gr2->SetName("eff_all"); - gr2->BayesDivide(found2,all2); + gr2->BayesDivide(found2, all2); - for(int j = 0; jSetPointError(j, 0., 0., gr->GetErrorYlow(j),gr->GetErrorYhigh(j) ); - gr2->SetPointError(j, 0., 0., gr2->GetErrorYlow(j),gr2->GetErrorYhigh(j) ); + for (int j = 0; j < nLayers + 1; j++) { + gr->SetPointError(j, 0., 0., gr->GetErrorYlow(j), gr->GetErrorYhigh(j)); + gr2->SetPointError(j, 0., 0., gr2->GetErrorYlow(j), gr2->GetErrorYhigh(j)); } - gr->GetXaxis()->SetLimits(0,nLayers); + gr->GetXaxis()->SetLimits(0, nLayers); gr->SetMarkerColor(2); gr->SetMarkerSize(1.2); gr->SetLineColor(2); @@ -1314,7 +1385,7 @@ void SiStripHitResolFromCalibTree::makeSummary() { gStyle->SetTitleBorderSize(0); gr->SetTitle(_title); - gr2->GetXaxis()->SetLimits(0,nLayers); + gr2->GetXaxis()->SetLimits(0, nLayers); gr2->SetMarkerColor(1); gr2->SetMarkerSize(1.2); gr2->SetLineColor(1); @@ -1325,322 +1396,327 @@ void SiStripHitResolFromCalibTree::makeSummary() { gr2->GetYaxis()->SetTitle("Efficiency"); gr2->SetTitle(_title); - for ( Long_t k=1; kGetXaxis()->SetBinLabel(((k+1)*100+2)/(nLayers)-4,label); - gr2->GetXaxis()->SetBinLabel(((k+1)*100+2)/(nLayers)-4,label); - } - else { - gr->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-6,label); - gr2->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-6,label); - } - } - else { - if(_showEndcapSides) { - gr->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-4,label); - gr2->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-4,label); - } - else { - gr->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-7,label); - gr2->GetXaxis()->SetBinLabel((k+1)*100/(nLayers)-7,label); - } - } + for (Long_t k = 1; k < nLayers + 1; k++) { + TString label; + if (_showEndcapSides) + label = GetLayerSideName(k); + else + label = GetLayerName(k); + if (!_showTOB6TEC9) { + if (k == 10) + label = ""; + if (!_showRings && k == nLayers) + label = ""; + if (!_showRings && _showEndcapSides && k == 25) + label = ""; + } + if (!_showRings) { + if (_showEndcapSides) { + gr->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label); + gr2->GetXaxis()->SetBinLabel(((k + 1) * 100 + 2) / (nLayers)-4, label); + } else { + gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label); + gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-6, label); + } + } else { + if (_showEndcapSides) { + gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label); + gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-4, label); + } else { + gr->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label); + gr2->GetXaxis()->SetBinLabel((k + 1) * 100 / (nLayers)-7, label); + } + } } - + gr->Draw("AP"); gr->GetXaxis()->SetNdivisions(36); c7->cd(); - TPad *overlay = new TPad("overlay","",0,0,1,1); + TPad* overlay = new TPad("overlay", "", 0, 0, 1, 1); overlay->SetFillStyle(4000); overlay->SetFillColor(0); overlay->SetFrameFillStyle(4000); overlay->Draw("same"); overlay->cd(); - if(!_showOnlyGoodModules) gr2->Draw("AP"); + if (!_showOnlyGoodModules) + gr2->Draw("AP"); - TLegend *leg = new TLegend(0.70,0.27,0.88,0.40); - leg->AddEntry(gr,"Good Modules","p"); - if(!_showOnlyGoodModules) leg->AddEntry(gr2,"All Modules","p"); + TLegend* leg = new TLegend(0.70, 0.27, 0.88, 0.40); + leg->AddEntry(gr, "Good Modules", "p"); + if (!_showOnlyGoodModules) + leg->AddEntry(gr2, "All Modules", "p"); leg->SetTextSize(0.020); leg->SetFillColor(0); leg->Draw("same"); - + c7->SaveAs("Summary.png"); } - void SiStripHitResolFromCalibTree::makeSummaryVsBx() { - cout<<"Computing efficiency vs bx"<make(Form("foundVsBx_layer%i", ilayer),Form("layer %i", ilayer),3565,0,3565); - TH1F *htotal = fs->make(Form("totalVsBx_layer%i", ilayer),Form("layer %i", ilayer),3565,0,3565); - - for(unsigned int ibx=0; ibx<3566; ibx++){ - hfound->SetBinContent(ibx, 1e-6); - htotal->SetBinContent(ibx, 1); - } - map >::iterator iterMapvsBx; - for(iterMapvsBx=layerfound_perBx.begin(); iterMapvsBx!=layerfound_perBx.end(); ++iterMapvsBx) - hfound->SetBinContent( iterMapvsBx->first, iterMapvsBx->second[ilayer]); - for(iterMapvsBx=layertotal_perBx.begin(); iterMapvsBx!=layertotal_perBx.end(); ++iterMapvsBx) - if(iterMapvsBx->second[ilayer]>0) htotal->SetBinContent( iterMapvsBx->first, iterMapvsBx->second[ilayer]); - - hfound->Sumw2(); - htotal->Sumw2(); - - TGraphAsymmErrors *geff = fs->make(3564); - geff->SetName(Form("effVsBx_layer%i", ilayer)); - geff->SetTitle("Hit Efficiency vs bx - "+GetLayerName(ilayer)); - geff->BayesDivide(hfound,htotal); - - //Average over trains - TGraphAsymmErrors *geff_avg = fs->make(); - geff_avg->SetName(Form("effVsBxAvg_layer%i", ilayer)); - geff_avg->SetTitle("Hit Efficiency vs bx - "+GetLayerName(ilayer)); - geff_avg->SetMarkerStyle(20); - int ibx=0; - int previous_bx=-80; - int delta_bx=0; - int nbx=0; - int found=0; - int total=0; - double sum_bx=0; - int ipt=0; - float low, up, eff; - int firstbx=0; - for(iterMapvsBx=layertotal_perBx.begin(); iterMapvsBx!=layertotal_perBx.end(); ++iterMapvsBx){ - ibx=iterMapvsBx->first; - delta_bx=ibx-previous_bx; - // consider a new train - if(delta_bx>(int)_spaceBetweenTrains && nbx>0 && total>0){ - eff=found/(float)total; - //cout<<"new train "<SetPoint(ipt, sum_bx/nbx, eff); - low = TEfficiency::Bayesian(total, found, .683, 1, 1, false); - up = TEfficiency::Bayesian(total, found, .683, 1, 1, true); - geff_avg->SetPointError(ipt, sum_bx/nbx-firstbx, previous_bx-sum_bx/nbx, eff-low, up-eff); - ipt++; - sum_bx=0; - found=0; - total=0; - nbx=0; - firstbx=ibx; - } - sum_bx+=ibx; - found+=hfound->GetBinContent(ibx); - total+=htotal->GetBinContent(ibx); - nbx++; - - previous_bx=ibx; - } - //last train - eff=found/(float)total; - //cout<<"new train "<SetPoint(ipt, sum_bx/nbx, eff); - low = TEfficiency::Bayesian(total, found, .683, 1, 1, false); - up = TEfficiency::Bayesian(total, found, .683, 1, 1, true); - geff_avg->SetPointError(ipt, sum_bx/nbx-firstbx, previous_bx-sum_bx/nbx, eff-low, up-eff); + if (_showRings) + nLayers = 20; + + for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) { + TH1F* hfound = fs->make(Form("foundVsBx_layer%i", ilayer), Form("layer %i", ilayer), 3565, 0, 3565); + TH1F* htotal = fs->make(Form("totalVsBx_layer%i", ilayer), Form("layer %i", ilayer), 3565, 0, 3565); + + for (unsigned int ibx = 0; ibx < 3566; ibx++) { + hfound->SetBinContent(ibx, 1e-6); + htotal->SetBinContent(ibx, 1); + } + map >::iterator iterMapvsBx; + for (iterMapvsBx = layerfound_perBx.begin(); iterMapvsBx != layerfound_perBx.end(); ++iterMapvsBx) + hfound->SetBinContent(iterMapvsBx->first, iterMapvsBx->second[ilayer]); + for (iterMapvsBx = layertotal_perBx.begin(); iterMapvsBx != layertotal_perBx.end(); ++iterMapvsBx) + if (iterMapvsBx->second[ilayer] > 0) + htotal->SetBinContent(iterMapvsBx->first, iterMapvsBx->second[ilayer]); + + hfound->Sumw2(); + htotal->Sumw2(); + + TGraphAsymmErrors* geff = fs->make(3564); + geff->SetName(Form("effVsBx_layer%i", ilayer)); + geff->SetTitle("Hit Efficiency vs bx - " + GetLayerName(ilayer)); + geff->BayesDivide(hfound, htotal); + + //Average over trains + TGraphAsymmErrors* geff_avg = fs->make(); + geff_avg->SetName(Form("effVsBxAvg_layer%i", ilayer)); + geff_avg->SetTitle("Hit Efficiency vs bx - " + GetLayerName(ilayer)); + geff_avg->SetMarkerStyle(20); + int ibx = 0; + int previous_bx = -80; + int delta_bx = 0; + int nbx = 0; + int found = 0; + int total = 0; + double sum_bx = 0; + int ipt = 0; + float low, up, eff; + int firstbx = 0; + for (iterMapvsBx = layertotal_perBx.begin(); iterMapvsBx != layertotal_perBx.end(); ++iterMapvsBx) { + ibx = iterMapvsBx->first; + delta_bx = ibx - previous_bx; + // consider a new train + if (delta_bx > (int)_spaceBetweenTrains && nbx > 0 && total > 0) { + eff = found / (float)total; + //cout<<"new train "<SetPoint(ipt, sum_bx / nbx, eff); + low = TEfficiency::Bayesian(total, found, .683, 1, 1, false); + up = TEfficiency::Bayesian(total, found, .683, 1, 1, true); + geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff); + ipt++; + sum_bx = 0; + found = 0; + total = 0; + nbx = 0; + firstbx = ibx; + } + sum_bx += ibx; + found += hfound->GetBinContent(ibx); + total += htotal->GetBinContent(ibx); + nbx++; + + previous_bx = ibx; + } + //last train + eff = found / (float)total; + //cout<<"new train "<SetPoint(ipt, sum_bx / nbx, eff); + low = TEfficiency::Bayesian(total, found, .683, 1, 1, false); + up = TEfficiency::Bayesian(total, found, .683, 1, 1, true); + geff_avg->SetPointError(ipt, sum_bx / nbx - firstbx, previous_bx - sum_bx / nbx, eff - low, up - eff); } } - TString SiStripHitResolFromCalibTree::GetLayerName(Long_t k) { + TString layername = ""; + TString ringlabel = "D"; + if (_showRings) + ringlabel = "R"; + if (k > 0 && k < 5) { + layername = TString("TIB L") + k; + } else if (k > 4 && k < 11) { + layername = TString("TOB L") + (k - 4); + } else if (k > 10 && k < 14) { + layername = TString("TID ") + ringlabel + (k - 10); + } else if (k > 13 && k < 14 + nTEClayers) { + layername = TString("TEC ") + ringlabel + (k - 13); + } - TString layername=""; - TString ringlabel="D"; - if(_showRings) ringlabel="R"; - if (k>0 && k<5) { - layername = TString("TIB L") + k; - } else if (k>4 && k<11) { - layername = TString("TOB L")+(k-4); - } else if (k>10 && k<14) { - layername = TString("TID ")+ringlabel+(k-10); - } else if (k>13 && k<14+nTEClayers) { - layername = TString("TEC ")+ringlabel+(k-13); - } - - return layername; + return layername; } -void SiStripHitResolFromCalibTree::ComputeEff(vector< TH1F* > &vhfound, vector< TH1F* > &vhtotal, string name) { - +void SiStripHitResolFromCalibTree::ComputeEff(vector& vhfound, vector& vhtotal, string name) { unsigned int nLayers = 22; - if(_showRings) nLayers = 20; - + if (_showRings) + nLayers = 20; + TH1F* hfound; TH1F* htotal; - - for(unsigned int ilayer=1; ilayerSumw2(); - htotal->Sumw2(); - - // new ROOT version: TGraph::Divide don't handle null or negative values - for (Long_t i=0; i< hfound->GetNbinsX()+1; ++i) { - if( hfound->GetBinContent(i)==0) hfound->SetBinContent(i,1e-6); - if( htotal->GetBinContent(i)==0) htotal->SetBinContent(i,1); - } - - TGraphAsymmErrors *geff = fs->make(hfound->GetNbinsX()); - geff->SetName(Form("%s_layer%i", name.c_str(), ilayer)); - geff->BayesDivide(hfound, htotal); - if(name=="effVsLumi") geff->SetTitle("Hit Efficiency vs inst. lumi. - "+GetLayerName(ilayer)); - if(name=="effVsPU") geff->SetTitle("Hit Efficiency vs pileup - "+GetLayerName(ilayer)); - if(name=="effVsCM") geff->SetTitle("Hit Efficiency vs common Mode - "+GetLayerName(ilayer)); - geff->SetMarkerStyle(20); - - } + for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) { + hfound = vhfound[ilayer]; + htotal = vhtotal[ilayer]; + + hfound->Sumw2(); + htotal->Sumw2(); + + // new ROOT version: TGraph::Divide don't handle null or negative values + for (Long_t i = 0; i < hfound->GetNbinsX() + 1; ++i) { + if (hfound->GetBinContent(i) == 0) + hfound->SetBinContent(i, 1e-6); + if (htotal->GetBinContent(i) == 0) + htotal->SetBinContent(i, 1); + } + + TGraphAsymmErrors* geff = fs->make(hfound->GetNbinsX()); + geff->SetName(Form("%s_layer%i", name.c_str(), ilayer)); + geff->BayesDivide(hfound, htotal); + if (name == "effVsLumi") + geff->SetTitle("Hit Efficiency vs inst. lumi. - " + GetLayerName(ilayer)); + if (name == "effVsPU") + geff->SetTitle("Hit Efficiency vs pileup - " + GetLayerName(ilayer)); + if (name == "effVsCM") + geff->SetTitle("Hit Efficiency vs common Mode - " + GetLayerName(ilayer)); + geff->SetMarkerStyle(20); + } } void SiStripHitResolFromCalibTree::makeSummaryVsLumi() { - cout<<"Computing efficiency vs lumi"<GetEntries()) // from infos per event - cout<<"Avg conditions (avg+/-rms): lumi :"<GetMean()<<"+/-"<GetRMS() - <<" pu: "<GetMean()<<"+/-"<GetRMS()<GetMean(); - layerPU=layertotal_vsPU[ilayer]->GetMean(); - //cout<<" layer "<GetEntries()) // from infos per event + cout << "Avg conditions (avg+/-rms): lumi :" << instLumiHisto->GetMean() << "+/-" << instLumiHisto->GetRMS() + << " pu: " << PUHisto->GetMean() << "+/-" << PUHisto->GetRMS() << endl; + + else { // from infos per hit + + unsigned int nLayers = 22; + if (_showRings) + nLayers = 20; + unsigned int nLayersForAvg = 0; + float layerLumi = 0; + float layerPU = 0; + float avgLumi = 0; + float avgPU = 0; + + cout << "Lumi summary: (avg over trajectory measurements)" << endl; + for (unsigned int ilayer = 1; ilayer < nLayers; ilayer++) { + layerLumi = layertotal_vsLumi[ilayer]->GetMean(); + layerPU = layertotal_vsPU[ilayer]->GetMean(); + //cout<<" layer "< 0 && k < 5) { + layername = TString("TIB L") + k; + } else if (k > 4 && k < 11) { + layername = TString("TOB L") + (k - 4); + } else if (k > 10 && k < 14) { + layername = TString("TID- ") + ringlabel + (k - 10); + } else if (k > 13 && k < 17) { + layername = TString("TID+ ") + ringlabel + (k - 13); + } else if (k > 16 && k < 17 + nTEClayers) { + layername = TString("TEC- ") + ringlabel + (k - 16); + } else if (k > 16 + nTEClayers) { + layername = TString("TEC+ ") + ringlabel + (k - 16 - nTEClayers); + } - TString layername=""; - TString ringlabel="D"; - if(_showRings) ringlabel="R"; - if (k>0 && k<5) { - layername = TString("TIB L") + k; - } else if (k>4&&k<11) { - layername = TString("TOB L")+(k-4); - } else if (k>10&&k<14) { - layername = TString("TID- ")+ringlabel+(k-10); - } else if (k>13&&k<17) { - layername = TString("TID+ ")+ringlabel+(k-13); - } else if (k>16&&k<17+nTEClayers) { - layername = TString("TEC- ")+ringlabel+(k-16); - } else if (k>16+nTEClayers) { - layername = TString("TEC+ ")+ringlabel+(k-16-nTEClayers); - } - - return layername; + return layername; } std::unique_ptr SiStripHitResolFromCalibTree::getNewObject() { //Need this for a Condition DB Writer //Initialize a return variable - auto obj= std::make_unique(); - - SiStripBadStrip::RegistryIterator rIter=quality_->getRegistryVectorBegin(); - SiStripBadStrip::RegistryIterator rIterEnd=quality_->getRegistryVectorEnd(); - - for(;rIter!=rIterEnd;++rIter){ - SiStripBadStrip::Range range(quality_->getDataVectorBegin()+rIter->ibegin,quality_->getDataVectorBegin()+rIter->iend); - if ( ! obj->put(rIter->detid,range) ) - edm::LogError("SiStripHitResolFromCalibTree")<<"[SiStripHitResolFromCalibTree::getNewObject] detid already exists"<(); + + SiStripBadStrip::RegistryIterator rIter = quality_->getRegistryVectorBegin(); + SiStripBadStrip::RegistryIterator rIterEnd = quality_->getRegistryVectorEnd(); + + for (; rIter != rIterEnd; ++rIter) { + SiStripBadStrip::Range range(quality_->getDataVectorBegin() + rIter->ibegin, + quality_->getDataVectorBegin() + rIter->iend); + if (!obj->put(rIter->detid, range)) + edm::LogError("SiStripHitResolFromCalibTree") + << "[SiStripHitResolFromCalibTree::getNewObject] detid already exists" << std::endl; } - - return std::move(obj); + + return obj; } float SiStripHitResolFromCalibTree::calcPhi(float x, float y) { float phi = 0; float Pi = 3.14159; - if((x>=0)&&(y>=0)) phi = atan(y/x); - else if((x>=0)&&(y<=0)) phi = atan(y/x) + 2*Pi; - else if((x<=0)&&(y>=0)) phi = atan(y/x) + Pi; - else phi = atan(y/x) + Pi; - phi = phi*180.0/Pi; + if ((x >= 0) && (y >= 0)) + phi = atan(y / x); + else if ((x >= 0) && (y <= 0)) + phi = atan(y / x) + 2 * Pi; + else if ((x <= 0) && (y >= 0)) + phi = atan(y / x) + Pi; + else + phi = atan(y / x) + Pi; + phi = phi * 180.0 / Pi; return phi; -} - -void SiStripHitResolFromCalibTree::SetBadComponents(int i, int component,SiStripQuality::BadComponent& BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4]){ - - int napv=reader->getNumberOfApvsAndStripLength(BC.detid).first; - - ssV[i][component] << "\n\t\t " - << BC.detid - << " \t " << BC.BadModule << " \t " - << ( (BC.BadFibers)&0x1 ) << " "; - if (napv==4) - ssV[i][component] << "x " <<( (BC.BadFibers>>1)&0x1 ); - - if (napv==6) - ssV[i][component] << ( (BC.BadFibers>>1)&0x1 ) << " " - << ( (BC.BadFibers>>2)&0x1 ); - ssV[i][component] << " \t " - << ( (BC.BadApvs)&0x1 ) << " " - << ( (BC.BadApvs>>1)&0x1 ) << " "; - if (napv==4) - ssV[i][component] << "x x " << ( (BC.BadApvs>>2)&0x1 ) << " " - << ( (BC.BadApvs>>3)&0x1 ); - if (napv==6) - ssV[i][component] << ( (BC.BadApvs>>2)&0x1 ) << " " - << ( (BC.BadApvs>>3)&0x1 ) << " " - << ( (BC.BadApvs>>4)&0x1 ) << " " - << ( (BC.BadApvs>>5)&0x1 ) << " "; - - if (BC.BadApvs){ - NBadComponent[i][0][2]+= ( (BC.BadApvs>>5)&0x1 )+ ( (BC.BadApvs>>4)&0x1 ) + ( (BC.BadApvs>>3)&0x1 ) + - ( (BC.BadApvs>>2)&0x1 )+ ( (BC.BadApvs>>1)&0x1 ) + ( (BC.BadApvs)&0x1 ); - NBadComponent[i][component][2]+= ( (BC.BadApvs>>5)&0x1 )+ ( (BC.BadApvs>>4)&0x1 ) + ( (BC.BadApvs>>3)&0x1 ) + - ( (BC.BadApvs>>2)&0x1 )+ ( (BC.BadApvs>>1)&0x1 ) + ( (BC.BadApvs)&0x1 ); +} + +void SiStripHitResolFromCalibTree::SetBadComponents( + int i, int component, SiStripQuality::BadComponent& BC, std::stringstream ssV[4][19], int NBadComponent[4][19][4]) { + int napv = _detInfo.getNumberOfApvsAndStripLength(BC.detid).first; + + ssV[i][component] << "\n\t\t " << BC.detid << " \t " << BC.BadModule << " \t " << ((BC.BadFibers) & 0x1) << " "; + if (napv == 4) + ssV[i][component] << "x " << ((BC.BadFibers >> 1) & 0x1); + + if (napv == 6) + ssV[i][component] << ((BC.BadFibers >> 1) & 0x1) << " " << ((BC.BadFibers >> 2) & 0x1); + ssV[i][component] << " \t " << ((BC.BadApvs) & 0x1) << " " << ((BC.BadApvs >> 1) & 0x1) << " "; + if (napv == 4) + ssV[i][component] << "x x " << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1); + if (napv == 6) + ssV[i][component] << ((BC.BadApvs >> 2) & 0x1) << " " << ((BC.BadApvs >> 3) & 0x1) << " " + << ((BC.BadApvs >> 4) & 0x1) << " " << ((BC.BadApvs >> 5) & 0x1) << " "; + + if (BC.BadApvs) { + NBadComponent[i][0][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) + ((BC.BadApvs >> 3) & 0x1) + + ((BC.BadApvs >> 2) & 0x1) + ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1); + NBadComponent[i][component][2] += ((BC.BadApvs >> 5) & 0x1) + ((BC.BadApvs >> 4) & 0x1) + + ((BC.BadApvs >> 3) & 0x1) + ((BC.BadApvs >> 2) & 0x1) + + ((BC.BadApvs >> 1) & 0x1) + ((BC.BadApvs) & 0x1); + } + if (BC.BadFibers) { + NBadComponent[i][0][1] += ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1); + NBadComponent[i][component][1] += + ((BC.BadFibers >> 2) & 0x1) + ((BC.BadFibers >> 1) & 0x1) + ((BC.BadFibers) & 0x1); } - if (BC.BadFibers){ - NBadComponent[i][0][1]+= ( (BC.BadFibers>>2)&0x1 )+ ( (BC.BadFibers>>1)&0x1 ) + ( (BC.BadFibers)&0x1 ); - NBadComponent[i][component][1]+= ( (BC.BadFibers>>2)&0x1 )+ ( (BC.BadFibers>>1)&0x1 ) + ( (BC.BadFibers)&0x1 ); - } - if (BC.BadModule){ + if (BC.BadModule) { NBadComponent[i][0][0]++; NBadComponent[i][component][0]++; } diff --git a/CalibTracker/SiStripHitResolution/python/SiStripHitResol_cff.py b/CalibTracker/SiStripHitResolution/python/SiStripHitResol_cff.py index 6c91d97723ff6..3b480e7406864 100644 --- a/CalibTracker/SiStripHitResolution/python/SiStripHitResol_cff.py +++ b/CalibTracker/SiStripHitResolution/python/SiStripHitResol_cff.py @@ -9,35 +9,24 @@ compressionSettings = 201 anResol = cms.EDAnalyzer("HitResol", - #CompressionSettings = cms.untracked.int32(compressionSettings), - Debug = cms.bool(False), - Layer = cms.int32(0), # =0 means do all layers - #combinatorialTracks = cms.InputTag("ctfWithMaterialTracksP5"), - #combinatorialTracks = cms.InputTag("TrackRefitterP5"), - #combinatorialTracks = cms.InputTag("ALCARECOTkAlCosmicsCTF0T"), - combinatorialTracks = cms.InputTag("generalTracks"), - #trajectories = cms.InputTag("ctfWithMaterialTracksP5"), - #trajectories = cms.InputTag("TrackRefitterP5"), - #trajectories = cms.InputTag("CalibrationTracksRefit"), - trajectories = cms.InputTag("generalTracks"), - siStripClusters = cms.InputTag("siStripClusters"), - siStripDigis = cms.InputTag("siStripDigis"), - trackerEvent = cms.InputTag("MeasurementTrackerEvent"), - lumiScalers = cms.InputTag("scalersRawToDigi"), - addLumi = cms.untracked.bool(False), - commonMode = cms.InputTag("siStripDigis", "CommonMode"), - addCommonMode = cms.untracked.bool(False), - # do not cut on the total number of tracks - cutOnTracks = cms.untracked.bool(True), - # compatibility - trackMultiplicity = cms.untracked.uint32(100), - # use or not first and last measurement of a trajectory (biases), default is false - useFirstMeas = cms.untracked.bool(False), - useLastMeas = cms.untracked.bool(False), - # use or not all hits when some missing hits in the trajectory (bias), default is false - useAllHitsFromTracksWithMissingHits = cms.untracked.bool(False), - MomentumCut = cms.untracked.double(3.), - UsePairsOnly = cms.untracked.uint32(1) - ) + CompressionSettings = cms.untracked.int32(compressionSettings), + Debug = cms.bool(False), + Layer = cms.int32(0), # = 0 means do all layers + #combinatorialTracks = cms.InputTag("ctfWithMaterialTracksP5"), + #combinatorialTracks = cms.InputTag("TrackRefitterP5"), + #combinatorialTracks = cms.InputTag("ALCARECOTkAlCosmicsCTF0T"), + combinatorialTracks = cms.InputTag("generalTracks"), + #trajectories = cms.InputTag("ctfWithMaterialTracksP5"), + #trajectories = cms.InputTag("TrackRefitterP5"), + #trajectories = cms.InputTag("CalibrationTracksRefit"), + trajectories = cms.InputTag("generalTracks"), + lumiScalers = cms.InputTag("scalersRawToDigi"), + addLumi = cms.untracked.bool(False), + # do not cut on the total number of tracks + cutOnTracks = cms.untracked.bool(True), + # compatibility + trackMultiplicity = cms.untracked.uint32(100), + MomentumCut = cms.untracked.double(3.), + UsePairsOnly = cms.untracked.uint32(1)) hitresol = cms.Sequence( anResol ) diff --git a/CalibTracker/SiStripHitResolution/src/HitResol.cc b/CalibTracker/SiStripHitResolution/src/HitResol.cc index 3e39107f7d555..2acb8a268a570 100644 --- a/CalibTracker/SiStripHitResolution/src/HitResol.cc +++ b/CalibTracker/SiStripHitResolution/src/HitResol.cc @@ -1,8 +1,9 @@ //////////////////////////////////////////////////////////////////////////////// // Package: CalibTracker/SiStripHitResolution // Class: HitResol -// Original Author: DG -// adapted from HitEff +// Original Authors: Denis Gele and Kathryn Coldham (adapted from HitEff) +// modified by Khawla Jaffel for CPE studies +// ported to cmssw by M. Musich // /////////////////////////////////////////////////////////////////////////////// @@ -11,386 +12,353 @@ #include #include -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "UserCode/SiStripHitResolution/interface/HitResol.h" -#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" +#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h" +#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" +#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h" +#include "CalibTracker/Records/interface/SiStripQualityRcd.h" +#include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h" +#include "CalibTracker/SiStripHitResolution/interface/HitResol.h" +#include "CommonTools/Statistics/interface/ChiSquaredProbability.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" #include "DataFormats/Common/interface/Handle.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" +#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "DataFormats/GeometryVector/interface/GlobalVector.h" #include "DataFormats/GeometryVector/interface/LocalVector.h" -#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" -#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" -#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDetType.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/SiStripDetId/interface/StripSubdetector.h" #include "DataFormats/TrackReco/interface/Track.h" -#include "DataFormats/TrackReco/interface/TrackFwd.h" -#include "DataFormats/TrackReco/interface/TrackExtra.h" #include "DataFormats/TrackReco/interface/TrackBase.h" -#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" -#include "TrackingTools/Records/interface/TransientRecHitRecord.h" -#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" -#include "UserCode/SiStripHitResolution/interface/TrajectoryAtInvalidHit.h" +#include "DataFormats/TrackReco/interface/TrackExtra.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" +#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" +#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonDetUnit/interface/GeomDetType.h" +#include "Geometry/CommonDetUnit/interface/GluedGeomDet.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" -#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h" -#include "DataFormats/TrackReco/interface/DeDxData.h" -#include "TrackingTools/DetLayers/interface/DetLayer.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" - -#include "RecoTracker/Record/interface/CkfComponentsRecord.h" -#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h" -#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h" -#include "CalibTracker/Records/interface/SiStripQualityRcd.h" -#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" -#include "Geometry/CommonDetUnit/interface/GluedGeomDet.h" -#include "DataFormats/Common/interface/DetSetVector.h" -#include "DataFormats/Common/interface/DetSetVectorNew.h" -#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" - -#include "DataFormats/MuonReco/interface/Muon.h" -#include "DataFormats/MuonReco/interface/MuonFwd.h" - +#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h" +#include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h" +#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h" +#include "TrackingTools/DetLayers/interface/DetLayer.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h" +#include "TrackingTools/Records/interface/TransientRecHitRecord.h" +#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" #include "TMath.h" #include "TH1F.h" -//ModifDG -#include "CommonTools/Statistics/interface/ChiSquaredProbability.h" -#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" -#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" -#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" -#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" -#include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h" -#include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h" -#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h" +#include "vdt/vdtMath.h" + // // constructors and destructor // - using namespace std; -HitResol::HitResol(const edm::ParameterSet& conf) : - scalerToken_( consumes< LumiScalersCollection >(conf.getParameter("lumiScalers")) ), - commonModeToken_( mayConsume< edm::DetSetVector >(conf.getParameter("commonMode")) ), - combinatorialTracks_token_( consumes< reco::TrackCollection >(conf.getParameter("combinatorialTracks")) ), - trajectories_token_( consumes< std::vector >(conf.getParameter("trajectories")) ), - trajTrackAsso_token_( consumes< TrajTrackAssociationCollection >(conf.getParameter("trajectories")) ), -//ModifDG - tjToken_( consumes< std::vector >( conf.getParameter("trajectories") )), -// tkToken_( consumes< reco::TrackCollection >( conf.getParameter("tracks") )), -//EndModifDG - clusters_token_( consumes< edmNew::DetSetVector >(conf.getParameter("siStripClusters")) ), - digis_token_( consumes< DetIdCollection >(conf.getParameter("siStripDigis")) ), -//ModifDG -// tracks_token_( consumes >( conf.getParameter("Tracks") )), -//EndModifDG - trackerEvent_token_( consumes< MeasurementTrackerEvent>(conf.getParameter("trackerEvent")) ), - conf_(conf) -{ - compSettings=conf_.getUntrackedParameter("CompressionSettings",-1); - layers =conf_.getParameter("Layer"); - DEBUG = conf_.getParameter("Debug"); - addLumi_ = conf_.getUntrackedParameter("addLumi", false); - addCommonMode_ = conf_.getUntrackedParameter("addCommonMode", false); - cutOnTracks_ = conf_.getUntrackedParameter("cutOnTracks", false); - trackMultiplicityCut_ = conf.getUntrackedParameter("trackMultiplicity",100); - useFirstMeas_ = conf_.getUntrackedParameter("useFirstMeas", false); - useLastMeas_ = conf_.getUntrackedParameter("useLastMeas", false); - useAllHitsFromTracksWithMissingHits_ = conf_.getUntrackedParameter("useAllHitsFromTracksWithMissingHits", false); - MomentumCut_ = conf_.getUntrackedParameter("MomentumCut", 3.); - UsePairsOnly_ = conf.getUntrackedParameter("UsePairsOnly",1); +HitResol::HitResol(const edm::ParameterSet& conf) + : scalerToken_(consumes(conf.getParameter("lumiScalers"))), + combinatorialTracks_token_( + consumes(conf.getParameter("combinatorialTracks"))), + tjToken_(consumes >(conf.getParameter("trajectories"))), + topoToken_(esConsumes()), + geomToken_(esConsumes()), + cpeToken_(esConsumes(edm::ESInputTag("", "StripCPEfromTrackAngle"))), + siStripQualityToken_(esConsumes()), + magFieldToken_(esConsumes()), + addLumi_(conf.getUntrackedParameter("addLumi", false)), + DEBUG_(conf.getParameter("Debug")), + cutOnTracks_(conf.getUntrackedParameter("cutOnTracks", false)), + momentumCut_(conf.getUntrackedParameter("MomentumCut", 3.)), + compSettings_(conf.getUntrackedParameter("CompressionSettings", -1)), + usePairsOnly_(conf.getUntrackedParameter("UsePairsOnly", 1)), + layers_(conf.getParameter("Layer")), + trackMultiplicityCut_(conf.getUntrackedParameter("trackMultiplicity", 100)) { + usesResource(TFileService::kSharedResource); } -// Virtual destructor needed. -HitResol::~HitResol() { } - -void HitResol::beginJob(){ - +void HitResol::beginJob() { edm::Service fs; - if(compSettings>0){ - edm::LogInfo("SiStripHitResolution:HitResol")<<"the compressions settings are:"<< compSettings << std::endl; - fs->file().SetCompressionSettings(compSettings); + if (compSettings_ > 0) { + edm::LogInfo("SiStripHitResolution:HitResol") << "the compressions settings are:" << compSettings_ << std::endl; + fs->file().SetCompressionSettings(compSettings_); } - reso = fs->make("reso","tree hit pairs for resolution studies"); - reso->Branch("momentum",&mymom,"momentum/F"); - reso->Branch("numHits",&numHits,"numHits/I"); - reso->Branch("trackChi2",&ProbTrackChi2,"trackChi2/F"); - reso->Branch("detID1",&iidd1,"detID1/I"); - reso->Branch("pitch1",&mypitch1,"pitch1/F"); - reso->Branch("clusterW1",&clusterWidth,"clusterW1/I"); - reso->Branch("expectedW1",&expWidth,"expectedW1/F"); - reso->Branch("atEdge1",&atEdge,"atEdge1/F"); - reso->Branch("simpleRes",&simpleRes,"simpleRes/F"); - reso->Branch("detID2",&iidd2,"detID2/I"); - reso->Branch("clusterW2",&clusterWidth_2,"clusterW2/I"); - reso->Branch("expectedW2",&expWidth_2,"expectedW2/F"); - reso->Branch("atEdge2",&atEdge_2,"atEdge2/F"); - reso->Branch("pairPath",&pairPath,"pairPath/F"); - reso->Branch("hitDX",&hitDX,"hitDX/F"); - reso->Branch("trackDX",&trackDX,"trackDX/F"); - reso->Branch("trackDXE",&trackDXE,"trackDXE/F"); - reso->Branch("trackParamX",&trackParamX,"trackParamX/F"); - reso->Branch("trackParamY",&trackParamY,"trackParamY/F"); - reso->Branch("trackParamDXDZ",&trackParamDXDZ,"trackParamDXDZ/F"); - reso->Branch("trackParamDYDZ",&trackParamDYDZ,"trackParamDYDZ/F"); - reso->Branch("trackParamXE",&trackParamXE,"trackParamXE/F"); - reso->Branch("trackParamYE",&trackParamYE,"trackParamYE/F"); - reso->Branch("trackParamDXDZE",&trackParamDXDZE,"trackParamDXDZE/F"); - reso->Branch("trackParamDYDZE",&trackParamDYDZE,"trackParamDYDZE/F"); - reso->Branch("pairsOnly",&pairsOnly,"pairsOnly/I"); - treso = fs->make("treso","tree tracks for resolution studies"); - treso->Branch("track_momentum",&track_momentum,"track_momentum/F"); - treso->Branch("track_eta",&track_eta,"track_eta/F"); - treso->Branch("track_trackChi2",&track_trackChi2,"track_trackChi2/F"); + reso = fs->make("reso", "tree hit pairs for resolution studies"); + reso->Branch("momentum", &mymom, "momentum/F"); + reso->Branch("numHits", &numHits, "numHits/I"); + reso->Branch("trackChi2", &ProbTrackChi2, "trackChi2/F"); + reso->Branch("detID1", &iidd1, "detID1/I"); + reso->Branch("pitch1", &mypitch1, "pitch1/F"); + reso->Branch("clusterW1", &clusterWidth, "clusterW1/I"); + reso->Branch("expectedW1", &expWidth, "expectedW1/F"); + reso->Branch("atEdge1", &atEdge, "atEdge1/F"); + reso->Branch("simpleRes", &simpleRes, "simpleRes/F"); + reso->Branch("detID2", &iidd2, "detID2/I"); + reso->Branch("clusterW2", &clusterWidth_2, "clusterW2/I"); + reso->Branch("expectedW2", &expWidth_2, "expectedW2/F"); + reso->Branch("atEdge2", &atEdge_2, "atEdge2/F"); + reso->Branch("pairPath", &pairPath, "pairPath/F"); + reso->Branch("hitDX", &hitDX, "hitDX/F"); + reso->Branch("trackDX", &trackDX, "trackDX/F"); + reso->Branch("trackDXE", &trackDXE, "trackDXE/F"); + reso->Branch("trackParamX", &trackParamX, "trackParamX/F"); + reso->Branch("trackParamY", &trackParamY, "trackParamY/F"); + reso->Branch("trackParamDXDZ", &trackParamDXDZ, "trackParamDXDZ/F"); + reso->Branch("trackParamDYDZ", &trackParamDYDZ, "trackParamDYDZ/F"); + reso->Branch("trackParamXE", &trackParamXE, "trackParamXE/F"); + reso->Branch("trackParamYE", &trackParamYE, "trackParamYE/F"); + reso->Branch("trackParamDXDZE", &trackParamDXDZE, "trackParamDXDZE/F"); + reso->Branch("trackParamDYDZE", &trackParamDYDZE, "trackParamDYDZE/F"); + reso->Branch("pairsOnly", &pairsOnly, "pairsOnly/I"); + treso = fs->make("treso", "tree tracks for resolution studies"); + treso->Branch("track_momentum", &track_momentum, "track_momentum/F"); + treso->Branch("track_pt", &track_pt, "track_pt/F"); + treso->Branch("track_eta", &track_eta, "track_eta/F"); + treso->Branch("track_phi", &track_phi, "track_phi/F"); + treso->Branch("track_trackChi2", &track_trackChi2, "track_trackChi2/F"); + treso->Branch("track_width", &expWidth, "track_width/F"); // from 1D HIT + treso->Branch("NumberOf_tracks", &NumberOf_tracks, "NumberOf_tracks/I"); events = 0; EventTrackCKF = 0; -} + histos2d_["track_phi_vs_eta"] = new TH2F("track_phi_vs_eta", ";track phi;track eta", 60, -3.5, 3.5, 60, -3., 3.); + histos2d_["residual_vs_trackMomentum"] = new TH2F("residual_vs_trackMomentum", + ";track momentum [GeV]; x_{pred_track} - x_{reco_hit} [#mum]", + 60, + 0., + 10., + 60, + 0., + 200.); + histos2d_["residual_vs_trackPt"] = new TH2F( + "residual_vs_trackPt", ";track p_{T}[GeV];x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 10., 60, 0., 200.); + histos2d_["residual_vs_trackEta"] = + new TH2F("residual_vs_trackEta", ";track #eta;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3., 60, 0., 200.); + histos2d_["residual_vs_trackPhi"] = + new TH2F("residual_vs_trackPhi", ";track #phi;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3.5, 60, 0., 200.); + histos2d_["residual_vs_expectedWidth"] = new TH2F( + "residual_vs_expectedWidth", ";track Width;x_{pred_track} - x_{reco_hit} [#mum]", 3, 0., 3., 60, 0., 200.); + histos2d_["numHits_vs_residual"] = + new TH2F("numHits_vs_residual", ";x_{pred_track} - x_{reco_hit} [#mum];N Hits", 60, 0., 200., 15, 0., 15.); +} -void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es){ +void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es) { //Retrieve tracker topology from geometry - edm::ESHandle tTopoHandle; - es.get().get(tTopoHandle); - const TrackerTopology* const tTopo = tTopoHandle.product(); - - // bool DEBUG = false; + const TrackerTopology* tTopo = &es.getData(topoToken_); - LogDebug("SiStripHitResolution:HitResol") << "beginning analyze from HitResol" << endl; + LogDebug("SiStripHitResolution:HitResol") << "beginning analyze from HitResol" << endl; using namespace edm; using namespace reco; + // Step A: Get Inputs int run_nr = e.id().run(); int ev_nr = e.id().event(); - // - edm::Handle > commonModeDigis; - if(addCommonMode_) e.getByToken(commonModeToken_, commonModeDigis); - - //CombinatoriaTrack + // get the tracks edm::Handle trackCollectionCKF; - e.getByToken(combinatorialTracks_token_,trackCollectionCKF); + e.getByToken(combinatorialTracks_token_, trackCollectionCKF); + const reco::TrackCollection* tracksCKF = trackCollectionCKF.product(); - edm::Handle > TrajectoryCollectionCKF; - e.getByToken(trajectories_token_,TrajectoryCollectionCKF); - - edm::Handle trajTrackAssociationHandle; - e.getByToken(trajTrackAsso_token_, trajTrackAssociationHandle); - - // Clusters - // get the SiStripClusters from the event - edm::Handle< edmNew::DetSetVector > theClusters; - e.getByToken(clusters_token_, theClusters); + // get the trajectory collection + edm::Handle > trajectoryCollectionHandle; + e.getByToken(tjToken_, trajectoryCollectionHandle); + const TrajectoryCollection* trajectoryCollection = trajectoryCollectionHandle.product(); //get tracker geometry - edm::ESHandle tracker; - es.get().get(tracker); - const TrackerGeometry * tkgeom=&(* tracker); + edm::ESHandle tracker = es.getHandle(geomToken_); + const TrackerGeometry* tkgeom = &(*tracker); //get Cluster Parameter Estimator - edm::ESHandle parameterestimator; - es.get().get("StripCPEfromTrackAngle", parameterestimator); - const StripClusterParameterEstimator &stripcpe(*parameterestimator); + edm::ESHandle parameterestimator = es.getHandle(cpeToken_); + const StripClusterParameterEstimator& stripcpe(*parameterestimator); // get the SiStripQuality records - edm::ESHandle SiStripQuality_; - es.get().get(SiStripQuality_); + edm::ESHandle SiStripQuality_ = es.getHandle(siStripQualityToken_); - edm::ESHandle magFieldHandle; - es.get().get(magFieldHandle); - const MagneticField* magField_ = magFieldHandle.product(); + // get the magnetic field + const MagneticField* magField_ = &es.getData(magFieldToken_); - // get the list of module IDs with FED-detected errors - edm::Handle< DetIdCollection > fedErrorIds; - e.getByToken(digis_token_, fedErrorIds ); - - ESHandle measurementTrackerHandle; - es.get().get(measurementTrackerHandle); - - edm::Handle measurementTrackerEvent; - e.getByToken(trackerEvent_token_,measurementTrackerEvent); - - edm::ESHandle est; - es.get().get("Chi2",est); - -// edm::ESHandle prop; -// es.get().get("PropagatorWithMaterial",prop); -// const Propagator* thePropagator = prop.product(); + events++; - edm::Handle > trajectoryCollectionHandle; - e.getByToken(tjToken_,trajectoryCollectionHandle); - const TrajectoryCollection* trajectoryCollection = trajectoryCollectionHandle.product(); + // List of variables for SiStripHitResolution ntuple + mymom = 0; + numHits = 0; + ProbTrackChi2 = 0; + iidd1 = 0; + mypitch1 = 0; + clusterWidth = 0; + expWidth = 0; + atEdge = 0; + simpleRes = 0; + iidd2 = 0; + clusterWidth_2 = 0; + expWidth_2 = 0; + atEdge_2 = 0; + pairPath = 0; + hitDX = 0; + trackDX = 0; + trackDXE = 0; + trackParamX = 0; + trackParamY = 0; + trackParamDXDZ = 0; + trackParamDYDZ = 0; + trackParamXE = 0; + trackParamYE = 0; + trackParamDXDZE = 0; + trackParamDYDZE = 0; + pairsOnly = 0; - events++; + LogDebug("HitResol") << "Starting analysis, nrun nevent, tracksCKF->size(): " << run_nr << " " << ev_nr << " " + << tracksCKF->size() << std::endl; -// // *************** SiStripCluster Collection -// const edmNew::DetSetVector& input = *theClusters; - - -// List of variables for SiStripHitResolution ntuple - mymom = 0; - numHits = 0; - ProbTrackChi2 = 0; - iidd1 = 0; - mypitch1 = 0; - clusterWidth = 0; - expWidth = 0; - atEdge = 0; - simpleRes = 0; - iidd2 = 0; - clusterWidth_2 = 0; - expWidth_2 = 0; - atEdge_2 = 0; - pairPath = 0; - hitDX = 0; - trackDX = 0; - trackDXE = 0; - trackParamX = 0; - trackParamY = 0; - trackParamDXDZ = 0; - trackParamDYDZ = 0; - trackParamXE = 0; - trackParamYE = 0; - trackParamDXDZE = 0; - trackParamDYDZE = 0; - pairsOnly = 0; - - // Tracking - const reco::TrackCollection *tracksCKF=trackCollectionCKF.product(); - -////// Plugin of Nico code: - - std::cout<<"Starting analysis, nrun nevent, tracksCKF->size(): "<size() <size(); ++iT){ - track_momentum = tracksCKF->at(iT).pt(); + for (unsigned int iT = 0; iT < tracksCKF->size(); ++iT) { + track_momentum = tracksCKF->at(iT).p(); + track_pt = tracksCKF->at(iT).p(); track_eta = tracksCKF->at(iT).eta(); - track_trackChi2 = ChiSquaredProbability((double)( tracksCKF->at(iT).chi2() ),(double)( tracksCKF->at(iT).ndof() )); + track_phi = tracksCKF->at(iT).phi(); + track_trackChi2 = ChiSquaredProbability((double)(tracksCKF->at(iT).chi2()), (double)(tracksCKF->at(iT).ndof())); treso->Fill(); } + histos2d_["track_phi_vs_eta"]->Fill(track_phi, track_eta); + // loop over trajectories from refit - for ( const auto& traj : *trajectoryCollection ) { + for (const auto& traj : *trajectoryCollection) { const auto& TMeas = traj.measurements(); // Loop on each measurement and take it into consideration //-------------------------------------------------------- - for ( auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm ) { + for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) { if (!itm->updatedState().isValid()) { - std::cout<<"NONVALIDE"<recHit(); - const TrackingRecHit* myhit = (*itm->recHit()).hit(); + const TrackingRecHit* myhit = (*itm->recHit()).hit(); ProbTrackChi2 = 0; numHits = 0; -//// std::cout<<"TrackChi2 = "<< ChiSquaredProbability((double)( itm->chiSquared() ),(double)( itm->ndof(false) )) <updatedState().globalMomentum().perp(): "<< itm->updatedState().globalMomentum().perp() <foundHits() <updatedState().globalMomentum().perp(): " + << itm->updatedState().globalMomentum().perp() << "\n" + << "numhits " << traj.foundHits() << std::endl; - mymom = itm->updatedState().globalMomentum().perp(); + numHits = traj.foundHits(); + ProbTrackChi2 = ChiSquaredProbability((double)(traj.chiSquared()), (double)(traj.ndof(false))); - std::cout<<"mymom "<updatedState().globalMomentum().perp(); //Now for the first hit TrajectoryStateOnSurface mytsos = itm->updatedState(); const auto hit1 = itm->recHit(); DetId id1 = hit1->geographicalId(); -// if(id1.subdetId() < StripSubdetector::TIB || id1.subdetId() > StripSubdetector::TEC) continue; -// if(id1.subdetId() < StripSubdetector::TIB || id1.subdetId() > StripSubdetector::TEC) std::cout<<"AUTRE"<isValid() && mymom > MomentumCut_) std::cout<<"mymom: "<< itm->updatedState().globalMomentum().perp() <isValid() && - mymom > MomentumCut_ && - ( id1.subdetId() >= StripSubdetector::TIB && id1.subdetId() <= StripSubdetector::TEC) ) { - - const auto stripdet = dynamic_cast(tkgeom->idToDetUnit( hit1->geographicalId())); - const StripTopology& Topo = stripdet->specificTopology(); - int Nstrips = Topo.nstrips(); - mypitch1 = stripdet->surface().bounds().width() / Topo.nstrips(); - - const auto det = dynamic_cast(tkgeom->idToDetUnit(mypointhit->geographicalId())); - - TrajectoryStateOnSurface mytsos = itm->updatedState(); - LocalVector trackDirection = mytsos.localDirection(); - LocalVector drift = stripcpe.driftDirection(stripdet); - - const auto hit1d = dynamic_cast(myhit); - - if (hit1d) { -// float myres = getSimHitRes(det,trackDirection,*hit1d,expWidth,&mypitch1,drift); - getSimHitRes(det,trackDirection,*hit1d,expWidth,&mypitch1,drift); - clusterWidth = hit1d->cluster()->amplitudes().size(); - uint16_t firstStrip = hit1d->cluster()->firstStrip(); - uint16_t lastStrip = firstStrip + (hit1d->cluster()->amplitudes()).size() -1; - atEdge = (firstStrip == 0 || lastStrip == (Nstrips-1) ); - } - - - const auto hit2d = dynamic_cast(myhit); - - if (hit2d) { -// float myres = getSimHitRes(det,trackDirection,*hit2d, expWidth,&mypitch1,drift); - getSimHitRes(det,trackDirection,*hit2d, expWidth,&mypitch1,drift); - clusterWidth = hit2d->cluster()->amplitudes().size(); - uint16_t firstStrip = hit2d->cluster()->firstStrip(); - uint16_t lastStrip = firstStrip + (hit2d->cluster()->amplitudes()).size() -1; - atEdge = (firstStrip == 0 || lastStrip == (Nstrips-1) ); - } - - simpleRes = getSimpleRes(&(*itm)); // simple resolution by using the track re-fit forward and backward predicted state - - // Now to see if there is a match - pair method - hit in overlapping sensors - vector < TrajectoryMeasurement >::const_iterator itTraj2 = TMeas.end(); // last hit along the fitted track - - for ( auto itmCompare = itm-1; - // start to compare from the 5th hit - itmCompare >= TMeas.cbegin() && itmCompare > itm-4; - --itmCompare ) { - const auto hit2 = itmCompare->recHit(); - if (!hit2->isValid()) continue; - DetId id2 = hit2->geographicalId(); - - //must be from the same detector and layer - iidd1 = hit1->geographicalId().rawId(); - iidd2 = hit2->geographicalId().rawId(); - if ( id1.subdetId() != id2.subdetId() || checkLayer(iidd1, tTopo) != checkLayer(iidd2, tTopo) ) break; - //must both be stereo if one is - if(tTopo->isStereo(id1) != tTopo->isStereo(id2) ) continue; - //A check i dont completely understand but might as well keep there - if (tTopo->glued(id1) == id1.rawId() ) cout << "BAD GLUED: Have glued layer with id = " << id1.rawId() << " and glued id = " << tTopo->glued(id1) << " and stereo = " << tTopo->isStereo(id1) << endl; - if (tTopo->glued(id2) == id2.rawId() ) cout << "BAD GLUED: Have glued layer with id = " << id2.rawId() << " and glued id = " << tTopo->glued(id2) << " and stereo = " << tTopo->isStereo(id2) << endl; - - itTraj2 = itmCompare; - break; - } - - if ( itTraj2 == TMeas.cend() ) { - } else { -// std::cout<<"Found overlapping sensors "< StripSubdetector::TEC) + continue; -// pairsOnly = 1; - pairsOnly = UsePairsOnly_; + if (hit1->isValid() && mymom > momentumCut_ && + (id1.subdetId() >= StripSubdetector::TIB && id1.subdetId() <= StripSubdetector::TEC)) { + const auto stripdet = dynamic_cast(tkgeom->idToDetUnit(hit1->geographicalId())); + const StripTopology& Topo = stripdet->specificTopology(); + int Nstrips = Topo.nstrips(); + mypitch1 = stripdet->surface().bounds().width() / Topo.nstrips(); + + const auto det = dynamic_cast(tkgeom->idToDetUnit(mypointhit->geographicalId())); + + TrajectoryStateOnSurface mytsos = itm->updatedState(); + LocalVector trackDirection = mytsos.localDirection(); + LocalVector drift = stripcpe.driftDirection(stripdet); + + const auto hit1d = dynamic_cast(myhit); + + if (hit1d) { + getSimHitRes(det, trackDirection, *hit1d, expWidth, &mypitch1, drift); + clusterWidth = hit1d->cluster()->amplitudes().size(); + uint16_t firstStrip = hit1d->cluster()->firstStrip(); + uint16_t lastStrip = firstStrip + (hit1d->cluster()->amplitudes()).size() - 1; + atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1)); + } + + const auto hit2d = dynamic_cast(myhit); + + if (hit2d) { + getSimHitRes(det, trackDirection, *hit2d, expWidth, &mypitch1, drift); + clusterWidth = hit2d->cluster()->amplitudes().size(); + uint16_t firstStrip = hit2d->cluster()->firstStrip(); + uint16_t lastStrip = firstStrip + (hit2d->cluster()->amplitudes()).size() - 1; + atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1)); + } + + simpleRes = + getSimpleRes(&(*itm)); // simple resolution by using the track re-fit forward and backward predicted state + + histos2d_["residual_vs_trackMomentum"]->Fill(itm->updatedState().globalMomentum().mag(), + simpleRes * 10000); // reso in cm *10000 == micro-meter + histos2d_["residual_vs_trackPt"]->Fill(mymom, simpleRes * 10000); // reso in cm *10000 == micro-meter + histos2d_["residual_vs_trackEta"]->Fill(itm->updatedState().globalMomentum().eta(), simpleRes * 10000); + histos2d_["residual_vs_trackPhi"]->Fill(itm->updatedState().globalMomentum().phi(), simpleRes * 10000); + histos2d_["residual_vs_expectedWidth"]->Fill(expWidth, simpleRes * 10000); + histos2d_["numHits_vs_residual"]->Fill(simpleRes * 10000, numHits); + + // Now to see if there is a match - pair method - hit in overlapping sensors + vector::const_iterator itTraj2 = TMeas.end(); // last hit along the fitted track + + for (auto itmCompare = itm - 1; + // start to compare from the 5th hit + itmCompare >= TMeas.cbegin() && itmCompare > itm - 4; + --itmCompare) { + const auto hit2 = itmCompare->recHit(); + if (!hit2->isValid()) + continue; + DetId id2 = hit2->geographicalId(); + + //must be from the same detector and layer + iidd1 = hit1->geographicalId().rawId(); + iidd2 = hit2->geographicalId().rawId(); + if (id1.subdetId() != id2.subdetId() || checkLayer(iidd1, tTopo) != checkLayer(iidd2, tTopo)) + break; + //must both be stereo if one is + if (tTopo->isStereo(id1) != tTopo->isStereo(id2)) + continue; + //A check i dont completely understand but might as well keep there + if (tTopo->glued(id1) == id1.rawId()) + LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id1.rawId() + << " and glued id = " << tTopo->glued(id1) << " and stereo = " << tTopo->isStereo(id1) + << endl; + if (tTopo->glued(id2) == id2.rawId()) + LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id2.rawId() + << " and glued id = " << tTopo->glued(id2) << " and stereo = " << tTopo->isStereo(id2) + << endl; + + itTraj2 = itmCompare; + break; + } + + if (itTraj2 == TMeas.cend()) { + } else { + LogDebug("HitResol") << "Found overlapping sensors " << std::endl; + pairsOnly = usePairsOnly_; //We found one....let's fill in the truth info! TrajectoryStateOnSurface tsos_2 = itTraj2->updatedState(); @@ -398,12 +366,12 @@ void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es){ const auto myhit2 = itTraj2->recHit(); const auto myhit_2 = (*itTraj2->recHit()).hit(); const auto stripdet_2 = dynamic_cast(tkgeom->idToDetUnit(myhit2->geographicalId())); - const StripTopology& Topo_2 = stripdet_2->specificTopology(); + const StripTopology& Topo_2 = stripdet_2->specificTopology(); int Nstrips_2 = Topo_2.nstrips(); float mypitch_2 = stripdet_2->surface().bounds().width() / Topo_2.nstrips(); - if ( mypitch1 != mypitch_2 ) return; // for PairsOnly - + if (mypitch1 != mypitch_2) + return; // for PairsOnly const auto det_2 = dynamic_cast(tkgeom->idToDetUnit(myhit2->geographicalId())); @@ -411,179 +379,158 @@ void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es){ const auto hit1d_2 = dynamic_cast(myhit_2); if (hit1d_2) { - getSimHitRes(det_2,trackDirection_2,*hit1d_2,expWidth_2,&mypitch_2,drift_2); - clusterWidth_2 = hit1d_2->cluster()->amplitudes().size(); - uint16_t firstStrip_2 = hit1d_2->cluster()->firstStrip(); - uint16_t lastStrip_2 = firstStrip_2 + (hit1d_2->cluster()->amplitudes()).size() -1; - atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2-1) ); + getSimHitRes(det_2, trackDirection_2, *hit1d_2, expWidth_2, &mypitch_2, drift_2); + clusterWidth_2 = hit1d_2->cluster()->amplitudes().size(); + uint16_t firstStrip_2 = hit1d_2->cluster()->firstStrip(); + uint16_t lastStrip_2 = firstStrip_2 + (hit1d_2->cluster()->amplitudes()).size() - 1; + atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1)); } const auto hit2d_2 = dynamic_cast(myhit_2); if (hit2d_2) { -// float myres_2 = getSimHitRes(det_2,trackDirection_2,*hit2d_2, expWidth_2,&mypitch_2,drift_2); - getSimHitRes(det_2,trackDirection_2,*hit2d_2, expWidth_2,&mypitch_2,drift_2); - clusterWidth_2 = hit2d_2->cluster()->amplitudes().size(); - uint16_t firstStrip_2 = hit2d_2->cluster()->firstStrip(); - uint16_t lastStrip_2 = firstStrip_2 + (hit2d_2->cluster()->amplitudes()).size() -1; - atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2-1) ); + getSimHitRes(det_2, trackDirection_2, *hit2d_2, expWidth_2, &mypitch_2, drift_2); + clusterWidth_2 = hit2d_2->cluster()->amplitudes().size(); + uint16_t firstStrip_2 = hit2d_2->cluster()->firstStrip(); + uint16_t lastStrip_2 = firstStrip_2 + (hit2d_2->cluster()->amplitudes()).size() - 1; + atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1)); } - -// if(pairsOnly && (pitch != pitch2) ) fill = false; - - - // Make AnalyticalPropagator to use in getPairParameters - AnalyticalPropagator mypropagator(magField_,anyDirection); - - if(!getPairParameters(&(*magField_),mypropagator,&(*itTraj2),&(*itm),pairPath,hitDX,trackDX,trackDXE,trackParamX,trackParamY,trackParamDXDZ,trackParamDYDZ,trackParamXE,trackParamYE,trackParamDXDZE,trackParamDYDZE) ){ - - - }else { - -// std::cout<<" "<Fill(); - } - }//itTraj2 != TMeas.end() - - - }//hit1->isValid().... - - } // itm - } // it - - - - - - - - - -////// Endof Plugin of Nico code: - - + // if(pairsOnly && (pitch != pitch2) ) fill = false; + + // Make AnalyticalPropagator to use in getPairParameters + AnalyticalPropagator mypropagator(magField_, anyDirection); + + if (!getPairParameters(&(*magField_), + mypropagator, + &(*itTraj2), + &(*itm), + pairPath, + hitDX, + trackDX, + trackDXE, + trackParamX, + trackParamY, + trackParamDXDZ, + trackParamDYDZ, + trackParamXE, + trackParamYE, + trackParamDXDZE, + trackParamDYDZE)) { + } else { + LogDebug("HitResol") << " \n\n\n" + << " momentum " << mymom << "\n" + << " numHits " << numHits << "\n" + << " trackChi2 " << ProbTrackChi2 << "\n" + << " detID1 " << iidd1 << "\n" + << " pitch1 " << mypitch1 << "\n" + << " clusterW1 " << clusterWidth << "\n" + << " expectedW1 " << expWidth << "\n" + << " atEdge1 " << atEdge << "\n" + << " simpleRes " << simpleRes << "\n" + << " detID2 " << iidd2 << "\n" + << " clusterW2 " << clusterWidth_2 << "\n" + << " expectedW2 " << expWidth_2 << "\n" + << " atEdge2 " << atEdge_2 << "\n" + << " pairPath " << pairPath << "\n" + << " hitDX " << hitDX << "\n" + << " trackDX " << trackDX << "\n" + << " trackDXE " << trackDXE << "\n" + << " trackParamX " << trackParamX << "\n" + << " trackParamY " << trackParamY << "\n" + << " trackParamDXDZ " << trackParamDXDZ << "\n" + << " trackParamDYDZ " << trackParamDYDZ << "\n" + << " trackParamXE " << trackParamXE << "\n" + << " trackParamYE " << trackParamYE << "\n" + << " trackParamDXDZE" << trackParamDXDZE << "\n" + << " trackParamDYDZE" << trackParamDYDZE << std::endl; + reso->Fill(); + } + } //itTraj2 != TMeas.end() + } //hit1->isValid().... + } // itm + } // it } - - - - - -void HitResol::endJob(){ -// traj->GetDirectory()->cd(); -// traj->Write(); - - LogDebug("SiStripHitResolution:HitResol") << " Events Analysed " << events << endl; - LogDebug("SiStripHitResolution:HitResol") << " Number Of Tracked events " << EventTrackCKF << endl; +void HitResol::endJob() { + LogDebug("SiStripHitResolution:HitResol") << " Events Analysed " << events << endl; + LogDebug("SiStripHitResolution:HitResol") << " Number Of Tracked events " << EventTrackCKF << endl; reso->GetDirectory()->cd(); reso->Write(); treso->Write(); } -double HitResol::checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters, double xx, double xerr) { - double error = sqrt(parameters.second.xx() + xerr*xerr); +double HitResol::checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters, + double xx, + double xerr) { + double error = sqrt(parameters.second.xx() + xerr * xerr); double separation = abs(parameters.first.x() - xx); - double consistency = separation/error; + double consistency = separation / error; return consistency; } bool HitResol::isDoubleSided(unsigned int iidd, const TrackerTopology* tTopo) const { - StripSubdetector strip=StripSubdetector(iidd); - unsigned int subid=strip.subdetId(); + StripSubdetector strip = StripSubdetector(iidd); + unsigned int subid = strip.subdetId(); unsigned int layer = 0; - if (subid == StripSubdetector::TIB) { + if (subid == StripSubdetector::TIB) { layer = tTopo->tibLayer(iidd); - if (layer == 1 || layer == 2) return true; - else return false; - } - else if (subid == StripSubdetector::TOB) { - layer = tTopo->tobLayer(iidd) + 4 ; - if (layer == 5 || layer == 6) return true; - else return false; - } - else if (subid == StripSubdetector::TID) { + if (layer == 1 || layer == 2) + return true; + else + return false; + } else if (subid == StripSubdetector::TOB) { + layer = tTopo->tobLayer(iidd) + 4; + if (layer == 5 || layer == 6) + return true; + else + return false; + } else if (subid == StripSubdetector::TID) { layer = tTopo->tidRing(iidd) + 10; - if (layer == 11 || layer == 12) return true; - else return false; - } - else if (subid == StripSubdetector::TEC) { - layer = tTopo->tecRing(iidd) + 13 ; - if (layer == 14 || layer == 15 || layer == 18) return true; - else return false; - } - else + if (layer == 11 || layer == 12) + return true; + else + return false; + } else if (subid == StripSubdetector::TEC) { + layer = tTopo->tecRing(iidd) + 13; + if (layer == 14 || layer == 15 || layer == 18) + return true; + else + return false; + } else return false; } -//float HitResol::getSimHitRes(const GeomDetUnit * det, const LocalVector& trackdirection, const TrackingRecHit& recHit, float& trackWidth, float* pitch, LocalVector& drift){ -void HitResol::getSimHitRes(const GeomDetUnit * det, const LocalVector& trackdirection, const TrackingRecHit& recHit, float& trackWidth, float* pitch, LocalVector& drift){ +void HitResol::getSimHitRes(const GeomDetUnit* det, + const LocalVector& trackdirection, + const TrackingRecHit& recHit, + float& trackWidth, + float* pitch, + LocalVector& drift) { const auto stripdet = dynamic_cast(det); const auto& topol = dynamic_cast(stripdet->topology()); LocalPoint position = recHit.localPosition(); (*pitch) = topol.localPitch(position); -// float rechitrphiresMF = -1; - float anglealpha = 0; if (trackdirection.z() != 0) { anglealpha = atan(trackdirection.x() / trackdirection.z()) * TMath::RadToDeg(); } - -// LocalVector drift = stripcpe.driftDirection(stripdet); + // LocalVector drift = stripcpe.driftDirection(stripdet); float thickness = stripdet->surface().bounds().thickness(); float tanalpha = tan(anglealpha * TMath::DegToRad()); float tanalphaL = drift.x() / drift.z(); (trackWidth) = fabs((thickness / (*pitch)) * tanalpha - (thickness / (*pitch)) * tanalphaL); - - -// return rechitrphiresMF; - } -double HitResol::getSimpleRes(const TrajectoryMeasurement* traj1){ +double HitResol::getSimpleRes(const TrajectoryMeasurement* traj1) { TrajectoryStateOnSurface theCombinedPredictedState; - if ( traj1->backwardPredictedState().isValid() ) - theCombinedPredictedState = TrajectoryStateCombiner().combine( traj1->forwardPredictedState(), - traj1->backwardPredictedState()); + if (traj1->backwardPredictedState().isValid()) + theCombinedPredictedState = + TrajectoryStateCombiner().combine(traj1->forwardPredictedState(), traj1->backwardPredictedState()); else theCombinedPredictedState = traj1->forwardPredictedState(); @@ -591,55 +538,74 @@ double HitResol::getSimpleRes(const TrajectoryMeasurement* traj1){ return -100; } - const TransientTrackingRecHit::ConstRecHitPointer firstRecHit = traj1->recHit(); + const TransientTrackingRecHit::ConstRecHitPointer& firstRecHit = traj1->recHit(); double recHitX_1 = firstRecHit->localPosition().x(); return (theCombinedPredictedState.localPosition().x() - recHitX_1); } - //traj1 is the matched trajectory...traj2 is the original -bool HitResol::getPairParameters(const MagneticField* magField_, AnalyticalPropagator& propagator, const TrajectoryMeasurement* traj1, const TrajectoryMeasurement* traj2, float & pairPath,float & hitDX, float & trackDX, float & trackDXE, float & trackParamX, float &trackParamY , float & trackParamDXDZ, float &trackParamDYDZ , float & trackParamXE, float &trackParamYE, float & trackParamDXDZE, float &trackParamDYDZE){ - pairPath = 0; - hitDX = 0; - trackDX = 0; - trackDXE = 0; - - trackParamX = 0; - trackParamY = 0; - trackParamDXDZ = 0; - trackParamDYDZ = 0; - trackParamXE = 0; - trackParamYE = 0; +bool HitResol::getPairParameters(const MagneticField* magField_, + AnalyticalPropagator& propagator, + const TrajectoryMeasurement* traj1, + const TrajectoryMeasurement* traj2, + float& pairPath, + float& hitDX, + float& trackDX, + float& trackDXE, + float& trackParamX, + float& trackParamY, + float& trackParamDXDZ, + float& trackParamDYDZ, + float& trackParamXE, + float& trackParamYE, + float& trackParamDXDZE, + float& trackParamDYDZE) { + pairPath = 0; + hitDX = 0; + trackDX = 0; + trackDXE = 0; + + trackParamX = 0; + trackParamY = 0; + trackParamDXDZ = 0; + trackParamDYDZ = 0; + trackParamXE = 0; + trackParamYE = 0; trackParamDXDZE = 0; trackParamDYDZE = 0; TrajectoryStateCombiner combiner_; // backward predicted state at module 1 - TrajectoryStateOnSurface bwdPred1 = traj1->backwardPredictedState(); - if ( !bwdPred1.isValid() ) return false; - //cout << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl; + const TrajectoryStateOnSurface& bwdPred1 = traj1->backwardPredictedState(); + if (!bwdPred1.isValid()) + return false; + LogDebug("HitResol") << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl; // forward predicted state at module 2 - TrajectoryStateOnSurface fwdPred2 = traj2->forwardPredictedState(); - //cout << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl; - if ( !fwdPred2.isValid() ) return false; + const TrajectoryStateOnSurface& fwdPred2 = traj2->forwardPredictedState(); + LogDebug("HitResol") << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl; + if (!fwdPred2.isValid()) + return false; // extrapolate fwdPred2 to module 1 - TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2,bwdPred1.surface()); - if ( !fwdPred2At1.isValid() ) return false; + TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2, bwdPred1.surface()); + if (!fwdPred2At1.isValid()) + return false; // combine fwdPred2At1 with bwdPred1 (ref. state, best estimate without hits 1 and 2) - TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1,fwdPred2At1); - if ( !comb1.isValid() ) return false; + TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1, fwdPred2At1); + if (!comb1.isValid()) + return false; // // propagation of reference parameters to module 2 // - std::pair tsosWithS = - propagator.propagateWithPath(comb1,fwdPred2.surface()); + std::pair tsosWithS = propagator.propagateWithPath(comb1, fwdPred2.surface()); TrajectoryStateOnSurface comb1At2 = tsosWithS.first; - if ( !comb1At2.isValid() ) return false; + if (!comb1At2.isValid()) + return false; //distance of propagation from one surface to the next==could cut here pairPath = tsosWithS.second; - if (TMath::Abs(pairPath) > 15 ) return false; //cut to remove hit pairs > 15 cm apart + if (TMath::Abs(pairPath) > 15) + return false; //cut to remove hit pairs > 15 cm apart // local parameters and errors on module 1 AlgebraicVector5 pars = comb1.localParameters().vector(); @@ -647,14 +613,14 @@ bool HitResol::getPairParameters(const MagneticField* magField_, AnalyticalPropa //number 3 is predX double predX1 = pars[3]; //track fitted parameters in local coordinates for position 0 - (trackParamX ) = pars[3]; - (trackParamY ) = pars[4]; - (trackParamDXDZ ) = pars[1]; - (trackParamDYDZ ) = pars[2]; - (trackParamXE ) = TMath::Sqrt(errs(3,3)); - (trackParamYE ) = TMath::Sqrt(errs(4,4)); - (trackParamDXDZE) = TMath::Sqrt(errs(1,1)); - (trackParamDYDZE) = TMath::Sqrt(errs(2,2)); + (trackParamX) = pars[3]; + (trackParamY) = pars[4]; + (trackParamDXDZ) = pars[1]; + (trackParamDYDZ) = pars[2]; + (trackParamXE) = TMath::Sqrt(errs(3, 3)); + (trackParamYE) = TMath::Sqrt(errs(4, 4)); + (trackParamDXDZE) = TMath::Sqrt(errs(1, 1)); + (trackParamDYDZE) = TMath::Sqrt(errs(2, 2)); // local parameters and errors on module 2 pars = comb1At2.localParameters().vector(); @@ -664,41 +630,35 @@ bool HitResol::getPairParameters(const MagneticField* magField_, AnalyticalPropa //// //// jacobians (local-to-global@1,global 1-2,global-to-local@2) //// - JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(), - comb1.localParameters(), - *magField_); - AnalyticalCurvilinearJacobian jacCurvToCurv(comb1.globalParameters(), - comb1At2.globalPosition(), - comb1At2.globalMomentum(), - tsosWithS.second); - JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(), - comb1At2.localParameters(), - *magField_); + JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(), comb1.localParameters(), *magField_); + AnalyticalCurvilinearJacobian jacCurvToCurv( + comb1.globalParameters(), comb1At2.globalPosition(), comb1At2.globalMomentum(), tsosWithS.second); + JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(), comb1At2.localParameters(), *magField_); // combined jacobian local-1-to-local-2 - AlgebraicMatrix55 jacobian = - jacLocToCurv.jacobian()*jacCurvToCurv.jacobian()*jacCurvToLoc.jacobian(); + AlgebraicMatrix55 jacobian = jacLocToCurv.jacobian() * jacCurvToCurv.jacobian() * jacCurvToLoc.jacobian(); // covariance on module 1 AlgebraicSymMatrix55 covComb1 = comb1.localError().matrix(); // variance and correlations for predicted local_x on modules 1 and 2 - double c00 = covComb1(3,3); + double c00 = covComb1(3, 3); double c10(0.); double c11(0.); - for ( int i=1; i<5; ++i ) { - c10 += jacobian(3,i)*covComb1(i,3); - for ( int j=1; j<5; ++j ) c11 += jacobian(3,i)*covComb1(i,j)*jacobian(3,j); + for (int i = 1; i < 5; ++i) { + c10 += jacobian(3, i) * covComb1(i, 3); + for (int j = 1; j < 5; ++j) + c11 += jacobian(3, i) * covComb1(i, j) * jacobian(3, j); } // choose relative sign in order to minimize error on difference - double diff = c00 - 2*fabs(c10) + c11; - diff = diff>0 ? sqrt(diff) : -sqrt(-diff); + double diff = c00 - 2 * fabs(c10) + c11; + diff = diff > 0 ? sqrt(diff) : -sqrt(-diff); (trackDXE) = diff; - double relativeXSign_ = c10>0 ? -1 : 1; + double relativeXSign_ = c10 > 0 ? -1 : 1; - (trackDX) = predX1 + relativeXSign_*predX2; + (trackDX) = predX1 + relativeXSign_ * predX2; double recHitX_1 = traj1->recHit()->localPosition().x(); double recHitX_2 = traj2->recHit()->localPosition().x(); - (hitDX) = recHitX_1 + relativeXSign_*recHitX_2; + (hitDX) = recHitX_1 + relativeXSign_ * recHitX_2; return true; } @@ -707,35 +667,53 @@ bool HitResol::check2DPartner(unsigned int iidd, const std::vectorgeographicalId().rawId()==partner_iidd) { + for (const auto& tm : traj) { + if (tm.recHit()->geographicalId().rawId() == partner_iidd) { found2DPartner = true; } } return found2DPartner; } -unsigned int HitResol::checkLayer( unsigned int iidd, const TrackerTopology* tTopo) { - StripSubdetector strip=StripSubdetector(iidd); - unsigned int subid=strip.subdetId(); - if (subid == StripSubdetector::TIB) { +unsigned int HitResol::checkLayer(unsigned int iidd, const TrackerTopology* tTopo) { + StripSubdetector strip = StripSubdetector(iidd); + unsigned int subid = strip.subdetId(); + if (subid == StripSubdetector::TIB) { return tTopo->tibLayer(iidd); } - if (subid == StripSubdetector::TOB) { - return tTopo->tobLayer(iidd) + 4 ; + if (subid == StripSubdetector::TOB) { + return tTopo->tobLayer(iidd) + 4; } - if (subid == StripSubdetector::TID) { + if (subid == StripSubdetector::TID) { return tTopo->tidWheel(iidd) + 10; } - if (subid == StripSubdetector::TEC) { - return tTopo->tecWheel(iidd) + 13 ; + if (subid == StripSubdetector::TEC) { + return tTopo->tecWheel(iidd) + 13; } return 0; } +void HitResol::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("lumiScalers", edm::InputTag("scalersRawToDigi")); + desc.add("combinatorialTracks", edm::InputTag("generalTracks")); + desc.add("trajectories", edm::InputTag("generalTracks")); + desc.addUntracked("CompressionSettings", -1); + desc.add("Layer", 0); + desc.add("Debug", false); + desc.addUntracked("addLumi", false); + desc.addUntracked("cutOnTracks", false); + desc.addUntracked("trackMultiplicity", 100); + desc.addUntracked("MomentumCut", 3.); + desc.addUntracked("UsePairsOnly", 1); + descriptions.addWithDefaultLabel(desc); +} + //define this as a plug-in DEFINE_FWK_MODULE(HitResol); diff --git a/CalibTracker/SiStripHitResolution/src/SiStripOverlapHit.cc b/CalibTracker/SiStripHitResolution/src/SiStripOverlapHit.cc new file mode 100644 index 0000000000000..c14b0b5dbdea5 --- /dev/null +++ b/CalibTracker/SiStripHitResolution/src/SiStripOverlapHit.cc @@ -0,0 +1,61 @@ +#include "CalibTracker/SiStripHitResolution/interface/SiStripOverlapHit.h" +#include + +SiStripOverlapHit::SiStripOverlapHit(TrajectoryMeasurement const& measA, TrajectoryMeasurement const& measB) { + // check which hit is closer to the IP + // assign it to hitA_, the other to hitB_ + double rA = measA.recHit()->globalPosition().perp(); + double rB = measB.recHit()->globalPosition().perp(); + if (rA < rB) { + measA_ = measA; + measB_ = measB; + } else { + measA_ = measB; + measB_ = measA; + } +} + +TrajectoryStateOnSurface const& SiStripOverlapHit::trajectoryStateOnSurface(unsigned int hit, bool updated) const { + assert(hit < 2); + switch (hit) { + case 0: + return updated ? measA_.updatedState() : measA_.predictedState(); + case 1: + return updated ? measB_.updatedState() : measB_.predictedState(); + default: + return measA_.updatedState(); + } +} + +double SiStripOverlapHit::getTrackLocalAngle(unsigned int hit) const { + // since x is the precise coordinate and z is pointing out, we want the angle between z and x + return hit ? atan(trajectoryStateOnSurface(hit - 1).localDirection().x() / + trajectoryStateOnSurface(hit - 1).localDirection().z()) + : (getTrackLocalAngle(1) + getTrackLocalAngle(2)) / 2.; +} + +double SiStripOverlapHit::offset(unsigned int hit) const { + assert(hit < 2); + // x is the precise coordinate + return this->hit(hit)->localPosition().x() - trajectoryStateOnSurface(hit, false).localPosition().x(); +} + +double SiStripOverlapHit::shift() const { + // so this is the double difference + return offset(0) - offset(1); +} + +double SiStripOverlapHit::distance(bool fromTrajectory) const { + if (fromTrajectory) { + return (trajectoryStateOnSurface(0, true).globalPosition().basicVector() - + trajectoryStateOnSurface(1, true).globalPosition().basicVector()) + .mag(); + } else { + return (hitA()->globalPosition().basicVector() - hitB()->globalPosition().basicVector()).mag(); + } +} + +GlobalPoint SiStripOverlapHit::position() const { + auto vector = (hitA()->globalPosition().basicVector() + hitB()->globalPosition().basicVector()) / 2.; + return GlobalPoint(vector); +} diff --git a/CalibTracker/SiStripHitResolution/test/BuildFile.xml b/CalibTracker/SiStripHitResolution/test/BuildFile.xml new file mode 100644 index 0000000000000..fa4cf908e6453 --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/BuildFile.xml @@ -0,0 +1,4 @@ + + + + diff --git a/CalibTracker/SiStripHitResolution/test/CPEanconfig.py b/CalibTracker/SiStripHitResolution/test/CPEanconfig.py new file mode 100644 index 0000000000000..03a43d3a518d6 --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/CPEanconfig.py @@ -0,0 +1,92 @@ +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +options = VarParsing.VarParsing() +options.register("isUnitTest", + False, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.bool, + "are we running the unit test") +options.parseArguments() + +process = cms.Process("CPEana") + +### Standard Configurations +process.load('Configuration.StandardSequences.Services_cff') +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.RawToDigi_cff') +process.load('Configuration.StandardSequences.L1Reco_cff') +process.load('Configuration.StandardSequences.Reconstruction_cff') + +### global tag +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run3_data', '') + +### initialize MessageLogger and output report +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.enable = False +process.MessageLogger.SiStripCPEAnalyzer =dict() +process.MessageLogger.cout = cms.untracked.PSet( + enable = cms.untracked.bool(True), + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), + FwkReport = cms.untracked.PSet(limit = cms.untracked.int32(-1), + reportEvery = cms.untracked.int32((10 if options.isUnitTest else 100 )) + ), + SiStripCPEAnalyzer = cms.untracked.PSet( limit = cms.untracked.int32(-1)), + enableStatistics = cms.untracked.bool(True) + ) + +process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) ) + +### Events and data source + +if(options.isUnitTest): + process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10) ) +else: + process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) + +process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/ef6009e4-6857-40a1-9a55-0c702021caad.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/ef6cdbda-400c-4813-b4c7-9dfacd070e08.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f1a88b5f-8573-403e-aa35-0ad6b57125c0.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f1c537d0-2265-403b-84f9-dacb5a63c03f.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f30d1b57-eda6-4836-9ed6-cd683945a1e0.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f3c7f61d-7f6b-4021-b6c2-a15b66e3f375.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f3e11a67-7a78-4f6e-9b9b-b7687ce16c68.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4250ffa-e73e-4e42-baa7-aebd8b169105.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4c4cb8e-5c92-49b4-9fd3-40e09c4cf48a.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4f5e6bd-0a16-4937-a3eb-5b76333d9c4d.root")) + +### Track refitter specific stuff +process.load("RecoVertex.BeamSpotProducer.BeamSpot_cff") +import RecoTracker.TrackProducer.TrackRefitter_cfi +import CommonTools.RecoAlgos.recoTrackRefSelector_cfi +process.mytkselector = CommonTools.RecoAlgos.recoTrackRefSelector_cfi.recoTrackRefSelector.clone() +process.mytkselector.src = 'ALCARECOSiStripCalMinBias' +process.mytkselector.quality = ['highPurity'] +process.mytkselector.min3DLayer = 2 +process.mytkselector.ptMin = 0.5 +process.mytkselector.tip = 1.0 +process.myRefittedTracks = RecoTracker.TrackProducer.TrackRefitter_cfi.TrackRefitter.clone() +process.myRefittedTracks.src= 'mytkselector' +process.myRefittedTracks.NavigationSchool = '' +process.myRefittedTracks.Fitter = 'FlexibleKFFittingSmoother' + +### Analyzer +process.SiStripCPEAnalyzer = cms.EDAnalyzer('SiStripCPEAnalyzer', + tracks = cms.untracked.InputTag("ALCARECOSiStripCalMinBias",""), + trajectories = cms.untracked.InputTag('myRefittedTracks'), + association = cms.untracked.InputTag('myRefittedTracks'), + clusters = cms.untracked.InputTag('ALCARECOSiStripCalMinBias'), + StripCPE = cms.ESInputTag('StripCPEfromTrackAngleESProducer:StripCPEfromTrackAngle')) + +### TFileService: output histogram or ntuple +process.TFileService = cms.Service("TFileService", + fileName = cms.string('histodemo.root')) + +### Finally, put together the sequence +process.p = cms.Path(process.offlineBeamSpot*process.mytkselector+process.myRefittedTracks+process.SiStripCPEAnalyzer) diff --git a/CalibTracker/SiStripHitResolution/test/SiStripHitResol_test.py b/CalibTracker/SiStripHitResolution/test/SiStripHitResol_test.py new file mode 100644 index 0000000000000..965695ab80fa2 --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/SiStripHitResol_test.py @@ -0,0 +1,55 @@ +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +options = VarParsing.VarParsing() +options.register("isUnitTest", + False, + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.bool, + "are we running the unit test") +options.parseArguments() + +process = cms.Process("HitResol") + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.MessageLogger.cerr.FwkReport.reportEvery = 1000 +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run3_data', '') + +InputTagName = "ALCARECOSiStripCalMinBias" +OutputRootFile = "hitresol_ALCARECO_2022F.root" +fileNames=cms.untracked.vstring("/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/ef6009e4-6857-40a1-9a55-0c702021caad.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/ef6cdbda-400c-4813-b4c7-9dfacd070e08.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f1a88b5f-8573-403e-aa35-0ad6b57125c0.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f1c537d0-2265-403b-84f9-dacb5a63c03f.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f30d1b57-eda6-4836-9ed6-cd683945a1e0.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f3c7f61d-7f6b-4021-b6c2-a15b66e3f375.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f3e11a67-7a78-4f6e-9b9b-b7687ce16c68.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4250ffa-e73e-4e42-baa7-aebd8b169105.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4c4cb8e-5c92-49b4-9fd3-40e09c4cf48a.root", + "/store/express/Run2022F/StreamExpress/ALCARECO/SiStripCalMinBias-Express-v1/000/362/167/00000/f4f5e6bd-0a16-4937-a3eb-5b76333d9c4d.root") + +process.source = cms.Source("PoolSource", fileNames=fileNames) + +if(options.isUnitTest): + process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(20)) +else: + process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(10000)) + +process.load("RecoVertex.BeamSpotProducer.BeamSpot_cfi") +#process.load("RecoLocalTracker.SiStripRecHitConverter.StripCPEfromTrackAngle_cfi") +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") +process.refitTracks = process.TrackRefitterP5.clone(src=cms.InputTag(InputTagName)) +process.load("CalibTracker.SiStripHitResolution.SiStripHitResol_cff") +process.anResol.combinatorialTracks = cms.InputTag("refitTracks") +process.anResol.trajectories = cms.InputTag("refitTracks") + +process.TFileService = cms.Service("TFileService", + fileName = cms.string(OutputRootFile) + ) + +process.allPath = cms.Path(process.MeasurementTrackerEvent*process.offlineBeamSpot*process.refitTracks*process.hitresol) diff --git a/CalibTracker/SiStripHitResolution/test/SiStripHitResolutionFromCalibTree_cfg.py b/CalibTracker/SiStripHitResolution/test/SiStripHitResolutionFromCalibTree_cfg.py new file mode 100644 index 0000000000000..bbfc99817d68b --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/SiStripHitResolutionFromCalibTree_cfg.py @@ -0,0 +1,84 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("HitEff") +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_data', '') + +RunNumberBegin = 315264 +RunNumberEnd = 315264 + +process.source = cms.Source("EmptyIOVSource", + firstValue = cms.uint64(RunNumberBegin), + lastValue = cms.uint64(RunNumberEnd), + timetype = cms.string('runnumber'), + interval = cms.uint64(1) +) + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(1)) + +InputFilePath = 'root://cms-xrd-global.cern.ch///eos/cms/store/group/dpg_tracker_strip/comm_tracker/Strip/Calibration/calibrationtree/GR18/calibTree_' +InputFilePathEnd = '.root' + +FileName1 = InputFilePath + str(RunNumberBegin) + InputFilePathEnd + +process.SiStripHitEff = cms.EDAnalyzer("SiStripHitResolFromCalibTree", + CalibTreeFilenames = cms.untracked.vstring(FileName1), + Threshold = cms.double(0.2), + nModsMin = cms.int32(25), + doSummary = cms.int32(0), + #ResXSig = cms.untracked.double(5), + SinceAppendMode = cms.bool(True), + IOVMode = cms.string('Run'), + Record = cms.string('SiStripBadStrip'), + doStoreOnDB = cms.bool(True), + BadModulesFile = cms.untracked.string("BadModules_input.txt"), # default "" no input + AutoIneffModTagging = cms.untracked.bool(True), # default true, automatic limit for each layer to identify inefficient modules + ClusterMatchingMethod = cms.untracked.int32(4), # default 0 case0,1,2,3,4 + ClusterTrajDist = cms.untracked.double(64), # default 64 + StripsApvEdge = cms.untracked.double(10), # default 10 + UseOnlyHighPurityTracks = cms.untracked.bool(True), # default True + SpaceBetweenTrains = cms.untracked.int32(25), # default 25 + UseCommonMode = cms.untracked.bool(False), # default False + ShowEndcapSides = cms.untracked.bool(True), # default True + ShowRings = cms.untracked.bool(True), # default False + ShowTOB6TEC9 = cms.untracked.bool(False), # default False + ShowOnlyGoodModules = cms.untracked.bool(False), # default False + TkMapMin = cms.untracked.double(0.95), # default 0.90 + EffPlotMin = cms.untracked.double(0.90), # default 0.90 + Title = cms.string(' Hit Efficiency ') +) + +process.PoolDBOutputService = cms.Service("PoolDBOutputService", + BlobStreamerName = cms.untracked.string('TBufferBlobStreamingService'), + DBParameters = cms.PSet( + authenticationPath = cms.untracked.string('/afs/cern.ch/cms/DB/conddb') + ), + timetype = cms.untracked.string('runnumber'), + connect = cms.string('sqlite_file:dbfile.db'), + toPut = cms.VPSet(cms.PSet( + record = cms.string('SiStripBadStrip'), + tag = cms.string('SiStripHitEffBadModules') + )) +) + +UnitString = "strip unit" +#UnitString = "cm" + +if UnitString == "cm": + RootFileEnding = "_CM.root" +elif UnitString == "strip unit": + RootFileEnding = "_StripUnit.root" +else: + print('ERROR: Unit must be cm or strip unit') + +RootFileBeginning = "SiStripHitEffHistos_run_" +RootFileName = RootFileBeginning + str(RunNumberBegin) + "_" + str(RunNumberEnd) + RootFileEnding + +process.TFileService = cms.Service("TFileService", + fileName = cms.string(RootFileName) +) + +process.allPath = cms.Path(process.SiStripHitEff) diff --git a/CalibTracker/SiStripHitResolution/test/TestDriver.cpp b/CalibTracker/SiStripHitResolution/test/TestDriver.cpp new file mode 100644 index 0000000000000..2f0e0c40064da --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/TestDriver.cpp @@ -0,0 +1,2 @@ +#include "FWCore/Utilities/interface/TestHelper.h" +RUNTEST() diff --git a/CalibTracker/SiStripHitResolution/test/test_SiStripHitResolution.sh b/CalibTracker/SiStripHitResolution/test/test_SiStripHitResolution.sh new file mode 100755 index 0000000000000..b3a147f04da9c --- /dev/null +++ b/CalibTracker/SiStripHitResolution/test/test_SiStripHitResolution.sh @@ -0,0 +1,26 @@ +#!/bin/bash +function die { echo $1: status $2; exit $2; } + +if [ "${SCRAM_TEST_NAME}" != "" ] ; then + mkdir ${SCRAM_TEST_NAME} + cd ${SCRAM_TEST_NAME} +fi + +echo -e "Testing SiStripHitEfficencyWorker \n\n" +cmsRun ${LOCAL_TEST_DIR}/SiStripHitResol_test.py isUnitTest=True || die 'failed running SiStripHitResol_test.py' $? + +echo -e "Testing CPEanconfig.py \n\n" +cmsRun ${LOCAL_TEST_DIR}/CPEanconfig.py isUnitTest=True || die 'failed running CPEanconfig.py' $? + +echo -e "Testing SiStripHitResolutionFromCalibTree_cfg.py \n\n" +cmsRun ${LOCAL_TEST_DIR}/SiStripHitResolutionFromCalibTree_cfg.py || die 'failed running SiStripHitResolutionFromCalibTree_cfg.py' $? + +### To be implemented + +#echo -e "Testing SiStripHitEfficencyHarvester \n\n" +#cmsRun ${LOCAL_TEST_DIR}/testHitEffHarvester.py isUnitTest=True || die 'failed running testHitEffHarvester.py' $? + +#echo -e " testing tSiStripHitEffFromCalibTree \n\n" +#cmsRun ${LOCAL_TEST_DIR}/testSiStripHitEffFromCalibTree_cfg.py inputFiles=HitEffTree.root runNumber=325172 || die 'failed running testSiStripHitEffFromCalibTree_cfg.py' $? + +echo -e "Done with the tests!"