diff --git a/offline/packages/mbd/MbdCalib.cc b/offline/packages/mbd/MbdCalib.cc index 3ebe139ffb..e8d10ba00a 100644 --- a/offline/packages/mbd/MbdCalib.cc +++ b/offline/packages/mbd/MbdCalib.cc @@ -2497,7 +2497,7 @@ void MbdCalib::Reset_Pileup() _pileup_p0err.fill(std::numeric_limits::quiet_NaN()); _pileup_p1err.fill(std::numeric_limits::quiet_NaN()); _pileup_p2err.fill(std::numeric_limits::quiet_NaN()); - _qfit_chi2ndf.fill(std::numeric_limits::quiet_NaN()); + _pileup_chi2ndf.fill(std::numeric_limits::quiet_NaN()); } void MbdCalib::Reset_Thresholds() @@ -2590,3 +2590,31 @@ TGraph *MbdCalib::get_lut_graph(const int pmtch, std::string_view type) return g; } + +void MbdCalib::set_pileup(const int ifeech, const int ipar, const float pval) +{ + int chtype = (ifeech / 8) % 2; // 0=T-ch, 1=Q-ch + + if (ipar==0) + { + _pileup_p0[ifeech] = pval; + } + else if (ipar==1) + { + _pileup_p1[ifeech] = pval; + } + else if (ipar==2) + { + _pileup_p2[ifeech] = pval; + } + else if (ipar==3 && chtype==0) + { + _pileup_p1err[ifeech] = pval; + } + else if (ipar==4 && chtype==0) + { + _pileup_p2err[ifeech] = pval; + } +} + + diff --git a/offline/packages/mbd/MbdCalib.h b/offline/packages/mbd/MbdCalib.h index 49129ac749..6f0a6e579a 100644 --- a/offline/packages/mbd/MbdCalib.h +++ b/offline/packages/mbd/MbdCalib.h @@ -38,6 +38,36 @@ class MbdCalib float get_pedrms(const int ifeech) const { return _pedsigma[ifeech]; } int get_sampmax(const int ifeech) const { return _sampmax[ifeech]; } int get_status(const int ifeech) const { return _mbdstatus[ifeech]; } + + float get_pileup(const int ifeech, const int ipar) const + { + int chtype = (ifeech / 8) % 2; // 0=T-ch, 1=Q-ch + + if (ipar==0) + { + return _pileup_p0[ifeech]; + } + else if (ipar==1) + { + return _pileup_p1[ifeech]; + } + else if (ipar==2) + { + return _pileup_p2[ifeech]; + } + else if (ipar==3 && chtype==0) + { + return _pileup_p1err[ifeech]; + } + else if (ipar==4 && chtype==0) + { + return _pileup_p2err[ifeech]; + } + + return std::numeric_limits::quiet_NaN(); + } + + float get_tcorr(const int ifeech, const int tdc) const { if (tdc<0) { @@ -86,24 +116,6 @@ class MbdCalib std::vector get_shape(const int ifeech) const { return _shape_y[ifeech]; } std::vector get_sherr(const int ifeech) const { return _sherr_yerr[ifeech]; } - float get_pileup(const int ifeech, const int ipar) const { - - if (ipar==0) - { - return _pileup_p0[ifeech]; - } - else if (ipar==1) - { - return _pileup_p1[ifeech]; - } - else if (ipar==2) - { - return _pileup_p2[ifeech]; - } - - return std::numeric_limits::quiet_NaN(); - } - float get_threshold(const int pmtch, const int rel_or_abs = 0); TGraph *get_lut_graph(const int pmtch, std::string_view type); @@ -111,6 +123,7 @@ class MbdCalib void set_sampmax(const int ifeech, const int val) { _sampmax[ifeech] = val; } void set_status(const int ifeech, const int val) { _mbdstatus[ifeech] = val; } void set_ped(const int ifeech, const float m, const float merr, const float s, const float serr); + void set_pileup(const int ifeech, const int ipar, const float val); void set_tt0(const int ipmt, const float t0) { _ttfit_t0mean[ipmt] = t0; } void set_tq0(const int ipmt, const float t0) { _tqfit_t0mean[ipmt] = t0; } diff --git a/offline/packages/mbd/MbdEvent.cc b/offline/packages/mbd/MbdEvent.cc index 72cecb212d..62a6c96fe7 100644 --- a/offline/packages/mbd/MbdEvent.cc +++ b/offline/packages/mbd/MbdEvent.cc @@ -380,6 +380,7 @@ int MbdEvent::End() for (auto & sig : _mbdsig) { + sig.WritePedvsEvent(); sig.WriteChi2Hist(); } @@ -1127,11 +1128,6 @@ int MbdEvent::Calculate(MbdPmtContainer *bbcpmts, MbdOut *bbcout, PHCompositeNod gausfit[iarm]->SetRange(hevt_bbct[iarm]->GetMean() - 5, hevt_bbct[iarm]->GetMean() + 5); */ - if ( hevt_bbct[iarm]->GetEntries()==0 )//chiu - { - std::cout << PHWHERE << " hevt_bbct EMPTY" << std::endl; - } - hevt_bbct[iarm]->Fit(gausfit[iarm], "BNQLR"); // m_bbct[iarm] = m_bbct[iarm] / m_bbcn[iarm]; @@ -1462,10 +1458,6 @@ int MbdEvent::CalcPedCalib() pedgaus->SetParameters(ampl,mean,sigma); pedgaus->SetRange(mean-(4*sigma), mean+(4*sigma)); - if ( hped0->GetEntries()==0 ) //chiu - { - std::cout << "HPED0 EMPTY" << std::endl; - } hped0->Fit(pedgaus,"RNQ"); mean = pedgaus->GetParameter(1); diff --git a/offline/packages/mbd/MbdSig.cc b/offline/packages/mbd/MbdSig.cc index 1118ba2e60..9b270e8ee8 100644 --- a/offline/packages/mbd/MbdSig.cc +++ b/offline/packages/mbd/MbdSig.cc @@ -49,12 +49,13 @@ void MbdSig::Init() gSubPulse->SetName(name); gSubPulse->GetHistogram()->SetXTitle("sample"); gSubPulse->GetHistogram()->SetYTitle("ADC"); + gSubPulse->GetHistogram()->SetTitle(name); hpulse = hRawPulse; // hpulse,gpulse point to raw by default gpulse = gRawPulse; // we switch to sub for default if ped is applied - //ped0stats = std::make_unique(100); // use the last 100 events for running pedestal - ped0stats = new MbdRunningStats(100); // use the last 100 events for running pedestal + ped0stats = new MbdRunningStats(8); // use the last 8 samples for running pedestal + name = "hPed0_"; name += _ch; hPed0 = new TH1F(name, name, 3000, -0.5, 2999.5); @@ -62,6 +63,15 @@ void MbdSig::Init() name = "hPedEvt_"; name += _ch; hPedEvt = new TH1F(name, name, 3000, -0.5, 2999.5); + if ( _pedstudyflag ) + { + gPedvsEvent = new TGraphErrors(); + name = "gpedvsevent"; + name += _ch; + gPedvsEvent->SetName(name); + gPedvsEvent->GetHistogram()->SetXTitle("evtnum"); + gPedvsEvent->GetHistogram()->SetYTitle("ped"); + } SetTemplateSize(900, 1000, -10., 20.); // SetTemplateSize(300,300,0.,15.); @@ -137,9 +147,9 @@ MbdSig::~MbdSig() delete hSubPulse; delete gRawPulse; delete gSubPulse; - delete hPed0; delete ped0stats; - // h2Template->Write(); + delete hPed0; + delete hPedEvt; delete h2Template; delete h2Residuals; delete hAmpl; @@ -149,6 +159,10 @@ MbdSig::~MbdSig() delete ped_fcn; delete ped_tail; delete h_chi2ndf; + if ( _pedstudyflag ) + { + delete gPedvsEvent; + } } void MbdSig::SetEventPed0PreSamp(const Int_t presample, const Int_t nsamps, const int max_samp) @@ -327,54 +341,134 @@ void MbdSig::Remove_Pileup() if ( (_ch/8)%2 == 0 ) // time ch { - float offset = _pileup_p0*gSubPulse->GetPointY(0); + double x_at_max = TMath::LocMax( 5, gSubPulse->GetY() ); - for (int isamp = 0; isamp < _nsamples; isamp++) + if ( x_at_max != 0 ) { - double x = gSubPulse->GetPointX(isamp); - double y = gSubPulse->GetPointY(isamp); + // time hit in prev crossing + if ( fit_pileup == nullptr ) + { + TString name = "fit_pileup"; name += _ch; + fit_pileup = new TF1(name,"pol3",0,16000); + for (int ipar=0; ipar<4; ipar++) + { + fit_pileup->SetParameter( ipar, _mbdcal->get_pileup(_ch,ipar+1) ); + } + } - hSubPulse->SetBinContent( isamp + 1, y - offset ); - gSubPulse->SetPoint( isamp, x, y - offset ); + int sampmax = _mbdcal->get_sampmax(_ch); + if ( (sampmax-6) > 0 ) + { + double x_sampmax = gSubPulse->GetPointX(sampmax); + double y_sampmax = gSubPulse->GetPointY(sampmax); + double y_min6 = gSubPulse->GetPointY(sampmax-6); + + double offset = y_min6*fit_pileup->Eval(y_min6); + + hSubPulse->SetBinContent( sampmax + 1, y_sampmax - offset ); + gSubPulse->SetPoint( sampmax, x_sampmax, y_sampmax - offset ); + } + else + { + static int ctr = 0; + if ( ctr<10 ) + { + std::cout << PHWHERE << " WARNING, sampmax too early for time pileup corr" << std::endl; + ctr++; + } + } } - } - else - { - if ( fit_pileup == nullptr ) + else { - TString name = "fit_pileup"; name += _ch; - fit_pileup = new TF1(name,"gaus",-0.1,4.1); - fit_pileup->SetLineColor(2); - } + // time hit in 2 crossings before + float offset = _pileup_p0*gSubPulse->GetPointY(0); - fit_pileup->SetRange(-0.1,4.1); - fit_pileup->SetParameters( _pileup_p0*gSubPulse->GetPointY(0), _pileup_p1, _pileup_p2 ); - - // fix par limits - double plow{0.}; - double phigh{0.}; - fit_pileup->GetParLimits(2,plow,phigh); - if ( phigh < _pileup_p2 ) - { - phigh = 2*_pileup_p2; - fit_pileup->SetParLimits(2,plow,phigh); + for (int isamp = 0; isamp < _nsamples; isamp++) + { + double x = gSubPulse->GetPointX(isamp); + double y = gSubPulse->GetPointY(isamp); + + hSubPulse->SetBinContent( isamp + 1, y - offset ); + gSubPulse->SetPoint( isamp, x, y - offset ); + } } + } + else // charge ch + { + double ymax = TMath::MaxElement( 5, gSubPulse->GetY() ); + double x_at_max = TMath::LocMax( 5, gSubPulse->GetY() ); - if ( _verbose ) + if ( x_at_max != 0 ) { - gSubPulse->Fit( fit_pileup, "R" ); - gSubPulse->Draw("ap"); - PadUpdate(); + // Fit a pulse in prev crossing + template_fcn->SetParameters(ymax, x_at_max); + template_fcn->SetRange(0, x_at_max+2.1); + + if (_verbose == 0) + { + //std::cout << PHWHERE << std::endl; + gSubPulse->Fit(template_fcn, "RNQ"); + } + else + { + std::cout << "pre-pileup " << _ch << "\t" << x_at_max << "\t" << ymax << std::endl; + gSubPulse->Fit(template_fcn, "R"); + gSubPulse->Draw("ap"); + gSubPulse->GetHistogram()->SetTitle(gSubPulse->GetName()); + gPad->SetGridy(1); + PadUpdate(); + //gSubPulse->Print("ALL"); + } + } else { - gSubPulse->Fit( fit_pileup, "RNQ" ); + // Fit the tail + if ( fit_pileup == nullptr ) + { + TString name = "fit_pileup"; name += _ch; + fit_pileup = new TF1(name,"gaus",-0.1,4.1); + fit_pileup->SetLineColor(2); + } + + fit_pileup->SetRange(-0.1,4.1); + fit_pileup->SetParameters( _pileup_p0*gSubPulse->GetPointY(0), _pileup_p1, _pileup_p2 ); + + // fix par limits + double plow{0.}; + double phigh{0.}; + fit_pileup->GetParLimits(2,plow,phigh); + if ( phigh < _pileup_p2 ) + { + phigh = 2*_pileup_p2; + fit_pileup->SetParLimits(2,plow,phigh); + } + + if ( _verbose ) + { + gSubPulse->Fit( fit_pileup, "R" ); + gSubPulse->Draw("ap"); + PadUpdate(); + } + else + { + gSubPulse->Fit( fit_pileup, "RNQ" ); + } } + // subtract pre-pulse for (int isamp = 0; isamp < _nsamples; isamp++) { - double bkg = fit_pileup->Eval(isamp); + double bkg = 0.; + if ( x_at_max != 0 ) + { + bkg = template_fcn->Eval(isamp); + } + else + { + bkg = fit_pileup->Eval(isamp); + } double x = gSubPulse->GetPointX(isamp); double y = gSubPulse->GetPointY(isamp); @@ -388,6 +482,7 @@ void MbdSig::Remove_Pileup() if ( _verbose ) { + std::cout << "pileup sub " << _ch << std::endl; gSubPulse->Draw("ap"); PadUpdate(); } @@ -428,6 +523,14 @@ void MbdSig::WritePedHist() hPed0->Write(); } +void MbdSig::WritePedvsEvent() +{ + if ( _pedstudyflag ) + { + gPedvsEvent->Write(); + } +} + void MbdSig::FillPed0(const Int_t sampmin, const Int_t sampmax) { Double_t x; @@ -438,13 +541,6 @@ void MbdSig::FillPed0(const Int_t sampmin, const Int_t sampmax) // gRawPulse->Print("all"); hPed0->Fill(y); - /* - // chiu taken out - ped0stats->Push( y ); - ped0 = ped0stats->Mean(); - ped0rms = ped0stats->RMS(); - */ - // std::cout << "ped0 " << _ch << " " << n << "\t" << ped0 << std::endl; // std::cout << "ped0 " << _ch << "\t" << ped0 << std::endl; } @@ -463,10 +559,10 @@ void MbdSig::FillPed0(const Double_t begin, const Double_t end) hPed0->Fill(y); /* - ped0stats->Push( y ); - ped0 = ped0stats->Mean(); - ped0rms = ped0stats->RMS(); - */ + ped0stats->Push( y ); + ped0 = ped0stats->Mean(); + ped0rms = ped0stats->RMS(); + */ // std::cout << "ped0 " << _ch << " " << n << "\t" << x << "\t" << y << std::endl; } @@ -548,9 +644,7 @@ void MbdSig::CalcEventPed0(const Double_t minpedx, const Double_t maxpedx) // If a prev event pileup is detected, return 1, otherwise, return 0 int MbdSig::CalcEventPed0_PreSamp(const int presample, const int nsamps) { - //std::cout << PHWHERE << std::endl; //chiu //_verbose = 100; - //ped0stats->Clear(); int status = 0; // assume no pileup @@ -608,11 +702,6 @@ int MbdSig::CalcEventPed0_PreSamp(const int presample, const int nsamps) ped_fcn->SetRange(minsamp-0.1,maxsamp+0.1); ped_fcn->SetParameter(0,1500.); - if ( gRawPulse->GetN()==0 )//chiu - { - std::cout << PHWHERE << " gRawPulse 0" << std::endl; - } - if ( _verbose ) { gRawPulse->Fit( ped_fcn, "RQ" ); @@ -630,18 +719,6 @@ int MbdSig::CalcEventPed0_PreSamp(const int presample, const int nsamps) //std::cout << PHWHERE << std::endl; gRawPulse->Fit( ped_fcn, "RNQ" ); - /* - double chi2ndf = ped_fcn->GetChisquare()/ped_fcn->GetNDF(); - if ( _pileupfile != nullptr && chi2ndf > 4.0 ) - { - *_pileupfile << "ped " << _ch << " mean " << mean << "\t"; - for ( int i=0; iGetN(); i++) - { - *_pileupfile << std::setw(6) << gRawPulse->GetPointY(i); - } - *_pileupfile << std::endl; - } - */ } double chi2 = ped_fcn->GetChisquare(); @@ -670,6 +747,28 @@ int MbdSig::CalcEventPed0_PreSamp(const int presample, const int nsamps) << "isamp " << isamp << "\t" << x << "\t" << y << std::endl; } } + + // study pedestal vs event + if ( _pedstudyflag ) + { + // running pedestal (replace with mean and meanerr for evt-by-evt) + double ped_evtnum = _evt_counter; + double ped_mean = ped0stats->Mean(); + double ped_meanerr = rms; + if ( ped0stats->Size()>1 ) + { + ped_meanerr = ped0stats->RMS()/std::sqrt(ped0stats->Size()); + } + else + { + ped_meanerr = _mbdcal->get_pedrms(_ch); + } + + int n = gPedvsEvent->GetN(); + gPedvsEvent->SetPoint(n,ped_evtnum,ped_mean); + gPedvsEvent->SetPointError(n,0,ped_meanerr); + } + } else { @@ -700,11 +799,18 @@ int MbdSig::CalcEventPed0_PreSamp(const int presample, const int nsamps) } } - // use straight mean for pedestal - // Could consider using fit to hPed0 to remove outliers - //rms = ped0stats->RMS(); - //Double_t mean = hPed0->GetMean(); - //Double_t rms = hPed0->GetRMS(); + // uncomment this to write out file with pileup waveforms + /* + if ( _pileupfile != nullptr ) + { + *_pileupfile << "ped " << _ch << " mean " << mean << "\t"; + for ( int i=0; iGetN(); i++) + { + *_pileupfile << std::setw(6) << gRawPulse->GetPointY(i); + } + *_pileupfile << std::endl; + } + */ } SetPed0(mean, rms); @@ -982,9 +1088,10 @@ void MbdSig::PadUpdate() const std::cout << PHWHERE << " PadUpdate\t_verbose = " << _verbose << std::endl; if ( _verbose>5 ) { + gPad->SetGridy(1); gPad->Modified(); gPad->Update(); - std::cout << _ch << " ? "; + std::cout << _evt_counter << ": " << _ch << " ? "; if ( _verbose>10 ) { std::string junk; @@ -1124,7 +1231,7 @@ Double_t MbdSig::TemplateFcn(const Double_t* x, const Double_t* par) // sampmax>0 means fit to the peak near sampmax int MbdSig::FitTemplate( const Int_t sampmax ) { - //std::cout << PHWHERE << std::endl; //chiu + //std::cout << PHWHERE << std::endl; /* if ( _evt_counter==2142 && _ch==92 ) { @@ -1265,7 +1372,7 @@ int MbdSig::FitTemplate( const Int_t sampmax ) // fit was out of time, likely from pileup, try two waveforms if ( (f_time<(sampmax-2.5) || f_time>sampmax) && (nsaturated<=3) ) { - //_verbose = 100; //chiu + //_verbose = 100; if ( _verbose ) { @@ -1310,15 +1417,6 @@ int MbdSig::FitTemplate( const Int_t sampmax ) } } - - /* chiu - if ( f_time<0. || f_time>9 ) - { - _verbose = 100; - f_time = _nsamples*0.5; // bad fit last time - } - */ - // refit with new range to exclude after-pulses template_fcn->SetParameters(ymax, x_at_max); //template_fcn->SetParameters( f_ampl, f_time ); @@ -1369,14 +1467,6 @@ int MbdSig::FitTemplate( const Int_t sampmax ) h_chi2ndf->Fill( f_chi2/f_ndf ); - /* - if ( (f_chi2/f_ndf) > 100. ) //chiu - { - std::cout << "very bad chi2ndf after refit " << f_ampl << "\t" << f_time << "\t" << f_chi2/f_ndf << std::endl; - //_verbose = 100; - } - */ - //if ( f_time<0 || f_time>30 ) //if ( (_ch==185||_ch==155||_ch==249) && (fabs(f_ampl) > 44000.) ) //double chi2 = template_fcn->GetChisquare(); diff --git a/offline/packages/mbd/MbdSig.h b/offline/packages/mbd/MbdSig.h index a3ee986c33..bafdd9339c 100644 --- a/offline/packages/mbd/MbdSig.h +++ b/offline/packages/mbd/MbdSig.h @@ -116,6 +116,7 @@ class MbdSig void SetMinMaxFitTime(const Double_t mintime, const Double_t maxtime); void WritePedHist(); + void WritePedvsEvent(); void WriteChi2Hist(); void DrawWaveform(); /// Draw Subtracted Waveform @@ -157,22 +158,22 @@ class MbdSig TGraphErrors *gpulse{nullptr}; //! /** for CalcPed0 */ - //std::unique_ptr ped0stats{nullptr}; //! - MbdRunningStats *ped0stats{nullptr}; //! - TH1 *hPed0{nullptr}; //! all events - TH1 *hPedEvt{nullptr}; //! evt-by-event pedestal + MbdRunningStats *ped0stats{nullptr}; //! running pedestal + TH1 *hPed0{nullptr}; //! all events + TH1 *hPedEvt{nullptr}; //! evt-by-event pedestal + TGraphErrors *gPedvsEvent{nullptr}; //! Keep track of pedestal vs evtnum TF1 *ped_fcn{nullptr}; - TF1 *ped_tail{nullptr}; //! tail of prev signal + TF1 *ped_tail{nullptr}; //! tail of prev signal Double_t ped0{0.}; //! Double_t ped0rms{0.}; //! - int use_ped0{0}; //! whether to apply ped0 - Int_t minped0samp{-9999}; //! min sample for event-by-event ped, inclusive - Int_t maxped0samp{-9999}; //! max sample for event-by-event ped, inclusive + int use_ped0{0}; //! whether to apply ped0 + Int_t minped0samp{-9999}; //! min sample for event-by-event ped, inclusive + Int_t maxped0samp{-9999}; //! max sample for event-by-event ped, inclusive Double_t minped0x{0.}; //! min x for event-by-event ped, inclusive Double_t maxped0x{0.}; //! max x for event-by-event ped, inclusive - Double_t ped_presamp{}; //! presamples for ped calculation - Double_t ped_presamp_nsamps{}; //! num of presamples for ped calculation - Double_t ped_presamp_maxsamp{-1}; //! a peak sample for ped calc (-1 = use max) + Double_t ped_presamp{}; //! presamples for ped calculation + Double_t ped_presamp_nsamps{}; //! num of presamples for ped calculation + Double_t ped_presamp_maxsamp{-1}; //! a peak sample for ped calc (-1 = use max) /** for time calibration */ // Double_t time_calib; @@ -200,6 +201,7 @@ class MbdSig TH1 *h_chi2ndf{nullptr}; //! for eval int _verbose{0}; + bool _pedstudyflag{false}; }; #endif // __MBDSIG_H__