21 #include "TFitResult.h" 22 #include "TFitResultPtr.h" 44 #include "Utilities/func/MathUtil.h" 45 #include "NovaDAQConventions/DAQConventions.h" 122 gStyle->SetOptStat(0);
187 it->second += *profsmap;
217 it->second += *profsmap;
231 TFile*
outfile =
new TFile(
"Atten_Fit_Results.root",
"RECREATE");
239 const TString mcStr =
fIsMC ?
"_mc" :
"";
240 FILE* fConstsCSV = fopen(
"calib_atten_consts"+mcStr+
".csv",
"w");
241 FILE* fPointsCSV = fopen(
"calib_atten_points"+mcStr+
".csv",
"w");
243 fprintf(fConstsCSV,
"#coeff_exp,atten_length,background,edge_low,edge_high,coeff_low,coeff_high,chisq,channel,tv\n");
244 fprintf(fPointsCSV,
"#w,factor,channel,tv\n");
248 <<
plane <<
"/" << plane_max;
262 <<
cell <<
"/" << cell_max;
292 if(prof->GetXaxis()->GetXmin() == 0 ||
293 prof->GetXaxis()->GetXmax() == 0 ||
309 if(prof->GetXaxis()->GetXmin() == 0 ||
310 prof->GetXaxis()->GetXmax() == 0 ||
338 if(!dir) dir = outfile->mkdir(Form(
"plane_%03d",
plane));
340 prof->Write(Form(
"Atten_Fit_%03d_%03d",
plane,
cell));
371 int isMuonCatcher=firstPlane <=
plane;
372 view = geom->
Plane(plane)->
View()+(2*isMuonCatcher);
385 bool isConsolidated =
true;
388 const std::vector<geo::OfflineChan> chans =
fChannelMapProf.find(0)->second.Channels();
390 if((
chan.Plane() > 3 && isND) || (
chan.Plane() > 1 && isND==
false)){
391 isConsolidated =
false;
398 mf::LogInfo(
"AttenFit") <<
"Input file appears to already be consolidated by view. Skipping that step.";
403 mf::LogInfo(
"AttenFit") <<
"Input file is not consolidated by view. Consolidating";
412 const std::vector<geo::OfflineChan> chans = profmap.
Channels();
439 fit->GetRange(xmin, xmax);
443 for(
int i = 0;
i < prof->GetNbinsX()+2; ++
i){
444 const double x = prof->GetXaxis()->GetBinCenter(
i);
446 if(x < xmin || x > xmax) {
451 diff +=
util::sqr((fit->Eval(x)-prof->GetBinContent(
i))/fit->Eval(x));
454 const double ret =
sqrt(diff/tot);
475 TProfile*
p = h->ProfileX();
477 TAxis* xax = h->GetXaxis();
478 TH1*
ret =
new TH1D(
"",
"", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
479 ret->GetXaxis()->SetTitle(xax->GetTitle());
480 for(
int ix = 0; ix < xax->GetNbins()+2; ++ix){
481 ret->SetBinContent(ix, p->GetBinContent(ix));
482 ret->SetBinError(ix, p->GetBinError(ix));
491 TAxis* xax = h->GetXaxis();
492 TAxis* yax = h->GetYaxis();
493 TH1*
ret =
new TH1D(
"",
"", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
494 ret->GetXaxis()->SetTitle(xax->GetTitle());
495 for(
int ix = 0; ix < xax->GetNbins()+2; ++ix){
498 for(
int iy = 0; iy < yax->GetNbins()+2; ++iy){
499 tot += h->GetBinContent(ix, iy);
502 for(
int iy = 0; iy < yax->GetNbins()+2; ++iy){
503 seen += h->GetBinContent(ix, iy);
505 const double z = h->GetBinContent(ix, iy);
507 const double frac = 1-(seen-tot/2)/z;
508 ret->SetBinContent(ix, yax->GetBinLowEdge(iy)+frac*yax->GetBinWidth(iy));
510 ret->SetBinError(ix,
sqrt(tot)*yax->GetBinWidth(iy)/(2*
z));
521 TAxis* xax = h->GetXaxis();
522 TAxis* yax = h->GetYaxis();
523 TH1*
ret =
new TH1D(
"",
"", xax->GetNbins(), xax->GetXmin(), xax->GetXmax());
524 ret->GetXaxis()->SetTitle(xax->GetTitle());
525 for(
int ix = 0; ix < xax->GetNbins()+2; ++ix){
528 for(
int iy = 0; iy < yax->GetNbins()+2; ++iy){
529 tot += h->GetBinContent(ix, iy);
536 for(
int iy = 0; iy < yax->GetNbins()+2; ++iy){
537 const double z = h->GetBinContent(ix, iy);
543 double weight_lo = (seen-(.25*
tot))/
z;
544 double weight_hi = 1-(seen-(.75*
tot))/
z;
545 if(weight_lo > 1) weight_lo = 1;
546 if(weight_hi > 1) weight_hi = 1;
547 double weight = weight_lo+weight_hi-1;
548 if(weight < 0) weight = 0;
551 const double c = yax->GetBinCenter(iy);
566 ret->SetBinContent(ix, avg);
567 ret->SetBinError(ix,
sqrt(rms)/
sqrt(tot/2));
570 ret->SetBinContent(ix, 0);
571 ret->SetBinError(ix, 1e10);
582 for(
int i = 0;
i < prof->GetNbinsX()+2; ++
i){
583 const double w = prof->GetXaxis()->GetBinCenter(
i);
584 const double y = prof->GetBinContent(
i);
587 prof->SetBinContent(
i, y/corr);
588 prof->SetBinError(
i, prof->GetBinError(
i)/
corr);
602 view = (chan.
Plane());
620 double operator()(
double* xs,
double*
ps)
622 fRes->coeff_exp = ps[0];
623 fRes->atten_length = ps[1];
624 fRes->background = ps[2];
625 return fRes->MeanPEPerCmAt(xs[0]);
629 ExpFitAdaptor* adapt =
new ExpFitAdaptor(res);
631 TString fitStrArr[4]={
"fitX",
"fitY",
"fitMuX",
"fitMuY"};
632 TString fit_name = fitStrArr[
view];
633 double rangeup=
range;
638 TF1* fit =
new TF1(fit_name, adapt, -range, +rangeup, 3);
640 fit->SetParameter(0, 11);
642 fit->SetParameter(0, 43);
643 fit->SetParLimits(0, 0, 1000);
644 fit->SetParLimits(1, 0, 1000);
645 fit->SetParLimits(2, 0, 1000);
648 fit->SetParameter(0, 70);
650 fit->SetParLimits(0, 0, 1000);
651 fit->SetParLimits(1, 0, 1000);
652 fit->SetParLimits(2, 0, 1000);
654 fit->SetParameter(1, 400);
655 fit->SetParameter(2, 7.5);
656 fit->SetParName(0,
"norm");
657 fit->SetParName(1,
"atten");
658 fit->SetParName(2,
"bkgd");
660 fit->SetRange(-range, +rangeup);
663 TFitResultPtr fr = prof->Fit(fit,
"RS"+opt);
670 fit->SetParameter(0, 43);
671 fit->SetParameter(1, 452);
672 fit->SetParameter(2, 7.6);
698 struct RolloffFitAdaptor
701 double operator()(
double* xs,
double*
ps)
703 fRes->coeff_low = ps[0];
704 fRes->coeff_high = ps[1];
705 fRes->edge_low = ps[2];
706 fRes->edge_high = ps[3];
707 return fRes->MeanPEPerCmAt(xs[0]);
711 RolloffFitAdaptor* adapt =
new RolloffFitAdaptor(res);
720 double rangeup=
range;
726 TF1* fit =
new TF1(view ==
geo::kX ?
"rollFitX" :
"rollFitY",
727 adapt, -range, +rangeup, 4);
728 fit->SetParName(0,
"coeffL");
729 fit->SetParName(1,
"coeffR");
730 fit->SetParName(2,
"edgeL");
731 fit->SetParName(3,
"edgeR");
732 fit->SetParameter(0, 0);
733 fit->SetParameter(1, 0);
734 fit->SetParameter(2, -edge);
735 fit->SetParameter(3, edgeup);
736 fit->SetRange(-range, +rangeup);
740 TFitResultPtr fr = prof->Fit(fit,
"RS"+opt);
744 res->
edge_low = fit->GetParameter(2);
761 TGraph*
ret =
new TGraph;
764 const int kNumControlPts = 20;
770 const double kRange = 1.5*(x1-x0)/kNumControlPts;
772 for(
int n = 0;
n < kNumControlPts; ++
n){
773 const double x = x0+double(
n)/(kNumControlPts-1)*(x1-x0);
775 std::vector<double> xs, ys, ws;
777 for(
int i = 0;
i < prof->GetNbinsX()+2; ++
i){
778 const double xi = prof->GetBinCenter(
i);
779 if(xi < x0 || xi > x1)
continue;
780 if(
fabs(xi-x) > kRange)
continue;
783 ys.push_back(prof->GetBinContent(
i));
792 const double y = m*x+
c;
794 ret->SetPoint(
n, x, y);
798 ret->SetLineWidth(2);
799 ret->SetMarkerColor(
kRed);
800 ret->SetLineColor(
kRed);
801 ret->SetMarkerStyle(kFullCircle);
809 TGraph*
g =
new TGraph;
810 g->SetPoint(0,
sign*range, 0);
811 g->SetPoint(1,
sign*range, 1e4);
820 TGraph*
g =
new TGraph;
821 g->SetPoint(0, rangedown, 0);
822 g->SetPoint(1, rangedown, 1e4);
825 TGraph* g2 =
new TGraph;
826 g2->SetPoint(0, rangeup, 0);
827 g2->SetPoint(1, rangeup, 1e4);
844 const TString dataMCStr =
fIsMC ?
"_mc" :
"_data";
845 const TString viewStrArr[4]={
"fitX",
"fitY",
"fitMuX",
"fitMuY"};
846 const TString viewStr =viewStrArr[
view];
856 TString suffix = dataMCStr+viewStr;
861 if(cell >= 0) quiet =
"Q";
868 if (
fIsMC) title +=
"Monte Carlo";
else title +=
"data";
870 TString viewStrArrVertHor[4]={
"vertical",
"horizontal",
"vertical",
"horizontal"};
871 TString viewStrVertHor =viewStrArrVertHor[
view];
875 viewStrVertHor.Data(),
883 atten->SetTitle(title+
";Distance from center (cm);Mean PE / cm");
892 const double attenErr = fit->GetParError(0);
894 atten->GetYaxis()->SetRangeUser(0, 120);
896 TLegend*
leg =
new TLegend(.3, .15, .7, .35);
897 leg->SetFillStyle(0);
906 gSystem->mkdir(plotsDirectory.c_str());
907 gPad->Print(plotsDirectory+
"/atten"+suffix+
".eps");
911 TH1D* profCorr = (TH1D*)atten->Clone();
912 profCorr->GetYaxis()->SetTitle(
"Ratio to exponential fit");
913 fit->SetRange(-1.1*cell_length, +1.1*cell_length);
914 profCorr->Divide(fit);
922 profCorr->GetYaxis()->SetRangeUser(.4, 1.1);
924 leg =
new TLegend(.3, .15, .7, .35);
925 leg->SetFillStyle(0);
938 if(!plotsDirectory.empty()){
939 gPad->Print(plotsDirectory+
"/rolloff"+suffix+
".eps");
944 TH1* profFullCorr = (TH1*)atten->Clone();
945 profFullCorr->GetYaxis()->SetTitle(
"Ratio to full fit");
946 rollFit->SetRange(-800, +800);
948 rollFit->SetRange(-250, +250);
950 profFullCorr->Divide(rollFit);
951 profFullCorr->GetYaxis()->SetRangeUser(0.9, 1.1);
952 profFullCorr->Draw();
954 TGraph*
one =
new TGraph;
955 one->SetPoint(0, -1e4, 1);
956 one->SetPoint(1, +1e4, 1);
957 one->SetLineStyle(7);
968 TGraph* lowess =
lowessFit(profFullCorr, plane, res);
969 lowess->Draw(
"lp same");
974 if(!plotsDirectory.empty()){
975 gPad->Print(plotsDirectory+
"/residual"+suffix+
".eps");
980 TH1* plotProf = (TH1*)atten->Clone();
983 struct FullFitAdaptor
986 double operator()(
double* xs,
double*
ps)
988 return fRes->MeanPEPerCmAt(xs[0]);
993 const double range = cell_length/2;
994 double rangeup=
range;
998 TF1* fullFit =
new TF1(
"fullFit", adapt, -range, +rangeup, 0);
999 fullFit->SetLineColor(
kBlue);
1000 fullFit->SetRange(-range, +rangeup);
1001 fullFit->Draw(
"same");
1003 plotProf->SetLineColor(kBlack);
1004 plotProf->GetYaxis()->SetTitle(
"Mean PE / cm");
1005 plotProf->GetYaxis()->SetRangeUser(0, 150);
1008 plotProf->GetXaxis()->CenterTitle();
1009 plotProf->GetYaxis()->CenterTitle();
1010 plotProf->GetYaxis()->SetTitleSize(plotProf->GetXaxis()->GetTitleSize());
1011 plotProf->GetYaxis()->SetTitleOffset(plotProf->GetXaxis()->GetTitleOffset());
1020 if(!plotsDirectory.empty()){
1021 gPad->Print(plotsDirectory+
"/totfit"+suffix+
".eps");
1022 gPad->Print(plotsDirectory+
"/totfit"+suffix+
".png");
1026 TH1* resid =
new TH1F(
"",
";W (cm);Ratio to total fit", atten->GetNbinsX(), atten->GetXaxis()->GetXmin(), atten->GetXaxis()->GetXmax());
1028 for(
int n = 0;
n < atten->GetNbinsX()+2; ++
n){
1029 resid->SetBinContent(
n, atten->GetBinContent(
n)/res->
MeanPEPerCmAt(atten->GetBinCenter(
n)));
1031 resid->GetYaxis()->SetRangeUser(.9, 1.1);
1033 one->Draw(
"l same");
1035 if(!plotsDirectory.empty()){
1036 gPad->Print(plotsDirectory+
"/residual_final"+suffix+
".eps");
1042 delete profFullCorr;
1053 const int vldTime = 1;
1058 assert(plane == 0 || plane == 1 || plane == 2 || plane == 3);
1062 fprintf(f,
"-1,-1,-1,-1,-1,-1,-1,1000,%d,%d\n",
caldp::AttenHists GetChannelHists(bool isData, int plane, int cell)
void fit_channel_prof(TH1 *atten, int plane, int cell, calib::AttenCurve *&res)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
int isinf(const stan::math::var &a)
bool HasProfiles(const geo::OfflineChan &c) const
std::map< std::string, double > xmax
fvar< T > fabs(const fvar< T > &x)
void DrawDetectorEdges(double range)
virtual void beginRun(const art::Run &run)
double ThresholdCorrAt(int view, int cell, double w) const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void reconfigure(const fhicl::ParameterSet &pset)
Float_t x1[n_points_granero]
TF1 * expFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
static AttenCurve * Uninitialized(int det, geo::OfflineChan chan)
Return a new AttenCurve objects with fields uninitialized.
TH2 * GetBestHistogram(caldp::AttenHists v, caldp::AttenProfiles u)
Vertical planes which measure X.
unsigned int Ncells() const
Number of cells in this plane.
TGraph * lowessFit(TH1 *prof, int plane, calib::AttenCurve *&res)
AttenProfiles for many channels.
virtual void endSubRun(const art::SubRun &sr)
TH1 * MedianProfile(TH2 *f)
double corr(TGraph *g, double thresh)
std::vector< geo::OfflineChan > Channels() const
List of channels with profiles on.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
TH1 * TruncatedMeanProfile(TH2 *h)
const PlaneGeo * Plane(unsigned int i) const
std::string fThresholdCorrMapFile
TH1 * MeanProfile(TH2 *f)
T cube(T x)
More efficient cube function than pow(x,3)
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
int isnan(const stan::math::var &a)
AttenFit(const fhicl::ParameterSet &pset)
int ToDBValidityChan() const
View_t View() const
Which coordinate does this plane measure.
Far Detector at Ash River, MN.
block
print "ROW IS " print row
std::map< int, caldp::AttenProfilesMap > fChannelMapProf
bool fConsolidateViewsMode
void WriteDummyCSVRow(FILE *f, int plane, int cell) const
Ensure every cell has a CSV entry even if it couldn't be fit.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
AttenHists & GetHists(const geo::OfflineChan &c)
Near Detector in the NuMI cavern.
virtual void analyze(const art::Event &evt)
void ApplyThresholdCorrection(int view, int cell, TH1 *prof)
Histograms used by attenuation calibration.
std::string fAttenHistsMapLabel
double frac(double x)
Fractional part.
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
AttenHists for many channels.
unsigned short Plane() const
Profiles used by attenuation calibration.
unsigned short Cell() const
caldp::AttenHistsMap fChannelMap
This is where GetChannelHist stores its histograms. Map is from block.
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
virtual void endRun(const art::Run &run)
assert(nhit_max >=nhit_nbins)
int GetView(int plane, art::ServiceHandle< geo::Geometry > geom)
float MeanPEPerCmAt(double w) const
Mean response of this channel at this distance from detector centre.
caldp::AttenProfiles GetChannelProfiles(bool isData, int plane, int cell)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
bool removeCachedProduct(Handle< PROD > &h) const
unsigned int NPlanes() const
nova::dbi::DataType GetDataType()
std::string fPlotsDirectory
TF1 * rolloffFit(TH1 *prof, TString opt, geo::OfflineChan chan, calib::AttenCurve *&res)
AttenProfiles & GetProfiles(const geo::OfflineChan &c)
double FitQuality(TH1 *prof, TF1 *fit) const
Encapsulate the geometry of one entire detector (near, far, ndos)
void WriteToCSVs(FILE *fConsts, FILE *fPoints, bool mc) const
ThresholdCorrMap * fThreshCorrMap
void AddInterpPoint(float w, float factor)
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.