diff --git a/[refs] b/[refs] index 36aeb268edf68..b03ed3dd58822 100644 --- a/[refs] +++ b/[refs] @@ -1,2 +1,2 @@ --- -"refs/heads/CMSSW_8_1_X": 12b33af54478ad3148d97ade71d9ae645c0dd461 +"refs/heads/CMSSW_8_1_X": b4bb9749585578dd928efe7ebdee3ca8ba75d423 diff --git a/trunk/FWCore/Services/src/InitRootHandlers.cc b/trunk/FWCore/Services/src/InitRootHandlers.cc index e24313f84032d..5c7467fc2b1a3 100644 --- a/trunk/FWCore/Services/src/InitRootHandlers.cc +++ b/trunk/FWCore/Services/src/InitRootHandlers.cc @@ -116,7 +116,8 @@ namespace { (el_location.find("TDecompChol::Solve") != std::string::npos) || (el_location.find("THistPainter::PaintInit") != std::string::npos) || (el_location.find("TUnixSystem::SetDisplay") != std::string::npos) || - (el_location.find("TGClient::GetFontByName") != std::string::npos)) { + (el_location.find("TGClient::GetFontByName") != std::string::npos) || + (el_message.find("nbins is <=0 - set to nbins = 1") != std::string::npos)) { el_severity = SeverityLevel::kInfo; } diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_CosmicData_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_CosmicData_cfi.py index 6e63110184a2c..e715e66bc5681 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_CosmicData_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_CosmicData_cfi.py @@ -13,14 +13,16 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftFromDBAlgo_CosmicData, - segmCleanerMode = cms.int32(1), AlphaMaxPhi = cms.double(100.0), + AlphaMaxTheta = cms.double(100.0), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(100.0), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator and performT0SegCorrection = cms.bool(False), performT0_vdriftSegCorrection = cms.bool(False), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_cfi.py index 3125421d42b6f..55108e457b2a6 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDriftFromDB_cfi.py @@ -13,14 +13,16 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftFromDBAlgo, - segmCleanerMode = cms.int32(2), AlphaMaxPhi = cms.double(1.0), + AlphaMaxTheta = cms.double(0.9), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(2), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(0.9), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator performT0_vdriftSegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_CosmicData_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_CosmicData_cfi.py index 6e6d98e0dde5b..81c9ec10edd10 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_CosmicData_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_CosmicData_cfi.py @@ -13,14 +13,16 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftAlgo_CosmicData, - segmCleanerMode = cms.int32(1), AlphaMaxPhi = cms.double(100.0), + AlphaMaxTheta = cms.double(100.0), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(100.0), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # Parameters for the Updator and T0 fit segment performT0_vdriftSegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_cfi.py index 4380dd0640202..e320fe00297d8 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_LinearDrift_cfi.py @@ -13,14 +13,16 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftAlgo, - segmCleanerMode = cms.int32(2), AlphaMaxPhi = cms.double(1.0), + AlphaMaxTheta = cms.double(0.9), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(2), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(0.9), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator and hit_afterT0_resolution = cms.double(0.03), performT0_vdriftSegCorrection = cms.bool(False), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_ParamDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_ParamDrift_cfi.py index 17ad5ce1e4d41..630724324adf8 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_ParamDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco2DAlgo_ParamDrift_cfi.py @@ -13,14 +13,16 @@ # Parameters for the updator # this is the RecHit1D algo!! DTParametrizedDriftAlgo, - segmCleanerMode = cms.int32(2), AlphaMaxPhi = cms.double(1.0), + AlphaMaxTheta = cms.double(0.9), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(2), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(0.9), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator and performT0SegCorrection = cms.bool(False), performT0_vdriftSegCorrection = cms.bool(False), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_CosmicData_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_CosmicData_cfi.py index 0bae89a5841b0..e9d9bb198022d 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_CosmicData_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_CosmicData_cfi.py @@ -18,17 +18,20 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftFromDBAlgo_CosmicData, - segmCleanerMode = cms.int32(1), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator and performT0SegCorrection = cms.bool(False), performT0_vdriftSegCorrection = cms.bool(False), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_cfi.py index 158c34ba5a191..a0e69e5ecc882 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDriftFromDB_cfi.py @@ -19,14 +19,17 @@ # this is the RecHit1D algo!! DTLinearDriftFromDBAlgo, debug = cms.untracked.bool(False), + # Parameters for the cleaner nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator performT0SegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_CosmicData_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_CosmicData_cfi.py index 46aaab785c05f..ef4cd92a12953 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_CosmicData_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_CosmicData_cfi.py @@ -18,17 +18,20 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftAlgo_CosmicData, - segmCleanerMode = cms.int32(1), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for the Updator and T0 fit segment performT0SegCorrection = cms.bool(False), performT0_vdriftSegCorrection = cms.bool(False), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_cfi.py index 2e7ae632967c2..5af50502cd4d5 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_LinearDrift_cfi.py @@ -20,12 +20,14 @@ DTCombinatorialPatternReco2DAlgo_LinearDrift, debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator performT0SegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_ParamDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_ParamDrift_cfi.py index aa330617b083d..ae8f235f55633 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_ParamDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTCombinatorialPatternReco4DAlgo_ParamDrift_cfi.py @@ -20,12 +20,14 @@ DTParametrizedDriftAlgo, debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator performT0SegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDBLoose_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDBLoose_cfi.py index 76940a1dffb8c..3803e7599b174 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDBLoose_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDBLoose_cfi.py @@ -15,20 +15,20 @@ DTLinearDriftFromDBAlgo, AlphaMaxPhi = cms.double(100.0), AlphaMaxTheta = cms.double(100.), - MaxChi2 = cms.double(4.0), - MaxT0 = cms.double(100.0), - MinT0 = cms.double(-100.0), + MaxChi2 = cms.double(8.0), MaxAllowedHits = cms.uint32(50), debug = cms.untracked.bool(False), + # Parameters for the cleaner segmCleanerMode = cms.int32(2), nSharedHitsMax = cms.int32(2), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator performT0_vdriftSegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), performT0SegCorrection = cms.bool(False), - perform_delta_rejecting = cms.bool(True) + perform_delta_rejecting = cms.bool(False) ), Reco2DAlgoName = cms.string('DTMeantimerPatternReco') ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDB_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDB_cfi.py index 7c96cd043247b..c6c27a10b5784 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDB_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDriftFromDB_cfi.py @@ -15,20 +15,20 @@ DTLinearDriftFromDBAlgo, AlphaMaxPhi = cms.double(1.0), AlphaMaxTheta = cms.double(0.9), - MaxChi2 = cms.double(4.0), - MaxT0 = cms.double(100.0), - MinT0 = cms.double(-100.0), + MaxChi2 = cms.double(8.0), MaxAllowedHits = cms.uint32(50), debug = cms.untracked.bool(False), + # Parameters for the cleaner segmCleanerMode = cms.int32(2), nSharedHitsMax = cms.int32(2), nUnSharedHitsMin = cms.int32(2), + # Parameters for T0 fit segment in the Updator performT0_vdriftSegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), performT0SegCorrection = cms.bool(False), - perform_delta_rejecting = cms.bool(True) + perform_delta_rejecting = cms.bool(False) ), Reco2DAlgoName = cms.string('DTMeantimerPatternReco') ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDrift_cfi.py index 56751e37269ce..b1743701cd55e 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_LinearDrift_cfi.py @@ -13,17 +13,22 @@ # Parameters for the updator # this is the RecHit1D algo!! DTLinearDriftAlgo, - segmCleanerMode = cms.int32(1), AlphaMaxPhi = cms.double(1.0), + AlphaMaxTheta = cms.double(0.1), MaxChi2 = cms.double(8.0), - MaxT0 = cms.double(50.0), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(0.1), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), - MinT0 = cms.double(-10.0) + + # Parameters for T0 fit segment in the Updator + performT0_vdriftSegCorrection = cms.bool(False), + hit_afterT0_resolution = cms.double(0.03), + performT0SegCorrection = cms.bool(False), + perform_delta_rejecting = cms.bool(False) ), Reco2DAlgoName = cms.string('DTMeantimerPatternReco') ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_ParamDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_ParamDrift_cfi.py index ec7aac0da2531..73ffa65d8bdc1 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_ParamDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco2DAlgo_ParamDrift_cfi.py @@ -13,17 +13,22 @@ # Parameters for the updator # this is the RecHit1D algo!! DTParametrizedDriftAlgo, - segmCleanerMode = cms.int32(1), AlphaMaxPhi = cms.double(1.0), + AlphaMaxTheta = cms.double(0.1), MaxChi2 = cms.double(8.0), - MaxT0 = cms.double(50.0), MaxAllowedHits = cms.uint32(50), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - AlphaMaxTheta = cms.double(0.1), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), - MinT0 = cms.double(-10.0) + + # Parameters for T0 fit segment in the Updator + performT0_vdriftSegCorrection = cms.bool(False), + hit_afterT0_resolution = cms.double(0.03), + performT0SegCorrection = cms.bool(False), + perform_delta_rejecting = cms.bool(False) ), Reco2DAlgoName = cms.string('DTMeantimerPatternReco') ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDBLoose_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDBLoose_cfi.py index 652efe4f2a722..a651e643735de 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDBLoose_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDBLoose_cfi.py @@ -19,19 +19,22 @@ # this are the RecSegment2D algo parameters! DTMeantimerPatternReco2DAlgo_LinearDriftFromDBLoose, debug = cms.untracked.bool(False), + # Parameters for the cleaner nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator performT0SegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), performT0_vdriftSegCorrection = cms.bool(False), - perform_delta_rejecting = cms.bool(True) + perform_delta_rejecting = cms.bool(False) ) ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDB_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDB_cfi.py index 3a14a7b4d6758..aa9177dfc9831 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDB_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDriftFromDB_cfi.py @@ -19,19 +19,22 @@ # this are the RecSegment2D algo parameters! DTMeantimerPatternReco2DAlgo_LinearDriftFromDB, debug = cms.untracked.bool(False), + # Parameters for the cleaner nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits # If false the theta segment will be taken from the Event. Caveat: in this case the # event must contain the 2D segments! AllDTRecHits = cms.bool(True), + # Parameters for T0 fit segment in the Updator performT0SegCorrection = cms.bool(False), hit_afterT0_resolution = cms.double(0.03), performT0_vdriftSegCorrection = cms.bool(False), - perform_delta_rejecting = cms.bool(True) + perform_delta_rejecting = cms.bool(False) ) ) diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDrift_cfi.py index 5f4b051c1bd7a..424b95025de61 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_LinearDrift_cfi.py @@ -18,11 +18,13 @@ DTLinearDriftAlgo, # this are the RecSegment2D algo parameters! DTMeantimerPatternReco2DAlgo_LinearDrift, - segmCleanerMode = cms.int32(1), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits diff --git a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_ParamDrift_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_ParamDrift_cfi.py index d502fd4bc0798..909145306e921 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_ParamDrift_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/DTMeantimerPatternReco4DAlgo_ParamDrift_cfi.py @@ -18,11 +18,13 @@ # Parameters for the updator # this is the RecHit1D algo!! DTParametrizedDriftAlgo, - segmCleanerMode = cms.int32(1), + debug = cms.untracked.bool(False), + # Parameters for the cleaner + segmCleanerMode = cms.int32(1), nSharedHitsMax = cms.int32(2), - debug = cms.untracked.bool(False), nUnSharedHitsMin = cms.int32(2), + # the input type. # If true the instructions in setDTRecSegment2DContainer will be schipped and the # theta segment will be recomputed from the 1D rechits diff --git a/trunk/RecoLocalMuon/DTSegment/python/dt4DSegments_ApplyT0Correction_cfi.py b/trunk/RecoLocalMuon/DTSegment/python/dt4DSegments_ApplyT0Correction_cfi.py index 26908ba746698..7517ee23ea703 100644 --- a/trunk/RecoLocalMuon/DTSegment/python/dt4DSegments_ApplyT0Correction_cfi.py +++ b/trunk/RecoLocalMuon/DTSegment/python/dt4DSegments_ApplyT0Correction_cfi.py @@ -8,5 +8,5 @@ debug = cms.untracked.bool(False), performT0_vdriftSegCorrection = cms.bool(True), hit_afterT0_resolution = cms.double(0.03), - perform_delta_rejecting = cms.bool(False) + perform_delta_rejecting = cms.bool(False) ) diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.cc b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.cc index b100a983de4bd..a0d0de9d3d7a0 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.cc @@ -57,7 +57,7 @@ DTCombinatorialExtendedPatternReco::reconstruct(const DTSuperLayer* sl, if(debug) cout << "DTCombinatorialExtendedPatternReco::reconstruct" << endl; theTriedPattern.clear(); edm::OwnVector result; - vector hitsForFit = initHits(sl, pairs); + vector> hitsForFit = initHits(sl, pairs); vector candidates = buildSegments(sl, hitsForFit); @@ -77,9 +77,6 @@ DTCombinatorialExtendedPatternReco::reconstruct(const DTSuperLayer* sl, delete *(cand++); // delete the candidate! } - for (vector::iterator it = hitsForFit.begin(), ed = hitsForFit.end(); - it != ed; ++it) delete *it; - return result; } @@ -93,29 +90,29 @@ void DTCombinatorialExtendedPatternReco::setClusters(const vector +vector> DTCombinatorialExtendedPatternReco::initHits(const DTSuperLayer* sl, const std::vector& hits){ - vector result; + vector> result; for (vector::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) { - result.push_back(new DTHitPairForFit(*hit, *sl, theDTGeometry)); + result.push_back(std::make_shared(*hit, *sl, theDTGeometry)); } return result; } vector DTCombinatorialExtendedPatternReco::buildSegments(const DTSuperLayer* sl, - const std::vector& hits){ + const std::vector>& hits){ - typedef vector hitCont; + typedef vector> hitCont; typedef hitCont::const_iterator hitIter; vector result; if(debug) { cout << "DTCombinatorialExtendedPatternReco::buildSegments: " << sl->id() << " nHits " << hits.size() << endl; - for (vector::const_iterator hit=hits.begin(); + for (vector>::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) cout << **hit<< endl; } @@ -176,7 +173,7 @@ DTCombinatorialExtendedPatternReco::buildSegments(const DTSuperLayer* sl, ((*lastHit)->localPosition(codes[lastLR])-posIni).unit(); // search for other compatible hits, with or without the L/R solved - vector assHits = findCompatibleHits(posIni, dirIni, hits); + vector assHits = findCompatibleHits(posIni, dirIni, hits); if(debug) cout << "compatible hits " << assHits.size() << endl; @@ -234,18 +231,18 @@ DTCombinatorialExtendedPatternReco::buildSegments(const DTSuperLayer* sl, } -vector +vector DTCombinatorialExtendedPatternReco::findCompatibleHits(const LocalPoint& posIni, const LocalVector& dirIni, - const vector& hits) { + const vector>& hits) { if (debug) cout << "Pos: " << posIni << " Dir: "<< dirIni << endl; - vector result; + vector result; // counter to early-avoid double counting in hits pattern vector tried; int nCompatibleHits=0; - typedef vector hitCont; + typedef vector> hitCont; typedef hitCont::const_iterator hitIter; for (hitIter hit=hits.begin(); hit!=hits.end(); ++hit) { pair isCompatible = (*hit)->isCompatible(posIni, dirIni); @@ -276,7 +273,7 @@ DTCombinatorialExtendedPatternReco::findCompatibleHits(const LocalPoint& posIni, tried.push_back(0); continue; // neither is compatible } - result.push_back(AssPoint(*hit, lrcode)); + result.push_back(DTSegmentCand::AssPoint(*hit, lrcode)); } @@ -299,7 +296,7 @@ DTCombinatorialExtendedPatternReco::findCompatibleHits(const LocalPoint& posIni, } DTSegmentExtendedCand* -DTCombinatorialExtendedPatternReco::buildBestSegment(std::vector& hits, +DTCombinatorialExtendedPatternReco::buildBestSegment(std::vector& hits, const DTSuperLayer* sl) { if (debug) cout << "DTCombinatorialExtendedPatternReco::buildBestSegment " << hits.size() << endl; @@ -309,14 +306,14 @@ DTCombinatorialExtendedPatternReco::buildBestSegment(std::vector& hits } // hits with defined LR - vector points; + vector points; // without: I store both L and R, a deque since I need front insertion and // deletion - deque pointsNoLR; + deque> pointsNoLR; // first add only the hits with LR assigned - for (vector::const_iterator hit=hits.begin(); + for (vector::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) { if ((*hit).second != DTEnums::undefLR) { points.push_back(*hit); @@ -371,8 +368,8 @@ DTCombinatorialExtendedPatternReco::buildBestSegment(std::vector& hits } void -DTCombinatorialExtendedPatternReco::buildPointsCollection(vector& points, - deque& pointsNoLR, +DTCombinatorialExtendedPatternReco::buildPointsCollection(vector& points, + deque>& pointsNoLR, vector& candidates, const DTSuperLayer* sl) { @@ -381,11 +378,11 @@ DTCombinatorialExtendedPatternReco::buildPointsCollection(vector& poin cout << "points: " << points.size() << " NOLR: " << pointsNoLR.size()<< endl; } if (pointsNoLR.size()>0) { // still unassociated points! - DTHitPairForFit* unassHit = pointsNoLR.front(); + std::shared_ptr unassHit = pointsNoLR.front(); // try with the right if(debug) cout << "Right hit" << endl; - points.push_back(AssPoint(unassHit, DTEnums::Right)); + points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Right)); pointsNoLR.pop_front(); buildPointsCollection(points, pointsNoLR, candidates, sl); pointsNoLR.push_front((unassHit)); @@ -394,7 +391,7 @@ DTCombinatorialExtendedPatternReco::buildPointsCollection(vector& poin // try with the left if(debug) cout << "Left hit" << endl; - points.push_back(AssPoint(unassHit, DTEnums::Left)); + points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Left)); pointsNoLR.pop_front(); buildPointsCollection(points, pointsNoLR, candidates, sl); pointsNoLR.push_front((unassHit)); @@ -410,7 +407,7 @@ DTCombinatorialExtendedPatternReco::buildPointsCollection(vector& poin } DTSegmentCand::AssPointCont pointsSet; - // for (vector::const_iterator point=points.begin(); + // for (vector::const_iterator point=points.begin(); // point!=points.end(); ++point) pointsSet.insert(points.begin(),points.end()); diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.h b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.h index d253263cb2922..ed22c82f7355a 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.h +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialExtendedPatternReco.h @@ -33,6 +33,7 @@ class DTSegmentExtendedCand; #include "Geometry/DTGeometry/interface/DTGeometry.h" #include "FWCore/Framework/interface/ESHandle.h" +#include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h" /* ====================================================================== */ @@ -68,23 +69,22 @@ class DTCombinatorialExtendedPatternReco : private DTRecSegment2DBaseAlgo { protected: private: - typedef std::pair AssPoint; // create the DTHitPairForFit from the pairs for easy use - std::vector initHits(const DTSuperLayer* sl, + std::vector> initHits(const DTSuperLayer* sl, const std::vector& hits); // search for candidate, starting from pairs of hits in different layers std::vector buildSegments(const DTSuperLayer* sl, - const std::vector& hits); + const std::vector>& hits); // find all the hits compatible with the candidate - std::vector findCompatibleHits(const LocalPoint& pos, + std::vector findCompatibleHits(const LocalPoint& pos, const LocalVector& dir, - const std::vector& hits); + const std::vector>& hits); // build segments from hits collection - DTSegmentExtendedCand* buildBestSegment(std::vector& assHits, + DTSegmentExtendedCand* buildBestSegment(std::vector& assHits, const DTSuperLayer* sl) ; bool checkDoubleCandidates(std::vector& segs, @@ -92,8 +92,8 @@ class DTCombinatorialExtendedPatternReco : private DTRecSegment2DBaseAlgo { /** build collection of compatible hits for L/R hits: the candidates is * updated with the segment candidates found */ - void buildPointsCollection(std::vector& points, - std::deque& pointsNoLR, + void buildPointsCollection(std::vector& points, + std::deque>& pointsNoLR, std::vector& candidates, const DTSuperLayer* sl); diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.cc b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.cc index 038c56346a890..a25ea504e4b5e 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.cc @@ -56,7 +56,7 @@ DTCombinatorialPatternReco::reconstruct(const DTSuperLayer* sl, const std::vector& pairs){ edm::OwnVector result; - vector hitsForFit = initHits(sl, pairs); + vector> hitsForFit = initHits(sl, pairs); vector candidates = buildSegments(sl, hitsForFit); @@ -76,9 +76,6 @@ DTCombinatorialPatternReco::reconstruct(const DTSuperLayer* sl, delete *(cand++); // delete the candidate! } - for (vector::iterator it = hitsForFit.begin(), ed = hitsForFit.end(); - it != ed; ++it) delete *it; - return result; } @@ -88,34 +85,34 @@ void DTCombinatorialPatternReco::setES(const edm::EventSetup& setup){ theUpdator->setES(setup); } -vector +vector> DTCombinatorialPatternReco::initHits(const DTSuperLayer* sl, const std::vector& hits){ - vector result; + vector> result; for (vector::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) { - result.push_back(new DTHitPairForFit(*hit, *sl, theDTGeometry)); + result.push_back(std::make_shared(*hit, *sl, theDTGeometry)); } return result; } vector DTCombinatorialPatternReco::buildSegments(const DTSuperLayer* sl, - const std::vector& hits){ + const std::vector>& hits){ // clear the patterns tried if (debug) { cout << "theTriedPattern.size is " << theTriedPattern.size() << "\n"; } theTriedPattern.clear(); - typedef vector hitCont; + typedef vector> hitCont; typedef hitCont::const_iterator hitIter; vector result; if(debug) { cout << "buildSegments: " << sl->id() << " nHits " << hits.size() << endl; - for (vector::const_iterator hit=hits.begin(); + for (vector>::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) cout << **hit<< endl; } @@ -177,7 +174,7 @@ DTCombinatorialPatternReco::buildSegments(const DTSuperLayer* sl, ((*lastHit)->localPosition(codes[lastLR])-posIni).unit(); // search for other compatible hits, with or without the L/R solved - vector assHits = findCompatibleHits(posIni, dirIni, hits); + vector assHits = findCompatibleHits(posIni, dirIni, hits); if(debug) cout << "compatible hits " << assHits.size() << endl; @@ -232,18 +229,18 @@ DTCombinatorialPatternReco::buildSegments(const DTSuperLayer* sl, } -vector +vector DTCombinatorialPatternReco::findCompatibleHits(const LocalPoint& posIni, const LocalVector& dirIni, - const vector& hits) { + const vector>& hits) { if (debug) cout << "Pos: " << posIni << " Dir: "<< dirIni << endl; - vector result; + vector result; // counter to early-avoid double counting in hits pattern TriedPattern tried; int nCompatibleHits=0; - typedef vector hitCont; + typedef vector> hitCont; typedef hitCont::const_iterator hitIter; for (hitIter hit=hits.begin(); hit!=hits.end(); ++hit) { pair isCompatible = (*hit)->isCompatible(posIni, dirIni); @@ -274,7 +271,7 @@ DTCombinatorialPatternReco::findCompatibleHits(const LocalPoint& posIni, tried.push_back(0); continue; // neither is compatible } - result.push_back(AssPoint(*hit, lrcode)); + result.push_back(DTSegmentCand::AssPoint(*hit, lrcode)); } @@ -313,7 +310,7 @@ DTCombinatorialPatternReco::findCompatibleHits(const LocalPoint& posIni, } DTSegmentCand* -DTCombinatorialPatternReco::buildBestSegment(std::vector& hits, +DTCombinatorialPatternReco::buildBestSegment(std::vector& hits, const DTSuperLayer* sl) { if (hits.size()<3) { //cout << "buildBestSegment: hits " << hits.size()<< endl; @@ -321,14 +318,14 @@ DTCombinatorialPatternReco::buildBestSegment(std::vector& hits, } // hits with defined LR - vector points; + vector points; // without: I store both L and R, a deque since I need front insertion and // deletion - deque pointsNoLR; + deque > pointsNoLR; // first add only the hits with LR assigned - for (vector::const_iterator hit=hits.begin(); + for (vector::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) { if ((*hit).second != DTEnums::undefLR) { points.push_back(*hit); @@ -379,8 +376,8 @@ DTCombinatorialPatternReco::buildBestSegment(std::vector& hits, } void -DTCombinatorialPatternReco::buildPointsCollection(vector& points, - deque& pointsNoLR, +DTCombinatorialPatternReco::buildPointsCollection(vector& points, + deque>& pointsNoLR, vector& candidates, const DTSuperLayer* sl) { @@ -389,11 +386,11 @@ DTCombinatorialPatternReco::buildPointsCollection(vector& points, cout << "points: " << points.size() << " NOLR: " << pointsNoLR.size()<< endl; } if (pointsNoLR.size()>0) { // still unassociated points! - DTHitPairForFit* unassHit = pointsNoLR.front(); + std::shared_ptr unassHit = pointsNoLR.front(); // try with the right if(debug) cout << "Right hit" << endl; - points.push_back(AssPoint(unassHit, DTEnums::Right)); + points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Right)); pointsNoLR.pop_front(); buildPointsCollection(points, pointsNoLR, candidates, sl); pointsNoLR.push_front((unassHit)); @@ -402,7 +399,7 @@ DTCombinatorialPatternReco::buildPointsCollection(vector& points, // try with the left if(debug) cout << "Left hit" << endl; - points.push_back(AssPoint(unassHit, DTEnums::Left)); + points.push_back(DTSegmentCand::AssPoint(unassHit, DTEnums::Left)); pointsNoLR.pop_front(); buildPointsCollection(points, pointsNoLR, candidates, sl); pointsNoLR.push_front((unassHit)); @@ -418,7 +415,7 @@ DTCombinatorialPatternReco::buildPointsCollection(vector& points, } DTSegmentCand::AssPointCont pointsSet; - // for (vector::const_iterator point=points.begin(); + // for (vector::const_iterator point=points.begin(); // point!=points.end(); ++point) pointsSet.insert(points.begin(),points.end()); diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h index 25c42ef76d69f..9ccced7cd89f9 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco.h @@ -32,6 +32,7 @@ class DTSegmentCand; #include "Geometry/DTGeometry/interface/DTGeometry.h" #include "FWCore/Framework/interface/ESHandle.h" +#include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h" /* ====================================================================== */ @@ -66,23 +67,23 @@ class DTCombinatorialPatternReco : public DTRecSegment2DBaseAlgo { private: friend class DTCombinatorialPatternReco4D; - typedef std::pair AssPoint; + // typedef std::pair AssPoint; // create the DTHitPairForFit from the pairs for easy use - std::vector initHits(const DTSuperLayer* sl, + std::vector> initHits(const DTSuperLayer* sl, const std::vector& hits); // search for candidate, starting from pairs of hits in different layers std::vector buildSegments(const DTSuperLayer* sl, - const std::vector& hits); + const std::vector>& hits); // find all the hits compatible with the candidate - std::vector findCompatibleHits(const LocalPoint& pos, + std::vector findCompatibleHits(const LocalPoint& pos, const LocalVector& dir, - const std::vector& hits); + const std::vector>& hits); // build segments from hits collection - DTSegmentCand* buildBestSegment(std::vector& assHits, + DTSegmentCand* buildBestSegment(std::vector& assHits, const DTSuperLayer* sl) ; bool checkDoubleCandidates(std::vector& segs, @@ -90,8 +91,8 @@ class DTCombinatorialPatternReco : public DTRecSegment2DBaseAlgo { /** build collection of compatible hits for L/R hits: the candidates is * updated with the segment candidates found */ - void buildPointsCollection(std::vector& points, - std::deque& pointsNoLR, + void buildPointsCollection(std::vector& points, + std::deque >& pointsNoLR, std::vector& candidates, const DTSuperLayer* sl); private: diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.cc b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.cc index 220171d4cd46e..2b83e9115cdcc 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.cc @@ -37,7 +37,9 @@ DTCombinatorialPatternReco4D::DTCombinatorialPatternReco4D(const ParameterSet& p //do you want the T0 correction? applyT0corr = pset.getParameter("performT0SegCorrection"); - computeT0corr = pset.getUntrackedParameter("computeT0Seg",true); + + computeT0corr = pset.existsAs("computeT0Seg") ? + pset.getParameter("computeT0Seg") : true; // the updator theUpdator = new DTSegmentUpdator(pset); @@ -131,7 +133,7 @@ DTCombinatorialPatternReco4D::reconstruct() { cout << "Reconstructing of the Phi segments" << endl; } - vector pairPhiOwned; + vector> pairPhiOwned; vector resultPhi = buildPhiSuperSegmentsCandidates(pairPhiOwned); if (debug) cout << "There are " << resultPhi.size() << " Phi cand" << endl; @@ -250,15 +252,13 @@ DTCombinatorialPatternReco4D::reconstruct() { // finally delete the candidates! for (vector::iterator phi=resultPhi.begin(); phi!=resultPhi.end(); ++phi) delete *phi; - for (vector::iterator phiPair = pairPhiOwned.begin(); - phiPair!=pairPhiOwned.end(); ++phiPair) delete *phiPair; return result; } -vector DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates(vector &pairPhiOwned){ +vector DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandidates(vector> &pairPhiOwned){ DTSuperLayerId slId; @@ -275,9 +275,9 @@ vector DTCombinatorialPatternReco4D::buildPhiSuperSegmentsCandid const DTSuperLayer *sl = theDTGeometry->superLayer(slId); - vector pairPhi1 = the2DAlgo->initHits(sl,theHitsFromPhi1); + vector> pairPhi1 = the2DAlgo->initHits(sl,theHitsFromPhi1); // same sl!! Since the fit will be in the sl phi 1! - vector pairPhi2 = the2DAlgo->initHits(sl,theHitsFromPhi2); + vector> pairPhi2 = the2DAlgo->initHits(sl,theHitsFromPhi2); // copy the pairPhi2 in the pairPhi1 vector copy(pairPhi2.begin(),pairPhi2.end(),back_inserter(pairPhi1)); diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h index 52b19d2b14471..39deeb2b99c0d 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h +++ b/trunk/RecoLocalMuon/DTSegment/src/DTCombinatorialPatternReco4D.h @@ -61,7 +61,7 @@ class DTCombinatorialPatternReco4D : public DTRecSegment4DBaseAlgo { protected: private: - std::vector buildPhiSuperSegmentsCandidates(std::vector &pairPhiOwned); + std::vector buildPhiSuperSegmentsCandidates(std::vector> &pairPhiOwned); DTRecSegment4D* segmentSpecialZed(const DTRecSegment4D* seg); std::string theAlgoName; diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTLinearFit.cc b/trunk/RecoLocalMuon/DTSegment/src/DTLinearFit.cc index 1adbb669c33b7..5a056de610393 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTLinearFit.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTLinearFit.cc @@ -9,6 +9,7 @@ /* Collaborating Class Header */ /* C++ Headers */ +#include using namespace std; /* ====================================================================== */ @@ -25,32 +26,355 @@ DTLinearFit::~DTLinearFit() { /* Operations */ void DTLinearFit::fit(const vector & x, const vector & y, - int ndat, + int nptfit, const vector & sigy, float& slope, - float& intercept, + float& intercept, + double& chi2, float& covss, float& covii, float& covsi) const { - - float g1 = 0, g2 = 0; - float s11 = 0, s12 = 0, s22 = 0; - for (int i = 0; i != ndat; i++) { - float sy2 = sigy[i] * sigy[i]; - g1 += y[i] / sy2; - g2 += x[i]*y[i] / sy2; + chi2=0; + float sy = 0, sx = 0; + float s11 = 0, sxy = 0, sxx = 0; + for (int i = 0; i != nptfit; i++) { + float sy2 = sigy[i]*sigy[i]; + sy += y[i] / sy2; + sxy += x[i]*y[i] / sy2; s11 += 1. / sy2; - s12 += x[i] / sy2; - s22 += x[i]*x[i] / sy2; + sx += x[i] / sy2; + sxx += x[i]*x[i] / sy2; + } + + float delta = s11*sxx - sx*sx; + + if (delta==0) return; + + intercept = (sy*sxx - sxy*sx) / delta; + slope = (sxy*s11 - sy*sx) / delta; + + covii = sxx / delta; + covss = s11 / delta; + covsi = -sx / delta; + + for (int j=0; j& xfit, + const vector& yfit, + const vector& lfit, + const vector& tfit, + const vector & sigy, + float& aminf, + float& bminf, + float& cminf, + float& vminf, + double& chi2fit, + const bool debug=0) const { + + int nptfit=xfit.size(); + if (nptfitmargin) && (fabs(ycorr-xwire)>margin) && ((yfit[j]-xwire)*(ycorr-xwire)<0)) { +// cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl; +// cout << " XXX hit moved across wire!!!" << endl; + chi2fit=-1.; + return; + } + + if (fabs(ypred-xwire)>2.1 + margin) { +// cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl; +// cout << " XXX segment outside chamber!!!" << endl; + chi2fit=-1.; + return; + } + + if (fabs(ycorr-xwire)>2.1 + margin) { +// cout << " segmentpred: " << ypred << " hit: " << ycorr << " xwire: " << xwire << " yhit: " << yfit[j] << " t0: " << -cminf/0.00543 << endl; +// cout << " XXX hit outside chamber!!!" << endl; + chi2fit=-1.; + return; + } + if (fabs(chi2fit<0.0001)) chi2fit=0; + } //end loop chi2 + } + } else if (npar==4) { + const double det = (a1*a1*(b2*v6 - b6*b6) - a1*(a2*a2*v6 - 2*a2*a6*b6 + a6*a6*b2 + b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)) + + a2*a2*c6*c6 + 2*a2*(a3*(b3*v6 - b6*c6) - a6*b3*c6) + a3*a3*(b6*b6 - b2*v6) + + a6*(2*a3*(b2*c6 - b3*b6) + a6*b3*b3)); + + if (det != 0) { + // computation of a, b, c and v + bminf = (a1*(a2*(b6*d6 - v6*d2) + a6*(b6*d2 - b2*d6) + d1*(b2*v6 - b6*b6)) - a2*(b3*(c6*d6 - v6*d3) + + c6*(b6*d3 - c6*d2)) + a3*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) + + a6*(b2*c6*d3 + b3*(b3*d6 - b6*d3 - c6*d2)) - d1*(b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)))/det; + aminf = - (a1*a1*(b6*d6 - v6*d2) - a1*(a2*(a6*d6 - v6*d1) - a6*a6*d2 + a6*b6*d1 + b3*(c6*d6 - v6*d3) + + c6*(b6*d3 - c6*d2)) + a2*(a3*(c6*d6 - v6*d3) + c6*(a6*d3 - c6*d1)) + a3*a3*(v6*d2 - b6*d6) + + a3*(a6*(b3*d6 + b6*d3 - 2*c6*d2) - d1*(b3*v6 - b6*c6)) - a6*b3*(a6*d3 - c6*d1))/det; + cminf = -(a1*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) + a2*a2*(v6*d3 - c6*d6) + + a2*(a3*(b6*d6 - v6*d2) + a6*(b3*d6 - 2*b6*d3 + c6*d2) - d1*(b3*v6 - b6*c6)) + + a3*(d1*(b2*v6 - b6*b6) - a6*(b2*d6 - b6*d2)) + a6*(a6*(b2*d3 - b3*d2) - d1*(b2*c6 - b3*b6)))/det; + vminf = - (a1*a1*(b2*d6 - b6*d2) - a1*(a2*a2*d6 - a2*(a6*d2 + b6*d1) + a6*b2*d1 + b2*c6*d3 + + b3*(b3*d6 - b6*d3 - c6*d2)) + a2*a2*c6*d3 + a2*(a3*(2*b3*d6 - b6*d3 - c6*d2) - b3*(a6*d3 + c6*d1)) + + a3*a3*(b6*d2 - b2*d6) + a3*(a6*(b2*d3 - b3*d2) + d1*(b2*c6 - b3*b6)) + a6*b3*b3*d1)/det; + + chi2fit = 0.; + for (int j=0; j& xfit, + const vector& yfit, + const vector& lfit, + const int nptfit, + const vector & sigy, + float& aminf, + float& bminf, + float& cminf, + double& chi2fit, + const bool debug=0) const { + + float vminf; + vector tfit( nptfit, 0.); + + fitNpar(3,xfit,yfit,lfit,tfit,sigy,aminf,bminf,cminf,vminf,chi2fit,debug); +} - covii = s22 / d; - covss = s11 / d; - covsi = -s12 / d; +void DTLinearFit::fit4Var( const vector& xfit, + const vector& yfit, + const vector& lfit, + const vector& tfit, + const int nptfit, + float& aminf, + float& bminf, + float& cminf, + float& vminf, + double& chi2fit, + const bool vdrift_4parfit=0, + const bool debug=0) const { + + if (debug) cout << "Entering Fit4Var" << endl; + + const double sigma = 0.0295;// FIXME errors can be inserted .just load them/that is the usual TB resolution value for DT chambers + vector sigy; + + aminf = 0.; + bminf = 0.; + cminf = -999.; + vminf = 0.; + + int nppar = 0; + double chi2fitN2 = -1. ; + double chi2fit3 = -1.; + double chi2fitN3 = -1. ; + double chi2fitN4 = -1.; + float bminf3 = bminf; + float aminf3 = aminf; + float cminf3 = cminf; + int nppar2 = 0; + int nppar3 = 0; + int nppar4 = 0; + + cminf = -999.; + vminf = 0.; + + for (int j=0; j=5) { + fitNpar(4,xfit,yfit,lfit,tfit,sigy,aminf,bminf,cminf,vminf,chi2fit,debug); + + if (vminf!=0) { + nppar = 4; + if (nptfit<=nppar){ + chi2fitN4=-1; + } else{ + chi2fitN4= chi2fit / (nptfit-nppar); + } + } else { + if (nptfit <= nppar) chi2fitN4=-1; + else chi2fitN4 = chi2fit / (nptfit-nppar); + } + + if (fabs(vminf) >= 0.29) { + // for safety and for code construction..dont accept correction on dv/vdrift greater then 0.09 + vminf = 0.; + cminf = cminf3; + aminf = aminf3; + bminf = bminf3; + nppar = 3; + chi2fit = chi2fit3; + } + + } + + if (!vdrift_4parfit){ //if not required explicitly leave the t0 and track step as at step 3 + // just update vdrift value vmin for storing in the segments for monitoring + cminf = cminf3; + aminf = aminf3; + bminf = bminf3; + nppar = 3; + chi2fit = chi2fit3; + } + + nppar4 = nppar; + + } //end nptfit >=3 + + if (debug) { + cout << " dt0= 0 : slope 4 = " << bminf << " pos out = " << aminf <<" chi2fitN4 = " << chi2fitN4 + << " nppar= " << nppar4 << " T0_ev ns= " << cminf/0.00543 <<" delta v = " << vminf < @@ -30,17 +31,19 @@ using namespace std; /* ====================================================================== */ + typedef std::vector> hitCont; + typedef hitCont::const_iterator hitIter; + /// Constructor DTMeantimerPatternReco::DTMeantimerPatternReco(const edm::ParameterSet& pset) : -DTRecSegment2DBaseAlgo(pset), theAlgoName("DTMeantimerPatternReco") + DTRecSegment2DBaseAlgo(pset), + theFitter(new DTLinearFit()), + theAlgoName("DTMeantimerPatternReco") { - theMaxAllowedHits = pset.getParameter("MaxAllowedHits"); // 100 theAlphaMaxTheta = pset.getParameter("AlphaMaxTheta");// 0.1 ; theAlphaMaxPhi = pset.getParameter("AlphaMaxPhi");// 1.0 ; - theMaxT0 = pset.getParameter("MaxT0"); - theMinT0 = pset.getParameter("MinT0"); theMaxChi2 = pset.getParameter("MaxChi2");// 8.0 ; debug = pset.getUntrackedParameter("debug"); theUpdator = new DTSegmentUpdator(pset); @@ -51,6 +54,7 @@ DTRecSegment2DBaseAlgo(pset), theAlgoName("DTMeantimerPatternReco") DTMeantimerPatternReco::~DTMeantimerPatternReco() { delete theUpdator; delete theCleaner; + delete theFitter; } /* Operations */ @@ -59,7 +63,7 @@ DTMeantimerPatternReco::reconstruct(const DTSuperLayer* sl, const std::vector& pairs){ edm::OwnVector result; - vector hitsForFit = initHits(sl, pairs); + std::vector> hitsForFit = initHits(sl, pairs); vector candidates = buildSegments(sl, hitsForFit); @@ -67,19 +71,15 @@ DTMeantimerPatternReco::reconstruct(const DTSuperLayer* sl, while (candupdate(segment); if (debug) cout<<"Reconstructed 2D segments "<<*segment<::iterator it = hitsForFit.begin(), ed = hitsForFit.end(); - it != ed; ++it) delete *it; - return result; } @@ -89,24 +89,22 @@ void DTMeantimerPatternReco::setES(const edm::EventSetup& setup){ theUpdator->setES(setup); } -vector +vector> DTMeantimerPatternReco::initHits(const DTSuperLayer* sl, const std::vector& hits){ - vector result; + hitCont result; for (vector::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) { - result.push_back(new DTHitPairForFit(*hit, *sl, theDTGeometry)); + result.push_back(std::make_shared(*hit, *sl, theDTGeometry)); } return result; } vector DTMeantimerPatternReco::buildSegments(const DTSuperLayer* sl, - const std::vector& hits){ + const std::vector>& hits){ - typedef vector hitCont; - typedef hitCont::const_iterator hitIter; vector result; DTEnums::DTCellSide codes[2]={DTEnums::Left, DTEnums::Right}; @@ -115,9 +113,6 @@ DTMeantimerPatternReco::buildSegments(const DTSuperLayer* sl, for (hitIter hit=hits.begin(); hit!=hits.end(); ++hit) cout << **hit<< " wire: " << (*hit)->id() << " DigiTime: " << (*hit)->digiTime() << endl; } - // 10-Mar-2004 SL - // put a protection against heavily populated chambers, for which the segment - // building could lead to infinite memory usage... if (hits.size() > theMaxAllowedHits ) { if(debug) { cout << "Warning: this SuperLayer " << sl->id() << " has too many hits : " @@ -134,8 +129,6 @@ DTMeantimerPatternReco::buildSegments(const DTSuperLayer* sl, else // Phi SL DAlphaMax=theAlphaMaxPhi; - vector usedHits; - // get two hits in different layers and see if there are other hits // compatible with them for (hitCont::const_iterator firstHit=hits.begin(); firstHit!=hits.end(); @@ -166,37 +159,36 @@ DTMeantimerPatternReco::buildSegments(const DTSuperLayer* sl, float DAlpha=fabs(gvec.theta()-gvecIP.theta()); if (DAlpha>DAlphaMax) continue; - if(debug) { - cout << "Selected hit pair:" << endl; - cout << " First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << " Side: " << firstLR << " DigiTime: " << (*firstHit)->digiTime() << endl; - cout << " Last " << *(*lastHit) << " Layer Id: " << (*lastHit)->id().layerId() << " Side: " << lastLR << " DigiTime: " << (*lastHit)->digiTime() << endl; - } +// if(debug) { +// cout << "Selected hit pair:" << endl; +// cout << " First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << " Side: " << firstLR << " DigiTime: " << (*firstHit)->digiTime() << endl; +// cout << " Last " << *(*lastHit) << " Layer Id: " << (*lastHit)->id().layerId() << " Side: " << lastLR << " DigiTime: " << (*lastHit)->digiTime() << endl; +// } - vector assHits; + vector assHits; // create a candidate hit list - assHits.push_back(AssPoint(*firstHit,codes[firstLR])); - assHits.push_back(AssPoint(*lastHit,codes[lastLR])); + assHits.push_back(DTSegmentCand::AssPoint(*firstHit,codes[firstLR])); + assHits.push_back(DTSegmentCand::AssPoint(*lastHit,codes[lastLR])); // run hit adding/segment building maxfound = 3; - addHits(sl,assHits,hitsForFit,result, usedHits); + addHits(sl,assHits,hitsForFit,result); } } } - } + } - // now I have a couple of segment hypotheses, should check for ghost + // now I have a couple of segment hypotheses, should check for ghosts if (debug) { cout << "Result (before cleaning): " << result.size() << endl; - for (vector::const_iterator seg=result.begin(); seg!=result.end(); ++seg) - cout << *(*seg) << endl; + for (vector::const_iterator seg=result.begin(); seg!=result.end(); ++seg) cout << *(*seg) << endl; } + result = theCleaner->clean(result); + if (debug) { cout << "Result (after cleaning): " << result.size() << endl; - for (vector::const_iterator seg=result.begin(); - seg!=result.end(); ++seg) - cout << *(*seg) << endl; + for (vector::const_iterator seg=result.begin(); seg!=result.end(); ++seg) cout << *(*seg) << endl; } return result; @@ -205,37 +197,41 @@ DTMeantimerPatternReco::buildSegments(const DTSuperLayer* sl, void -DTMeantimerPatternReco::addHits(const DTSuperLayer* sl, vector& assHits, const vector& hits, - vector &result, vector& usedHits) { +DTMeantimerPatternReco::addHits(const DTSuperLayer* sl, vector& assHits, const vector>& hits, + vector &result) { - typedef vector hitCont; double chi2l,chi2r,t0_corrl,t0_corrr; bool foundSomething = false; - if (debug) - cout << "DTMeantimerPatternReco::addHit " << endl << " Picked " << assHits.size() << " hits, " << hits.size() << " left." << endl; +// if (debug) +// cout << "DTMeantimerPatternReco::addHit " << endl << " Picked " << assHits.size() << " hits, " << hits.size() << " left." << endl; if (assHits.size()+hits.size()id() << endl; + +// if (debug) { +// int nHits=assHits.size()+1; +// cout << " Trying B: " << **hit<< " wire: " << (*hit)->id() << endl; +// printPattern(assHits,*hit); +// } - assHits.push_back(AssPoint(*hit, DTEnums::Left)); - bool left_ok=fitWithT0(assHits, chi2l, t0_corrl,0); + assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Left)); + std::unique_ptr left_seg = fitWithT0(sl,assHits, chi2l, t0_corrl,0); assHits.pop_back(); +// if (debug) +// cout << " Left: t0= " << t0_corrl << " chi2/nHits= " << chi2l << "/" << nHits << " ok: " << left_ok << endl; - assHits.push_back(AssPoint(*hit, DTEnums::Right)); - bool right_ok=fitWithT0(assHits, chi2r, t0_corrr,0); + assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Right)); + std::unique_ptr right_seg = fitWithT0(sl,assHits, chi2r, t0_corrr,0); assHits.pop_back(); +// if (debug) +// cout << " Right: t0= " << t0_corrr << " chi2/nHits= " << chi2r << "/" << nHits << " ok: " << right_ok << endl; - if (debug) { - int nHits=assHits.size()+1; - cout << " Left: t0_corr = " << t0_corrl << "ns chi2/nHits = " << chi2l << "/" << nHits << " ok: " << left_ok << endl; - cout << " Right: t0_corr = " << t0_corrr << "ns chi2/nHits = " << chi2r << "/" << nHits << " ok: " << right_ok << endl; - } - + bool left_ok=(left_seg)?true:false; + bool right_ok=(right_seg)?true:false; + if (!left_ok && !right_ok) continue; foundSomething = true; @@ -253,46 +249,40 @@ DTMeantimerPatternReco::addHits(const DTSuperLayer* sl, vector& assHit if (chi2r seg = fitWithT0(sl,assHits, chi2l, t0_corrl,debug); + if (!seg) return; + + if (!seg->good()) { +// if (debug) cout << " Segment not good() - skipping" << endl; + return; + } - DTSegmentCand::AssPointCont pointsSet; - pointsSet.insert(assHits.begin(),assHits.end()); - DTSegmentCand* seg = new DTSegmentCand(pointsSet,sl); - theUpdator->fit(seg); - - if (seg) { - for (vector::const_iterator hiti = assHits.begin()+1; hiti != assHits.end()-1; hiti++) - usedHits.push_back(*hiti); - - if (assHits.size()>maxfound) maxfound = assHits.size(); - seg->setChi2(chi2l); // we need to set the chi^2 so that the cleaner can pick the best segments - if (debug) cout << "Segment built: " << *seg<< endl; - if (checkDoubleCandidates(result,seg)) { - result.push_back(seg); - if (debug) cout << " Result is now " << result.size() << endl; - } else { - if (debug) cout << " Exists - skipping" << endl; - delete seg; - } + if (assHits.size()>maxfound) maxfound = assHits.size(); + if (debug) cout << endl << " Seg t0= " << t0_corrl << *seg<< endl; + + if (checkDoubleCandidates(result,seg.get())) { + result.push_back(seg.release()); + if (debug) cout << " Result is now " << result.size() << endl; + } else { + if (debug) cout << " Exists - skipping" << endl; } } @@ -333,147 +323,44 @@ DTMeantimerPatternReco::geometryFilter( const DTWireId first, const DTWireId sec } -bool -DTMeantimerPatternReco::fitWithT0(const vector &assHits, double &chi2, double &t0_corr, const bool fitdebug) +std::unique_ptr +DTMeantimerPatternReco::fitWithT0(const DTSuperLayer* sl, const vector &assHits, double &chi2, double &t0_corr, const bool fitdebug) { - typedef vector < pair > hitCoord; - double a,b,coordError,x,y; - double sx=0,sy=0,sxy=0,sxx=0,ssx=0,ssy=0,s=0,ss=0; - int leftHits=0,rightHits=0; - - if (assHits.size()<3) return false; - - // I'm assuming the single hit error is the same for all hits... - coordError=((*(assHits.begin())).first)->localPositionError().xx(); - - for (vector::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) { - if (coordError!=((*(assHits.begin())).first)->localPositionError().xx()) - cout << " DTMeantimerPatternReco: Warning! Hit errors different within one segment!" << endl; - - x=((*assHit).first)->localPosition((*assHit).second).z(); - y=((*assHit).first)->localPosition((*assHit).second).x(); - - sx+=x; - sy+=y; - sxy+=x*y; - sxx+=x*x; - s++; - - if ((*assHit).second==DTEnums::Left) { - leftHits++; - ssx+=x; - ssy+=y; - ss++; - } else { - rightHits++; - ssx-=x; - ssy-=y; - ss--; - } - } - - if (fitdebug && debug) - cout << " DTMeantimerPatternReco::fitWithT0 Left hits: " << leftHits << " Right hits: " << rightHits << endl; - - if (leftHits && rightHits) { - - double delta = ss*ss*sxx+s*sx*sx+s*ssx*ssx-s*s*sxx-2*ss*sx*ssx; - t0_corr=0.; - - if (delta) { - a=(ssy*s*ssx+sxy*ss*ss+sy*sx*s-sy*ss*ssx-ssy*sx*ss-sxy*s*s)/delta; - b=(ssx*sy*ssx+sxx*ssy*ss+sx*sxy*s-sxx*sy*s-ssx*sxy*ss-sx*ssy*ssx)/delta; - t0_corr=(ssx*s*sxy+sxx*ss*sy+sx*sx*ssy-sxx*s*ssy-sx*ss*sxy-ssx*sx*sy)/delta; - } else return false; - } else { - double d = s*sxx - sx*sx; - b = (sxx*sy- sx*sxy)/ d; - a = (s*sxy - sx*sy) / d; - t0_corr=0; - } - -// Calculate the chi^2 of the hits AFTER t0 correction - - double chi,chi2_not0; - chi2=0; - chi2_not0=0; + // create a DTSegmentCand + DTSegmentCand::AssPointCont pointsSet; + pointsSet.insert(assHits.begin(),assHits.end()); + std::unique_ptr seg(new DTSegmentCand(pointsSet,sl)); - for (vector::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) { - x=((*assHit).first)->localPosition((*assHit).second).z(); - y=((*assHit).first)->localPosition((*assHit).second).x(); - - chi=y-a*x-b; - chi2_not0+=chi*chi/coordError; + // perform the 3 parameter fit on the segment candidate + theUpdator->fit(seg.get(),1,fitdebug); - if ((*assHit).second==DTEnums::Left) chi-=t0_corr; - else chi+=t0_corr; - - chi2+=chi*chi/coordError; - } - - // For 3-hit segments ignore timing information - if (assHits.size()<4) { - chi2=chi2_not0; -// if (chi2theMinT0)) || (assHits.size()==4)) { - cout << " a = " << a << " b = " << b << endl; - for (vector::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) { - x=((*assHit).first)->localPosition((*assHit).second).z(); - y=((*assHit).first)->localPosition((*assHit).second).x(); - cout << " z= " << x << " x= " << y; - if ((*assHit).second==DTEnums::Left) cout << " x_corr= " << y+t0_corr*0.00543; - else cout << " x_corr= " << y-t0_corr*0.00543; - cout << " seg= " << a*x+b << endl; - } - } + chi2 = seg->chi2(); + t0_corr = seg->t0(); + + // sanity check - drop segment candidates with a failed fit + // for a 3-par fit this includes segments with hits after the calculated t0 correction ending up + // beyond the chamber walls or on the other side of the wire + if (chi2==-1.) return nullptr; + + // at this point we keep all 3-hit segments that passed the above check + if (assHits.size()==3) return seg; + + // for segments with no t0 information we impose a looser chi2 cut + if (t0_corr==0) { + if (chi2<200.) return seg; + else return nullptr; } - if ((chi2/(assHits.size()-2)theMinT0)) return true; - else return false; + // cut on chi2/ndof of the segment candidate + if ((chi2/(assHits.size()-3) > &hits) { - double s=0,sx=0,sy=0,x,y; - double sxx=0,sxy=0; - - a=b=0; - if (hits.size()==0) return; - - if (hits.size()==1) { - b=(*(hits.begin())).second; - } else { - for (unsigned int i = 0; i != hits.size(); i++) { - x=hits[i].first; - y=hits[i].second; - sy += y; - sxy+= x*y; - s += 1.; - sx += x; - sxx += x*x; - } - // protect against a vertical line??? - - double d = s*sxx - sx*sx; - b = (sxx*sy- sx*sxy)/ d; - a = (s*sxy - sx*sy) / d; - } -} bool DTMeantimerPatternReco::checkDoubleCandidates(vector& cands, - DTSegmentCand* seg) { + DTSegmentCand* seg) { for (vector::iterator cand=cands.begin(); cand!=cands.end(); ++cand) { if (*(*cand)==*seg) return false; @@ -482,3 +369,26 @@ DTMeantimerPatternReco::checkDoubleCandidates(vector& cands, } return true; } + + + +void +DTMeantimerPatternReco::printPattern( vector& assHits, const DTHitPairForFit* hit) { + + char mark[26]={". . . . . . . . . . . . "}; + int wire[12]={0,0,0,0,0,0,0,0,0,0,0,0}; + + for (vector::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) { + int lay = (((*assHit).first)->id().superlayerId().superLayer()-1)*4 + ((*assHit).first)->id().layerId().layer()-1; + wire[lay]= ((*assHit).first)->id().wire(); + if ((*assHit).second==DTEnums::Left) mark[lay*2]='L'; else mark[lay*2]='R'; + } + + int lay = ((*hit).id().superlayerId().superLayer()-1)*4 + (*hit).id().layerId().layer() - 1; + wire[lay]= (*hit).id().wire(); + mark[lay*2]='*'; + + cout << " " << mark << endl << " "; + for (int i=0; i<12; i++) if (wire[i]) cout << setw(2) << wire[i]; else cout << " "; + cout << endl; +} diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco.h b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco.h index 572e296a53824..8ed65cc9c805d 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco.h +++ b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco.h @@ -25,6 +25,7 @@ class DTSegmentUpdator; class DTSegmentCleaner; class DTHitPairForFit; class DTSegmentCand; +class DTLinearFit; /* C++ Headers */ #include @@ -33,11 +34,13 @@ class DTSegmentCand; #include "Geometry/DTGeometry/interface/DTGeometry.h" #include "FWCore/Framework/interface/ESHandle.h" +#include "RecoLocalMuon/DTSegment/src/DTSegmentCand.h" /* ====================================================================== */ /* Class DTMeantimerPatternReco Interface */ + class DTMeantimerPatternReco : public DTRecSegment2DBaseAlgo { public: @@ -65,37 +68,42 @@ class DTMeantimerPatternReco : public DTRecSegment2DBaseAlgo { protected: private: + DTLinearFit* theFitter; // the linear fitter + friend class DTMeantimerPatternReco4D; - typedef std::pair AssPoint; + // typedef std::pair AssPoint; // create the DTHitPairForFit from the pairs for easy use - std::vector initHits(const DTSuperLayer* sl, - const std::vector& hits); + std::vector> initHits(const DTSuperLayer* sl, + const std::vector& hits); // search for candidate, starting from pairs of hits in different layers std::vector buildSegments(const DTSuperLayer* sl, - const std::vector& hits); + const std::vector>& hits); // try adding more hits to a candidate void addHits(const DTSuperLayer* sl, - std::vector& assHits, - const std::vector& hits, - std::vector &result, - std::vector& usedHits); + std::vector& assHits, + const std::vector>& hits, + std::vector &result); // fit a set of left/right hits, calculate t0 and chi^2 - bool fitWithT0(const std::vector &assHits, double &chi2, double &t0_corr, const bool fitdebug); + std::unique_ptr fitWithT0(const DTSuperLayer* sl, + const std::vector &assHits, + double &chi2, + double &t0_corr, + const bool fitdebug); // check if two hist can be considered in one segment (come from different layers, not too far away etc.) bool geometryFilter( const DTWireId first, const DTWireId second ) const; - // a generic least-square fit to a set of points - void rawFit(double &a, double &b, const std::vector< std::pair > &hits); - bool checkDoubleCandidates(std::vector& segs, DTSegmentCand* seg); + void printPattern( std::vector& assHits, const DTHitPairForFit* hit); + + private: std::string theAlgoName; @@ -103,8 +111,6 @@ class DTMeantimerPatternReco : public DTRecSegment2DBaseAlgo { double theAlphaMaxTheta; double theAlphaMaxPhi; double theMaxChi2; - double theMaxT0; - double theMinT0; bool debug; DTSegmentUpdator* theUpdator; // the updator and fitter DTSegmentCleaner* theCleaner; // the cleaner diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.cc b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.cc index 16f11a35bbc73..751cb19e100ad 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.cc @@ -134,7 +134,7 @@ DTMeantimerPatternReco4D::reconstruct(){ cout << "Reconstructing Phi segments"< pairPhiOwned; + vector> pairPhiOwned; vector resultPhi = buildPhiSuperSegmentsCandidates(pairPhiOwned); if (debug) cout << "There are " << resultPhi.size() << " Phi cand" << endl; @@ -163,7 +163,7 @@ DTMeantimerPatternReco4D::reconstruct(){ } // Now I want to build the concrete DTRecSegment4D. - if(debug) cout<<"Building the concrete DTRecSegment4D"<::const_iterator phi=resultPhi.begin(); phi!=resultPhi.end(); ++phi) { @@ -248,15 +248,13 @@ DTMeantimerPatternReco4D::reconstruct(){ // finally delete the candidates! for (vector::iterator phi=resultPhi.begin(); phi!=resultPhi.end(); ++phi) delete *phi; - for (vector::iterator phiPair = pairPhiOwned.begin(); - phiPair!=pairPhiOwned.end(); ++phiPair) delete *phiPair; return result; } -vector DTMeantimerPatternReco4D::buildPhiSuperSegmentsCandidates(vector &pairPhiOwned){ +vector DTMeantimerPatternReco4D::buildPhiSuperSegmentsCandidates(vector> &pairPhiOwned){ DTSuperLayerId slId; @@ -273,9 +271,9 @@ vector DTMeantimerPatternReco4D::buildPhiSuperSegmentsCandidates const DTSuperLayer *sl = theDTGeometry->superLayer(slId); - vector pairPhi1 = the2DAlgo->initHits(sl,theHitsFromPhi1); + vector> pairPhi1 = the2DAlgo->initHits(sl,theHitsFromPhi1); // same sl!! Since the fit will be in the sl phi 1! - vector pairPhi2 = the2DAlgo->initHits(sl,theHitsFromPhi2); + vector> pairPhi2 = the2DAlgo->initHits(sl,theHitsFromPhi2); // copy the pairPhi2 in the pairPhi1 vector copy(pairPhi2.begin(),pairPhi2.end(),back_inserter(pairPhi1)); diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.h b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.h index 7ab28a8b568d8..5f35d0dbea71c 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.h +++ b/trunk/RecoLocalMuon/DTSegment/src/DTMeantimerPatternReco4D.h @@ -61,10 +61,9 @@ class DTMeantimerPatternReco4D : public DTRecSegment4DBaseAlgo { protected: private: - std::vector buildPhiSuperSegmentsCandidates(std::vector &pairPhiOwned); + std::vector buildPhiSuperSegmentsCandidates(std::vector> &pairPhiOwned); DTRecSegment4D* segmentSpecialZed(DTRecSegment4D* seg); - std::string theAlgoName; bool debug; diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.cc b/trunk/RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.cc index aa1d8e699eadd..3b2217a8c8469 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTRefitAndCombineReco4D.cc @@ -186,7 +186,7 @@ vector DTRefitAndCombineReco4D::refitSuperSegments(){ DTChamberRecSegment2D superPhi(chId,recHitsSeg2DPhi1); // refit it! - theUpdator->fit(&superPhi); + theUpdator->fit(&superPhi,0); // cut on the chi^2 if (superPhi.chi2() > theMaxChi2forPhi) diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCand.cc b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCand.cc index cfc8244682604..090cc9b65db50 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCand.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCand.cc @@ -58,7 +58,7 @@ bool DTSegmentCand::operator<(const DTSegmentCand& seg){ return (nHits() hit, DTEnums::DTCellSide code) { theHits.insert(AssPoint(hit,code)); } @@ -107,6 +107,7 @@ DTSegmentCand::conflictingHitPairs(const DTSegmentCand& seg) const{ bool DTSegmentCand::good() const { + // std::cout << NDOF() << " " << chi2()/NDOF() << " " << nHits() << std::endl; if(NDOF() == 0) return false; if(chi2()/NDOF() > chi2max || nHits() < nHitsMin) return false; @@ -117,17 +118,21 @@ bool DTSegmentCand::good() const bool DTSegmentCand::hitsShareLayer() const { - std::vector layerN; + int layerN[20]; + int i=0; + + // we don't expect so many 1D hits, if such a segment arrives just drop it + if (theHits.size()>20) return false; for(DTSegmentCand::AssPointCont::iterator assHit=theHits.begin(); assHit!=theHits.end(); ++assHit) { - layerN.push_back((*assHit).first->id().layerId().layer()); - - //std::cout << (*assHit).first->id().layerId().layer() << std::endl; + layerN[i]=(*assHit).first->id().layerId().layer()+10*(*assHit).first->id().superlayerId().superlayer(); + i++; +// std::cout << (*assHit).first->id().layerId().layer()+10*(*assHit).first->id().superlayerId().superlayer()) << std::endl; } - for(int i=0;i<(int)layerN.size();i++){ - for(int j=i+1;j<(int)layerN.size();j++){ + for(int i=0;i<(int)theHits.size();i++){ + for(int j=0;j AssPoint; + typedef std::pair, DTEnums::DTCellSide> AssPoint; typedef std::set AssPointCont; /// Constructor @@ -60,10 +61,13 @@ class DTSegmentCand{ virtual unsigned int nHits() const { return theHits.size(); } /// the chi2 (NOT chi2/NDOF) of the fit - virtual double chi2() const {return theChi2; } + virtual double chi2() const { return theChi2; } /// the chi2/NDOF of the fit - virtual double chi2ndof() const {return theChi2/(nHits()-2.); } + virtual double chi2ndof() const { return theChi2/(nHits()-2.); } + + /// the t0 of the segment + virtual double t0() const { return thet0; } /// equality operator based on position, direction, chi2 and nHits virtual bool operator==(const DTSegmentCand& seg); @@ -89,16 +93,19 @@ class DTSegmentCand{ virtual void setPosition(LocalPoint& pos) { thePosition=pos; } /// set direction - virtual void setDirection(LocalVector& dir) { theDirection = dir ; } + virtual void setDirection(LocalVector& dir) { theDirection = dir; } /// add hits to the hit list. - virtual void add(DTHitPairForFit* hit, DTEnums::DTCellSide code) ; + virtual void add(std::shared_ptr hit, DTEnums::DTCellSide code) ; /// remove hit from the candidate virtual void removeHit(AssPoint hit) ; /// set chi2 - virtual void setChi2(double& chi2) { theChi2 = chi2 ;} + virtual void setChi2(double& chi2) { theChi2 = chi2; } + + /// set t0 + virtual void sett0(double& t0) { thet0 = t0; } /// number of shared hit pair with other segment candidate virtual int nSharedHitPairs(const DTSegmentCand& seg) const; @@ -124,25 +131,25 @@ class DTSegmentCand{ /// convert this DTSegmentCand into a DTChamberRecSegment2D operator DTChamberRecSegment2D*() const; - struct AssPointLessZ : public std::binary_function { public: bool operator()(const AssPoint& pt1, const AssPoint& pt2) const ; }; + private: const DTSuperLayer* theSL; // the SL LocalPoint thePosition; // in SL frame LocalVector theDirection; // in SL frame double theChi2; // chi2 of the fit + double thet0; // the t0 offset /// mat[1][1]=sigma (dx/dz) /// mat[2][2]=sigma (x) /// mat[1][2]=cov(dx/dz,x) AlgebraicSymMatrix theCovMatrix; // the covariance matrix - AssPointCont theHits; // the used hits protected: diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCleaner.cc b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCleaner.cc index d4f762f69de72..14a2bcb314dc2 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCleaner.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentCleaner.cc @@ -36,9 +36,10 @@ DTSegmentCleaner::~DTSegmentCleaner() { /* Operations */ vector DTSegmentCleaner::clean(const std::vector& inputCands) const { if (inputCands.size()<2) return inputCands; - // cout << "[DTSegmentCleaner] # of candidates: " << inputCands.size() << endl; + //cout << "[DTSegmentCleaner] # of candidates: " << inputCands.size() << endl; vector result = solveConflict(inputCands); + //cout << "[DTSegmentCleaner] to ghostbuster: " << result.size() << endl; result = ghostBuster(result); return result; @@ -55,7 +56,7 @@ vector DTSegmentCleaner::solveConflict(const std::vector::const_iterator cand2 = cand+1 ; cand2!=inputCands.end() ; ++cand2) { DTSegmentCand::AssPointCont confHits=(*cand)->conflictingHitPairs(*(*cand2)); - + if (confHits.size()) { ///treatment of LR ambiguity cases: 1 chooses the best chi2 /// 2 chooses the smaller angle @@ -91,7 +92,6 @@ vector DTSegmentCleaner::solveConflict(const std::vector DAlpha2) ? (*cand) : (*cand2); - } for (DTSegmentCand::AssPointCont::const_iterator cHit=confHits.begin() ; @@ -113,10 +113,9 @@ vector DTSegmentCleaner::solveConflict(const std::vector::const_iterator cand=inputCands.begin(); while ( cand < inputCands.end() ) { - if ((*cand)->good()) result.push_back(*cand); + if ((*cand)->good()) result.push_back(*cand); else { vector::const_iterator badCand=cand; delete *badCand; @@ -134,14 +133,20 @@ DTSegmentCleaner::ghostBuster(const std::vector& inputCands) con for (vector::const_iterator cand2=cand+1; cand2!=inputCands.end(); ++cand2) { unsigned int nSharedHits=(*cand)->nSharedHitPairs(*(*cand2)); - // cout << "Sharing " << (**cand) << " " << (**cand2) << " " << nSharedHits - // << " (first or second) " << ((**cand)<(**cand2)) << endl; + //cout << "Sharing " << (**cand) << " " << (**cand2) << " " << nSharedHits + // << " (first or second) " << ((**cand)<(**cand2)) << endl; if ((nSharedHits==((*cand)->nHits())) && (nSharedHits==((*cand2)->nHits())) &&(fabs((*cand)->chi2()-(*cand2)->chi2())<0.1) &&(segmCleanerMode==3)) { continue; } + + if (((*cand2)->nHits()==3 || (*cand2)->nHits()==3) + &&(fabs((*cand)->chi2()-(*cand2)->chi2())<0.0001)) + { + continue; + } // remove the worst segment if too many shared hits or too few unshared if ((int)nSharedHits >= nSharedHitsMax || diff --git a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentExtendedCand.cc b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentExtendedCand.cc index 033773697ae71..fddd65b308a84 100644 --- a/trunk/RecoLocalMuon/DTSegment/src/DTSegmentExtendedCand.cc +++ b/trunk/RecoLocalMuon/DTSegment/src/DTSegmentExtendedCand.cc @@ -30,8 +30,8 @@ bool DTSegmentExtendedCand::isCompatible(const DTSegmentExtendedCand::DTSLRecClu LocalPoint posAtSL = position()+direction()*(clus.pos.z()-position().z())/cos(direction().theta()); // cout << "pos :" << clus.pos << " posAtSL " << posAtSL << endl; - constexpr float errScaleFact=10.; - constexpr float minError=25.; // (cm) + static float errScaleFact=10.; + static float minError=25.; // (cm) // cout << "clus.err.xx() " << clus.err << endl; return fabs((posAtSL-clus.pos).x())hasPhi(); const bool hasZed = seg->hasZed(); @@ -98,18 +98,41 @@ void DTSegmentUpdator::update(DTRecSegment4D* seg, const bool calcT0) const { } void DTSegmentUpdator::update(DTRecSegment2D* seg) const { + if(debug) cout << "[DTSegmentUpdator] Starting to update the DTRecSegment2D" << endl; GlobalPoint pos = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localPosition()); GlobalVector dir = (theGeom->idToDet(seg->geographicalId()))->toGlobal(seg->localDirection()); updateHits(seg,pos,dir); - fit(seg); + fit(seg,1); } void DTSegmentUpdator::fit(DTRecSegment4D* seg) const { + if(debug) cout << "[DTSegmentUpdator] Fit DTRecSegment4D:" << endl; // after the update must refit the segments - if(seg->hasPhi()) fit(seg->phiSegment()); - if(seg->hasZed()) fit(seg->zSegment()); + + if (debug) { + if(seg->hasPhi()) cout << " 4D Segment contains a Phi segment. t0= " << seg->phiSegment()->t0() << " chi2= " << seg->phiSegment()->chi2() << endl; + if(seg->hasZed()) cout << " 4D Segment contains a Zed segment. t0= " << seg->zSegment()->t0() << " chi2= " << seg->zSegment()->chi2() << endl; + } + + + // use the 2-par fit UNLESS the 4D segment only has the Phi projection or is out of time by more than 20ns + // - in these cases use the 3-par fit + if(seg->hasPhi()) { + if(seg->hasZed()) { + + // fit in-time Phi segments with the 2par fit and out-of-time segments with the 3par fit + if (fabs(seg->phiSegment()->t0())<40.) { + fit(seg->phiSegment(),0); + fit(seg->zSegment(),0); + } else { + fit(seg->phiSegment(),1); + fit(seg->zSegment(),1); + } + } else fit(seg->phiSegment(),1); + } else fit(seg->zSegment(),1); + const DTChamber* theChamber = theGeom->chamber(seg->chamberId()); if(seg->hasPhi() && seg->hasZed() ) { @@ -149,9 +172,7 @@ void DTSegmentUpdator::fit(DTRecSegment4D* seg) const { seg->setCovMatrix(mat); seg->setCovMatrixForZed(posZInCh); - } - else if (seg->hasPhi()) { DTChamberRecSegment2D *segPhi=seg->phiSegment(); @@ -166,7 +187,6 @@ void DTSegmentUpdator::fit(DTRecSegment4D* seg) const { seg->setCovMatrix(mat); } - else if (seg->hasZed()) { DTSLRecSegment2D *segZed = seg->zSegment(); @@ -190,49 +210,91 @@ void DTSegmentUpdator::fit(DTRecSegment4D* seg) const { } } -bool DTSegmentUpdator::fit(DTSegmentCand* seg) const { + + +bool DTSegmentUpdator::fit(DTSegmentCand* seg, bool allow3par, const bool fitdebug) const { + +// if (debug && fitdebug) cout << "[DTSegmentUpdator] Fit DTRecSegment2D" << endl; if (!seg->good()) return false; +// DTSuperLayerId DTid = (DTSuperLayerId)seg->superLayer()->id(); +// if (DTid.superlayer()==2) +// allow3par = 0; + vector x; vector y; vector sigy; + vector lfit; + vector dist; + int i=0; + + x.reserve(8); + y.reserve(8); + sigy.reserve(8); + lfit.reserve(8); + dist.reserve(8); DTSegmentCand::AssPointCont hits=seg->hits(); - for (DTSegmentCand::AssPointCont::const_iterator iter=hits.begin(); - iter!=hits.end(); ++iter) { + for (DTSegmentCand::AssPointCont::const_iterator iter=hits.begin(); iter!=hits.end(); ++iter) { LocalPoint pos = (*iter).first->localPosition((*iter).second); + float xwire = (((*iter).first)->localPosition(DTEnums::Left).x() + ((*iter).first)->localPosition(DTEnums::Right).x()) /2.; + float distance = pos.x() - xwire; + + if ((*iter).second==DTEnums::Left) lfit.push_back(1); + else lfit.push_back(-1); + + dist.push_back(distance); + sigy.push_back(sqrt((*iter).first->localPositionError().xx())); x.push_back(pos.z()); y.push_back(pos.x()); - sigy.push_back(sqrt((*iter).first->localPositionError().xx())); + i++; } LocalPoint pos; LocalVector dir; AlgebraicSymMatrix covMat(2); - + float cminf =0.; + float vminf =0.; double chi2 = 0.; - fit(x,y,sigy,pos,dir,covMat,chi2); + double t0_corr = 0.; + + fit(x,y,lfit,dist,sigy,pos,dir,cminf,vminf,covMat,chi2,allow3par); + if (cminf!=0) t0_corr=-cminf/0.00543; // convert drift distance to time + + if (debug && fitdebug) cout << " DTcand chi2: " << chi2 << endl; + if (debug && fitdebug) cout << " DTcand t0: " << t0_corr << endl; seg->setPosition(pos); seg->setDirection(dir); - - //cout << "pos " << segPosition<< endl; - //cout << "dir " << segDirection<< endl; - + seg->sett0(t0_corr); seg->setCovMatrix(covMat); + seg->setChi2(chi2); + + // cout << "pos " << pos << endl; + // cout << "dir " << dir << endl; // cout << "Mat " << covMat << endl; - seg->setChi2(chi2); return true; } -void DTSegmentUpdator::fit(DTRecSegment2D* seg) const { - // WARNING: since this method is called both with a 2D and a 2DPhi as argument - // seg->geographicalId() can be a superLayerId or a chamberId +void DTSegmentUpdator::fit(DTRecSegment2D* seg, bool allow3par) const { + + if(debug) cout << "[DTSegmentUpdator] Fit DTRecSegment2D - 3par: " << allow3par << endl; vector x; vector y; vector sigy; + vector lfit; + vector dist; + x.reserve(8); + y.reserve(8); + sigy.reserve(8); + lfit.reserve(8); + dist.reserve(8); + +// DTSuperLayerId DTid = (DTSuperLayerId)seg->geographicalId(); +// if (DTid.superlayer()==2) +// allow3par = 0; vector hits=seg->specificRecHits(); for (vector::const_iterator hit=hits.begin(); @@ -241,18 +303,22 @@ void DTSegmentUpdator::fit(DTRecSegment2D* seg) const { // I have to get the hits position (the hit is in the layer rf) in SL frame... GlobalPoint glbPos = ( theGeom->layer( hit->wireId().layerId() ) )->toGlobal(hit->localPosition()); LocalPoint pos = ( theGeom->idToDet(seg->geographicalId()) )->toLocal(glbPos); - x.push_back(pos.z()); y.push_back(pos.x()); + const DTLayer* layer = theGeom->layer( hit->wireId().layerId() ); + float xwire = layer->specificTopology().wirePosition(hit->wireId().wire()); + float distance = fabs(hit->localPosition().x() - xwire); + dist.push_back(distance); + + int ilc = ( hit->lrSide() == DTEnums::Left ) ? 1 : -1; + lfit.push_back(ilc); + // Get local error in SL frame //RB: is it right in this way? ErrorFrameTransformer tran; - GlobalError glbErr = - tran.transform( hit->localPositionError(),(theGeom->layer( hit->wireId().layerId() ))->surface()); - LocalError slErr = - tran.transform( glbErr, (theGeom->idToDet(seg->geographicalId()))->surface()); - + GlobalError glbErr = tran.transform( hit->localPositionError(),(theGeom->layer( hit->wireId().layerId() ))->surface()); + LocalError slErr = tran.transform( glbErr, (theGeom->idToDet(seg->geographicalId()))->surface()); sigy.push_back(sqrt(slErr.xx())); } @@ -260,28 +326,38 @@ void DTSegmentUpdator::fit(DTRecSegment2D* seg) const { LocalVector dir; AlgebraicSymMatrix covMat(2); double chi2 = 0.; + float cminf=0.; + float vminf=0.; + double t0_corr = 0.; - fit(x,y,sigy,pos,dir,covMat,chi2); + fit(x,y,lfit,dist,sigy,pos,dir,cminf,vminf,covMat,chi2,allow3par); + if (cminf!=0) t0_corr=-cminf/0.00543; // convert drift distance to time + + if (debug) cout << " DTSeg2d chi2: " << chi2 << endl; + if (debug) cout << " DTSeg2d Fit t0: " << t0_corr << endl; + // cout << "pos " << segPosition << endl; + // cout << "dir " << segDirection << endl; + // cout << "Mat " << mat << endl; seg->setPosition(pos); seg->setDirection(dir); - - //cout << "pos " << segPosition << endl; - //cout << "dir " << segDirection << endl; - seg->setCovMatrix(covMat); - // cout << "Mat " << mat << endl; - seg->setChi2(chi2); + seg->setT0(t0_corr); } void DTSegmentUpdator::fit(const vector& x, const vector& y, + const vector& lfit, + const vector& dist, const vector& sigy, LocalPoint& pos, LocalVector& dir, + float& cminf, + float& vminf, AlgebraicSymMatrix& covMatrix, - double& chi2) const { + double& chi2, + const bool allow3par) const { float slope = 0.; float intercept = 0.; @@ -289,8 +365,26 @@ void DTSegmentUpdator::fit(const vector& x, float covii = 0.; float covsi = 0.; - // do the fit - theFitter->fit(x,y,x.size(),sigy,slope,intercept,covss,covii,covsi); + cminf=0; + vminf=0; + + int leftHits=0, rightHits=0; + for (unsigned int i=0; ifit(x,y,x.size(),sigy,slope,intercept,chi2,covss,covii,covsi); + + // If we have at least one left and one right hit we can try the 3 parameter fit (if it is switched on) + // FIXME: currently the covariance matrix from the 2-par fit is kept + if (leftHits && rightHits && (leftHits+rightHits>3)) { + theFitter->fitNpar(3,x,y,lfit,dist,sigy,slope,intercept,cminf,vminf,chi2,debug); + double t0_corr=-cminf/0.00543; + if (fabs(t0_corr)<20. && !allow3par) { + theFitter->fit(x,y,x.size(),sigy,slope,intercept,chi2,covss,covii,covsi); + cminf=0; + } + } + // cout << "slope " << slope << endl; // cout << "intercept " << intercept << endl; @@ -307,13 +401,6 @@ void DTSegmentUpdator::fit(const vector& x, covMatrix[0][0] = covss; // this is var(dy/dz) covMatrix[1][1] = covii; // this is var(y) covMatrix[1][0] = covsi; // this is cov(dy/dz,y) - - /* Calculate chi2. */ - chi2 = 0.; - for(unsigned int i=0; ilocalPosition().x(),+segPosAtLayer.y(),0.); - GlobalPoint glbpos= theGeom->layer( hit->wireId().layerId() )->toGlobal(hitPos); - newHit1D.setPosition(hitPos); - ok = theAlgo->compute(layer,*hit,angle,glbpos,newHit1D); } else if (step == 4) { @@ -373,9 +456,7 @@ void DTSegmentUpdator::updateHits(DTRecSegment2D* seg, GlobalPoint &gpos, const float xwire = layer->specificTopology().wirePosition(hit->wireId().wire()); const float distance = fabs(hit->localPosition().x() - xwire); - const int ilc = ( hit->lrSide() == DTEnums::Left ) ? 1 : -1; - const double dy_corr = (vminf*ilc*distance-cminf*ilc ); LocalPoint point(hit->localPosition().x() + dy_corr, +segPosAtLayer.y(), 0.); @@ -383,9 +464,7 @@ void DTSegmentUpdator::updateHits(DTRecSegment2D* seg, GlobalPoint &gpos, //double final_hit_resol = T0_hit_resolution; //if(newHit1D.wireId().layerId().superlayerId().superLayer() != 2) final_hit_resol = final_hit_resol * 0.8; //LocalError error(final_hit_resol * final_hit_resol,0.,0.); - LocalError error(T0_hit_resolution*T0_hit_resolution,0.,0.); - newHit1D.setPositionAndError(point, error); //FIXME: check that the hit is still inside the cell @@ -472,8 +551,7 @@ void DTSegmentUpdator::rejectBadHits(DTChamberRecSegment2D* phiSeg) const { } if(debug) cout << " Residuals computed! "<< endl; - - + // Perform bad hit rejecting -- update hits vector updatedRecHits; @@ -529,8 +607,6 @@ void DTSegmentUpdator::rejectBadHits(DTChamberRecSegment2D* phiSeg) const { y_upd.push_back(pos.x()); cout << " x_upd: "<< pos.z() << " y_upd: "<< pos.x() << endl; - - } cout << " end of segment! " << endl; @@ -549,6 +625,7 @@ void DTSegmentUpdator::calculateT0corr(DTRecSegment4D* seg) const { void DTSegmentUpdator::calculateT0corr(DTRecSegment2D* seg) const { // WARNING: since this method is called both with a 2D and a 2DPhi as argument // seg->geographicalId() can be a superLayerId or a chamberId + if(debug) cout << "[DTSegmentUpdator] CalculateT0corr DTRecSegment4D" << endl; vector d_drift; vector x; @@ -584,11 +661,12 @@ void DTSegmentUpdator::calculateT0corr(DTRecSegment2D* seg) const { double chi2fit = 0.; float cminf = 0.; - double vminf = 0.; + float vminf = 0.; + float a,b; if ( nptfit > 2 ) { //NB chi2fit is normalized - Fit4Var(x,y,lc,d_drift,nptfit,cminf,vminf,chi2fit); + theFitter->fit4Var(x,y,lc,d_drift,nptfit,a,b,cminf,vminf,chi2fit,vdrift_4parfit,debug); double t0cor = -999.; if(cminf > -998.) t0cor = - cminf/0.00543 ; // in ns @@ -602,246 +680,3 @@ void DTSegmentUpdator::calculateT0corr(DTRecSegment2D* seg) const { seg->setVdrift(vminf); // vdrift correction are recorded in the segment } } - - -void DTSegmentUpdator::Fit4Var(const vector& xfit, - const vector& yfit, - const vector& lfit, - const vector& tfit, - const int nptfit, - float& cminf, - double& vminf, - double& chi2fit) const { - - const double sigma = 0.0295;// errors can be inserted .just load them/that is the usual TB resolution value for DT chambers - double aminf = 0.; - double bminf = 0.; - int nppar = 0; - double sx = 0.; - double sx2 = 0.; - double sy = 0.; - double sxy = 0.; - double sl = 0.; - double sl2 = 0.; - double sly = 0.; - double slx = 0.; - double st = 0.; - double st2 = 0.; - double slt = 0.; - double sltx = 0.; - double slty = 0.; - double chi2fitN2 = -1. ; - double chi2fit3 = -1.; - double chi2fitN3 = -1. ; - double chi2fitN4 = -1.; - float bminf3 = bminf; - float aminf3 = aminf; - float cminf3 = cminf; - int nppar2 = 0; - int nppar3 = 0; - int nppar4 = 0; - - cminf = -999.; - vminf = 0.; - - for (int j=0; j=5) { - const double det = (a1*a1*(b2*v6 - b6*b6) - a1*(a2*a2*v6 - 2*a2*a6*b6 + a6*a6*b2 + b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)) - + a2*a2*c6*c6 + 2*a2*(a3*(b3*v6 - b6*c6) - a6*b3*c6) + a3*a3*(b6*b6 - b2*v6) - + a6*(2*a3*(b2*c6 - b3*b6) + a6*b3*b3)); - - // the dv/vdrift correction may be computed under vdrift_4parfit request; - if (det != 0) { - nppar = 4; - chi2fit = 0.; - // computation of a, b, c e v - aminf = (a1*(a2*(b6*d6 - v6*d2) + a6*(b6*d2 - b2*d6) + d1*(b2*v6 - b6*b6)) - a2*(b3*(c6*d6 - v6*d3) - + c6*(b6*d3 - c6*d2)) + a3*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) - + a6*(b2*c6*d3 + b3*(b3*d6 - b6*d3 - c6*d2)) - d1*(b2*c6*c6 + b3*(b3*v6 - 2*b6*c6)))/det; - bminf = - (a1*a1*(b6*d6 - v6*d2) - a1*(a2*(a6*d6 - v6*d1) - a6*a6*d2 + a6*b6*d1 + b3*(c6*d6 - v6*d3) - + c6*(b6*d3 - c6*d2)) + a2*(a3*(c6*d6 - v6*d3) + c6*(a6*d3 - c6*d1)) + a3*a3*(v6*d2 - b6*d6) - + a3*(a6*(b3*d6 + b6*d3 - 2*c6*d2) - d1*(b3*v6 - b6*c6)) - a6*b3*(a6*d3 - c6*d1))/det; - cminf = -(a1*(b2*(c6*d6 - v6*d3) + b3*(v6*d2 - b6*d6) + b6*(b6*d3 - c6*d2)) + a2*a2*(v6*d3 - c6*d6) - + a2*(a3*(b6*d6 - v6*d2) + a6*(b3*d6 - 2*b6*d3 + c6*d2) - d1*(b3*v6 - b6*c6)) - + a3*(d1*(b2*v6 - b6*b6) - a6*(b2*d6 - b6*d2)) + a6*(a6*(b2*d3 - b3*d2) - d1*(b2*c6 - b3*b6)))/det; - vminf = - (a1*a1*(b2*d6 - b6*d2) - a1*(a2*a2*d6 - a2*(a6*d2 + b6*d1) + a6*b2*d1 + b2*c6*d3 - + b3*(b3*d6 - b6*d3 - c6*d2)) + a2*a2*c6*d3 + a2*(a3*(2*b3*d6 - b6*d3 - c6*d2) - b3*(a6*d3 + c6*d1)) - + a3*a3*(b6*d2 - b2*d6) + a3*(a6*(b2*d3 - b3*d2) + d1*(b2*c6 - b3*b6)) + a6*b3*b3*d1)/det; - - // chi 2 - for (int j=0; j= 0.29) { - // for safety and for code construction..dont accept correction on dv/vdrift greater then 0.09 - vminf = 0.; - cminf = cminf3; - aminf = aminf3; - bminf = bminf3; - nppar = 3; - chi2fit = chi2fit3; - } - - } //end if vdrift - - if(!vdrift_4parfit){ //if not required explicitly leave the t0 and track step as at step 3 - // just update vdrift value vmin for storing in the segments for monitoring - cminf = cminf3; - aminf = aminf3; - bminf = bminf3; - nppar = 3; - chi2fit = chi2fit3; - } - - nppar4 = nppar; - - } //end nptfit >=3 - - if (debug) { - cout << " dt0= 0 : slope 4 = " << bminf << " pos out = " << aminf <<" chi2fitN4 = " << chi2fitN4 - << " nppar= " << nppar4 << " T0_ev ns= " << cminf/0.00543 <<" delta v = " << vminf <& x, const std::vector& y, + const std::vector& lfit, + const std::vector& dist, const std::vector& sigy, LocalPoint& pos, LocalVector& dir, + float& cminf, + float& vminf, AlgebraicSymMatrix& covMat, - double& chi2) const; - - void Fit4Var(const std::vector& xfit, - const std::vector& yfit, - const std::vector& lfit, - const std::vector& tfit, - const int nptfit, - float& cminf, double& vminf, - double& chi2fit) const; + double& chi2, + const bool allow3par = false) const; bool vdrift_4parfit; double T0_hit_resolution; diff --git a/trunk/RecoMuon/MuonIdentification/src/DTTimingExtractor.cc b/trunk/RecoMuon/MuonIdentification/src/DTTimingExtractor.cc index 43987e9c94627..a5c651d16cf85 100644 --- a/trunk/RecoMuon/MuonIdentification/src/DTTimingExtractor.cc +++ b/trunk/RecoMuon/MuonIdentification/src/DTTimingExtractor.cc @@ -141,6 +141,10 @@ DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRe // get the DT segments that were used to construct the muon std::vector range = theMatcher->matchDT(*muonTrack,iEvent); + + if (debug) + std::cout << " The muon track matches " << range.size() << " segments." << std::endl; + // create a collection on TimeMeasurements for the track for (std::vector::iterator rechit = range.begin(); rechit!=range.end();++rechit) { @@ -149,6 +153,7 @@ DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRe DetId id = (*rechit)->geographicalId(); DTChamberId chamberId(id.rawId()); int station = chamberId.station(); + if (debug) std::cout << "Matched DT segment in station " << station << std::endl; // use only segments with both phi and theta projections present (optional) bool bothProjections = ( ((*rechit)->hasPhi()) && ((*rechit)->hasZed()) ); @@ -218,6 +223,8 @@ DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRe thisHit.posInLayer += slCorr; } + if (debug) std::cout << " dist: " << dist << " t0: " << thisHit.posInLayer << std::endl; + tms.push_back(thisHit); } } // phi = (0,1) @@ -284,7 +291,7 @@ DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRe std::cout << " t0 = zero, Left hits: " << hitxl.size() << " Right hits: " << hitxr.size() << std::endl; continue; } - + // a segment must have at least one left and one right hit if ((!hitxl.size()) || (!hityl.size())) continue; @@ -301,6 +308,8 @@ DTTimingExtractor::fillTiming(TimeMeasurementSequence &tmSequence, reco::TrackRe double hitLocalPos = tm->posInLayer; int hitSide = -tm->isLeft*2+1; double t0_segm = (-(hitSide*segmLocalPos)+(hitSide*hitLocalPos))/0.00543+tm->timeCorr; + + if (debug) std::cout << " Segm hit. dstnc: " << tm->distIP << " t0: " << t0_segm << std::endl; dstnc.push_back(tm->distIP); dsegm.push_back(t0_segm); diff --git a/trunk/RecoMuon/TrackingTools/src/MuonSegmentMatcher.cc b/trunk/RecoMuon/TrackingTools/src/MuonSegmentMatcher.cc index 4139c046e3e89..e9f3852aa0268 100644 --- a/trunk/RecoMuon/TrackingTools/src/MuonSegmentMatcher.cc +++ b/trunk/RecoMuon/TrackingTools/src/MuonSegmentMatcher.cc @@ -76,14 +76,16 @@ vector MuonSegmentMatcher::matchDT(const reco::Track &muo // Loop and select DT recHits for(trackingRecHit_iterator hit = muon.recHitsBegin(); hit != muon.recHitsEnd(); ++hit) { - if ( !(*hit)->isValid()) continue; - if ( (*hit)->geographicalId().det() != DetId::Muon ) continue; - if ( (*hit)->geographicalId().subdetId() != MuonSubdetId::DT ) continue; if (!(*hit)->isValid()) continue; - if ((*hit)->recHits().size()>1) segments = true; + if ( (*hit)->geographicalId().det() != DetId::Muon ) continue; + if ( (*hit)->geographicalId().subdetId() != MuonSubdetId::DT ) continue; + if ((*hit)->recHits().size()) + if ((*(*hit)->recHits().begin())->recHits().size()>1) segments = true; dtHits.push_back(*hit); } + // cout << "Muon DT hits found: " << dtHits.size() << " segments " << segments << endl; + double PhiCutParameter=dtRadius_; double ZCutParameter=dtRadius_; double matchRatioZ=0; @@ -91,13 +93,11 @@ vector MuonSegmentMatcher::matchDT(const reco::Track &muo for (DTRecSegment4DCollection::const_iterator rechit = dtRecHits->begin(); rechit!=dtRecHits->end();++rechit) { - if ( !rechit->isValid()) continue; LocalPoint pointLocal = rechit->localPosition(); if (segments) { // Loop over muon recHits for(trackingRecHit_iterator hit = dtHits.begin(); hit != dtHits.end(); ++hit) { - if ( !(*hit)->isValid()) continue; // Pick the one in the same DT Chamber as the muon DetId idT = (*hit)->geographicalId(); diff --git a/trunk/Validation/DTRecHits/interface/DTHitQualityUtils.h b/trunk/Validation/DTRecHits/interface/DTHitQualityUtils.h index a4411dc447b64..966b7d127b9bd 100644 --- a/trunk/Validation/DTRecHits/interface/DTHitQualityUtils.h +++ b/trunk/Validation/DTRecHits/interface/DTHitQualityUtils.h @@ -43,11 +43,15 @@ class DTHitQualityUtils { static std::pair findMuSimSegmentDirAndPos(const std::pair& inAndOutSimHit, const DetId detId, const DTGeometry *muonGeom); /// Find the angles from a segment direction: + /// atan(dx/dz) = "phi" angle in the chamber RF + /// atan(dy/dz) = "theta" angle in the chamber RF (note: this has opposite sign in the SLZ RF!) static std::pair findSegmentAlphaAndBeta(const LocalVector& direction); // Set the verbosity level static bool debug; + //Find angle error + static double sigmaAngle(double Angle, double sigma2TanAngle); protected: diff --git a/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.cc b/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.cc index 191abadac0785..62ea09b1bc7c3 100644 --- a/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.cc +++ b/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.cc @@ -30,8 +30,15 @@ using namespace std; using namespace edm; - - +// In phi SLs, The dependency on X and angle is specular in positive +// and negative wheels. Since positive and negative wheels are filled +// together into the same plots, it is useful to mirror negative wheels +// so that the actual dependency can be observerd instead of an artificially +// simmetrized one. +// Set mirrorMinusWheels to avoid this. +namespace { + bool mirrorMinusWheels = true; +} // Constructor DTRecHitQuality::DTRecHitQuality(const ParameterSet& pset){ @@ -52,11 +59,9 @@ DTRecHitQuality::DTRecHitQuality(const ParameterSet& pset){ doStep3 = pset.getUntrackedParameter("doStep3", false); doall = pset.getUntrackedParameter("doall", false); local = pset.getUntrackedParameter("local", true); - // if(doall) doStep1 - // Create the root file - //theFile = new TFile(rootFileName.c_str(), "RECREATE"); - //theFile->cd(); +} +void DTRecHitQuality::beginRun(const edm::Run& iRun, const edm::EventSetup &setup) { // ---------------------- // get hold of back-end interface @@ -350,7 +355,7 @@ float DTRecHitQuality::simHitDistFromWire(const DTLayer* layer, return fabs(xEntry - (entryP.z()*(xExit-xEntry))/(exitP.z()-entryP.z()));//FIXME: check... } -// Compute SimHit impact angle (in direction perp to wire) +// Compute SimHit impact angle (in direction perp to wire), in the SL RF float DTRecHitQuality::simHitImpactAngle(const DTLayer* layer, DTWireId wireId, const PSimHit& hit) { @@ -367,7 +372,10 @@ float DTRecHitQuality::simHitDistFromFE(const DTLayer* layer, LocalPoint entryP = hit.entryPoint(); LocalPoint exitP = hit.exitPoint(); float wireLenght=layer->specificTopology().cellLenght(); - return (entryP.y()+exitP.y())/2.+wireLenght; + // FIXME: should take only wireLenght/2.; + // moreover, pos+cellLenght/2. is shorter than the distance from FE. + // In fact it would make more sense to make plots vs y. + return (entryP.y()+exitP.y())/2.+wireLenght; } @@ -422,6 +430,9 @@ void DTRecHitQuality::compute(const DTGeometry *dtGeom, wireAndSHits != simHitsPerWire.end(); wireAndSHits++) { DTWireId wireId = (*wireAndSHits).first; + int wheel = wireId.wheel(); + int sl = wireId.superLayer(); + vector simHitsInCell = (*wireAndSHits).second; // Get the layer @@ -472,77 +483,83 @@ void DTRecHitQuality::compute(const DTGeometry *dtGeom, float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer); if(debug) cout << " SimHit distance from wire: " << simHitWireDist << endl - << " SimHit distance from FE: " << simHitFEDist << endl - << " SimHit distance angle " << simHitTheta << endl - << " RecHit distance from wire: " << recHitWireDist << endl; + << " SimHit distance from FE: " << simHitFEDist << endl + << " SimHit angle in layer RF: " << simHitTheta << endl + << " RecHit distance from wire: " << recHitWireDist << endl; float recHitErr = recHitPositionError(*theBestRecHit); HRes1DHit *hRes = 0; HRes1DHit *hResTot = 0; + // Mirror angle in phi so that + and - wheels can be plotted together + if (mirrorMinusWheels && wheel<0 && sl!=2){ + simHitTheta *= -1.; + // Note: local X, if used, would have to be mirrored as well + } + // Fill residuals and pulls // Select the histo to be filled if(step == 1) { // Step 1 - if(wireId.superLayer() != 2) { + if(sl != 2) { hResTot = hRes_S1RPhi; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S1RPhi_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S1RPhi_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S1RPhi_W2; } else { hResTot = hRes_S1RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S1RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S1RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S1RZ_W2; } } else if(step == 2) { // Step 2 - if(wireId.superlayer() != 2) { + if(sl != 2) { hRes = hRes_S2RPhi; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S2RPhi_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S2RPhi_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S2RPhi_W2; } else { hResTot = hRes_S2RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S2RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S2RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S2RZ_W2; } } else if(step == 3) { // Step 3 - if(wireId.superlayer() != 2) { + if(sl != 2) { hResTot = hRes_S3RPhi; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S3RPhi_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S3RPhi_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S3RPhi_W2; - if (local) hRes_S3RPhiWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitTheta, simHitFEDist, recHitWireDist, simHitGlobalPos.eta(),simHitGlobalPos.phi(),recHitErr,wireId.station()); + if (local) hRes_S3RPhiWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitTheta, simHitFEDist, recHitWireDist, simHitGlobalPos.eta(),simHitGlobalPos.phi(),recHitErr,wireId.station()); } else { hResTot = hRes_S3RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hRes = hRes_S3RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hRes = hRes_S3RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hRes = hRes_S3RZ_W2; - if (local) hRes_S3RZWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitTheta, simHitFEDist, recHitWireDist, simHitGlobalPos.eta(),simHitGlobalPos.phi(),recHitErr,wireId.station()); + if (local) hRes_S3RZWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitTheta, simHitFEDist, recHitWireDist, simHitGlobalPos.eta(),simHitGlobalPos.phi(),recHitErr,wireId.station()); } } // Fill @@ -559,48 +576,48 @@ void DTRecHitQuality::compute(const DTGeometry *dtGeom, HEff1DHit *hEffTot = 0; if(step == 1) { // Step 1 - if(wireId.superlayer() != 2) { + if(sl != 2) { hEff = hEff_S1RPhi; - if (local) hEff_S1RPhiWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); + if (local) hEff_S1RPhiWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); } else { hEffTot = hEff_S1RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hEff = hEff_S1RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hEff = hEff_S1RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hEff = hEff_S1RZ_W2; - if (local) hEff_S1RZWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); + if (local) hEff_S1RZWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); } } else if(step == 2) { // Step 2 - if(wireId.superlayer() != 2) { + if(sl != 2) { hEff = hEff_S2RPhi; } else { hEffTot = hEff_S2RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hEff = hEff_S2RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hEff = hEff_S2RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hEff = hEff_S2RZ_W2; } } else if(step == 3) { // Step 3 - if(wireId.superlayer() != 2) { + if(sl != 2) { hEff = hEff_S3RPhi; - if (local) hEff_S3RPhiWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); + if (local) hEff_S3RPhiWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); } else { hEffTot = hEff_S3RZ; - if(wireId.wheel() == 0) + if(wheel == 0) hEff = hEff_S3RZ_W0; - if(abs(wireId.wheel()) == 1) + if(abs(wheel) == 1) hEff = hEff_S3RZ_W1; - if(abs(wireId.wheel()) == 2) + if(abs(wheel) == 2) hEff = hEff_S3RZ_W2; - if (local) hEff_S3RZWS[abs(wireId.wheel())][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); + if (local) hEff_S3RZWS[abs(wheel)][wireId.station()-1]->Fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed); } } diff --git a/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.h b/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.h index 8098996ac9990..d208db238947e 100644 --- a/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.h +++ b/trunk/Validation/DTRecHits/plugins/DTRecHitQuality.h @@ -60,8 +60,10 @@ class DTRecHitQuality : public edm::EDAnalyzer { void analyze(const edm::Event & event, const edm::EventSetup& eventSetup); // Write the histos to file + virtual void beginRun(const edm::Run& iRun, const edm::EventSetup &setup); + void endJob(); -void endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, + void endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, edm::EventSetup const& c); protected: diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.cc b/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.cc index 985e58a373cf4..48b70340f27ea 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.cc +++ b/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.cc @@ -48,13 +48,12 @@ DTSegment2DQuality::DTSegment2DQuality(const ParameterSet& pset) { //sigma resolution on angle sigmaResAngle = pset.getParameter("sigmaResAngle"); - // Create the root file - //theFile = new TFile(rootFileName.c_str(), "RECREATE"); - //theFile->cd(); - if(debug) cout << "[DTSegment2DQuality] Constructor called " << endl; - +} + +void DTSegment2DQuality::beginRun(const edm::Run& iRun, const edm::EventSetup &setup) { + // ---------------------- // get hold of back-end interface dbe_ = 0; diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.h b/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.h index c2aa5ea256927..df36d84bc8d14 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.h +++ b/trunk/Validation/DTRecHits/plugins/DTSegment2DQuality.h @@ -39,6 +39,9 @@ class DTSegment2DQuality : public edm::EDAnalyzer { /// Perform the real analysis void analyze(const edm::Event & event, const edm::EventSetup& eventSetup); + + virtual void beginRun(const edm::Run& iRun, const edm::EventSetup &setup); + // Write the histos to file void endJob(); diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.cc b/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.cc index 07c433d0ac441..fb930bdd20449 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.cc +++ b/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.cc @@ -53,11 +53,11 @@ DTSegment2DSLPhiQuality::DTSegment2DSLPhiQuality(const ParameterSet& pset) { sigmaResAngle = pset.getParameter("sigmaResAngle"); doall = pset.getUntrackedParameter("doall", false); local = pset.getUntrackedParameter("local", false); +} + + +void DTSegment2DSLPhiQuality::beginRun(const edm::Run& iRun, const edm::EventSetup &setup) { - // Create the root file - //theFile = new TFile(rootFileName.c_str(), "RECREATE"); - //theFile->cd(); - // ---------------------- // get hold of back-end interface dbe_ = 0; dbe_ = Service().operator->(); diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.h b/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.h index fe457f0715941..dc1f6c09d22bc 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.h +++ b/trunk/Validation/DTRecHits/plugins/DTSegment2DSLPhiQuality.h @@ -39,6 +39,9 @@ class DTSegment2DSLPhiQuality : public edm::EDAnalyzer { /// Perform the real analysis void analyze(const edm::Event & event, const edm::EventSetup& eventSetup); + + virtual void beginRun(const edm::Run& iRun, const edm::EventSetup &setup); + // Write the histos to file void endJob(); void endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.cc b/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.cc index 5e6e092f9bfcc..646aeb6aca871 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.cc +++ b/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.cc @@ -33,7 +33,15 @@ using namespace std; using namespace edm; -// Book the histos +// In phi SLs, The dependency on X and angle is specular in positive +// and negative wheels. Since positive and negative wheels are filled +// together into the same plots, it is useful to mirror negative wheels +// so that the actual dependency can be observerd instead of an artificially +// simmetrized one. +// Set mirrorMinusWheels to avoid this. +namespace { + bool mirrorMinusWheels = true; +} // Constructor DTSegment4DQuality::DTSegment4DQuality(const ParameterSet& pset) { @@ -55,12 +63,10 @@ DTSegment4DQuality::DTSegment4DQuality(const ParameterSet& pset) { sigmaResBeta = pset.getParameter("sigmaResBeta"); doall = pset.getUntrackedParameter("doall", false); local = pset.getUntrackedParameter("local", false); +} +void DTSegment4DQuality::beginRun(const edm::Run& iRun, const edm::EventSetup &setup) { - // Create the root file - //theFile = new TFile(rootFileName.c_str(), "RECREATE"); - //theFile->cd(); -// ---------------------- // get hold of back-end interface dbe_ = 0; dbe_ = Service().operator->(); @@ -96,13 +102,11 @@ DTSegment4DQuality::DTSegment4DQuality(const ParameterSet& pset) { TString nameWS=(name+w+"_St"+s); h4DHitWS[w][s-1] = new HRes4DHit(nameWS.Data(),dbe_,doall,local); hEffWS[w][s-1] = new HEff4DHit (nameWS.Data(),dbe_); - dbe_->setCurrentFolder("DT/4DSegments/"); - hHitMult[w][s-1] = dbe_->book2D("4D_" +nameWS+ "_hNHits", "NHits",12,0,12, 6,0,6); - ht0[w][s-1] = dbe_->book2D("4D_" +nameWS+ "_ht0", "t0",200,-25,25,200,-25,25); } } } -} +}; + // Destructor DTSegment4DQuality::~DTSegment4DQuality(){ @@ -110,8 +114,6 @@ DTSegment4DQuality::~DTSegment4DQuality(){ } void DTSegment4DQuality::endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, edm::EventSetup const& c){ - - } void DTSegment4DQuality::endJob() { @@ -142,6 +144,8 @@ void DTSegment4DQuality::endJob() { void DTSegment4DQuality::analyze(const Event & event, const EventSetup& eventSetup){ //theFile->cd(); + const float epsilon=5e-5; // numerical accuracy on angles [rad} + // Get the DT Geometry ESHandle dtGeom; eventSetup.get().get(dtGeom); @@ -154,10 +158,14 @@ void DTSegment4DQuality::endJob() { map simHitsPerCh; for(PSimHitContainer::const_iterator simHit = simHits->begin(); simHit != simHits->end(); simHit++){ - // Create the id of the chamber (the simHits in the DT known their wireId) - DTChamberId chamberId = (((DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId(); - // Fill the map - simHitsPerCh[chamberId].push_back(*simHit); + + // Consider only muon simhits; the others are not considered elsewhere in this class! + if (abs((*simHit).particleType())==13) { + // Create the id of the chamber (the simHits in the DT known their wireId) + DTChamberId chamberId = (((DTWireId(simHit->detUnitId())).layerId()).superlayerId()).chamberId(); + // Fill the map + simHitsPerCh[chamberId].push_back(*simHit); + } } // Get the 4D rechits from the event @@ -170,41 +178,41 @@ void DTSegment4DQuality::endJob() { return; } - // Loop over all chambers containing a segment - DTRecSegment4DCollection::id_iterator chamberId; - for (chamberId = segment4Ds->id_begin(); - chamberId != segment4Ds->id_end(); - ++chamberId){ + // Loop over all chambers containing a (muon) simhit + for (map::const_iterator simHitsInChamber = simHitsPerCh.begin(); simHitsInChamber != simHitsPerCh.end(); ++simHitsInChamber) { - if((*chamberId).station() == 4) - continue; //use DTSegment2DSLPhiQuality to analyze MB4 performaces + DTChamberId chamberId = simHitsInChamber->first; + int station = chamberId.station(); + if (station == 4 && !(local) ) continue; //use DTSegment2DSLPhiQuality to analyze MB4 performaces in DQM + int wheel = chamberId.wheel(); //------------------------- simHits ---------------------------// - //Get simHits of each chamber - PSimHitContainer simHits = simHitsPerCh[(*chamberId)]; + //Get simHits of this chamber + const PSimHitContainer& simHits = simHitsInChamber->second; // Map simhits per wire map simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(simHits); map muSimHitPerWire = DTHitQualityUtils::mapMuSimHitsPerWire(simHitsPerWire); int nMuSimHit = muSimHitPerWire.size(); - if(nMuSimHit == 0 || nMuSimHit == 1) { - if(debug && nMuSimHit == 1) - cout << "[DTSegment4DQuality] Only " << nMuSimHit << " mu SimHit in this chamber, skipping!" << endl; - continue; // If no or only one mu SimHit is found skip this chamber + if(nMuSimHit <2) { // Skip chamber with less than 2 cells with mu hits + continue; } if(debug) - cout << "=== Chamber " << (*chamberId) << " has " << nMuSimHit << " SimHits" << endl; + cout << "=== Chamber " << chamberId << " has " << nMuSimHit << " SimHits" << endl; //Find outer and inner mu SimHit to build a segment - pair inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire); + pair inAndOutSimHit = DTHitQualityUtils::findMuSimSegment(muSimHitPerWire); + + // Consider only sim segments crossing at least 2 SLs + if ((DTWireId(inAndOutSimHit.first->detUnitId())).superlayer() == (DTWireId(inAndOutSimHit.second->detUnitId())).superLayer()) continue; //Find direction and position of the sim Segment in Chamber RF pair dirAndPosSimSegm = DTHitQualityUtils::findMuSimSegmentDirAndPos(inAndOutSimHit, - (*chamberId),&(*dtGeom)); + chamberId,&(*dtGeom)); LocalVector simSegmLocalDir = dirAndPosSimSegm.first; LocalPoint simSegmLocalPos = dirAndPosSimSegm.second; - const DTChamber* chamber = dtGeom->chamber(*chamberId); + const DTChamber* chamber = dtGeom->chamber(chamberId); GlobalPoint simSegmGlobalPos = chamber->toGlobal(simSegmLocalPos); GlobalVector simSegmGlobalDir = chamber->toGlobal(simSegmLocalDir); @@ -217,6 +225,8 @@ void DTSegment4DQuality::endJob() { //Position (in eta,phi coordinates) in lobal RF float etaSimSeg = simSegmGlobalPos.eta(); float phiSimSeg = simSegmGlobalPos.phi(); + + double count_seg = 0; if(debug) cout<<" Simulated segment: local direction "<get(*chamberId); + DTRecSegment4DCollection::range range = segment4Ds->get(chamberId); int nsegm = distance(range.first, range.second); if(debug) - cout << " Chamber: " << *chamberId << " has " << nsegm + cout << " Chamber: " << chamberId << " has " << nsegm << " 4D segments" << endl; if (nsegm!=0) { @@ -240,63 +250,62 @@ void DTSegment4DQuality::endJob() { // the residual distribution (we are looking for residuals of segments // usefull to the track fit) for efficency purpose const DTRecSegment4D* bestRecHit = 0; - bool bestRecHitFound = false; double deltaAlpha = 99999; + double deltaBeta = 99999; // Loop over the recHits of this chamberId for (DTRecSegment4DCollection::const_iterator segment4D = range.first; segment4D!=range.second; ++segment4D){ - if (local) { - const DTChamberRecSegment2D* phiSeg = (*segment4D).phiSegment(); - const DTSLRecSegment2D* zSeg = (*segment4D).zSegment(); - - float t0phi = -999; - float t0z = -999; - int nHitPhi=0; - int nHitZ=0; - if (phiSeg) { - t0phi = phiSeg->t0(); - nHitPhi = phiSeg->recHits().size(); - } - - if (zSeg) { - t0z = zSeg->t0(); - nHitZ = zSeg->recHits().size(); - } - - hHitMult[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(nHitPhi,nHitZ); - ht0[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(t0phi,t0z); - - // Uncomment to skip segments w/o t0 computed - // if (!(t0phi>-998&&t0z>-998)) continue - } - - // Check the dimension - if((*segment4D).dimension() != 4) { - if(debug)cout << "[DTSegment4DQuality]***Error: This is not 4D segment!!!" << endl; + // Consider only segments with both projections + if(station!=4 && (*segment4D).dimension() != 4) { continue; } // Segment Local Direction and position (in Chamber RF) LocalVector recSegDirection = (*segment4D).localDirection(); - //LocalPoint recSegPosition = (*segment4D).localPosition(); + LocalPoint recSegPosition = (*segment4D).localPosition(); + + pair ab = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection); + float recSegAlpha = ab.first; + float recSegBeta = ab.second; - float recSegAlpha = DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).first; if(debug) - cout << " RecSegment direction: " << recSegDirection << endl - << " position : " << (*segment4D).localPosition() << endl - << " alpha : " << recSegAlpha << endl - << " beta : " << DTHitQualityUtils::findSegmentAlphaAndBeta(recSegDirection).second << endl; - - if(fabs(recSegAlpha - alphaSimSeg) < deltaAlpha) { - deltaAlpha = fabs(recSegAlpha - alphaSimSeg); - bestRecHit = &(*segment4D); - bestRecHitFound = true; - } + cout << &(*segment4D) + << " RecSegment direction: " << recSegDirection << endl + << " position : " << (*segment4D).localPosition() << endl + << " alpha : " << recSegAlpha << endl + << " beta : " << recSegBeta << endl + << " nhits : " << (*segment4D).phiSegment()->recHits().size() << " " << (((*segment4D).zSegment()!=0)?(*segment4D).zSegment()->recHits().size():0) + << endl; + + + float dAlphaRecSim=fabs(recSegAlpha - alphaSimSeg); + float dBetaRecSim =fabs(recSegBeta - betaSimSeg); + + + if ((fabs(recSegPosition.x()-simSegmLocalPos.x())<4) // require rec and sim segments to be ~in the same cell in x + && ((fabs(recSegPosition.y()-simSegmLocalPos.y())<4)||(*segment4D).dimension()<4)) { // ~in the same cell in y, if segment has the theta view + ++count_seg; + + if (fabs(dAlphaRecSim-deltaAlpha) < epsilon) { // Numerically equivalent alphas, choose based on beta + if (dBetaRecSim < deltaBeta) { + deltaAlpha = dAlphaRecSim; + deltaBeta = dBetaRecSim; + bestRecHit = &(*segment4D); + } + + } else if (dAlphaRecSim < deltaAlpha) { + deltaAlpha = dAlphaRecSim; + deltaBeta = dBetaRecSim; + bestRecHit = &(*segment4D); + } + } + } // End of Loop over all 4D RecHits - if(bestRecHitFound) { + if(bestRecHit) { + if (debug) cout << endl << "Chosen: " << bestRecHit << endl; // Best rechit direction and position in Chamber RF LocalPoint bestRecHitLocalPos = bestRecHit->localPosition(); LocalVector bestRecHitLocalDir = bestRecHit->localDirection(); @@ -304,8 +313,9 @@ void DTSegment4DQuality::endJob() { LocalError bestRecHitLocalPosErr = bestRecHit->localPositionError(); LocalError bestRecHitLocalDirErr = bestRecHit->localDirectionError(); - float alphaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).first; - float betaBestRHit = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir).second; + pair ab = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDir); + float alphaBestRHit = ab.first; + float betaBestRHit = ab.second; // Errors on alpha and beta // Get position and direction using the rx projection (so in SL @@ -313,94 +323,136 @@ void DTSegment4DQuality::endJob() { // frame //if (bestRecHit->hasZed() ) { const DTSLRecSegment2D * zedRecSeg = bestRecHit->zSegment(); - LocalPoint bestRecHitLocalPosRZ = zedRecSeg->localPosition(); - LocalVector bestRecHitLocalDirRZ = zedRecSeg->localDirection(); - // Errors on x and y - LocalError bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError(); - LocalError bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError(); - - float alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first; - - // Get SimSeg position and Direction in rZ SL frame - const DTSuperLayer* sl = dtGeom->superLayer(zedRecSeg->superLayerId()); - LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos); - LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir); - LocalPoint simSegLocalPosRZ = - simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.z()/(cos(simSegLocalDirRZ.theta()))); - float alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first; - - if (debug) cout << - "RZ SL: recPos " << bestRecHitLocalPosRZ << - "recDir " << bestRecHitLocalDirRZ << - "recAlpha " << alphaBestRHitRZ << endl << - "RZ SL: simPos " << simSegLocalPosRZ << - "simDir " << simSegLocalDirRZ << - "simAlpha " << alphaSimSegRZ << endl ; - //} - - - if(fabs(alphaBestRHit - alphaSimSeg) < 5*sigmaResAlpha && - fabs(betaBestRHit - betaSimSeg) < 5*sigmaResBeta && - fabs(bestRecHitLocalPos.x() - xSimSeg) < 5*sigmaResX && - fabs(bestRecHitLocalPos.y() - ySimSeg) < 5*sigmaResY) { - recHitFound = true; - } + LocalPoint bestRecHitLocalPosRZ; + LocalVector bestRecHitLocalDirRZ; + LocalError bestRecHitLocalPosErrRZ; + LocalError bestRecHitLocalDirErrRZ; + LocalPoint simSegLocalPosRZ; // FIXME: this is not set for segments with only the phi view. + float alphaBestRHitRZ = 0; // angle measured in the RZ SL, in its own frame + float alphaSimSegRZ = betaSimSeg; + if (zedRecSeg) { + bestRecHitLocalPosRZ = zedRecSeg->localPosition(); + bestRecHitLocalDirRZ = zedRecSeg->localDirection(); + // Errors on x and y + bestRecHitLocalPosErrRZ = zedRecSeg->localPositionError(); + bestRecHitLocalDirErrRZ = zedRecSeg->localDirectionError(); + + // angle measured in the RZ SL, in its own frame + alphaBestRHitRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(bestRecHitLocalDirRZ).first; + + // Get SimSeg position and Direction in rZ SL frame + const DTSuperLayer* sl = dtGeom->superLayer(zedRecSeg->superLayerId()); + LocalPoint simSegLocalPosRZTmp = sl->toLocal(simSegmGlobalPos); + LocalVector simSegLocalDirRZ = sl->toLocal(simSegmGlobalDir); + simSegLocalPosRZ = + simSegLocalPosRZTmp + simSegLocalDirRZ*(-simSegLocalPosRZTmp.z()/(cos(simSegLocalDirRZ.theta()))); + alphaSimSegRZ = DTHitQualityUtils::findSegmentAlphaAndBeta(simSegLocalDirRZ).first; + + if (debug) cout << + "RZ SL: recPos " << bestRecHitLocalPosRZ << + "recDir " << bestRecHitLocalDirRZ << + "recAlpha " << alphaBestRHitRZ << endl << + "RZ SL: simPos " << simSegLocalPosRZ << + "simDir " << simSegLocalDirRZ << + "simAlpha " << alphaSimSegRZ << endl ; + } + + // get nhits and t0 + const DTChamberRecSegment2D* phiSeg = bestRecHit->phiSegment(); + + float t0phi = -999; + float t0theta = -999; + int nHitPhi=0; + int nHitTheta=0; + + if (phiSeg) { + t0phi = phiSeg->t0(); + nHitPhi = phiSeg->recHits().size(); + } + + if (zedRecSeg) { + t0theta = zedRecSeg->t0(); + nHitTheta = zedRecSeg->recHits().size(); + } + // Nominal cuts for efficiency - consider all matched segments for the time being +// if(fabs(alphaBestRHit - alphaSimSeg) < 5*sigmaResAlpha && +// fabs(betaBestRHit - betaSimSeg) < 5*sigmaResBeta && +// fabs(bestRecHitLocalPos.x() - xSimSeg) < 5*sigmaResX && +// fabs(bestRecHitLocalPos.y() - ySimSeg) < 5*sigmaResY) { +// recHitFound = true; +// } + recHitFound = true; + + // Mirror alpha in phi SLs so that + and - wheels can be plotted together + if (mirrorMinusWheels && wheel<0){ + alphaSimSeg*=-1.; + alphaBestRHit*=-1.; + // Note: local X (xSimSeg, bestRecHitLocalPos.x() would have to be mirrored as well; + // but at the moment there is no plot of dependency vs X, except for efficiency. + } + // Fill Residual histos HRes4DHit *histo=0; - if((*chamberId).wheel() == 0) + if(wheel == 0) histo = h4DHit_W0; - else if(abs((*chamberId).wheel()) == 1) + else if(abs(wheel) == 1) histo = h4DHit_W1; - else if(abs((*chamberId).wheel()) == 2) + else if(abs(wheel) == 2) histo = h4DHit_W2; - histo->Fill(alphaSimSeg, - alphaBestRHit, - betaSimSeg, - betaBestRHit, - xSimSeg, - bestRecHitLocalPos.x(), - ySimSeg, - bestRecHitLocalPos.y(), - etaSimSeg, - phiSimSeg, - bestRecHitLocalPosRZ.x(), - simSegLocalPosRZ.x(), - alphaBestRHitRZ, - alphaSimSegRZ, - sqrt(bestRecHitLocalDirErr.xx()), - sqrt(bestRecHitLocalDirErr.yy()), - sqrt(bestRecHitLocalPosErr.xx()), - sqrt(bestRecHitLocalPosErr.yy()), - sqrt(bestRecHitLocalDirErrRZ.xx()), - sqrt(bestRecHitLocalPosErrRZ.xx()) - ); + float sigmaAlphaBestRhit = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHit,bestRecHitLocalDirErr.xx())); + float sigmaBetaBestRhit = sqrt(DTHitQualityUtils::sigmaAngle(betaBestRHit,bestRecHitLocalDirErr.yy())); // FIXME this misses the contribution from uncertainty in extrapolation! + float sigmaAlphaBestRhitRZ = sqrt(DTHitQualityUtils::sigmaAngle(alphaBestRHitRZ,bestRecHitLocalDirErrRZ.xx())); + + histo->Fill(alphaSimSeg, + alphaBestRHit, + betaSimSeg, + betaBestRHit, + xSimSeg, + bestRecHitLocalPos.x(), + ySimSeg, + bestRecHitLocalPos.y(), + etaSimSeg, + phiSimSeg, + bestRecHitLocalPosRZ.x(), + simSegLocalPosRZ.x(), + alphaBestRHitRZ, + alphaSimSegRZ, + sigmaAlphaBestRhit, + sigmaBetaBestRhit, + sqrt(bestRecHitLocalPosErr.xx()), + sqrt(bestRecHitLocalPosErr.yy()), + sigmaAlphaBestRhitRZ, + sqrt(bestRecHitLocalPosErrRZ.xx()), + nHitPhi,nHitTheta,t0phi,t0theta + ); h4DHit->Fill(alphaSimSeg, - alphaBestRHit, - betaSimSeg, - betaBestRHit, - xSimSeg, - bestRecHitLocalPos.x(), - ySimSeg, - bestRecHitLocalPos.y(), - etaSimSeg, - phiSimSeg, - bestRecHitLocalPosRZ.x(), - simSegLocalPosRZ.x(), - alphaBestRHitRZ, - alphaSimSegRZ, - sqrt(bestRecHitLocalDirErr.xx()), - sqrt(bestRecHitLocalDirErr.yy()), - sqrt(bestRecHitLocalPosErr.xx()), - sqrt(bestRecHitLocalPosErr.yy()), - sqrt(bestRecHitLocalDirErrRZ.xx()), - sqrt(bestRecHitLocalPosErrRZ.xx()) + alphaBestRHit, + betaSimSeg, + betaBestRHit, + xSimSeg, + bestRecHitLocalPos.x(), + ySimSeg, + bestRecHitLocalPos.y(), + etaSimSeg, + phiSimSeg, + bestRecHitLocalPosRZ.x(), + simSegLocalPosRZ.x(), + alphaBestRHitRZ, + alphaSimSegRZ, + sigmaAlphaBestRhit, + sigmaBetaBestRhit, + sqrt(bestRecHitLocalPosErr.xx()), + sqrt(bestRecHitLocalPosErr.yy()), + sigmaAlphaBestRhitRZ, + sqrt(bestRecHitLocalPosErrRZ.xx()), + nHitPhi,nHitTheta,t0phi,t0theta ); - if (local) h4DHitWS[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(alphaSimSeg, + if (local) h4DHitWS[abs(wheel)][station-1]->Fill(alphaSimSeg, alphaBestRHit, betaSimSeg, betaBestRHit, @@ -414,16 +466,16 @@ void DTSegment4DQuality::endJob() { simSegLocalPosRZ.x(), alphaBestRHitRZ, alphaSimSegRZ, - sqrt(bestRecHitLocalDirErr.xx()), - sqrt(bestRecHitLocalDirErr.yy()), - sqrt(bestRecHitLocalPosErr.xx()), + sigmaAlphaBestRhit, + sigmaBetaBestRhit, + sqrt(bestRecHitLocalPosErr.xx()), sqrt(bestRecHitLocalPosErr.yy()), - sqrt(bestRecHitLocalDirErrRZ.xx()), - sqrt(bestRecHitLocalPosErrRZ.xx()) + sigmaAlphaBestRhitRZ, + sqrt(bestRecHitLocalPosErrRZ.xx()), + nHitPhi,nHitTheta,t0phi,t0theta ); - - } //end of if(bestRecHitFound) + } //end of if(bestRecHit) } //end of if(nsegm!=0) @@ -431,15 +483,15 @@ void DTSegment4DQuality::endJob() { if(doall){ HEff4DHit *heff = 0; - if((*chamberId).wheel() == 0) + if(wheel == 0) heff = hEff_W0; - else if(abs((*chamberId).wheel()) == 1) + else if(abs(wheel) == 1) heff = hEff_W1; - else if(abs((*chamberId).wheel()) == 2) + else if(abs(wheel) == 2) heff = hEff_W2; - heff->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound); - hEff_All->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound); - if (local) hEffWS[abs((*chamberId).wheel())][(*chamberId).station()-1]->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound); + heff->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg); + hEff_All->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg); + if (local) hEffWS[abs(wheel)][station-1]->Fill(etaSimSeg, phiSimSeg, xSimSeg, ySimSeg, alphaSimSeg, betaSimSeg, recHitFound,count_seg); } } // End of loop over chambers diff --git a/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.h b/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.h index eb62491862da1..0226d51657c0b 100644 --- a/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.h +++ b/trunk/Validation/DTRecHits/plugins/DTSegment4DQuality.h @@ -23,7 +23,6 @@ #include "FWCore/Framework/interface/EDAnalyzer.h" #include "Histograms.h" #include "DQMServices/Core/interface/DQMStore.h" -#include "DQMServices/Core/interface/MonitorElement.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Utilities/interface/InputTag.h" @@ -38,7 +37,6 @@ namespace edm { } class TFile; -class MonitorElement; class DTSegment4DQuality : public edm::EDAnalyzer { public: @@ -52,6 +50,9 @@ class DTSegment4DQuality : public edm::EDAnalyzer { /// Perform the real analysis void analyze(const edm::Event & event, const edm::EventSetup& eventSetup); + + virtual void beginRun(const edm::Run& iRun, const edm::EventSetup &setup); + // Write the histos to file void endJob(); void endLuminosityBlock(edm::LuminosityBlock const& lumiSeg, @@ -89,9 +90,6 @@ class DTSegment4DQuality : public edm::EDAnalyzer { HEff4DHit *hEff_W2; HEff4DHit *hEffWS[3][4]; - MonitorElement* hHitMult[3][4]; - MonitorElement* ht0[3][4]; - DQMStore* dbe_; bool doall; bool local; diff --git a/trunk/Validation/DTRecHits/plugins/Histograms.h b/trunk/Validation/DTRecHits/plugins/Histograms.h index 195d3a1f3fc34..17a8117ce0787 100644 --- a/trunk/Validation/DTRecHits/plugins/Histograms.h +++ b/trunk/Validation/DTRecHits/plugins/Histograms.h @@ -37,11 +37,11 @@ class HRes1DHit{ if(doall){ hDist=0; hDist = dbe_->book1D(pre + "_hDist" ,"1D RHit distance from wire", 100, 0,2.5); //hDist = new TH1F ("1D_"+N+"_hDist", "1D RHit distance from wire", 100, 0,2.5); - hResVsAngle = 0; hResVsAngle = dbe_->book2D(pre+"_hResVsAngle", "1D RHit residual vs impact angle",100, 0.,1.2, 150, -0.5,0.5); + hResVsAngle = 0; hResVsAngle = dbe_->book2D(pre+"_hResVsAngle", "1D RHit residual vs impact angle",100, -1.2,1.2, 100, -0.2,0.2); hResVsDistFE = 0; hResVsDistFE = dbe_->book2D(pre+"_hResVsDistFE", "1D RHit residual vs FE distance", 100, 0.,400., 150, -0.5,0.5); dbe_->setCurrentFolder("DT/1DRecHits/Pull/"); hPullVsPos= 0; hPullVsPos = dbe_->book2D (pre+"_hPullVsPos", "1D RHit pull vs position", 100, 0,2.5, 100, -5,5); - hPullVsAngle = 0; hPullVsAngle = dbe_->book2D (pre+"_hPullVsAngle", "1D RHit pull vs impact angle",100, 0.,+1.2, 100, -5,5); + hPullVsAngle = 0; hPullVsAngle = dbe_->book2D (pre+"_hPullVsAngle", "1D RHit pull vs impact angle",100, -1.2,1.2, 100, -5,5); hPullVsDistFE = 0; hPullVsDistFE = dbe_->book2D (pre+"_hPullVsDistFE", "1D RHit pull vs FE distance", 100, 0., 400., 100, -5,5); } dbe_->setCurrentFolder("DT/1DRecHits/Res/"); @@ -335,9 +335,9 @@ class HRes2DHit{ pre += name_; dbe_->setCurrentFolder("DT/2DSegments/Res/"); if(doall){ - hRecAngle=0;hRecAngle = dbe_->book1D (pre+"_hRecAngle", "Distribution of Rec segment angles;angle (rad)",100, -3.5, 3.5); - hSimAngle=0;hSimAngle = dbe_->book1D (pre+"_hSimAngle", "Distribution of segment angles from SimHits;angle (rad)",100, -3.5, 3.5); - hRecVsSimAngle=0;hRecVsSimAngle = dbe_->book2D (pre+"_hRecVsSimAngle", "Rec angle vs sim angle;angle (rad)",100, -3.5, 3.5, 100, -3.5, 3.5); + hRecAngle=0;hRecAngle = dbe_->book1D (pre+"_hRecAngle", "Distribution of Rec segment angles;angle (rad)",100, -1.5, 1.5); + hSimAngle=0;hSimAngle = dbe_->book1D (pre+"_hSimAngle", "Distribution of segment angles from SimHits;angle (rad)",100, -1.5, 1.5); + hRecVsSimAngle=0;hRecVsSimAngle = dbe_->book2D (pre+"_hRecVsSimAngle", "Rec angle vs sim angle;angle (rad)",100, -1.5, 1.5, 100, -1.5, 1.5); hResAngleVsEta=0;hResAngleVsEta = dbe_->book2D (pre+"_hResAngleVsEta", "Residual on 2D segment angle vs Eta; #eta; res (rad)",100, -2.5, 2.5, 200, -0.2, 0.2); hResAngleVsPhi=0;hResAngleVsPhi = dbe_->book2D (pre+"_hResAngleVsPhi", "Residual on 2D segment angle vs Phi; #phi (rad);res (rad)", 100, -3.2, 3.2, 150, -0.2, 0.2); @@ -665,24 +665,24 @@ class HEff2DHit{ // Histos of residuals for 4D rechits class HRes4DHit{ public: - HRes4DHit(std::string name_,DQMStore *dbe_,bool doall=true,bool local=true){ + HRes4DHit(std::string name_,DQMStore *dbe_,bool doall=true,bool local=true) : isLocal(local){ std::string pre ="4D_"; pre += name_; _doall = doall; dbe_->setCurrentFolder("DT/4DSegments/Res/"); if(doall){ - hRecAlpha=0;hRecAlpha = dbe_->book1D (pre+"_hRecAlpha", "4D RecHit alpha (RPhi) distribution;#alpha^{x} (rad)", 100, -3.5, 3.5); - hRecBeta=0;hRecBeta = dbe_->book1D (pre+"_hRecBeta", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -3.5, 3.5); + hRecAlpha=0;hRecAlpha = dbe_->book1D (pre+"_hRecAlpha", "4D RecHit alpha (RPhi) distribution;#alpha^{x} (rad)", 100, -1.5, 1.5); + hRecBeta=0;hRecBeta = dbe_->book1D (pre+"_hRecBeta", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -1.5, 1.5); hSimAlpha=0;hSimAlpha = dbe_->book1D(pre+"_hSimAlpha", "4D segment from SimHit alpha (RPhi) distribution;i#alpha^{x} (rad)", - 100, -3.5, 3.5); + 100, -1.5, 1.5); hSimBeta=0;hSimBeta = dbe_->book1D(pre+"_hSimBeta", "4D segment from SimHit beta distribution;#alpha^{y} (rad)", - 100, -3.5, 3.5); + 100, -1.5, 1.5); hRecVsSimAlpha=0;hRecVsSimAlpha = dbe_->book2D(pre+"_hRecVsSimAlpha", "4D segment rec alpha {v}s sim alpha (RPhi);#alpha^{x} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); + 100, -1.5, 1.5, 100, -1.5, 1.5); hRecVsSimBeta=0;hRecVsSimBeta = dbe_->book2D(pre+"_hRecVsSimBeta", "4D segment rec beta vs sim beta (RZ);#alpha^{y} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); + 100, -1.5, 1.5, 100, -1.5, 1.5); hResAlphaVsEta=0;hResAlphaVsEta = dbe_->book2D (pre+"_hResAlphaVsEta", "4D RecHit residual on #alpha_x direction vs eta;#eta;#alpha^{x}_{rec}-#alpha^{x}_{sim} (rad)", @@ -717,12 +717,12 @@ class HRes4DHit{ hResAlphaVsResY=0;hResAlphaVsResY = dbe_->book2D(pre+"_hResAlphaVsResY", "4D RecHit residual on alpha vs residual on y", 150, -0.6, 0.6, 500, -0.15, 0.15); - hRecBetaRZ=0;hRecBetaRZ = dbe_->book1D (pre+"_hRecBetaRZ", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -3.5, 3.5); + hRecBetaRZ=0;hRecBetaRZ = dbe_->book1D (pre+"_hRecBetaRZ", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -1.5, 1.5); hSimBetaRZ=0;hSimBetaRZ = dbe_->book1D(pre+"_hSimBetaRZ", "4D segment from SimHit beta distribution in RZ SL;#alpha^{y} (rad)", - 100, -3.5, 3.5); + 100, -1.5, 1.5); hRecVsSimBetaRZ=0;hRecVsSimBetaRZ = dbe_->book2D(pre+"_hRecVsSimBetaRZ", "4D segment rec beta vs sim beta (RZ) in RZ SL;#alpha^{y} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); + 100, -1.5, 1.5, 100, -1.5, 1.5); hResBetaVsEtaRZ=0;hResBetaVsEtaRZ = dbe_->book2D (pre+"_hResBetaVsEtaRZ", "4D RecHit residual on beta direction vs eta;#eta in RZ SL;#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", @@ -798,7 +798,6 @@ class HRes4DHit{ 150, -0.15, 0.15); // Pulls - // Pulls dbe_->setCurrentFolder("DT/4DSegments/Pull/"); hPullAlpha=0;hPullAlpha = dbe_->book1D (pre+"_hPullAlpha", @@ -822,84 +821,16 @@ class HRes4DHit{ hPullYRZ=0;hPullYRZ = dbe_->book1D (pre+"_hPullYRZ", "4D RecHit pull on position (y) in chamber in RZ SL;(y_{rec}-y_{sim})/#sigma", 150, -5, 5); - } -/* - HRes4DHit (TString name_, TFile* file){ - name=name_; - - hRecAlpha = (TH1F *) file->Get("DQMData/4D_"+name+"_hRecAlpha"); - hRecBeta = (TH1F *) file->Get("DQMData/4D_"+name+"_hRecBeta"); - - hSimAlpha = (TH1F *) file->Get("DQMData/4D_"+name+"_hSimAlpha"); - hSimBeta = (TH1F *) file->Get("DQMData/4D_"+name+"_hSimBeta"); - - hRecVsSimAlpha = (TH2F *) file->Get("DQMData/4D_"+name+"_hRecVsSimAlpha"); - hRecVsSimBeta = (TH2F *) file->Get("DQMData/4D_"+name+"_hRecVsSimBeta"); - - hResAlpha = (TH1F *) file->Get("DQMData/4D_"+name+"_hResAlpha"); - hResAlphaVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hResAlphaVsEta"); - hResAlphaVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hResAlphaVsPhi"); - - hResBeta = (TH1F *) file->Get("DQMData/4D_"+name+"_hResBeta"); - hResBetaVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hResBetaVsEta"); - hResBetaVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hResBetaVsPhi"); - - hResX = (TH1F *) file->Get("DQMData/4D_"+name+"_hResX"); - hResXVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hResXVsEta"); - hResXVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hResXVsPhi"); - - hResY = (TH1F *) file->Get("DQMData/4D_"+name+"_hResY"); - hResYVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hResYVsEta"); - hResYVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hResYVsPhi"); - - hResAlphaVsResBeta = (TH2F *) file->Get("DQMData/4D_"+name+"_hResAlphaVsResBeta"); - hResXVsResY = (TH2F *) file->Get("DQMData/4D_"+name+"_hResXVsResY"); - hResAlphaVsResX = (TH2F *) file->Get("DQMData/4D_"+name+"_hResAlphaVsResX"); - hResAlphaVsResY = (TH2F *) file->Get("DQMData/4D_"+name+"_hResAlphaVsResY"); - - hPullAlpha = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullAlpha"); - hPullAlphaVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullAlphaVsEta"); - hPullAlphaVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullAlphaVsPhi"); - - hPullBeta = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullBeta"); - hPullBetaVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullBetaVsEta"); - hPullBetaVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullBetaVsPhi"); - hPullX = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullX"); - hPullXVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullXVsEta"); - hPullXVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullXVsPhi"); - - hPullY = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullY"); - hPullYVsEta = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullYVsEta"); - hPullYVsPhi = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullYVsPhi"); - - // RX SL frame - hRecBetaRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hRecBetaRZ"); - - hSimBetaRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hSimBetaRZ"); - - hRecVsSimBetaRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hRecVsSimBetaRZ"); - - hResBetaRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hResBetaRZ"); - hResBetaVsEtaRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hResBetaVsEtaRZ"); - hResBetaVsPhiRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hResBetaVsPhiRZ"); - - hResYRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hResYRZ"); - hResYVsEtaRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hResYVsEtaRZ"); - hResYVsPhiRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hResYVsPhiRZ"); - - hPullBetaRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullBetaRZ"); - hPullBetaVsEtaRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullBetaVsEtaRZ"); - hPullBetaVsPhiRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullBetaVsPhiRZ"); - - hPullYRZ = (TH1F *) file->Get("DQMData/4D_"+name+"_hPullYRZ"); - hPullYVsEtaRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullYVsEtaRZ"); - hPullYVsPhiRZ = (TH2F *) file->Get("DQMData/4D_"+name+"_hPullYVsPhiRZ"); + // NHits, t0 + if (isLocal) { + dbe_->setCurrentFolder("DT/4DSegments/"); + hHitMult = dbe_->book2D(pre+"_hNHits", "NHits", 12,0,12, 6,0,6); + ht0 = dbe_->book2D(pre+"_ht0", "t0", 200,-25,25,200,-25,25); } - ~HRes4DHit(){ } -*/ + void Fill(float simDirectionAlpha, float recDirectionAlpha, float simDirectionBeta, @@ -919,7 +850,11 @@ class HRes4DHit{ float sigmaX, float sigmaY, float sigmaBetaRZ, - float sigmaYRZ + float sigmaYRZ, + int nHitsPhi, + int nHitsTheta, + float t0Phi, + float t0Theta ) { float resAlpha = recDirectionAlpha - simDirectionAlpha; hResAlpha->Fill(resAlpha); @@ -981,62 +916,12 @@ class HRes4DHit{ hPullYVsEtaRZ->Fill(simEta, resYRZ/sigmaYRZ); hPullYVsPhiRZ->Fill(simPhi, resYRZ/sigmaYRZ); } + if (isLocal){ + hHitMult->Fill(nHitsPhi, nHitsTheta); + ht0->Fill(t0Phi,t0Theta); + } } - /*void Write() { - hRecAlpha->Write(); - hRecBeta->Write(); - hSimAlpha->Write(); - hSimBeta->Write(); - hRecVsSimAlpha->Write(); - hRecVsSimBeta->Write(); - hResAlpha->Write(); - hResAlphaVsEta->Write(); - hResAlphaVsPhi->Write(); - hResBeta->Write(); - hResBetaVsEta->Write(); - hResBetaVsPhi->Write(); - hResX->Write(); - hResXVsEta->Write(); - hResXVsPhi->Write(); - hResY->Write(); - hResYVsEta->Write(); - hResYVsPhi->Write(); - hResAlphaVsResBeta->Write(); - hResXVsResY->Write(); - hResAlphaVsResX->Write(); - hResAlphaVsResY->Write(); - hPullAlpha->Write(); - hPullAlphaVsEta->Write(); - hPullAlphaVsPhi->Write(); - hPullBeta->Write(); - hPullBetaVsEta->Write(); - hPullBetaVsPhi->Write(); - hPullX->Write(); - hPullXVsEta->Write(); - hPullXVsPhi->Write(); - hPullY->Write(); - hPullYVsEta->Write(); - hPullYVsPhi->Write(); - - - hRecBetaRZ->Write(); - hSimBetaRZ->Write(); - hRecVsSimBetaRZ->Write(); - hResBetaRZ->Write(); - hResBetaVsEtaRZ->Write(); - hResBetaVsPhiRZ->Write(); - hResYRZ->Write(); - hResYVsEtaRZ->Write(); - hResYVsPhiRZ->Write(); - hPullBetaRZ->Write(); - hPullBetaVsEtaRZ->Write(); - hPullBetaVsPhiRZ->Write(); - hPullYRZ->Write(); - hPullYVsEtaRZ->Write(); - hPullYVsPhiRZ->Write(); - }*/ - public: MonitorElement *hRecAlpha; @@ -1107,8 +992,13 @@ class HRes4DHit{ MonitorElement *hPullYRZ; MonitorElement *hPullYVsEtaRZ; MonitorElement *hPullYVsPhiRZ; + + MonitorElement *hHitMult; + MonitorElement *ht0; + bool _doall; - TString name; + bool isLocal; + TString name; }; //--------------------------------------------------------------------------------------- @@ -1156,56 +1046,13 @@ class HEff4DHit{ 100, -2, 2); hEffVsBeta = 0; - } -/* - HEff4DHit (TString name_, TFile* file){ - name=name_; - hEtaSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hEtaSimSegm"); - hEtaRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hEtaRecHit"); - hEffVsEta = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsEta"); + hNSeg =0; hNSeg = dbe_->book1D(pre+"_hNSeg", "Number of rec segment per sim seg", + 20, 0, 20); - hPhiSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hPhiSimSegm"); - hPhiRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hPhiRecHit"); - hEffVsPhi = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsPhi"); - - hXSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hXSimSegm"); - hXRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hXRecHit"); - hEffVsX = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsX"); - - hYSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hYSimSegm"); - hYRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hYRecHit"); - hEffVsY = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsY"); - - hAlphaSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hAlphaSimSegm"); - hAlphaRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hAlphaRecHit"); - hEffVsAlpha = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsAlpha"); - - hBetaSimSegm = (TH1F *) file->Get("DQMData/4D_"+name+"_hBetaSimSegm"); - hBetaRecHit = (TH1F *) file->Get("DQMData/4D_"+name+"_hBetaRecHit"); - hEffVsBeta = (TH1F *) file->Get("DQMData/4D_"+name+"_hEffVsBeta"); - } -*/ + } ~HEff4DHit(){ - /*delete hEtaSimSegm; - delete hEtaRecHit; - delete hEffVsEta; - delete hPhiSimSegm; - delete hPhiRecHit; - delete hEffVsPhi; - delete hXSimSegm; - delete hXRecHit; - delete hEffVsX; - delete hYSimSegm; - delete hYRecHit; - delete hEffVsY; - delete hAlphaSimSegm; - delete hAlphaRecHit; - delete hEffVsAlpha; - delete hBetaSimSegm; - delete hBetaRecHit; - delete hEffVsBeta;*/ } void Fill(float etaSimSegm, @@ -1214,7 +1061,8 @@ class HEff4DHit{ float ySimSegm, float alphaSimSegm, float betaSimSegm, - bool fillRecHit) { + bool fillRecHit, + int nSeg) { hEtaSimSegm->Fill(etaSimSegm); hPhiSimSegm->Fill(phiSimSegm); @@ -1222,6 +1070,7 @@ class HEff4DHit{ hYSimSegm->Fill(ySimSegm); hAlphaSimSegm->Fill(alphaSimSegm); hBetaSimSegm->Fill(betaSimSegm); + hNSeg->Fill(nSeg); if(fillRecHit) { hEtaRecHit->Fill(etaSimSegm); @@ -1342,34 +1191,6 @@ class HEff4DHit{ } } - /* void Write() { - hEtaSimSegm->Write(); - hEtaRecHit->Write(); - if(hEffVsEta != 0) - hEffVsEta->Write(); - hPhiSimSegm->Write(); - hPhiRecHit->Write(); - if(hEffVsPhi != 0) - hEffVsPhi->Write(); - hXSimSegm->Write(); - hXRecHit->Write(); - if(hEffVsX != 0) - hEffVsX->Write(); - hYSimSegm->Write(); - hYRecHit->Write(); - if(hEffVsY != 0) - hEffVsY->Write(); - hAlphaSimSegm->Write(); - hAlphaRecHit->Write(); - if(hEffVsAlpha != 0) - hEffVsAlpha->Write(); - hBetaSimSegm->Write(); - hBetaRecHit->Write(); - if(hEffVsBeta != 0) - hEffVsBeta->Write(); - - }*/ - public: MonitorElement *hEtaSimSegm; MonitorElement *hEtaRecHit; @@ -1389,6 +1210,7 @@ class HEff4DHit{ MonitorElement *hBetaSimSegm; MonitorElement *hBetaRecHit; TH1F *hEffVsBeta; + MonitorElement *hNSeg; TString name; diff --git a/trunk/Validation/DTRecHits/src/DTHitQualityUtils.cc b/trunk/Validation/DTRecHits/src/DTHitQualityUtils.cc index 7b2b195903494..32bcee4afd85d 100644 --- a/trunk/Validation/DTRecHits/src/DTHitQualityUtils.cc +++ b/trunk/Validation/DTRecHits/src/DTHitQualityUtils.cc @@ -182,15 +182,21 @@ DTHitQualityUtils::findMuSimSegmentDirAndPos(const pair DTHitQualityUtils::findSegmentAlphaAndBeta(const LocalVector& direction) { - //return make_pair(atan(direction.x()/direction.z()), atan(direction.y()/direction.z())); - return make_pair((direction.x()/direction.z()), (direction.y()/direction.z())); + return make_pair(atan(direction.x()/direction.z()), atan(direction.y()/direction.z())); } +//Find error on angle (squared) from localDirectionError, which is the error on tan(Angle) +double DTHitQualityUtils::sigmaAngle(double Angle, double sigma2TanAngle) { + + double XdivZ = tan(Angle); + double sigma2Angle = 1/(1+XdivZ*XdivZ); + sigma2Angle *= sigma2Angle*sigma2TanAngle; + + return sigma2Angle; +} diff --git a/trunk/Validation/DTRecHits/test/DTValidation_RelVal_fromRECO_local_cfg.py b/trunk/Validation/DTRecHits/test/DTValidation_RelVal_fromRECO_local_cfg.py index 21c24cacc8f35..362966b67d580 100644 --- a/trunk/Validation/DTRecHits/test/DTValidation_RelVal_fromRECO_local_cfg.py +++ b/trunk/Validation/DTRecHits/test/DTValidation_RelVal_fromRECO_local_cfg.py @@ -6,8 +6,18 @@ # # Configurable options: -reReco = False # Set this to True to re-reconstruct hits -skipDeltaSuppr = False # Skip DRR (only when reReco=True) +reReco = True # Set this to True to re-reconstruct hits +skipDeltaSuppr = True # Skip DRR (only when reReco=True) + + +doAngleCorr = False + +#SAMPLE = "536" +#SAMPLE = "5312" +#SAMPLE = "620p6" +#SAMPLE = "700p4" +SAMPLE = "700p7" + #### @@ -18,8 +28,6 @@ ## Conditions process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") -process.GlobalTag.globaltag = "START53_V7G::All" -#process.GlobalTag.globaltag = "PRE_ST62_V6::All" process.load("Configuration.StandardSequences.MagneticField_cff") process.load("Configuration.Geometry.GeometryIdeal_cff") @@ -46,6 +54,11 @@ process.seg4dvalidation.doall = True process.seg4dvalidation.local = True +# Debug validation +#process.rechivalidation.debug = True +#process.seg4dvalidation.debug = True + + process.options = cms.untracked.PSet( #FailPath = cms.untracked.vstring('ProductNotFound'), makeTriggerResults = cms.untracked.bool(True), @@ -55,7 +68,11 @@ process.load("FWCore.MessageService.MessageLogger_cfi") process.MessageLogger.cerr.FwkReport.reportEvery = 1000 -process.source = cms.Source("PoolSource", + + +if SAMPLE == "536" : + process.GlobalTag.globaltag = "START53_V7G::All" + process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( '/store/relval/CMSSW_5_3_6-START53_V14/RelValZMM/GEN-SIM-RECO/v2/00000/08C1D822-F629-E211-A6B1-003048679188.root', '/store/relval/CMSSW_5_3_6-START53_V14/RelValZMM/GEN-SIM-RECO/v2/00000/76156813-F529-E211-917B-003048678FA6.root', @@ -66,28 +83,76 @@ '/store/relval/CMSSW_5_3_6-START53_V14/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/v2/00000/7A9F10B4-EB29-E211-88F1-003048FFCBA8.root', '/store/relval/CMSSW_5_3_6-START53_V14/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/v2/00000/64ECECBC-ED29-E211-AB98-002618943939.root' ) -) - -# process.source = cms.Source("PoolSource", -# fileNames = cms.untracked.vstring( -# # '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValSingleMuPt100/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/E4C71BBB-EDBE-E211-8CAF-002590593920.root' + ) -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/1A1EDFF1-D5BE-E211-AE75-003048FFCB9E.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/3E430421-D9BE-E211-B2EB-0026189438A2.root' -# ), -# secondaryFileNames = cms.untracked.vstring( -# # '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValSingleMuPt100/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/005D0C6A-D9BE-E211-A130-0026189438D6.root' - -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/4A3DC569-C9BE-E211-9B0B-003048678ED4.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/5446A469-C3BE-E211-9E1C-00259059642E.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/7025B923-D1BE-E211-8FE8-0026189438AA.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/861065B9-C5BE-E211-B254-003048678A7E.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/B085805B-C3BE-E211-9207-0026189437E8.root', -# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/C637F0BC-C5BE-E211-831A-00261894397B.root' +elif SAMPLE == "5312" : + process.GlobalTag.globaltag = "START53_LV2::All" + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-RECO/START53_LV2-v1/00000/28D552C4-A82B-E311-952C-002590596486.root', + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-RECO/START53_LV2-v1/00000/EA8973F4-A92B-E311-A0A6-00261894396D.root', + ), + secondaryFileNames = cms.untracked.vstring( + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/START53_LV2-v1/00000/1AC2E0AF-912B-E311-B9DE-003048679080.root', + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/START53_LV2-v1/00000/6247C454-912B-E311-BCCC-0025905822B6.root', + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/START53_LV2-v1/00000/CE2C4A01-9F2B-E311-8A2A-00261894395B.root', + '/store/relval/CMSSW_5_3_12_patch2/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/START53_LV2-v1/00000/F6429A09-9F2B-E311-A76F-003048FFCC18.root' + ) + ) -# ) -# ) +elif SAMPLE == "620p6" : + process.GlobalTag.globaltag = "PRE_ST62_V6::All" + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( +# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValSingleMuPt100/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/E4C71BBB-EDBE-E211-8CAF-002590593920.root' + + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/1A1EDFF1-D5BE-E211-AE75-003048FFCB9E.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-RECO/PRE_ST62_V6-v1/00000/3E430421-D9BE-E211-B2EB-0026189438A2.root' + ), + secondaryFileNames = cms.untracked.vstring( +# '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValSingleMuPt100/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/005D0C6A-D9BE-E211-A130-0026189438D6.root' + + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/4A3DC569-C9BE-E211-9B0B-003048678ED4.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/5446A469-C3BE-E211-9E1C-00259059642E.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/7025B923-D1BE-E211-8FE8-0026189438AA.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/861065B9-C5BE-E211-B254-003048678A7E.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/B085805B-C3BE-E211-9207-0026189437E8.root', + '/store/relval/CMSSW_6_2_0_pre6_patch1/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V6-v1/00000/C637F0BC-C5BE-E211-831A-00261894397B.root' + + ) + ) + +elif SAMPLE == "700p4" : + process.GlobalTag.globaltag = "PRE_ST62_V8::All" + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-RECO/PRE_ST62_V8-v1/00000/086DEE6A-1325-E311-BEB2-003048FFD752.root', + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-RECO/PRE_ST62_V8-v1/00000/A8B465D7-1025-E311-B656-003048D3C010.root' + ), + secondaryFileNames = cms.untracked.vstring( + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v1/00000/04C8F739-0225-E311-B778-00261894385A.root', + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v1/00000/429368D4-0125-E311-B80D-002618943882.root', + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v1/00000/8E337AD2-FB24-E311-8208-00261894397E.root', + '/store/relval/CMSSW_7_0_0_pre4/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v1/00000/E47362B8-0025-E311-AA6C-003048FFCBA4.root' + ) + ) + +elif SAMPLE == "700p7" : + #RelValZMM/CMSSW_7_0_0_pre7-PRE_ST62_V8-v2 + process.GlobalTag.globaltag = "PRE_ST62_V8::All" + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-RECO/PRE_ST62_V8-v2/00000/1A18DD97-1946-E311-9689-003048FFD7A2.root', + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-RECO/PRE_ST62_V8-v2/00000/7CEB59F4-0F46-E311-BE5C-0025905964B2.root' + ), + secondaryFileNames = cms.untracked.vstring( + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v2/00000/022794E1-0346-E311-8A49-0026189438FD.root', + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v2/00000/4C9B84CF-0246-E311-BBAD-003048FFCB84.root', + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v2/00000/BE8B527A-0E46-E311-92DC-0026189438F7.root', + '/store/relval/CMSSW_7_0_0_pre7/RelValZMM/GEN-SIM-DIGI-RAW-HLTDEBUG/PRE_ST62_V8-v2/00000/C65E52BB-FC45-E311-B207-00304867916E.root' + ) + ) process.source.inputCommands = cms.untracked.vstring("drop *", "keep PSimHits_g4SimHits_MuonDTHits_SIM", @@ -110,6 +175,9 @@ process.dt4DSegments.Reco4DAlgoConfig.perform_delta_rejecting = False; process.dt4DSegments.Reco4DAlgoConfig.Reco2DAlgoConfig.perform_delta_rejecting = False; +if (doAngleCorr) : + process.dt4DSegments.Reco4DAlgoConfig.recAlgoConfig.doAngleCorr = True; + process.dt4DSegments.Reco4DAlgoConfig.Reco2DAlgoConfig.doAngleCorr = True; #Note: hit recomputation @step2 is not activated anyhow! if (reReco) : #add process.dt2DSegments if needed diff --git a/trunk/Validation/DTRecHits/test/Histograms.h b/trunk/Validation/DTRecHits/test/Histograms.h index 6140d2ee9f7d0..4cd04a91fc8a7 100644 --- a/trunk/Validation/DTRecHits/test/Histograms.h +++ b/trunk/Validation/DTRecHits/test/Histograms.h @@ -3,6 +3,7 @@ /** \class Histograms * Collection of histograms for DT RecHit and Segment test. + * This interface is intended only for reading histogram sets from root files. cf ../plugins/Histograms.h * * \author S. Bolognesi and G. Cerminara - INFN Torino */ @@ -18,84 +19,78 @@ #include #include + +// Retrieve standard name string for histogram sets. Note that validation plots are currently done per abs(wheel). +TString buildName(int wheel, int station, int sl) { + TString name_; + if(sl == 0) { + name_+="W"; + } else if(sl == 2) { + name_+="RZ_W"; + } else { + name_+="RPhi_W"; + } + if (station==0) { + name_+=long(abs(wheel)); + } else { + name_=name_+long(abs(wheel))+"_St"+long(station); + } + return name_; +} + + + + //--------------------------------------------------------------------------------------- /// A set of histograms of residuals and pulls for 1D RecHits class HRes1DHit{ public: - HRes1DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - cout << "constructor Histo"<FIXME: move into a new set - hResSt[0] = new TH1F("1D_"+N + "_hResMB1","1D RHit residual", 300, -0.5,0.5); - hResSt[1] = new TH1F("1D_"+N + "_hResMB2","1D RHit residual", 300, -0.5,0.5); - hResSt[2] = new TH1F("1D_"+N + "_hResMB3","1D RHit residual", 300, -0.5,0.5); - hResSt[3] = new TH1F("1D_"+N + "_hResMB4","1D RHit residual", 300, -0.5,0.5); - //<--- - hResVsEta = new TH2F("1D_"+N+"_hResVsEta", "1D RHit residual vs eta", - 50, -1.25,1.25, 150, -1.5,1.5); - hResVsPhi = new TH2F("1D_"+N+"_hResVsPhi", "1D RHit residual vs phi", - 100, -3.2, 3.2, 150, -1.5,1.5); - hResVsPos = new TH2F("1D_"+N+"_hResVsPos", "1D RHit residual vs position", - 100, 0, 2.5, 150, -1.5,1.5); - hResVsAngle = new TH2F("1D_"+N+"_hResVsAngle", "1D RHit residual vs impact angle", - 100, 0.,1.2, 150, -1.5,1.5); - hResVsDistFE = new TH2F("1D_"+N+"_hResVsDistFE", "1D RHit residual vs FE distance", - 100, 0.,400., 150, -1.5,1.5); - hPull = new TH1F ("1D_"+N+"_hPull", "1D RHit pull", 100, -5,5); - //--->FIXME: move into a new set - hPullSt[0] = new TH1F("1D_"+N + "_hPullMB1","1D RHit residual", 300, -5,5); - hPullSt[1] = new TH1F("1D_"+N + "_hPullMB2","1D RHit residual", 300, -5,5); - hPullSt[2] = new TH1F("1D_"+N + "_hPullMB3","1D RHit residual", 300, -5,5); - hPullSt[3] = new TH1F("1D_"+N + "_hPullMB4","1D RHit residual", 300, -5,5); - //<--- - hPullVsPos = new TH2F ("1D_"+N+"_hPullVsPos", "1D RHit pull vs position", 100, 0,2.5, 100, -5,5); - hPullVsAngle = new TH2F ("1D_"+N+"_hPullVsAngle", "1D RHit pull vs impact angle", - 100, 0.,+1.2, 100, -5,5); - hPullVsDistFE = new TH2F ("1D_"+N+"_hPullVsDistFE", "1D RHit pull vs FE distance",100, 0., 400., 100, -5,5); - hPullVsEta = new TH2F("1D_"+N+"_hPullVsEta", "1D RHit residual vs eta", - 50, -1.25,1.25, 150, -5,5); - hPullVsPhi = new TH2F("1D_"+N+"_hResVsPhi", "1D RHit residual vs phi", - 100, -3.2, 3.2, 150, -5,5); - - } - - HRes1DHit(TString name_, TFile* file){ - name=name_; - hDist = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hDist"); - hRes = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hRes"); - hResSt[0] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB1"); - hResSt[1] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB2"); - hResSt[2] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB3"); - hResSt[3] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB4"); - hResVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsEta"); - hResVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsPhi"); - hResVsPos = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsPos"); - hResVsAngle = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsAngle"); - hResVsDistFE = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsDistFE"); - hPull = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPull"); - hPullSt[0] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB1"); - hPullSt[1] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB2"); - hPullSt[2] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB3"); - hPullSt[3] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB4"); - hPullVsPos = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsPos"); - hPullVsAngle = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsAngle"); - hPullVsDistFE = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsDistFE"); - } - - ~HRes1DHit(){ - // delete hDist; - // delete hRes; - // delete hResVsEta; - // delete hResVsPhi; - // delete hResVsPos; - // delete hPull; - } + HRes1DHit(TFile* file, int wheel, int station, int sl, const TString& step){ + TString name_=step; + name_+=buildName(wheel,station,sl); + initFromFile(name_,file); + } + + + HRes1DHit(TString name_, TFile* file){ + initFromFile(name_,file); + } + + + void initFromFile(TString name_, TFile* file){ + name=name_; + hDist = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hDist"); + hRes = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hRes"); + hResSt[0] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB1"); + hResSt[1] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB2"); + hResSt[2] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB3"); + hResSt[3] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResMB4"); + hResVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsEta"); + hResVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsPhi"); + hResVsPos = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsPos"); + hResVsAngle = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsAngle"); + hResVsDistFE = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Res/1D_"+name+"_hResVsDistFE"); + hPull = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPull"); + hPullSt[0] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB1"); + hPullSt[1] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB2"); + hPullSt[2] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB3"); + hPullSt[3] = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullMB4"); + hPullVsPos = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsPos"); + hPullVsAngle = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsAngle"); + hPullVsDistFE = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/Pull/1D_"+name+"_hPullVsDistFE"); + + if (hRes) { + hRes->SetXTitle("|d_{hit}|-|d_{true}| (cm)"); + hResVsPos->SetXTitle("|d_{true}| (cm)"); + hResVsPos->SetYTitle("|d_{hit}|-|d_{true}| (cm)"); + hResVsAngle->SetXTitle("#alpha_{true} (rad)"); + hResVsAngle->SetYTitle("|d_{hit}|-|d_{true}| (cm)"); + } + } + + + ~HRes1DHit(){} void Fill(float distSimHit, float thetaSimHit, @@ -123,31 +118,6 @@ class HRes1DHit{ else std::cout<<"Error: RecHit error = 0" << std::endl; } - void Write() { - hDist->Write(); - hRes->Write(); - hResSt[0]->Write(); - hResSt[1]->Write(); - hResSt[2]->Write(); - hResSt[3]->Write(); - hResVsEta->Write(); - hResVsPhi->Write(); - hResVsPos->Write(); - hResVsAngle->Write(); - hResVsDistFE->Write(); - hPull->Write(); - hPullSt[0]->Write(); - hPullSt[1]->Write(); - hPullSt[2]->Write(); - hPullSt[3]->Write(); - hPullVsPos->Write(); - hPullVsPhi->Write(); - hPullVsPos->Write(); - hPullVsAngle->Write(); - hPullVsDistFE->Write(); - } - - public: TH1F* hDist; TH1F* hRes; @@ -173,37 +143,21 @@ class HRes1DHit{ /// A set of histograms for efficiency 1D DT RecHits class HEff1DHit{ public: - HEff1DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - - hEtaMuSimHit = new TH1F("1D_"+N+"_hEtaMuSimHit", "SimHit Eta distribution", - 100, -1.5, 1.5); - hEtaRecHit = new TH1F("1D_"+N+"_hEtaRecHit", "SimHit Eta distribution with 1D RecHit", - 100, -1.5, 1.5); - hEffVsEta = 0; - - - hPhiMuSimHit = new TH1F("1D_"+N+"_hPhiMuSimHit", "SimHit Phi distribution", - 100, -TMath::Pi(),TMath::Pi()); - hPhiRecHit = new TH1F("1D_"+N+"_hPhiRecHit", "SimHit Phi distribution with 1D RecHit", - 100, -TMath::Pi(),TMath::Pi()); - hEffVsPhi = 0; - - - hDistMuSimHit = new TH1F("1D_"+N+"_hDistMuSimHit", "SimHit Distance from wire distribution", - 100, 0, 2.5); - hDistRecHit = new TH1F("1D_"+N+"_hDistRecHit", "SimHit Distance from wire distribution with 1D RecHit", - 100, 0, 2.5); - hEffVsDist = 0; - + + HEff1DHit (TFile* file, int wheel, int station, int sl, const TString& step){ + TString name_=step; + name_+=buildName(wheel,station,sl); + initFromFile(name_,file); } + HEff1DHit (TString name_, TFile* file){ + initFromFile(name_,file); + } + + void initFromFile (TString name_, TFile* file){ name=name_; -% cout << "HEff1DHit: 1D_" << name_ << endl; - hEtaMuSimHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/1D_"+name+"_hEtaMuSimHit"); hEtaRecHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/1D_"+name+"_hEtaRecHit"); hEffVsEta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/1DRecHits/1D_"+name+"_hEffVsEta"); @@ -221,24 +175,7 @@ class HEff1DHit{ } - ~HEff1DHit(){ - - // delete hEtaMuSimHit; - // delete hEtaRecHit; - // if(hEffVsEta != 0) - // delete hEffVsEta; - - // delete hPhiMuSimHit; - // delete hPhiRecHit; - // if(hEffVsPhi != 0) - // delete hEffVsPhi; - - // delete hDistMuSimHit; - // delete hDistRecHit; - // if(hEffVsDist != 0) - // delete hEffVsDist; - - } + ~HEff1DHit(){} void Fill(float distSimHit, float etaSimHit, @@ -311,24 +248,8 @@ class HEff1DHit{ } hEffVsDist->SetBinError(bin, error); } - - } - void Write() { - hEtaMuSimHit->Write(); - hEtaRecHit->Write(); - if(hEffVsEta != 0) - hEffVsEta->Write(); - hPhiMuSimHit->Write(); - hPhiRecHit->Write(); - if(hEffVsPhi != 0) - hEffVsPhi->Write(); - hDistMuSimHit->Write(); - hDistRecHit->Write(); - if(hEffVsDist != 0) - hEffVsDist->Write(); - } public: TH1F* hEtaMuSimHit; @@ -352,40 +273,6 @@ class HEff1DHit{ // Histos of residuals for 2D rechits class HRes2DHit{ public: - HRes2DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - - hRecAngle = new TH1F ("2D_"+N+"_hRecAngle", "Distribution of Rec segment angles;angle (rad)", - 100, -3.5, 3.5); - hSimAngle = new TH1F ("2D_"+N+"_hSimAngle", "Distribution of segment angles from SimHits;angle (rad)", - 100, -3.5, 3.5); - hRecVsSimAngle = new TH2F ("2D_"+N+"_hRecVsSimAngle", "Rec angle vs sim angle;angle (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); - - - hResAngle = new TH1F ("2D_"+N+"_hResAngle", "Residual on 2D segment angle;angle_{rec}-angle_{sim} (rad)", 150, -0.15, 0.15); - hResAngleVsEta = new TH2F ("2D_"+N+"_hResAngleVsEta", "Residual on 2D segment angle vs Eta; #eta; res (rad)", - 100, -2.5, 2.5, 200, -0.2, 0.2); - hResAngleVsPhi = new TH2F ("2D_"+N+"_hResAngleVsPhi", "Residual on 2D segment angle vs Phi; #phi (rad);res (rad)", - 100, -3.2, 3.2, 150, -0.2, 0.2); - - hResPos = new TH1F ("2D_"+N+"_hResPos", "Residual on 2D segment position (x at SL center);x_{rec}-x_{sim} (cm)", - 150, -0.2, 0.2); - hResPosVsEta = new TH2F ("2D_"+N+"_hResPosVsEta", "Residual on 2D segment position vs Eta;#eta;res (cm)", - 100, -2.5, 2.5, 150, -0.2, 0.2); - hResPosVsPhi = new TH2F ("2D_"+N+"_hResPosVsPhi", "Residual on 2D segment position vs Phi;#phi (rad);res (cm)", - 100, -3.2, 3.2, 150, -0.2, 0.2); - - hResPosVsResAngle = new TH2F("2D_"+N+"_hResPosVsResAngle", - "Residual on 2D segment position vs Residual on 2D segment angle;angle (rad);res (cm)", - 100, -0.3, 0.3, 150, -0.2, 0.2); - - hPullAngle = new TH1F ("2D_"+N+"_hPullAngle", "Pull on 2D segment angle;(angle_{rec}-angle_{sim})/#sigma (rad)", 150, -5, 5); - hPullPos = new TH1F ("2D_"+N+"_hPullPos", "Pull on 2D segment position (x at SL center);(x_{rec}-x_{sim} (cm))/#sigma", - 150, -5, 5); - } - HRes2DHit (TString name_, TFile* file){ name=name_; @@ -405,8 +292,7 @@ class HRes2DHit{ } - ~HRes2DHit(){ - } + ~HRes2DHit(){} void Fill(float angleSimSegment, float angleRecSegment, @@ -433,23 +319,6 @@ class HRes2DHit{ hPullPos->Fill(resPos/sigmaPos); } - void Write() { - - hRecAngle->Write(); - hSimAngle->Write(); - hRecVsSimAngle->Write(); - hResAngle->Write(); - hResAngleVsEta->Write(); - hResAngleVsPhi->Write(); - hResPos->Write(); - hResPosVsEta->Write(); - hResPosVsPhi->Write(); - hResPosVsResAngle->Write(); - hPullAngle->Write(); - hPullPos->Write(); - } - - public: TH1F *hRecAngle; TH1F *hSimAngle; @@ -472,36 +341,6 @@ class HRes2DHit{ // Histos for 2D RecHit efficiency class HEff2DHit{ public: - HEff2DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - - hEtaSimSegm = new TH1F("2D_"+N+"_hEtaSimSegm", "Eta of SimHit segment", 100, -1.5, 1.5); - hEtaRecHit = new TH1F("2D_"+N+"_hEtaRecHit", "Eta distribution of SimHit segment with 2D RecHit", - 100, -1.5, 1.5); - hEffVsEta = 0; - - hPhiSimSegm = new TH1F("2D_"+N+"_hPhiSimSegm", "Phi of SimHit segment", - 100, -TMath::Pi(),TMath::Pi()); - hPhiRecHit = new TH1F("2D_"+N+"_hPhiRecHit", "Phi distribution of SimHit segment with 2D RecHit", - 100, -TMath::Pi(),TMath::Pi()); - hEffVsPhi = 0; - - - hPosSimSegm = new TH1F("2D_"+N+"_hPosSimSegm", "Position in SL of SimHit segment (cm)", - 100, -250, 250); - hPosRecHit = new TH1F("2D_"+N+"_hPosRecHit", "Position in SL of SimHit segment with 2D RecHit (cm)", - 100, -250, 250); - hEffVsPos = 0; - - - hAngleSimSegm = new TH1F("2D_"+N+"_hAngleSimSegm", "Angle of SimHit segment (rad)", - 100, -2, 2); - hAngleRecHit = new TH1F("2D_"+N+"_hAngleRecHit", "Angle of SimHit segment with 2D RecHit (rad)", - 100, -2, 2); - hEffVsAngle = 0; - - } HEff2DHit (TString name_, TFile* file){ name=name_; @@ -523,27 +362,7 @@ class HEff2DHit{ } - ~HEff2DHit(){ - //delete hEtaSimSegm; - //delete hEtaRecHit; - //if(hEffVsEta != 0) - // delete hEffVsEta; - - //delete hPhiSimSegm; - //delete hPhiRecHit; - //if(hEffVsPhi != 0) - // delete hEffVsPhi; - - //delete hPosSimSegm; - //delete hPosRecHit; - //if(hEffVsPos != 0) - // delete hEffVsPos; - - //delete hAngleSimSegm; - //delete hAngleRecHit; - //if(hEffVsAngle != 0) - // delete hEffVsAngle; - } + ~HEff2DHit(){} void Fill(float etaSimSegm, float phiSimSegm, @@ -634,25 +453,6 @@ class HEff2DHit{ } - void Write() { - hEtaSimSegm->Write(); - hEtaRecHit->Write(); - if(hEffVsEta != 0) - hEffVsEta->Write(); - hPhiSimSegm->Write(); - hPhiRecHit->Write(); - if(hEffVsPhi != 0) - hEffVsPhi->Write(); - hPosSimSegm->Write(); - hPosRecHit->Write(); - if(hEffVsPos != 0) - hEffVsPos->Write(); - hAngleSimSegm->Write(); - hAngleRecHit->Write(); - if(hEffVsAngle != 0) - hEffVsAngle->Write(); - - } public: @@ -676,158 +476,18 @@ class HEff2DHit{ // Histos of residuals for 4D rechits class HRes4DHit{ public: - HRes4DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - - hRecAlpha = new TH1F ("4D_"+N+"_hRecAlpha", "4D RecHit alpha (RPhi) distribution;#alpha^{x} (rad)", 100, -3.5, 3.5); - hRecBeta = new TH1F ("4D_"+N+"_hRecBeta", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -3.5, 3.5); - - hSimAlpha = new TH1F("4D_"+N+"_hSimAlpha", "4D segment from SimHit alpha (RPhi) distribution;i#alpha^{x} (rad)", - 100, -3.5, 3.5); - hSimBeta = new TH1F("4D_"+N+"_hSimBeta", "4D segment from SimHit beta distribution;#alpha^{y} (rad)", - 100, -3.5, 3.5); - hRecVsSimAlpha = new TH2F("4D_"+N+"_hRecVsSimAlpha", "4D segment rec alpha {v}s sim alpha (RPhi);#alpha^{x} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); - hRecVsSimBeta = new TH2F("4D_"+N+"_hRecVsSimBeta", "4D segment rec beta vs sim beta (RZ);#alpha^{y} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); - - hResAlpha = new TH1F ("4D_"+N+"_hResAlpha", - "4D RecHit residual on #alpha_x direction;#alpha^{x}_{rec}-#alpha^{x}_{sim} (rad)", - 200, -0.015, 0.015); - hResAlphaVsEta = new TH2F ("4D_"+N+"_hResAlphaVsEta", - "4D RecHit residual on #alpha_x direction vs eta;#eta;#alpha^{x}_{rec}-#alpha^{x}_{sim} (rad)", - 100, -2.5, 2.5, 100, -0.025, 0.025); - hResAlphaVsPhi = new TH2F ("4D_"+N+"_hResAlphaVsPhi", - "4D RecHit residual on #alpha_x direction vs phi (rad);#phi (rad);#alpha^{x}_{rec}-#alpha^{x}_{sim} (rad)", - 100, -3.2, 3.2, 100, -0.025, 0.025); - - hResBeta = new TH1F ("4D_"+N+"_hResBeta", - "4D RecHit residual on beta direction;#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 200, -0.1, 0.1); - hResBetaVsEta = new TH2F ("4D_"+N+"_hResBetaVsEta", - "4D RecHit residual on beta direction vs eta;#eta;#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 100, -2.5, 2.5, 200, -0.2, 0.2); - hResBetaVsPhi = new TH2F ("4D_"+N+"_hResBetaVsPhi", - "4D RecHit residual on beta direction vs phi;#phi (rad);#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 100, -3.2, 3.2, 200, -0.2, 0.2); - - hResX = new TH1F ("4D_"+N+"_hResX", "4D RecHit residual on position (x) in chamber;x_{rec}-x_{sim} (cm)", - 150, -0.15, 0.15); - hResXVsEta = new TH2F ("4D_"+N+"_hResXVsEta", "4D RecHit residual on position (x) in chamber vs eta;#eta;x_{rec}-x_{sim} (cm)", - 100, -2.5, 2.5, 150, -0.3, 0.3); - hResXVsPhi = new TH2F ("4D_"+N+"_hResXVsPhi", "4D RecHit residual on position (x) in chamber vs phi;#phi (rad);x_{rec}-x_{sim} (cm)", - 100, -3.2, 3.2, 150, -0.3, 0.3); - - hResY = new TH1F ("4D_"+N+"_hResY", "4D RecHit residual on position (y) in chamber;y_{rec}-y_{sim} (cm)", 150, -0.6, 0.6); - hResYVsEta = new TH2F ("4D_"+N+"_hResYVsEta", "4D RecHit residual on position (y) in chamber vs eta;#eta;y_{rec}-y_{sim} (cm)", - 100, -2.5, 2.5, 150, -0.6, 0.6); - hResYVsPhi = new TH2F ("4D_"+N+"_hResYVsPhi", "4D RecHit residual on position (y) in chamber vs phi;#phi (rad);y_{rec}-y_{sim} (cm)", - 100, -3.2, 3.2, 150, -0.6, 0.6); - - hResAlphaVsResBeta = new TH2F("4D_"+N+"_hResAlphaVsResBeta", "4D RecHit residual on alpha vs residual on beta", - 200, -0.3, 0.3, 500, -0.15, 0.15); - hResXVsResY = new TH2F("4D_"+N+"_hResXVsResY", "4D RecHit residual on X vs residual on Y", - 150, -0.6, 0.6, 50, -0.3, 0.3); - hResAlphaVsResX = new TH2F("4D_"+N+"_hResAlphaVsResX", "4D RecHit residual on alpha vs residual on x", - 150, -0.3, 0.3, 500, -0.15, 0.15); - - hResAlphaVsResY = new TH2F("4D_"+N+"_hResAlphaVsResY", "4D RecHit residual on alpha vs residual on y", - 150, -0.6, 0.6, 500, -0.15, 0.15); - - // Pulls - - hPullAlpha = new TH1F ("4D_"+N+"_hPullAlpha", - "4D RecHit pull on #alpha_x direction;(#alpha^{x}_{rec}-#alpha^{x}_{sim})/#sigma", - 200, -5, 5); - hPullAlphaVsEta = new TH2F ("4D_"+N+"_hPullAlphaVsEta", - "4D RecHit pull on #alpha_x direction vs eta;#eta;(#alpha^{x}_{rec}-#alpha^{x}_{sim})/#sigma", - 100, -2.5, 2.5, 100, -5, 5); - hPullAlphaVsPhi = new TH2F ("4D_"+N+"_hPullAlphaVsPhi", - "4D RecHit pull on #alpha_x direction vs phi (rad);#phi (rad);(#alpha^{x}_{rec}-#alpha^{x}_{sim})/#sigma", - 100, -3.2, 3.2, 100, -5, 5); - - hPullBeta = new TH1F ("4D_"+N+"_hPullBeta", - "4D RecHit pull on beta direction;(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 200, -5, 5); - hPullBetaVsEta = new TH2F ("4D_"+N+"_hPullBetaVsEta", - "4D RecHit pull on beta direction vs eta;#eta;(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 100, -2.5, 2.5, 200, -5, 5); - hPullBetaVsPhi = new TH2F ("4D_"+N+"_hPullBetaVsPhi", - "4D RecHit pull on beta direction vs phi;#phi (rad);(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 100, -3.2, 3.2, 200, -5, 5); - - hPullX = new TH1F ("4D_"+N+"_hPullX", - "4D RecHit pull on position (x) in chamber;(x_{rec}-x_{sim})#sigma", - 150, -5, 5); - hPullXVsEta = new TH2F ("4D_"+N+"_hPullXVsEta", - "4D RecHit pull on position (x) in chamber vs eta;#eta;(x_{rec}-x_{sim})#sigma", - 100, -2.5, 2.5, 150, -5, 5); - hPullXVsPhi = new TH2F ("4D_"+N+"_hPullXVsPhi", - "4D RecHit pull on position (x) in chamber vs phi;#phi (rad);(x_{rec}-x_{sim})/#sigma", - 100, -3.2, 3.2, 150, -5, 5); - - hPullY = new TH1F ("4D_"+N+"_hPullY", - "4D RecHit pull on position (y) in chamber;(y_{rec}-y_{sim})/#sigma", 150, -5, 5); - hPullYVsEta = new TH2F ("4D_"+N+"_hPullYVsEta", - "4D RecHit pull on position (y) in chamber vs eta;#eta;(y_{rec}-y_{sim})/#sigma", - 100, -2.5, 2.5, 150, -5, 5); - hPullYVsPhi = new TH2F ("4D_"+N+"_hPullYVsPhi", - "4D RecHit pull on position (y) in chamber vs phi;#phi (rad);(y_{rec}-y_{sim})/#sigma", - 100, -3.2, 3.2, 150, -5, 5); - - // histo in rz SL reference frame. - - hRecBetaRZ = new TH1F ("4D_"+N+"_hRecBetaRZ", "4D RecHit beta distribution:#alpha^{y} (rad)", 100, -3.5, 3.5); - - hSimBetaRZ = new TH1F("4D_"+N+"_hSimBetaRZ", "4D segment from SimHit beta distribution in RZ SL;#alpha^{y} (rad)", - 100, -3.5, 3.5); - hRecVsSimBetaRZ = new TH2F("4D_"+N+"_hRecVsSimBetaRZ", "4D segment rec beta vs sim beta (RZ) in RZ SL;#alpha^{y} (rad)", - 100, -3.5, 3.5, 100, -3.5, 3.5); - - hResBetaRZ = new TH1F ("4D_"+N+"_hResBetaRZ", - "4D RecHit residual on beta direction in RZ SL;#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 200, -0.1, 0.1); - hResBetaVsEtaRZ = new TH2F ("4D_"+N+"_hResBetaVsEtaRZ", - "4D RecHit residual on beta direction vs eta;#eta in RZ SL;#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 100, -2.5, 2.5, 200, -0.2, 0.2); - hResBetaVsPhiRZ = new TH2F ("4D_"+N+"_hResBetaVsPhiRZ", - "4D RecHit residual on beta direction vs phi in RZ SL;#phi (rad);#alpha^{y}_{rec}-#alpha^{y}_{sim} (rad)", - 100, -3.2, 3.2, 200, -0.2, 0.2); - - hResYRZ = new TH1F ("4D_"+N+"_hResYRZ", - "4D RecHit residual on position (y) in chamber in RZ SL;y_{rec}-y_{sim} (cm)", - 150, -0.15, 0.15); - hResYVsEtaRZ = new TH2F ("4D_"+N+"_hResYVsEtaRZ", - "4D RecHit residual on position (y) in chamber vs eta in RZ SL;#eta;y_{rec}-y_{sim} (cm)", - 100, -2.5, 2.5, 150, -0.6, 0.6); - hResYVsPhiRZ = new TH2F ("4D_"+N+"_hResYVsPhiRZ", - "4D RecHit residual on position (y) in chamber vs phi in RZ SL;#phi (rad);y_{rec}-y_{sim} (cm)", - 100, -3.2, 3.2, 150, -0.6, 0.6); - - // Pulls - hPullBetaRZ = new TH1F ("4D_"+N+"_hPullBetaRZ", - "4D RecHit pull on beta direction in RZ SL;(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 200, -5, 5); - hPullBetaVsEtaRZ = new TH2F ("4D_"+N+"_hPullBetaVsEtaRZ", - "4D RecHit pull on beta direction vs eta;#eta in RZ SL;(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 100, -2.5, 2.5, 200, -5, 5); - hPullBetaVsPhiRZ = new TH2F ("4D_"+N+"_hPullBetaVsPhiRZ", - "4D RecHit pull on beta direction vs phi in RZ SL;#phi (rad);(#alpha^{y}_{rec}-#alpha^{y}_{sim})/#sigma", - 100, -3.2, 3.2, 200, -5, 5); - - hPullYRZ = new TH1F ("4D_"+N+"_hPullYRZ", - "4D RecHit pull on position (y) in chamber in RZ SL;(y_{rec}-y_{sim})/#sigma", - 150, -5, 5); - hPullYVsEtaRZ = new TH2F ("4D_"+N+"_hPullYVsEtaRZ", - "4D RecHit pull on position (y) in chamber vs eta in RZ SL;#eta;(y_{rec}-y_{sim})/#sigma", - 100, -2.5, 2.5, 150, -5, 5); - hPullYVsPhiRZ = new TH2F ("4D_"+N+"_hPullYVsPhiRZ", - "4D RecHit pull on position (y) in chamber vs phi in RZ SL;#phi (rad);(y_{rec}-y_{sim})/#sigma", - 100, -3.2, 3.2, 150, -5, 5); + + HRes4DHit (TFile* file, int wheel, int station, int sl){ + initFromFile(buildName(wheel,station,sl),file); } + HRes4DHit (TString name_, TFile* file){ + initFromFile(name_,file); + } + + + void initFromFile (TString name_, TFile* file){ name=name_; hRecAlpha = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hRecAlpha"); @@ -860,21 +520,21 @@ class HRes4DHit{ hResAlphaVsResX = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hResAlphaVsResX"); hResAlphaVsResY = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hResAlphaVsResY"); - hPullAlpha = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullAlpha"); - hPullAlphaVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullAlphaVsEta"); - hPullAlphaVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullAlphaVsPhi"); + hPullAlpha = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullAlpha"); + hPullAlphaVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullAlphaVsEta"); + hPullAlphaVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullAlphaVsPhi"); - hPullBeta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBeta"); - hPullBetaVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBetaVsEta"); - hPullBetaVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBetaVsPhi"); + hPullBeta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBeta"); + hPullBetaVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBetaVsEta"); + hPullBetaVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBetaVsPhi"); - hPullX = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullX"); - hPullXVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullXVsEta"); - hPullXVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullXVsPhi"); + hPullX = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullX"); + hPullXVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullXVsEta"); + hPullXVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullXVsPhi"); - hPullY = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullY"); - hPullYVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullYVsEta"); - hPullYVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullYVsPhi"); + hPullY = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullY"); + hPullYVsEta = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullYVsEta"); + hPullYVsPhi = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullYVsPhi"); // RX SL frame hRecBetaRZ = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hRecBetaRZ"); @@ -891,157 +551,20 @@ class HRes4DHit{ hResYVsEtaRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hResYVsEtaRZ"); hResYVsPhiRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hResYVsPhiRZ"); - hPullBetaRZ = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBetaRZ"); - hPullBetaVsEtaRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBetaVsEtaRZ"); - hPullBetaVsPhiRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullBetaVsPhiRZ"); + hPullBetaRZ = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBetaRZ"); + hPullBetaVsEtaRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBetaVsEtaRZ"); + hPullBetaVsPhiRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullBetaVsPhiRZ"); - hPullYRZ = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullYRZ"); - hPullYVsEtaRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullYVsEtaRZ"); - hPullYVsPhiRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Res/4D_"+name+"_hPullYVsPhiRZ"); - } + hPullYRZ = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullYRZ"); + hPullYVsEtaRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullYVsEtaRZ"); + hPullYVsPhiRZ = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/Pull/4D_"+name+"_hPullYVsPhiRZ"); - ~HRes4DHit(){ + hHitMult = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hNHits"); + ht0 = (TH2F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_ht0"); + } - void Fill(float simDirectionAlpha, - float recDirectionAlpha, - float simDirectionBeta, - float recDirectionBeta, - float simX, - float recX, - float simY, - float recY, - float simEta, - float simPhi, - float recYRZ, - float simYRZ, - float recBetaRZ, - float simBetaRZ, - float sigmaAlpha, - float sigmaBeta, - float sigmaX, - float sigmaY, - float sigmaBetaRZ, - float sigmaYRZ - ) { - - hRecAlpha->Fill(recDirectionAlpha); - hRecBeta->Fill(recDirectionBeta); - hSimAlpha->Fill(simDirectionAlpha); - hSimBeta->Fill(simDirectionBeta); - - hRecVsSimAlpha->Fill(simDirectionAlpha, recDirectionAlpha); - hRecVsSimBeta->Fill(simDirectionBeta, recDirectionBeta); - - float resAlpha = recDirectionAlpha - simDirectionAlpha; - hResAlpha->Fill(resAlpha); - hResAlphaVsEta->Fill(simEta, resAlpha); - hResAlphaVsPhi->Fill(simPhi, resAlpha); - hPullAlpha->Fill(resAlpha/sigmaAlpha); - hPullAlphaVsEta->Fill(simEta, resAlpha/sigmaAlpha); - hPullAlphaVsPhi->Fill(simPhi, resAlpha/sigmaAlpha); - float resBeta = recDirectionBeta - simDirectionBeta; - hResBeta->Fill(resBeta); - hResBetaVsEta->Fill(simEta, resBeta); - hResBetaVsPhi->Fill(simPhi, resBeta); - hPullBeta->Fill(resBeta/sigmaBeta); - hPullBetaVsEta->Fill(simEta, resBeta/sigmaBeta); - hPullBetaVsPhi->Fill(simPhi, resBeta/sigmaBeta); - float resX = recX - simX; - hResX->Fill(resX); - hResXVsEta->Fill(simEta, resX); - hResXVsPhi->Fill(simPhi, resX); - hPullX->Fill(resX/sigmaX); - hPullXVsEta->Fill(simEta, resX/sigmaX); - hPullXVsPhi->Fill(simPhi, resX/sigmaX); - float resY = recY - simY; - hResY->Fill(resY); - hResYVsEta->Fill(simEta, resY); - hResYVsPhi->Fill(simPhi, resY); - hPullY->Fill(resY/sigmaY); - hPullYVsEta->Fill(simEta, resY/sigmaY); - hPullYVsPhi->Fill(simPhi, resY/sigmaY); - - hResAlphaVsResBeta->Fill(resBeta, resAlpha); - hResXVsResY->Fill(resY, resX); - hResAlphaVsResX->Fill(resX, resAlpha); - hResAlphaVsResY->Fill(resY, resAlpha); - - // RZ SuperLayer - hRecBetaRZ->Fill(recBetaRZ); - hSimBetaRZ->Fill(simBetaRZ); - - hRecVsSimBetaRZ->Fill(simBetaRZ, recBetaRZ); - - float resBetaRZ = recBetaRZ - simBetaRZ; - hResBetaRZ->Fill(resBetaRZ); - hResBetaVsEtaRZ->Fill(simEta, resBetaRZ); - hResBetaVsPhiRZ->Fill(simPhi, resBetaRZ); - hPullBetaRZ->Fill(resBetaRZ/sigmaBetaRZ); - hPullBetaVsEtaRZ->Fill(simEta, resBetaRZ/sigmaBetaRZ); - hPullBetaVsPhiRZ->Fill(simPhi, resBetaRZ/sigmaBetaRZ); - float resYRZ = recYRZ - simYRZ; - hResYRZ->Fill(resYRZ); - hResYVsEtaRZ->Fill(simEta, resYRZ); - hResYVsPhiRZ->Fill(simPhi, resYRZ); - hPullYRZ->Fill(resYRZ/sigmaYRZ); - hPullYVsEtaRZ->Fill(simEta, resYRZ/sigmaYRZ); - hPullYVsPhiRZ->Fill(simPhi, resYRZ/sigmaYRZ); - } - - void Write() { - hRecAlpha->Write(); - hRecBeta->Write(); - hSimAlpha->Write(); - hSimBeta->Write(); - hRecVsSimAlpha->Write(); - hRecVsSimBeta->Write(); - hResAlpha->Write(); - hResAlphaVsEta->Write(); - hResAlphaVsPhi->Write(); - hResBeta->Write(); - hResBetaVsEta->Write(); - hResBetaVsPhi->Write(); - hResX->Write(); - hResXVsEta->Write(); - hResXVsPhi->Write(); - hResY->Write(); - hResYVsEta->Write(); - hResYVsPhi->Write(); - hResAlphaVsResBeta->Write(); - hResXVsResY->Write(); - hResAlphaVsResX->Write(); - hResAlphaVsResY->Write(); - hPullAlpha->Write(); - hPullAlphaVsEta->Write(); - hPullAlphaVsPhi->Write(); - hPullBeta->Write(); - hPullBetaVsEta->Write(); - hPullBetaVsPhi->Write(); - hPullX->Write(); - hPullXVsEta->Write(); - hPullXVsPhi->Write(); - hPullY->Write(); - hPullYVsEta->Write(); - hPullYVsPhi->Write(); - - - hRecBetaRZ->Write(); - hSimBetaRZ->Write(); - hRecVsSimBetaRZ->Write(); - hResBetaRZ->Write(); - hResBetaVsEtaRZ->Write(); - hResBetaVsPhiRZ->Write(); - hResYRZ->Write(); - hResYVsEtaRZ->Write(); - hResYVsPhiRZ->Write(); - hPullBetaRZ->Write(); - hPullBetaVsEtaRZ->Write(); - hPullBetaVsPhiRZ->Write(); - hPullYRZ->Write(); - hPullYVsEtaRZ->Write(); - hPullYVsPhiRZ->Write(); - } + ~HRes4DHit(){} public: @@ -1114,60 +637,80 @@ class HRes4DHit{ TH2F *hPullYVsEtaRZ; TH2F *hPullYVsPhiRZ; + TH2F* hHitMult; + TH2F *ht0; + TString name; }; //--------------------------------------------------------------------------------------- /// A set of histograms for efficiency 4D RecHits class HEff4DHit{ - public: - HEff4DHit(std::string name_){ - TString N = name_.c_str(); - name=N; - - hEtaSimSegm = new TH1F("4D_"+N+"_hEtaSimSegm", "Eta of SimHit segment", 100, -1.5, 1.5); - hEtaRecHit = new TH1F("4D_"+N+"_hEtaRecHit", "Eta distribution of SimHit segment with 4D RecHit", - 100, -1.5, 1.5); - hEffVsEta = 0; - - hPhiSimSegm = new TH1F("4D_"+N+"_hPhiSimSegm", "Phi of SimHit segment", - 100, -TMath::Pi(),TMath::Pi()); - hPhiRecHit = new TH1F("4D_"+N+"_hPhiRecHit", "Phi distribution of SimHit segment with 4D RecHit", - 100, -TMath::Pi(),TMath::Pi()); - hEffVsPhi = 0; - - - hXSimSegm = new TH1F("4D_"+N+"_hXSimSegm", "X position in Chamber of SimHit segment (cm)", - 100, -200, 200); - hXRecHit = new TH1F("4D_"+N+"_hXRecHit", "X position in Chamber of SimHit segment with 4D RecHit (cm)", - 100, -200, 200); - hEffVsX = 0; - - hYSimSegm = new TH1F("4D_"+N+"_hYSimSegm", "Y position in Chamber of SimHit segment (cm)", - 100, -200, 200); - hYRecHit = new TH1F("4D_"+N+"_hYRecHit", "Y position in Chamber of SimHit segment with 4D RecHit (cm)", - 100, -200, 200); - hEffVsY = 0; - - hAlphaSimSegm = new TH1F("4D_"+N+"_hAlphaSimSegm", "Alpha of SimHit segment (rad)", - 100, -1.5, 1.5); - hAlphaRecHit = new TH1F("4D_"+N+"_hAlphaRecHit", "Alpha of SimHit segment with 4D RecHit (rad)", - 100, -1.5, 1.5); - hEffVsAlpha = 0; - - hBetaSimSegm = new TH1F("4D_"+N+"_hBetaSimSegm", "Beta of SimHit segment (rad)", - 100, -2, 2); - hBetaRecHit = new TH1F("4D_"+N+"_hBetaRecHit", "Beta of SimHit segment with 4D RecHit (rad)", - 100, -2, 2); - hEffVsBeta = 0; - - } - - HEff4DHit (TString name_, TFile* file){ - name=name_; - hEtaSimSegm = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEtaSimSegm"); - hEtaRecHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEtaRecHit"); - hEffVsEta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEffVsEta"); + + public: + + HEff4DHit (TFile* file, int wheel, int station, int sl){ + initFromFile(buildName(wheel,station,sl),file); + } + + + HEff4DHit (TString name_, TFile* file){ + initFromFile(name_,file); + } + + + + HEff4DHit(std::string name_){ + TString N = name_.c_str(); + name=N; + + hEtaSimSegm = new TH1F("4D_"+N+"_hEtaSimSegm", "Eta of SimHit segment", 100, -1.5, 1.5); + hEtaRecHit = new TH1F("4D_"+N+"_hEtaRecHit", "Eta distribution of SimHit segment with 4D RecHit", + 100, -1.5, 1.5); + hEffVsEta = 0; + + hPhiSimSegm = new TH1F("4D_"+N+"_hPhiSimSegm", "Phi of SimHit segment", + 100, -TMath::Pi(),TMath::Pi()); + hPhiRecHit = new TH1F("4D_"+N+"_hPhiRecHit", "Phi distribution of SimHit segment with 4D RecHit", + 100, -TMath::Pi(),TMath::Pi()); + hEffVsPhi = 0; + + + hXSimSegm = new TH1F("4D_"+N+"_hXSimSegm", "X position in Chamber of SimHit segment (cm)", + 100, -200, 200); + hXRecHit = new TH1F("4D_"+N+"_hXRecHit", "X position in Chamber of SimHit segment with 4D RecHit (cm)", + 100, -200, 200); + hEffVsX = 0; + + hYSimSegm = new TH1F("4D_"+N+"_hYSimSegm", "Y position in Chamber of SimHit segment (cm)", + 100, -200, 200); + hYRecHit = new TH1F("4D_"+N+"_hYRecHit", "Y position in Chamber of SimHit segment with 4D RecHit (cm)", + 100, -200, 200); + hEffVsY = 0; + + hAlphaSimSegm = new TH1F("4D_"+N+"_hAlphaSimSegm", "Alpha of SimHit segment (rad)", + 100, -1.5, 1.5); + hAlphaRecHit = new TH1F("4D_"+N+"_hAlphaRecHit", "Alpha of SimHit segment with 4D RecHit (rad)", + 100, -1.5, 1.5); + hEffVsAlpha = 0; + + hBetaSimSegm = new TH1F("4D_"+N+"_hBetaSimSegm", "Beta of SimHit segment (rad)", + 100, -2, 2); + hBetaRecHit = new TH1F("4D_"+N+"_hBetaRecHit", "Beta of SimHit segment with 4D RecHit (rad)", + 100, -2, 2); + hEffVsBeta = 0; + + hNSeg = new TH1F(" 4D_"+N+"_hNSeg", "Number of rec segment per sim seg", + 20, 0, 20); + + } + + void initFromFile (TString name_, TFile* file){ + + name=name_; + hEtaSimSegm = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEtaSimSegm"); + hEtaRecHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEtaRecHit"); + hEffVsEta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEffVsEta"); hPhiSimSegm = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hPhiSimSegm"); hPhiRecHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hPhiRecHit"); @@ -1188,38 +731,24 @@ class HEff4DHit{ hBetaSimSegm = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hBetaSimSegm"); hBetaRecHit = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hBetaRecHit"); hEffVsBeta = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hEffVsBeta"); - } + hNSeg = (TH1F *) file->Get("DQMData/Run 1/DT/Run summary/4DSegments/4D_"+name+"_hNSeg"); + + ComputeEfficiency(); - ~HEff4DHit(){ - - /*delete hEtaSimSegm; - delete hEtaRecHit; - delete hEffVsEta; - delete hPhiSimSegm; - delete hPhiRecHit; - delete hEffVsPhi; - delete hXSimSegm; - delete hXRecHit; - delete hEffVsX; - delete hYSimSegm; - delete hYRecHit; - delete hEffVsY; - delete hAlphaSimSegm; - delete hAlphaRecHit; - delete hEffVsAlpha; - delete hBetaSimSegm; - delete hBetaRecHit; - delete hEffVsBeta;*/ } + + ~HEff4DHit(){} + void Fill(float etaSimSegm, float phiSimSegm, float xSimSegm, float ySimSegm, float alphaSimSegm, float betaSimSegm, - bool fillRecHit) { + bool fillRecHit, + int nSeg) { hEtaSimSegm->Fill(etaSimSegm); hPhiSimSegm->Fill(phiSimSegm); @@ -1227,6 +756,7 @@ class HEff4DHit{ hYSimSegm->Fill(ySimSegm); hAlphaSimSegm->Fill(alphaSimSegm); hBetaSimSegm->Fill(betaSimSegm); + hNSeg->Fill(nSeg); if(fillRecHit) { hEtaRecHit->Fill(etaSimSegm); @@ -1341,33 +871,6 @@ class HEff4DHit{ } } - void Write() { - hEtaSimSegm->Write(); - hEtaRecHit->Write(); - if(hEffVsEta != 0) - hEffVsEta->Write(); - hPhiSimSegm->Write(); - hPhiRecHit->Write(); - if(hEffVsPhi != 0) - hEffVsPhi->Write(); - hXSimSegm->Write(); - hXRecHit->Write(); - if(hEffVsX != 0) - hEffVsX->Write(); - hYSimSegm->Write(); - hYRecHit->Write(); - if(hEffVsY != 0) - hEffVsY->Write(); - hAlphaSimSegm->Write(); - hAlphaRecHit->Write(); - if(hEffVsAlpha != 0) - hEffVsAlpha->Write(); - hBetaSimSegm->Write(); - hBetaRecHit->Write(); - if(hEffVsBeta != 0) - hEffVsBeta->Write(); - - } public: TH1F *hEtaSimSegm; @@ -1388,7 +891,8 @@ class HEff4DHit{ TH1F *hBetaSimSegm; TH1F *hBetaRecHit; TH1F *hEffVsBeta; - + TH1F *hNSeg; + TString name; }; diff --git a/trunk/Validation/DTRecHits/test/macros.C b/trunk/Validation/DTRecHits/test/macros.C index d77f07295f98f..669d8275062f7 100644 --- a/trunk/Validation/DTRecHits/test/macros.C +++ b/trunk/Validation/DTRecHits/test/macros.C @@ -3,12 +3,6 @@ * * N. Amapane 2002-2004 */ -/* - * Draw 2-D plots superimposed to their profiles - * - * 2003 NCA - */ -// Draw a 2-D plot within the specified Y range and superimpose its X profile #include #include @@ -69,34 +63,133 @@ void setStyle(TH2 *histo) { histo->GetYaxis()->SetLabelSize(gStyle->GetLabelSize()); } -void plotAndProfileX (TH2* h2, float min, float max, bool profile=false) { - setStyle(h2); + +bool addProfile=false; +bool addSlice=true; +bool addMedian = false; +TString opt2Dplot = "col"; + +// Plot a TH2 + add profiles on top of it +// minY, maxY: Y range for plotting and for computing profile if addProfile==true. +// Note that the simple profile is very sensitive to the Y range used! +TH1F* plotAndProfileX (TH2* theh, int rebinX, int rebinY, int rebinProfile, float minY, float maxY, float minX=0, float maxX=0) { + TH2* h2=theh->Clone(); + + // setStyle(h2); + if (h2==0) { + cout << "plotAndProfileX: null histo ptr" << endl; + return; + } + gPad->SetGrid(1,1); gStyle->SetGridColor(15); - h2->GetYaxis()->SetRangeUser(min,max); - h2->Draw(); - if (profile) { - TProfile* prof = h2->ProfileX(); + h2->Rebin2D(rebinX,rebinY); + // h2->GetYaxis()->SetRangeUser(minY,maxY); + + TLine * l = new TLine(h2->GetXaxis()->GetXmin(),0,h2->GetXaxis()->GetXmax(),0); + if (maxX>minX) { + h2->GetXaxis()->SetRangeUser(minX,maxX); + l->SetX1(minX); + l->SetX2(maxX); + } + + h2->SetMarkerStyle(1); + h2->Draw(opt2Dplot); + l->SetLineColor(3); + l->SetLineWidth(2); + l->Draw(); + if (addProfile) { + TAxis* yaxis = h2->GetYaxis(); + //Add option "s" to draw RMS as error instead than RMS/sqrt(N) + TProfile* prof = h2->ProfileX("_pfx", + TMath::Max(1,yaxis->FindBin(minY)), + TMath::Min(yaxis->GetNbins(),yaxis->FindBin(maxY))); +// cout << yaxis->FindBin(minY) << " " << yaxis->FindBin(maxY) << endl; +// cout << yaxis->GetNbins(); + //TProfile* prof = h2->ProfileX("_pfx"); prof->SetMarkerColor(2); + prof->SetMarkerStyle(20); + prof->SetMarkerSize(0.4); prof->SetLineColor(2); + prof->Rebin(rebinProfile); prof->Draw("same"); } - TLine * l = new TLine(h2->GetXaxis()->GetXmin(),0,h2->GetXaxis()->GetXmax(),0); - l->SetLineColor(3); - l->Draw(); + + TH1F* ht=0; + + if (addSlice) { + TObjArray aSlices; + // TF1 fff("a", "gaus", -0.1, 0.1); + h2->FitSlicesY(0, 0, -1, 0, "QNR", &aSlices); // add "G2" to merge 2 consecutive bins + + TH1F* ht = (TH1F*) aSlices[1]->Clone(); + // Remove bins with failed fits, based on fit errors + float thr = (maxY-minY)/4.; + for (int bin=1; bin<=ht->GetNbinsX();++bin){ + if (ht->GetBinError(bin)>thr) { + ht->SetBinContent(bin,0); + ht->SetBinError(bin,0); + } + } + ht->SetMarkerColor(4); + ht->Draw("same"); + } + + if (addMedian) { + double xq[1] = {0.5}; + double median[1]; + + TAxis* axis = h2->GetXaxis(); + TH1F* medprof = new TH1F(h2->GetName()+TString("medians"),"medians", axis->GetNbins(), axis->GetXmin(), axis->GetXmax()); + float bw = h2->GetYaxis()->GetBinLowEdge(2)-h2->GetYaxis()->GetBinLowEdge(1); + + + TString projname = h2->GetName()+TString("_pmedian"); + for (int bin=1; bin<=h2->GetNbinsX(); ++bin){ + TH1D * proj = h2->ProjectionY(projname, bin, bin); + double integral = proj->Integral(); + if (integral==0) continue; + // Take overflow and underflow into account + int nbins = proj->GetNbinsX(); + proj->SetBinContent(1, proj->GetBinContent(0)+proj->GetBinContent(1)); + proj->SetBinContent(0,0); + proj->SetBinContent(nbins, proj->GetBinContent(nbins)+proj->GetBinContent(nbins+1)); + proj->SetBinContent(nbins+1,0); + proj->GetQuantiles(1,median,xq); + medprof->SetBinContent(bin,median[0]); + // Approximated uncertainty on median, probably underestimated. + medprof->SetBinError(bin,bw*sqrt(integral/2.)/2./TMath::Max(1.,proj->GetBinContent(proj->FindBin(median[0])))); + } + medprof->SetMarkerColor(2); + medprof->SetMarkerStyle(20); + medprof->SetMarkerSize(0.4); + medprof->Draw("Esame"); + } + + h2->GetYaxis()->SetRangeUser(minY,maxY); + + return ht; } -// void plotAndProfileY (TH2* h2, float min, float max) { + +// void plotAndProfileX (TH2* h2, float min, float max, bool profile=false) { +// setStyle(h2); +// gPad->SetGrid(1,1); +// gStyle->SetGridColor(15); // h2->GetYaxis()->SetRangeUser(min,max); // h2->Draw(); -// TProfile* prof = h2->ProfileY(); -// prof->SetMarkerStyle(8); -// prof->SetMarkerSize(0.7); -// prof->SetMarkerColor(2); -// prof->SetLineColor(2); -// prof->Draw("same"); +// if (profile) { +// TProfile* prof = h2->ProfileX(); +// prof->SetMarkerColor(2); +// prof->SetLineColor(2); +// prof->Draw("same"); +// } +// TLine * l = new TLine(h2->GetXaxis()->GetXmin(),0,h2->GetXaxis()->GetXmax(),0); +// l->SetLineColor(3); +// l->Draw(); // } + // Draw a 2-D plot within the specified Y range and superimpose its X profile, // setting as sigmas that of the fit (and not the error of the mean) void plotAndProfileXSpread (TH2* h2, float min, float max, bool profile=false, float ymin=-5., float ymax=5.) { @@ -128,44 +221,47 @@ void plotAndProfileXSpread (TH2* h2, float min, float max, bool profile=false, f l->SetLineColor(3); l->Draw(); } -/* - * Draw and format a fitted histogram - * - * 2003 NCA - */ -// Fit a histogram with a gaussian and draw it in the range -// mean+-nsigmas*RMS. -void drawGFit(TH1 * h1, float nsigmas, float min, float max){ + + + +// Fit the gaussian core of an histogram with in the range mean+-nsigmas. +TF1* drawGFit(TH1 * h1, float nsigmas, float min, float max){ + + gPad->SetGrid(1,1); + gStyle->SetGridColor(15); + h1->GetXaxis()->SetRangeUser(min,max); float minfit = h1->GetMean() - h1->GetRMS(); float maxfit = h1->GetMean() + h1->GetRMS(); - drawGFit(h1, min, max, minfit, maxfit); - gPad->Draw(); -} -// Fit a histogram with a gaussian and draw it in the specified range. -void drawGFit(TH1 * h1, float min, float max){ - drawGFit(h1, min, max, min, max); - gPad->Draw(); -} -// Fit a histogram in the range (minfit, maxfit) with a gaussian and -// draw it in the range (min, max) -void drawGFit(TH1 * h1, float min, float max, float minfit, float maxfit) { - setStyle(h1); + + TLine * l = new TLine(0,0,0,h1->GetMaximum()*1.05); + + l->SetLineColor(3); + l->SetLineWidth(2); + static int i = 0; + TString nameF1 = TString("g") + (Long_t)i; i++; - gPad->SetGrid(1,1); - gStyle->SetGridColor(15); - h1->GetXaxis()->SetRangeUser(min,max); - TString fitName = "g"; - fitName += i; - TF1* g1 = new TF1(fitName.Data(),"gaus",minfit,maxfit); + TF1* g1 = new TF1(nameF1,"gaus",minfit,maxfit); + g1->SetLineColor(2); g1->SetLineWidth(2); - h1->Fit(g1,"R"); - h1->Draw(); -// TPaveStats *st = (TPaveStats*)h1->GetListOfFunctions()->FindObject("stats"); -// st->SetX2NDC(0.905); -// st->SetY2NDC(0.905); + h1->Fit(g1,"RQ"); + + minfit = g1->GetParameter("Mean") - nsigmas*g1->GetParameter("Sigma"); + maxfit = g1->GetParameter("Mean") + nsigmas*g1->GetParameter("Sigma"); + g1->SetRange(minfit,maxfit); + + h1->Fit(g1,"RQ"); + TF1* fh=h1->GetFunction(nameF1); + if (fh) fh->FixParameter(0,g1->GetParameter(0)); // so that it is not shown in legend + + gPad->Draw(); + l->Draw(); + h1->Draw("same"); //redraw on top of the line + return g1; } + + /* * Create a new TCanvas setting its properties * @@ -243,15 +339,17 @@ void printCanvasesEps(){ void printCanvasesEps2() { gROOT->GetListOfCanvases()->Print(".eps"); } -// Print all canvases in separate EPS files -void printCanvases(TString type=".eps"){ + +// Print all canvases in separate files +void printCanvases(TString type="png", TString path="."){ TIter iter(gROOT->GetListOfCanvases()); TCanvas *c; while( (c = (TCanvas *)iter()) ) { - c->cd(); - c->Print(0,type); + TString name = c->GetTitle(); + c->Print(path+"/"+name+"."+type,type); } } + /* * Define different TStyles; use them with: * getStyle->cd(); @@ -261,6 +359,9 @@ void printCanvases(TString type=".eps"){ TStyle * getStyle(TString name="myStyle") { TStyle *theStyle; + + gROOT->ForceStyle(); + if ( name == "myStyle" ) { theStyle = new TStyle("myStyle", "myStyle"); // theStyle->SetOptStat(0); @@ -283,8 +384,11 @@ TStyle * getStyle(TString name="myStyle") // For the canvas: theStyle->SetCanvasBorderMode(0); theStyle->SetCanvasColor(kWhite); - theStyle->SetCanvasDefH(600); //Height of canvas - theStyle->SetCanvasDefW(600); //Width of canvas +// theStyle->SetCanvasDefH(600); //Height of canvas +// theStyle->SetCanvasDefW(800); //Width of canvas + theStyle->SetCanvasDefH(750); //Height of canvas + theStyle->SetCanvasDefW(1000); //Width of canvas + theStyle->SetCanvasDefX(0); //POsition on screen theStyle->SetCanvasDefY(0); @@ -310,17 +414,21 @@ TStyle * getStyle(TString name="myStyle") // For the histo: // theStyle->SetHistFillColor(1); // theStyle->SetHistFillStyle(0); - theStyle->SetHistLineColor(1); - theStyle->SetHistLineStyle(0); - theStyle->SetHistLineWidth(1); + theStyle->SetHistLineColor(kBlue); + theStyle->SetMarkerColor(kBlue); + // theStyle->SetHistLineStyle(0); + // theStyle->SetHistLineWidth(1); // theStyle->SetLegoInnerR(Float_t rad = 0.5); // theStyle->SetNumberContours(Int_t number = 20); - theStyle->SetEndErrorSize(2); + + theStyle->SetEndErrorSize(2); // theStyle->SetErrorMarker(20); - theStyle->SetErrorX(0.); - +// theStyle->SetErrorX(0.); + theStyle->SetMarkerStyle(20); + theStyle->SetMarkerSize(0.5); + //For the fit/function: theStyle->SetOptFit(1); @@ -337,57 +445,60 @@ TStyle * getStyle(TString name="myStyle") // For the statistics box: theStyle->SetOptFile(0); // theStyle->SetOptStat(0); // To display the mean and RMS: SetOptStat("mr"); - theStyle->SetOptStat(10); + + theStyle->SetOptStat("e"); theStyle->SetStatColor(kWhite); - theStyle->SetStatFont(42); - theStyle->SetStatFontSize(0.07); + // theStyle->SetStatFont(42); + // theStyle->SetStatFontSize(0.05); theStyle->SetStatTextColor(1); theStyle->SetStatFormat("6.4g"); theStyle->SetStatBorderSize(1); - theStyle->SetStatH(0.3); - theStyle->SetStatW(0.2); +// theStyle->SetStatH(0.02); +// theStyle->SetStatW(0.2); // theStyle->SetStatStyle(Style_t style = 1001); - // theStyle->SetStatX(Float_t x = 0); - // theStyle->SetStatY(Float_t y = 0); + theStyle->SetStatX(0.94); + theStyle->SetStatY(0.96); // Margins: - theStyle->SetPadTopMargin(0.05); - theStyle->SetPadBottomMargin(0.13); - theStyle->SetPadLeftMargin(0.16); - theStyle->SetPadRightMargin(0.02); +// theStyle->SetPadTopMargin(0.1); + theStyle->SetPadBottomMargin(0.11); +// theStyle->SetPadLeftMargin(0.1); +// theStyle->SetPadRightMargin(0.05); + theStyle->SetPadLeftMargin(0.15); // For the Global title: - - theStyle->SetOptTitle(0); - theStyle->SetTitleFont(42); - theStyle->SetTitleColor(1); - theStyle->SetTitleTextColor(1); - theStyle->SetTitleFillColor(10); - theStyle->SetTitleFontSize(0.05); + + // theStyle->SetOptTitle(0); // Uncomment to remove title +// theStyle->SetTitleFont(42); +// theStyle->SetTitleColor(1); +// theStyle->SetTitleTextColor(1); + theStyle->SetTitleFillColor(0); +// theStyle->SetTitleFontSize(0.05); // theStyle->SetTitleH(0); // Set the height of the title box // theStyle->SetTitleW(0); // Set the width of the title box // theStyle->SetTitleX(0); // Set the position of the title box - // theStyle->SetTitleY(0.985); // Set the position of the title box - // theStyle->SetTitleStyle(Style_t style = 1001); - // theStyle->SetTitleBorderSize(2); + theStyle->SetTitleY(0.96); // Set the position of the title box + theStyle->SetTitleStyle(0); + theStyle->SetTitleBorderSize(0); + // For the axis titles: - theStyle->SetTitleColor(1, "XYZ"); - theStyle->SetTitleFont(42, "XYZ"); - theStyle->SetTitleSize(0.06, "XYZ"); +// theStyle->SetTitleColor(1, "XYZ"); +// theStyle->SetTitleFont(42, "XYZ"); + // theStyle->SetTitleSize(0.05, "XYZ"); // theStyle->SetTitleXSize(Float_t size = 0.02); // Another way to set the size? // theStyle->SetTitleYSize(Float_t size = 0.02); - theStyle->SetTitleXOffset(0.9); - theStyle->SetTitleYOffset(1.25); +// theStyle->SetTitleXOffset(0.9); +// theStyle->SetTitleYOffset(1.25); // theStyle->SetTitleOffset(1.1, "Y"); // Another way to set the Offset // For the axis labels: - theStyle->SetLabelColor(1, "XYZ"); - theStyle->SetLabelFont(42, "XYZ"); - theStyle->SetLabelOffset(0.007, "XYZ"); - theStyle->SetLabelSize(0.045, "XYZ"); +// theStyle->SetLabelColor(1, "XYZ"); +// theStyle->SetLabelFont(42, "XYZ"); +// theStyle->SetLabelOffset(0.007, "XYZ"); +// theStyle->SetLabelSize(0.045, "XYZ"); // For the axis: @@ -416,10 +527,11 @@ TStyle * getStyle(TString name="myStyle") // theStyle->SetPalette(Int_t ncolors = 0, Int_t* colors = 0); // theStyle->SetTimeOffset(Double_t toffset); // theStyle->SetHistMinimumZero(kTRUE); - - + theStyle->SetTextSize(0.045); + // theStyle->SetTextFont(42); + // style->SetOptFit(101); - // style->SetOptStat(1111111); + // style->SetOptStat(1111111); } else { // Avoid modifying the default style! @@ -427,3 +539,88 @@ TStyle * getStyle(TString name="myStyle") } return theStyle; } + + +setPalette() +{ + const Int_t NRGBs = 5; + const Int_t NCont = 255; + +// { // Fine rainbow +// Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; +// Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; +// Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; +// Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; +// } + +// { // blues +// Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; +// Double_t red[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 }; +// Double_t green[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 }; +// Double_t blue[NRGBs] = { 1.00, 1.00, 1.00, 1.00, 1.00 }; +// } + + +// { // Gray (white->black) +// Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; +// Double_t red[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 }; +// Double_t green[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 }; +// Double_t blue[NRGBs] = { 1.00, 0.84, 0.61, 0.34, 0.00 }; +// } + + + { // Gray (white->gray) + // similar to gStyle->SetPalette(5); + float max = 0.3; + float step=(1-max)/(NRGBs-1); + Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; + Double_t red[NRGBs] = { 1.00, 1-step, 1-2*step, 1-3*step, 1-4*step }; + Double_t green[NRGBs] = { 1.00, 1-step, 1-2*step, 1-3*step, 1-4*step }; + Double_t blue[NRGBs] = { 1.00, 1-step, 1-2*step, 1-3*step, 1-4*step }; + } + + + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + gStyle->SetNumberContours(NCont); +} + + +// Overlay efficiency plots for S1, S2, S3 +void plotEff(TH1* h1, TH1* h2=0, TH1* h3=0) { + float minY=0.6; + float maxY=1.05; + h1->GetYaxis()->SetRangeUser(minY,maxY); + h1->GetXaxis()->SetRangeUser(0.,2.1); + h1->SetMarkerStyle(20); + h1->SetMarkerSize(0.5); + h1->SetStats(0); + h1->Draw(); + + if (h2) { + h2->SetLineColor(kRed); + h2->SetMarkerColor(kRed); + h2->SetMarkerStyle(20); + h2->SetMarkerSize(0.5); + h2->SetStats(0); + h2->Draw("same"); + } + + if (h3) { + h3->SetLineColor(kBlack); + h3->SetMarkerColor(kBlack); + h3->SetMarkerStyle(20); + h3->SetMarkerSize(0.5); + h3->SetStats(0); + h3->Draw("same"); + } +} + +TH1F* getEffPlot(TH1* hnum, TH1* hden, int rebin=1) { + hnum->Rebin(rebin); + hden->Rebin(rebin); + TH1F* h = hnum->Clone(); + h->Sumw2(); + h->Divide(hnum,hden,1.,1.,"B"); + h->SetStats(0); + return h; +} diff --git a/trunk/Validation/DTRecHits/test/plotValidation.r b/trunk/Validation/DTRecHits/test/plotValidation.r new file mode 100644 index 0000000000000..a5a99e2624ad6 --- /dev/null +++ b/trunk/Validation/DTRecHits/test/plotValidation.r @@ -0,0 +1,376 @@ +//------------------------------ +// +// A macro to decorate validation plots produced with "local" option. +// +// Usage: +// +// .x plotValidation.r("file",wheel,station) +// +// Some drawing options are set in the code (see below). +// +// +// Author: N. Amapane +// +//------------------------------ + + +void plotValidation(){ + cout << endl << "Usage: .x plotValidation.r(\"inputFile.root\",,)" << endl << endl; +} + + +void plotValidation(TString filename, int wheel, int station) { + + if (! TString(gSystem->GetLibraries()).Contains("Histograms_h")) { + gROOT->LoadMacro("$CMSSW_BASE/src/Validation/DTRecHits/test/Histograms.h+"); + gROOT->LoadMacro("macros.C"); + } + + //---------------------------------------------------------------------- + // Configurable options + addProfile = false; + addSlice = true; + + int rbx =2; // rebin x in scatter plots + int rby =1; // rebin y in scatter plots + int rbp = 1; // rebin profiles + float nsigma = 2; // interval for the fit of residual distributions + + // Canvases to plot + bool doPhiAndThetaS3 =true; + bool doHitPull = true; + bool doEff = true; + bool doT0= true; + bool doSegRes=true; + bool doSegPull = true; + bool doAngularDeps = true; + bool doEff4D = true; + + + //---------------------------------------------------------------------- + + TStyle * style = getStyle("tdr"); + style->cd(); + setPalette(); + gStyle->SetTitleSize(0.05,"XYZ"); // Set larger axis titles + gStyle->SetTitleOffset(1.3,"Y"); + // gStyle->SetOptTitle(0); // remove histogram titles + + // gStyle->SetOptStat(111); + + float cmToMicron = 10000.; + float vdrift = 54.3; + + + TFile *file = new TFile(filename); + + HRes1DHit *hResPhi1 = new HRes1DHit(file, wheel, station, 1, "S3"); + HRes1DHit *hResTheta = new HRes1DHit(file, wheel, station, 2, "S3"); + // HRes1DHit *hResPhi2 = new HRes1DHit(file, wheel, station, 4, "S3"); + HRes1DHit *hResPhi = hResPhi1; + + HEff1DHit* hEffS1RPhi= new HEff1DHit(file, wheel, station, 1, "S1"); + HEff1DHit* hEffS3RPhi= new HEff1DHit(file, wheel, station, 1, "S3"); + HEff1DHit* hEffS1RZ=0; + HEff1DHit* hEffS3RZ=0; + if (station!=4) { + hEffS1RZ= new HEff1DHit(file, wheel, station, 2, "S1"); + hEffS3RZ= new HEff1DHit(file, wheel, station, 2, "S3"); + } + + + HRes4DHit* hRes4D= new HRes4DHit(file, wheel, station, 0); + HEff4DHit* hEff4D = new HEff4DHit(file, wheel, station, 0); + + + // Result of fits + float m_phi = 0.; + float s_phi = 0.; + float m_theta = 0.; + float s_theta = 0.; + float m_phi1 = 0.; + float s_phi1 = 0.; + float m_phi2 = 0.; + float s_phi2 = 0.; + float m_phiS1 = 0.; + float s_phiS1 = 0.; + float m_phiS2 = 0.; + float s_phiS2 = 0.; + float m_thetaS1 = 0.; + float s_thetaS1 = 0.; + float m_thetaS2 = 0.; + float s_thetaS2 = 0.; + float t0phi =0.; + float t0theta =0.; + + + TString canvbasename = filename; + canvbasename = canvbasename.Replace(canvbasename.Length()-5,5,"") + TString("_W") + (long) wheel + "_St" + (long) station ; + + + //-------------------- Hit Residuals at step 3 in phi and theta (full distrib and vs distance from wire) + if (doPhiAndThetaS3) { + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_ResPhiTheta"); + c1->SetName(canvbasename+"_ResPhiTheta"); + c1->Divide(2,2); + + c1->cd(1); + + TH1F* hRes; + hRes=hResPhi->hRes; + + hResPhi->hRes->Rebin(2); + TF1* fphi=drawGFit(hRes, nsigma, -0.4, 0.4); + + c1->cd(2); + + plotAndProfileX(hResPhi->hResVsPos,rbx,rby,rbp,-.1, .1, 0, 2.1); + + + m_phi = fphi->GetParameter("Mean")*cmToMicron; + s_phi = fphi->GetParameter("Sigma")*cmToMicron; + + if (hResTheta->hRes) { + + c1->cd(3); + hRes=hResTheta->hRes; + + hResTheta->hRes->Rebin(2); + TF1* ftheta=drawGFit(hRes, nsigma, -0.4, 0.4); + + c1->cd(4); + plotAndProfileX(hResTheta->hResVsPos,rbx,rby,rbp,-.1, .1, 0, 2.1); + + m_theta = ftheta->GetParameter("Mean")*cmToMicron; + s_theta = ftheta->GetParameter("Sigma")*cmToMicron; + } + + + cout << canvbasename << " Step3 W" << wheel << " St" << station << endl + << " Res: Phi: M= " << int(floor(m_phi+0.5)) + << " S= " << int(s_phi+0.5) + << "; Theta: M= " << int(floor(m_theta+0.5)) + << " S= " << int(s_theta+0.5) << endl; + } + + + //-------------------- Hit pulls + if (doHitPull){ + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_PullPhiTheta"); + c1->SetName(canvbasename+"_PullPhiTheta"); + c1->Divide(2,2); + + c1->cd(1); + + TH1F* hPull; + hPull=hResPhi->hPull; + + // hResPhi->hPull->Rebin(2); + TF1* fphi=drawGFit(hPull, nsigma, -5, 5); + + c1->cd(2); + + plotAndProfileX(hResPhi->hPullVsPos,rbx,rby,rbp,-5, 5, 0, 2.1); + + + m_phi = fphi->GetParameter("Mean")*cmToMicron; + s_phi = fphi->GetParameter("Sigma")*cmToMicron; + + if (hResTheta->hPull) { + + c1->cd(3); + hPull=hResTheta->hPull; + + // hResTheta->hPull->Rebin(2); + TF1* ftheta=drawGFit(hPull, nsigma, -5, 5); + + c1->cd(4); + plotAndProfileX(hResTheta->hPullVsPos,rbx,rby,rbp,-5, 5, 0, 2.1); + + m_theta = ftheta->GetParameter("Mean"); + s_theta = ftheta->GetParameter("Sigma"); + } + + + cout << canvbasename << " Step3 W" << wheel << " St" << station << endl + << " Pulls: Phi: M= " << int(floor(m_phi+0.5)) + << " S= " << int(s_phi+0.5) + << "; Theta: M= " << int(floor(m_theta+0.5)) + << " S= " << int(s_theta+0.5) << endl; + } + + + //-------------------- Hit efficiencies as a function of distance from wire + if (doEff) { + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_EffPhiTheta"); + c1->SetName(canvbasename+"_EffPhiTheta"); + c1->SetWindowSize(325,750); + c1->Divide(1,2); + c1->cd(1); + plotEff(hEffS1RPhi->hEffVsDist, hEffS3RPhi->hEffVsDist); + c1->cd(2); + if (station!=4) plotEff(hEffS1RZ->hEffVsDist, hEffS3RZ->hEffVsDist); +// c1->cd(2); +// //plotEff(hEffS1RPhi->hEffVsPhi, hEffS3RPhi->hEffVsPhi); +// plotEff(hEffS1RPhi->hEffVsEta, hEffS3RPhi->hEffVsEta); +// c1->cd(4); +// // if (station!=4) plotEff(hEffS1RZ->hEffVsPhi, hEffS3RZ->hEffVsPhi); +// if (station!=4) plotEff(hEffS1RZ->hEffVsEta, hEffS3RZ->hEffVsEta); + + } + + + //-------------------- #hits, t0s + if (doT0) { + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_NHitsT0"); + c1->SetName(canvbasename+"_NHitsT0"); + c1->Divide(2,2); + c1->cd(1); + + TH2F* hNh = hRes4D->hHitMult; + + hNh->SetXTitle("#phi hits"); + hNh->SetYTitle("#theta hits"); + hNh->Draw("BOX"); + c1->cd(2); + hNh->ProjectionY()->Draw(); + c1->cd(3); + hNh->ProjectionX()->Draw(); + c1->cd(4); + + TH2F* ht = hRes4D->ht0; + ht->SetXTitle("t0 #phi"); + ht->SetYTitle("t0 #theta"); + ht->Draw("BOX"); + + + } + + //-------------------- Segment x, y, alpha, beta resolutions + if (doSegRes){ + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_ResSeg"); + c1->SetName(canvbasename+"_ResSeg"); + c1->Divide(2,2); + c1->cd(1); + // hRes4D->hResX->Rebin(2); + drawGFit(hRes4D->hResX, nsigma, -0.1, 0.1); + c1->cd(2); + // hRes4D->hResAlpha->Rebin(2); + drawGFit(hRes4D->hResAlpha, nsigma, -0.01, 0.01); + c1->cd(3); + hRes4D->hResYRZ->Rebin(2); + drawGFit(hRes4D->hResYRZ, nsigma, -0.4, 0.4); + c1->cd(4); + hRes4D->hResBeta->Rebin(2); + drawGFit(hRes4D->hResBeta, nsigma, -0.4, 0.4); + } + + + //-------------------- Angular dependencies + if (doAngularDeps) { + TCanvas* c1= new TCanvas; + + float min; + float max; + c1->SetTitle(canvbasename+"_MeanVsAngles"); + c1->SetName(canvbasename+"_MeanVsAngles"); + c1->SetTitle(canvbasename+" Angles"); + + c1->Divide(2,2); + + c1->cd(1); + hRes4D->hSimAlpha->SetLineColor(kGreen); + hRes4D->hSimAlpha->Draw(); + hRes4D->hRecAlpha->Draw("same"); + + c1->cd(2); + plotAndProfileX(hResPhi->hResVsAngle,1,1,1,-.04, .04, -1., 1.); + + c1->cd(3); + hRes4D->hSimBetaRZ->SetLineColor(kGreen); + hRes4D->hSimBetaRZ->Draw(); + hRes4D->hRecBetaRZ->Draw("same"); + + c1->cd(4); + plotAndProfileX(hResTheta->hResVsAngle,1,1,1,-.04, .04, -1.2.,1.2); + + + } + + + + //------------------- Efficiencies Vs X, Y. alpha, beta + + if(doEff4D){ + // FIXME: should rebin histograms. + + TH1F* hEffX; + TH1F* hEffY; + TH1F* hEffalpha; + TH1F* hEffbeta; + + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_Efficiencies"); + c1->SetName(canvbasename+"_Efficiencies"); + c1->SetTitle(canvbasename+"_Efficiencies"); + + c1->Divide(2,2); + + c1->cd(1); + hEffX = getEffPlot(hEff4D->hXRecHit, hEff4D->hXSimSegm,2); + hEffX->Draw(); + + c1->cd(2); + hEffalpha = getEffPlot(hEff4D->hAlphaRecHit, hEff4D->hAlphaSimSegm,2); + hEffalpha->Draw(); + + c1->cd(3); + hEffY = getEffPlot(hEff4D->hYRecHit, hEff4D->hYSimSegm,2); + hEffY->Draw(); + + c1->cd(4); + hEffbeta = getEffPlot(hEff4D->hBetaRecHit, hEff4D->hBetaSimSegm,2); + hEffbeta->Draw(); + + + + } + + //-------------------- Segment x, y, alpha, beta pull + if (doSegPull){ + TCanvas* c1= new TCanvas; + c1->SetTitle(canvbasename+"_PullSeg"); + c1->SetName(canvbasename+"_PullSeg"); + c1->Divide(2,2); + c1->cd(1); + + hRes4D->hPullX->Rebin(2); + drawGFit(hRes4D->hPullX, nsigma, -10.,10.); + c1->cd(2); + hRes4D->hPullAlpha->Rebin(2); + drawGFit(hRes4D->hPullAlpha, nsigma, -10., 10.); + c1->cd(3); + hRes4D->hPullYRZ->Rebin(2); + drawGFit(hRes4D->hPullYRZ, nsigma, -10., 10.); + c1->cd(4); + hRes4D->hPullBetaRZ->Rebin(2); + drawGFit(hRes4D->hPullBetaRZ, nsigma, -10.,10.); + + //Fixme: Move these to another canvas. Note that the error used for hPullY is not computed correctly.q +// c1->cd(3); +// hRes4D->hPullY->Rebin(2); +// drawGFit(hRes4D->hPullY, nsigma, -10., 10.); +// c1->cd(4); +// hRes4D->hPullBeta->Rebin(2); +// drawGFit(hRes4D->hPullBeta, nsigma, -10.,10.); + + + } + + +} diff --git a/trunk/Validation/DTRecHits/test/summaryPlot.gnu b/trunk/Validation/DTRecHits/test/summaryPlot.gnu new file mode 100644 index 0000000000000..1c8b7f24ff535 --- /dev/null +++ b/trunk/Validation/DTRecHits/test/summaryPlot.gnu @@ -0,0 +1,154 @@ +# Produce summary plots from the tables generated with writeSummaryTable.r. +# +# Requires gnuplot 4.6. The font path can be set with the environment variable GDFONTPATH. +# + +reset +set macro + +#Set the input file here +myfile = '" + + +void writeSummaryTable(){ + cout << endl << "Usage: .x writeSummaryTable.r(\"inputFile.root\",doEff,doEffab,doRes,doResab,doPull)" << endl << endl; +} + + +void writeSummaryTable(TString filename, bool doEff=true, bool doEffab = true, bool doRes=true, bool doResab = true, bool doPull=true,bool doPullabxy=true, bool doSegRes=true) { + + //---------------------------------------------------------------------- + // Configurable options + + // Interval for the fit of residual distributions + float nsigma = 2.; + + //---------------------------------------------------------------------- + + if (! TString(gSystem->GetLibraries()).Contains("Histograms_h")) { + gROOT->LoadMacro("$CMSSW_BASE/src/Validation/DTRecHits/test/Histograms.h+"); + gROOT->LoadMacro("macros.C"); + } + + float cmToMicron = 10000.; + + TFile *file = new TFile(filename); + + TString table = filename.Replace(filename.Length()-5,5,"_summary.txt"); + ofstream f(table,ios_base::out); + f << fixed; + f << "# W St sec SL effS1RPhi effS3RPhi effSeg resHit pullHit meanAngle sigmaAngle meanSegPos sigmaSegPos" << endl; + + // All sectors are collapsed together as of now + int smin = 0; + int smax = 0; + + for (int wheel = -2; wheel<3; ++wheel) { + for (int station =1; station<=4; ++station) { + for (int sector = smin; sector<=smax; ++sector) { + + if (station!=4 && sector>12) continue; + + double stheta = -1.; + double sphi = -1.; + double ptheta = -1.; + double pphi = -1.; + double salpha_res = -1.; + double sbeta_res = -1.; + double malpha_res = -1.; + double mbeta_res = -1.; + double salpha_pull = -1.; + double sbeta_pull = -1.; + double malpha_pull = -1.; + double mbeta_pull = -1.; + double sx = -1.; + double sy = -1.; + double mx = -1.; + double my = -1.; + double sx_pull = -1.; + double sy_pull = -1.; + double mx_pull = -1.; + double my_pull = -1.; + + float effS1RPhi=0.; + float effS3RPhi=0.; + float effS1RZ=0.; + float effS3RZ=0.; + float effSeg = 0; + + + + HRes1DHit *hResPhi1 = new HRes1DHit(file, wheel, station, 1, "S3"); + HRes1DHit *hResTheta = new HRes1DHit(file, wheel, station, 2, "S3"); + + HEff1DHit* hEffS1RPhi= new HEff1DHit(file, wheel, station, 1, "S1"); + HEff1DHit* hEffS3RPhi= new HEff1DHit(file, wheel, station, 1, "S3"); + HEff1DHit* hEffS1RZ=0; + HEff1DHit* hEffS3RZ=0; + + HRes4DHit* hRes4D= new HRes4DHit(file, wheel, station, 0); + HEff4DHit* hEff4D = new HEff4DHit(file, wheel, station, 0); + + + if (station!=4) { + hEffS1RZ= new HEff1DHit(file, wheel, station, 2, "S1"); + hEffS3RZ= new HEff1DHit(file, wheel, station, 2, "S3"); + } + + + + TCanvas c; + + //--- Hit efficiency + if (doEff) { + if (station!=4) { + effS1RZ = hEffS1RZ->hDistRecHit->Integral()/hEffS1RZ->hDistMuSimHit->Integral(); + effS3RZ = hEffS3RZ->hDistRecHit->Integral()/hEffS3RZ->hDistMuSimHit->Integral(); + } + effS1RPhi = hEffS1RPhi->hDistRecHit->Integral()/hEffS1RPhi->hDistMuSimHit->Integral(); + effS3RPhi = hEffS3RPhi->hDistRecHit->Integral()/hEffS3RPhi->hDistMuSimHit->Integral(); + + cout << " " << wheel << " " << station << " " << sector << " effPhi: " << effS1RPhi << " " << effS3RPhi; + if (station!=4) cout << " effZ: " << effS1RZ << " " << effS3RZ; + cout << endl; + } + + + //--- Segment efficiency + if (doEffab) { + if (station!=4) { + + effSeg = hEff4D->hBetaRecHit->Integral()/hEff4D->hBetaSimSegm->Integral(); + + cout << " " << wheel << " " << station << " " << sector << " effSeg: " << effSeg << endl; + + + } + } + + + //--- Hit resolution + if (doRes) { + if (station!=4) { + TH1F* tmpTheta = (TH1F*) hResTheta->hRes->Clone("tmpTheta"); + tmpTheta->Rebin(2); + TF1* ftheta = drawGFit(tmpTheta, nsigma, -2. , 2.); + stheta = ftheta->GetParameter("Sigma"); + } + TH1F* tmpPhi = (TH1F*) hResPhi1->hRes->Clone("tmpPhi"); + tmpPhi->Rebin(2); + TF1* fphi = drawGFit(tmpPhi, nsigma, -2. , 2. ); + sphi = fphi->GetParameter("Sigma"); + + cout << " " << wheel << " " << station << " " << sector << " sphi: " << sphi; + if (station!=4) cout << " stheta: " << stheta ; + cout << endl; + } + + + //--- Hit pull + if (doPull) { + if (station!=4) { + TH1F* tmpTheta = (TH1F*) hResTheta->hPull->Clone("tmpTheta"); + TF1* ftheta = drawGFit(tmpTheta, nsigma, -2. , 2.); + ptheta = ftheta->GetParameter("Sigma"); + } + TH1F* tmpPhi = (TH1F*) hResPhi1->hPull->Clone("tmpPhi"); + TF1* fphi = drawGFit(tmpPhi, nsigma, -2. , 2. ); + pphi = fphi->GetParameter("Sigma"); + + cout << " " << wheel << " " << station << " " << sector + << " pphi: " << pphi; + if (station!=4) { + cout << " ptheta: " << ptheta ; + } + cout << endl; + } + + + //--- Segment position resolution + if (doSegRes) { + if (station!=4) { + TH1F* tmpY = (TH1F*) hRes4D->hResYRZ->Clone("tmpY"); + TF1* fy = drawGFit(tmpY, nsigma, -2. , 2. ); + sy = fy->GetParameter("Sigma"); + my = fy->GetParameter("Mean"); + } + + TH1F* tmpX = (TH1F*) hRes4D->hResX->Clone("tmpX"); + TF1* fx = drawGFit(tmpX, nsigma, -2. , 2.); + sx = fx->GetParameter("Sigma"); + mx = fx->GetParameter("Mean"); + } + + //--- Segment angle resolution + if (doResab) { + if (station!=4) { + TH1F* tmpBeta = (TH1F*) hRes4D->hResBeta->Clone("tmpBeta"); + tmpBeta->Rebin(2); + TF1* fbeta = drawGFit(tmpBeta, nsigma, -2. , 2.); + sbeta_res = fbeta->GetParameter("Sigma"); + mbeta_res = fbeta->GetParameter("Mean"); + } + + TH1F* tmpAlpha = (TH1F*) hRes4D->hResAlpha->Clone("tmpAlpha"); + // tmpAlpha->Rebin(2); + TF1* falpha = drawGFit(tmpAlpha, nsigma, -2. , 2. ); + salpha_res = falpha->GetParameter("Sigma"); + malpha_res = falpha->GetParameter("Mean"); + + cout << " " << wheel << " " << station << " " << sector << " salpha_res: " << salpha_res << " malpha_res: " << malpha_res; + cout << " sbeta_res: " << sbeta_res << " mbeta_res: " << mbeta_res; + cout << endl; + + } + + // Segment angle pull + if (doPullabxy) { + if (station!=4) { + TH1F* tmpBeta = (TH1F*) hRes4D->hPullBetaRZ->Clone("tmpBeta"); + tmpBeta->Rebin(2); + TF1* fbeta = drawGFit(tmpBeta, nsigma, -2. , 2.); + sbeta_pull = fbeta->GetParameter("Sigma"); + mbeta_pull = fbeta->GetParameter("Mean"); + + TH1F* tmpY = (TH1F*) hRes4D->hPullYRZ->Clone("tmpY"); + tmpY->Rebin(2); + TF1* fy = drawGFit(tmpY, nsigma, -2. , 2. ); + sy_pull = fy->GetParameter("Sigma"); + my_pull = fy->GetParameter("Mean"); + } + + TH1F* tmpAlpha = (TH1F*) hRes4D->hPullAlpha->Clone("tmpAlpha"); + tmpAlpha->Rebin(2); + TF1* falpha = drawGFit(tmpAlpha, nsigma, -2. , 2. ); + salpha_pull = falpha->GetParameter("Sigma"); + malpha_pull = falpha->GetParameter("Mean"); + + TH1F* tmpX = (TH1F*) hRes4D->hPullX->Clone("tmpX"); + tmpX->Rebin(2); + TF1* fx = drawGFit(tmpX, nsigma, -2. , 2.); + sx_pull = fx->GetParameter("Sigma"); + mx_pull = fx->GetParameter("Mean"); + + + cout << " " << wheel << " " << station << " " << sector << " salpha_pull: " << salpha_pull << " malpha_pull: " << malpha_pull; + cout << " sbeta_pull: " << sbeta_pull << " mbeta_pull: " << mbeta_pull; + cout << " sx_pull: " << sx_pull << " mx_pull: " << mx_pull; + cout << " sy_pull: " << sy_pull << " my_pull: " << my_pull; + cout << endl; + + } + + + // Write summary table (one row per SL) + int secmin=sector; + int secmax=sector; + if (sector ==0) { + secmin=1; + secmax=14; + } + for (int sec = secmin; sec<=secmax; sec++) { + if (station!=4 && sec>12) continue; + f << wheel << " " << station << " " << sec << " " << 1 << " " << effS1RPhi << " " << effS3RPhi << " " << effSeg << " " << sphi << " " << pphi << " " << malpha_res << " " <