diff --git a/Makefile b/Makefile index d9924ad1..a260b9c0 100755 --- a/Makefile +++ b/Makefile @@ -38,15 +38,33 @@ MAKEFLAGS += --no-builtin-rules #.NOTPARALLEL: megalib .SILENT: + +#---------------------------------------------------------------- +# External libraries +# + +H5CXXFLAGS = +H5LIBS = +ifeq ("$(shell pkg-config --exists hdf5 1>&2 2> /dev/null; echo $$?)", "0") + H5CXXFLAGS += $(shell pkg-config --cflags hdf5) + H5LIBS += $(shell pkg-config --libs hdf5) + H5LIBS += -lhdf5_cpp +else + $(error "Unable to find HDF5 headers and libraries") +endif + + #---------------------------------------------------------------- # Definitions based on architecture and user options # CMD="" -CXXFLAGS += -I$(IN) -I$(MEGALIB)/include -I/opt/local/include +CXXFLAGS += -I$(IN) -I$(MEGALIB)/include -I/opt/local/include $(H5CXXFLAGS) # Comment this line out if you want to accept warnings #CXXFLAGS += -Werror -Wno-unused-variable +LIBS += $(H5LIBS) + # Names of the program NUCLEARIZER_PRG = $(BN)/nuclearizer NUCLEARIZER_CXX_MAIN = src/MNuclearizerMain.cxx @@ -62,6 +80,7 @@ $(LB)/MAspectReconstruction.o \ $(LB)/MHit.o \ $(LB)/MTimeAndCoordinate.o \ $(LB)/MStripHit.o \ +$(LB)/MStripMap.o \ $(LB)/MGuardringHit.o \ $(LB)/MDetectorEffectsEngineBalloon.o \ $(LB)/MModuleLoaderSimulationsBalloon.o \ @@ -71,6 +90,8 @@ $(LB)/MGUIOptionsLoaderSimulations.o \ $(LB)/MModuleLoaderMeasurements.o \ $(LB)/MModuleLoaderMeasurementsROA.o \ $(LB)/MGUIOptionsLoaderMeasurements.o \ +$(LB)/MModuleLoaderMeasurementsHDF.o \ +$(LB)/MGUIOptionsLoaderMeasurementsHDF.o \ $(LB)/MBinaryFlightDataParser.o \ $(LB)/MModuleReceiverBalloon.o \ $(LB)/MGUIOptionsReceiverBalloon.o \ @@ -116,6 +137,11 @@ $(LB)/MModuleDiagnostics.o \ $(LB)/MModuleDiagnosticsEnergyPerStrip.o \ $(LB)/MGUIExpoDiagnosticsEnergyPerStrip.o \ $(LB)/MGUIExpoDiagnostics.o \ +$(LB)/MModuleTACcut.o \ +$(LB)/MGUIExpoTACcut.o \ +$(LB)/MGUIOptionsTACcut.o \ +$(LB)/MModuleNearestNeighbor.o \ + @@ -190,6 +216,7 @@ $(FRETALON_DEP_FILES): $(LB)/%.d: $(FRETALON_DIR)/src/%.cxx $(FRETALON_LIBS): $(LB)/%.o: $(FRETALON_DIR)/src/%.cxx $(FRETALON_DIR)/inc/%.h $(LB)/%.d @echo "Compiling $(subst $(FRETALON_DIR)/src/,,$<) ..." @$(CXX) $(CXXFLAGS) -c $< -o $@ + @echo "$(CXX) $(CXXFLAGS) -c $< -o $@" $(NUCLEARIZER_DEP_FILES): $(LB)/%.d: src/%.cxx @echo "Creating dependencies for $(subst src/,,$<) ..." @@ -203,7 +230,7 @@ $(NUCLEARIZER_DICT): $(FRETALON_H_FILES) $(NUCLEARIZER_H_FILES) @echo "Generating LinkDef ..." @$(BN)/generatelinkdef -o $(NUCLEARIZER_LINKDEF) -i $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) @echo "Generating dictionary ..." - @rootcling -f $@ -I$(IN) -I$(MEGALIB)/include -D___CLING___ -rmf $(NUCLEARIZER_ROOTMAP) -s libNuclearizer -c $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) $(NUCLEARIZER_LINKDEF) + @rootcling -f $@ -I$(IN) -I$(MEGALIB)/include $(H5CXXFLAGS) -D___CLING___ -rmf $(NUCLEARIZER_ROOTMAP) -s libNuclearizer -c $(NUCLEARIZER_H_FILES) $(FRETALON_H_FILES) $(NUCLEARIZER_LINKDEF) @mv $(NUCLEARIZER_ROOTPCM) $(LB) $(NUCLEARIZER_DICT_LIB): $(NUCLEARIZER_DICT) @@ -218,6 +245,7 @@ $(NUCLEARIZER_PRG): $(NUCLEARIZER_SHARED_LIB) $(NUCLEARIZER_CXX_MAIN) @echo "Linking and compiling $(subst $(BN)/,,$(NUCLEARIZER_PRG)) ... Please stand by ... " @$(CXX) $(CXXFLAGS) $(LDFLAGS) $(NUCLEARIZER_CXX_MAIN) $(NUCLEARIZER_SHARED_LIB) $(ALLLIBS) $(GLIBS) $(LIBS) -o $(NUCLEARIZER_PRG) + ifneq ($(MAKECMDGOALS),clean) -include $(NUCLEARIZER_DEP_FILES) -include $(FRETALON_DEP_FILES) diff --git a/apps/TrappingCorrection.cxx b/apps/TrappingCorrection.cxx new file mode 100644 index 00000000..306a4c04 --- /dev/null +++ b/apps/TrappingCorrection.cxx @@ -0,0 +1,762 @@ +/* + * TrappingCorrection.cxx + * + * + * Copyright (C) by Sean Pike. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Sean Pike. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + +// Standard +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +// ROOT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace ROOT::Minuit2; + +// MEGAlib +#include "MGlobal.h" +#include "MFile.h" +#include "MReadOutElementDoubleStrip.h" +#include "MFileReadOuts.h" +#include "MReadOutAssembly.h" +#include "MStripHit.h" +#include "MReadOutSequence.h" +#include "MSupervisor.h" +#include "MModuleLoaderMeasurements.h" +#include "MModuleLoaderMeasurementsROA.h" +#include "MModuleLoaderMeasurementsHDF.h" +#include "MModuleEnergyCalibrationUniversal.h" +#include "MModuleEventFilter.h" +#include "MModuleStripPairingGreedy.h" +#include "MModuleStripPairingChiSquare.h" +#include "MModuleTACcut.h" +#include "MAssembly.h" + + +int g_HistBins = 75; +double g_MinCTD = -250; +double g_MaxCTD = 250; +double g_MinRatio = 0.96; +double g_MaxRatio = 1.04; + + +//////////////////////////////////////////////////////////////////////////////// + + +class SymmetryFCN : public FCNBase +{ +public: + //! Operator which returns the symmetry of m_CTDHistogram given the parameters passed + double operator()(vector const &v) const override; + double Up() const override { return 1; } + + void AddCTD(double CTD){ m_CTD.push_back(CTD); } + void AddHVEnergy(double HVEnergy, double HVEnergyResolution){ m_HVEnergy.push_back(HVEnergy); m_HVEnergyResolution.push_back(HVEnergyResolution); } + void AddLVEnergy(double LVEnergy, double LVEnergyResolution){ m_LVEnergy.push_back(LVEnergy); m_LVEnergyResolution.push_back(LVEnergyResolution); } + + vector GetCTD(){ return m_CTD; } + vector GetHVEnergy(){ return m_HVEnergy; } + vector GetLVEnergy(){ return m_LVEnergy; } + vector GetHVEnergyResolution(){ return m_HVEnergyResolution; } + vector GetLVEnergyResolution(){ return m_LVEnergyResolution; } + + void SetCTD(vector CTD){ m_CTD = CTD; } + void SetHVEnergy(vector HVEnergy, vector HVEnergyResolution){ m_HVEnergy = HVEnergy; m_HVEnergyResolution = HVEnergyResolution; } + void SetLVEnergy(vector LVEnergy, vector LVEnergyResolution){ m_LVEnergy = LVEnergy; m_LVEnergyResolution = LVEnergyResolution; } + + //! The measured CTD and HV/LV energies + vector m_CTD; + vector m_HVEnergy; + vector m_LVEnergy; + vector m_HVEnergyResolution; + vector m_LVEnergyResolution; + +}; + + +//////////////////////////////////////////////////////////////////////////////// + + +class ChiSquaredFCN : public SymmetryFCN +{ +public: + //! Operator which returns the symmetry of m_CTDHistogram given the parameters passed + double operator()(vector const &v) const override; + double Up() const override { return 1; } + +}; + + +//////////////////////////////////////////////////////////////////////////////// + + + +//! A standalone program based on MEGAlib and ROOT +class TrappingCorrection +{ +public: + //! Default constructor + TrappingCorrection(); + //! Default destructor + ~TrappingCorrection(); + + //! Parse the command line + bool ParseCommandLine(int argc, char** argv); + //! Analyze what ever needs to be analyzed... + bool Analyze(); + //! Interrupt the analysis + void Interrupt() { m_Interrupt = true; } + + void dummy_func() { return; } + + MStripHit* GetDominantStrip(vector& Strips, double& EnergyFraction); + +private: + //! True, if the analysis needs to be interrupted + bool m_Interrupt; + //! The input file name + MString m_FileName; + MString m_EcalFile; + MString m_TACCalFile; + MString m_TACCutFile; + MString m_StripMapFile; + //! output file names + MString m_OutFile; + //! option to do a pixel-by-pixel calibration (instead of detector-by-detector) + bool m_PixelCorrect; + bool m_GreedyPairing; + bool m_ExcludeNN; + + double m_MinEnergy; + double m_MaxEnergy; + +}; + +//////////////////////////////////////////////////////////////////////////////// + + +double SymmetryFCN::operator()(vector const &v) const +{ + double HVSlope = v[0]; + double LVSlope = v[1]; + + char name[64]; sprintf(name,"name"); + int HistBins = g_HistBins; + if (HistBins%2 == 0) { + HistBins += 1; + } + TH2D CorrectedHistogram(name, name, HistBins, g_MinCTD, g_MaxCTD, HistBins, g_MinRatio, g_MaxRatio); + + for (unsigned int i = 0; i < m_CTD.size(); ++i) { + double CTDHVShift = m_CTD[i] + g_MaxCTD; + double CTDLVShift = m_CTD[i] + g_MinCTD; + // Correct the HV and LV energies by dividing by the CCE. DeltaCCE is defined as a linear function with units percentage energy lost. + double CorrectedHVEnergy = m_HVEnergy[i]/(1 - (HVSlope*CTDHVShift)/100); + double CorrectedLVEnergy = m_LVEnergy[i]/(1 - (LVSlope*CTDLVShift)/100); + CorrectedHistogram.Fill(m_CTD[i], CorrectedHVEnergy/CorrectedLVEnergy); + } + + vector> BinValues; + vector> ReflectedBinValues; + double Asymmetry = 0; + + for (unsigned int y = 0; y < CorrectedHistogram.GetNbinsY(); ++y) { + + vector XValues; + + for (unsigned int x = 0; x < CorrectedHistogram.GetNbinsX(); ++x) { + XValues.push_back(CorrectedHistogram.GetBinContent(x,y)); + } + + // BinValues.push_back(XValues); + vector ReflectedXValues = XValues; + reverse(ReflectedXValues.begin(), ReflectedXValues.end()); + + for (unsigned int x = 0; x < XValues.size(); ++x) { + Asymmetry += pow(XValues[x] - ReflectedXValues[x], 2)/(XValues[x] + ReflectedXValues[x]); + } + } + + return Asymmetry*Asymmetry; + +} + + +//////////////////////////////////////////////////////////////////////////////// + + +double ChiSquaredFCN::operator()(vector const &v) const +{ + double HVLVFactor = v[0]; + double HVSlope = v[1]; + double LVSlope = v[2]; + + double ChiSquare = 0; + + for (unsigned int i = 0; i < m_CTD.size(); ++i) { + double CTDHVShift = m_CTD[i] + g_MaxCTD; + double CTDLVShift = m_CTD[i] + g_MinCTD; + // Correct the HV and LV energies by dividing by the CCE. DeltaCCE is defined as a linear function with units percentage energy lost. + double CorrectedHVEnergy = m_HVEnergy[i]/((1 - (HVSlope*CTDHVShift)/100)*HVLVFactor); + double CorrectedLVEnergy = m_LVEnergy[i]/(1 - (LVSlope*CTDLVShift)/100); + + ChiSquare += pow(CorrectedHVEnergy - CorrectedLVEnergy, 2)/(m_HVEnergyResolution[i] + m_LVEnergyResolution[i]); + } + + // ChiSquare /= m_CTD.size(); + + return ChiSquare; + +} + + +//! Default constructor +TrappingCorrection::TrappingCorrection() : m_Interrupt(false) +{ + gStyle->SetPalette(1, 0); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Default destructor +TrappingCorrection::~TrappingCorrection() +{ + // Intentionally left blank +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Parse the command line +bool TrappingCorrection::ParseCommandLine(int argc, char** argv) +{ + ostringstream Usage; + Usage<"< i+1) && (argv[i+1][0] != '-' || isalpha(argv[i+1][1]) == 0))){ + cout<<"Error: Option "< FileNames; + + if ((m_FileName.GetSubString(m_FileName.Length() - 4)) == "hdf5") { + FileNames.push_back(m_FileName); + } else if ((m_FileName.GetSubString(m_FileName.Length() - 3)) == "txt") { + MFile F; + if (F.Open(m_FileName)==false) { + cout<<"Error: Failed to open input file."< Histograms; + map FCNs; + + for (unsigned int f = 0; fSetFileNameStripMap(m_StripMapFile); + Loader->SetFileName(InFile); + S->SetModule(Loader, MNumber); + ++MNumber; + + cout<<"Creating TAC calibrator"<SetTACCalFileName(m_TACCalFile); + TACCalibrator->SetTACCutFileName(m_TACCutFile); + S->SetModule(TACCalibrator, MNumber); + ++MNumber; + + cout<<"Creating energy calibrator"<SetFileName(m_EcalFile); + EnergyCalibrator->EnablePreampTempCorrection(false); + S->SetModule(EnergyCalibrator, MNumber); + ++MNumber; + + cout<<"Creating Event filter"<SetMinimumLVStrips(1); + EventFilter->SetMaximumLVStrips(3); + EventFilter->SetMinimumHVStrips(1); + EventFilter->SetMaximumHVStrips(3); + EventFilter->SetMinimumHits(0); + EventFilter->SetMaximumHits(100); + EventFilter->SetMinimumTotalEnergy(m_MinEnergy); + EventFilter->SetMaximumTotalEnergy(m_MaxEnergy); + S->SetModule(EventFilter, MNumber); + ++MNumber; + + cout<<"Creating strip pairing"<SetModule(Pairing, MNumber); + + cout<<"Initializing Loader"<Initialize() == false) return false; + cout<<"Initializing TAC calibrator"<Initialize() == false) return false; + cout<<"Initializing Energy calibrator"<Initialize() == false) return false; + cout<<"Initializing Event filter"<Initialize() == false) return false; + cout<<"Initializing Pairing"<Initialize() == false) return false; + + bool IsFinished = false; + MReadOutAssembly* Event = new MReadOutAssembly(); + + cout<<"Analyzing..."<Clear(); + + if (Loader->IsReady()){ + + Loader->AnalyzeEvent(Event); + TACCalibrator->AnalyzeEvent(Event); + EnergyCalibrator->AnalyzeEvent(Event); + bool Unfiltered = EventFilter->AnalyzeEvent(Event); + + if (Unfiltered == true) { + + Pairing->AnalyzeEvent(Event); + + if (Event->HasAnalysisProgress(MAssembly::c_StripPairing) == true) { + for (unsigned int h = 0; h < Event->GetNHits(); ++h) { + double HVEnergy = 0.0; + double LVEnergy = 0.0; + double HVEnergyResolution = 0.0; + double LVEnergyResolution = 0.0; + vector HVStrips; + vector LVStrips; + + MHit* H = Event->GetHit(h); + + int DetID = H->GetStripHit(0)->GetDetectorID(); + TH2D* Hist = Histograms[DetID]; + SymmetryFCN* FCN = FCNs[DetID]; + + if (Hist == nullptr) { + char name[64]; sprintf(name,"Detector %d (Uncorrected)",DetID); + Hist = new TH2D(name, name, g_HistBins, g_MinCTD, g_MaxCTD, g_HistBins, g_MinRatio, g_MaxRatio); + Hist->SetXTitle("CTD (ns)"); + Hist->SetYTitle("HV/LV Energy Ratio"); + Histograms[DetID] = Hist; + } + + if (FCN == nullptr) { + FCN = new SymmetryFCN(); + FCNs[DetID] = FCN; + } + + for (unsigned int sh = 0; sh < H->GetNStripHits(); ++sh) { + MStripHit* SH = H->GetStripHit(sh); + + if ((m_ExcludeNN==false) || ((m_ExcludeNN==true) && (SH->IsNearestNeighbor()==false))) { + if (SH->IsLowVoltageStrip()==true) { + LVEnergy += SH->GetEnergy(); + LVEnergyResolution += (SH->GetEnergyResolution())*(SH->GetEnergyResolution()); + LVStrips.push_back(SH); + } else { + HVEnergy += SH->GetEnergy(); + HVEnergyResolution += (SH->GetEnergyResolution())*(SH->GetEnergyResolution()); + HVStrips.push_back(SH); + } + } + } + + if ((HVStrips.size()>0) && (LVStrips.size()>0)) { + double HVEnergyFraction = 0; + double LVEnergyFraction = 0; + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + double EnergyFraction = HVEnergy/LVEnergy; + if ((LVSH->HasCalibratedTiming()==true) && (HVSH->HasCalibratedTiming()==true)) { + double CTD = LVSH->GetTiming() - HVSH->GetTiming(); + Hist->Fill(CTD, EnergyFraction); + FCN->AddCTD(CTD); + FCN->AddHVEnergy(HVEnergy, HVEnergyResolution); + FCN->AddLVEnergy(LVEnergy, LVEnergyResolution); + } + } + } + } + } + } + IsFinished = Loader->IsFinished(); + } + } + + //setup output file + ofstream OutputCalFile; + OutputCalFile.open(m_OutFile+MString("_parameters.txt")); + OutputCalFile<<"Det"<<'\t'<<"HV Slope"<<'\t'<<"LV Slope"<<'\t'<<"HV/LV Factor"<SetLogz(); + C->cd(); + H.second->Draw("colz"); + + H.second->Write(); + f.Close(); + + } + + for (auto F: FCNs) { + + int DetID = F.first; + MnUserParameters* InitialStateSym = new MnUserParameters(); + InitialStateSym->Add("HVSlope", 2e-3, 1e-4, 0, 3e-1); + InitialStateSym->Add("LVSlope", 0, 1e-4, -3e-2, 0); + + InitialStateSym->Fix("LVSlope"); + + MnMigrad migradSym(*F.second, *InitialStateSym); + // Minimize + FunctionMinimum MinimumSym = migradSym(); + + MnUserParameters ParametersSym = MinimumSym.UserParameters(); + double HVSlope = ParametersSym.Value("HVSlope"); + double LVSlope = ParametersSym.Value("LVSlope"); + + // output + cout<Add("HVLVFactor", 1.0, 0.01, 0.95, 1.05); + // InitialStateChi->Add("HVSlope", ParametersSym.Value("HVSlope"), ParametersSym.Error("HVSlope")); + // InitialStateChi->Add("LVSlope", ParametersSym.Value("LVSlope"), ParametersSym.Error("LVSlope")); + + // InitialStateChi->Fix("HVSlope"); + // InitialStateChi->Fix("LVSlope"); + + + // ChiSquaredFCN* ChiSquaredF = new ChiSquaredFCN(); + // ChiSquaredF->SetCTD(F.second->GetCTD()); + // ChiSquaredF->SetHVEnergy(F.second->GetHVEnergy(), F.second->GetHVEnergyResolution()); + // ChiSquaredF->SetLVEnergy(F.second->GetLVEnergy(), F.second->GetLVEnergyResolution()); + // MnMigrad migradChi(*ChiSquaredF, *InitialStateChi); + // // Minimize + // FunctionMinimum MinimumChi = migradChi(); + + // MnUserParameters ParametersChi = MinimumChi.UserParameters(); + // double HVLVFactor = ParametersChi.Value("HVLVFactor"); + double HVLVFactor = 1; + // // output + // cout<SetXTitle("CTD (ns)"); + Hist->SetYTitle("HV/LV Energy Ratio"); + + vector CTDList = F.second->GetCTD(); + vector HVEnergyList = F.second->GetHVEnergy(); + vector LVEnergyList = F.second->GetLVEnergy(); + for (unsigned int i=0; iFill(CTDList[i], CorrectedHVEnergy/CorrectedLVEnergy); + } + + TCanvas* C = new TCanvas(); + C->SetLogz(); + C->cd(); + Hist->Draw("colz"); + + TFile f(m_OutFile+MString("_Det")+DetID+MString("_Hist_Corr.root"),"recreate"); + Hist->Write(); + f.Close(); + + OutputCalFile<& Strips, double& EnergyFraction) +{ + double MaxEnergy = -numeric_limits::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one + double TotalEnergy = 0.0; + MStripHit* MaxStrip = nullptr; + + // Iterate through strip hits and get the strip with highest energy + for (const auto SH : Strips) { + double Energy = SH->GetEnergy(); + TotalEnergy += Energy; + if (Energy > MaxEnergy) { + MaxStrip = SH; + MaxEnergy = Energy; + } + } + if (TotalEnergy == 0) { + EnergyFraction = 0; + } else { + EnergyFraction = MaxEnergy/TotalEnergy; + } + return MaxStrip; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Called when an interrupt signal is flagged +//! All catched signals lead to a well defined exit of the program +void CatchSignal(int a) +{ + if (g_Prg != 0 && g_NInterruptCatches-- > 0) { + cout<<"Catched signal Ctrl-C (ID="<Interrupt(); + } else { + abort(); + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Main program +int main(int argc, char** argv) +{ + // Catch a user interupt for graceful shutdown + signal(SIGINT, CatchSignal); + + // Initialize global MEGALIB variables, especially mgui, etc. + MGlobal::Initialize("Standalone", "a standalone example program"); + + TApplication TrappingCorrectionApp("TrappingCorrectionApp", 0, 0); + + g_Prg = new TrappingCorrection(); + + if (g_Prg->ParseCommandLine(argc, argv) == false) { + cerr<<"Error during parsing of command line!"<Analyze() == false) { + cerr<<"Error during analysis!"< +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +// ROOT +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib +#include "MGlobal.h" +#include "MFile.h" +#include "MReadOutElementDoubleStrip.h" +#include "MFileReadOuts.h" +#include "MReadOutAssembly.h" +#include "MStripHit.h" +#include "MReadOutSequence.h" +#include "MSupervisor.h" +#include "MModuleLoaderMeasurementsHDF.h" +#include "MModuleEnergyCalibrationUniversal.h" +#include "MModuleEventFilter.h" +#include "MModuleStripPairingGreedy.h" +#include "MModuleStripPairingChiSquare.h" +#include "MModuleTACcut.h" +#include "MAssembly.h" + + +double g_MinCTD = -400; +double g_MaxCTD = 400; +int g_MinCounts = 250; + +int g_HVStrips = 64; +int g_LVStrips = 64; + +double g_AmPhotopeak = 59.54; + +//////////////////////////////////////////////////////////////////////////////// + + +//! A standalone program based on MEGAlib and ROOT +class TrappingCorrectionAm241 +{ +public: + //! Default constructor + TrappingCorrectionAm241(); + //! Default destructor + ~TrappingCorrectionAm241(); + + //! Parse the command line + bool ParseCommandLine(int argc, char** argv); + //! Analyze what ever needs to be analyzed... + bool Analyze(); + //! Interrupt the analysis + void Interrupt() { m_Interrupt = true; } + + //! Produce functions for fitting + TF1* GenerateCTDFunction(double CTDFitMin, double CTDFitMax, double CTDGuess, double FlipSwitch); + TF1* GeneratePhotopeakFunction(); + + MStripHit* GetDominantStrip(vector& Strips, double& EnergyFraction); + +private: + //! True, if the analysis needs to be interrupted + bool m_Interrupt; + //! The input file name + MString m_HVFileName; + MString m_LVFileName; + MString m_EcalFile; + MString m_TACCalFile; + MString m_TACCutFile; + MString m_StripMapFile; + //! output file names + MString m_OutFile; + //! option to do a pixel-by-pixel calibration (instead of detector-by-detector) + bool m_PixelCorrect; + bool m_GreedyPairing; + bool m_ExcludeNN; + bool m_ContinueHDF5; + + double m_MinEnergy; + double m_MaxEnergy; + +}; + +//////////////////////////////////////////////////////////////////////////////// + + +//! Default constructor +TrappingCorrectionAm241::TrappingCorrectionAm241() : m_Interrupt(false) +{ + gStyle->SetPalette(1, 0); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Default destructor +TrappingCorrectionAm241::~TrappingCorrectionAm241() +{ + // Intentionally left blank +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Parse the command line +bool TrappingCorrectionAm241::ParseCommandLine(int argc, char** argv) +{ + ostringstream Usage; + Usage<"< i+1) && (argv[i+1][0] != '-' || isalpha(argv[i+1][1]) == 0))){ + cout<<"Error: Option "<>> Endpoints; + map>> FullDetEndpoints; + + // Store the HV and LV input files + vector FileNames; + FileNames.push_back(m_HVFileName); + FileNames.push_back(m_LVFileName); + + // Map the side integer to HV and LV labels + // i.e. 0=HV, 1=LV + vector IllumSide; + IllumSide.push_back(MString("HV")); + IllumSide.push_back(MString("LV")); + + // Make a directory in which to store the pixel-level data + MString PixelDir = m_OutFile + MString("_pixeldata"); + if (m_PixelCorrect==true) { + if (std::filesystem::create_directories(PixelDir.Data())==false) { + cout<<"Directory '"< HDFNames; + + TH1::SetDefaultSumw2(); + + map CTDHistograms; + map HVEnergyHistograms; + map LVEnergyHistograms; + + // Read in the input files and make a list of hdf5 files to calibrate + if ((InputFile.GetSubString(InputFile.Length() - 4)) == "hdf5") { + HDFNames.push_back(InputFile); + } else if ((InputFile.GetSubString(InputFile.Length() - 3)) == "txt") { + cout<<"Reading input file "<SetFileNameStripMap(m_StripMapFile); + Loader->SetFileName(File); + Loader->SetLoadContinuationFiles(m_ContinueHDF5); + S->SetModule(Loader, MNumber); + ++MNumber; + + cout<<"Creating TAC calibrator"<SetTACCalFileName(m_TACCalFile); + TACCalibrator->SetTACCutFileName(m_TACCutFile); + S->SetModule(TACCalibrator, MNumber); + ++MNumber; + + cout<<"Creating energy calibrator"<SetFileName(m_EcalFile); + EnergyCalibrator->EnablePreampTempCorrection(false); + S->SetModule(EnergyCalibrator, MNumber); + ++MNumber; + + cout<<"Creating Event filter"<SetMinimumLVStrips(1); + EventFilter->SetMaximumLVStrips(3); + EventFilter->SetMinimumHVStrips(1); + EventFilter->SetMaximumHVStrips(3); + EventFilter->SetMinimumHits(0); + EventFilter->SetMaximumHits(100); + EventFilter->SetMinimumTotalEnergy(m_MinEnergy); + EventFilter->SetMaximumTotalEnergy(m_MaxEnergy*2); // Multiply by 2 because this is the event-level energy, i.e. sum over both sides + S->SetModule(EventFilter, MNumber); + ++MNumber; + + cout<<"Creating strip pairing"<SetModule(Pairing, MNumber); + + cout<<"Initializing Loader"<Initialize() == false) return false; + cout<<"Initializing TAC calibrator"<Initialize() == false) return false; + cout<<"Initializing Energy calibrator"<Initialize() == false) return false; + cout<<"Initializing Event filter"<Initialize() == false) return false; + cout<<"Initializing Pairing"<Initialize() == false) return false; + + bool IsFinished = false; + MReadOutAssembly* Event = new MReadOutAssembly(); + + // Pass Events through each module. Once calibrated, add the Event to the histograms + cout<<"Analyzing..."<Clear(); + + if (Loader->IsReady()) { + + // Load Event from HDF5 file, then do TAC cut/calibration, energy calibration, and filtering + Loader->AnalyzeEvent(Event); + TACCalibrator->AnalyzeEvent(Event); + EnergyCalibrator->AnalyzeEvent(Event); + bool Unfiltered = EventFilter->AnalyzeEvent(Event); + + if (Unfiltered == true) { + + // Pair strips + Pairing->AnalyzeEvent(Event); + + // Only look at events with 1 Hit since we're interested in the 60keV photopeak. + if ((Event->HasAnalysisProgress(MAssembly::c_StripPairing) == true) && (Unfiltered==true) && (Event->GetNHits()==1)) { + + // for each Hit + for (unsigned int h = 0; h < Event->GetNHits(); ++h) { + double HVEnergy = 0.0; + double LVEnergy = 0.0; + double HVEnergyResolution = 0.0; + double LVEnergyResolution = 0.0; + vector HVStrips; + vector LVStrips; + + MHit* H = Event->GetHit(h); + + int DetID = H->GetStripHit(0)->GetDetectorID(); + + // for each Strip Hit in the Hit + for (unsigned int sh = 0; sh < H->GetNStripHits(); ++sh) { + MStripHit* SH = H->GetStripHit(sh); + + // Sum up the HV and LV energies, excluding nearest neighbors if m_ExcludeNN==true + if ((m_ExcludeNN==false) || ((m_ExcludeNN==true) && (SH->IsNearestNeighbor()==false))) { + if (SH->IsLowVoltageStrip()==true) { + LVEnergy += SH->GetEnergy(); + LVEnergyResolution += (SH->GetEnergyResolution())*(SH->GetEnergyResolution()); + LVStrips.push_back(SH); + } else { + HVEnergy += SH->GetEnergy(); + HVEnergyResolution += (SH->GetEnergyResolution())*(SH->GetEnergyResolution()); + HVStrips.push_back(SH); + } + } + } + + // As long as there are Strip Hits remaining on each side, add the energies and CTDs to histograms + if ((HVStrips.size()>0) && (LVStrips.size()>0)) { + + double HVEnergyFraction = 0; + double LVEnergyFraction = 0; + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + + if ((LVSH->HasCalibratedTiming()==true) && (HVSH->HasCalibratedTiming()==true)) { + + double CTD = LVSH->GetTiming() - HVSH->GetTiming(); + + int PixelID = (10000*DetID) + (100*LVSH->GetStripID()) + (HVSH->GetStripID()); + + // If this PixelID hasn't been encountered yet, make an entry for it in Endpoints + if (Endpoints.find(PixelID)==Endpoints.end()) { + vector tempHVvec; + vector tempLVvec; + vector> tempvec; + Endpoints[PixelID] = tempvec; + Endpoints[PixelID].push_back(tempHVvec); + Endpoints[PixelID].push_back(tempLVvec); + } + + TH1D* CTDHist = CTDHistograms[PixelID]; + TH1D* HVHist = HVEnergyHistograms[PixelID]; + TH1D* LVHist = LVEnergyHistograms[PixelID]; + + // If we didnt find an entry for the CTD and energy histograms, make new ones + if (CTDHist == nullptr) { + char name[64]; sprintf(name,"CTD: PixelID %d %s Illumination",PixelID,IllumSide[s].Data()); + CTDHist = new TH1D(name, name, (g_MaxCTD - g_MinCTD)/2, g_MinCTD, g_MaxCTD); + CTDHist->SetXTitle("CTD (ns)"); + CTDHist->SetYTitle("Hits"); + CTDHistograms[PixelID] = CTDHist; + } + if (HVHist == nullptr) { + char name[64]; sprintf(name,"HV energy: PixelID %d %s Illumination",PixelID,IllumSide[s].Data()); + HVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy)*2, m_MinEnergy, m_MaxEnergy); + HVHist->SetXTitle("HV Energy (keV)"); + HVHist->SetYTitle("Hits"); + HVEnergyHistograms[PixelID] = HVHist; + } + if (LVHist == nullptr) { + char name[64]; sprintf(name,"LV energy: PixelID %d %s Illumination",PixelID,IllumSide[s].Data()); + LVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy)*2, m_MinEnergy, m_MaxEnergy); + LVHist->SetXTitle("LV Energy (keV)"); + LVHist->SetYTitle("Hits"); + LVEnergyHistograms[PixelID] = LVHist; + } + + // Add the data to the histograms + CTDHist->Fill(CTD); + HVHist->Fill(HVEnergy); + LVHist->Fill(LVEnergy); + } + } + } + } + } + } + IsFinished = Loader->IsFinished(); + } + } + + map FullDetCTDHistograms; + map FullDetHVEnergyHistograms; + map FullDetLVEnergyHistograms; + + // Loop over all pixels for which there are pixel-level histograms and sum to get full detector histograms + // Fit to the pixel-level histograms and record the results + for (auto H: CTDHistograms) { + + int PixelID = H.first; + int DetID = (PixelID-(PixelID%10000))/10000; + + // If this DetID hasn't been encountered yet, make an entry for it in Endpoints + if (FullDetEndpoints.find(DetID)==FullDetEndpoints.end()) { + vector tempHVvec; + vector tempLVvec; + vector> tempvec; + FullDetEndpoints[DetID] = tempvec; + FullDetEndpoints[DetID].push_back(tempHVvec); + FullDetEndpoints[DetID].push_back(tempLVvec); + } + + TH1D* CTDHist = FullDetCTDHistograms[DetID]; + TH1D* HVHist = FullDetHVEnergyHistograms[DetID]; + TH1D* LVHist = FullDetLVEnergyHistograms[DetID]; + + // If we didnt find an entry for the CTD and energy histograms, make new ones + if (CTDHist == nullptr) { + char name[64]; sprintf(name,"CTD: DetID %d %s Illumination",DetID,IllumSide[s].Data()); + CTDHist = new TH1D(name, name, (g_MaxCTD - g_MinCTD)/2, g_MinCTD, g_MaxCTD); + CTDHist->SetXTitle("CTD (ns)"); + CTDHist->SetYTitle("Hits"); + FullDetCTDHistograms[DetID] = CTDHist; + } + if (HVHist == nullptr) { + char name[64]; sprintf(name,"HV energy: DetID %d %s Illumination",DetID,IllumSide[s].Data()); + HVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy)*2, m_MinEnergy, m_MaxEnergy); + HVHist->SetXTitle("HV Energy (keV)"); + HVHist->SetYTitle("Hits"); + FullDetHVEnergyHistograms[DetID] = HVHist; + } + if (LVHist == nullptr) { + char name[64]; sprintf(name,"LV energy: DetID %d %s Illumination",DetID,IllumSide[s].Data()); + LVHist = new TH1D(name, name, (m_MaxEnergy - m_MinEnergy)*2, m_MinEnergy, m_MaxEnergy); + LVHist->SetXTitle("LV Energy (keV)"); + LVHist->SetYTitle("Hits"); + FullDetLVEnergyHistograms[DetID] = LVHist; + } + + CTDHist->Add(CTDHist, H.second); + HVHist->Add(HVHist, HVEnergyHistograms[PixelID]); + LVHist->Add(LVHist, LVEnergyHistograms[PixelID]); + + if (m_PixelCorrect==true) { + + // Fit histograms and save the results + // Only perform fits if the total counts in each is above the threshold given by g_MinCounts + if ((H.second->Integral() > g_MinCounts) && (HVEnergyHistograms[PixelID]->Integral() > g_MinCounts) && (LVEnergyHistograms[PixelID]->Integral() > g_MinCounts)) { + + // Initialize and run the CTD fit + double CTDGuess = H.second->GetBinCenter(H.second->GetMaximumBin()); + TF1* CTDFunction = GenerateCTDFunction(CTDFitMin, CTDFitMax, CTDGuess, FlipSwitch); + TFitResultPtr CTDFit = H.second->Fit(CTDFunction, "SQ", "", CTDFitMin, CTDFitMax); + + // Save the results of the CTD fit + ofstream CTDFitFile(PixelDir +MString("/") + PixelID+MString("_CTDFitResult_")+IllumSide[s]+ MString("Illum.txt")); + streambuf* coutbuf = cout.rdbuf(); + cout.rdbuf(CTDFitFile.rdbuf()); + if (CTDFit.Get()) { + CTDFit->Print(); + } + cout.rdbuf(coutbuf); + CTDFitFile.close(); + + // Initialize and run the HV photopeak fit + TF1* PhotopeakFunctionHV = GeneratePhotopeakFunction(); + TFitResultPtr HVFit = HVEnergyHistograms[PixelID]->Fit(PhotopeakFunctionHV, "SQ", "", 55, 70); + + // Save the results of the HV photopeak fit + ofstream HVFitFile(PixelDir +MString("/") + PixelID+MString("_HVEnergyFitResult_")+IllumSide[s]+ MString("Illum.txt")); + coutbuf = cout.rdbuf(); + cout.rdbuf(HVFitFile.rdbuf()); + if (HVFit.Get()) { + HVFit->Print(); + } + cout.rdbuf(coutbuf); + HVFitFile.close(); + + // Initialize and run the LV photopeak fit + TF1* PhotopeakFunctionLV = GeneratePhotopeakFunction(); + TFitResultPtr LVFit = LVEnergyHistograms[PixelID]->Fit(PhotopeakFunctionLV, "SQ", "", 55, 70); + + // Save the results of the LV photopeak fit + ofstream LVFitFile(PixelDir +MString("/") + PixelID+MString("_LVEnergyFitResult_")+IllumSide[s]+ MString("Illum.txt")); + coutbuf = cout.rdbuf(); + cout.rdbuf(LVFitFile.rdbuf()); + if (LVFit.Get()) { + LVFit->Print(); + } + cout.rdbuf(coutbuf); + LVFitFile.close(); + + // Save the resulting HV and LV photopeak centroids and CTD centroids to the Endpoints map. + // Only save these values if the fits ran successfully and if their reduced chi-square is less than 5 + if ((!(CTDFit->IsEmpty())) && (!(HVFit->IsEmpty())) && (!(LVFit->IsEmpty())) && ((CTDFit->Chi2()/CTDFit->Ndf()) < 5) && ((HVFit->Chi2()/HVFit->Ndf()) < 5) && ((LVFit->Chi2()/LVFit->Ndf()) < 5)) { + Endpoints[PixelID][s].push_back(CTDFit->Parameter(2)); + Endpoints[PixelID][s].push_back(HVFit->Parameter(1)); + Endpoints[PixelID][s].push_back(LVFit->Parameter(1)); + } else { + cout<<"Fits failed for Pixel "<Write(); + CTDHistFile.Close(); + + TFile HVHistFile(PixelDir+MString("/")+PixelID+MString("_HVEnergyHist_") + IllumSide[s]+ MString("Illum.root"),"recreate"); + HVEnergyHistograms[PixelID]->Write(); + HVHistFile.Close(); + + TFile LVHistFile(PixelDir+MString("/")+PixelID+MString("_LVEnergyHist_") + IllumSide[s]+ MString("Illum.root"),"recreate"); + LVEnergyHistograms[PixelID]->Write(); + LVHistFile.Close(); + } + } + + // Do the same function fitting and recording as above, but for the full detectors rather than pixel-by-pixel + for (auto H: FullDetCTDHistograms) { + + int DetID = H.first; + + if ((H.second->Integral() > g_MinCounts) && (FullDetHVEnergyHistograms[DetID]->Integral() > g_MinCounts) && (FullDetLVEnergyHistograms[DetID]->Integral() > g_MinCounts)) { + + double CTDGuess = H.second->GetBinCenter(H.second->GetMaximumBin()); + TF1* CTDFunction = GenerateCTDFunction(CTDFitMin, CTDFitMax, CTDGuess, FlipSwitch); + TFitResultPtr CTDFit = H.second->Fit(CTDFunction, "SQ", "", CTDFitMin, CTDFitMax); + + TF1* PhotopeakFunctionHV = GeneratePhotopeakFunction(); + TFitResultPtr HVFit = FullDetHVEnergyHistograms[DetID]->Fit(PhotopeakFunctionHV, "SQ", "", 55, 70); + + TF1* PhotopeakFunctionLV = GeneratePhotopeakFunction(); + TFitResultPtr LVFit = FullDetLVEnergyHistograms[DetID]->Fit(PhotopeakFunctionLV, "SQ", "", 55, 70); + + if ((!(CTDFit->IsEmpty())) && (!(HVFit->IsEmpty())) && (!(LVFit->IsEmpty()))) { + FullDetEndpoints[DetID][s].push_back(CTDFit->Parameter(2)); + FullDetEndpoints[DetID][s].push_back(HVFit->Parameter(1)); + FullDetEndpoints[DetID][s].push_back(LVFit->Parameter(1)); + + ofstream CTDFitFile(DetID+MString("_CTDFitResult_")+IllumSide[s]+ MString("Illum.txt")); + streambuf* coutbuf = cout.rdbuf(); + cout.rdbuf(CTDFitFile.rdbuf()); + if (CTDFit.Get()) { + CTDFit->Print(); + } + cout.rdbuf(coutbuf); + CTDFitFile.close(); + + ofstream HVFitFile(DetID+MString("_HVEnergyFitResult_")+IllumSide[s]+ MString("Illum.txt")); + coutbuf = cout.rdbuf(); + cout.rdbuf(HVFitFile.rdbuf()); + if (HVFit.Get()) { + HVFit->Print(); + } + cout.rdbuf(coutbuf); + HVFitFile.close(); + + ofstream LVFitFile(DetID+MString("_LVEnergyFitResult_")+IllumSide[s]+ MString("Illum.txt")); + coutbuf = cout.rdbuf(); + cout.rdbuf(LVFitFile.rdbuf()); + if (LVFit.Get()) { + LVFit->Print(); + } + cout.rdbuf(coutbuf); + LVFitFile.close(); + + TFile CTDFile(m_OutFile+MString("_Det")+DetID+MString("_CTDHist_") + IllumSide[s]+ MString("Illum.root"),"recreate"); + + TCanvas* CTDCanvas = new TCanvas(); + CTDCanvas->cd(); + H.second->Draw("hist"); + CTDFunction->Draw("same"); + + H.second->Write(); + CTDFile.Close(); + + TFile HVHistFile(m_OutFile+MString("_Det")+DetID+MString("_HVEnergyHist_") + IllumSide[s]+ MString("Illum.root"),"recreate"); + + TCanvas* HVHistCanvas = new TCanvas(); + HVHistCanvas->cd(); + FullDetHVEnergyHistograms[DetID]->Draw("hist"); + PhotopeakFunctionHV->Draw("same"); + + FullDetHVEnergyHistograms[DetID]->Write(); + HVHistFile.Close(); + + TFile LVHistFile(m_OutFile+MString("_Det")+DetID+MString("_LVEnergyHist_") + IllumSide[s]+ MString("Illum.root"),"recreate"); + + TCanvas* LVHistCanvas = new TCanvas(); + LVHistCanvas->cd(); + FullDetLVEnergyHistograms[DetID]->Draw("hist"); + PhotopeakFunctionLV->Draw("same"); + + FullDetLVEnergyHistograms[DetID]->Write(); + LVHistFile.Close(); + + } else { + cout<<"Fits failed for Det "< DeltaHVMap; + map DeltaLVMap; + + map DeltaHVHist; + map DeltaLVHist; + + // Write the results of good fits to the parameter file and produce summary plots. + // If an individual pixel did not produce a good fit, default to the detector-level parameters. + // Loop over each detector + for (auto E: FullDetEndpoints) { + + int DetID = E.first; + + double HVIllumCTD = 0; + double LVIllumCTD = 0; + double HVIllumHVCentroid = 0; + double LVIllumHVCentroid = 0; + double HVIllumLVCentroid = 0; + double LVIllumLVCentroid = 0; + + if ((E.second[0].size() > 0) && (E.second[1].size() > 0)) { + HVIllumCTD = E.second[0][0]; + LVIllumCTD = E.second[1][0]; + HVIllumHVCentroid = E.second[0][1]; + LVIllumHVCentroid = E.second[1][1]; + HVIllumLVCentroid = E.second[0][2]; + LVIllumLVCentroid = E.second[1][2]; + } + + // Loop over each pixel + for (int hv=0; hvSetXTitle("HV Strip"); + TempHVMap->SetYTitle("LV Strip"); + TempHVMap->SetZTitle("Delta HV Energy"); + + char LVname[64]; sprintf(LVname,"Delta LV Map: Det %d",DetID); + TempLVMap = new TH2D(LVname, LVname, g_HVStrips, -0.5, g_HVStrips-0.5, g_LVStrips, -0.5, g_LVStrips-0.5); + TempLVMap->SetXTitle("HV Strip"); + TempLVMap->SetYTitle("LV Strip"); + TempLVMap->SetZTitle("Delta LV Energy"); + + DeltaHVMap[DetID] = TempHVMap; + DeltaLVMap[DetID] = TempLVMap; + } + + if (TempHVHist == nullptr) { + + char HVname[64]; sprintf(HVname,"Delta HV Hist: Det %d",DetID); + TempHVHist = new TH1D(HVname, HVname, 50, -3.0, 3.0); + TempHVHist->SetXTitle("Delta HV Energy (keV)"); + TempHVHist->SetYTitle("Number of pixels"); + + char LVname[64]; sprintf(LVname,"Delta LV Hist: Det %d",DetID); + TempLVHist = new TH1D(LVname, LVname, 50, -3.0, 3.0); + TempLVHist->SetXTitle("Delta LV Energy (keV)"); + TempLVHist->SetYTitle("Number of pixels"); + + DeltaHVHist[DetID] = TempHVHist; + DeltaLVHist[DetID] = TempLVHist; + } + + DeltaHVMap[DetID]->SetBinContent(hv+1, lv+1, -100); + DeltaLVMap[DetID]->SetBinContent(hv+1, lv+1, -100); + + // Retrieve the stored HV and LV energies and CTDs + // Plot the energy shifts in the HV and LV Maps, and fill in the histograms + if (Endpoints.find(PixelID)!=Endpoints.end()) { + if ((Endpoints[PixelID][0].size() > 0) && (Endpoints[PixelID][1].size() > 0)) { + + HVIllumCTD = Endpoints[PixelID][0][0]; + LVIllumCTD = Endpoints[PixelID][1][0]; + HVIllumHVCentroid = Endpoints[PixelID][0][1]; + LVIllumHVCentroid = Endpoints[PixelID][1][1]; + HVIllumLVCentroid = Endpoints[PixelID][0][2]; + LVIllumLVCentroid = Endpoints[PixelID][1][2]; + + // Normalize the differences in Energy to account for calibration differences + double HVDiff = (HVIllumHVCentroid - LVIllumHVCentroid)*(g_AmPhotopeak/HVIllumHVCentroid); + double LVDiff = (HVIllumLVCentroid - LVIllumLVCentroid)*(g_AmPhotopeak/HVIllumLVCentroid); + + DeltaHVMap[DetID]->SetBinContent(hv+1, lv+1, HVDiff); + DeltaLVMap[DetID]->SetBinContent(hv+1, lv+1, LVDiff); + + DeltaHVHist[DetID]->Fill(HVDiff); + DeltaLVHist[DetID]->Fill(LVDiff); + } + } + } + + // Record the Energies and CTDs in the parameter file + // Note that we default to the detector value if the pixel value isn't available and to 0 if the detector value isn't available. + OutputCalFile<cd(); + H.second->SetMinimum(-5.); + H.second->Draw("colz"); + H.second->Write(); + HVFile.Close(); + + TFile LVFile(m_OutFile+MString("_Det")+DetID+MString("_DeltaLVMap.root"),"recreate"); + TCanvas* LVCanvas = new TCanvas(); + LVCanvas->cd(); + DeltaLVMap[DetID]->SetMinimum(-5.); + DeltaLVMap[DetID]->Draw("colz"); + DeltaLVMap[DetID]->Write(); + LVFile.Close(); + + TFile HVHistFile(m_OutFile+MString("_Det")+DetID+MString("_DeltaHVHist.root"),"recreate"); + TCanvas* HVHistCanvas = new TCanvas(); + HVHistCanvas->cd(); + DeltaHVHist[DetID]->Draw("hist"); + DeltaHVHist[DetID]->Write(); + HVHistFile.Close(); + + TFile LVHistFile(m_OutFile+MString("_Det")+DetID+MString("_DeltaLVHist.root"),"recreate"); + TCanvas* LVHistCanvas = new TCanvas(); + LVHistCanvas->cd(); + DeltaLVHist[DetID]->Draw("hist"); + DeltaLVHist[DetID]->Write(); + LVHistFile.Close(); + } + } + + watch.Stop(); + cout<<"total time (s): "<SetParName(0, "Gauss norm"); + PhotopeakFunction->SetParName(1, "Mu"); + PhotopeakFunction->SetParName(2, "Sigma"); + PhotopeakFunction->SetParName(3, "Shelf norm"); + + PhotopeakFunction->SetParameter("Gauss norm", 1000); + PhotopeakFunction->SetParameter("Mu", 60); + PhotopeakFunction->SetParameter("Sigma", 2); + PhotopeakFunction->SetParameter("Shelf norm", 0.05); + + PhotopeakFunction->SetParLimits(0, 10, 1e8); + PhotopeakFunction->SetParLimits(1, 55, 65); + PhotopeakFunction->SetParLimits(2, 1.0, 10); + PhotopeakFunction->SetParLimits(3, 0, 0.1); + + return PhotopeakFunction; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +TF1* TrappingCorrectionAm241::GenerateCTDFunction(double CTDFitMin, double CTDFitMax, double CTDGuess, double FlipSwitch) +{ + // Exponentially modified gaussian + TF1* CTDFunction = new TF1("CTDFunction", "[0]*([1]/2)*exp(([1]/2)*(([1]*[3]*[3]) - 2*[4]*(x-[2])))*erfc((([1]*[3]*[3]) - [4]*(x-[2]))/([3]*sqrt(2)))", CTDFitMin, CTDFitMax); + + CTDFunction->SetParName(0, "Norm"); + CTDFunction->SetParName(1, "Lambda"); + CTDFunction->SetParName(2, "Mu"); + CTDFunction->SetParName(3, "Sigma"); + CTDFunction->SetParName(4, "Flip"); + + CTDFunction->SetParameter("Norm", 1000); + CTDFunction->SetParameter("Lambda", 0.05); + CTDFunction->SetParameter("Sigma", 12); + CTDFunction->SetParameter("Mu", CTDGuess); + CTDFunction->FixParameter(4, FlipSwitch); + + CTDFunction->SetParLimits(0, 0, 1e8); + CTDFunction->SetParLimits(1, 0.01, 1); + CTDFunction->SetParLimits(2, CTDFitMin, CTDFitMax); + CTDFunction->SetParLimits(3, 6, 30); + + return CTDFunction; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +TrappingCorrectionAm241* g_Prg = 0; +int g_NInterruptCatches = 1; + +MStripHit* TrappingCorrectionAm241::GetDominantStrip(vector& Strips, double& EnergyFraction) +{ + double MaxEnergy = -numeric_limits::max(); // AZ: When both energies are zero (which shouldn't happen) we still pick one + double TotalEnergy = 0.0; + MStripHit* MaxStrip = nullptr; + + // Iterate through strip hits and get the strip with highest energy + for (const auto SH : Strips) { + double Energy = SH->GetEnergy(); + TotalEnergy += Energy; + if (Energy > MaxEnergy) { + MaxStrip = SH; + MaxEnergy = Energy; + } + } + if (TotalEnergy == 0) { + EnergyFraction = 0; + } else { + EnergyFraction = MaxEnergy/TotalEnergy; + } + return MaxStrip; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Called when an interrupt signal is flagged +//! All catched signals lead to a well defined exit of the program +void CatchSignal(int a) +{ + if (g_Prg != 0 && g_NInterruptCatches-- > 0) { + cout<<"Catched signal Ctrl-C (ID="<Interrupt(); + } else { + abort(); + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Main program +int main(int argc, char** argv) +{ + // Catch a user interupt for graceful shutdown + signal(SIGINT, CatchSignal); + + // Initialize global MEGALIB variables, especially mgui, etc. + MGlobal::Initialize("Standalone", "a standalone example program"); + + TApplication TrappingCorrectionApp("TrappingCorrectionApp", 0, 0); + + g_Prg = new TrappingCorrectionAm241(); + + if (g_Prg->ParseCommandLine(argc, argv) == false) { + cerr<<"Error during parsing of command line!"<Analyze() == false) { + cerr<<"Error during analysis!"< +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIERBList.h" + +// NuSTAR libs +#include "MGUIExpo.h" + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +class MGUIExpoTACcut : public MGUIExpo +{ + // public Session: + public: + //! Default constructor + MGUIExpoTACcut(MModule* Module); + //! Default destructor + virtual ~MGUIExpoTACcut(); + + //! The creation part which gets overwritten + virtual void Create(); + + //! Update the frame + virtual void Update(); + + //! Reset the data in the UI + virtual void Reset(); + + //! Export the data in the UI + virtual void Export(const MString& FileName); + + //! Set the arrangment of the TAC histogram + //! 0 1 2 3 + //! 4 5 6 7 + //! 8 9 10 11 + void SetTACHistogramArrangement(const vector DetIDs); + + //! Set the energy histogram parameters + void SetTACHistogramParameters(unsigned int DetID, unsigned int NBins, double TACMin, double TACMax); + + //! Set the energy histogram parameters + void SetTACHistogramName(unsigned int DetID, MString Name); + + //! Add data to the TAC histogram + //! 0 1 2 3 + //! 4 5 6 7 + //! 8 9 10 11 + void AddTAC(unsigned int DetID, double TAC); + + // protected methods: + protected: + + + // protected members: + protected: + + // private members: + private: + //! TAC canvas + unordered_map m_TACCanvases; + //! TAC vs detector ID histogram + unordered_map m_TACHistograms; + + //! Detectors in x direction + unsigned int m_NColumns; + //! Detectors in y direction + unsigned int m_NRows; + + // Map the detector ID to the x,y position of histograms. + vector> m_DetectorMap; + + //! The number of bins of the histogram + unordered_map m_NBins; + //! The minimum TAC + unordered_map m_Min; + //! The maximum TAC + unordered_map m_Max; + + +#ifdef ___CLING___ + public: + ClassDef(MGUIExpoTACcut, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsDepthCalibration2024.h b/include/MGUIOptionsDepthCalibration2024.h index a9853d86..15c72ed6 100644 --- a/include/MGUIOptionsDepthCalibration2024.h +++ b/include/MGUIOptionsDepthCalibration2024.h @@ -74,9 +74,6 @@ class MGUIOptionsDepthCalibration2024 : public MGUIOptions //! Select spline file to load, splines will convert CTD->Depth MGUIEFileSelector* m_SplinesFileSelector; - //! Select TAC Calibration file to load, converts readout timing to nanoseconds - MGUIEFileSelector* m_TACCalFileSelector; - //! Check button if working with the Card Cage at UCSD TGCheckButton* m_UCSDOverride; diff --git a/include/MGUIOptionsEventSaver.h b/include/MGUIOptionsEventSaver.h index 9d2d405c..47eff234 100644 --- a/include/MGUIOptionsEventSaver.h +++ b/include/MGUIOptionsEventSaver.h @@ -43,6 +43,7 @@ //////////////////////////////////////////////////////////////////////////////// +//! UI for the event saver class MGUIOptionsEventSaver : public MGUIOptions { // public Session: @@ -78,6 +79,9 @@ class MGUIOptionsEventSaver : public MGUIOptions //! Checkbutton to save or reject bad events TGCheckButton* m_SaveBadEvents; + //! Checkbutton to save veto events + TGCheckButton* m_SaveVetoEvents; + //! Checkbutton to add a time tag TGCheckButton* m_AddTimeTag; @@ -85,6 +89,23 @@ class MGUIOptionsEventSaver : public MGUIOptions TGCheckButton* m_SplitFile; //! Entry field for the time after which to split the file MGUIEEntry* m_SplitFileTime; + + //! Checkbutton to include or exclude ADCs in the roa file + TGCheckButton* m_RoaWithADCs; + //! Checkbutton to include or exclude TACs in the roa file + TGCheckButton* m_RoaWithTACs; + //! Checkbutton to include or exclude energies in the roa file + TGCheckButton* m_RoaWithEnergies; + //! Checkbutton to include or exclude timings in the roa file + TGCheckButton* m_RoaWithTimings; + //! Checkbutton to include or exclude temperatures in the roa file + TGCheckButton* m_RoaWithTemperatures; + //! Checkbutton to include or exclude flags in the roa file + TGCheckButton* m_RoaWithFlags; + //! Checkbutton to include or exclude origins in the roa file + TGCheckButton* m_RoaWithOrigins; + //! Checkbutton to include or exclude nearest neighbor hits in the roa file + TGCheckButton* m_RoaWithNearestNeighbors; #ifdef ___CLING___ public: diff --git a/include/MGUIOptionsLoaderMeasurements.h b/include/MGUIOptionsLoaderMeasurements.h index 2e9082b4..bbcbd5e8 100644 --- a/include/MGUIOptionsLoaderMeasurements.h +++ b/include/MGUIOptionsLoaderMeasurements.h @@ -46,7 +46,7 @@ class MGUIOptionsLoaderMeasurements : public MGUIOptions // public Session: public: //! Default constructor - MGUIOptionsLoaderMeasurements(MModule* Module); + MGUIOptionsLoaderMeasurements(MModule* Module, MString FileType); //! Default destructor virtual ~MGUIOptionsLoaderMeasurements(); @@ -71,6 +71,8 @@ class MGUIOptionsLoaderMeasurements : public MGUIOptions //! Select which file to load MGUIEFileSelector* m_FileSelector; + //! The file type to load + MString m_FileType; #ifdef ___CLING___ public: diff --git a/include/MGUIOptionsLoaderMeasurementsHDF.h b/include/MGUIOptionsLoaderMeasurementsHDF.h new file mode 100644 index 00000000..b453cb2a --- /dev/null +++ b/include/MGUIOptionsLoaderMeasurementsHDF.h @@ -0,0 +1,93 @@ +/* + * MGUIOptionsLoaderMeasurementsHDF.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIOptionsLoaderMeasurementsHDF__ +#define __MGUIOptionsLoaderMeasurementsHDF__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIEFileSelector.h" +#include "MGUIOptions.h" + +// Nuclearizer libs: +#include "MModule.h" + + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! UI settings for the HDF measurements loader +class MGUIOptionsLoaderMeasurementsHDF : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsLoaderMeasurementsHDF(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsLoaderMeasurementsHDF(); + + //! Process all button, etc. messages + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + + //! The creation part which gets overwritten + virtual void Create(); + + // protected methods: + protected: + + //! Actions after the Apply or OK button has been pressed + virtual bool OnApply(); + + + // protected members: + protected: + + // private members: + private: + //! Select which file to load + MGUIEFileSelector* m_FileSelectorHDF; + + //! Check button for switch between loading continuation files or not + TGCheckButton* m_LoadContinuationFiles; + + //! Select which file to load + MGUIEFileSelector* m_FileSelectorStripMap; + + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsLoaderMeasurementsHDF, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MGUIOptionsTACcut.h b/include/MGUIOptionsTACcut.h new file mode 100644 index 00000000..d83df4cd --- /dev/null +++ b/include/MGUIOptionsTACcut.h @@ -0,0 +1,94 @@ +/* + * MGUIOptionsTACcut.h + * + * Copyright (C) 2008-2010 by Jau-Shian Liang. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MGUIOptionsTACcut__ +#define __MGUIOptionsTACcut__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// ROOT libs: +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MGlobal.h" +#include "MGUIERBList.h" +#include "MModule.h" +#include "MGUIOptions.h" + +// Forward declarations: +class MGUIEFileSelector; +class MGUIEMinMaxEntry; +class MGUIEEntry; + +//////////////////////////////////////////////////////////////////////////////// + + +class MGUIOptionsTACcut : public MGUIOptions +{ + // public Session: + public: + //! Default constructor + MGUIOptionsTACcut(MModule* Module); + //! Default destructor + virtual ~MGUIOptionsTACcut(); + + //! Process all button, etc. messages + virtual bool ProcessMessage(long Message, long Parameter1, long Parameter2); + + //! The creation part which gets overwritten + virtual void Create(); + + // protected methods: + protected: + + //! Actions after the Apply or OK button has been pressed + virtual bool OnApply(); + + + // protected members: + protected: + //! The detector IDs as a string + TGTextEntry* m_Detectors; + + //! Select TAC Calibration file to load, converts readout timing to nanoseconds + MGUIEFileSelector* m_TACCalFileSelector; + + //! Select TAC Cut file to load, which specifies the parameters for removing strip hits + MGUIEFileSelector* m_TACCutFileSelector; + + + // private members: + private: + + +#ifdef ___CLING___ + public: + ClassDef(MGUIOptionsTACcut, 1) // basic class for dialog windows +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleDepthCalibration2024.h b/include/MModuleDepthCalibration2024.h index 5e03bb0a..21f48291 100644 --- a/include/MModuleDepthCalibration2024.h +++ b/include/MModuleDepthCalibration2024.h @@ -72,11 +72,6 @@ class MModuleDepthCalibration2024 : public MModule //! Get filename for CTD->Depth splines MString GetSplinesFileName() const {return m_SplinesFile;} - //! Set filename for TAC Calibration - void SetTACCalFileName( const MString& FileName) {m_TACCalFile = FileName;} - //! Get filename for TAC Calibration - MString GetTACCalFileName() const {return m_TACCalFile;} - //! Set whether the data came from the card cage at UCSD void SetUCSDOverride( bool Override ) {m_UCSDOverride = Override;} //! Get whether the data came from the card cage at UCSD @@ -109,13 +104,11 @@ class MModuleDepthCalibration2024 : public MModule //! Load in the specified coefficients file bool LoadCoeffsFile(MString FName); //! Return the coefficients for a pixel - vector* GetPixelCoeffs(int pixel_code); + vector* GetPixelCoeffs(int PixelCode); //! Load the splines file bool LoadSplinesFile(MString FName); - //! Load the TAC Calibration file - bool LoadTACCalFile(MString FName); //! Get the timing FWHM noise for the specified pixel and Energy - double GetTimingNoiseFWHM(int pixel_code, double Energy); + double GetTimingNoiseFWHM(int PixelCode, double Energy); // private methods @@ -127,12 +120,8 @@ class MModuleDepthCalibration2024 : public MModule unordered_map> m_Coeffs; double m_Coeffs_Energy; - unordered_map>> m_HVTACCal; - unordered_map>> m_LVTACCal; MString m_CoeffsFile; MString m_SplinesFile; - MString m_TACCalFile; - unordered_map m_DetectorNames; unordered_map m_Thicknesses; unordered_map m_NXStrips; unordered_map m_NYStrips; @@ -146,7 +135,9 @@ class MModuleDepthCalibration2024 : public MModule uint64_t m_Error5; uint64_t m_Error6; uint64_t m_ErrorSH; - vector m_Detectors; + uint64_t m_ErrorNullSH; + uint64_t m_ErrorNoE; + unordered_map m_Detectors; vector m_DetectorIDs; MModuleEnergyCalibrationUniversal* m_EnergyCalibration; MGUIExpoDepthCalibration2024* m_ExpoDepthCalibration; @@ -156,7 +147,6 @@ class MModuleDepthCalibration2024 : public MModule unordered_map> m_DepthGrid; bool m_SplinesFileIsLoaded; bool m_CoeffsFileIsLoaded; - bool m_TACCalFileIsLoaded; // boolean for use with the card cage at UCSD since it tags all events as detector 11 bool m_UCSDOverride; diff --git a/include/MModuleEventSaver.h b/include/MModuleEventSaver.h index 28f1e5ad..7211f27d 100644 --- a/include/MModuleEventSaver.h +++ b/include/MModuleEventSaver.h @@ -58,10 +58,15 @@ class MModuleEventSaver : public MModule //! Set the file name void SetFileName(const MString& Name) { m_FileName = Name; } - //! Get the file name + //! Return true if the Bad events should be saved bool GetSaveBadEvents() const { return m_SaveBadEvents; } - //! Set the file name + //! Set whether the Bad events should be saved void SetSaveBadEvents(bool SaveBadEvents) { m_SaveBadEvents = SaveBadEvents; } + + //! Return true if the Veto events should be saved + bool GetSaveVetoEvents() const { return m_SaveVetoEvents; } + //! Set whether the Veto events should be saved + void SetSaveVetoEvents(bool SaveVetoEvents) {m_SaveVetoEvents = SaveVetoEvents; } //! Get the add time tag flag bool GetAddTimeTag() const { return m_AddTimeTag; } @@ -77,7 +82,47 @@ class MModuleEventSaver : public MModule MTime GetSplitFileTime() const { return m_SplitFileTime; } //! Set the time after which the file should be split void SetSplitFileTime(MTime SplitFileTime) { m_SplitFileTime = SplitFileTime; } - + + //! Return whether ADCs should be included in the roa file + bool GetRoaWithADCs() const { return m_RoaWithADCs; } + //! Set whether ADCs should be included in the roa file + void SetRoaWithADCs(bool Flag) { m_RoaWithADCs = Flag; } + + //! Return whether TACs should be included in the roa file + bool GetRoaWithTACs() const { return m_RoaWithTACs; } + //! Set whether TACs should be included in the roa file + void SetRoaWithTACs(bool Flag) { m_RoaWithTACs = Flag; } + + //! Return whether energies should be included in the roa file + bool GetRoaWithEnergies() const { return m_RoaWithEnergies; } + //! Set whether energies should be included in the roa file + void SetRoaWithEnergies(bool Flag) { m_RoaWithEnergies = Flag; } + + //! Return whether timings should be included in the roa file + bool GetRoaWithTimings() const { return m_RoaWithTimings; } + //! Set whether timings should be included in the roa file + void SetRoaWithTimings(bool Flag) { m_RoaWithTimings = Flag; } + + //! Return whether temperatures should be included in the roa file + bool GetRoaWithTemperatures() const { return m_RoaWithTemperatures; } + //! Set whether temperatures should be included in the roa file + void SetRoaWithTemperatures(bool Flag) { m_RoaWithTemperatures = Flag; } + + //! Return whether flags should be included in the roa file + bool GetRoaWithFlags() const { return m_RoaWithFlags; } + //! Set whether flags should be included in the roa file + void SetRoaWithFlags(bool Flag) { m_RoaWithFlags = Flag; } + + //! Return whether origins should be included in the roa file + bool GetRoaWithOrigins() const { return m_RoaWithOrigins; } + //! Set whether origins should be included in the roa file + void SetRoaWithOrigins(bool Flag) { m_RoaWithOrigins = Flag; } + + //! Return whether nearest neighbors should be included in the roa file + bool GetRoaWithNearestNeighbors() const { return m_RoaWithNearestNeighbors; } + //! Set whether nearest neighbors should be included in the roa file + void SetRoaWithNearestNeighbors(bool Flag) { m_RoaWithNearestNeighbors = Flag; } + //! Set the start area of the far field simulation if there was any void SetStartAreaFarField(double Area) { m_StartAreaFarField = Area; } //! Set the number if simulated events @@ -148,6 +193,9 @@ class MModuleEventSaver : public MModule //! Save bad events bool m_SaveBadEvents; + //! Save Veto events + bool m_SaveVetoEvents; + //! Add a time tag to the file bool m_AddTimeTag; @@ -158,6 +206,25 @@ class MModuleEventSaver : public MModule bool m_SplitFile; //! If we split the file, this is the time in seconds after which we split MTime m_SplitFileTime; + + // Roa options + + //! If true write ADC values to the roa file + bool m_RoaWithADCs; + //! If true write TAC values to the roa file + bool m_RoaWithTACs; + //! If true write energy values to the roa file + bool m_RoaWithEnergies; + //! If true write timing values to the roa file + bool m_RoaWithTimings; + //! If true write temperature values to the roa file + bool m_RoaWithTemperatures; + //! If true write flags to the roa file + bool m_RoaWithFlags; + //! If true write origins to the roa file + bool m_RoaWithOrigins; + //! True if we should include next neighbors in the data stream + bool m_RoaWithNearestNeighbors; //! Main output stream for file MFile m_Out; diff --git a/include/MModuleLoaderMeasurements.h b/include/MModuleLoaderMeasurements.h index a1c31e39..bf9b59e0 100644 --- a/include/MModuleLoaderMeasurements.h +++ b/include/MModuleLoaderMeasurements.h @@ -51,9 +51,6 @@ class MModuleLoaderMeasurements : public MModule, public MFileEvents //! Main data analysis routine, which updates the event to a new level virtual bool AnalyzeEvent(MReadOutAssembly* Event); - //! Show the options GUI - virtual void ShowOptionsGUI(); - //! Read the configuration data from an XML node virtual bool ReadXmlConfiguration(MXmlNode* Node) = 0; //! Create an XML node tree from the configuration diff --git a/include/MModuleLoaderMeasurementsHDF.h b/include/MModuleLoaderMeasurementsHDF.h new file mode 100644 index 00000000..a926e449 --- /dev/null +++ b/include/MModuleLoaderMeasurementsHDF.h @@ -0,0 +1,239 @@ +/* + * MModuleLoaderMeasurementsHDF.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleLoaderMeasurementsHDF__ +#define __MModuleLoaderMeasurementsHDF__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" +#include "MFileReadOuts.h" + +// Nuclearizer libs: +#include "MStripMap.h" +#include "MModuleLoaderMeasurements.h" + +// H5 libs +#include "H5Cpp.h" +using namespace H5; + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! The versions of the strip hits stored in the HDF file +enum class MHDFStripHitVersion { V1_0, V1_2 }; + +//! And a streamer for it +inline ostream& operator<<(ostream& os, MHDFStripHitVersion Version) { + switch (Version) { + case MHDFStripHitVersion::V1_0: return os<<"1.0"; + case MHDFStripHitVersion::V1_2: return os<<"1.2"; + default: return os<<"Unknown"; + } +} + +//! Version 1.0 & 1.1 of the HDF5 hit info +struct MHDFStripHit_V1_0 { + uint16_t m_EventID; + uint32_t m_TimeCode; + uint8_t m_HitType; + uint8_t m_TimingType; + uint16_t m_StripID; + uint8_t m_CrystalID; + uint8_t m_Gain; + uint8_t m_Overflow; + uint16_t m_CurrentMaximum; + uint16_t m_HighCurrentSamples; + uint16_t m_EnergyData; + uint16_t m_EnergyDataLowGain; + uint16_t m_EnergyDataHighGain; + uint16_t m_TimingData; + uint8_t m_Pad; + uint8_t m_Hits; + uint8_t m_EventType; + uint8_t m_CRC; +}; + +//! Version 1.2 of the HDF5 hit info +struct MHDFStripHit_V1_2 { + uint16_t m_EventID; + uint64_t m_TimeCode; + double m_GSETimeCode; + uint8_t m_HitType; + uint8_t m_TimingType; + uint16_t m_StripID; + uint8_t m_CrystalID; + uint8_t m_Gain; + uint8_t m_Overflow; + uint16_t m_CurrentMaximum; + uint16_t m_HighCurrentSamples; + uint16_t m_EnergyData; + uint16_t m_EnergyDataLowGain; + uint16_t m_EnergyDataHighGain; + uint16_t m_TimingData; + uint8_t m_Pad; + uint8_t m_Hits; + uint16_t m_Bytes; + uint8_t m_EventType; + uint8_t m_CRC; +}; + +//! The version string +struct MHDFStripHitVersionString { + char string_col[256]; +}; + +//////////////////////////////////////////////////////////////////////////////// + + +//! A module to load HDF5 data files +class MModuleLoaderMeasurementsHDF : public MModuleLoaderMeasurements +{ + // public interface: + public: + //! Default constructor + MModuleLoaderMeasurementsHDF(); + //! Default destructor + virtual ~MModuleLoaderMeasurementsHDF(); + + //! Create a new object of this class + virtual MModuleLoaderMeasurementsHDF* Clone() { return new MModuleLoaderMeasurementsHDF(); } + + //! Get the file name of the strip map + MString GetFileNameStripMap() const { return m_FileNameStripMap; } + //! Set the file name of the strip map + void SetFileNameStripMap(const MString& Name) { m_FileNameStripMap = Name; } + + + //! Enable/Disable loading continuation files + bool GetLoadContinuationFiles() const { return m_LoadContinuationFiles; } + //! Set loading continuation files + void SetLoadContinuationFiles(bool LoadContinuationFiles) { m_LoadContinuationFiles = LoadContinuationFiles; } + + //! Initialize the module + virtual bool Initialize(); + + //! Initialize the module + virtual void Finalize(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + //! Read the configuration data from an XML node + virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration + virtual MXmlNode* CreateXmlConfiguration(); + + + // protected methods: + protected: + //! Convert more data from raw to intermediate format - return false if no more data can be converted + bool OpenHDF5File(MString FileName); + //! Read a batch of hits using a hyperslab + bool ReadBatchHits(); + + // private methods: + private: + + + + // protected members: + protected: + + + // private members: + private: + /* + //! Start of the observation time + MTime m_StartObservationTime; + //! Clock time belonging to the start of the observation time + unsigned long m_StartClock; + //! End of the observation time + MTime m_EndObservationTime; + //! Clock time belonging to the end of the observation time + unsigned long m_EndClock; + */ + + //! The HDF5 file + H5File m_HDFFile; + + //! True, if we want to load continuation files + bool m_LoadContinuationFiles; + + //! The ID of the currently loaded continuation file + unsigned int m_ContinuationFileID; + + //! The HDF5 data set + DataSet m_HDFDataSet; + + //! The HDF5 compond data type + CompType m_HDFCompoundDataType; + + //! The version of the strip hit structure + MHDFStripHitVersion m_HDFStripHitVersion; + + //! The default batch size + static constexpr unsigned int m_DefaultBatchSize = 10000; + + //! The current batch size + unsigned int m_CurrentBatchSize; + + //! The current index in the batch + unsigned int m_CurrentBatchIndex; + + // The various batches: + //! The MHDFStripHit_V1_0 batch: + vector m_Buffer_1_0; + //! The MHDFStripHit_V1_2 batch: + vector m_Buffer_1_2; + + //! Total number of hits in a file + hsize_t m_TotalHits; + + //! Current hit number in the file + hsize_t m_CurrentHit; + + //! Number of event ID roll-overs: + unsigned int m_NumberOfEventIDRollOvers; + + //! The last handled event ID + unsigned int m_LastEventID; + + //! The file name of the strip map + MString m_FileNameStripMap; + + //! The strip map + MStripMap m_StripMap; + +#ifdef ___CLING___ + public: + ClassDef(MModuleLoaderMeasurementsHDF, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleLoaderMeasurementsROA.h b/include/MModuleLoaderMeasurementsROA.h index 9db50119..3b302782 100644 --- a/include/MModuleLoaderMeasurementsROA.h +++ b/include/MModuleLoaderMeasurementsROA.h @@ -57,6 +57,9 @@ class MModuleLoaderMeasurementsROA : public MModuleLoaderMeasurements //! Main data analysis routine, which updates the event to a new level virtual bool AnalyzeEvent(MReadOutAssembly* Event); + //! Show the options GUI + virtual void ShowOptionsGUI(); + //! Read the configuration data from an XML node virtual bool ReadXmlConfiguration(MXmlNode* Node); //! Create an XML node tree from the configuration diff --git a/include/MModuleNearestNeighbor.h b/include/MModuleNearestNeighbor.h new file mode 100644 index 00000000..c46f05c5 --- /dev/null +++ b/include/MModuleNearestNeighbor.h @@ -0,0 +1,90 @@ +/* + * MModuleNearestNeighbor.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleNearestNeighbor__ +#define __MModuleNearestNeighbor__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" +#include "MModule.h" + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +class MModuleNearestNeighbor : public MModule +{ + // public interface: + public: + //! Default constructor + MModuleNearestNeighbor(); + //! Default destructor + virtual ~MModuleNearestNeighbor(); + + //! Create a new object of this class + virtual MModuleNearestNeighbor* Clone() { return new MModuleNearestNeighbor(); } + + //! Initialize the module + virtual bool Initialize(); + + //! Finalize the module + virtual void Finalize(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + //! Read the configuration data from an XML node + virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration + virtual MXmlNode* CreateXmlConfiguration(); + + // protected methods: + protected: + + // private methods: + private: + + + + // protected members: + protected: + + + // private members: + private: + + + + +#ifdef ___CLING___ + public: + ClassDef(MModuleNearestNeighbor, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MModuleTACcut.h b/include/MModuleTACcut.h new file mode 100644 index 00000000..7bf3836a --- /dev/null +++ b/include/MModuleTACcut.h @@ -0,0 +1,134 @@ +/* + * MModuleTACcut.h + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MModuleTACcut__ +#define __MModuleTACcut__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include + +// ROOT libs: +#include "TGClient.h" +#include "TH1.h" + +// MEGAlib libs: +#include "MGlobal.h" +#include "MModule.h" +#include "MGUIExpoTACcut.h" + + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +class MModuleTACcut : public MModule +{ + // public interface: + public: + //! Default constructor + MModuleTACcut(); + //! Default destructor + virtual ~MModuleTACcut(); + + //! Create a new object of this class + virtual MModuleTACcut* Clone() { return new MModuleTACcut(); } + + //! Initialize the module + virtual bool Initialize(); + + //! Create expos + virtual void CreateExpos(); + + //! Finalize the module + virtual void Finalize(); + + //! Main data analysis routine, which updates the event to a new level + virtual bool AnalyzeEvent(MReadOutAssembly* Event); + + //! Show the options GUI + virtual void ShowOptionsGUI(); + + + //! Read the configuration data from an XML node + virtual bool ReadXmlConfiguration(MXmlNode* Node); + //! Create an XML node tree from the configuration + virtual MXmlNode* CreateXmlConfiguration(); + + ///////////// Creating functions that will update and get the min/max TAC values ////////////////////////// + + //! Set filename for TAC Calibration + void SetTACCalFileName( const MString& FileName) {m_TACCalFile = FileName;} + //! Get filename for TAC Calibration + MString GetTACCalFileName() const {return m_TACCalFile;} + + //! Set filename for TAC Cut + void SetTACCutFileName( const MString& FileName) {m_TACCutFile = FileName;} + //! Get filename for TAC Cut + MString GetTACCutFileName() const {return m_TACCutFile;} + + //! Load the TAC calibration file + bool LoadTACCalFile(MString FName); + + //! Load the TAC cut file + bool LoadTACCutFile(MString FName); + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + + // protected methods: + protected: + + + // private methods: + private: + + + + // protected members: + protected: + + // private members: + private: + +//! TAC cut and TAC calibration parameter files +MString m_TACCalFile; +MString m_TACCutFile; + +//! Map DetID -> Side (LV=0, HV=1) -> Strip ID -> TAC calibration/cut parameters +unordered_map>>> m_TACCal; +unordered_map>>> m_TACCut; + +//! Map characters representing sides of the detectors indices to avoid mistakes +unordered_map m_SideToIndex; + +vector m_DetectorIDs; + +MGUIExpoTACcut* m_ExpoTACcut; + + +#ifdef ___CLING___ + public: + ClassDef(MModuleTACcut, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/include/MReadOutAssembly.h b/include/MReadOutAssembly.h index 8f5e0fa3..a332f860 100644 --- a/include/MReadOutAssembly.h +++ b/include/MReadOutAssembly.h @@ -109,25 +109,15 @@ class MReadOutAssembly : public MReadOutSequence //! Find out if the event contains strip hits in a given detector bool InDetector(int DetectorID); - //! Set the vetoed flag - void SetVeto(bool Veto = true) { m_Veto = Veto; } - //! Return the veto flag - bool GetVeto() const { return m_Veto; } - - //! Set the guard ring 0 veto flag - void SetGR0Veto(bool Veto = true) { m_VetoGR0 = Veto; } - //! Get the guard ring 0 veto flag - bool GetGR0Veto() const { return m_VetoGR0; } - - //! Set the guard ring 1 veto flag - void SetGR1Veto(bool Veto = true) { m_VetoGR1 = Veto; } - //! Get the guard ring 1 veto flag - bool GetGR1Veto() const { return m_VetoGR1; } + //! Set the guard ring veto flag + void SetGuardRingVeto(bool Veto = true) { m_GuardRingVeto = Veto; } + //! Get the guard ring veto flag + bool GetGuardRingVeto() const { return m_GuardRingVeto; } //! Set the shield veto flag - void SetShieldVeto(bool Veto = true) { m_VetoShield = Veto; } + void SetShieldVeto(bool Veto = true) { m_ShieldVeto = Veto; } //! Get the shield veto flag - bool GetShieldVeto() const { return m_VetoShield; } + bool GetShieldVeto() const { return m_ShieldVeto; } //! Set the triggered flag void SetTrigger(bool Trigger = true) { m_Trigger = Trigger; } @@ -254,6 +244,9 @@ class MReadOutAssembly : public MReadOutSequence //! Get the filgtered-out flag bool IsFilteredOut() const { return m_FilteredOut; } + //! Returns true if any of the "veto" flags have been set + bool IsVeto() const; + //! Returns true if none of the "bad" or "incomplete" flags has been set and the event has not been filtered out or rejected bool IsGood() const; //! Returns true if any of the "bad" or "incomplete" flags has been set @@ -276,7 +269,7 @@ class MReadOutAssembly : public MReadOutSequence //! Stream the content in MEGAlib's evta format void StreamEvta(ostream& S); //! Stream the content in MEGAlib's roa format - void StreamRoa(ostream& S, bool WithDescriptor = true); + void StreamRoa(ostream& S, bool WithADCs = true, bool WithTACs = true, bool WithEnergies = false, bool WithTimings = false, bool WithTemperatures = false, bool WithFlags = false, bool WithOrigins = false, bool WithNearestNeighbors = false); //! Build the next MReadoutAssemply from a .dat file bool GetNextFromDatFile(MFile &F); //! Use the info in m_Aspect to turn m_CL into an absolute UTC time @@ -335,17 +328,11 @@ class MReadOutAssembly : public MReadOutSequence //! Quality of this event double m_EventQuality; - //! Veto flag of this event - bool m_Veto; - - //! Guard ring 0 veto flag - bool m_VetoGR0; - - //! Guard ring 1 veto flag - bool m_VetoGR1; + //! Guard ring veto flag + bool m_GuardRingVeto; //! Shield veto flag - bool m_VetoShield; + bool m_ShieldVeto; //! Trigger flag of this event bool m_Trigger; diff --git a/include/MStripHit.h b/include/MStripHit.h index d9aa131a..bb72b792 100644 --- a/include/MStripHit.h +++ b/include/MStripHit.h @@ -101,11 +101,26 @@ class MStripHit //! Return the calibrated energy double GetEnergyResolution() const { return m_EnergyResolution; } - //! Set the Timing of the top side + //! Set the TAC + void SetTAC(double TAC) { m_TAC = TAC; } + //! Return the TAC + double GetTAC() const { return m_TAC; } + + //! Set the TAC resolution + void SetTACResolution(double TACResolution) { m_TACResolution = TACResolution; } + //! Return the TAC resolution + double GetTACResolution() const { return m_TACResolution; } + + //! Set the Timing in nanoseconds void SetTiming(double Timing) { m_Timing = Timing; } - //! Return the Timing of the top side + //! Return the Timing in nanoseconds double GetTiming() const { return m_Timing; } + //! Set the Timing resolution + void SetTimingResolution(double TimingResolution) { m_TimingResolution = TimingResolution; } + //! Return the Timing resolution + double GetTimingResolution() const { return m_TimingResolution; } + //! Set the Temperature of the relavent preamp (in degrees C) void SetPreampTemp(double PreampTemp) { m_PreampTemp = PreampTemp; } //! Return the Temperature of the relavent preamp (in degrees C) @@ -115,15 +130,37 @@ class MStripHit void AddOrigins(vector Origins); //! Get the origins from the simulation vector GetOrigins() const { return m_Origins; } - - - + + //! Set the Guard Ring flag + void IsGuardRing(bool GuardRing) { m_IsGuardRing = GuardRing; } + //! Return a boolean indicating whether the strip is a Guard Ring + bool IsGuardRing() const { return m_IsGuardRing; } + //! Set the Nearest Neighbor flag + void IsNearestNeighbor(bool NearestNeighbor) { m_IsNearestNeighbor = NearestNeighbor; } + //! Return a boolean indicating whether the strip is a Nearest Neighbor + bool IsNearestNeighbor() const { return m_IsNearestNeighbor; } + + //! Set the Fast Timing flag + void HasFastTiming(bool FastTiming) { m_HasFastTiming = FastTiming; } + //! Return a boolean indicating whether the strip timing is fast; + bool HasFastTiming() const { return m_HasFastTiming; } + + //! Set the Calibrated Timing flag + void HasCalibratedTiming(bool CalibratedTiming) { m_HasCalibratedTiming = CalibratedTiming; } + //! Return a boolean indicating whether the strip timing has been calibrated; + bool HasCalibratedTiming() const { return m_HasCalibratedTiming; } + + //! Produce an unsigned int with bitwise values representing flags + unsigned int MakeFlags(); + //! Read in unsigned int with bitwise values representing flags and update boolean flags + void ParseFlags(unsigned int Flags); + //! Parse some content from a line bool Parse(MString& Line, int Version = 1); //! Dump the content into a file stream bool StreamDat(ostream& S, int Version = 1); //! Stream the content in MEGAlib's roa format - void StreamRoa(ostream& S); + void StreamRoa(ostream& S, bool WithADC = true, bool WithTAC = true, bool WithEnergy = false, bool WithTiming = false, bool WithTemperature = false, bool WithFlags = false, bool WithOrigins = false); // protected methods: @@ -152,11 +189,26 @@ class MStripHit double m_Energy; //! The energy resolution double m_EnergyResolution; - //! Timing of the top side + //! TAC timing + double m_TAC; + //! TAC timing resolution + double m_TACResolution; + //! Timing in ns double m_Timing; + //! Timing resolution in ns + double m_TimingResolution; //! Temperature of Preamp double m_PreampTemp; - + + //! Flags denoting the type of strip hit + bool m_IsGuardRing; + bool m_IsNearestNeighbor; + + //! Flag indicating whether the hit has fast timing + bool m_HasFastTiming; + //! Flag indicating whether the hit has calibrated timing + bool m_HasCalibratedTiming; + //! Origin IAs from simulations vector m_Origins; diff --git a/include/MStripMap.h b/include/MStripMap.h new file mode 100644 index 00000000..cce86254 --- /dev/null +++ b/include/MStripMap.h @@ -0,0 +1,105 @@ +/* + * MStripMap.h + * + * Copyright (C) by Andreas Zoglauer + * All rights reserved. + * + * Please see the source-file for the copyright-notice. + * + */ + + +#ifndef __MStripMap__ +#define __MStripMap__ + + +//////////////////////////////////////////////////////////////////////////////// + + +// Standard libs: +#include +#include +using namespace std; + +// ROOT libs: + +// MEGAlib libs: +#include "MGlobal.h" + +// Forward declarations: + + +//////////////////////////////////////////////////////////////////////////////// + + +//! This class represents the mapping from asic channels to detector, side, and strip ID +class MStripMap +{ + // public interface: + public: + //! Default constructor + MStripMap(); + //! Default destructor + virtual ~MStripMap(); + + //! Load a strip map - return false on error + bool Open(MString FileName); + + //! Check if we have a certain read-out ID + bool HasReadOutID(unsigned int ROI) const; + + //! Get detector by read out ID - check with HasReadOutID(ROI) first + unsigned int GetDetectorID(unsigned int ROI) const; + + //! Get detector side by read out ID - check with HasReadOutID(ROI) first + bool IsLowVoltage(unsigned int ROI) const; + + //! Get strip ID by read out ID - check with HasReadOutID(ROI) first + unsigned int GetStripNumber(unsigned int ROI) const; + + + // protected methods: + protected: + //! Return the index of the read-out ID or throw an exception + unsigned int GetReadOutIDIndex(unsigned int ROI) const; + + + // private methods: + private: + + + + // protected members: + protected: + + + // private members: + private: + //! The internal struct for the map + struct MSingleStripMapping { + unsigned int m_ReadOutID; + unsigned int m_RTB; + unsigned int m_DRM; + bool m_IsPrimary; + unsigned int m_ASICID; + unsigned int m_ChannelID; + unsigned int m_DetectorID; + bool m_IsLowVoltage; + unsigned int m_StripNumber; + }; + + //! The strip mapping data + vector m_StripMappings; + + +#ifdef ___CLING___ + public: + ClassDef(MStripMap, 0) // no description +#endif + +}; + +#endif + + +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MAssembly.cxx b/src/MAssembly.cxx index 14f6e193..2c91b546 100644 --- a/src/MAssembly.cxx +++ b/src/MAssembly.cxx @@ -56,10 +56,10 @@ using namespace std; #include "MModule.h" #include "MGUIExpoCombinedViewer.h" #include "MModuleTransmitterRealta.h" - #include "MModuleLoaderSimulationsBalloon.h" #include "MModuleLoaderSimulationsSMEX.h" #include "MModuleLoaderMeasurementsROA.h" +#include "MModuleLoaderMeasurementsHDF.h" #include "MModuleReceiverBalloon.h" #include "MModuleLoaderMeasurementsBinary.h" #include "MModuleEnergyCalibration.h" @@ -74,6 +74,8 @@ using namespace std; #include "MModuleEventFilter.h" #include "MModuleEventSaver.h" #include "MModuleResponseGenerator.h" +#include "MModuleTACcut.h" +#include "MModuleNearestNeighbor.h" #include "MModuleDiagnostics.h" #include "MModuleDiagnosticsEnergyPerStrip.h" @@ -112,6 +114,7 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleLoaderSimulationsBalloon()); m_Supervisor->AddAvailableModule(new MModuleLoaderSimulationsSMEX()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsROA()); + m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsHDF()); m_Supervisor->AddAvailableModule(new MModuleReceiverBalloon()); m_Supervisor->AddAvailableModule(new MModuleLoaderMeasurementsBinary()); @@ -130,6 +133,8 @@ MAssembly::MAssembly() m_Supervisor->AddAvailableModule(new MModuleEventSaver()); m_Supervisor->AddAvailableModule(new MModuleTransmitterRealta()); m_Supervisor->AddAvailableModule(new MModuleResponseGenerator()); + m_Supervisor->AddAvailableModule(new MModuleTACcut()); + m_Supervisor->AddAvailableModule(new MModuleNearestNeighbor()); m_Supervisor->AddAvailableModule(new MModuleDiagnostics()); m_Supervisor->AddAvailableModule(new MModuleDiagnosticsEnergyPerStrip()); diff --git a/src/MBinaryFlightDataParser.cxx b/src/MBinaryFlightDataParser.cxx index 4be64951..ddc86fac 100644 --- a/src/MBinaryFlightDataParser.cxx +++ b/src/MBinaryFlightDataParser.cxx @@ -1242,16 +1242,16 @@ bool MBinaryFlightDataParser::ConvertToMReadOutAssemblys( dataframe * DataIn, ve NewEvent->SetTime( NewTime ); if( E.TrigAndVetoInfo & 0x70 ){ - NewEvent->SetVeto(); + NewEvent->SetShieldVeto(); // Was ->SetVeto() } if( E.TrigAndVetoInfo & 0x10 ){ - //had a guard ring 0 veto - NewEvent->SetGR0Veto(true); + //had a guard ring veto + NewEvent->SetGuardRingVeto(true); } if( E.TrigAndVetoInfo & 0x20 ){ - NewEvent->SetGR1Veto(true); + NewEvent->SetGuardRingVeto(true); } if( E.TrigAndVetoInfo & 0x40 ){ diff --git a/src/MGUIExpoTACcut.cxx b/src/MGUIExpoTACcut.cxx new file mode 100644 index 00000000..aa7ed9f4 --- /dev/null +++ b/src/MGUIExpoTACcut.cxx @@ -0,0 +1,282 @@ +/* + * MGUIExpoTACcut.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIExpoTACcut.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" + + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIExpoTACcut) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIExpoTACcut::MGUIExpoTACcut(MModule* Module) : MGUIExpo(Module) +{ + // standard constructor + + // Set the new title of the tab here: + m_TabTitle = "TAC Calibration"; + + // Set the histogram arrangment + // SetTACHistogramArrangement(1, 1); + + // use hierarchical cleaning + SetCleanup(kDeepCleanup); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIExpoTACcut::~MGUIExpoTACcut() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::Reset() +{ + //! Reset the data in the UI + + m_Mutex.Lock(); + for (auto H: m_TACHistograms) { + (H.second)->Reset(); + } + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::SetTACHistogramArrangement(const vector DetIDs) +{ + // Take in the list of detector IDs and determine the number in X and number in Y + // Update the variable m_DetectorMap. + m_Mutex.Lock(); + + unsigned int column = 0; + unsigned int row = 0; + + unsigned int max_columns = 4; + + unsigned int NDetectors = DetIDs.size(); + cout<<"MGUIExpoTACcut::SetTACHistogramArrangement: Number of detectors:"<< NDetectors< new_row; + m_DetectorMap.push_back(new_row); + column = 1; + } + + unsigned int DetID = DetIDs.at(i); + m_DetectorMap[row-1].push_back(DetID); + + TH1D* TAC = new TH1D("", "TAC", m_NBins[DetID], m_Min[DetID], m_Max[DetID]); + TAC->SetXTitle("Calibrated TAC (ns)"); + TAC->SetYTitle("counts"); + TAC->SetFillColor(kAzure+7); + + m_TACHistograms[DetID] = TAC; + // m_TACCanvases[DetID] = 0; + + ++column; + } + + if (NDetectors < max_columns) { + m_NColumns = NDetectors; + } else { + m_NColumns=max_columns; + } + + m_NRows = (NDetectors/max_columns) + 1; + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::SetTACHistogramParameters(unsigned int DetID, unsigned int NBins, double MinimumTAC, double MaximumTAC) +{ + // Set the energy histogram parameters + + m_Mutex.Lock(); + + m_NBins[DetID] = NBins; + m_Min[DetID] = MinimumTAC; + m_Max[DetID] = MaximumTAC; + TH1D* H = m_TACHistograms[DetID]; + H->SetBins(NBins, MinimumTAC, MaximumTAC); + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::SetTACHistogramName(unsigned int DetID, MString Name) +{ + // Set the title of the histogram + + m_Mutex.Lock(); + + if (m_TACHistograms.find(DetID) != m_TACHistograms.end()) { + m_TACHistograms[DetID]->SetTitle(Name); + } + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::AddTAC(unsigned int DetID, double TAC) +{ + // Add data to the energy histogram + + m_Mutex.Lock(); + + if (m_TACHistograms.find(DetID) != m_TACHistograms.end()) { + m_TACHistograms[DetID]->Fill(TAC); + } + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::Create() +{ + // Add the GUI options here + + // Do not create it twice! + if (m_IsCreated == true) return; + + m_Mutex.Lock(); + + TGLayoutHints* CanvasLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX | kLHintsExpandY, 2, 2, 2, 2); + + for (unsigned int y = 0; y < m_DetectorMap.size(); ++y) { + TGHorizontalFrame* HFrame = new TGHorizontalFrame(this); + AddFrame(HFrame, CanvasLayout); + + for (unsigned int x = 0; x < m_DetectorMap[y].size(); ++x) { + unsigned int DetID = m_DetectorMap[y][x]; + TRootEmbeddedCanvas* TACCanvas = new TRootEmbeddedCanvas("TAC", HFrame, 100, 100); + HFrame->AddFrame(TACCanvas, CanvasLayout); + m_TACCanvases[DetID] = TACCanvas; + + TACCanvas->GetCanvas()->cd(); + m_TACHistograms[DetID]->Draw("colz"); + TACCanvas->GetCanvas()->Update(); + } + } + + m_IsCreated = true; + + m_Mutex.UnLock(); +} + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::Update() +{ + //! Update the frame + + m_Mutex.Lock(); + + double Max = 0; + // for (auto H : m_TACHistograms) { + for (const auto dethistpair : m_TACHistograms) { + TH1D* H = dethistpair.second; + for (int bx = 2; bx < H->GetNbinsX(); ++bx) { // Skip first and last + if (Max < H->GetBinContent(bx)) { + Max = H->GetBinContent(bx); + } + } + } + Max *= 1.1; + for (const auto dethistpair : m_TACHistograms) { + TH1D* H = dethistpair.second; + H->SetMaximum(Max); + } + + for (auto C : m_TACCanvases) { + + (C.second)->GetCanvas()->Modified(); + (C.second)->GetCanvas()->Update(); + } + + m_Mutex.UnLock(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIExpoTACcut::Export(const MString& FileName) +{ + // Add data to the energy histogram + + m_Mutex.Lock(); + + TCanvas* P = new TCanvas(); + P->Divide(m_NColumns, m_NRows); + for (unsigned int y = 0; y < m_DetectorMap.size(); ++y) { + for (unsigned int x = 0; x < m_DetectorMap[y].size(); ++x) { + unsigned int DetID = m_DetectorMap[y][x]; + P->cd((x+1) + m_NColumns*y); + m_TACHistograms[DetID]->DrawCopy("colz"); + } + } + P->SaveAs(FileName); + delete P; + + m_Mutex.UnLock(); +} diff --git a/src/MGUIOptionsDepthCalibration2024.cxx b/src/MGUIOptionsDepthCalibration2024.cxx index 53688321..cd6c07ec 100644 --- a/src/MGUIOptionsDepthCalibration2024.cxx +++ b/src/MGUIOptionsDepthCalibration2024.cxx @@ -78,12 +78,6 @@ void MGUIOptionsDepthCalibration2024::Create() TGLayoutHints* Label2Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); m_OptionsFrame->AddFrame(m_SplinesFileSelector, Label2Layout); - m_TACCalFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a TAC Calibration file:", - dynamic_cast(m_Module)->GetTACCalFileName()); - m_TACCalFileSelector->SetFileType("TAC", "*.csv"); - TGLayoutHints* Label3Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - m_OptionsFrame->AddFrame(m_TACCalFileSelector, Label3Layout); - 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* Label4Layout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); @@ -133,7 +127,6 @@ bool MGUIOptionsDepthCalibration2024::OnApply() dynamic_cast(m_Module)->SetCoeffsFileName(m_CoeffsFileSelector->GetFileName()); dynamic_cast(m_Module)->SetSplinesFileName(m_SplinesFileSelector->GetFileName()); - dynamic_cast(m_Module)->SetTACCalFileName(m_TACCalFileSelector->GetFileName()); dynamic_cast(m_Module)->SetUCSDOverride(m_UCSDOverride->IsOn()); return true; diff --git a/src/MGUIOptionsEnergyCalibrationUniversal.cxx b/src/MGUIOptionsEnergyCalibrationUniversal.cxx index 6aafead3..d2ad6d72 100644 --- a/src/MGUIOptionsEnergyCalibrationUniversal.cxx +++ b/src/MGUIOptionsEnergyCalibrationUniversal.cxx @@ -105,8 +105,8 @@ bool MGUIOptionsEnergyCalibrationUniversal::ProcessMessage(long Message, long Pa { // Modify here if you have more buttons - bool Status = true; - + bool Status = true; + switch (GET_MSG(Message)) { case kC_COMMAND: switch (GET_SUBMSG(Message)) { @@ -123,7 +123,7 @@ bool MGUIOptionsEnergyCalibrationUniversal::ProcessMessage(long Message, long Pa m_TempFile->SetEnabled(false); } break; - } + } default: break; } @@ -146,17 +146,14 @@ bool MGUIOptionsEnergyCalibrationUniversal::ProcessMessage(long Message, long Pa bool MGUIOptionsEnergyCalibrationUniversal::OnApply() { - // Modify this to store the data in the module! + // Modify this to store the data in the module! dynamic_cast(m_Module)->SetFileName(m_FileSelector->GetFileName()); dynamic_cast(m_Module)->SetTempFileName(m_TempFile->GetFileName()); - if (dynamic_cast(m_Module)->GetPreampTempCorrection() != m_UseTempCal) dynamic_cast(m_Module)->EnablePreampTempCorrection(m_UseTempCal); - - - return true; + return true; } diff --git a/src/MGUIOptionsEventSaver.cxx b/src/MGUIOptionsEventSaver.cxx index 76020f67..dd821c0c 100644 --- a/src/MGUIOptionsEventSaver.cxx +++ b/src/MGUIOptionsEventSaver.cxx @@ -66,8 +66,9 @@ void MGUIOptionsEventSaver::Create() { PreCreate(); - TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); - + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft, 10, 10, 10, 10); + TGLayoutHints* RoaCheckButtonLayout = new TGLayoutHints(kLHintsTop | kLHintsLeft, 10, 10, 2, 2); + m_Mode = new MGUIERBList(m_OptionsFrame, "Please select an output mode:"); m_Mode->Add("*.roa file to use with melinator"); m_Mode->Add("*.dat file containing all information"); @@ -88,6 +89,10 @@ void MGUIOptionsEventSaver::Create() m_SaveBadEvents->SetOn(dynamic_cast(m_Module)->GetSaveBadEvents()); m_OptionsFrame->AddFrame(m_SaveBadEvents, LabelLayout); + m_SaveVetoEvents = new TGCheckButton(m_OptionsFrame, "Save guard ring and shield veto events (Veto)", 1); + m_SaveVetoEvents->SetOn(dynamic_cast(m_Module)->GetSaveVetoEvents()); + m_OptionsFrame->AddFrame(m_SaveVetoEvents, LabelLayout); + m_AddTimeTag = new TGCheckButton(m_OptionsFrame, "Add a unique time tag", 3); m_AddTimeTag->SetOn(dynamic_cast(m_Module)->GetAddTimeTag()); m_OptionsFrame->AddFrame(m_AddTimeTag, LabelLayout); @@ -103,8 +108,50 @@ void MGUIOptionsEventSaver::Create() dynamic_cast(m_Module)->GetSplitFileTime().GetAsSystemSeconds(), true, 0l); if (m_SplitFile->IsOn() == false) m_SplitFileTime->SetEnabled(false); m_OptionsFrame->AddFrame(m_SplitFileTime, SplitFileTimeLayout); - - + + TGLabel* ROAOptionsLabel = new TGLabel(m_OptionsFrame, "Special options for roa files:"); + m_OptionsFrame->AddFrame(ROAOptionsLabel, LabelLayout); + + m_RoaWithADCs = new TGCheckButton(m_OptionsFrame, "Include ADCs", 3); + m_RoaWithADCs->Associate(this); + m_RoaWithADCs->SetOn(dynamic_cast(m_Module)->GetRoaWithADCs()); + m_OptionsFrame->AddFrame(m_RoaWithADCs, RoaCheckButtonLayout); + + m_RoaWithTACs = new TGCheckButton(m_OptionsFrame, "Include TACs", 3); + m_RoaWithTACs->Associate(this); + m_RoaWithTACs->SetOn(dynamic_cast(m_Module)->GetRoaWithTACs()); + m_OptionsFrame->AddFrame(m_RoaWithTACs, RoaCheckButtonLayout); + + m_RoaWithEnergies = new TGCheckButton(m_OptionsFrame, "Include energies", 3); + m_RoaWithEnergies->Associate(this); + m_RoaWithEnergies->SetOn(dynamic_cast(m_Module)->GetRoaWithEnergies()); + m_OptionsFrame->AddFrame(m_RoaWithEnergies, RoaCheckButtonLayout); + + m_RoaWithTimings = new TGCheckButton(m_OptionsFrame, "Include timings", 3); + m_RoaWithTimings->Associate(this); + m_RoaWithTimings->SetOn(dynamic_cast(m_Module)->GetRoaWithTimings()); + m_OptionsFrame->AddFrame(m_RoaWithTimings, RoaCheckButtonLayout); + + m_RoaWithTemperatures = new TGCheckButton(m_OptionsFrame, "Include temperatures", 3); + m_RoaWithTemperatures->Associate(this); + m_RoaWithTemperatures->SetOn(dynamic_cast(m_Module)->GetRoaWithTemperatures()); + m_OptionsFrame->AddFrame(m_RoaWithTemperatures, RoaCheckButtonLayout); + + m_RoaWithFlags = new TGCheckButton(m_OptionsFrame, "Include flags", 3); + m_RoaWithFlags->Associate(this); + m_RoaWithFlags->SetOn(dynamic_cast(m_Module)->GetRoaWithFlags()); + m_OptionsFrame->AddFrame(m_RoaWithFlags, RoaCheckButtonLayout); + + m_RoaWithOrigins = new TGCheckButton(m_OptionsFrame, "Include origins", 3); + m_RoaWithOrigins->Associate(this); + m_RoaWithOrigins->SetOn(dynamic_cast(m_Module)->GetRoaWithOrigins()); + m_OptionsFrame->AddFrame(m_RoaWithOrigins, RoaCheckButtonLayout); + + m_RoaWithNearestNeighbors = new TGCheckButton(m_OptionsFrame, "Include nearest neighbor Hits", 3); + m_RoaWithNearestNeighbors->Associate(this); + m_RoaWithNearestNeighbors->SetOn(dynamic_cast(m_Module)->GetRoaWithNearestNeighbors()); + m_OptionsFrame->AddFrame(m_RoaWithNearestNeighbors, RoaCheckButtonLayout); + PostCreate(); } @@ -153,12 +200,24 @@ bool MGUIOptionsEventSaver::OnApply() // Modify this to store the data in the module! dynamic_cast(m_Module)->SetMode(m_Mode->GetSelected()); + dynamic_cast(m_Module)->SetFileName(m_FileSelector->GetFileName()); + dynamic_cast(m_Module)->SetSaveBadEvents(m_SaveBadEvents->IsOn()); + dynamic_cast(m_Module)->SetSaveVetoEvents(m_SaveVetoEvents->IsOn()); dynamic_cast(m_Module)->SetAddTimeTag(m_AddTimeTag->IsOn()); dynamic_cast(m_Module)->SetSplitFile(m_SplitFile->IsOn()); dynamic_cast(m_Module)->SetSplitFileTime(MTime(m_SplitFileTime->GetAsInt())); - + + dynamic_cast(m_Module)->SetRoaWithADCs(m_RoaWithADCs->IsOn()); + dynamic_cast(m_Module)->SetRoaWithTACs(m_RoaWithTACs->IsOn()); + dynamic_cast(m_Module)->SetRoaWithEnergies(m_RoaWithEnergies->IsOn()); + dynamic_cast(m_Module)->SetRoaWithTimings(m_RoaWithTimings->IsOn()); + dynamic_cast(m_Module)->SetRoaWithTemperatures(m_RoaWithTemperatures->IsOn()); + dynamic_cast(m_Module)->SetRoaWithFlags(m_RoaWithFlags->IsOn()); + dynamic_cast(m_Module)->SetRoaWithOrigins(m_RoaWithOrigins->IsOn()); + dynamic_cast(m_Module)->SetRoaWithNearestNeighbors(m_RoaWithNearestNeighbors->IsOn()); + return true; } diff --git a/src/MGUIOptionsLoaderMeasurements.cxx b/src/MGUIOptionsLoaderMeasurements.cxx index c9dfe019..79f38dde 100644 --- a/src/MGUIOptionsLoaderMeasurements.cxx +++ b/src/MGUIOptionsLoaderMeasurements.cxx @@ -43,8 +43,8 @@ ClassImp(MGUIOptionsLoaderMeasurements) //////////////////////////////////////////////////////////////////////////////// -MGUIOptionsLoaderMeasurements::MGUIOptionsLoaderMeasurements(MModule* Module) - : MGUIOptions(Module) +MGUIOptionsLoaderMeasurements::MGUIOptionsLoaderMeasurements(MModule* Module, MString FileType) + : MGUIOptions(Module), m_FileType(FileType) { // standard constructor } @@ -66,12 +66,12 @@ void MGUIOptionsLoaderMeasurements::Create() { PreCreate(); - m_FileSelector = new MGUIEFileSelector(m_OptionsFrame, "Please select a data file:", + m_FileSelector = new MGUIEFileSelector(m_OptionsFrame, MString("Please select a ") + m_FileType + " file:", dynamic_cast(m_Module)->GetFileName()); - m_FileSelector->SetFileType("Roa file", "*.roa"); - m_FileSelector->SetFileType("Roa file", "*.roa.gz"); - m_FileSelector->SetFileType("Data file", "*.dat"); - m_FileSelector->SetFileType("Data file", "*.dat.gz"); + m_FileSelector->SetFileType(m_FileType + " file", MString("*.") + m_FileType); + if (m_FileType == "roa") { + m_FileSelector->SetFileType(m_FileType + " file", MString("*.") + m_FileType + ".gz"); + } TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); m_OptionsFrame->AddFrame(m_FileSelector, LabelLayout); diff --git a/src/MGUIOptionsLoaderMeasurementsHDF.cxx b/src/MGUIOptionsLoaderMeasurementsHDF.cxx new file mode 100644 index 00000000..5126b5b8 --- /dev/null +++ b/src/MGUIOptionsLoaderMeasurementsHDF.cxx @@ -0,0 +1,141 @@ +/* + * MGUIOptionsLoaderMeasurementsHDF.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIOptionsLoaderMeasurementsHDF.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MModuleLoaderMeasurementsHDF.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsLoaderMeasurementsHDF) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsHDF::MGUIOptionsLoaderMeasurementsHDF(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsLoaderMeasurementsHDF::~MGUIOptionsLoaderMeasurementsHDF() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsLoaderMeasurementsHDF::Create() +{ + PreCreate(); + + TGLayoutHints* LabelLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + + m_FileSelectorHDF = new MGUIEFileSelector(m_OptionsFrame, "Please select a HDF5 file:", + dynamic_cast(m_Module)->GetFileName()); + m_FileSelectorHDF->SetFileType("HDF5 file", "*.hdf5"); + m_FileSelectorHDF->SetFileType("HDF5 file", "*.hdf"); + m_OptionsFrame->AddFrame(m_FileSelectorHDF, LabelLayout); + + + m_LoadContinuationFiles = new TGCheckButton(m_OptionsFrame, "Enable loading continuation HDF5 files", 1); + m_LoadContinuationFiles->SetOn(dynamic_cast(m_Module)->GetLoadContinuationFiles()); + m_LoadContinuationFiles->Associate(this); + m_OptionsFrame->AddFrame(m_LoadContinuationFiles, LabelLayout); + + + m_FileSelectorStripMap = new MGUIEFileSelector(m_OptionsFrame, "Please select a strip map file:", + dynamic_cast(m_Module)->GetFileNameStripMap()); + m_FileSelectorStripMap->SetFileType("Strip map file", "*.map"); + m_OptionsFrame->AddFrame(m_FileSelectorStripMap, LabelLayout); + + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsLoaderMeasurementsHDF::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + // Modify here if you have more buttons + + bool Status = true; + + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + default: + break; + } + break; + default: + break; + } + + if (Status == false) { + return false; + } + + // Call also base class + return MGUIOptions::ProcessMessage(Message, Parameter1, Parameter2); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsLoaderMeasurementsHDF::OnApply() +{ + // Modify this to store the data in the module! + + dynamic_cast(m_Module)->SetFileName(m_FileSelectorHDF->GetFileName()); + dynamic_cast(m_Module)->SetLoadContinuationFiles(m_LoadContinuationFiles->IsOn()); + dynamic_cast(m_Module)->SetFileNameStripMap(m_FileSelectorStripMap->GetFileName()); + + return true; +} + + +// MGUIOptionsLoaderMeasurementsHDF: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MGUIOptionsTACcut.cxx b/src/MGUIOptionsTACcut.cxx new file mode 100644 index 00000000..e410d65b --- /dev/null +++ b/src/MGUIOptionsTACcut.cxx @@ -0,0 +1,137 @@ +/* + * MGUIOptionsTACcut +.cxx + * + * + * Copyright (C) by Andreas Zoglauer + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Jau-Shian Liang. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +// Include the header: +#include "MGUIOptionsTACcut.h" + +// Standard libs: + +// ROOT libs: +#include +#include +#include +#include +#include + +// MEGAlib libs: +#include "MStreams.h" +#include "MString.h" +#include "MGUIEFileSelector.h" +#include "MGUIEMinMaxEntry.h" +#include "MGUIEEntry.h" + +// Nuclearizer libs: +#include "MModuleTACcut.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MGUIOptionsTACcut +) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsTACcut::MGUIOptionsTACcut(MModule* Module) + : MGUIOptions(Module) +{ + // standard constructor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MGUIOptionsTACcut::~MGUIOptionsTACcut() +{ + // kDeepCleanup is activated +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MGUIOptionsTACcut::Create() +{ + PreCreate(); + + m_TACCalFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a TAC Calibration file:", dynamic_cast(m_Module)->GetTACCalFileName()); + m_TACCalFileSelector->SetFileType("TAC", "*.csv"); + TGLayoutHints* TACCalLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_TACCalFileSelector, TACCalLayout); + + m_TACCutFileSelector = new MGUIEFileSelector(m_OptionsFrame, "Select a TAC Cut file:", dynamic_cast(m_Module)->GetTACCutFileName()); + m_TACCutFileSelector->SetFileType("TAC", "*.csv"); + TGLayoutHints* TACCutLayout = new TGLayoutHints(kLHintsTop | kLHintsCenterX | kLHintsExpandX, 10, 10, 10, 10); + m_OptionsFrame->AddFrame(m_TACCutFileSelector, TACCutLayout); + + PostCreate(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsTACcut::ProcessMessage(long Message, long Parameter1, long Parameter2) +{ + // Modify here if you have more buttons + + bool Status = true; + + switch (GET_MSG(Message)) { + case kC_COMMAND: + switch (GET_SUBMSG(Message)) { + case kCM_BUTTON: + break; + default: + break; + } + break; + default: + break; + } + + if (Status == false) { + return false; + } + + // Call also base class + return MGUIOptions::ProcessMessage(Message, Parameter1, Parameter2); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MGUIOptionsTACcut::OnApply() +{ + // Store the data in the module + dynamic_cast(m_Module)->SetTACCalFileName(m_TACCalFileSelector->GetFileName()); + dynamic_cast(m_Module)->SetTACCutFileName(m_TACCutFileSelector->GetFileName()); + + return true; +} + + +// MGUIOptionsTACcut: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleDepthCalibration2024.cxx b/src/MModuleDepthCalibration2024.cxx index 3dd095ee..135db6dc 100644 --- a/src/MModuleDepthCalibration2024.cxx +++ b/src/MModuleDepthCalibration2024.cxx @@ -63,6 +63,7 @@ MModuleDepthCalibration2024::MModuleDepthCalibration2024() : MModule() AddPreceedingModuleType(MAssembly::c_EventLoader, true); AddPreceedingModuleType(MAssembly::c_EnergyCalibration, true); AddPreceedingModuleType(MAssembly::c_StripPairing, true); + AddPreceedingModuleType(MAssembly::c_TACcut, true); // AddPreceedingModuleType(MAssembly::c_CrosstalkCorrection, false); // Soft requirement // Set all types this modules handles @@ -82,6 +83,8 @@ MModuleDepthCalibration2024::MModuleDepthCalibration2024() : MModule() m_AllowMultiThreading = true; m_AllowMultipleInstances = false; + m_Coeffs_Energy = 0; + m_NoError = 0; m_Error1 = 0; m_Error2 = 0; @@ -90,6 +93,8 @@ MModuleDepthCalibration2024::MModuleDepthCalibration2024() : MModule() m_Error5 = 0; m_Error6 = 0; m_ErrorSH = 0; + m_ErrorNullSH=0; + m_ErrorNoE=0; } @@ -108,58 +113,64 @@ MModuleDepthCalibration2024::~MModuleDepthCalibration2024() bool MModuleDepthCalibration2024::Initialize() { - if( LoadCoeffsFile(m_CoeffsFile) == false ){ - return false; - } - if( LoadSplinesFile(m_SplinesFile) == false ){ - return false; - } - // The detectors need to be in the same order as DetIDs. // ie DetID=0 should be the 0th detector in m_Detectors, DetID=1 should the 1st, etc. - m_Detectors = m_Geometry->GetDetectorList(); - - if( LoadTACCalFile(m_TACCalFile) == false ){ - cout << "No TAC Calibration file loaded. Proceeding without TAC Calibration." << endl; - } + vector DetList = m_Geometry->GetDetectorList(); // Look through the Geometry and get the names and thicknesses of all the detectors. - for(unsigned int i = 0; i < m_Detectors.size(); ++i){ + for(unsigned int i = 0; i < DetList.size(); ++i){ // For now, DetID is in order of detectors, which puts contraints on how the geometry file should be written. // If using the card cage at UCSD, default to DetID=11. unsigned int DetID = i; if ( m_UCSDOverride ){ DetID = 11; } - MDDetector* det = m_Detectors[i]; - // MString det_name = (det->GetDetectorVolume())->GetNamedDetectorName(0); - if (det->GetNNamedDetectors() > 0){ - // TODO: determine thickness of each detector using the geometry file - // cout << "Trying to get thickness from the geometry file..." << endl; - // cout << "step 1" << endl; - // MDVolume* vol = det->GetSensitiveVolume(0); - // cout << "step 2" << endl; - // MDShapeBRIK* shape = dynamic_cast(vol->GetShape()); - // cout << "step 3" << endl; - // double thickness = (shape->GetSize()).GetZ(); - // cout << "Success, the thickness is " << thickness << " cm" << endl; - // m_Thicknesses[DetID] = thickness; - MString det_name = det->GetNamedDetectorName(0); - m_DetectorNames[DetID] = det_name; - MDStrip3D* strip = dynamic_cast(det); - m_XPitches[DetID] = strip->GetPitchX(); - m_YPitches[DetID] = strip->GetPitchY(); - m_NXStrips[DetID] = strip->GetNStripsX(); - m_NYStrips[DetID] = strip->GetNStripsY(); - cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl; - cout << "Number of X strips: " << m_NXStrips[DetID] << endl; - cout << "Number of Y strips: " << m_NYStrips[DetID] << endl; - cout << "X strip pitch: " << m_XPitches[DetID] << endl; - cout << "Y strip pitch: " << m_YPitches[DetID] << endl; - m_DetectorIDs.push_back(DetID); + + MDDetector* det = DetList[i]; + vector DetectorNames; + if (det->GetTypeName() == "Strip3D") { + if (det->GetNSensitiveVolumes() == 1) { + MDVolume* vol = det->GetSensitiveVolume(0); + string det_name = vol->GetName().GetString(); + if (find(DetectorNames.begin(), DetectorNames.end(), det_name) == DetectorNames.end()) { + DetectorNames.push_back(det_name); + m_Thicknesses[DetID] = 2*(det->GetStructuralSize().GetZ()); + MDStrip3D* strip = dynamic_cast(det); + m_XPitches[DetID] = strip->GetPitchX(); + m_YPitches[DetID] = strip->GetPitchY(); + m_NXStrips[DetID] = strip->GetNStripsX(); + m_NYStrips[DetID] = strip->GetNStripsY(); + cout << "Found detector " << det_name << " corresponding to DetID=" << DetID << "." << endl; + cout << "Detector thickness: " << m_Thicknesses[DetID] << endl; + cout << "Number of X strips: " << m_NXStrips[DetID] << endl; + cout << "Number of Y strips: " << m_NYStrips[DetID] << endl; + cout << "X strip pitch: " << m_XPitches[DetID] << endl; + cout << "Y strip pitch: " << m_YPitches[DetID] << endl; + m_DetectorIDs.push_back(DetID); + m_Detectors[DetID] = det; + } else { + cout<<"ERROR in MModuleDepthCalibration2024::Initialize: Found a duplicate detector: "<GetNSensitiveVolumes()<<" Sensitive Volumes."<GetAvailableModuleByXmlTag("EnergyCalibrationUniversal"); if (m_EnergyCalibration == nullptr) { @@ -194,237 +205,183 @@ void MModuleDepthCalibration2024::CreateExpos() bool MModuleDepthCalibration2024::AnalyzeEvent(MReadOutAssembly* Event) { - for( unsigned int i = 0; i < Event->GetNHits(); ++i ){ - // Each event represents one photon. It contains Hits, representing interaction sites. - // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by - // instances of the MStripHit class. - MHit* H = Event->GetHit(i); - - int Grade = GetHitGrade(H); - - // cout << "got a hit with grade " << Grade << endl; - - // Handle different grades differently - // GRADE=-1 is an error. Break from the loop and continue. - if ( Grade == -1 ){ - H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); - ++m_ErrorSH; - } - - // GRADE=5 is some complicated geometry with multiple hits on a single strip. - // GRADE=4 means there are more than 2 strip hits on one or both sides. - else if( Grade > 3 ){ - H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); - if(Grade==4){ - ++m_Error4; - } - else if(Grade==5){ - ++m_Error5; - } - else if(Grade==6){ - ++m_Error6; - } - } + if (Event->GetGuardRingVeto()==true) { + + Event->SetDepthCalibrationIncomplete(); + return false; + + } else { + + for( unsigned int i = 0; i < Event->GetNHits(); ++i ){ + // Each event represents one photon. It contains Hits, representing interaction sites. + // H is a pointer to an instance of the MHit class. Each Hit has activated strips, represented by + // instances of the MStripHit class. + MHit* H = Event->GetHit(i); - // If the Grade is 0-3, we can handle it. - else { + int Grade = GetHitGrade(H); - MVector LocalPosition, PositionResolution, GlobalPosition, GlobalResolution, LocalOrigin; - - // 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 XStrips; - std::vector YStrips; - // cout << "looping over strip hits..." << endl; - for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ - // cout << "strip hit " << j << endl; - MStripHit* SH = H->GetStripHit(j); - if( SH->IsLowVoltageStrip() ) XStrips.push_back(SH); else YStrips.push_back(SH); - } - - // cout << "finished looping over strip hits" << endl; - double XEnergyFraction; - double YEnergyFraction; - MStripHit* XSH = GetDominantStrip(XStrips, XEnergyFraction); - MStripHit* YSH = GetDominantStrip(YStrips, YEnergyFraction); - - // cout << "found the dominant strips" << endl; - - double CTD_s = 0.0; - - //now try and get z position - int DetID = XSH->GetDetectorID(); - int XStripID = XSH->GetStripID(); - int YStripID = YSH->GetStripID(); - int pixel_code = 10000*DetID + 100*XStripID + YStripID; - - // TODO: Calculate X and Y positions more rigorously using charge sharing. - // Somewhat confusing notation: XStrips run parallel to X-axis, so we calculate X position with YStrips. - double Xpos = m_YPitches[DetID]*((m_NYStrips[DetID]/2.0) - ((double)YStripID)); - double Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)XStripID)); - // cout << "X position " << Xpos << endl; - // cout << "Y position " << Ypos << endl; - double Zpos = 0.0; - - double Xsigma = m_YPitches[DetID]/sqrt(12.0); - double Ysigma = m_XPitches[DetID]/sqrt(12.0); - double Zsigma = m_Thicknesses[DetID]/sqrt(12.0); - - // cout << "looking up the coefficients" << endl; - vector* Coeffs = GetPixelCoeffs(pixel_code); - - // TODO: For Card Cage, may need to add noise - double XTiming = XSH->GetTiming(); - double YTiming = YSH->GetTiming(); - if ( !m_UCSDOverride ) { - if ( m_TACCalFileIsLoaded ) { - if ( XSH->IsLowVoltageStrip() ){ - XTiming = XTiming*m_LVTACCal[DetID][XStripID][0] + m_LVTACCal[DetID][XStripID][1]; - YTiming = YTiming*m_HVTACCal[DetID][YStripID][0] + m_HVTACCal[DetID][YStripID][1]; - } - else { - XTiming = XTiming*m_HVTACCal[DetID][XStripID][0] + m_HVTACCal[DetID][XStripID][1]; - YTiming = YTiming*m_LVTACCal[DetID][YStripID][0] + m_LVTACCal[DetID][YStripID][1]; - } - } - else { - if ( XSH->IsLowVoltageStrip() ){ - XTiming = XTiming*0.405 - 525.; - YTiming = YTiming*0.43 - 500.; - } - else { - XTiming = XTiming*0.43 - 500.; - YTiming = YTiming*0.405 -525.; - } + // Handle different grades differently + // GRADE=-1 is an error. Break from the loop and continue. + if ( Grade < 0 ){ + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + if (Grade == -1) { + ++m_ErrorSH; + } else if (Grade == -2) { + ++m_ErrorNullSH; + } else if (Grade == -3) { + ++m_ErrorNoE; } - } - - // cout << "Got the coefficients: " << Coeffs << endl; - - // If there aren't coefficients loaded, then calibration is incomplete. - if( Coeffs == nullptr ){ - //set the bad flag for depth + } else if (Grade > 4) { // GRADE=5 is some complicated geometry with multiple hits on a single strip. GRADE=6 means not all strips are adjacent. 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( (XTiming < 1.0E-6) || (YTiming < 1.0E-6) ){ - // cout << "no timing info" << endl; - ++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 (Grade==5) { + ++m_Error5; + } else if (Grade==6) { + ++m_Error6; + } + } else { // If the Grade is 0-4, we can handle it. - if ( ctdvec.size() == 0){ - cout << "Empty CTD vector" << endl; - H->SetNoDepth(); - Event->SetDepthCalibrationIncomplete(); - } + MVector LocalPosition, PositionResolution, GlobalPosition, GlobalResolution, LocalOrigin; - double CTD; - if ( XSH->IsLowVoltageStrip() ){ - CTD = (YTiming - XTiming); - } - else { - CTD = (XTiming - YTiming); + // 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; + 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); } - // cout << "Got the CTD: " << CTD << endl; + double LVEnergyFraction; + double HVEnergyFraction; + MStripHit* LVSH = GetDominantStrip(LVStrips, LVEnergyFraction); + MStripHit* HVSH = GetDominantStrip(HVStrips, HVEnergyFraction); - // Confirmed that this matches SP's python code. - CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset + double CTD_s = 0.0; - // cout << "Transformed CTD: " << CTD_s << endl; + //now try and get z position + int DetID = LVSH->GetDetectorID(); + int LVStripID = LVSH->GetStripID(); + int HVStripID = HVSH->GetStripID(); + int PixelCode = 10000*DetID + 100*LVStripID + HVStripID; - double Xmin = * std::min_element(ctdvec.begin(), ctdvec.end()); - double Xmax = * std::max_element(ctdvec.begin(), ctdvec.end()); + // 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)); + double Ypos = m_XPitches[DetID]*((m_NXStrips[DetID]/2.0) - ((double)HVStripID)); + double Zpos = 0.0; - // cout << "Got the min and max ctd values: " << Xmin << "; " << Xmax << endl; + double Xsigma = m_YPitches[DetID]/sqrt(12.0); + double Ysigma = m_XPitches[DetID]/sqrt(12.0); + double Zsigma = m_Thicknesses[DetID]/sqrt(12.0); - double noise = GetTimingNoiseFWHM(pixel_code, H->GetEnergy()); + vector* Coeffs = GetPixelCoeffs(PixelCode); - // cout << "Got the timing noise: " << noise << endl; + // TODO: For Card Cage, may need to add noise + double LVTiming = LVSH->GetTiming(); + double HVTiming = HVSH->GetTiming(); - //if the CTD is out of range, check if we should reject the event. - if( (CTD_s < (Xmin - 2.0*noise)) || (CTD_s > (Xmax + 2.0*noise)) ){ + // If there aren't coefficients loaded, then calibration is incomplete. + if( Coeffs == nullptr ){ + //set the bad flag for depth H->SetNoDepth(); Event->SetDepthCalibrationIncomplete(); - ++m_Error2; + ++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) ){ + ++m_Error3; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); } - // If the CTD is in range, calculate the depth + // If there are coefficients and timing information is loaded, try calculating the CTD and depth else { - // cout << "Calculating depth" << endl; - // 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); - - // 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); - //cout << "summed probability: " << prob_sum << endl; - double weighted_depth = 0.0; - 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; k ctdvec = GetCTD(DetID, Grade); + vector depthvec = GetDepth(DetID); - Zsigma = sqrt(depth_var/prob_sum); - Zpos = mean_depth - (m_Thicknesses[DetID]/2.0); - // Zpos = mean_depth; - // cout << "calculated depth: " << Zpos << endl; + if ( ctdvec.size() == 0){ + cout << "Empty CTD vector" << endl; + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + } - // Add the depth to the GUI histogram. - m_ExpoDepthCalibration->AddDepth(DetID, Zpos); + double CTD = (HVTiming - LVTiming); - m_NoError+=1; - } - } + // Confirmed that this matches SP's python code. + CTD_s = (CTD - Coeffs->at(1))/(Coeffs->at(0)); //apply inverse stretch and offset - LocalPosition.SetXYZ(Xpos, Ypos, Zpos); - LocalOrigin.SetXYZ(0.0,0.0,0.0); - // cout << m_DetectorNames[DetID] << endl; - GlobalPosition = m_Geometry->GetGlobalPosition(LocalPosition, m_DetectorNames[DetID]); - // cout << "Found the GlobalPosition" << endl; + double Xmin = * std::min_element(ctdvec.begin(), ctdvec.end()); + double Xmax = * std::max_element(ctdvec.begin(), ctdvec.end()); - // Make sure XYZ resolution are correctly mapped to the global coord system. - PositionResolution.SetXYZ(Xsigma, Ysigma, Zsigma); - GlobalResolution = (m_Geometry->GetGlobalPosition(PositionResolution, m_DetectorNames[DetID]) - m_Geometry->GetGlobalPosition(LocalOrigin, m_DetectorNames[DetID])).Abs(); - - // cout << "Set the PositionResolution vector" << endl; + double noise = GetTimingNoiseFWHM(PixelCode, H->GetEnergy()); - H->SetPosition(GlobalPosition); + //if the CTD is out of range, check if we should reject the event. + if( (CTD_s < (Xmin - 2.0*noise)) || (CTD_s > (Xmax + 2.0*noise)) ){ + H->SetNoDepth(); + Event->SetDepthCalibrationIncomplete(); + ++m_Error2; + } - // cout << "Set the global position for the strip hit" << endl; + // If the CTD is in range, calculate the depth + else { + // 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); + + // 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]; + } + // 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) { + if (HasExpos() == true) { + m_ExpoDepthCalibration->AddDepth(DetID, Zpos); + } + } + m_NoError+=1; + } + } + + LocalPosition.SetXYZ(Xpos, Ypos, Zpos); + LocalOrigin.SetXYZ(0.0,0.0,0.0); + GlobalPosition = m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalPosition); - H->SetPositionResolution(GlobalResolution); + // Make sure XYZ resolution are correctly mapped to the global coord system. + PositionResolution.SetXYZ(Xsigma, Ysigma, Zsigma); + GlobalResolution = ((m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(PositionResolution)) - (m_Detectors[DetID]->GetSensitiveVolume(0)->GetPositionInWorldVolume(LocalOrigin))).Abs(); + + H->SetPosition(GlobalPosition); - // cout << "Set the position resolution for the strip hit" << endl; + H->SetPositionResolution(GlobalResolution); + + } } } - + Event->SetAnalysisProgress(MAssembly::c_DepthCorrection | MAssembly::c_PositionDetermiation); return true; @@ -453,14 +410,14 @@ MStripHit* MModuleDepthCalibration2024::GetDominantStrip(vector& Str return MaxStrip; } -double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int pixel_code, double Energy) +double MModuleDepthCalibration2024::GetTimingNoiseFWHM(int PixelCode, double Energy) { // Placeholder for determining the timing noise with energy, and possibly even on a pixel-by-pixel basis. // Should follow 1/E relation // TODO: Determine real energy dependence and implement it here. double noiseFWHM = 0.0; - if ( m_Coeffs_Energy != NULL ){ - noiseFWHM = m_Coeffs[pixel_code][2] * m_Coeffs_Energy/Energy; + if ( m_Coeffs_Energy != 0 ){ + noiseFWHM = m_Coeffs[PixelCode][2] * m_Coeffs_Energy/Energy; if ( noiseFWHM < 3.0*2.355 ){ noiseFWHM = 3.0*2.355; } @@ -476,10 +433,10 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) // Read in the stretch and offset file, which should have a header line with information on the measurements: // ### 800 V 80 K 59.5 keV // And which should contain for each pixel: - // Pixel code (10000*det + 100*Xchannel + Ychannel), Stretch, Offset, Timing/CTD noise, Chi2 for the CTD fit (for diagnostics mainly) + // Pixel code (10000*det + 100*LVStrip + HVStrip), Stretch, Offset, Timing/CTD noise, Chi2 for the CTD fit (for diagnostics mainly) MFile F; if( F.Open(FName) == false ){ - cout << "MModuleDepthCalibration2024: failed to open coefficients file..." << endl; + cout << "ERROR in MModuleDepthCalibration2024::LoadCoeffsFile: failed to open coefficients file." << endl; return false; } else { MString Line; @@ -492,7 +449,7 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) else { std::vector Tokens = Line.Tokenize(","); if( Tokens.size() == 5 ){ - int pixel_code = Tokens[0].ToInt(); + int PixelCode = Tokens[0].ToInt(); double Stretch = Tokens[1].ToDouble(); double Offset = Tokens[2].ToDouble(); double CTD_FWHM = Tokens[3].ToDouble() * 2.355; @@ -500,75 +457,25 @@ bool MModuleDepthCalibration2024::LoadCoeffsFile(MString FName) // Previous iteration of depth calibration read in "Scale" instead of ctd resolution. vector coeffs; coeffs.push_back(Stretch); coeffs.push_back(Offset); coeffs.push_back(CTD_FWHM); coeffs.push_back(Chi2); - m_Coeffs[pixel_code] = coeffs; + m_Coeffs[PixelCode] = coeffs; } } } F.Close(); } - m_CoeffsFileIsLoaded = true; - return true; } -bool MModuleDepthCalibration2024::LoadTACCalFile(MString FName) -{ - // Read in the TAC Calibration file, which should contain for each strip: - // DetID, h or l for high or low voltage, TAC cal, TAC cal error, TAC cal offset, TAC offset error - MFile F; - if( F.Open(FName) == false ){ - cout << "MModuleDepthCalibration2024: failed to open TAC Calibration file." << endl; - m_TACCalFileIsLoaded = false; - return false; - } else { - for(unsigned int i = 0; i < m_Detectors.size(); ++i){ - unordered_map> temp_map_HV; - m_HVTACCal[i] = temp_map_HV; - unordered_map> temp_map_LV; - m_LVTACCal[i] = temp_map_LV; - } - MString Line; - while( F.ReadLine( Line ) ){ - if( !Line.BeginsWith("#") ){ - std::vector Tokens = Line.Tokenize(","); - if( Tokens.size() == 7 ){ - int DetID = Tokens[0].ToInt(); - int StripID = Tokens[2].ToInt(); - double taccal = Tokens[3].ToDouble(); - double taccal_err = Tokens[4].ToDouble(); - double offset = Tokens[5].ToDouble(); - double offset_err = Tokens[6].ToDouble(); - vector cal_vals; - cal_vals.push_back(taccal); cal_vals.push_back(offset); cal_vals.push_back(taccal_err); cal_vals.push_back(offset_err); - if ( Tokens[1] == "l" ){ - m_LVTACCal[DetID][StripID] = cal_vals; - } - else if ( Tokens[1] == "h" ){ - m_HVTACCal[DetID][StripID] = cal_vals; - } - } - } - } - F.Close(); - } - - m_TACCalFileIsLoaded = true; - - return true; - -} - - -std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int pixel_code) +std::vector* MModuleDepthCalibration2024::GetPixelCoeffs(int PixelCode) { // Check to see if the stretch and offset have been loaded. If so, try to get the coefficients for the specified pixel. if( m_CoeffsFileIsLoaded ){ - if( m_Coeffs.count(pixel_code) > 0 ){ - return &m_Coeffs[pixel_code]; + if( m_Coeffs.count(PixelCode) > 0 ){ + return &m_Coeffs[PixelCode]; } else { - cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << pixel_code << " not found." << endl; + cout << "MModuleDepthCalibration2024::GetPixelCoeffs: cannot get stretch and offset; pixel code " << PixelCode << " not found." << endl; return nullptr; } } else { @@ -599,6 +506,7 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) // '' '' '' 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; @@ -608,45 +516,44 @@ bool MModuleDepthCalibration2024::LoadSplinesFile(MString FName) vector temp_vec; ctdarr.push_back(temp_vec); } + bool Result = true; MString line; - int DetID, NewDetID; - while( F.ReadLine(line) ){ - if( line.Length() != 0 ){ - if( line.BeginsWith("#") ){ + int DetID = 0; + while (F.ReadLine(line)) { + if (line.Length() != 0) { + 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(" "); - NewDetID = tokens[1].ToInt(); - if( depthvec.size() > 0 ) { - AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); + 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 ){ + for (unsigned int i=0; i < 5; ++i) { vector temp_vec; ctdarr.push_back(temp_vec); } - DetID = NewDetID; + DetID = tokens[1].ToInt(); } else { vector tokens = line.Tokenize(","); depthvec.push_back(tokens[0].ToDouble()); // Multiple CTDs allowed. - for( unsigned int i = 0; i < (tokens.size() - 1); ++i ){ + 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){ + for (unsigned int i=tokens.size()-1; i<5; ++i) { ctdarr[i].push_back(tokens[1].ToDouble()); } } } } //make last spline - if( depthvec.size() > 0 ){ - AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); + if (depthvec.size() > 0) { + Result &= AddDepthCTD(depthvec, ctdarr, DetID, m_DepthGrid, m_CTDMap); } - m_SplinesFileIsLoaded = true; - return true; + return Result; } @@ -666,21 +573,21 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ } // Take a Hit and separate its activated p and n strips into separate vectors. - std::vector PStrips; - std::vector NStrips; - vector PStripIDs; - vector NStripIDs; + std::vector LVStrips; + std::vector HVStrips; + vector LVStripIDs; + vector HVStripIDs; for( unsigned int j = 0; j < H->GetNStripHits(); ++j){ MStripHit* SH = H->GetStripHit(j); if( SH == NULL ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got NULL strip hit :( " << endl; return -1;} if( SH->GetEnergy() == 0 ) { cout << "ERROR in MModuleDepthCalibration2024: Depth Calibration: got strip without energy :( " << endl; return -1;} if( SH->IsLowVoltageStrip() ){ - PStrips.push_back(SH); - PStripIDs.push_back(SH->GetStripID()); + LVStrips.push_back(SH); + LVStripIDs.push_back(SH->GetStripID()); } else { - NStrips.push_back(SH); - NStripIDs.push_back(SH->GetStripID()); + HVStrips.push_back(SH); + HVStripIDs.push_back(SH->GetStripID()); } } @@ -691,15 +598,15 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ return 5; } - if( PStrips.size()>0 && NStrips.size()>0 ){ - int Nmin = * std::min_element(NStripIDs.begin(), NStripIDs.end()); - int Nmax = * std::max_element(NStripIDs.begin(), NStripIDs.end()); + if( LVStrips.size()>0 && HVStrips.size()>0 ){ + int HVmin = * std::min_element(HVStripIDs.begin(), HVStripIDs.end()); + int HVmax = * std::max_element(HVStripIDs.begin(), HVStripIDs.end()); - int Pmin = * std::min_element(PStripIDs.begin(), PStripIDs.end()); - int Pmax = * std::max_element(PStripIDs.begin(), PStripIDs.end()); + int LVmin = * std::min_element(LVStripIDs.begin(), LVStripIDs.end()); + int LVmax = * std::max_element(LVStripIDs.begin(), LVStripIDs.end()); // If the strip hits are not all adjacent, it's a bad grade. - if ( ((Nmax - Nmin) >= (NStrips.size())) || ((Pmax - Pmin) >= (PStrips.size())) ){ + if ( ((HVmax - HVmin) >= (HVStrips.size())) || ((LVmax - LVmin) >= (LVStrips.size())) ){ return 6; } } @@ -711,48 +618,48 @@ int MModuleDepthCalibration2024::GetHitGrade(MHit* H){ int return_value; // If 1 strip on each side, GRADE=0 // This represents the center of the pixel - if( (PStrips.size() == 1) && (NStrips.size() == 1) ){ + if( ((LVStrips.size() == 1) && (HVStrips.size() == 1)) || ((LVStrips.size() == 3) && (HVStrips.size() == 3)) ){ return_value = 0; } // If 2 hits on N side and 1 on P, GRADE=1 // This represents the middle of the edges of the pixel - else if( (PStrips.size() == 1) && (NStrips.size() == 2) ){ + else if( (LVStrips.size() == 1) && (HVStrips.size() == 2) ){ return_value = 1; } // If 2 hits on P and 1 on N, GRADE=2 // This represents the middle of the edges of the pixel - else if( (PStrips.size() == 2) && (NStrips.size() == 1) ){ + else if( (LVStrips.size() == 2) && (HVStrips.size() == 1) ){ return_value = 2; } // If 2 strip hits on both sides, GRADE=3 // This represents the corners the pixel - else if( (PStrips.size() == 2) && (NStrips.size() == 2) ){ + else if( (LVStrips.size() == 2) && (HVStrips.size() == 2) ){ return_value = 3; } // If 3 hits on N side and 1 on P, GRADE=0 // This represents the middle of the pixel, near the p (LV) side of the detector. - else if( (PStrips.size() == 1) && (NStrips.size() == 3) ){ + else if( (LVStrips.size() == 1) && (HVStrips.size() == 3) ){ return_value = 0; } // If 3 hits on P and 1 on N, GRADE=0 // This represents the middle of the pixel, near the n (HV) side of the detector. - else if( (PStrips.size() == 3) && (NStrips.size() == 1) ){ + else if( (LVStrips.size() == 3) && (HVStrips.size() == 1) ){ return_value = 0; } // If 3 hits on N side and 2 on P, GRADE=0 // This represents the middle of the edge of the pixel, near the p (LV) side of the detector. - else if( (PStrips.size() == 2) && (NStrips.size() == 3) ){ + else if( (LVStrips.size() == 2) && (HVStrips.size() == 3) ){ return_value = 2; } // If 3 hits on P and 2 on N, GRADE=0 // This represents the middle of the edge of the pixel, near the n (HV) side of the detector. - else if( (PStrips.size() == 3) && (NStrips.size() == 2) ){ + else if( (LVStrips.size() == 3) && (HVStrips.size() == 2) ){ return_value = 1; } @@ -778,21 +685,24 @@ bool MModuleDepthCalibration2024::AddDepthCTD(vector depthvec, vector 0) ){ - cout << "MModuleDepthCalibration2024::AddDepthCTD: The number of values in the CTD list is not equal to the number of depth values." << endl; + for (unsigned int i = 0; i < ctdarr.size(); ++i) { + if ((ctdarr[i].size() != depthvec.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) { + cout<<"ERROR in MModuleDepthCalibration2024::AddDepthCTD: The thickness of detector "<GetValue(); } - MXmlNode* TACCalFileNameNode = Node->GetNode("TACCalFileName"); - if (TACCalFileNameNode != 0) { - m_TACCalFile = TACCalFileNameNode->GetValue(); + MXmlNode* UCSDOverrideNode = Node->GetNode("UCSDOverride"); + if( UCSDOverrideNode != NULL ){ + m_UCSDOverride = (bool) UCSDOverrideNode->GetValueAsInt(); } return true; @@ -888,8 +798,8 @@ MXmlNode* MModuleDepthCalibration2024::CreateXmlConfiguration() MXmlNode* Node = new MXmlNode(0,m_XmlTag); new MXmlNode(Node, "CoeffsFileName", m_CoeffsFile); new MXmlNode(Node, "SplinesFileName", m_SplinesFile); - new MXmlNode(Node, "TACCalFileName", m_TACCalFile); - + new MXmlNode(Node, "UCSDOverride",(unsigned int) m_UCSDOverride); + return Node; } @@ -908,6 +818,8 @@ void MModuleDepthCalibration2024::Finalize() cout << "Number of hits with non-adjacent strip hits: " << m_Error6 << endl; cout << "Number of hits with too many strip hits: " << m_Error4 << endl; 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 ); diff --git a/src/MModuleEnergyCalibrationUniversal.cxx b/src/MModuleEnergyCalibrationUniversal.cxx index 3180cba1..8d705851 100644 --- a/src/MModuleEnergyCalibrationUniversal.cxx +++ b/src/MModuleEnergyCalibrationUniversal.cxx @@ -74,12 +74,13 @@ MModuleEnergyCalibrationUniversal::MModuleEnergyCalibrationUniversal() : MModule // Set all modules, which have to be done before this module AddPreceedingModuleType(MAssembly::c_EventLoader); + // AddPreceedingModuleType(MAssembly::c_TACcut); // Set all types this modules handles AddModuleType(MAssembly::c_EnergyCalibration); // Set all modules, which can follow this module - AddSucceedingModuleType(MAssembly::c_NoRestriction); + AddSucceedingModuleType(MAssembly::c_TACcut); // Set if this module has an options GUI m_HasOptionsGUI = true; @@ -394,18 +395,6 @@ bool MModuleEnergyCalibrationUniversal::AnalyzeEvent(MReadOutAssembly* Event) } } - for (unsigned int i = 0; i < Event->GetNStripHits(); ) { - MStripHit* SH = Event->GetStripHit(i); - if (SH->GetEnergy() < 8 || SH->GetTiming() < 8700 || SH->GetTiming() > 12000) { - // cout<<"HACK: Removing strip hit due to TAC "<GetTiming()<<" cut or energy "<GetEnergy()<RemoveStripHit(i); - delete SH; - } else { - ++i; - } - } - - Event->SetAnalysisProgress(MAssembly::c_EnergyCalibration); return true; diff --git a/src/MModuleEventSaver.cxx b/src/MModuleEventSaver.cxx index 8bace658..b7070b7a 100644 --- a/src/MModuleEventSaver.cxx +++ b/src/MModuleEventSaver.cxx @@ -75,8 +75,18 @@ MModuleEventSaver::MModuleEventSaver() : MModule() m_InternalFileName = ""; m_Zip = false; m_SaveBadEvents = true; + m_SaveVetoEvents = true; m_AddTimeTag = false; - + + m_RoaWithADCs = true; + m_RoaWithTACs = true; + m_RoaWithEnergies = false; + m_RoaWithTimings = false; + m_RoaWithTemperatures = false; + m_RoaWithFlags = false; + m_RoaWithOrigins = false; + m_RoaWithNearestNeighbors = true; + m_SplitFile = true; m_SplitFileTime.Set(60*10); // seconds m_SubFileStart.Set(0); @@ -163,14 +173,49 @@ bool MModuleEventSaver::Initialize() Header<= c_Error) cout<IsBad() == true) return true; } + if (m_SaveVetoEvents == false) { + if (Event->IsVeto() == true) return true; + } + + MFile* Choosen = 0; // Wish C++ would allow unassigned references... if (m_SplitFile == true) { MTime Current = Event->GetTime(); @@ -302,14 +352,14 @@ bool MModuleEventSaver::AnalyzeEvent(MReadOutAssembly* Event) } else { Choosen = &m_Out; } - + ostringstream Out; if (m_Mode == c_EvtaFile) { Event->StreamEvta(Out); } else if (m_Mode == c_DatFile) { Event->StreamDat(Out, 1); } else if (m_Mode == c_RoaFile) { - Event->StreamRoa(Out); + Event->StreamRoa(Out, m_RoaWithADCs, m_RoaWithTACs, m_RoaWithEnergies, m_RoaWithTimings, m_RoaWithTemperatures, m_RoaWithFlags, m_RoaWithOrigins, m_RoaWithNearestNeighbors); } Choosen->Write(Out); @@ -351,6 +401,10 @@ bool MModuleEventSaver::ReadXmlConfiguration(MXmlNode* Node) if (SaveBadEventsNode != 0) { m_SaveBadEvents = SaveBadEventsNode->GetValueAsBoolean(); } + MXmlNode* SaveVetoEventsNode = Node->GetNode("SaveVetoEvents"); + if (SaveVetoEventsNode != 0) { + m_SaveVetoEvents = SaveVetoEventsNode->GetValueAsBoolean(); + } MXmlNode* AddTimeTagNode = Node->GetNode("AddTimeTag"); if (AddTimeTagNode != 0) { m_AddTimeTag = AddTimeTagNode->GetValueAsBoolean(); @@ -364,6 +418,39 @@ bool MModuleEventSaver::ReadXmlConfiguration(MXmlNode* Node) m_SplitFileTime.Set(SplitFileTimeNode->GetValueAsInt()); } + MXmlNode* RoaWithADCsNode = Node->GetNode("RoaWithADCs"); + if (RoaWithADCsNode != nullptr) { + m_RoaWithADCs = RoaWithADCsNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithTACsNode = Node->GetNode("RoaWithTACs"); + if (RoaWithTACsNode != nullptr) { + m_RoaWithTACs = RoaWithTACsNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithEnergiesNode = Node->GetNode("RoaWithEnergies"); + if (RoaWithEnergiesNode != nullptr) { + m_RoaWithEnergies = RoaWithEnergiesNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithTimingsNode = Node->GetNode("RoaWithTimings"); + if (RoaWithTimingsNode != nullptr) { + m_RoaWithTimings = RoaWithTimingsNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithTemperaturesNode = Node->GetNode("RoaWithTemperatures"); + if (RoaWithTemperaturesNode != nullptr) { + m_RoaWithTemperatures = RoaWithTemperaturesNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithFlagsNode = Node->GetNode("RoaWithFlags"); + if (RoaWithFlagsNode != nullptr) { + m_RoaWithFlags = RoaWithFlagsNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithOriginsNode = Node->GetNode("RoaWithOrigins"); + if (RoaWithOriginsNode != nullptr) { + m_RoaWithOrigins = RoaWithOriginsNode->GetValueAsBoolean(); + } + MXmlNode* RoaWithNearestNeighborsNode = Node->GetNode("RoaWithNearestNeighbors"); + if (RoaWithNearestNeighborsNode != nullptr) { + m_RoaWithNearestNeighbors = RoaWithNearestNeighborsNode->GetValueAsBoolean(); + } + return true; } @@ -379,9 +466,17 @@ MXmlNode* MModuleEventSaver::CreateXmlConfiguration() new MXmlNode(Node, "FileName", m_FileName); new MXmlNode(Node, "Mode", m_Mode); new MXmlNode(Node, "SaveBadEvents", m_SaveBadEvents); + new MXmlNode(Node, "SaveVetoEvents", m_SaveVetoEvents); new MXmlNode(Node, "AddTimeTag", m_AddTimeTag); new MXmlNode(Node, "SplitFile", m_SplitFile); new MXmlNode(Node, "SplitFileTime", m_SplitFileTime.GetAsSystemSeconds()); + new MXmlNode(Node, "RoaWithADCs", m_RoaWithADCs); + new MXmlNode(Node, "RoaWithTACs", m_RoaWithTACs); + new MXmlNode(Node, "RoaWithEnergies", m_RoaWithEnergies); + new MXmlNode(Node, "RoaWithTimings", m_RoaWithTimings); + new MXmlNode(Node, "RoaWithFlags", m_RoaWithFlags); + new MXmlNode(Node, "RoaWithOrigins", m_RoaWithOrigins); + new MXmlNode(Node, "RoaWithNearestNeighbors", m_RoaWithNearestNeighbors); return Node; } diff --git a/src/MModuleLoaderMeasurements.cxx b/src/MModuleLoaderMeasurements.cxx index 5c6b890c..acc133bb 100644 --- a/src/MModuleLoaderMeasurements.cxx +++ b/src/MModuleLoaderMeasurements.cxx @@ -95,18 +95,5 @@ bool MModuleLoaderMeasurements::AnalyzeEvent(MReadOutAssembly* Event) } -/////////////////////////////////////////////////////////////////////////////// - - -void MModuleLoaderMeasurements::ShowOptionsGUI() -{ - //! Show the options GUI - - MGUIOptionsLoaderMeasurements* Options = new MGUIOptionsLoaderMeasurements(this); - Options->Create(); - gClient->WaitForUnmap(Options); -} - - // MModuleLoaderMeasurements.cxx: the end... //////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleLoaderMeasurementsBinary.cxx b/src/MModuleLoaderMeasurementsBinary.cxx index af547b39..4d0e0343 100644 --- a/src/MModuleLoaderMeasurementsBinary.cxx +++ b/src/MModuleLoaderMeasurementsBinary.cxx @@ -1,19 +1,19 @@ /* - * MModuleLoaderMeasurementsBinary.cxx - * - * - * Copyright (C) by Alex Lowell & Andreas Zoglauer. - * All rights reserved. - * - * - * This code implementation is the intellectual property of - * Andreas Zoglauer. - * - * By copying, distributing or modifying the Program (or any work - * based on the Program) you indicate your acceptance of this statement, - * and all its terms. - * - */ +* MModuleLoaderMeasurementsBinary.cxx +* +* +* Copyright (C) by Alex Lowell & Andreas Zoglauer. +* All rights reserved. +* +* +* This code implementation is the intellectual property of +* Andreas Zoglauer. +* +* By copying, distributing or modifying the Program (or any work +* based on the Program) you indicate your acceptance of this statement, +* and all its terms. +* +*/ @@ -84,11 +84,11 @@ MModuleLoaderMeasurementsBinary::MModuleLoaderMeasurementsBinary() : MModule(), m_IgnoreAspect = false; //this was set to true and was causing events to be pushed through the pipeline before aspect info was available for them AWL Sep 20 2016 m_FileIsDone = false; - - m_IsZipped = false; + + m_IsZipped = false; m_ZipFile = NULL; - m_ExpoAspectViewer = nullptr; + m_ExpoAspectViewer = nullptr; } @@ -106,13 +106,13 @@ MModuleLoaderMeasurementsBinary::~MModuleLoaderMeasurementsBinary() void MModuleLoaderMeasurementsBinary::CreateExpos() { - // Create all expos - - if (HasExpos() == true) return; - - // Set the histogram display - m_ExpoAspectViewer = new MGUIExpoAspectViewer(this); - m_Expos.push_back(m_ExpoAspectViewer); + // Create all expos + + if (HasExpos() == true) return; + + // Set the histogram display + m_ExpoAspectViewer = new MGUIExpoAspectViewer(this); + m_Expos.push_back(m_ExpoAspectViewer); } @@ -122,35 +122,35 @@ FILE * f_TOnly; bool MModuleLoaderMeasurementsBinary::OpenNextFile() { - //! Open next file, return false on error - - ++m_OpenFileID; - if (m_OpenFileID >= (int) m_BinaryFileNames.size()) return false; - - m_IsZipped = m_BinaryFileNames[m_OpenFileID].EndsWith(".gz"); - - if (m_IsZipped == false) { - if (m_In.is_open()) m_In.close(); - m_In.clear(); - - m_In.open(m_BinaryFileNames[m_OpenFileID], ios::binary); - if (m_In.is_open() == false) { - if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Info) cout<= (int) m_BinaryFileNames.size()) return false; + + m_IsZipped = m_BinaryFileNames[m_OpenFileID].EndsWith(".gz"); + + if (m_IsZipped == false) { + if (m_In.is_open()) m_In.close(); + m_In.clear(); + + m_In.open(m_BinaryFileNames[m_OpenFileID], ios::binary); + if (m_In.is_open() == false) { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Info) cout<= c_Error) cout<= c_Error) cout<= c_Info) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Info) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout< Stream(Size); // Check if we reached the end of the file, if yes, truncate, and set the OK flag to false - // when the end of the file is reached, we want to + // when the end of the file is reached, we want to - //AWL restructured this so that we don't allocate/fill a 1MB array when there is nothing to read. + //AWL restructured this so that we don't allocate/fill a 1MB array when there is nothing to read. vector Stream; unsigned int Size = 1000000; @@ -261,55 +261,55 @@ bool MModuleLoaderMeasurementsBinary::IsReady() Read = 0; } else { Stream.reserve(Size); - if (m_IsZipped == false) { - m_In.read(&Stream[0], Size); - Read = m_In.gcount(); - } else { - Read = 0; - for (unsigned int i = 0; i < Size; ++i) { - int c = gzgetc(m_ZipFile); - if (c == -1) { - break; - } - Stream[i] = (char) c; - Read = i+1; - } - } + if (m_IsZipped == false) { + m_In.read(&Stream[0], Size); + Read = m_In.gcount(); + } else { + Read = 0; + for (unsigned int i = 0; i < Size; ++i) { + int c = gzgetc(m_ZipFile); + if (c == -1) { + break; + } + Stream[i] = (char) c; + Read = i+1; + } + } } // If we do not read anything, try again with the next file - if (Read == 0) { - if (OpenNextFile() == true) { - Stream.reserve(Size); - if (m_IsZipped == false) { - m_In.read(&Stream[0], Size); - Read = m_In.gcount(); - } else { - Read = 0; - for (unsigned int i = 0; i < Size; ++i) { - int c = gzgetc(m_ZipFile); - if (c == -1) { - break; - } - Stream[i] = (char) c; - Read = i+1; - } - } - } - } + if (Read == 0) { + if (OpenNextFile() == true) { + Stream.reserve(Size); + if (m_IsZipped == false) { + m_In.read(&Stream[0], Size); + Read = m_In.gcount(); + } else { + Read = 0; + for (unsigned int i = 0; i < Size; ++i) { + int c = gzgetc(m_ZipFile); + if (c == -1) { + break; + } + Stream[i] = (char) c; + Read = i+1; + } + } + } + } /* - if (Read < Size) { - m_IsOK = false; - } - */ + * i f (Read < Size) { * + * m_IsOK = false; +} +*/ /* - if (Read < Size) { - m_FileIsDone = true; - SetIsDone(true); - } - */ + * i f (Read < Size) { * + * m_FileIsDone = true; + * SetIsDone(true); +} +*/ if (Read == 0) { m_FileIsDone = true; @@ -330,7 +330,7 @@ bool MModuleLoaderMeasurementsBinary::IsReady() } //cout<<"Received: "<GetCL() < LastCL){ - cout << LastCL << "--->" << NewEvent->GetCL() << endl; - } - LastCL = NewEvent->GetCL(); - */ + * i f(NewEvent->GetCL() < LastCL){ * + * cout << LastCL << "--->" << NewEvent->GetCL() << endl; +} +LastCL = NewEvent->GetCL(); +*/ //print TOnly info for these events /* @@ -424,7 +424,7 @@ bool MModuleLoaderMeasurementsBinary::AnalyzeEvent(MReadOutAssembly* Event) } //cout<<"Adding: "<GetTime()<<":"<GetHeading()<AddHeading(NewEvent->GetTime(), A->GetHeading(), A->GetGPS_or_magnetometer(), A->GetBRMS(), A->GetAttFlag()); + m_ExpoAspectViewer->AddHeading(NewEvent->GetTime(), A->GetHeading(), A->GetGPS_or_magnetometer(), A->GetBRMS(), A->GetAttFlag()); } Event->SetAnalysisProgress(MAssembly::c_Aspect); } else { @@ -432,9 +432,9 @@ bool MModuleLoaderMeasurementsBinary::AnalyzeEvent(MReadOutAssembly* Event) Event->SetAspectIncomplete(true); } } - Event->SetAnalysisProgress(MAssembly::c_EventLoader | - MAssembly::c_EventLoaderMeasurement | - MAssembly::c_EventOrdering); + Event->SetAnalysisProgress(MAssembly::c_EventLoader | + MAssembly::c_EventLoaderMeasurement | + MAssembly::c_EventOrdering); if (Event->GetTime().GetAsSystemSeconds() == 0) { Event->SetTimeIncomplete(true); diff --git a/src/MModuleLoaderMeasurementsHDF.cxx b/src/MModuleLoaderMeasurementsHDF.cxx new file mode 100644 index 00000000..24f1f918 --- /dev/null +++ b/src/MModuleLoaderMeasurementsHDF.cxx @@ -0,0 +1,568 @@ +/* + * MModuleLoaderMeasurementsHDF.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleLoaderMeasurementsHDF +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleLoaderMeasurementsHDF.h" + +// Standard libs: +#include +using namespace std; + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MGUIOptionsLoaderMeasurementsHDF.h" +#include "MReadOut.h" +#include "MReadOutSequence.h" +#include "MReadOutElementDoubleStrip.h" +#include "MReadOutDataADCValue.h" +#include "MReadOutDataTiming.h" +#include "MReadOutDataOrigins.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleLoaderMeasurementsHDF) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsHDF::MModuleLoaderMeasurementsHDF() : MModuleLoaderMeasurements() +{ + // Construct an instance of MModuleLoaderMeasurementsHDF + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "Measurement loader for HDF files"; + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "XmlTagMeasurementLoaderHDF"; + + // This is a special start module which can generate its own events + m_IsStartModule = true; + + // Allow the use of multiple threads and instances + m_AllowMultiThreading = true; + m_AllowMultipleInstances = false; + + m_LoadContinuationFiles = false; + m_FileNameStripMap = ""; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleLoaderMeasurementsHDF::~MModuleLoaderMeasurementsHDF() +{ + // Delete this instance of MModuleLoaderMeasurementsHDF +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleLoaderMeasurementsHDF::Initialize() +{ + // Initialize the module + + // Clean: + m_FileType = "Unknown"; + m_Detector = "Unknown"; + m_Version = -1; + /* + m_StartObservationTime = MTime(0); + m_EndObservationTime = MTime(0); + m_StartClock = numeric_limits::max(); + m_EndClock = numeric_limits::max(); + */ + + m_TotalHits = 0; + m_CurrentHit = 0; + + m_NumberOfEventIDRollOvers = 0; + m_LastEventID = 0; + + if (MFile::Exists(m_FileName) == false) { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout< 0) { + DataSet VersionDataset = m_HDFFile.openDataSet("/HDFVersion"); + + // Create compound type for version + CompType VersionType(sizeof(MHDFStripHitVersionString)); + StrType StringType(PredType::C_S1, 256); + VersionType.insertMember("string_col", HOFFSET(MHDFStripHitVersionString, string_col), StringType); + + MHDFStripHitVersionString VS; + VersionDataset.read(&VS, VersionType); + + if (string(VS.string_col) == "1.2") { + m_HDFStripHitVersion = MHDFStripHitVersion::V1_2; + } else { + if (g_Verbosity >= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<= c_Error) cout<= m_TotalHits) { + if (m_LoadContinuationFiles == true) { + // Check if we have more files to load: + MString FileName = m_FileName; + MString NextSuffix = MString("_") + (m_ContinuationFileID+1) + ".hdf5"; + FileName.ReplaceAllInPlace(".hdf5", NextSuffix); + if (MFile::Exists(FileName) == true) { + if (OpenHDF5File(FileName) == false) { + cout<= m_CurrentBatchSize) { + if (ReadBatchHits() == false) { + return false; + } + } + + // Extract the data we need + uint16_t EventID; + uint64_t TimeCode; + uint16_t StripID; + uint16_t ADCs; + uint16_t TACs; + uint8_t NumberOfHits; + uint8_t HitType; + uint8_t TimingType; + + if (m_HDFStripHitVersion == MHDFStripHitVersion::V1_0) { + MHDFStripHit_V1_0& Hit = m_Buffer_1_0[m_CurrentBatchIndex]; + ++m_CurrentBatchIndex; + ++m_CurrentHit; + + EventID = Hit.m_EventID; + TimeCode = Hit.m_TimeCode; + StripID = Hit.m_StripID; + ADCs = Hit.m_EnergyData; + TACs = Hit.m_TimingData; + NumberOfHits = Hit.m_Hits; + HitType = Hit.m_HitType; + TimingType = Hit.m_TimingType; + + } else if (m_HDFStripHitVersion == MHDFStripHitVersion::V1_2) { + MHDFStripHit_V1_2& Hit = m_Buffer_1_2[m_CurrentBatchIndex]; + ++m_CurrentBatchIndex; + ++m_CurrentHit; + + EventID = Hit.m_EventID; + TimeCode = Hit.m_TimeCode; + StripID = Hit.m_StripID; + ADCs = Hit.m_EnergyData; + TACs = Hit.m_TimingData; + NumberOfHits = Hit.m_Hits; + HitType = Hit.m_HitType; + TimingType = Hit.m_TimingType; + + } else { + if (g_Verbosity >= c_Error) cout<= c_Info) { + cout<= c_Error) cout<GetNStripHits() > 0) { + MStripHit* H = Event->GetStripHit(0); + delete H; + Event->RemoveStripHit(0); + } + } + IsZeroDataBug = false; + } + + if (EventID < m_LastEventID) { + m_NumberOfEventIDRollOvers++; + } + m_LastEventID = EventID; + + unsigned long LongEventID = EventID + m_NumberOfEventIDRollOvers*(numeric_limits::max() + 1); + + Event->SetID(LongEventID); + if (m_HDFStripHitVersion == MHDFStripHitVersion::V1_0) { + Event->SetCL(TimeCode); + } else { + Event->SetTI(TimeCode); + } + + if (m_StripMap.HasReadOutID(StripID) == true) { + MStripHit* H = new MStripHit(); + H->SetDetectorID(m_StripMap.GetDetectorID(StripID)); + H->SetStripID(m_StripMap.GetStripNumber(StripID)); + H->IsLowVoltageStrip(m_StripMap.IsLowVoltage(StripID)); + H->SetADCUnits(ADCs); + H->SetTAC(TACs); + + + // Set boolean flags based on HitType and TimingType + H->IsGuardRing(HitType == 2); + if (H->IsGuardRing() == true) { + Event->SetGuardRingVeto(true); + } + H->IsNearestNeighbor(HitType == 1); + H->HasFastTiming(TimingType == 1); + + Event->AddStripHit(H); + + } else { + if (g_Verbosity >= c_Error) cout<(NumberOfHits); + StripHitIndex++; + } + + Event->SetAnalysisProgress(MAssembly::c_EventLoader | MAssembly::c_EventLoaderMeasurement); + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleLoaderMeasurementsHDF::Finalize() +{ + // Initialize the module + + MModule::Finalize(); + + cout<<"MModuleLoaderMeasurementsHDF: "<GetNode("FileNameHDF5"); + if (FileNameHDF5Node != nullptr) { + m_FileName = FileNameHDF5Node->GetValue(); + } + + MXmlNode* LoadContinuationFilesNode = Node->GetNode("LoadContinuationFiles"); + if (LoadContinuationFilesNode != nullptr) { + m_LoadContinuationFiles = LoadContinuationFilesNode->GetValueAsBoolean(); + } + + MXmlNode* FileNameStripMapNode = Node->GetNode("FileNameStripMap"); + if (FileNameStripMapNode != nullptr) { + m_FileNameStripMap = FileNameStripMapNode->GetValue(); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleLoaderMeasurementsHDF::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + new MXmlNode(Node, "FileNameHDF5", m_FileName); + new MXmlNode(Node, "LoadContinuationFiles", m_LoadContinuationFiles); + new MXmlNode(Node, "FileNameStripMap", m_FileNameStripMap); + + return Node; +} + + +/////////////////////////////////////////////////////////////////////////////// + + +void MModuleLoaderMeasurementsHDF::ShowOptionsGUI() +{ + //! Show the options GUI + + MGUIOptionsLoaderMeasurementsHDF* Options = new MGUIOptionsLoaderMeasurementsHDF(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +// MModuleLoaderMeasurementsHDF.cxx: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleLoaderMeasurementsROA.cxx b/src/MModuleLoaderMeasurementsROA.cxx index ba7da260..5d17749a 100644 --- a/src/MModuleLoaderMeasurementsROA.cxx +++ b/src/MModuleLoaderMeasurementsROA.cxx @@ -33,6 +33,7 @@ // MEGAlib libs: #include "MGUIOptionsTemplate.h" +#include "MGUIOptionsLoaderMeasurements.h" #include "MReadOut.h" #include "MReadOutSequence.h" #include "MReadOutElementDoubleStrip.h" @@ -40,6 +41,7 @@ #include "MReadOutDataTiming.h" #include "MReadOutDataOrigins.h" + //////////////////////////////////////////////////////////////////////////////// @@ -195,7 +197,7 @@ bool MModuleLoaderMeasurementsROA::ReadNextEvent(MReadOutAssembly* Event) SH->SetStripID(Strip->GetStripID()); if (Timing != nullptr) { - SH->SetTiming(Timing->GetTiming()); + SH->SetTAC(Timing->GetTiming()); } else { cout<Create(); + gClient->WaitForUnmap(Options); +} + + // MModuleLoaderMeasurementsROA.cxx: the end... //////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleNearestNeighbor.cxx b/src/MModuleNearestNeighbor.cxx new file mode 100644 index 00000000..6d79b7ee --- /dev/null +++ b/src/MModuleNearestNeighbor.cxx @@ -0,0 +1,189 @@ +/* + * MModuleNearestNeighbor.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleNearestNeighbor +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleNearestNeighbor.h" + +// Standard libs: + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MModule.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleNearestNeighbor) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleNearestNeighbor::MModuleNearestNeighbor() : MModule() +{ + // Construct an instance of MModuleNearestNeighbor + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "NearestNeighbor"; + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "XmlTagNearestNeighbor"; + + // Set all modules, which have to be done before this module + AddPreceedingModuleType(MAssembly::c_EventLoader); + + // Set all types this modules handles + AddModuleType(MAssembly::c_NearestNeighbor); + + // Set all modules, which can follow this module + AddSucceedingModuleType(MAssembly::c_NoRestriction); + + // Set if this module has an options GUI + // Overwrite ShowOptionsGUI() with the call to the GUI! + m_HasOptionsGUI = false; + // If true, you have to derive a class from MGUIOptions (use MGUIOptionsNearestNeighbor) + // and implement all your GUI options + + // Can the program be run multi-threaded + m_AllowMultiThreading = true; + + // Can we use multiple instances of this class + m_AllowMultipleInstances = true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleNearestNeighbor::~MModuleNearestNeighbor() +{ + // Delete this instance of MModuleNearestNeighbor +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleNearestNeighbor::Initialize() +{ + // Initialize the module + + return MModule::Initialize(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleNearestNeighbor::AnalyzeEvent(MReadOutAssembly* Event) +{ + // For the time being remove all next neighbor hits + + for (unsigned int i = 0; i < Event->GetNStripHits();) { + MStripHit* SH = Event->GetStripHit(i); + if (SH->IsNearestNeighbor() == true) { + Event->RemoveStripHit(i); + delete SH; + } else { + ++i; + } + } + + Event->SetAnalysisProgress(MAssembly::c_NearestNeighbor); + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleNearestNeighbor::Finalize() +{ + // Finalize the analysis - do all cleanup, i.e., undo Initialize() + + MModule::Finalize(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleNearestNeighbor::ShowOptionsGUI() +{ + //! Show the options GUI --- has to be overwritten! + + /* + MGUIOptionsNearestNeighbor* Options = new MGUIOptionsNearestNeighbor(this); + Options->Create(); + gClient->WaitForUnmap(Options); + */ +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleNearestNeighbor::ReadXmlConfiguration(MXmlNode* Node) +{ + //! Read the configuration data from an XML node + + /* + MXmlNode* SomeTagNode = Node->GetNode("SomeTag"); + if (SomeTagNode != 0) { + m_SomeTagValue = SomeTagNode->GetValue(); + } + */ + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleNearestNeighbor::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + + /* + MXmlNode* SomeTagNode = new MXmlNode(Node, "SomeTag", "SomeValue"); + */ + + return Node; +} + + +// MModuleNearestNeighbor.cxx: the end... +//////////////////////////////////////////////////////////////////////////////// diff --git a/src/MModuleStripPairingChiSquare.cxx b/src/MModuleStripPairingChiSquare.cxx index fe60b8d4..28c846d2 100644 --- a/src/MModuleStripPairingChiSquare.cxx +++ b/src/MModuleStripPairingChiSquare.cxx @@ -62,6 +62,7 @@ MModuleStripPairingChiSquare::MModuleStripPairingChiSquare() : MModule() // Set all modules, which have to be done before this module AddPreceedingModuleType(MAssembly::c_EventLoader); AddPreceedingModuleType(MAssembly::c_EnergyCalibration); + AddPreceedingModuleType(MAssembly::c_TACcut); // Set all types this modules handles AddModuleType(MAssembly::c_StripPairing); diff --git a/src/MModuleStripPairingGreedy.cxx b/src/MModuleStripPairingGreedy.cxx index 124aa061..a697a947 100644 --- a/src/MModuleStripPairingGreedy.cxx +++ b/src/MModuleStripPairingGreedy.cxx @@ -61,6 +61,7 @@ MModuleStripPairingGreedy::MModuleStripPairingGreedy() : MModule() // Set all modules, which have to be done before this module AddPreceedingModuleType(MAssembly::c_EventLoader); AddPreceedingModuleType(MAssembly::c_EnergyCalibration); + AddPreceedingModuleType(MAssembly::c_TACcut); //AddPreceedingModuleType(MAssembly::c_CrosstalkCorrection); // Set all types this modules handles diff --git a/src/MModuleTACcut.cxx b/src/MModuleTACcut.cxx new file mode 100644 index 00000000..07b2de7f --- /dev/null +++ b/src/MModuleTACcut.cxx @@ -0,0 +1,459 @@ +/* + * MModuleTACcut.cxx + * + * + * Copyright (C) by Andreas Zoglauer. + * All rights reserved. + * + * + * This code implementation is the intellectual property of + * Andreas Zoglauer. + * + * By copying, distributing or modifying the Program (or any work + * based on the Program) you indicate your acceptance of this statement, + * and all its terms. + * + */ + + +//////////////////////////////////////////////////////////////////////////////// +// +// MModuleTACcut +// +//////////////////////////////////////////////////////////////////////////////// + + +// Include the header: +#include "MModuleTACcut.h" + +// Standard libs: + +// ROOT libs: +#include "TGClient.h" + +// MEGAlib libs: +#include "MModule.h" +#include "MGUIOptionsTACcut.h" +#include "MGUIExpoTACcut.h" + + +//////////////////////////////////////////////////////////////////////////////// + + +#ifdef ___CLING___ +ClassImp(MModuleTACcut) +#endif + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleTACcut::MModuleTACcut() : MModule() +{ + // Construct an instance of MModuleTACcut + + // Set all module relevant information + + // Set the module name --- has to be unique + m_Name = "TAC Calibration"; + + // Set the XML tag --- has to be unique --- no spaces allowed + m_XmlTag = "XmlTagTACcut"; + + // Set all modules, which have to be done before this module + AddPreceedingModuleType(MAssembly::c_EventLoader); + //AddPreceedingModuleType(MAssembly::c_EnergyCalibration); + + // Set all types this modules handles + AddModuleType(MAssembly::c_TACcut); + + // Set all modules, which can follow this module + AddSucceedingModuleType(MAssembly::c_StripPairing); + AddSucceedingModuleType(MAssembly::c_DepthCorrection); + + // Set if this module has an options GUI + // Overwrite ShowOptionsGUI() with the call to the GUI! + m_HasOptionsGUI = true; + // If true, you have to derive a class from MGUIOptions (use MGUIOptionsTACcut) + // and implement all your GUI options + + // Can the program be run multi-threaded + m_AllowMultiThreading = true; + + // Can we use multiple instances of this class + m_AllowMultipleInstances = true; + + m_SideToIndex = {{'l', 0}, {'h', 1}, {'0', 0}, {'1', 1}, {'p', 0}, {'n', 1}}; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MModuleTACcut::~MModuleTACcut() +{ + // Delete this instance of MModuleTACcut +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTACcut::Initialize() +{ + // Initialize the module + + if (LoadTACCalFile(m_TACCalFile) == false) { + cout<SetTACHistogramArrangement(m_DetectorIDs); + for (unsigned int i = 0; i < m_DetectorIDs.size(); ++i) { + unsigned int DetID = m_DetectorIDs[i]; + m_ExpoTACcut->SetTACHistogramParameters(DetID, 200, 0, 6000); + } + m_Expos.push_back(m_ExpoTACcut); +} + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTACcut::AnalyzeEvent(MReadOutAssembly* Event) +{ + // Start with sanity checks: + for (unsigned int i = 0; i < Event->GetNStripHits(); ++i) { + MStripHit* SH = Event->GetStripHit(i); + + int DetID = SH->GetDetectorID(); + int StripID = SH->GetStripID(); + char Side = SH->IsLowVoltageStrip() ? 'l' : 'h'; + + if (DetID >= m_TACCal.size()) { + cout<IsGuardRing() == false) { + cout<= m_TACCal[DetID][m_SideToIndex[Side]].size()) && (SH->IsGuardRing()==false)) { + cout<= m_TACCut.size()) { + cout<IsGuardRing()==false) { + cout<= m_TACCut[DetID][m_SideToIndex[Side]].size()) && (SH->IsGuardRing()==false)) { + cout<::max(); + for (unsigned int i = 0; i < Event->GetNStripHits(); ++i) { + MStripHit* SH = Event->GetStripHit(i); + double TAC_timing = SH->GetTAC(); + double ns_timing = 0; + int DetID = SH->GetDetectorID(); + int StripID = SH->GetStripID(); + char Side = SH->IsLowVoltageStrip() ? 'l' : 'h'; + if ((TAC_timing == 0) || (SH->IsGuardRing() == true)) { + SH->HasCalibratedTiming(false); + } else { + ns_timing = TAC_timing*m_TACCal[DetID][m_SideToIndex[Side]][StripID][0] + m_TACCal[DetID][m_SideToIndex[Side]][StripID][1]; + SH->HasCalibratedTiming(true); + } + SH->SetTiming(ns_timing); + if ((SH->HasCalibratedTiming()==true) && (ns_timing > MaxTAC) && (SH->HasFastTiming()==true) && (SH->IsNearestNeighbor()==false)) { + MaxTAC = ns_timing; + } + } + + for (unsigned int i = 0; i < Event->GetNStripHits();) { + MStripHit* SH = Event->GetStripHit(i); + bool Passed = true; + if ((SH->HasCalibratedTiming() == true) && (SH->IsGuardRing()==false)) { + double SHTiming = SH->GetTiming(); + int DetID = SH->GetDetectorID(); + int StripID = SH->GetStripID(); + char Side = SH->IsLowVoltageStrip() ? 'l' : 'h'; + double FLNoiseCut = m_TACCut[DetID][m_SideToIndex[Side]][StripID][4]; + if ((SHTiming < FLNoiseCut)) { + Passed = false; + } else if (SH->HasFastTiming()==true) { + double ShapingOffset = m_TACCut[DetID][m_SideToIndex[Side]][StripID][0]; + double CoincidenceWindow = m_TACCut[DetID][m_SideToIndex[Side]][StripID][1]; + double DisableTime = m_TACCut[DetID][m_SideToIndex[Side]][StripID][2]; + double FlagToEnDelay = m_TACCut[DetID][m_SideToIndex[Side]][StripID][3]; + double FlagDelay = m_TACCut[DetID][m_SideToIndex[Side]][StripID][5]; + double TotalOffset = ShapingOffset + DisableTime + FlagToEnDelay + FlagDelay; + if ((SHTiming > TotalOffset + CoincidenceWindow) || (SHTiming < TotalOffset) || (SHTiming < MaxTAC - CoincidenceWindow)) { + Passed = false; + } else if (HasExpos()==true) { + m_ExpoTACcut->AddTAC(DetID, SHTiming); + } + } + } + if (Passed == true) { + ++i; + } else { + Event->RemoveStripHit(i); + delete SH; + } + } + + Event->SetAnalysisProgress(MAssembly::c_TACcut); + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleTACcut::Finalize() +{ + MModule::Finalize(); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +void MModuleTACcut::ShowOptionsGUI() +{ + //! Show the options GUI --- has to be overwritten! + + MGUIOptionsTACcut* Options = new MGUIOptionsTACcut(this); + Options->Create(); + gClient->WaitForUnmap(Options); +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTACcut::ReadXmlConfiguration(MXmlNode* Node) +{ + //! Read the configuration data from an XML node + + MXmlNode* TACCalFileNameNode = Node->GetNode("TACCalFileName"); + if (TACCalFileNameNode != nullptr) { + SetTACCalFileName(TACCalFileNameNode->GetValue()); + } + + MXmlNode* TACCutFileNameNode = Node->GetNode("TACCutFileName"); + if (TACCutFileNameNode != nullptr) { + SetTACCutFileName(TACCutFileNameNode->GetValue()); + } + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +MXmlNode* MModuleTACcut::CreateXmlConfiguration() +{ + //! Create an XML node tree from the configuration + + MXmlNode* Node = new MXmlNode(0, m_XmlTag); + + new MXmlNode(Node, "TACCalFileName", m_TACCalFile); + new MXmlNode(Node, "TACCutFileName", m_TACCutFile); + + return Node; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +bool MModuleTACcut::LoadTACCalFile(MString FName) +{ + // Read in the TAC Calibration file, which should contain for each strip: + // DetID, Side (h or l for high or low voltage), TAC cal, TAC cal error, TAC cal offset, TAC offset error + // OR: + // ReadOutID, Detector, Side, Strip, TAC cal, TAC cal error, TAC offset, TAC offset error + MFile F; + if (F.Open(FName) == false) { + cout< Tokens = Line.Tokenize(","); + if ((Tokens.size() == 7) || (Tokens.size() == 8)) { + int IndexOffset = Tokens.size() % 7; + int DetID = Tokens[0+IndexOffset].ToInt(); + MString SideString = Tokens[1+IndexOffset].Trim(); + char Side; + if (SideString.Length()!=1) { + cout< CalValues; + CalValues.push_back(TACCal); CalValues.push_back(Offset); CalValues.push_back(TACCalError); CalValues.push_back(OffsetError); + + if (m_TACCal.find(DetID) == m_TACCal.end()) { + vector>> TempVector; + unordered_map> TempMapLV; + unordered_map> TempMapHV; + m_TACCal[DetID] = TempVector; + m_TACCal[DetID].push_back(TempMapLV); + m_TACCal[DetID].push_back(TempMapHV); + } + + if (find(m_DetectorIDs.begin(), m_DetectorIDs.end(), DetID) == m_DetectorIDs.end()) { + m_DetectorIDs.push_back(DetID); + } + + if (m_SideToIndex.find(Side) != m_SideToIndex.end()) { + m_TACCal[DetID][m_SideToIndex[Side]][StripID] = CalValues; + } else { + cout< Tokens = Line.Tokenize(","); + if (Tokens.size() >= 5) { + int DetID = Tokens[0].ToInt(); + MString SideString = Tokens[1]; + char Side; + if (SideString.Length()!=1) { + cout< CutParams; + CutParams.push_back(ShapingOffset); CutParams.push_back(CoincidenceWindow); CutParams.push_back(DisableTime); CutParams.push_back(FlagToEnDelay); CutParams.push_back(FLNoiseCut); CutParams.push_back(FlagDelay); + + if (m_TACCut.find(DetID) == m_TACCut.end()) { + vector>> TempVector; + unordered_map> TempMapLV; + unordered_map> TempMapHV; + m_TACCut[DetID] = TempVector; + m_TACCut[DetID].push_back(TempMapLV); + m_TACCut[DetID].push_back(TempMapHV); + } + + if (find(m_DetectorIDs.begin(), m_DetectorIDs.end(), DetID) == m_DetectorIDs.end()) { + m_DetectorIDs.push_back(DetID); + } + + if (m_SideToIndex.find(Side) != m_SideToIndex.end()) { + m_TACCut[DetID][m_SideToIndex[Side]][StripID] = CutParams; + } else { + cout<StreamEvta(S); } @@ -640,8 +657,16 @@ void MReadOutAssembly::StreamRoa(ostream& S, bool) S<StreamRoa(S); + if (WithNearestNeighbors == false && m_StripHits[h]->IsNearestNeighbor() == true) { + continue; + } + m_StripHits[h]->StreamRoa(S, WithADCs, WithTACs, WithEnergies, WithTimings, WithTemperatures, WithFlags); + ++Counter; + } + if (Counter == 0) { + S<<"BD No strip hits"< tokens = Line.Tokenize(" "); - if( tokens.size() >= 10 ){ - SetDetectorID(tokens.at(1).ToInt()); - tokens.at(2) == "p" ? IsLowVoltageStrip(true) : IsLowVoltageStrip(false); - SetStripID(tokens.at(3).ToInt()); - tokens.at(4) == "0" ? HasTriggered(false) : HasTriggered(true); - SetTiming( tokens.at(5).ToDouble() ); - SetUncorrectedADCUnits( tokens.at(6).ToDouble() ); - SetADCUnits( tokens.at(7).ToDouble() ); - SetEnergy( tokens.at(8).ToDouble() ); - SetEnergyResolution( tokens.at(9).ToDouble() ); - int det_id, pos_strip, strip_id - return true; + const char* line = Line.Data(); + if (line[0] == 'S' && line[1] == 'H') { + int det_id, strip_id, has_triggered, timing, un_adc, adc; + float energy, energy_res; + char pos_strip; + unsigned int flags; + sscanf(&line[3],"%d %c %d %d %d %d %d %f %f %u", &det_id, &pos_strip, &strip_id, &has_triggered, &timing, &un_adc, &adc, &energy, &energy_res, &flags); + SetDetectorID(det_id); + pos_strip == 'l' ? IsLowVoltageStrip(true) : IsLowVoltageStrip(false); + SetStripID(strip_id); + has_triggered == 0 ? HasTriggered(false) : HasTriggered(true); + SetTiming((double)timing); + SetUncorrectedADCUnits((double)un_adc); + SetADCUnits((double)adc); + SetEnergy(energy); + SetEnergyResolution(energy_res); + ParseFlags(flags); + return true; } else { - return false; + return false; } - */ } @@ -170,7 +154,8 @@ bool MStripHit::StreamDat(ostream& S, int Version) <GetDetectorID()<<" " <GetStripID()<<" " - <<((m_ReadOutElement->IsLowVoltageStrip() == true) ? "l" : "h")<<" " - <IsLowVoltageStrip() == true) ? "l" : "h")<<" "; + if (WithADC == true) { + S<= c_Error) cout<<"MStripMap: Unable to load file: "<GetNTokens() == 0) continue; + if (Parser.GetTokenizerAt(i)->GetTokenAtAsString(0).BeginsWith("#") == true) continue; + if (Parser.GetTokenizerAt(i)->GetNTokens() == 9) { + MSingleStripMapping SM; + SM.m_ReadOutID = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(0); + SM.m_RTB = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(1); + SM.m_DRM = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(2); + SM.m_IsPrimary = Parser.GetTokenizerAt(i)->GetTokenAtAsBoolean(3); + SM.m_ASICID = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(4); + SM.m_ChannelID = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(5); + SM.m_DetectorID = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(6); + SM.m_IsLowVoltage = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(7) == 0 ? true : false; + SM.m_StripNumber = Parser.GetTokenizerAt(i)->GetTokenAtAsUnsignedInt(8); + m_StripMappings.push_back(SM); + } + } + + // Sort by m_ReadOutID: + sort(m_StripMappings.begin(), m_StripMappings.end(), [](const MSingleStripMapping& A, const MSingleStripMapping& B) { return A.m_ReadOutID < B.m_ReadOutID; }); + + return true; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Return true if the given read-out ID is on file +bool MStripMap::HasReadOutID(unsigned int ROI) const +{ + auto Iter = lower_bound(m_StripMappings.begin(), m_StripMappings.end(), ROI, [](const MSingleStripMapping& SSM, unsigned int ID) { return SSM.m_ReadOutID < ID; }); + return Iter != m_StripMappings.end() && Iter->m_ReadOutID == ROI; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Return the index of the ROI, throw an exception otherwise +unsigned int MStripMap::GetReadOutIDIndex(unsigned int ROI) const +{ + auto Iter = lower_bound(m_StripMappings.begin(), m_StripMappings.end(), ROI, [](const MSingleStripMapping& SSM, unsigned int ID) { return SSM.m_ReadOutID < ID; }); + + if (Iter != m_StripMappings.end() && Iter->m_ReadOutID == ROI) { + return distance(m_StripMappings.begin(), Iter); + } else { + throw MExceptionValueNotFound(ROI, "vector of read-out IDs"); + } +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Get detector by read-out ID +unsigned int MStripMap::GetDetectorID(unsigned int ROI) const +{ + return m_StripMappings[GetReadOutIDIndex(ROI)].m_DetectorID; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Get detector side by read-out ID +bool MStripMap::IsLowVoltage(unsigned int ROI) const +{ + return m_StripMappings[GetReadOutIDIndex(ROI)].m_IsLowVoltage; +} + + +//////////////////////////////////////////////////////////////////////////////// + + +//! Get strip ID by read-out ID +unsigned int MStripMap::GetStripNumber(unsigned int ROI) const +{ + return m_StripMappings[GetReadOutIDIndex(ROI)].m_StripNumber; +} + + +// MStripMap.cxx: the end... +////////////////////////////////////////////////////////////////////////////////