diff --git a/CMakeLists.txt b/CMakeLists.txt index 11df14a6b8..9dad2b30b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -123,6 +123,7 @@ add_subdirectory(Blinding) add_subdirectory(bin) add_subdirectory(boost_fix) add_subdirectory(BTrkLegacy) +add_subdirectory(CaloCalibration) add_subdirectory(CaloCluster) add_subdirectory(CaloConditions) add_subdirectory(CaloConfig) diff --git a/CaloCalibration/CMakeLists.txt b/CaloCalibration/CMakeLists.txt new file mode 100644 index 0000000000..6b5da16551 --- /dev/null +++ b/CaloCalibration/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(CosmicsCalib) +add_subdirectory(NoiseAnalysis) diff --git a/CaloCalibration/CosmicsCalib/CMakeLists.txt b/CaloCalibration/CosmicsCalib/CMakeLists.txt new file mode 100755 index 0000000000..6eecf48a87 --- /dev/null +++ b/CaloCalibration/CosmicsCalib/CMakeLists.txt @@ -0,0 +1,40 @@ +cet_build_plugin(CaloCosmicEnecalib art::module + REG_SOURCE src/CaloCosmicEnecalib_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::CalorimeterGeom + Offline::DataProducts + Offline::GeometryService + Offline::RecoDataProducts + ROOT::Hist + ROOT::Tree +) + +cet_build_plugin(CaloCosmicEnergy art::module + REG_SOURCE src/CaloCosmicEnergy_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::CalorimeterGeom + Offline::DataProducts + Offline::GeometryService + Offline::RecoDataProducts + ROOT::Hist + ROOT::Tree +) + +cet_build_plugin(caloT0alig art::module + REG_SOURCE src/caloT0alig_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::CalorimeterGeom + Offline::DataProducts + Offline::GeometryService + Offline::RecoDataProducts + ROOT::Hist + ROOT::Tree +) + + +#install_headers(SUBDIRS inc) +install_source(SUBDIRS src) +install_fhicl(SUBDIRS fcl SUBDIRNAME Offline/CaloCalibration/CosmicsCalib/fcl) diff --git a/CaloCalibration/CosmicsCalib/fcl/prolog.fcl b/CaloCalibration/CosmicsCalib/fcl/prolog.fcl new file mode 100644 index 0000000000..65ea0a3f75 --- /dev/null +++ b/CaloCalibration/CosmicsCalib/fcl/prolog.fcl @@ -0,0 +1,41 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/CaloMC/fcl/prolog.fcl" +#include "Offline/CaloReco/fcl/prolog.fcl" +BEGIN_PROLOG + +CaloCosmicEnecalib : { + module_type : CaloCosmicEnecalib + diagLevel : 0 + useMeV : false + CutNCryHit : 3 + CutEnergyDep : 250. + LYSOcut : 500. + CutChi2Norm : 2.5 + CaloClusterTag : CaloClusterMaker + CaloHitTag : CaloHitMaker + OutCalibFile : "MIPfile.dat" +} + +CaloCosmicEnergy : { + module_type : CaloCosmicEnergy + diagLevel : 0 + CutNCryHit : 3 + CutEnergyDep : 10. + CutChi2Norm : 2.5 + CaloClusterTag : CaloClusterMaker + CaloHitTag : CaloHitMaker + OutCalibFile : "calib_parameters.dat" + OutCalibFileAsym : "npe_noise.dat" +} + + +CosmicsCalib: { + analyzers: { + CaloCosmicEnecalib : @local::CaloCosmicEnecalib + } +} + + +END_PROLOG diff --git a/CaloCalibration/CosmicsCalib/src/CaloCosmicEnecalib_module.cc b/CaloCalibration/CosmicsCalib/src/CaloCosmicEnecalib_module.cc new file mode 100644 index 0000000000..c244b1ac38 --- /dev/null +++ b/CaloCalibration/CosmicsCalib/src/CaloCosmicEnecalib_module.cc @@ -0,0 +1,769 @@ +// ============================================================================= +// extrapolate the MIP MPV for offline calorimeter calibration with cosmics +// 09-Mar-2026 - S.Salamino & S.Giovannella +// +// Starting code for track selection and path normalization: +// CaloCosmicEnergy module from P.Fedeli & S.Giovannella +// +// ============================================================================= + +// includes + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Provenance.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/Selector.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileDirectory.h" +#include "art_root_io/TFileService.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Sequence.h" + +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CalorimeterGeom/inc/CrystalMapper.hh" +#include "Offline/DataProducts/inc/CaloConst.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" + +#include "messagefacility/MessageLogger/MessageLogger.h" + +// C++ includes +#include +#include +#include +#include +#include + +// ROOT includes +#include "TF1.h" +#include "TFitResult.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TH1F.h" +#include "TH2F.h" + +namespace mu2e { + +class CaloCosmicEnecalib : public art::EDAnalyzer { + +public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom diagLevel{Name("diagLevel"), Comment("Diagnostic level"), 0}; + fhicl::Atom useMeV{Name("useMeV"), Comment("Set to true for data in MeV, false for ADC"), + false}; + fhicl::Atom CutNCryHit{Name("CutNCryHit"), + Comment("Minimum number of crystals in the event"), 3}; + fhicl::Atom CutEnergyDep{Name("CutEnergyDep"), Comment("Minimum energy of the hit"), + 250.}; // ADC + fhicl::Atom LYSOcut{Name("LYSOcut"), Comment("Minimum energy of the hit for LYSO"), + 500.}; // ADC + fhicl::Atom CutChi2Norm{Name("CutChi2Norm"), + Comment("Maximum chi2 threshold for linear fits"), 2.5}; + fhicl::Atom CaloClusterTag{ + Name("CaloClusterTag"), Comment("Tag for Calorimeter cluster collection"), art::InputTag()}; + fhicl::Atom CaloHitTag{ + Name("CaloHitTag"), Comment("Tag for Calorimeter hit collection"), art::InputTag()}; + fhicl::Atom OutCalibFile{Name("OutCalibFile"), + Comment("Name for output .dat MIP calibration file"), + "calib_parameters.dat"}; + }; + + explicit CaloCosmicEnecalib(const art::EDAnalyzer::Table& config); + virtual ~CaloCosmicEnecalib() {}; + void beginJob() override; + void beginRun(art::Run const& run) override; + void analyze(art::Event const& event) override; + void endJob() override; + +private: + art::ServiceHandle tfs; + + art::TFileDirectory tfdir = tfs->mkdir("All_tracks"); + art::TFileDirectory tfdirv = tfs->mkdir("Vertical_tracks"); + art::TFileDirectory tfdird = tfs->mkdir("Diagonal_tracks"); + art::TFileDirectory tfdirg = tfs->mkdir("General_tracks"); + art::TFileDirectory tfdir_res = tfs->mkdir("MPV_distributions"); + + // calo parameters + static constexpr int nSiPMs = CaloConst::_nSiPMPerCrystal; + static constexpr int nCrystals = CaloConst::_nCrystalPerDisk; + static constexpr int nDisks = CaloConst::_nDisk; + static constexpr int nROchan = nDisks * nCrystals * nSiPMs; + const Calorimeter* cal; // For calorimeter geometry + float cryDim; // Total crystal dimension + float MaxDxVertical; // Max Dx value for vertical tracks + + // fcl parameters + int _diagLevel; + bool _useMeV; + int CutNCryHit; // > + float CutEnergyDep; // > + float LYSOcut; // > + float CutChi2Norm; // < + art::ProductToken _caloClusterToken; + art::ProductToken _caloHitToken; + std::ofstream _outputfile; + + // histogram variables + TH1F* hSiPM[nROchan] = {nullptr}; + TH1F* hSiPMv[nROchan] = {nullptr}; + TH1F* hSiPMd[nROchan] = {nullptr}; + TH1F* hSiPMg[nROchan] = {nullptr}; + + TH1F* MPV = nullptr; + TH1F* MPVv = nullptr; + TH1F* MPVd = nullptr; + TH1F* MPVg = nullptr; + TH1F* hWidth = nullptr; + TH1F* hSigma = nullptr; + + int hSiPMbins = 330; // bins for hSiPM histograms + int lowStatThreshold = 370; // Number of events threshold to rebin histos with low statistics + int midStatThreshold = 600; // Number of events threshold to rebin histos with medium statistics + double hSiPMxmax; + double hMPVxmax; + double hWxmax; + double hSxmax; + + TF1* myPoly = nullptr; + TF1* mylang = nullptr; + TF1* mygaus = nullptr; + + // helper variables for Langaus fit + struct LangausResult { // langaus fit parameters + float mpv = -99.; + float mpvErr = -99.; + float width = -99.; + float widthErr = -99.; + float sigma = -99.; + float sigmaErr = -99.; + float chi2 = -99.; + int ndf = 0; + int nev = 0; + void reset() { *this = LangausResult(); } + }; + + struct FitFlags { // flags to determine fit quality + bool badWidthLow = false; + bool badWidthHigh = false; + bool badSigmaLow = false; + bool badSigmaHigh = false; + bool isBad() const { return badWidthLow || badWidthHigh || badSigmaLow || badSigmaHigh; } + void fitQuality(float width, float sigma, TH1F* hW, TH1F* hS) { + if (!hW || !hS || width <= 0 || sigma <= 0) + return; + badWidthLow = (width < hW->GetMean() - 3 * hW->GetRMS()); + badWidthHigh = (width > hW->GetMean() + 4 * hW->GetRMS()); + badSigmaLow = (sigma < hS->GetMean() - 3 * hS->GetRMS()); + badSigmaHigh = (sigma > hS->GetMean() + 4 * hS->GetRMS()); + } + }; + + std::vector sipmParams; + std::vector sipmParamsv; + std::vector sipmParamsd; + std::vector sipmParamsg; + + static double langaus(Double_t* x, Double_t* par); + float findpath(float m, float q, float x, float y); + LangausResult fitLangaus(TH1F* histo, int refit = 0, int iter = 0, float oldwidth = 0, + float oldsigma = 0); + +}; // end class + +CaloCosmicEnecalib::CaloCosmicEnecalib(const art::EDAnalyzer::Table& config) : + art::EDAnalyzer(config), _diagLevel(config().diagLevel()), _useMeV(config().useMeV()), + CutNCryHit(config().CutNCryHit()), CutEnergyDep(config().CutEnergyDep()), + LYSOcut(config().LYSOcut()), CutChi2Norm(config().CutChi2Norm()), + _caloClusterToken(consumes(config().CaloClusterTag())), + _caloHitToken(consumes(config().CaloHitTag())), + _outputfile(config().OutCalibFile()) // opens file.dat +{ + if (!_outputfile.is_open()) { + throw cet::exception("CaloCosmicEnecal") + << "ERROR! Cannot open output file " << config().OutCalibFile() << std::endl; + } +} + +void CaloCosmicEnecalib::beginRun(art::Run const& run) { + + art::ServiceHandle geom; + if (!geom->hasElement()) { + throw cet::exception("NO-CALO-GEOM") << "CaloEnergy: Calorimeter geometry not found"; + } + + GeomHandle ch; + cal = ch.get(); + cryDim = cal->G4Info().get("crystalXYLength") + + 2. * cal->G4Info().get("wrapperThickness"); + MaxDxVertical = cryDim * 1.1; +} + +void CaloCosmicEnecalib::beginJob() { + if (_diagLevel > 0) + std::cout << "CaloCosmicEnecalib: Entering beginJob" << std::endl; + + if (!_useMeV) { + hSiPMxmax = 2000.; + hMPVxmax = 1000.; + hWxmax = 200.; + hSxmax = 400.; + } else { + hSiPMxmax = 100.; + hMPVxmax = 35.; + hWxmax = 10.; + hSxmax = 20.; + } + + for (int iRou = 0; iRou < nROchan; iRou++) { + hSiPM[iRou] = + tfdir.make(Form("Channel_%i", iRou), Form("Energy deposited in channel %i", iRou), + hSiPMbins, 0, hSiPMxmax); + hSiPMv[iRou] = tfdirv.make(Form("Channel_%i_v", iRou), + Form("Energy deposited in channel %i - vertical tracks", iRou), + 80, 0, hSiPMxmax); + hSiPMd[iRou] = tfdird.make(Form("Channel_%i_d", iRou), + Form("Energy deposited in channel %i - diagonal tracks", iRou), + 80, 0, hSiPMxmax); + hSiPMg[iRou] = tfdirg.make(Form("Channel_%i_g", iRou), + Form("Energy deposited in channel %i - general tracks", iRou), + 80, 0, hSiPMxmax); + } + + MPV = tfdir_res.make("MPV_distr", "MPV distribution", 400, 0., hMPVxmax); + MPVv = tfdir_res.make("MPV_v_distr", "MPV distribution only vertical tracks", 1000, 15., + hMPVxmax); + MPVd = tfdir_res.make("MPV_d_distr", "MPV distribution only diagonal tracks", 1000, 15., + hMPVxmax); + MPVg = tfdir_res.make("MPV_g_distr", "MPV distribution only general tracks", 1000, 15., + hMPVxmax); + hWidth = tfdir_res.make("Width_distr", "Width distribution", 150, 0., hWxmax); + hSigma = tfdir_res.make("Sigma_distr", "Sigma distribution", 300, 0., hSxmax); + + myPoly = new TF1("myPoly", "pol1", -600., 600.); // polinomial to fit the track + mylang = new TF1("mylang", langaus, 0., 2100., 4); // langaus function + mygaus = new TF1("mygaus", "gaus"); + +} // end begin job + +void CaloCosmicEnecalib::analyze(art::Event const& event) { + + auto const& caloClusters = event.getProduct(_caloClusterToken); + int nCluster = caloClusters.size(); + + // Loop over clusters in the event + for (int iClu = 0; iClu < nCluster; iClu++) { + if (_diagLevel > 0) + std::cout << "iClu: " << iClu << std::endl; + + int nhits = 0; // number of crystal hits passing the selection + float Dx = 1400.; // to calculate the spread along x + + // creating arrays to store information for hits passing the selection + std::array IDs; + std::array Vmax; + std::array path; + std::vector> PosXY; + + auto const& hitsPtrVector = caloClusters[iClu].caloHitsPtrVector(); + int diskID = caloClusters[iClu].diskID(); + + // Loop over crystals in the cluster + for (auto const& hitsPtr : hitsPtrVector) { + auto const& recoDigis = hitsPtr->recoCaloDigis(); + int cryID = hitsPtr->crystalID(); + if (cryID > nCrystals * nDisks) + continue; + + // loop over hits in the crystal + for (auto const& rdPtr : recoDigis) { + double VmaxHit = -1; + int sipmID = -1; + + if (_useMeV) { + VmaxHit = rdPtr->energyDep(); + sipmID = rdPtr->SiPMID(); + } + + else { + auto const& rawDigiPtr = rdPtr->caloDigiPtr(); + if (rawDigiPtr.isNonnull()) { + const std::vector& waveform = rawDigiPtr->waveform(); + + // evaluate baseline on first 5 samples in the waveform + double baseline = 0; + for (int i = 0; i < 5; ++i) + baseline += waveform.at(i); + baseline /= 5.0; + + // find maximum of the waveform (baseline subtracted) + for (int i = 0; i < (int)waveform.size(); ++i) { + if (waveform.at(i) > VmaxHit) + VmaxHit = waveform.at(i); + } + if (_diagLevel > 0) { + std::cout << "sipmID: " << rawDigiPtr->SiPMID() << " waveform peak: " << VmaxHit + << " baseline: " << baseline; + } + VmaxHit = VmaxHit - baseline; + sipmID = rawDigiPtr->SiPMID(); + } // caloDigi check + } // use MeV or ADC + + if (_diagLevel > 0) + std::cout << " Vmax: " << VmaxHit << std::endl; + + // cut on energy deposition + float enecut = 9999.; + if (cryID == 610 || cryID == 637 || cryID == 609 || cryID == 582) + enecut = LYSOcut; + else + enecut = CutEnergyDep; + if (VmaxHit > enecut && sipmID >= 0 && sipmID < nROchan) { + float PosX = cal->mu2eToDiskFF(diskID, cal->crystal(cryID).position()).getX(); + float PosY = cal->mu2eToDiskFF(diskID, cal->crystal(cryID).position()).getY(); + + PosXY.push_back({PosX, PosY}); + IDs[nhits] = sipmID; + Vmax[nhits] = VmaxHit; + nhits++; + } // energy cut + + } // loop over hits + } // loop over crystals + + if (_diagLevel > 0) + std::cout << "nhits: " << nhits << std::endl; + + if (nhits == 0) { + if (_diagLevel > 0) + std::cout << "Cluster " << iClu << " has no valid hit." << std::endl; + continue; + } + + // get rid of duplicate positions + std::vector> fitPosXY = PosXY; + std::sort(fitPosXY.begin(), fitPosXY.end()); + auto last = std::unique(fitPosXY.begin(), fitPosXY.end()); + fitPosXY.erase(last, fitPosXY.end()); + + // find min and max x coordinate + float min_x = fitPosXY.front().first; + float max_x = fitPosXY.back().first; + Dx = abs(max_x - min_x); + + if (_diagLevel > 0) { + std::cout << "dx: " << Dx << std::endl; + } + + // cut on ncrystalhit + if (int(fitPosXY.size()) > CutNCryHit) { + + // prepare arrays for track fit + std::vector fitX, fitY, fitErr; + for (auto const& p : fitPosXY) { + fitX.push_back(p.first); + fitY.push_back(p.second); + fitErr.push_back(9.81); + if (_diagLevel > 0) { + std::cout << "x: " << p.first << " y: " << p.second << std::endl; + } + } + + // geometry track + float TrkSlope = -999; + float TrkIntercept = -999; + float chi2norm = 999; + + if (Dx <= MaxDxVertical) { // track is vertical + TrkSlope = 144.; + TrkIntercept = 0.; + } else { + double m_init = (fitY.back() - fitY.front()) / (fitX.back() - fitX.front()); + double q_init = fitY.front() - m_init * fitX.front(); + + myPoly->SetParameters(q_init, m_init); + + TGraphErrors g(fitX.size(), fitX.data(), fitY.data(), fitErr.data(), fitErr.data()); + TFitResultPtr fitresult = g.Fit(myPoly, "SQRM"); + + if ((fitresult.Get() != nullptr) && (fitresult->IsValid()) && (fitresult->Ndf() > 0)) { + TrkSlope = fitresult->Parameter(1); + TrkIntercept = fitresult->Parameter(0); + chi2norm = fitresult->Chi2() / fitresult->Ndf(); + + if (_diagLevel > 0) + std::cout << "Track fit done " << std::endl; + } + } + + // renormalize diagonal or good quality tracks + if (Dx > MaxDxVertical) { + for (int i = 0; i < nhits; i++) { + if (_diagLevel > 0) { + std::cout << "using m:" << TrkSlope << " q: " << TrkIntercept + << " x: " << PosXY[i].first << " y: " << PosXY[i].second + << " Vmax: " << Vmax[i] << std::endl; + } + path[i] = findpath(TrkSlope, TrkIntercept, PosXY[i].first, PosXY[i].second); + } + } + + // fill histograms + for (int kk = 0; kk < nhits; kk++) { + + // vertical track: no normalization + if (Dx <= MaxDxVertical) { + hSiPMv[IDs[kk]]->Fill(Vmax[kk]); + hSiPM[IDs[kk]]->Fill(Vmax[kk]); + } + + // diagonal and general tracks + else if ((chi2norm < CutChi2Norm) && (path[kk] > 0)) { + + // diagonal track + if (chi2norm < 0.01) + hSiPMd[IDs[kk]]->Fill(Vmax[kk] * cryDim / path[kk]); + // general track + else + hSiPMg[IDs[kk]]->Fill(Vmax[kk] * cryDim / path[kk]); + + hSiPM[IDs[kk]]->Fill(Vmax[kk] * cryDim / path[kk]); + } + } + } // end ncrystal cut + } // end loop iClu in ncluster +} // end analyze + +void CaloCosmicEnecalib::endJob() { + std::cout << "CaloCosmicEnecalib: entering endjob" << std::endl; + + sipmParams.assign(nROchan, LangausResult()); + sipmParamsv.assign(nROchan, LangausResult()); + sipmParamsg.assign(nROchan, LangausResult()); + sipmParamsd.assign(nROchan, LangausResult()); + int nEntries; + + for (int irou = 0; irou < nROchan; irou++) { + nEntries = hSiPM[irou]->GetEntries(); + if (nEntries > 0) { + if (_diagLevel > 0) + std::cout << " >> Fitting channel " << irou << std::endl; + + // Rebin hitsograms to help with convercence + int rebinfactor = 1; + if (nEntries <= lowStatThreshold) + rebinfactor = 6; + else if (nEntries <= midStatThreshold) + rebinfactor = 5; + else + rebinfactor = 3; + hSiPM[irou] = (TH1F*)hSiPM[irou]->Rebin(rebinfactor); + + // First attempt for langaus fit + sipmParams[irou] = fitLangaus(hSiPM[irou]); + sipmParamsv[irou] = fitLangaus(hSiPMv[irou]); + sipmParamsd[irou] = fitLangaus(hSiPMd[irou]); + sipmParamsg[irou] = fitLangaus(hSiPMg[irou]); + + hWidth->Fill(sipmParams[irou].width); + hSigma->Fill(sipmParams[irou].sigma); + if (sipmParamsv[irou].mpv > 0) + MPVv->Fill(sipmParamsv[irou].mpv); + if (sipmParamsd[irou].mpv > 0) + MPVd->Fill(sipmParamsd[irou].mpv); + if (sipmParamsg[irou].mpv > 0) + MPVg->Fill(sipmParamsg[irou].mpv); + } + } + + // Distribution of Landau width and Gaussian sigma for all channels + mygaus->SetParameters(hWidth->GetEntries(), hWidth->GetMean(), hWidth->GetRMS()); + hWidth->Fit("mygaus", "Q0"); + mygaus->SetParameters(hSigma->GetEntries(), hSigma->GetMean(), hSigma->GetRMS()); + hSigma->Fit("mygaus", "Q0"); + + // Check fit quality + for (int irou = 0; irou < nROchan; irou++) { + if (hSiPM[irou]->GetEntries() > 0) { + FitFlags flags; + flags.fitQuality(sipmParams[irou].width, sipmParams[irou].sigma, hWidth, hSigma); + + // Refit up to 10 times if fit was not good + int iter = 0; + while (flags.isBad() && iter < 10) { + if (_diagLevel > 0) + std::cout << "Refitting channel " << irou << ", iteration " << iter << std::endl; + + sipmParams[irou] = + fitLangaus(hSiPM[irou], 1, iter, sipmParams[irou].width, sipmParams[irou].sigma); + flags.fitQuality(sipmParams[irou].width, sipmParams[irou].sigma, hWidth, hSigma); + iter++; + } + if (sipmParams[irou].mpv > 0) + MPV->Fill(sipmParams[irou].mpv); + } + } + + MPV->Fit("gaus", "Q"); + MPVv->Fit("gaus", "Q"); + MPVd->Fit("gaus", "Q"); + MPVg->Fit("gaus", "Q"); + + // Write output file + _outputfile << std::left << std::setw(6) << "#ROU" << std::setw(10) << "MPV" << std::setw(10) + << "MPVerr" << std::setw(10) << "Width" << std::setw(10) << "Werr" << std::setw(10) + << "Sigma" << std::setw(10) << "Serr" << std::setw(10) << "Chi2" << std::setw(8) + << "NDF" << std::setw(8) << "Nev" << std::endl; + + for (int irou = 0; irou < nROchan; irou++) { + const auto& res = sipmParams[irou]; + + _outputfile << std::left << std::setw(6) << irou << std::fixed << std::setprecision(2) + << std::setw(10) << res.mpv << std::setw(10) << res.mpvErr << std::setw(10) + << res.width << std::setw(10) << res.widthErr << std::setw(10) << res.sigma + << std::setw(10) << res.sigmaErr << std::setw(10) << res.chi2 << std::setw(8) + << res.ndf << std::setw(8) << res.nev << std::endl; + } + + _outputfile.close(); + delete mygaus; + delete myPoly; + delete mylang; + +} // end endjob + +double CaloCosmicEnecalib::langaus(double* x, double* par) { + // Fit parameters: + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + // + // In the Landau distribution (represented by the CERNLIB approximation), + // the maximum is located at x=-0.22278298 with the location parameter=0. + // This shift is corrected within this function, so that the actual + // maximum is identical to the MP parameter. + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // number of convolution steps + Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx; + Double_t mpc; + Double_t fland; + Double_t sum = 0.0; + Double_t xlow, xupp; + Double_t step; + Double_t i; + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp - xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for (i = 1.0; i <= np / 2; i++) { + xx = xlow + (i - .5) * step; + fland = TMath::Landau(xx, mpc, par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0], xx, par[3]); + + xx = xupp - (i - .5) * step; + fland = TMath::Landau(xx, mpc, par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0], xx, par[3]); + } + + return (par[2] * step * sum * invsq2pi / par[3]); +} // end langaus + +float CaloCosmicEnecalib::findpath(float m, float q, float x, float y) { + float xup = 0; + float yright = 0; + float path = 0; + float xlow = 0; + float yleft = 0; + + int diag = 0; + float halfcrysize = cryDim / 2.; + + if (m != 0) { + xup = ((y + halfcrysize) - q) / m; + xlow = ((y - halfcrysize) - q) / m; + yleft = m * (x - halfcrysize) + q; + yright = m * (x + halfcrysize) + q; + } + if (xup < (x + halfcrysize) && xup > (x - halfcrysize)) { // hit the top of crystal + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + path = halfcrysize * 2 * TMath::Sqrt(1 / (m * m) + 1); + if (diag == 1) + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << xlow << ";" + << y - halfcrysize << ")" << std::endl; + } + if (xlow < (x - halfcrysize)) { + path = TMath::Sqrt((xup - (x - halfcrysize)) * (xup - (x - halfcrysize)) + + ((y + halfcrysize) - yleft) * ((y + halfcrysize) - yleft)); + if (diag == 1) { + std::cout << "Track exit from left" << std::endl; + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << x - halfcrysize + << ";" << yleft << ")" << std::endl; + } + } + if (xlow > (x + halfcrysize)) { + path = TMath::Sqrt((xup - (x + halfcrysize)) * (xup - (x + halfcrysize)) + + ((y + halfcrysize) - yright) * ((y + halfcrysize) - yright)); + if (diag == 1) { + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << x + halfcrysize + << ";" << yright << ")" << std::endl; + std::cout << "Track exit from right" << std::endl; + } + } + } + // enter from left + else if (yleft < (y + halfcrysize) && yleft > (y - halfcrysize)) { + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + path = TMath::Sqrt((x - halfcrysize - xlow) * (x - halfcrysize - xlow) + + (yleft - (y - halfcrysize)) * (yleft - (y - halfcrysize))); + } else { // track exit from right + path = halfcrysize * 2 * TMath::Sqrt(1 + m * m); + } + } + + else if (yright < (y + halfcrysize) && yright > (y - halfcrysize)) { + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + path = TMath::Sqrt((x + halfcrysize - xlow) * (x + halfcrysize - xlow) + + (yright - (y - halfcrysize)) * (yright - (y - halfcrysize))); + } + } + + else { + if (diag == 1) + std::cout << "can't find a path" << std::endl; + } + + if (diag == 1) + std::cout << "path is: " << path << std::endl; + + return path; +} // end findpath + +CaloCosmicEnecalib::LangausResult CaloCosmicEnecalib::fitLangaus(TH1F* hist, int refit, int iter, + float oldwidth, float oldsigma) { + + LangausResult res; + + // discarding histograms with not enough entries + if (hist->GetEntries() < 10) + return res; + + // set initial parameters + int binMax = hist->GetMaximumBin(); + float maxVal = hist->GetBinContent(binMax); + float mean = hist->GetBinCenter(binMax); + float rms = hist->GetRMS(); + float norm = hist->GetEntries(); + + // Range for Gaussian fit + float x_lo = mean - rms / 2.; + float x_hi = mean + rms / 2.; + + // Find FWHM + int bin1 = binMax; + int bin2 = binMax; + while ((bin1 > 1) && (hist->GetBinContent(bin1) > maxVal / 2.)) + bin1--; + while ((bin2 < hist->GetNbinsX()) && (hist->GetBinContent(bin2) > maxVal / 2.)) + bin2++; + double FWHM = hist->GetBinCenter(bin2) - hist->GetBinCenter(bin1); + + // Perform Gaussian fit + mygaus->SetParameters(norm, mean, FWHM / 2.); + hist->Fit("mygaus", "Q0", "", x_lo, x_hi); + double mpvVal = mygaus->GetParameter(1); + + // Range for Langaus fit + double xmin = 0.; + double xmax; + if (!_useMeV) + xmax = 2100.; + else + xmax = 100.; + + mylang->SetRange(0., xmax); + + mylang->SetParName(0, "Width"); + mylang->SetParName(1, "MPV"); + mylang->SetParName(2, "Norm"); + mylang->SetParName(3, "Sigma"); + + // Initial guess for parameters in Langaus fit + if (refit == 0) + mylang->SetParameters(FWHM / 2., mpvVal, norm, FWHM / 2.); + else { + // Refit with different initial params + float par0, par3; + FitFlags flags; + flags.fitQuality(oldwidth, oldsigma, hWidth, hSigma); + + // check Landau width + if (flags.badWidthLow) + par0 = hWidth->GetMean() + hWidth->GetRMS() * iter * 0.5; + else if (flags.badWidthHigh) + par0 = hWidth->GetMean() - hWidth->GetRMS() * iter * 0.5; + else + par0 = hWidth->GetMean(); + + // check Gaus sigma + if (flags.badSigmaLow) + par3 = hSigma->GetMean() + hSigma->GetRMS() * iter * 0.5; + else if (flags.badSigmaHigh) + par3 = hSigma->GetMean() - hSigma->GetRMS() * iter * 0.5; + else + par3 = hSigma->GetMean(); + + // Use new values as initial parameters + mylang->SetParameters(par0, mpvVal, norm, par3); + } + + // Limit parameters + mylang->SetParLimits(0, xmin, hWxmax); + mylang->SetParLimits(1, std::max(xmin, mpvVal - FWHM / 2.), std::min(xmax, mpvVal + FWHM / 2.)); + mylang->SetParLimits(3, xmin, hSxmax); + + hist->Fit("mylang", "Q", "", 0, xmax); + + if (_diagLevel > 0) { + std::cout << "fit complete " << std::endl; + } + + res.mpv = mylang->GetParameter(1); + res.mpvErr = mylang->GetParError(1); + res.width = mylang->GetParameter(0); + res.widthErr = mylang->GetParError(0); + res.sigma = mylang->GetParameter(3); + res.sigmaErr = mylang->GetParError(3); + res.chi2 = mylang->GetChisquare(); + res.ndf = mylang->GetNDF(); + res.nev = norm; + + return res; +} + +} // namespace mu2e + +DEFINE_ART_MODULE(mu2e::CaloCosmicEnecalib) diff --git a/CaloCalibration/CosmicsCalib/src/CaloCosmicEnergy_module.cc b/CaloCalibration/CosmicsCalib/src/CaloCosmicEnergy_module.cc new file mode 100644 index 0000000000..c38be9134b --- /dev/null +++ b/CaloCalibration/CosmicsCalib/src/CaloCosmicEnergy_module.cc @@ -0,0 +1,861 @@ +// ============================================================================= +// extrapolate energy, position and time for offline calorimeter calibration +// 12-Dec-2023 - P.Fedeli & S.Giovannella +// +// Starting code for retrieving calo information: +// CaloEnergy & caloT0alig module from P.Fedeli & S.Giovannella +// +// ============================================================================= + +// includes + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Provenance.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/Selector.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileDirectory.h" +#include "art_root_io/TFileService.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Sequence.h" + +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CalorimeterGeom/inc/CrystalMapper.hh" +#include "Offline/DataProducts/inc/CaloConst.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" + +#include "messagefacility/MessageLogger/MessageLogger.h" + +// C++ includes +#include +#include +#include +#include + +// ROOT includes +#include "TF1.h" +#include "TFitResult.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TH1F.h" +#include "TH2F.h" + +namespace mu2e { + +class CaloCosmicEnergy : public art::EDAnalyzer { + +public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom diagLevel{Name("diagLevel"), Comment(""), 0}; + fhicl::Atom CutNCryHit{Name("CutNCryHit"), Comment(""), 3}; + fhicl::Atom CutEnergyDep{Name("CutEnergyDep"), Comment(""), 10.}; + fhicl::Atom CutChi2Norm{Name("CutChi2Norm"), Comment(""), 2.5}; + fhicl::Atom CaloClusterTag{ + Name("CaloClusterTag"), Comment("Tag for Calorimeter cluster collection"), art::InputTag()}; + fhicl::Atom CaloHitTag{ + Name("CaloHitTag"), Comment("Tag for Calorimeter hit collection"), art::InputTag()}; + fhicl::Atom OutCalibFile{Name("OutCalibFile"), + Comment("Name for output .dat MIP calibration file"), + "calib_parameters.dat"}; + fhicl::Atom OutCalibFileAsym{Name("OutCalibFileAsym"), + Comment("Name for output .dat MIP asymmetry file"), + "npe_noise.dat"}; + }; + + explicit CaloCosmicEnergy(const art::EDAnalyzer::Table& config); + virtual ~CaloCosmicEnergy() {}; + void beginJob() override; + void beginRun(art::Run const& run) override; + + void analyze(art::Event const& event) override; + + void endJob() override; + + static double langaufun(Double_t* x, Double_t* par); + + float findpath(float m, float q, float x, float y); + +private: + art::ServiceHandle tfs; + art::TFileDirectory tfdir = tfs->mkdir("All_tracks"); + art::TFileDirectory tfdirv = tfs->mkdir("Vertical_tracks"); + art::TFileDirectory tfdird = tfs->mkdir("Diagonal_tracks"); + art::TFileDirectory tfdirfp = tfs->mkdir("All_normalized_tracks"); + art::TFileDirectory tfdir_res = tfs->mkdir("MPV_distributions"); + art::TFileDirectory tfdir_cryalr = tfs->mkdir("Asymmetry_lr_Crystals"); + art::TFileDirectory tfdir_alr = tfs->mkdir("Analysis_asymmetry_lr"); + + // calo parameters + static constexpr int nSiPMs = CaloConst::_nSiPMPerCrystal; + static constexpr int nCrystals = CaloConst::_nCrystalPerDisk; + static constexpr int nDisks = CaloConst::_nDisk; + static constexpr int nROchan = nDisks * nCrystals * nSiPMs; + const Calorimeter* cal; // For calorimeter geometry + float cryDim; // Total crystal dimension + float MaxDxVertical; // Max Dx value for vertical tracks + static constexpr int Erange = 5; // width of energy range in MeV for the sigma asymmetry variable + static constexpr int Ebin = 11; // number of ranges for the sigma asymmetry variable + + // fcl parameters + int _diagLevel; + int CutNCryHit; //> + float CutEnergyDep; // > + float CutChi2Norm; // < + art::ProductToken _caloClusterToken; + art::ProductToken _caloHitToken; + + TH1F* hSiPM[nDisks * nCrystals][nSiPMs] = {nullptr}; + TH1F* hSiPMv[nDisks * nCrystals][nSiPMs] = {nullptr}; + TH1F* hSiPMd[nDisks * nCrystals][nSiPMs] = {nullptr}; + TH1F* hSiPMfp[nDisks * nCrystals][nSiPMs] = {nullptr}; + + TH1F* MPV = nullptr; + TH1F* MPVv = nullptr; + TH1F* MPVd = nullptr; + TH1F* MPVfp = nullptr; + + TH1F* LR_all = nullptr; + TH2F* ALR_all = nullptr; + TH1F* LR[Ebin] = {nullptr}; + TH1F* ALR[Ebin] = {nullptr}; + TH1F* CryALR[nDisks * nCrystals][Ebin] = {nullptr}; + + TH1F* CryNpe = nullptr; + TH1F* CrySigNoise = nullptr; + TH1F* CryNpeChi2 = nullptr; + + float Energy_band[Ebin], counter_energy_band[Ebin], Cry_Energy_band[nDisks * nCrystals][Ebin], + Cry_counter_energy_band[nDisks * nCrystals][Ebin]; + TF1* gaussianFit = new TF1("gaussianFit", "gaus"); + + std::ofstream _outputfile; + std::ofstream _outputfileAsym; + + // helper variables + TF1* myPoly = new TF1("myPoly", "pol1", -600., 600.); // polinomial to fit the track + TF1* mylang = new TF1("mylang", langaufun, 15., 50., 4); // langaus function + +}; // end class + +CaloCosmicEnergy::CaloCosmicEnergy(const art::EDAnalyzer::Table& config) : + art::EDAnalyzer(config), _diagLevel(config().diagLevel()), CutNCryHit(config().CutNCryHit()), + CutEnergyDep(config().CutEnergyDep()), CutChi2Norm(config().CutChi2Norm()), + _caloClusterToken(consumes(config().CaloClusterTag())), + _caloHitToken(consumes(config().CaloHitTag())), + _outputfile(config().OutCalibFile()), _outputfileAsym(config().OutCalibFileAsym()) { + + if (!_outputfile.is_open()) { + throw cet::exception("CaloCosmicEnergy") + << "ERROR! Cannot open output file " << config().OutCalibFile() << std::endl; + } + if (!_outputfileAsym.is_open()) { + throw cet::exception("CaloCosmicEnergy") + << "ERROR! Cannot open output file " << config().OutCalibFileAsym() << std::endl; + } +} + +void CaloCosmicEnergy::beginRun(art::Run const& run) { + + art::ServiceHandle geom; + if (!geom->hasElement()) { + throw cet::exception("NO-CALO-GEOM") << "CaloEnergy: Calorimeter geometry not found"; + } + + GeomHandle ch; + cal = ch.get(); + cryDim = cal->G4Info().get("crystalXYLength") + + 2. * cal->G4Info().get("wrapperThickness"); + MaxDxVertical = cryDim * 1.1; +} + +void CaloCosmicEnergy::beginJob() { + if (_diagLevel > 0) + std::cout << "CaloCosmicEnergy: Entering beginJob" << std::endl; + + for (int iCry = 0; iCry < nDisks * nCrystals; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + hSiPM[iCry][iSiPM] = tfdir.make( + Form("Crystal_%i_SiPM_%i", iCry, iSiPM), + Form("Energy deposited in SiPM %i crystal %d", iSiPM, iCry), 70, 15., 50.); + hSiPMv[iCry][iSiPM] = tfdirv.make( + Form("Crystal_%i_SiPM_%i_v", iCry, iSiPM), + Form("Energy deposited in SiPM %i crystal %d only vertical tracks", iSiPM, iCry), 70, 15., + 50.); + hSiPMd[iCry][iSiPM] = tfdird.make( + Form("Crystal_%i_SiPM_%i_d", iCry, iSiPM), + Form("Energy deposited in SiPM %i crystal %d only diagonal tracks", iSiPM, iCry), 70, 15., + 50.); + hSiPMfp[iCry][iSiPM] = tfdirfp.make( + Form("Crystal_%i_SiPM_%i_fp", iCry, iSiPM), + Form("Energy deposited in SiPM %i crystal %d normalized track length", iSiPM, iCry), 70, + 15., 50.); + } + for (int ibin = 0; ibin < Ebin; ibin++) { + CryALR[iCry][ibin] = new TH1F(Form("Asymmetry_Crystal_%i_band_%i", iCry, ibin), + Form("Asymmetry of crystal %i in band [%i, %i)", iCry, + (ibin + 1) * Erange, (ibin + 1) * Erange + Erange), + 200, -1., 1.); + } + } + MPV = tfdir_res.make("MPV_distr", "MPV distribution", 400, 10., 30.); + MPVv = + tfdir_res.make("MPV_v_distr", "MPV distribution only vertical tracks", 400, 10., 30.); + MPVd = + tfdir_res.make("MPV_d_distr", "MPV distribution only diagonal tracks", 400, 10., 30.); + MPVfp = + tfdir_res.make("MPV_fp_distr", "MPV distribution all normalized tracks", 400, 10., 30.); + ALR_all = + tfdir_alr.make("ALR", "Asymmetry left right all energies", 600, 0., 600., 100, -1., 1.); + LR_all = + tfdir_alr.make("ALR_All_Energies", "Asymmetry Left Right all energies", 200, -1, 1.); + for (int i = 0; i < Ebin; i++) { + LR[i] = tfdir_alr.make( + Form("LR_ene_band_%i", i), + Form("Ratio Left Right Energy band [%i - %i) MeV", i * Erange, (i + 1) * Erange), 100, 0.5, + 1.5); + ALR[i] = tfdir_alr.make( + Form("ALR_ene_band_%i", i), + Form("Asymmetry Left-Right Energy band [%i - %i) MeV", i * Erange, (i + 1) * Erange), 200, + -1., 1.); + } +} // end begin job + +void CaloCosmicEnergy::analyze(art::Event const& event) { + + // art::Handle caloClustersHandle; + // event.getByLabel(_caloClusterTag, caloClustersHandle); + // const CaloClusterCollection& caloClusters(*caloClustersHandle); + auto const& caloClusters = event.getProduct(_caloClusterToken); + int nCluster = caloClusters.size(); + + for (int iClu = 0; iClu < nCluster; iClu++) { + int ncry = 0; // number of crystals passing the selection + float Dx = 1400.; // to calculate the spread alond x and y axis + float Dy = 1400.; + if (_diagLevel > 0) { + std::cout << "iClu: " << iClu << std::endl; + } + + // Creating arrays to store the information for crystals passing the selection + std::array PosX; + std::array PosY; + std::array ErrPos; + ErrPos.fill(9.81f); // Error for a uniform distribution: 34/sqrt(12) ( 34 = crystal dimension ) + std::array IDs; + std::array path; + std::array whichHit; + + /* + float PosY [nCrystals] = {0.0f}; + float ErrPos [nCrystals] = {9.81f}; // Error for a uniform distribution: 34/sqrt(12) ( 34 = + crystal dimension ) int IDs [nCrystals] = {0}; float path [nCrystals] = {0.0f}; int + whichHit[nCrystals] = {0}; + */ + + // Loop over crystals + for (long unsigned int iCry = 0; iCry < caloClusters[iClu].caloHitsPtrVector().size(); iCry++) { + // cut on energy of crystal in cluster + if (caloClusters[iClu].caloHitsPtrVector()[iCry].get()->energyDep() > CutEnergyDep) { + // counter n crystal above 10 MeV + PosX[ncry] = + cal->mu2eToDiskFF( + caloClusters[iClu].diskID(), + cal->crystal(caloClusters[iClu].caloHitsPtrVector()[iCry].get()->crystalID()) + .position()) + .getX(); + PosY[ncry] = + cal->mu2eToDiskFF( + caloClusters[iClu].diskID(), + cal->crystal(caloClusters[iClu].caloHitsPtrVector()[iCry].get()->crystalID()) + .position()) + .getY(); + IDs[ncry] = caloClusters[iClu].caloHitsPtrVector()[iCry].get()->crystalID(); + whichHit[ncry] = iCry; + ncry++; + } + + } // end loop iCry in cluster + + // cut on ncrystalhit + if (ncry > CutNCryHit) { + + float min_x = PosX[0]; + float max_x = PosX[0]; // sbsitture xlictyposx->posX + + // bubble sort to find Xmax and Xmin + for (int h = 1; h < ncry; h++) { + if (PosX[h] > max_x) { + max_x = PosX[h]; + // cout << "Max: " << max << endl; + } + if (PosX[h] < min_x) { + min_x = PosX[h]; + // cout << "Min: " << min << endl; + } + Dx = abs(max_x - min_x); + } + + float min_y = PosY[0]; + float max_y = PosY[0]; + + // bubble sort to find Xmax and Xmin + for (int h = 1; h < ncry; h++) { + if (PosY[h] > max_y) { + max_y = PosX[h]; + // cout << "Max: " << max << endl; + } + if (PosY[h] < min_y) { + min_y = PosY[h]; + // cout << "Min: " << min << endl; + } + Dy = abs(max_y - min_y); + } + + if (_diagLevel > 0) { + std::cout << "dx: " << Dx << " dy: " << Dy << std::endl; + } + + // geometry track + TGraph* gfirst = new TGraph(ncry, PosX.data(), PosY.data()); + auto gfit = std::make_unique(ncry, PosX.data(), PosY.data(), ErrPos.data(), + ErrPos.data()); + gfit->SetMarkerSize(0.9); + gfit->SetMarkerStyle(8); + + float chi2norm = 0.; + + TFitResultPtr fitresult = gfirst->Fit("myPoly", "SQ"); + TFitResultPtr fit = gfit->Fit("myPoly", "SQ"); + // 2 fits help convergence + + chi2norm = fit->Chi2() / fit->Ndf(); + + // begin findpath + + float TrkSlope = -1; + TrkSlope = myPoly->GetParameter(1); + if (Dx < 33) { // if vertical tracks fit could be wrong + TrkSlope = 144.; // ca 89.6 deg + } + + for (int h = 0; h < ncry; h++) { + if (_diagLevel > 0) { + std::cout << "using m:" << TrkSlope << " q: " << myPoly->GetParameter(0) + << " x: " << PosX[h] << " y: " << PosY[h] << " ene: " + << caloClusters[iClu].caloHitsPtrVector()[whichHit[h]].get()->energyDep() + << std::endl; + } + path[h] = findpath(TrkSlope, myPoly->GetParameter(0), PosX[h], PosY[h]); + } // end findpath + + // only verical tracks + if (Dx < MaxDxVertical) { + for (int iCry = 0; iCry < ncry; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + hSiPMv[IDs[iCry]][iSiPM]->Fill(caloClusters[iClu] + .caloHitsPtrVector()[whichHit[iCry]] + .get() + ->recoCaloDigis()[iSiPM] + .get() + ->energyDep()); + } + } + } + + // cut on diagonal tracks + if (chi2norm < 0.01) { + for (int iCry = 0; iCry < ncry; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + hSiPMd[IDs[iCry]][iSiPM]->Fill(caloClusters[iClu] + .caloHitsPtrVector()[whichHit[iCry]] + .get() + ->recoCaloDigis()[iSiPM] + .get() + ->energyDep()); + } + } + } + + // cut on vertical tracks + good fit chi2 -> all tracks + if (Dx < 35 || chi2norm < CutChi2Norm) { + for (int iCry = 0; iCry < ncry; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + hSiPM[IDs[iCry]][iSiPM]->Fill(caloClusters[iClu] + .caloHitsPtrVector()[whichHit[iCry]] + .get() + ->recoCaloDigis()[iSiPM] + .get() + ->energyDep()); + // normalized tracks by fidpaths.h + hSiPMfp[IDs[iCry]][iSiPM]->Fill(caloClusters[iClu] + .caloHitsPtrVector()[whichHit[iCry]] + .get() + ->recoCaloDigis()[iSiPM] + .get() + ->energyDep() * + cryDim / path[iCry]); + } + } + } // chi2 cut + + // asymmetry left-right and ratio left/right-> we use only trigger selection + for (long unsigned int iCry = 0; iCry < caloClusters[iClu].caloHitsPtrVector().size(); + iCry++) { + float left_sipm = caloClusters[iClu] + .caloHitsPtrVector()[iCry] + .get() + ->recoCaloDigis()[0] + .get() + ->energyDep(); + float right_sipm = caloClusters[iClu] + .caloHitsPtrVector()[iCry] + .get() + ->recoCaloDigis()[1] + .get() + ->energyDep(); + float sipm_mean_e = caloClusters[iClu].caloHitsPtrVector()[iCry].get()->energyDep(); + int crystal_id = caloClusters[iClu].caloHitsPtrVector()[iCry].get()->crystalID(); + LR_all->Fill((left_sipm - right_sipm) / (left_sipm + right_sipm)); + ALR_all->Fill(sipm_mean_e, (left_sipm - right_sipm) / (left_sipm + right_sipm)); + if (sipm_mean_e >= 0. && sipm_mean_e < Ebin * Erange) { + int whichband = sipm_mean_e / Erange; + Energy_band[whichband] += sipm_mean_e; + counter_energy_band[whichband]++; + LR[whichband]->Fill(left_sipm / right_sipm); + ALR[whichband]->Fill((left_sipm - right_sipm) / (left_sipm + right_sipm)); + CryALR[crystal_id][whichband]->Fill((left_sipm - right_sipm) / (left_sipm + right_sipm)); + Cry_Energy_band[crystal_id][whichband] += sipm_mean_e; + Cry_counter_energy_band[crystal_id][whichband]++; + } + } + + delete gfirst; + + } // end ncrystal cut + + } // end loop iClu in ncluster +} // end enalyze + +void CaloCosmicEnergy::endJob() { + std::cout << "CaloCosmicEnergy: entering endjob" << std::endl; + TGraphErrors *NpeLR, *NpeALR, *sigmaELR, *sigmaEALR, *CryNpeALR[nDisks * nCrystals]; + + // fit params + int redo[nDisks * nCrystals][nSiPMs]; + int ibin = 0; + float Xmean = 0; + float Xsigma = 0; + float params[nDisks * nCrystals][nSiPMs][5]; + + for (int k = 0; k < nDisks * nCrystals; k++) { + for (int j = 0; j < nSiPMs; j++) { + redo[k][j] = 0; + for (int i = 0; i < 5; i++) { + params[k][j][i] = -1; + } + } + } + for (int iCry = 0; iCry < nDisks * nCrystals; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + + // set initial parameters + + ibin = hSiPM[iCry][iSiPM]->GetMaximumBin(); + Xmean = hSiPM[iCry][iSiPM]->GetBinCenter(ibin); + Xsigma = hSiPM[iCry][iSiPM]->GetRMS(); + + mylang->SetParameters(Xsigma, Xmean, hSiPM[iCry][iSiPM]->GetMaximum(), 1); + hSiPM[iCry][iSiPM]->Fit("mylang", "Q", "", 10, 50); + + if (_diagLevel > 0) { + std::cout << "fit complete " << std::endl; + } + if ((mylang->GetChisquare() / mylang->GetNDF()) < 1.8 && + (hSiPM[iCry][iSiPM]->GetMean() - mylang->GetParameter(1)) < 10.) { + if (_diagLevel > 0) { + std::cout << "Crystal: " << iCry << " SiPM: " << iSiPM + << " MPV: " << mylang->GetParameter(1) << std::endl; + } + + MPV->Fill(mylang->GetParameter(1)); + + params[iCry][iSiPM][0] = mylang->GetParameter(1); + params[iCry][iSiPM][1] = mylang->GetParError(1); + params[iCry][iSiPM][2] = mylang->GetParameter(0); + params[iCry][iSiPM][3] = mylang->GetParError(0); + params[iCry][iSiPM][4] = mylang->GetChisquare() / mylang->GetNDF(); + + } else { + redo[iCry][iSiPM] = 1; + } + + } // end sipms loop + } // end crystal loop + if (_diagLevel > 0) { + std::cout << " I am redoing fits" << std::endl; + } + + // redo wrong fit + for (int iCry = 0; iCry < nDisks * nCrystals; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + if (redo[iCry][iSiPM] == 1) { + ibin = hSiPM[iCry][iSiPM]->GetMaximumBin(); + Xmean = hSiPM[iCry][iSiPM]->GetBinCenter(ibin); + Xsigma = hSiPM[iCry][iSiPM]->GetRMS(); + + // set wider parameter + mylang->SetParameters(Xsigma * 1.2, Xmean, hSiPM[iCry][iSiPM]->GetMaximum(), 1.1); + + hSiPM[iCry][iSiPM]->Fit("mylang", "Q", "", 10, 50); + + MPV->Fill(mylang->GetParameter(1)); + + params[iCry][iSiPM][0] = mylang->GetParameter(1); + params[iCry][iSiPM][1] = mylang->GetParError(1); + params[iCry][iSiPM][2] = mylang->GetParameter(0); + params[iCry][iSiPM][3] = mylang->GetParError(0); + params[iCry][iSiPM][4] = mylang->GetChisquare() / mylang->GetNDF(); + + if (_diagLevel > 0) { + std::cout << "Crystal: " << iCry << " SiPM: " << iSiPM + << " MPV: " << mylang->GetParameter(1) << std::endl; + } + } + } + } + _outputfile << "#Cry #SiPM #MPV #MPV err #sigm #sigma err #chi2/ndf \n\n"; + for (int iCry = 0; iCry < nDisks * nCrystals; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + _outputfile << Form("%i %i %.3f %.3f %.3f %.3f %.3f\n", iCry, iSiPM, params[iCry][iSiPM][0], + params[iCry][iSiPM][1], params[iCry][iSiPM][2], params[iCry][iSiPM][3], + params[iCry][iSiPM][4]); + } + _outputfile << "\n"; + } + + _outputfile.close(); + + // fits for vertical diagonal and nromalized tracks + for (int iCry = 0; iCry < nDisks * nCrystals; iCry++) { + for (int iSiPM = 0; iSiPM < nSiPMs; iSiPM++) { + + // set initial parameters + + ibin = hSiPMv[iCry][iSiPM]->GetMaximumBin(); + Xmean = hSiPMv[iCry][iSiPM]->GetBinCenter(ibin); + Xsigma = hSiPMv[iCry][iSiPM]->GetRMS(); + + mylang->SetParameters(Xsigma, Xmean, hSiPMv[iCry][iSiPM]->GetMaximum(), 1); + hSiPMv[iCry][iSiPM]->Fit("mylang", "Q", "", 10, 50); + MPVv->Fill(mylang->GetParameter(1)); + + ibin = hSiPMd[iCry][iSiPM]->GetMaximumBin(); + Xmean = hSiPMd[iCry][iSiPM]->GetBinCenter(ibin); + Xsigma = hSiPMd[iCry][iSiPM]->GetRMS(); + + mylang->SetParameters(Xsigma, Xmean, hSiPMd[iCry][iSiPM]->GetMaximum(), 1); + hSiPMd[iCry][iSiPM]->Fit("mylang", "Q", "", 10, 50); + MPVd->Fill(mylang->GetParameter(1)); + + ibin = hSiPMfp[iCry][iSiPM]->GetMaximumBin(); + Xmean = hSiPMfp[iCry][iSiPM]->GetBinCenter(ibin); + Xsigma = hSiPMfp[iCry][iSiPM]->GetRMS(); + + mylang->SetParameters(Xsigma, Xmean, hSiPMfp[iCry][iSiPM]->GetMaximum(), 1); + hSiPMfp[iCry][iSiPM]->Fit("mylang", "Q", "", 10, 50); + MPVfp->Fill(mylang->GetParameter(1)); + } + } + MPV->Fit("gaus", "Q", "", 10., 30.); + MPVv->Fit("gaus", "Q", "", 10., 30.); + MPVd->Fit("gaus", "Q", "", 10., 30.); + MPVfp->Fit("gaus", "Q", "", 10., 30.); + + // start Asymmetry L-R analysis + + LR_all->Fit("gaus", "Q", "", -1., 1.); + + // helper variables + float sigmaLR[Ebin], sigmaLRerr[Ebin], sigmaALR[Ebin], sigmaALRerr[Ebin], NpeLRvar[Ebin], + NpeALRvar[Ebin], eband_mean[Ebin], errEne[Ebin], errLR[Ebin], errALR[Ebin], + CrySigmaALR[nDisks * nCrystals][Ebin], CrySigmaALRErr[nDisks * nCrystals][Ebin], + Cryeband_mean[nDisks * nCrystals][Ebin], CryEneErr[nDisks * nCrystals][Ebin]; + + // function to fit Asymmetry L-R(ALR) and L/R(LR) + TF1* functional_form_ALR = + new TF1("functional_form_ALR", "TMath::Sqrt(0.5*(1/([0]*x)+([1]/x)^2))", 0., 100.); + TF1* functional_form_LR = + new TF1("functional_form_LR", "TMath::Sqrt(2*(1/([0]*x))+([1]/x)^2)", 0., 100.); + + // peparing histogram for Npe e SigNoise distribution + CryNpe = tfdir_alr.make("Npe_distr", "Distribution of Npe of all crystalss", 120, 0., 60.); + CrySigNoise = tfdir_alr.make("SigNoise_distr", "Distribution of SigNoise of all crystals", + 80, 0., 0.4); + CryNpeChi2 = + tfdir_alr.make("CryNpeChi2", "Distribution of #chi2/ndof of all crystals", 100, 0., 50); + + // initializating string for npe-noise file dat + _outputfileAsym << "#Cry #Npe/MeV #Npe/MeVerr #Noise #Noiseerr \n\n"; + + // preparing variables to make tgraph for every crystal + for (int icry = 0; icry < nDisks * nCrystals; icry++) { + for (int ibin = 0; ibin < Ebin; ibin++) { + CryALR[icry][ibin]->Fit("gaussianFit", "Q", "", -1., 1.); + CrySigmaALR[icry][ibin] = gaussianFit->GetParameter(2); + CrySigmaALRErr[icry][ibin] = gaussianFit->GetParError(2); + + // avoid infinities due low statistics + if (Cry_counter_energy_band[icry][ibin] == 0) { + Cry_counter_energy_band[icry][ibin] = 1; + } + Cryeband_mean[icry][ibin] = Cry_Energy_band[icry][ibin] / Cry_counter_energy_band[icry][ibin]; + CryEneErr[icry][ibin] = 1 / TMath::Sqrt(Cry_counter_energy_band[icry][ibin]); + } // end ibin for + + // tgraph for single chanel + CryNpeALR[icry] = tfdir_cryalr.make(Ebin, Cryeband_mean[icry], CrySigmaALR[icry], + CryEneErr[icry], CrySigmaALRErr[icry]); + CryNpeALR[icry]->SetMarkerStyle(8); + CryNpeALR[icry]->SetMarkerSize(0.9); + functional_form_ALR->SetParameters(30., 0.2); + + // fit single channel + CryNpeALR[icry]->Fit("functional_form_ALR", "Q"); + + // fill npe-nois file dat with crsystal; npe/mev +/- err; noise +/- err + _outputfileAsym << Form("%i %.3f %.3f %.3f %.3f\n", icry, functional_form_ALR->GetParameter(0), + functional_form_ALR->GetParError(0), + functional_form_ALR->GetParameter(1), + functional_form_ALR->GetParError(1)); + // eventualy chi2/ndf cut + if (functional_form_ALR->GetChisquare() / functional_form_ALR->GetNDF() > 0) { + CryNpeALR[icry]->Write(Form("Cry_%i_NPe", icry)); + CryNpe->Fill(functional_form_ALR->GetParameter(0)); // fill Npe + CryNpe->Fit("gaus", "Q"); + CrySigNoise->Fill(functional_form_ALR->GetParameter(1)); // fill SigNoise + CrySigNoise->Fit("gaus", "Q"); + } + CryNpeChi2->Fill(functional_form_ALR->GetChisquare() / functional_form_ALR->GetNDF()); + } // end icry for + + _outputfileAsym.close(); + for (int j = 0; j < Ebin; j++) { + eband_mean[j] = Energy_band[j] / counter_energy_band[j]; + errEne[j] = 1 / TMath::Sqrt(counter_energy_band[j]); + + LR[j]->Fit("gaussianFit", "Q", "", 0.5, 1.5); + sigmaLR[j] = gaussianFit->GetParameter(2); + sigmaLRerr[j] = gaussianFit->GetParError(2); + NpeLRvar[j] = 2 / gaussianFit->GetParameter(2) / gaussianFit->GetParameter(2) / eband_mean[j]; + errLR[j] = 4 * gaussianFit->GetParError(2) / gaussianFit->GetParameter(2) / + gaussianFit->GetParameter(2) / gaussianFit->GetParameter(2) / eband_mean[j]; + + ALR[j]->Fit("gaussianFit", "Q", "", -1., 1.); + sigmaALR[j] = gaussianFit->GetParameter(2); + sigmaALRerr[j] = gaussianFit->GetParError(2); + NpeALRvar[j] = + 1 / gaussianFit->GetParameter(2) / gaussianFit->GetParameter(2) / eband_mean[j] / 2; + errALR[j] = 0.5 * gaussianFit->GetParError(2) / gaussianFit->GetParameter(2) / + gaussianFit->GetParameter(2) / gaussianFit->GetParameter(2) / eband_mean[j]; + } + + // analsis for all cristals togheter + NpeLR = tfdir_alr.make(Ebin, eband_mean, NpeLRvar, errEne, errLR); + NpeALR = tfdir_alr.make(Ebin, eband_mean, NpeALRvar, errEne, errALR); + sigmaELR = tfdir_alr.make(Ebin, eband_mean, sigmaLR, errEne, sigmaLRerr); + sigmaEALR = tfdir_alr.make(Ebin, eband_mean, sigmaALR, errEne, sigmaALRerr); + + NpeLR->SetMarkerSize(0.9); + NpeLR->SetMarkerStyle(8); + NpeALR->SetMarkerSize(0.9); + NpeALR->SetMarkerStyle(8); + sigmaELR->SetMarkerSize(0.9); + sigmaELR->SetMarkerStyle(8); + sigmaEALR->SetMarkerSize(0.9); + sigmaEALR->SetMarkerStyle(8); + NpeLR->Fit("pol0", "Q"); + NpeALR->Fit("pol0", "Q"); + functional_form_LR->SetParameters(30., 0.2); + functional_form_ALR->SetParameters(30., 0.2); + sigmaELR->Fit("functional_form_LR", "Q"); + sigmaEALR->Fit("functional_form_ALR", "Q"); + NpeLR->Write("NpeLR"); + NpeALR->Write("NpeALR"); + sigmaELR->Write("sigmaELR"); + sigmaEALR->Write("sigmaEALR"); + + if (_diagLevel > 0) { + std::cout << "[0]: " << functional_form_ALR->GetParameter(0) + << "[1]: " << functional_form_ALR->GetParameter(1); + std::cout << "[0]: " << functional_form_LR->GetParameter(0) + << "[1]: " << functional_form_LR->GetParameter(1); + } + + delete gaussianFit; + delete functional_form_ALR; + delete functional_form_LR; + for (int icry = 0; icry < nDisks * nCrystals; icry++) { + for (int ibin = 0; ibin < Ebin; ibin++) { + delete CryALR[icry][ibin]; + } + } + +} // end endjob + +double CaloCosmicEnergy::langaufun(double* x, double* par) { + + // Fit parameters: + // par[0]=Width (scale) parameter of Landau density + // par[1]=Most Probable (MP, location) parameter of Landau density + // par[2]=Total area (integral -inf to inf, normalization constant) + // par[3]=Width (sigma) of convoluted Gaussian function + // + // In the Landau distribution (represented by the CERNLIB approximation), + // the maximum is located at x=-0.22278298 with the location parameter=0. + // This shift is corrected within this function, so that the actual + // maximum is identical to the MP parameter. + + // Numeric constants + Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) + Double_t mpshift = -0.22278298; // Landau maximum location + + // Control constants + Double_t np = 100.0; // number of convolution steps + Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas + + // Variables + Double_t xx; + Double_t mpc; + Double_t fland; + Double_t sum = 0.0; + Double_t xlow, xupp; + Double_t step; + Double_t i; + + // MP shift correction + mpc = par[1] - mpshift * par[0]; + + // Range of convolution integral + xlow = x[0] - sc * par[3]; + xupp = x[0] + sc * par[3]; + + step = (xupp - xlow) / np; + + // Convolution integral of Landau and Gaussian by sum + for (i = 1.0; i <= np / 2; i++) { + xx = xlow + (i - .5) * step; + fland = TMath::Landau(xx, mpc, par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0], xx, par[3]); + + xx = xupp - (i - .5) * step; + fland = TMath::Landau(xx, mpc, par[0]) / par[0]; + sum += fland * TMath::Gaus(x[0], xx, par[3]); + } + + return (par[2] * step * sum * invsq2pi / par[3]); +} // end langaus + +float CaloCosmicEnergy::findpath(float m, float q, float x, float y) { + float xup = 0; + float yright = 0; + float path = 0; + float xlow = 0; + float yleft = 0; + + int diag = 0; + float halfcrysize = cryDim / 2.; + + if (m != 0) { + xup = ((y + halfcrysize) - q) / m; + xlow = ((y - halfcrysize) - q) / m; + yleft = m * (x - halfcrysize) + q; + yright = m * (x + halfcrysize) + q; + } + if (xup < (x + halfcrysize) && xup > (x - halfcrysize)) { // hit the top of crystal + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + path = halfcrysize * 2 * TMath::Sqrt(1 / (m * m) + 1); + if (diag == 1) { + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << xlow << ";" + << y - halfcrysize << ")" << std::endl; + } + } + if (xlow < (x - halfcrysize)) { + path = TMath::Sqrt((xup - (x - halfcrysize)) * (xup - (x - halfcrysize)) + + ((y + halfcrysize) - yleft) * ((y + halfcrysize) - yleft)); + if (diag == 1) { + std::cout << "Track exit from left" << std::endl; + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << x - halfcrysize + << ";" << yleft << ")" << std::endl; + } + } + if (xlow > (x + halfcrysize)) { + path = TMath::Sqrt((xup - (x + halfcrysize)) * (xup - (x + halfcrysize)) + + ((y + halfcrysize) - yright) * ((y + halfcrysize) - yright)); + if (diag == 1) { + std::cout << "x0,y0: (" << xup << ";" << y + halfcrysize << ") x1,y1: (" << x + halfcrysize + << ";" << yright << ")" << std::endl; + std::cout << "Track exit from right" << std::endl; + } + } + } + // enter from left + else if (yleft < (y + halfcrysize) && yleft > (y - halfcrysize)) { + if (diag == 1) { + std::cout << "track enter from left and "; + } + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + if (diag == 1) { + std::cout << " it exit bot " << std::endl; + std::cout << "x0,y0: (" << x - halfcrysize << ";" << yleft << ") x1,y1: (" << xlow << ";" + << y - halfcrysize << ")" << std::endl; + } + path = TMath::Sqrt((x - halfcrysize - xlow) * (x - halfcrysize - xlow) + + (yleft - (y - halfcrysize)) * (yleft - (y - halfcrysize))); + } else { // track exit from right + if (diag == 1) { + std::cout << " it exit from right" << std::endl; + std::cout << "x0,y0: (" << x - halfcrysize << ";" << yleft << ") x1,y1: (" + << x + halfcrysize << ";" << yright << ")" << std::endl; + } + path = halfcrysize * 2 * TMath::Sqrt(1 + m * m); + } + } + + else if (yright < (y + halfcrysize) && yright > (y - halfcrysize)) { + if (diag == 1) { + std::cout << "track enter from right and "; + } + if (xlow < (x + halfcrysize) && xlow > (x - halfcrysize)) { + if (diag == 1) { + std::cout << " it exit bot " << std::endl; + std::cout << "x0,y0: (" << x + halfcrysize << ";" << yright << ") x1,y1: (" << xlow << ";" + << y - halfcrysize << ")" << std::endl; + } + path = TMath::Sqrt((x + halfcrysize - xlow) * (x + halfcrysize - xlow) + + (yright - (y - halfcrysize)) * (yright - (y - halfcrysize))); + } + } + + else { + if (diag == 1) { + std::cout << "can't find a path" << std::endl; + } + } + if (diag == 1) { + std::cout << "path is: " << path << std::endl; + } + + return path; +} // end findpath + +} // namespace mu2e + +DEFINE_ART_MODULE(mu2e::CaloCosmicEnergy) diff --git a/CaloCalibration/CosmicsCalib/src/SConscript b/CaloCalibration/CosmicsCalib/src/SConscript new file mode 100644 index 0000000000..b1c656f973 --- /dev/null +++ b/CaloCalibration/CosmicsCalib/src/SConscript @@ -0,0 +1,63 @@ +#!/usr/bin/env python +# +# Original author Rob Kutschke. +# + +import os, re +Import('env') + +Import('mu2e_helper') + +helper=mu2e_helper(env) + +rootlibs = env['ROOTLIBS'] + +extrarootlibs = [ 'Geom', 'TMVA' , 'Minuit' , 'XMLIO' ] + +mainlib = helper.make_mainlib ( [] ) + +helper.make_plugins( [ mainlib, + 'mu2e_Mu2eUtilities', + 'mu2e_CaloCluster', + 'mu2e_ParticleID', + 'mu2e_RecoDataProducts', + 'mu2e_TrackerConditions', + 'mu2e_GeometryService', + 'mu2e_CalorimeterGeom', + 'mu2e_CosmicRayShieldGeom', + 'mu2e_ExtinctionMonitorFNAL_Geometry', + 'mu2e_ProtonBeamDumpGeom', + 'mu2e_StoppingTargetGeom', + 'mu2e_TrackerGeom', + 'mu2e_GlobalConstantsService_GlobalConstantsService_service', + 'mu2e_GlobalConstantsService', + 'mu2e_ConfigTools', + 'mu2e_DataProducts', + 'mu2e_GeneralUtilities', + 'art_Framework_Core', + 'art_Framework_Principal', + 'art_Framework_Services_Registry', + 'art_root_io_tfile_support', + 'art_root_io_TFileService', + 'art_Framework_Services_Optional_RandomNumberGenerator_service', + 'art_Persistency_Common', + 'art_Persistency_Provenance', + 'art_Utilities', + 'canvas', + 'MF_MessageLogger', + 'fhiclcpp', + 'fhiclcpp_types', + 'tbb', + 'cetlib', + 'cetlib_except', + 'CLHEP', + rootlibs, + extrarootlibs, + 'xerces-c', + 'boost_filesystem', + ] ) + +# this tells emacs to view this file in python mode. +# Local Variables: +# mode:python +# End: diff --git a/CaloCalibration/CosmicsCalib/src/caloT0alig_module.cc b/CaloCalibration/CosmicsCalib/src/caloT0alig_module.cc new file mode 100644 index 0000000000..ac70b9cf09 --- /dev/null +++ b/CaloCalibration/CosmicsCalib/src/caloT0alig_module.cc @@ -0,0 +1,503 @@ +// ============================================================================= +// Filter module for offline calorimeter time calibration +// 18-Mar-2021 - S.Giovannella +// +// Starting code for retrieving calo information: CaloExample module +// +// This module performs a single step of an iterative procedure running on +// calorimeter readout channels: +// 1. selection of straight cosmic ray events with mip-like energy deposit +// 2. fit to the 2d linear trajectory in the XY plane with least-square method +// 3. common energy-weighted T0 subtracted to all readout channel +// 4. speed of light imposed between timing of different readouts +// 5. residuals evaluated for each readout channel, to be used as input for +// the next iteration +// Iterations are controlled by an external script +// Good events can be filtered to speed up the procedure +// Last iteration creates an output file to be uploaded on Condition DB +// ============================================================================= + +#include "art/Framework/Core/EDFilter.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Provenance.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/Selector.h" +#include "art_root_io/TFileService.h" +#include "canvas/Utilities/InputTag.h" +#include "cetlib_except/exception.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" + +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CalorimeterGeom/inc/DiskCalorimeter.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" + +#include "Offline/DataProducts/inc/CaloConst.hh" +#include "Offline/DataProducts/inc/CaloSiPMId.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/CaloHit.hh" + +#include +#include +#include + +#include "TDirectory.h" +#include "TH1F.h" +#include +#include + +namespace mu2e { + +class caloT0alig : public art::EDFilter { + +public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom caloHitCollection{Name("caloHitCollection"), + Comment("Calo Hit collection name")}; + fhicl::Atom caloClusterCollection{Name("caloClusterCollection"), + Comment("Calo cluster collection name")}; + fhicl::Atom cluHits{Name("cluHits"), Comment("Minimum number of hits in a cluster"), 3}; + fhicl::Atom cryEmin{Name("cryEmin"), Comment("Minimum Energy for crystal hits"), 15.}; + fhicl::Atom ncryCut{Name("ncryCut"), Comment("Minimum number of mip-like crystals"), 7}; + fhicl::Atom iteration{Name("iteration"), Comment("Iteration"), "first"}; + fhicl::Atom fileT0{Name("fileT0"), Comment("T0 input file")}; + fhicl::Atom fileTcor{Name("fileTcor"), Comment("T0 corrections input file")}; + fhicl::Atom diagLevel{Name("diagLevel"), Comment("Diag Level"), 0}; + }; + + explicit caloT0alig(const art::EDFilter::Table& config); + virtual ~caloT0alig() {} + + virtual void beginJob(); + virtual void endJob(); + virtual bool filter(art::Event& e); + +private: + art::InputTag _caloHitTag; + art::InputTag _caloClusterTag; + int _cluHits; + float _cryEmin; + int _ncryCut; + std::string _iteration; + std::ifstream _fileT0; + std::string _fileT0Name; + std::string _fileTcorName; + int _diagLevel; + int _nProcessed; + int _nFiltered; + + static constexpr double cvel = 299.792458; + + static constexpr int nDisks = CaloConst::_nDisk; // Number of calorimeter disks + static constexpr int nCrystals = CaloConst::_nCrystalPerDisk; // Number of crystals for each disk + static constexpr int nSiPMs = CaloConst::_nSiPMPerCrystal; // Number of SIPMs for each crystal + static constexpr int nROchan = nDisks * nCrystals * nSiPMs; + + float Tcor[nROchan] = {0}; + float Toff[nROchan] = {0}; + + TH1F* hTres[nROchan]; + TH1F *hcosTh[nDisks], *hcosThCut[nDisks]; + TH1F *hcosFit[nDisks], *hChi2[nDisks], *hDcosth[nDisks], *hYdif[nDisks]; + TH1F *hcvel[nDisks], *hEcell[nDisks], *hEmean[nDisks]; + TH1F *hNphi[nDisks], *hNphiEcut[nDisks]; + TH1F *hNhit, *hNhitEcut; +}; + +caloT0alig::caloT0alig(const art::EDFilter::Table& config) : + EDFilter{config}, _caloHitTag(config().caloHitCollection()), + _caloClusterTag(config().caloClusterCollection()), _cluHits(config().cluHits()), + _cryEmin(config().cryEmin()), _ncryCut(config().ncryCut()), _iteration(config().iteration()), + _fileT0Name(config().fileT0()), _fileTcorName(config().fileTcor()), _diagLevel(config().diagLevel()), + _nProcessed(0), _nFiltered(0) { + + _fileT0.open(_fileT0Name); + if (!_fileT0.is_open()) { + throw cet::exception("caloT0alig") + << "ERROR! Cannot open input file " << config().fileT0() << std::endl; + } +} + +// =========================================================================== +// Begin job: +// - Get first step of T0 corrections, i.e. mean values of time distribution +// for cosmic ray events +// - Get T0 corrections from previous iterations, if applicable +// - Histogram booking +// =========================================================================== +void caloT0alig::beginJob() { + + if (_diagLevel > 0) + std::cout << "caloT0alig: Entering beginJob" << std::endl; + + // Read T0-step0 corrections + // *** TO BE IMPLEMENTED ... waiting for real data *** + + /////////////////////////////////////////////////////////////////////// + // TEMPORARY for running on MC data: + // T0 smearing (+-0.5 ns) to simulate real data and test the procedure + /////////////////////////////////////////////////////////////////////// + + int iChanT0, nValT0 = 0; + float TvalT0; + + if (_fileT0.is_open()) { + while (_fileT0 >> iChanT0 >> TvalT0) { + if (iChanT0 < 0 || iChanT0 >= nROchan) { + throw cet::exception("caloT0alig") << "ERROR! Read invalid channel " << iChanT0 + << "from file " << _fileT0Name << std::endl; + } + Toff[iChanT0] = TvalT0; + if (_diagLevel > 1) + std::cout << "IdxT0 " << iChanT0 << " Toff " << Toff[iChanT0] << " " << std::endl; + nValT0++; + } + _fileT0.close(); + } else { + mf::LogError("INPUT-NOT-FOUND") + << "T0 file from previous iteration not found: " << _fileT0Name << std::endl; + } + /////////////////////////////////////////////////////////////// + + // Read T0 corrections from previous iteration, if applicable + + if (_iteration == "middle" || _iteration == "last") { + if (_diagLevel > 0) + std::cout << "Iteration " << _iteration << ": Reading T0 corrections from previous step " + << std::endl; + + int iChan, nVal = 0; + int Nev; + float Tval, Tmea, Tres, Chi2; + + std::ifstream _fileTcor(_fileTcorName); + if (_fileTcor.is_open()) { + while (_fileTcor >> iChan >> Tval >> Tmea >> Tres >> Chi2 >> Nev) { + if (iChan < 0 || iChan >= nROchan) { + throw cet::exception("caloT0alig") << "ERROR! Read invalid channel " << iChan + << "from file " << _fileTcorName << std::endl; + } + + Tcor[iChan] = Tval; + if (_diagLevel > 1) + std::cout << "Idx " << iChan << " Tcor " << Tcor[iChan] << " " << Tmea << std::endl; + nVal++; + } + _fileTcor.close(); + if (nVal != nROchan) { + mf::LogError("WRONG-NCHAN") + << "Wrong number of readout channels: " << nVal << " " << nROchan << std::endl; + } + } else { + if (!_fileTcor.is_open()) { + throw cet::exception("caloT0alig") + << "ERROR! Cannot open output file " << _fileTcorName << std::endl; + } + } + } + + // Book histograms + + art::ServiceHandle tfs; + + for (int iChan = 0; iChan < nROchan; iChan++) { + hTres[iChan] = + tfs->make(Form("Tres_%04i", iChan), Form("Tres %04i", iChan), 200, -2., 2.); + } + + for (int iDisk = 0; iDisk < nDisks; iDisk++) { + hcosTh[iDisk] = tfs->make(Form("cosTh_%i", iDisk), Form("cosTh %i", iDisk), 200, -1., 1.); + hcosThCut[iDisk] = + tfs->make(Form("cosThCut_%i", iDisk), Form("cosThCut %i", iDisk), 200, -1., 1.); + hcosFit[iDisk] = + tfs->make(Form("cosFit_%i", iDisk), Form("cosFit %i", iDisk), 200, -1., 1.); + hDcosth[iDisk] = + tfs->make(Form("Dcosth_%i", iDisk), Form("Dcosth %i", iDisk), 2000, -1., 1.); + hYdif[iDisk] = + tfs->make(Form("Ydif_%i", iDisk), Form("Ydif %i", iDisk), 200, -1000., 1000.); + hChi2[iDisk] = tfs->make(Form("Chi2_%i", iDisk), Form("Chi2 %i", iDisk), 200, 0., 100.); + hcvel[iDisk] = tfs->make(Form("cvel_%i", iDisk), Form("cvel %i", iDisk), 200, 100., 500.); + hEcell[iDisk] = + tfs->make(Form("Ecell_%01i", iDisk), Form("Ecell %i", iDisk), 200, 0., 100.); + hEmean[iDisk] = + tfs->make(Form("Emean_%01i", iDisk), Form("Emean %i", iDisk), 200, 0., 100.); + + hNphi[iDisk] = + tfs->make(Form("Nphi_%01i", iDisk), Form("Nphi %i", iDisk), 180, -180., 180.); + hNphiEcut[iDisk] = + tfs->make(Form("NphiEcut_%01i", iDisk), Form("NphiEcut %i", iDisk), 180, -180., 180.); + } + hNhit = tfs->make("Nhit", "Nhit", 2 * nCrystals, 0., 2 * nCrystals); + hNhitEcut = tfs->make("NhitEcut", "NhitEcut", 2 * nCrystals, 0., 2 * nCrystals); + +} // End of beginJob + +// =========================================================================== +// Event: +// - Select straight cosmic ray tracks crossing calo disks and fill time +// histogram for each readout channel +// - Filters good events +// =========================================================================== +bool caloT0alig::filter(art::Event& event) { + + bool retval(false); // preset to fail + + ++_nProcessed; + + if (_nProcessed % 10000 == 0 && _diagLevel > 0) + std::cout << "caloT0alig: Processing run/event " << event.run() << " " << event.id().event() + << std::endl; + + // + // *** TO BE IMPLEMENTED *** + // Check trigger bits to run only on trigger selected cosmics + // + + // Handle to the calorimeter + art::ServiceHandle geom; + if (!geom->hasElement()) { + mf::LogError("NO-CALO-GEOM") << "Calorimeter geometry not found!" << std::endl; + } + const Calorimeter& caloGeom = *(GeomHandle()); + + // Calorimeter crystal hits (average from readouts) + art::Handle CaloHitsHandle; + event.getByLabel(_caloHitTag, CaloHitsHandle); + const CaloHitCollection& CaloHits(*CaloHitsHandle); + + // Calorimeter clusters + art::Handle caloClustersHandle; + event.getByLabel(_caloClusterTag, caloClustersHandle); + const CaloClusterCollection& caloClusters(*caloClustersHandle); + + float sx[nDisks] = {}; + float sy[nDisks] = {}; + float sxy[nDisks] = {}; + float sy2[nDisks] = {}; + + int Ival[nDisks][nROchan]; + float Xval[nDisks][nROchan]; + float Yval[nDisks][nROchan]; + float Eval[nDisks][nROchan]; + float Tval[nDisks][nROchan]; + int nCry[nDisks] = {}; // Good cells on calo disks + int nChan[nDisks] = {}; // Good readout channels on calo disks + + float discr, bb, costh, Tres; + + for (unsigned int iClu = 0; iClu < caloClusters.size(); ++iClu) { // Loop on clusters + const CaloCluster& cluster = caloClusters.at(iClu); + + if (cluster.size() >= _cluHits) { // Hits in a cluster + if (_diagLevel > 2) + std::cout << "Nproc cluId " << _nProcessed << " " << iClu << std::endl; + int diskId = cluster.diskID(); + + std::vector cryList; + for (auto cryPtr : cluster.caloHitsPtrVector()) + cryList.push_back(std::distance(&CaloHits.at(0), cryPtr.get())); + + for (int iCry = 0; iCry < cluster.size(); iCry++) { // Crystals connected to cluster + + const CaloHit& hit = CaloHits.at(cryList[iCry]); + + int cryId = hit.crystalID(); + int cryDisk = caloGeom.crystal(cryId).diskID(); + if (cryDisk != diskId) + mf::LogError("WRONG-DISK") + << "Different clu/cry disk: " << diskId << " " << cryDisk << std::endl; + float cryE = hit.energyDep(); + + // Hit distribution vs phi and cristal index before t0alig selection cuts + + CLHEP::Hep3Vector crystalPos = caloGeom.mu2eToDiskFF( + diskId, caloGeom.crystal(cryId).position()); // Crystal position in disk FF frame + float xval = crystalPos.x(); + float yval = crystalPos.y(); + + float phiDeg = atan2(yval, xval) * 180. / 3.1416; + hNphi[cryDisk]->Fill(phiDeg); + hNhit->Fill(cryId); + if (cryE > _cryEmin) { + hNphiEcut[cryDisk]->Fill(phiDeg); + hNhitEcut->Fill(cryId); + } + + if (cryDisk == diskId && cryE > _cryEmin && cryE < 30.) { + + sx[diskId] = sx[diskId] + xval; + sy[diskId] = sy[diskId] + yval; + sxy[diskId] = sxy[diskId] + xval * yval; + sy2[diskId] = sy2[diskId] + pow(yval, 2); + + // Get values for the two SiPMs from digis + for (unsigned int iCha = 0; iCha < hit.recoCaloDigis().size(); iCha++) { + + int idx = CaloSiPMId(hit.recoCaloDigis().at(iCha)->SiPMID()).SiPMLocalId(); + Ival[diskId][nChan[diskId]] = hit.recoCaloDigis().at(idx)->SiPMID(); + Xval[diskId][nChan[diskId]] = xval; + Yval[diskId][nChan[diskId]] = yval; + Eval[diskId][nChan[diskId]] = cryE; + Tval[diskId][nChan[diskId]] = hit.recoCaloDigis().at(idx)->time(); + Tval[diskId][nChan[diskId]] = + Tval[diskId][nChan[diskId]] + Toff[Ival[diskId][nChan[diskId]]]; + if (_diagLevel > 2) + std::cout << "Disk Cry x y idx E T : " << cryDisk << " " << cryId << " " << xval + << " " << yval << " " << Ival[diskId][nChan[diskId]] << " " << cryE << " " + << Tval[diskId][nChan[diskId]] << std::endl; + nChan[diskId]++; + } + nCry[diskId]++; + } // Disk + energy cuts + } // Loop on crystals + } // Minimum number of cells in a cluster + } // Loop on clusters + + for (int iDisk = 0; iDisk < nDisks; iDisk++) { + + int Nfit = 0; + float enetot = 0.; + double t0 = 0.; + float Xfit[nCrystals] = {}, Yfit[nCrystals] = {}, Efit[nCrystals] = {}; + + if (nCry[iDisk] > _ncryCut) { // Enough good cells + + discr = float(nCry[iDisk]) * sy2[iDisk] - pow(sy[iDisk], 2); + + if (discr != 0.) { + + bb = (float(nCry[iDisk]) * sxy[iDisk] - sx[iDisk] * sy[iDisk]) / discr; + costh = 1. / sqrt(1. + bb * bb); + hcosTh[iDisk]->Fill(costh); + + if (costh > 0.2) { + hcosThCut[iDisk]->Fill(costh); + + for (int iCha = 0; iCha < nChan[iDisk]; iCha++) { + t0 = t0 + Eval[iDisk][iCha] * (Tval[iDisk][iCha] - Tcor[Ival[iDisk][iCha]] + + Yval[iDisk][iCha] / (cvel * costh)); + enetot = enetot + Eval[iDisk][iCha]; + hEcell[iDisk]->Fill(Eval[iDisk][iCha]); + + Xfit[Nfit] = Xval[iDisk][iCha]; + Yfit[Nfit] = Yval[iDisk][iCha]; + Efit[Nfit] = 10.; + if (_diagLevel > 3) + std::cout << Nfit << " " << Xfit[Nfit] << " " << Yfit[Nfit] << " " << Efit[Nfit] + << std::endl; + Nfit++; + } + + // Linear fit to x,y cell positions + TGraphErrors* emcpos = new TGraphErrors(Nfit, Xfit, Yfit, Efit, Efit); + TF1* linfit = new TF1("linfit", "[0]+[1]*x", -640., 640.); + emcpos->Fit("linfit", "rboq"); + double Par0 = linfit->GetParameter(0); + double Par1 = linfit->GetParameter(1); + double Chi2 = linfit->GetChisquare(); + if (_diagLevel > 3) + std::cout << Par0 << " " << Par1 << std::endl; + if (_diagLevel > 3) + std::cout << "chi2: " << Chi2 << std::endl; + Chi2 = Chi2 / linfit->GetNDF(); + hChi2[iDisk]->Fill(Chi2); + float costhFit = cos(atan(Par1)); + hcosFit[iDisk]->Fill(costhFit); + + if (Chi2 < 2. && enetot != 0.) { + + t0 = t0 / enetot; + hEmean[iDisk]->Fill(enetot / nChan[iDisk]); + + // Good event + retval = true; + ++_nFiltered; + + for (int iCha = 0; iCha < nChan[iDisk]; iCha++) { + float Ylinfit = Par0 + Par1 * Xval[iDisk][iCha]; + hYdif[iDisk]->Fill(Ylinfit - Yval[iDisk][iCha]); + Tres = Tval[iDisk][iCha] - Tcor[Ival[iDisk][iCha]] + + Yval[iDisk][iCha] / (cvel * costh) - t0; + hTres[Ival[iDisk][iCha]]->Fill(Tres); + + } // Loop on readout channels + } // Chi2 cut + } // costh cut + } // discr cut + } // If enough good cells + } // Loop on disks + + if (_diagLevel > 0) + std::cout << "caloT0alig: end of event " << std::endl; + + return retval; + +} // End of filter + +// =========================================================================== +// End job: +// - Writing of temporary output file to be used from next iteration +// - Last iteration: Writing of calibration file to be uploaded in Condition DB +// =========================================================================== +void caloT0alig::endJob() { + + if (_diagLevel > 0) + std::cout << "caloT0alig: Entering endJob" << std::endl; + + // Skip output writing when not in a minimization loop (i.e. filter only) + if (_iteration == "filt") { + if (_diagLevel > 0) + std::cout << "Iteration " << _iteration << ": Skip output writing" << std::endl; + return; + } + + // Temporary output file for next iteration + std::ofstream _fileTcor(_fileTcorName); + if (!_fileTcor.is_open()) { + throw cet::exception("caloT0alig") + << "ERROR! Cannot open output file " << _fileTcorName << std::endl; + } + + int Nevt; + float Tmea, Tsig, Chi2, Tval; + + for (int iCha = 0; iCha < nROchan; iCha++) { + TF1* gfit = new TF1("Gaussian", "gaus"); + hTres[iCha]->Fit(gfit, "q"); + Tmea = gfit->GetParameter(1); + Tsig = gfit->GetParameter(2); + Chi2 = gfit->GetChisquare(); + // TEMPORARY: MEAN/RMS instead of Gaussian fit params for MDC2020r CORSIKA production + // Tmea = hTres[iCha]->GetMean(); + // Tsig = hTres[iCha]->GetRMS(); + Nevt = hTres[iCha]->GetEntries(); + Tval = Tmea + Tcor[iCha]; + + _fileTcor << iCha << " " << Tval << " " << Tmea << " " << Tsig << " " << Chi2 << " " + << Nevt << std::endl; + } + _fileTcor.close(); + + // LAST ITERATION: + // Write output file for DB: rouId Tcor Tsigma Chi2 Nev + // Link this calibration to the previous time calibration step from laser/cosmics + // (unique DB identifier, table name + starting range, something else?) + // What if different run numbers have different step-0 corrections? + if (_iteration == "last") { + if (_diagLevel > 0) + std::cout << "Iteration " << _iteration << ": Writing output file for DB" << std::endl; + } + +} // End of endJob + +} // namespace mu2e + +DEFINE_ART_MODULE(mu2e::caloT0alig); diff --git a/CaloCalibration/NoiseAnalysis/CMakeLists.txt b/CaloCalibration/NoiseAnalysis/CMakeLists.txt new file mode 100755 index 0000000000..4ead09bb39 --- /dev/null +++ b/CaloCalibration/NoiseAnalysis/CMakeLists.txt @@ -0,0 +1,18 @@ +cet_build_plugin(BaselineAnalyzer art::module + REG_SOURCE src/BaselineAnalyzer_module.cc + LIBRARIES REG + art_root_io::TFileService_service + Offline::CaloConditions + Offline::CaloVisualizer + Offline::CaloVisualizer_dict + Offline::DAQ + Offline::DataProducts + Offline::ProditionsService + Offline::RecoDataProducts + ROOT::Hist + ROOT::Tree +) + +#install_headers(SUBDIRS inc) +install_source(SUBDIRS src) +install_fhicl(SUBDIRS fcl SUBDIRNAME Offline/CaloCalibration/NoiseAnalysis/fcl) diff --git a/CaloCalibration/NoiseAnalysis/fcl/analyzeBaselines.fcl b/CaloCalibration/NoiseAnalysis/fcl/analyzeBaselines.fcl new file mode 100644 index 0000000000..c05c3b8836 --- /dev/null +++ b/CaloCalibration/NoiseAnalysis/fcl/analyzeBaselines.fcl @@ -0,0 +1,69 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardServices.fcl" + +process_name : BaselineAnalyzer + +source : { + module_type : RootInput + fileNames : @nil + maxEvents : -1 + inputCommands: [ + "keep *", + "drop *_*STM*_*_*" + ] +} + +services : { + + #@local::Services.Reco + GeometryService : { + inputFile: "Offline/Mu2eG4/geom/geom_common.txt" + bFieldFile: "Offline/Mu2eG4/geom/bfgeom_reco_v01.txt" + simulatedDetector : {tool_type: "Mu2e"} + KinKalMaterial : @local::KinKalGeom.KKMaterial + } + GlobalConstantsService : {inputFile : "Offline/GlobalConstantsService/data/globalConstants_01.txt"} + DbService : @local::DbEmpty + ProditionsService : @local::Proditions + TimeTracker : {} +} + +physics : { + + filters : { + } + + producers : { + + CaloDigis : { + module_type : CaloDigisFromDTCEvents + dataType : 0 + diagLevel : 0 + useOfflineID : true + } + + } + + analyzers : { + + Baselines : { + module_type : BaselineAnalyzer + caloDigiTag : CaloDigis + verbosity : 0 + writeTXT : false + writeCSV : true + writePDF : true + thresholdOffset : 50 + thresholdOffsetPin : 100 + } + + } + + t1 : [ CaloDigis ] + e1 : [ Baselines ] + + trigger_paths : [t1] + end_paths : [e1] +} + +services.scheduler.wantSummary: true diff --git a/CaloCalibration/NoiseAnalysis/fcl/prolog.fcl b/CaloCalibration/NoiseAnalysis/fcl/prolog.fcl new file mode 100644 index 0000000000..d9db4343a4 --- /dev/null +++ b/CaloCalibration/NoiseAnalysis/fcl/prolog.fcl @@ -0,0 +1,26 @@ +#include "Offline/fcl/minimalMessageService.fcl" +#include "Offline/fcl/standardProducers.fcl" +#include "Offline/fcl/standardServices.fcl" +#include "Offline/CaloMC/fcl/prolog.fcl" +#include "Offline/CaloReco/fcl/prolog.fcl" +BEGIN_PROLOG + +CaloBaselineAna : { + module_type : BaselineAnalyzer + caloDigiTag : CaloDigis + verbosity : 0 + writeTXT : false + writeCSV : true + writePDF : true + thresholdOffset : 100 + thresholdOffsetPin : 100 +} + +NoiseAnalysis: { + analyzers: { + CaloBaselineAna : @local::CaloBaselineAna + } +} + + +END_PROLOG diff --git a/CaloCalibration/NoiseAnalysis/src/BaselineAnalyzer_module.cc b/CaloCalibration/NoiseAnalysis/src/BaselineAnalyzer_module.cc new file mode 100644 index 0000000000..6eeabbfece --- /dev/null +++ b/CaloCalibration/NoiseAnalysis/src/BaselineAnalyzer_module.cc @@ -0,0 +1,582 @@ +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" + +#include "artdaq-core/Data/Fragment.hh" +#include + +#include "Offline/CaloVisualizer/inc/THMu2eCaloDisk.hh" +#include "Offline/DAQ/inc/CaloDAQUtilities.hh" +#include "Offline/DataProducts/inc/CaloConst.hh" +#include "Offline/RecoDataProducts/inc/CaloDigi.hh" +#include "artdaq-core-mu2e/Data/EventHeader.hh" +#include "artdaq-core-mu2e/Overlays/DTCEventFragment.hh" +#include "artdaq-core-mu2e/Overlays/Decoders/CalorimeterDataDecoder.hh" +#include "artdaq-core-mu2e/Overlays/FragmentType.hh" + +#include "cetlib_except/exception.h" + +//-- insert calls to proditions ..for calodmap----- +#include "Offline/CaloConditions/inc/CalDAQMap.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +//------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TCanvas.h" +#include "TEllipse.h" +#include "TF1.h" +#include "TFile.h" +#include "TFitResult.h" +#include "TFitResultPtr.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TH1.h" +#include "TH2.h" +#include "TLatex.h" +#include "TLegend.h" +#include "TLine.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TText.h" +#include "TTree.h" +#include "art_root_io/TFileService.h" +#include + +namespace mu2e { +class BaselineAnalyzer : public art::EDAnalyzer { +public: + // clang-format off + struct Config { + fhicl::Atom caloDigiTag {fhicl::Name("caloDigiTag" ) , fhicl::Comment("caloDigiTag"), ""}; + fhicl::Atom verbosity {fhicl::Name("verbosity" ) , fhicl::Comment("Verbosity [0-2]"), 1}; + fhicl::Atom writeTXT {fhicl::Name("writeTXT" ) , fhicl::Comment("Write per-board text files with thresholds"), false}; + fhicl::Atom TXTfoldername {fhicl::Name("TXTfoldername" ) , fhicl::Comment("Folder to write thresholds into"), ""}; + fhicl::Atom writeCSV {fhicl::Name("writeCSV" ) , fhicl::Comment("Write CSV file with thresholds"), false}; + fhicl::Atom CSVfilename {fhicl::Name("CSVfilename" ) , fhicl::Comment("CSV file to write thresholds into"), ""}; + fhicl::Atom writePDF {fhicl::Name("writePDF" ) , fhicl::Comment("Write PDF report"), false}; + fhicl::Atom PDFfilename {fhicl::Name("PDFfilename" ) , fhicl::Comment("PDF report file name"), ""}; + fhicl::Atom thresholdOffset {fhicl::Name("thresholdOffset" ) , fhicl::Comment("Offset with respect to gaussian mean"), 100}; + fhicl::Atom thresholdOffsetPin {fhicl::Name("thresholdOffsetPin" ) , fhicl::Comment("Offset with respect to gaussian mean (pin diodes)"), 50}; + fhicl::Atom hotStdDev {fhicl::Name("hotStdDev" ) , fhicl::Comment("StdDev limit for hot channels"), 6}; + fhicl::Atom coldStdDev {fhicl::Name("coldStdDev" ) , fhicl::Comment("StdDev limit for cold channels"), 3}; + }; + // clang-format on + + explicit BaselineAnalyzer(const art::EDAnalyzer::Table& config); + void beginRun(art::Run const& run) override; + void analyze(art::Event const& event) override; + void endJob() override; + void FitHistograms(); + void WriteReport(); + +private: + std::string caloDigiTag_; + int verbosity_; + bool writeTXT_; + std::string TXTfoldername_; + bool writeCSV_; + std::string CSVfilename_; + bool writePDF_; + std::string PDFfilename_; + int thresholdOffset_; + int thresholdOffsetPin_; + double hotStdDev_; + double coldStdDev_; + + mu2e::ProditionsHandle _calodaqconds_h; + + double xmin; + double xmax; + bool titlesSet; + + std::map> channelType; + std::map> h1_baseline_map; + + TH2D* h2_baselines; + TH1D* h1_means; + TH1D* h1_sigmas; + TH1D* h1_threshold; + TGraphErrors* g_baselines; + TGraph* g_fitsigmas; + TGraph* g_stddev; + TGraph* g_stddev_pin; + TGraph* g_stddev_lyso; + TGraph* g_stddev_empty; + TGraph* g_thresholds; + mu2e::THMu2eCaloDisk* h2_disk0; + mu2e::THMu2eCaloDisk* h2_disk1; +}; +} // namespace mu2e + +mu2e::BaselineAnalyzer::BaselineAnalyzer(const art::EDAnalyzer::Table& config) : + art::EDAnalyzer{config}, caloDigiTag_(config().caloDigiTag()), verbosity_(config().verbosity()), + writeTXT_(config().writeTXT()), TXTfoldername_(config().TXTfoldername()), + writeCSV_(config().writeCSV()), CSVfilename_(config().CSVfilename()), + writePDF_(config().writePDF()), PDFfilename_(config().PDFfilename()), + thresholdOffset_(config().thresholdOffset()), + thresholdOffsetPin_(config().thresholdOffsetPin()), hotStdDev_(config().hotStdDev()), + coldStdDev_(config().coldStdDev()) { + art::ServiceHandle tfs; + + xmin = 2048 - 150; + xmax = 2048 + 300; + h2_baselines = tfs->make("h2_baselines", "Baselines;BoardID*100 + ChannelID;Baseline [ADC]", + 16100, 0, 16100, xmax - xmin, xmin, xmax); + h1_means = tfs->make("h1_means", "All channels baselines;ADC", int(0.5 * (xmax - xmin)), + xmin, xmax); + h1_sigmas = tfs->make("h1_sigmas", "All channels baseline sigmas;ADC", 100, 0, 10); + h1_threshold = tfs->make("h1_threshold", "All channel thresholds;ADC", + int(0.5 * (xmax - xmin)), xmin, xmax); + g_baselines = tfs->makeAndRegister( + "g_baselines", "Baselines;BoardID*100 + ChannelID;Baseline [ADC]"); + g_baselines->SetMarkerStyle(20); + g_fitsigmas = tfs->makeAndRegister("g_fitsigmas", + "Gaussian fit sigmas;BoardID*100 + ChannelID;ADC"); + g_fitsigmas->SetMarkerStyle(20); + g_stddev = tfs->makeAndRegister("g_stddev", "Channel StdDev;BoardID*100 + ChannelID;ADC"); + g_stddev->SetMarkerStyle(20); + g_stddev_pin = + tfs->makeAndRegister("g_stddev_pin", "[pin-diode];BoardID*100 + ChannelID;ADC"); + g_stddev_pin->SetMarkerStyle(20); + g_stddev_pin->SetMarkerColor(kRed); + g_stddev_lyso = + tfs->makeAndRegister("g_stddev_lyso", "[LYSO];BoardID*100 + ChannelID;ADC"); + g_stddev_lyso->SetMarkerStyle(20); + g_stddev_lyso->SetMarkerColor(kGreen); + g_stddev_empty = + tfs->makeAndRegister("g_stddev_empty", "[empty];BoardID*100 + ChannelID;ADC"); + g_stddev_empty->SetMarkerStyle(20); + g_stddev_empty->SetMarkerColor(kGray); + g_thresholds = + tfs->makeAndRegister("g_thresholds", "Thresholds;BoardID*100 + ChannelID;ADC"); + g_thresholds->SetMarkerStyle(20); + g_thresholds->SetMarkerColor(kRed); + h2_disk0 = + tfs->makeAndRegister("h2_disk0", "Disk 0", "h2_disk0", "Disk 0", 0); + h2_disk1 = + tfs->makeAndRegister("h2_disk1", "Disk 1", "h2_disk1", "Disk 1", 1); + h2_disk0->SetCombineMode(mu2e::ECombineMode::kAverage); + h2_disk1->SetCombineMode(mu2e::ECombineMode::kAverage); + h2_disk0->SetDrawOption("COLZL"); + h2_disk1->SetDrawOption("COLZL"); + + auto baseBoardDir = tfs->mkdir("boards"); + for (int boardID = 0; boardID < CaloConst::_nDIRAC; boardID++) { + TString dirname(Form("Board%03d", boardID)); + auto boardDir = baseBoardDir.mkdir(dirname.Data()); + for (int chanID = 0; chanID < CaloConst::_nChPerDIRAC; chanID++) { + TString hname = Form("h1_baseline_b%03d_c%02d", boardID, chanID); + TString htitle = Form("Baseline of Board %03d channel %02d", boardID, chanID); + h1_baseline_map[boardID][chanID] = + boardDir.make(hname, htitle, xmax - xmin, xmin, xmax); + } + } + titlesSet = false; +} + +void mu2e::BaselineAnalyzer::beginRun(art::Run const& run) { + if (CSVfilename_ == "") { + CSVfilename_ = std::string("thresholds_run") + run.id().run() + "_offset" + thresholdOffset_ + + "_offsetPin" + thresholdOffsetPin_ + ".csv"; + } + if (PDFfilename_ == "") { + PDFfilename_ = std::string("thresholds_run") + run.id().run() + "_offset" + thresholdOffset_ + + "_offsetPin" + thresholdOffsetPin_ + ".pdf"; + } +} + +void mu2e::BaselineAnalyzer::analyze(art::Event const& event) { + const auto& caloDigis = *event.getValidHandle(consumes(caloDigiTag_)); + art::ServiceHandle tfs; + + mu2e::CalDAQMap const& calodaqconds = _calodaqconds_h.get(event.id()); + if (!titlesSet) { // Need to do this here in order to have calodaqconds + for (int sipmid = 0; sipmid < CaloConst::_nChannel; sipmid++) { + mu2e::CaloSiPMId SiPMID_(sipmid); + if (!SiPMID_.isValid()) + continue; + mu2e::CaloRawSiPMId rawId = calodaqconds.rawId(SiPMID_); + int boardID = rawId.dirac(); + int chanID = rawId.ROCchannel(); + channelType[boardID][chanID] = SiPMID_.detType(); + TString typeName = mu2e::CaloConst::detTypeName(channelType[boardID][chanID]); + TString htitle = + Form("Baseline of Board %03d channel %02d [%s]", boardID, chanID, typeName.Data()); + h1_baseline_map[boardID][chanID]->SetTitle(htitle); + } + titlesSet = true; + } + + // Loop over the calo digis of this event + for (uint ihit = 0; ihit < caloDigis.size(); ihit++) { + int SiPMID = caloDigis[ihit].SiPMID(); + std::vector waveform = caloDigis[ihit].waveform(); + + mu2e::CaloSiPMId SiPMID_(SiPMID); + if (!SiPMID_.isValid()) + continue; + mu2e::CaloRawSiPMId rawId = calodaqconds.rawId(SiPMID_); + int boardID = rawId.dirac(); + int chanID = rawId.ROCchannel(); + + // Fill hist + for (auto sample : waveform) { + h2_baselines->Fill(boardID * 100 + chanID, sample); + h1_baseline_map[boardID][chanID]->Fill(sample); + } + } +} + +void mu2e::BaselineAnalyzer::endJob() { + // Remove empty hists + for (int boardID = 0; boardID < CaloConst::_nDIRAC; boardID++) { + for (int chanID = 0; chanID < CaloConst::_nChPerDIRAC; chanID++) { + if (h1_baseline_map[boardID][chanID]->GetEntries() == 0) { + delete h1_baseline_map[boardID][chanID]; + h1_baseline_map[boardID].erase(chanID); + } + } + if (h1_baseline_map[boardID].empty()) { + h1_baseline_map.erase(boardID); + } + } + + FitHistograms(); + if (writePDF_) { + WriteReport(); + } +} + +void mu2e::BaselineAnalyzer::FitHistograms() { + // Perform gaussian fits + int failed_fits = 0; + int unprecise_fits = 0; + std::set> failed_map; + std::set> unprecise_map; + std::ofstream outputCSV; + if (writeCSV_) { + outputCSV.open(CSVfilename_); + if (!outputCSV.is_open()) { + std::cout << "Warning! Can't open file " << CSVfilename_ << "\n"; + writeCSV_ = false; + } + } + std::set all_boards; + std::map> all_baselines; + std::map> all_thresholds; + for (auto board_pair : h1_baseline_map) { + int board = board_pair.first; + all_boards.insert(board); + std::ofstream outputBaselineFile; + TString fname = Form("%s/dirac%03d.baseline", TXTfoldername_.c_str(), board); + + if (writeTXT_) { + outputBaselineFile.open(fname); + if (!outputBaselineFile.is_open()) { + std::cout << "Warning! Can't open file " << fname << "\n"; + writeTXT_ = false; + } + } + + for (auto channel_pair : board_pair.second) { + int channel = channel_pair.first; + + bool pindiode = (channelType[board][channel] == mu2e::CaloConst::detType::PINDiode); + + h1_baseline_map[board][channel]->Scale(1. / h1_baseline_map[board][channel]->Integral()); + + double stddev = h1_baseline_map[board][channel]->GetRMS(); + if (verbosity_ > 0) { + if (stddev > hotStdDev_) { + std::cout << "Hot channel: Board " << board << " Channel " << channel << " (" + << std::setprecision(4) << std::setw(5) << stddev << ") [" + << channelType[board][channel] << "]\n"; + } + if (stddev < coldStdDev_) { + std::cout << "Cold channel: Board " << board << " Channel " << channel << " (" + << std::setprecision(4) << std::setw(5) << stddev << ") [" + << channelType[board][channel] << "]\n"; + } + } + + switch (channelType[board][channel]) { + case mu2e::CaloConst::detType::CsI: + g_stddev->AddPoint(board * 100 + channel, stddev); + break; + case mu2e::CaloConst::detType::PINDiode: + g_stddev_pin->AddPoint(board * 100 + channel, stddev); + break; + case mu2e::CaloConst::detType::CAPHRI: + g_stddev_lyso->AddPoint(board * 100 + channel, stddev); + break; + case mu2e::CaloConst::detType::Invalid: + g_stddev_empty->AddPoint(board * 100 + channel, stddev); + break; + default: + std::cout << "Board " << board << " Channel " << channel << " : unknown type " + << channelType[board][channel] << "\n"; + } + + if (verbosity_ > 1) + std::cout << "Fitting board " << board << " channel " << channel << " ... "; + double guess_max = h1_baseline_map[board][channel]->GetBinContent( + h1_baseline_map[board][channel]->GetMaximumBin()); + double guess_mean = h1_baseline_map[board][channel]->GetMean(); + double guess_sigma = h1_baseline_map[board][channel]->GetRMS(); + + ((TF1*)gROOT->GetFunction("gaus"))->SetParameters(guess_max, guess_mean, guess_sigma); + TFitResultPtr fit_result = + h1_baseline_map[board][channel]->Fit("gaus", "SQW", "", xmin, xmax); + h1_baseline_map[board][channel]->GetFunction("gaus")->SetNpx(1000); + if (fit_result < 0) { + if (verbosity_ > 1) + std::cout << " <---- BAD FIT!!! Status " << fit_result << "\n"; + failed_fits++; + failed_map.insert(std::make_pair(board, channel)); + continue; + } + double fit_mean = fit_result->Parameter(1); + double fit_sigma = fit_result->Parameter(2); + double fit_chi2 = fit_result->Chi2(); + double fit_ndf = fit_result->Ndf(); + int threshold = int(round(fit_mean + thresholdOffset_)); + + // Pin diodes + if (pindiode) { + threshold = int(round(fit_mean + thresholdOffsetPin_)); + } + + if (verbosity_ > 0) + std::cout << fit_mean << " +- " << fit_sigma << " chi2/ndf=" << fit_chi2 << "/" << fit_ndf; + if (fit_result > 0) { + if (verbosity_ > 0) + std::cout << " <---- UNPRECISE! Status " << fit_result; + unprecise_fits++; + unprecise_map.insert(std::make_pair(board, channel)); + } + if (verbosity_ > 0) + std::cout << "\n"; + + h1_means->Fill(fit_mean); + h1_sigmas->Fill(fit_sigma); + h1_threshold->Fill(threshold); + g_baselines->SetPoint(g_baselines->GetN(), board * 100 + channel, fit_mean); + g_baselines->SetPointError(g_baselines->GetN() - 1, 0, fit_sigma); + g_fitsigmas->SetPoint(g_fitsigmas->GetN(), board * 100 + channel, fit_sigma); + g_thresholds->SetPoint(g_thresholds->GetN(), board * 100 + channel, threshold); + if (board < 80) { + h2_disk0->FillRaw(board, channel, stddev); + } else { + h2_disk1->FillRaw(board, channel, stddev); + } + + all_baselines[board][channel] = fit_mean; + all_thresholds[board][channel] = threshold; + + TString lineout = Form("%d\t%.2f\t%.2f\t%d", channel, fit_mean, fit_sigma, threshold); + if (writeTXT_) { + outputBaselineFile << lineout << "\n"; + } + } + if (writeTXT_) { + outputBaselineFile.close(); + } + } + + // Check for missing channels -> assign default values + for (auto b : all_boards) { + for (int c = 0; c < 20; c++) { + if (all_baselines[b].find(c) == all_baselines[b].end()) { + if (verbosity_ > 0) { + std::cout << "Setting default values for empty Board " << b << " Channel " << c << "\n"; + } + all_baselines[b][c] = 2048; + all_thresholds[b][c] = all_baselines[b][c] + thresholdOffset_; + } + } + } + + if (writeCSV_) { + for (auto b : all_boards) { + std::stringstream ss; + ss << b << ",[["; + for (int c = 0; c < 20; c++) { + if (c) + ss << ","; + ss << all_baselines[b][c]; + } + ss << "]],[["; + for (int c = 0; c < 20; c++) { + if (c) + ss << ","; + ss << all_thresholds[b][c]; + } + ss << "]],\"\",\"\",\"\""; + std::string csv_line = ss.str(); + outputCSV << csv_line << "\n"; + } + outputCSV.close(); + std::cout << "wrote file " << CSVfilename_ << "\n"; + } + + if (verbosity_ > 0) { + std::cout << "Failed fits: " << failed_fits << "\n"; + for (auto pair : failed_map) { + std::cout << "Board " << pair.first << " Channel " << pair.second << "\n"; + } + + std::cout << "Unprecise fits: " << unprecise_fits << "\n"; + for (auto pair : unprecise_map) { + std::cout << "Board " << pair.first << " Channel " << pair.second << "\n"; + } + } +} + +void mu2e::BaselineAnalyzer::WriteReport() { + TCanvas* can = new TCanvas("pdf_canvas", "Baseline Report", 2000, 1400); + can->SetBatch(kTRUE); + + TString pdfname(PDFfilename_.c_str()); + can->SaveAs(pdfname + "["); + + // Summary pages + can->Clear(); + can->Divide(2, 2); + can->cd(1); + h2_baselines->Draw("COLZ"); + gPad->SetLogz(); + can->cd(2); + h1_means->Draw(); + can->cd(3); + h1_sigmas->Draw(); + can->cd(4); + h1_threshold->Draw(); + can->SaveAs(pdfname); + + TText* text0 = new TText(0, 50, "DISK 0"); + TText* text1 = new TText(0, 50, "DISK 1"); + TText* textL = new TText(0, -50, "LEFT"); + TText* textR = new TText(0, -50, "RIGHT"); + text0->SetTextAlign(22); + text1->SetTextAlign(22); + textL->SetTextAlign(22); + textR->SetTextAlign(22); + + h2_disk0->GetZaxis()->SetRangeUser(coldStdDev_, hotStdDev_); + h2_disk1->GetZaxis()->SetRangeUser(coldStdDev_, hotStdDev_); + + gStyle->SetPalette(kViridis); + can->Clear(); + can->Divide(2, 2); + can->cd(1); + h2_disk0->SetCombineMode(mu2e::ECombineMode::kLeft); + h2_disk0->Draw("colz l"); + text0->Draw("same"); + textL->Draw("same"); + can->cd(2); + h2_disk0->SetCombineMode(mu2e::ECombineMode::kRight); + h2_disk0->Draw("colz l"); + text0->Draw("same"); + textR->Draw("same"); + can->cd(3); + h2_disk1->SetCombineMode(mu2e::ECombineMode::kLeft); + h2_disk1->Draw("colz l"); + text1->Draw("same"); + textL->Draw("same"); + can->cd(4); + h2_disk1->SetCombineMode(mu2e::ECombineMode::kRight); + h2_disk1->Draw("colz l"); + text1->Draw("same"); + textR->Draw("same"); + can->SaveAs(pdfname); + + // Stddev + can->Clear(); + g_stddev->GetXaxis()->SetLimits(0, 16100); + double xmin = g_stddev->GetXaxis()->GetXmin(); + double xmax = g_stddev->GetXaxis()->GetXmax(); + double ymax = g_stddev->GetYaxis()->GetXmax(); + + if (g_stddev_pin->GetYaxis()->GetXmax() > ymax) { + ymax = g_stddev_pin->GetYaxis()->GetXmax(); + } + if (g_stddev_lyso->GetYaxis()->GetXmax() > ymax) { + ymax = g_stddev_lyso->GetYaxis()->GetXmax(); + } + if (hotStdDev_ > ymax) { + ymax = hotStdDev_; + } + g_stddev->GetYaxis()->SetRangeUser(0, ymax + 0.1); + + g_stddev->Draw("AP"); + if (g_stddev_pin->GetN() > 0) + g_stddev_pin->Draw("P"); + if (g_stddev_lyso->GetN() > 0) + g_stddev_lyso->Draw("P"); + if (g_stddev_empty->GetN() > 0) + g_stddev_empty->Draw("P"); + TLegend* leg = new TLegend(0.65, 0.91, 0.9, 0.98); + leg->AddEntry(g_stddev_pin, g_stddev_pin->GetTitle(), "P"); + leg->AddEntry(g_stddev_lyso, g_stddev_lyso->GetTitle(), "P"); + leg->AddEntry(g_stddev_empty, g_stddev_empty->GetTitle(), "P"); + leg->Draw(); + + TLine* l_cold = new TLine(xmin, coldStdDev_, xmax, coldStdDev_); + TLine* l_hot = new TLine(xmin, hotStdDev_, xmax, hotStdDev_); + l_cold->SetLineColor(kBlue); + l_hot->SetLineColor(kRed); + l_cold->SetLineWidth(2); + l_hot->SetLineWidth(2); + l_cold->SetLineStyle(kDashed); + l_hot->SetLineStyle(kDashed); + l_cold->Draw("SAME"); + l_hot->Draw("SAME"); + TLatex* text_cold = new TLatex(xmax, coldStdDev_, "#downarrow COLD"); + TLatex* text_hot = new TLatex(xmax, hotStdDev_, "#uparrow HOT"); + text_cold->SetTextSize(0.04); + text_hot->SetTextSize(0.04); + text_cold->SetTextAlign(12); + text_hot->SetTextAlign(12); + text_cold->SetTextColor(kBlue); + text_hot->SetTextColor(kRed); + text_cold->Draw("SAME"); + text_hot->Draw("SAME"); + can->SaveAs(pdfname); + delete leg; + + // Per-board channel histograms, grouped into pages of 20 plots each + for (const auto& board_pair : h1_baseline_map) { + int board = board_pair.first; + can->Clear(); + can->Divide(5, 4); + for (int channel = 0; channel < 20; channel++) { + can->cd(channel + 1); + if (h1_baseline_map[board].find(channel) != h1_baseline_map[board].end()) { + h1_baseline_map[board][channel]->GetXaxis()->SetRangeUser(1950, 2150); + h1_baseline_map[board][channel]->Draw("HIST"); + gPad->SetLogy(); + } + } + can->SaveAs(pdfname); + } + + can->SaveAs(pdfname + "]"); + delete can; +} + +DEFINE_ART_MODULE(mu2e::BaselineAnalyzer) diff --git a/CaloCalibration/NoiseAnalysis/src/SConscript b/CaloCalibration/NoiseAnalysis/src/SConscript new file mode 100644 index 0000000000..edd24cdd17 --- /dev/null +++ b/CaloCalibration/NoiseAnalysis/src/SConscript @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import os, re +Import('env') + +Import('mu2e_helper') + +helper=mu2e_helper(env) + +rootlibs = env['ROOTLIBS'] + +extrarootlibs = [ 'Minuit' ] + +mainlib = helper.make_mainlib ( [] ) + +helper.make_plugins( [ mainlib, + 'mu2e_Mu2eUtilities', + 'mu2e_CaloCluster', + 'mu2e_DAQ', + 'mu2e_ParticleID', + 'mu2e_RecoDataProducts', + 'mu2e_TrackerConditions', + 'mu2e_GeometryService', + 'mu2e_CalorimeterGeom', + 'mu2e_CaloVisualizer', + 'mu2e_CaloConditions', + 'mu2e_DAQConditions', + 'mu2e_DbService', + 'mu2e_DbTables', + 'mu2e_CosmicRayShieldGeom', + 'mu2e_ExtinctionMonitorFNAL_Geometry', + 'mu2e_ProtonBeamDumpGeom', + 'mu2e_StoppingTargetGeom', + 'mu2e_TrackerGeom', + 'mu2e_GlobalConstantsService_GlobalConstantsService_service', + 'mu2e_GlobalConstantsService', + 'mu2e_ProditionsService', + 'mu2e_ConfigTools', + 'mu2e_DataProducts', + 'mu2e_GeneralUtilities', + 'art_Framework_Core', + 'art_Framework_Principal', + 'art_Framework_Services_Registry', + 'art_root_io_tfile_support', + 'art_root_io_TFileService', + 'art_Framework_Services_Optional_RandomNumberGenerator_service', + 'art_Persistency_Common', + 'art_Persistency_Provenance', + 'art_Utilities', + 'canvas', + 'MF_MessageLogger', + 'fhiclcpp', + 'fhiclcpp_types', + 'tbb', + 'cetlib', + 'cetlib_except', + 'CLHEP', + rootlibs, + extrarootlibs, + 'xerces-c', + 'boost_filesystem', + 'artdaq-core-mu2e_Overlays', + 'artdaq-core_Data', + 'artdaq-core-mu2e_Data', + 'artdaq-core-mu2e_Data_dict' + ] ) + +# this tells emacs to view this file in python mode. +# Local Variables: +# mode:python +# End: diff --git a/DAQ/fcl/prolog.fcl b/DAQ/fcl/prolog.fcl index 35e0c1395d..1a176b1e72 100644 --- a/DAQ/fcl/prolog.fcl +++ b/DAQ/fcl/prolog.fcl @@ -15,6 +15,13 @@ DAQ : { allowOfflineFallbackWhenPanelMapMissing : false forceOfflineAddressing : false } + + CaloDigis : { + module_type : CaloDigisFromDTCEvents + dataType : 0 + diagLevel : 0 + useOfflineID : true + } } filters : {