diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index 15c72ed6..363b7022 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -77,6 +77,8 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Check button if working with the Card Cage at UCSD TGCheckButton* m_UCSDOverride; + //! Check button if weighting X and Y by energy + TGCheckButton* m_WeightedXY; #ifdef ___CLING___ public: diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 21f48291..11582dd1 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -77,6 +77,10 @@ class MModuleDepthCalibration2024 : public MModule //! Get whether the data came from the card cage at UCSD bool GetUCSDOverride() const {return m_UCSDOverride;} + //! Set whether X and Y positions should be calculated by weighting strips by energy + void SetWeightedXY( bool WeightedXY ) {m_WeightedXY = WeightedXY;} + //! Get whether X and Y positions should be calculated by weighting strips by energy + bool GetWeightedXY() const {return m_WeightedXY;} //! Read the XML configuration bool ReadXmlConfiguration(MXmlNode* Node); @@ -95,10 +99,14 @@ class MModuleDepthCalibration2024 : public MModule vector GetDepth(int DetID); //! Retrieve the appropriate CTD values given the DetID and Grade vector GetCTD(int DetID, int Grade); + //! Retrieve the appropriate depth-to-CTD spline given the DetID and Grade + TSpline3* GetSpline(int DetID, int Grade); //! Normal distribution vector norm_pdf(vector x, double mu, double sigma); //! Adds a Depth-to-CTD relation - bool AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap); + bool AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints); + // Calculate the Energy-weighted X and Y strip position + bool CalculateEnergyWeightedPosition(vector LVStrips, vector HVStrips, double& WeightedLVStripID, double& WeightedHVStripID); //! Determine the Grade (geometry of charge sharing) of the Hit int GetHitGrade(MHit* H); //! Load in the specified coefficients file @@ -145,11 +153,13 @@ class MModuleDepthCalibration2024 : public MModule // The CTD Map maps each detector (int) to a 2D array of CTD values. unordered_map>> m_CTDMap; unordered_map> m_DepthGrid; + unordered_map> m_SplineMap; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; // boolean for use with the card cage at UCSD since it tags all events as detector 11 bool m_UCSDOverride; + bool m_WeightedXY; diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index cd6c07ec..e7c97f26 100644 --- a/src/MGUIOptionsDepthCalibration2024.cxx +++ b/src/MGUIOptionsDepthCalibration2024.cxx @@ -80,8 +80,13 @@ void MGUIOptionsDepthCalibration2024::Create() m_UCSDOverride = new TGCheckButton(m_OptionsFrame, "Check this box if you're using the card cage at UCSD", 1); m_UCSDOverride->SetOn(dynamic_cast(m_Module)->GetUCSDOverride()); + TGLayoutHints* Label3Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_UCSDOverride, Label3Layout); + + m_WeightedXY = new TGCheckButton(m_OptionsFrame, "Check this box to weight the X and Y positions by energy deposited.", 1); + m_WeightedXY->SetOn(dynamic_cast(m_Module)->GetWeightedXY()); TGLayoutHints* Label4Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - m_OptionsFrame->AddFrame(m_UCSDOverride, Label4Layout); + m_OptionsFrame->AddFrame(m_WeightedXY, Label4Layout); PostCreate(); } @@ -128,6 +133,7 @@ bool MGUIOptionsDepthCalibration2024::OnApply() dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); dynamic_cast(m_Module)->SetUCSDOverride(m_UCSDOverride->IsOn()); + dynamic_cast(m_Module)->SetWeightedXY(m_WeightedXY->IsOn()); return true; } diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 135db6dc..022ab1c8 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -246,8 +246,8 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) // Calculate the position. If error is thrown, record and no depth. // Take a Hit and separate its activated X- and Y-strips into separate vectors. - std::vector LVStrips; - std::vector HVStrips; + vector LVStrips; + vector HVStrips; for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ MStripHit* SH = H->GetStripHit(j); if( SH->IsLowVoltageStrip() ) LVStrips.push_back(SH); else HVStrips.push_back(SH); @@ -266,6 +266,19 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) int HVStripID = HVSH->GetStripID(); int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; + if (m_WeightedXY==true) { + + // Try to calculate the energy-weighted strip positions. If that doesn't work, place the Hit in the middle of the dominant strips. + double WeightedLVStripID = 0; + double WeightedHVStripID = 0; + bool WeightedPosResult = CalculateEnergyWeightedPosition(LVStrips, HVStrips, WeightedLVStripID, WeightedHVStripID); + if (WeightedPosResult == true) { + LVStripID = WeightedLVStripID; + HVStripID = WeightedHVStripID; + } + + } + // TODO: Calculate X and Y positions more rigorously using charge sharing. // Somewhat confusing notation: HVStrips run parallel to X-axis, so we calculate X position with LVStrips. double Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)LVStripID)); @@ -278,6 +291,9 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) vector* Coeffs = GetPixelCoeffs(PixelCode); + vector CTDVec = GetCTD(DetID, Grade); + vector DepthVec = GetDepth(DetID); + // TODO: For Card Cage, may need to add noise double LVTiming = LVSH->GetTiming(); double HVTiming = HVSH->GetTiming(); @@ -288,34 +304,29 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); ++m_Error1; - } - // If there isn't timing information, set no depth. - // Alex's old comments suggest assigning the event to the middle of the detector and the position resolution to be large. - else if( (LVTiming < 1.0E-6) || (HVTiming < 1.0E-6) ){ + } else if (CTDVec.size() == 0) { + cout << "Empty CTD vector" << endl; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + } else if (DepthVec.size() == 0) { + cout << "Empty Depth vector" << endl; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + } else if ((LVTiming < 1.0E-6) || (HVTiming < 1.0E-6)) { ++m_Error3; H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); - } - - // If there are coefficients and timing information is loaded, try calculating the CTD and depth - else { - - vector ctdvec = GetCTD(DetID, Grade); - vector depthvec = GetDepth(DetID); - - if ( ctdvec.size() == 0){ - cout << "Empty CTD vector" << endl; - H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); - } + } else { + + // If there are coefficients and timing information is loaded, try calculating the CTD and depth double CTD = (HVTiming - LVTiming); // Confirmed that this matches SP's python code. CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset - double Xmin = * std::min_element(ctdvec.begin(), ctdvec.end()); - double Xmax = * std::max_element(ctdvec.begin(), ctdvec.end()); + double Xmin = * std::min_element(CTDVec.begin(), CTDVec.end()); + double Xmax = * std::max_element(CTDVec.begin(), CTDVec.end()); double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); @@ -327,32 +338,34 @@ bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) } // If the CTD is in range, calculate the depth + // Rather than plugging CTD into a spline to get depth, use the depth-CTD relation to calculate a probability-weighted depth value. + // This way we can avoid problems like non-monotonicity or assigning depth to events "outside" the detector + // Note that this requires that we don't massively overestimate the timing noise else { - // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in depthvec + // Calculate the probability given timing noise of CTD_s corresponding to the values of depth in DepthVec // Utlize symmetry of the normal distribution. - vector prob_dist = norm_pdf(ctdvec, CTD_s, noise/2.355); + vector prob_dist = norm_pdf(CTDVec, CTD_s, noise/2.355); // Weight the depth by probability double prob_sum = 0.0; for( unsigned int k=0; k < prob_dist.size(); ++k ){ prob_sum += prob_dist[k]; } - //double prob_sum = std::accumulate(prob_dist.begin(), prob_dist.end(), 0); double weighted_depth = 0.0; - for( unsigned int k = 0; k < depthvec.size(); ++k ){ - weighted_depth += prob_dist[k]*depthvec[k]; + for( unsigned int k = 0; k < DepthVec.size(); ++k ){ + weighted_depth += prob_dist[k]*DepthVec[k]; } // Calculate the expectation value of the depth double mean_depth = weighted_depth/prob_sum; // Calculate the standard deviation of the depth double depth_var = 0.0; - for( unsigned int k=0; kIsStripPairingIncomplete()==false) { @@ -490,7 +503,6 @@ vector MModuleDepthCalibration2024::norm_pdf(vector x, double mu vector result; for( unsigned int i=0; i MModuleDepthCalibration2024::norm_pdf(vector x, double mu bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) { - //when invert flag is set to true, the splines returned are CTD->Depth - // Previously saved cathode and anode timing in addition to CTD. This may be redundant, commenting out for now. // Input spline files should have the following format: // ### DetID, HV, Temperature, Photopeak Energy (TODO: More? Fewer?) // depth, ctd0, ctd1, ctd2.... (Basically, allow for CTDs for different subpixel regions) // '' '' '' - MFile F; + MFile F; if( F.Open(FName) == false ){ cout << "ERROR in MModuleDepthCalibration2024::LoadSplinesFile: failed to open splines file." << endl; return false; } - // vector depthvec, ctdvec, anovec, catvec; - vector depthvec; - vector> ctdarr; - for( unsigned int i=0; i < 5; ++i ){ - vector temp_vec; - ctdarr.push_back(temp_vec); - } + vector DepthVec; + vector> CTDArr; bool Result = true; MString line; int DetID = 0; @@ -524,39 +529,71 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) if (line.BeginsWith("#")) { // If we've reached a new ctd spline then record the previous one in the m_SplineMaps and start a new one. vector tokens = line.Tokenize(" "); - if (depthvec.size() > 0) { - Result &= AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); - } - depthvec.clear(); ctdarr.clear(); - for (unsigned int i=0; i < 5; ++i) { - vector temp_vec; - ctdarr.push_back(temp_vec); + if (DepthVec.size() > 0) { + Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); } + DepthVec.clear(); CTDArr.clear(); DetID = tokens[1].ToInt(); } else { vector tokens = line.Tokenize(","); - depthvec.push_back(tokens[0].ToDouble()); - + DepthVec.push_back(tokens[0].ToDouble()); // Multiple CTDs allowed. for (unsigned int i = 0; i < (tokens.size() - 1); ++i) { - ctdarr[i].push_back(tokens[1+i].ToDouble()); - } - // Fill in the higher grades with the GRADE=0 CTD if there are none listed in the file. - for (unsigned int i=tokens.size()-1; i<5; ++i) { - ctdarr[i].push_back(tokens[1].ToDouble()); + while (i>=CTDArr.size()) { + vector TempVec; + CTDArr.push_back(TempVec); + } + CTDArr[i].push_back(tokens[1+i].ToDouble()); } } } } //make last spline - if (depthvec.size() > 0) { - Result &= AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); + if (DepthVec.size() > 0) { + Result &= AddDepthCTD(DepthVec, CTDArr, DetID, m_DepthGrid, m_CTDMap, m_SplineMap, 1000); } return Result; } +bool MModuleDepthCalibration2024::CalculateEnergyWeightedPosition(vector LVStrips, vector HVStrips, double& WeightedLVStripID, double& WeightedHVStripID) { + // Determine the weighted strip ID positions based on charge sharing. + WeightedLVStripID = 0; + WeightedHVStripID = 0; + double TotalLVEnergy = 0; + double TotalHVEnergy = 0; + + // If one side or another doesn't have strips, abort. + if ( (LVStrips.size() == 0) || (HVStrips.size()==0) ) { + return false; + } + + // Calculate the weighted values based on energy per strip. + for ( unsigned int i=0; iGetEnergy(); + int LVSHID = LVStrips[i]->GetStripID(); + TotalLVEnergy += LVSHEnergy; + WeightedLVStripID += LVSHEnergy*LVSHID; + } + for ( unsigned int i=0; iGetEnergy(); + int HVSHID = HVStrips[i]->GetStripID(); + TotalHVEnergy += HVSHEnergy; + WeightedHVStripID += HVSHEnergy*HVSHID; + } + + // If one side or another doesn't have energy, abort. + if ( (TotalLVEnergy == 0) || (TotalHVEnergy == 0) ) { + return false; + } + + WeightedLVStripID /= TotalLVEnergy; + WeightedHVStripID /= TotalHVEnergy; + + return true; +} + int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ // Function for choosing which Depth-to-CTD relation to use for a given event. // At time of writing, intention is to choose a CTD based on sub-pixel region determined via charge sharing (Event "grade"). @@ -611,7 +648,6 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ } } else{ - // cout << "ERROR in MModuleDepthCalibration2024: HIT with no strip hits on one side" << endl; return -1; } @@ -672,44 +708,101 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ return return_value; } -bool MModuleDepthCalibration2024::AddDepthCTD(vector depthvec, vector> ctdarr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap){ +bool MModuleDepthCalibration2024::AddDepthCTD(vector Depth, vector> CTDArr, int DetID, unordered_map>& DepthGrid, unordered_map>>& CTDMap, unordered_map>& SplineMap, unsigned int NPoints) +{ // Saves a CTD array, basically allowing for multiple CTDs as a function of depth - // depthvec: list of simulated depth values - // ctdarr: vector of vectors of simulated CTD values. Each vector of CTDs must be the same length as depthvec + // Depth: list of simulated depth values + // CTDArr: vector of vectors of simulated CTD values. Each vector of CTDs must be the same length as Depth // DetID: An integer which specifies which detector. // CTDMap: unordered map into which the array of CTDs should be placed + // SplineMap: unordered map to store splines + // NPoints: number of points in the depth/CTD grid to be produced // TODO: Possible energy dependence of CTD? - // TODO: Depth values need to be evenly spaced. Check this when reading the files in. // Check to make sure things look right. // First check that the CTDs all have the right length. - for (unsigned int i = 0; i < ctdarr.size(); ++i) { - if ((ctdarr[i].size() != depthvec.size()) && (ctdarr[i].size() > 0)) { + + for (unsigned int i = 0; i < CTDArr.size(); ++i) { + if ((CTDArr[i].size() != Depth.size()) && (CTDArr[i].size() > 0)) { cout<<"ERROR in MModuleDepthCalibration2024::AddDepthCTD: The number of values in the CTD list is not equal to the number of depth values."< 0.01) { + // Check that the geometry file and the depth-CTD curve match within 0.1mm tolerance + double MaxDepth = * std::max_element(Depth.begin(), Depth.end()); + double MinDepth = * std::min_element(Depth.begin(), Depth.end()); + if (fabs((MaxDepth-MinDepth) - m_Thicknesses[DetID]) > 0.01) { cout<<"ERROR in MModuleDepthCalibration2024::AddDepthCTD: The thickness of detector "<0) { + cout<<"MModuleDepthCalibration2024::AddDepthCTD: Splines already added for DetID "< TempVec; + SplineMap[DetID] = TempVec; + } + + // Add a new point to the depth vector on the bottom and top of the detector to make sure we're sampling the full depth + double dx_low = Depth[1] - Depth[0]; + double newx_low = Depth[0] - (dx_low/2.0); + Depth.insert(Depth.begin(), newx_low); + + size_t N = Depth.size(); + double dx_high = Depth[N-1] - Depth[N-2]; + double newx_high = Depth[N-1] + (dx_high/2.0); + Depth.push_back(newx_high); + + // If the depth vector is more dense than NPoints, use the length of the input vector + if (NPoints < N+1) { + NPoints = N+1; + } + + // Reconstitute the depth vector to ensure that the depth is evenly sampled + vector NewDepth; + for (unsigned int i=0; i> NewCTDArr; + for (unsigned int i=0; i CTD = CTDArr[i]; + vector NewCTD; - //Now make sure the values for the depth start with 0.0. - if (mindepth != 0.0) { - cout<<"MModuleDepthCalibration2024::AddDepthCTD: The minimum depth is not zero. Editing the depth vector."<Eval(NewDepth[d])); + } + NewCTDArr.push_back(NewCTD); } - CTDMap[DetID] = ctdarr; - DepthGrid[DetID] = depthvec; + CTDMap[DetID] = NewCTDArr; return true; } @@ -729,7 +822,6 @@ vector MModuleDepthCalibration2024::GetCTD(int DetID, int Grade) return (m_CTDMap[DetID][Grade]); } else { - cout << "MModuleDepthCalibration2024::GetCTD: No CTD map is loaded for Grade " << Grade << ". Returning Grade 0 CTD." << endl; return (m_CTDMap[DetID][0]); } } else { @@ -757,6 +849,31 @@ vector MModuleDepthCalibration2024::GetDepth(int DetID) } } + +TSpline3* MModuleDepthCalibration2024::GetSpline(int DetID, int Grade) +{ + // Retrieves the appropriate depth->CTD spline given the Detector ID and Event Grade passed + + if( !m_SplinesFileIsLoaded ){ + cout << "MModuleDepthCalibration2024::GetSpline: cannot return Depth to CTD spline because the file was not loaded." << endl; + return nullptr; + } + + // If there is a spline for the given detector, return it. + // If the Grade is larger than the number of splines stored, then just return Grade 0 spline. + if( m_SplineMap.count(DetID) > 0 ){ + if ( ((int)m_SplineMap[DetID].size()) > Grade) { + return (m_SplineMap[DetID][Grade]); + } + else { + return (m_SplineMap[DetID][0]); + } + } else { + cout << "MModuleDepthCalibration2024::GetSpline: No spline is loaded for Det " << DetID << "." << endl; + return nullptr; + } +} + void MModuleDepthCalibration2024::ShowOptionsGUI() { // Show the options GUI - or do nothing @@ -785,6 +902,11 @@ bool MModuleDepthCalibration2024::ReadXmlConfiguration(MXmlNode* Node) m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt(); } + MXmlNode* WeightedXYNode = Node->GetNode("WeightedXY"); + if( WeightedXYNode != NULL ){ + m_WeightedXY = (bool) WeightedXYNode->GetValueAsInt(); + } + return true; } @@ -799,6 +921,7 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration() new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride); + new MXmlNode(Node, "WeightedXY",(unsigned int) m_WeightedXY); return Node; } @@ -820,11 +943,19 @@ void MModuleDepthCalibration2024::Finalize() cout << "Number of hits with no strip hits on one or both sides: " << m_ErrorSH << endl; cout << "Number of hits with null strip hits: " << m_ErrorNullSH << endl; cout << "Number of hits 0 energy on a strip hit: " << m_ErrorNoE << endl; - /* - TFile* rootF = new TFile("EHist.root","recreate"); - rootF->WriteTObject( EHist ); - rootF->Close(); - */ + + // Clean up maps and vectors + m_Coeffs.clear(); + m_Thicknesses.clear(); + m_NXStrips.clear(); + m_NYStrips.clear(); + m_XPitches.clear(); + m_YPitches.clear(); + m_Detectors.clear(); + m_CTDMap.clear(); + m_DepthGrid.clear(); + m_SplineMap.clear(); + m_DetectorIDs.clear(); }