Commit cdb9745a authored by Tony Tong's avatar Tony Tong Committed by Graeme Stewart
Browse files

fix bug in phiMaximum over substracting (MuonLayerHough-00-00-33)

	* Fixed bug in phiMaximum over substracting
	* Tagged as MuonLayerHough-00-00-33

2016-09-17 Baojia Tong
	* Fixed bug in returning BEE chamber position
	* Add extrapolation function, with linear and parabolic options
	* Tagged as MuonLayerHough-00-00-32

2016-07-21  scott snyder  <snyder@bnl.gov>
	* Tagging MuonLayerHough-00-00-31.
	* Fix clang warnings.
	* Tagging MuonLayerHough-00-00-30.
	* Fix abs() calls.
	* Tagging MuonLayerHough-00-00-29.
	* Comply with ATLAS naming conventions.
parent 8917aaa5
......@@ -20,7 +20,8 @@ namespace MuonHough {
static const int UNINITIALIZED = -99999;
/// struct containing additional debug information on the hits that is not needed for the actual alg
/// but very useful for debugging
struct HitDebugInfo {
class HitDebugInfo {
public:
HitDebugInfo();
HitDebugInfo( int type_, int sector_, Muon::MuonStationIndex::DetectorRegionIndex region_,
......@@ -51,7 +52,8 @@ namespace MuonHough {
/// struct containing all hit information needed for the Hough transform
struct Hit {
class Hit {
public:
/// constructor, takes ownership of the HitDebugInfo pointer
Hit( int layer_, float x_, float ymin_, float ymax_, float w_, HitDebugInfo* d_ = 0,
......
......@@ -64,37 +64,37 @@ namespace MuonHough {
std::vector<MuonHough::MuonLayerHoughSelector> m_selectors;
std::vector<Plots*> m_hMaximaHeightPerChIndex;
bool DEBUG;
bool DEBUG_seg;
TH1F* h_dtheta;
TH1F* h_dtheta_truth;
TH1F* h_diff_MI_e;
TH1F* h_diff_MO_e;
TH1F* h_diff_MI_b;
TH1F* h_diff_MO_b;
TH1F* h_diff_MI_e_truth;
TH1F* h_diff_MO_e_truth;
TH1F* h_diff_MI_b_truth;
TH1F* h_diff_MO_b_truth;
TH1F* h_expt_MI_e;
TH1F* h_expt_MO_e;
TH1F* h_expt_MI_b;
TH1F* h_expt_MO_b;
TH1F* h_expt_MI_e_truth;
TH1F* h_expt_MO_e_truth;
TH1F* h_expt_MI_b_truth;
TH1F* h_expt_MO_b_truth;
TH1F* h_comp_MI_e;
TH1F* h_comp_MO_e;
TH1F* h_comp_MI_b;
TH1F* h_comp_MO_b;
TH1F* h_comp_MI_e_truth;
TH1F* h_comp_MO_e_truth;
TH1F* h_comp_MI_b_truth;
TH1F* h_comp_MO_b_truth;
bool m_DEBUG;
bool m_DEBUG_seg;
TH1F* m_h_dtheta;
TH1F* m_h_dtheta_truth;
TH1F* m_h_diff_MI_e;
TH1F* m_h_diff_MO_e;
TH1F* m_h_diff_MI_b;
TH1F* m_h_diff_MO_b;
TH1F* m_h_diff_MI_e_truth;
TH1F* m_h_diff_MO_e_truth;
TH1F* m_h_diff_MI_b_truth;
TH1F* m_h_diff_MO_b_truth;
TH1F* m_h_expt_MI_e;
TH1F* m_h_expt_MO_e;
TH1F* m_h_expt_MI_b;
TH1F* m_h_expt_MO_b;
TH1F* m_h_expt_MI_e_truth;
TH1F* m_h_expt_MO_e_truth;
TH1F* m_h_expt_MI_b_truth;
TH1F* m_h_expt_MO_b_truth;
TH1F* m_h_comp_MI_e;
TH1F* m_h_comp_MO_e;
TH1F* m_h_comp_MI_b;
TH1F* m_h_comp_MO_b;
TH1F* m_h_comp_MI_e_truth;
TH1F* m_h_comp_MO_e_truth;
TH1F* m_h_comp_MI_b_truth;
TH1F* m_h_comp_MO_b_truth;
MuonHough::MuonDetectorHough detectorHough;
MuonHough::MuonDetectorHough detectorHoughTruth;
MuonHough::MuonDetectorHough m_detectorHough;
MuonHough::MuonDetectorHough m_detectorHoughTruth;
};
......
......@@ -66,8 +66,11 @@ namespace MuonHough {
const MuonLayerHough* hough; // pointer to the corresponding hough
bool isEndcap() const{
Muon::MuonStationIndex::DetectorRegionIndex region = hough->m_descriptor.region;;
if (region != Muon::MuonStationIndex::Barrel) {return true;}
// refers to the chamber orientation!!!! so BEE is a barell in this def
Muon::MuonStationIndex::DetectorRegionIndex region = hough->m_descriptor.region;
Muon::MuonStationIndex::ChIndex chIndex = hough->m_descriptor.chIndex;
//need to make sure BEE's reference plane is the same as barrel
if (region != Muon::MuonStationIndex::Barrel && chIndex != Muon::MuonStationIndex::BEE) {return true;}
return false;
};
float getGlobalR() const{
......@@ -195,5 +198,7 @@ namespace MuonHough {
fill(hit.x,hit.ymin,hit.w);
}
//function for linear/parabolic extrapolation from maxima to maxima in a different station
float extrapolate(const MuonLayerHough::Maximum& ref, const MuonLayerHough::Maximum& ex, bool doparabolic=false);
}
#endif
......@@ -48,7 +48,7 @@ namespace MuonHough {
private:
std::vector< MuonLayerHough* > m_transforms; /// Hough transforms for all regions
int m_sector; /// sector number
//int m_sector; /// sector number
};
/** class managing all Hough transforms in the detector */
......
......@@ -522,5 +522,80 @@ namespace MuonHough {
return std::make_pair<int,int>(floor((zmin-m_descriptor.yMinRange)*m_invbinsize), floor((zmax-m_descriptor.yMinRange)*m_invbinsize));//convert the output to bins
}
float extrapolate(const MuonLayerHough::Maximum& ref, const MuonLayerHough::Maximum& ex, bool doparabolic){
//z is always the precision plane. r is the reference plane, simple and robust
double ref_z = ref.getGlobalZ();
double ref_r = ref.getGlobalR();
double ex_z = ex.getGlobalZ();
double ex_r = ex.getGlobalR();
float theta_ref = ref.getGlobalTheta();
if (!doparabolic || ref_z == 0 || theta_ref == 0){//do linear extrapolation
if (!ex.isEndcap()){//if extrapolate to barell
return ex_z - ex_r / ref_r * ref_z;
}
else{//if extrapolate to endcap
return ex_r - ex_z * ref_r / ref_z;
}
}
else{//do parabolic
float expected = 0;
float extrapolated_diff = 9999;
float tan_theta_ref = tan(theta_ref);
float invtan_theta_ref = 1./tan(theta_ref);
float r_start = ref.hough->m_descriptor.chIndex%2 > 0? 4900.:5200.; //start of barrel B field; values could be further optimized; 5500.:6500.
float z_start = 8500.; //start of endcap B field; used to be 6500; should start at 8500
float z_end = 12500.; //end of endcap B field; used to be 12000; should end at 12500
float r_SL = ref_r + (ex_z-ref_z) * tan_theta_ref;
float z_SL = ref_z + (ex_r-ref_r) * invtan_theta_ref;
//start extrapolation
if (!ref.isEndcap()){//this is starting with barrel chamber; BEE is 6
float rgeo = ref_r*ref_r - r_start*r_start;
float rhoInv = (invtan_theta_ref*ref_r - ref_z)/rgeo;
if (!ex.isEndcap()){//this ends with barrel chamber; BEE is 6
expected = ref_z + (ex_r-ref_r) * invtan_theta_ref + (ex_r-ref_r)*(ex_r-ref_r)*rhoInv;
extrapolated_diff = ex_z - expected;
}//end with barrel to barrel extrapolation
else{//this is a barrel to endcap extrapolation, mostly in the transition region
//recalculate z start
z_start = ref_z + (r_start-ref_r) * invtan_theta_ref + (r_start-ref_r)*(r_start-ref_r)*rhoInv;
float zgeo = ref_z*ref_z - z_start*z_start;
float rho = (tan_theta_ref*ref_z - ref_r)/zgeo;
expected = ref_r + (ex_z-ref_z) * tan_theta_ref + (ex_z-ref_z)*(ex_z-ref_z)*rho;
extrapolated_diff = ex_r - expected;
}
}
else{//this starts with endcap chamber;
//swap the starting position if on the other side
if(tan_theta_ref < 0){
z_start = - z_start;
z_end = - z_end;
}
if (ex.isEndcap()){//extrapolate to endcap
if (std::abs(ref_z) < std::abs(z_end)){//extrapolate from EI or EE, have to use linear
expected = r_SL;
}
else{// from EM or EO or EE
if (std::abs(ex_z) > std::abs(z_start)){//extrapolate to EM or EO
//extrapolate to outer layer, just using theta of the middle measurement; only works if the theta measurement is correct
//can extrapolate with either the outside theta or middle theta; outside theta is better; farther away from the B field
expected = r_SL;
}
else{//to EI
float r_end = ref_r + (z_end-ref_z)*tan_theta_ref;
float zgeo = z_start * z_start - z_end * z_end;
float rhoInv = (r_end - z_end*tan_theta_ref) / zgeo;
float tantheta = tan_theta_ref - 2*(z_end - z_start)*rhoInv;
expected = ex_z*tantheta + (ex_z-z_start)*(ex_z-z_start)*rhoInv;
}
}
extrapolated_diff = ex_r - expected;
}
else{//exrapolate to barrel; again verly likely to be in transition region, complicated B field; just use linear
expected = z_SL;
extrapolated_diff = ex_z - expected;
}
}
return extrapolated_diff;
}//end of parabolic extrapolation
}
}
......@@ -5,6 +5,7 @@
#include "MuonLayerHough/MuonLayerHoughSelector.h"
#include <iostream>
#include <algorithm>
#include <cmath>
namespace MuonHough {
......@@ -27,7 +28,7 @@ namespace MuonHough {
float MuonLayerHoughSelector::getCutValue(float position) const
{
const float pos = fabs(position);
const float pos = std::abs(position);
for (const auto& cut : m_cutValues)
{
if (cut.first < pos) return cut.second;//notice the array is already sorted
......
......@@ -54,7 +54,8 @@ namespace MuonHough {
// if we get to the next layer process the current one and fill the Hough space
if( prevlayer != (*it)->layer ) {
for( int i=0;i<m_nbins;++i ) {
m_histo[i] += layerCounts[i];
if( subtract && -layerCounts[i] >= static_cast<int>(m_histo[i]) ) m_histo[i] = 0;
else m_histo[i] += layerCounts[i];
//if( m_debug && layerCounts[i] != 0 ) std::cout << " filling layer " << prevlayer << " bin " << i << std::endl;
layerCounts[i] = 0; // reset bin
}
......
......@@ -7,7 +7,7 @@
namespace MuonHough {
MuonSectorHough::MuonSectorHough( int sector, const MuonDetectorDescription& regionDescriptions ) : m_sector(sector) {
MuonSectorHough::MuonSectorHough( int sector, const MuonDetectorDescription& regionDescriptions ) {
m_transforms.resize(Muon::MuonStationIndex::sectorLayerHashMax(),0);
// loop over regions and layers of the detector and build the transforms
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment