From 9515c5b0dc97ee9e182a1593e8f5caba9a87848f Mon Sep 17 00:00:00 2001 From: arabella Date: Wed, 19 Dec 2018 16:01:31 +0100 Subject: [PATCH] revisit MCP wf reconstruction --- RawToDigi/plugins/HGCalTBCAENProducer.cc | 133 +++++--- RawToDigi/plugins/HGCalTBCAENProducer.h | 93 ++++-- RawToDigi/plugins/HGCalTBCAENSource.cc | 165 ++++++---- RawToDigi/plugins/HGCalTBCAENSource.h | 2 +- Reco/interface/MCPwaveform.h | 95 ++++++ Reco/plugins/CAEN_NTupelizer.cc | 135 +++++--- Reco/plugins/MCPNtupler.cc | 107 ++++--- Reco/src/MCPwaveform.cc | 390 +++++++++++++++++++++++ 8 files changed, 890 insertions(+), 230 deletions(-) create mode 100644 Reco/interface/MCPwaveform.h create mode 100644 Reco/src/MCPwaveform.cc diff --git a/RawToDigi/plugins/HGCalTBCAENProducer.cc b/RawToDigi/plugins/HGCalTBCAENProducer.cc index c22f91c..6b52e01 100644 --- a/RawToDigi/plugins/HGCalTBCAENProducer.cc +++ b/RawToDigi/plugins/HGCalTBCAENProducer.cc @@ -61,21 +61,26 @@ void HGCalTBBeamWireChamberProducer::beginJob() { tree->SetBranchAddress("N_scintillator_veto_timestamps", &scintillator_vetos, &b_scintillator_vetos); tree->SetBranchAddress("valid_TS_MCP1", &valid_TS_MCP1, &b_valid_TS_MCP1); tree->SetBranchAddress("valid_TS_MCP2", &valid_TS_MCP2, &b_valid_TS_MCP2); - tree->SetBranchAddress("TS_15PercentRise_MCP1", &TS_15PercentRise_MCP1, &b_TS_15PercentRise_MCP1); - tree->SetBranchAddress("TS_15PercentRise_MCP2", &TS_15PercentRise_MCP2, &b_TS_15PercentRise_MCP2); - tree->SetBranchAddress("TS_30PercentRise_MCP1", &TS_30PercentRise_MCP1, &b_TS_30PercentRise_MCP1); - tree->SetBranchAddress("TS_30PercentRise_MCP2", &TS_30PercentRise_MCP2, &b_TS_30PercentRise_MCP2); - tree->SetBranchAddress("TS_45PercentRise_MCP1", &TS_45PercentRise_MCP1, &b_TS_45PercentRise_MCP1); - tree->SetBranchAddress("TS_45PercentRise_MCP2", &TS_45PercentRise_MCP2, &b_TS_45PercentRise_MCP2); - tree->SetBranchAddress("TS_60PercentRise_MCP1", &TS_60PercentRise_MCP1, &b_TS_60PercentRise_MCP1); - tree->SetBranchAddress("TS_60PercentRise_MCP2", &TS_60PercentRise_MCP2, &b_TS_60PercentRise_MCP2); + tree->SetBranchAddress("noise_MCP1", &noise_MCP1, &b_noise_MCP1); + tree->SetBranchAddress("noise_MCP2", &noise_MCP2, &b_noise_MCP1); + tree->SetBranchAddress("TSpeak_MCP1", &TSpeak_MCP1, &b_TSpeak_MCP1); + tree->SetBranchAddress("TSpeak_MCP2", &TSpeak_MCP2, &b_TSpeak_MCP2); tree->SetBranchAddress("amp_MCP1", &_MCP1, &b_amp_MCP1); tree->SetBranchAddress("amp_MCP2", &_MCP2, &b_amp_MCP2); - - tree->SetBranchAddress("TS_MCP1", &TS_MCP1, &b_TS_MCP1); - tree->SetBranchAddress("TS_MCP2", &TS_MCP2, &b_TS_MCP2); - tree->SetBranchAddress("TS_MCP1_to_last_falling_Edge", &TS_MCP1_to_last_falling_Edge, &b_TS_MCP1_to_last_falling_Edge); - tree->SetBranchAddress("TS_MCP2_to_last_falling_Edge", &TS_MCP2_to_last_falling_Edge, &b_TS_MCP2_to_last_falling_Edge); + tree->SetBranchAddress("ampFit_MCP1", &Fit_MCP1, &b_ampFit_MCP1); + tree->SetBranchAddress("ampFit_MCP2", &Fit_MCP2, &b_ampFit_MCP2); + tree->SetBranchAddress("TSfitPeak_MCP1", &TSfitPeak_MCP1, &b_TSfitPeak_MCP1); + tree->SetBranchAddress("TSfitPeak_MCP2", &TSfitPeak_MCP2, &b_TSfitPeak_MCP2); + tree->SetBranchAddress("TScf_MCP1", &TScf_MCP1, &b_TScf_MCP1); + tree->SetBranchAddress("TScf_MCP2", &TScf_MCP2, &b_TScf_MCP2); + tree->SetBranchAddress("charge5nsS_MCP1", &charge5nsS_MCP1, &b_charge5nsS_MCP1); + tree->SetBranchAddress("charge5nsS_MCP2", &charge5nsS_MCP2, &b_charge5nsS_MCP2); + tree->SetBranchAddress("charge5nsB_MCP1", &charge5nsB_MCP1, &b_charge5nsB_MCP1); + tree->SetBranchAddress("charge5nsB_MCP2", &charge5nsB_MCP2, &b_charge5nsB_MCP2); + tree->SetBranchAddress("TS_toClock_FE_MCP1", &TS_toClock_FE_MCP1, &b_TS_toClock_FE_MCP1); + tree->SetBranchAddress("TS_toClock_FE_MCP2", &TS_toClock_FE_MCP2, &b_TS_toClock_FE_MCP2); + tree->SetBranchAddress("meanClockFE", &meanClockFE, &b_meanClockFE); + tree->SetBranchAddress("rmsClockFE", &rmsClockFE, &b_rmsClockFE); loaded_run = -1; } @@ -219,22 +224,28 @@ void HGCalTBBeamWireChamberProducer::produce(edm::Event& event, const edm::Event rd_full->intUserRecords.add("scintillator_coincidence_timestamps", scintillator_coincidences_loaded[event_nr]); rd_full->intUserRecords.add("scintillator_coincidence_timestamps", scintillator_vetos_loaded[event_nr]); //MCPs - rd_full->intUserRecords.add("valid_TS_MCP1", valid_TS_MCP1_loaded[event_nr]); - rd_full->intUserRecords.add("valid_TS_MCP2", valid_TS_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_MCP1", TS_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_MCP2", TS_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_15PercentRise_MCP1", TS_15PercentRise_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_15PercentRise_MCP2", TS_15PercentRise_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_30PercentRise_MCP2", TS_30PercentRise_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_30PercentRise_MCP1", TS_30PercentRise_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_45PercentRise_MCP1", TS_45PercentRise_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_45PercentRise_MCP2", TS_45PercentRise_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_60PercentRise_MCP1", TS_60PercentRise_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_60PercentRise_MCP2", TS_60PercentRise_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("amp_MCP1", amp_MCP1_loaded[event_nr]); - rd_full->doubleUserRecords.add("amp_MCP2", amp_MCP2_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_MCP1_to_last_falling_Edge", TS_MCP1_to_last_falling_Edge_loaded[event_nr]); - rd_full->doubleUserRecords.add("TS_MCP2_to_last_falling_Edge", TS_MCP2_to_last_falling_Edge_loaded[event_nr]); + rd_full->intUserRecords.add("valid_TS_MCP1", valid_TS_MCP1_loaded[event_nr]); + rd_full->intUserRecords.add("valid_TS_MCP2", valid_TS_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("noise_MCP1", noise_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("noise_MCP2", noise_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("TSpeak_MCP1", TSpeak_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("TSpeak_MCP2", TSpeak_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("amp_MCP1", amp_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("amp_MCP2", amp_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("ampFit_MCP1", ampFit_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("ampFit_MCP2", ampFit_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("TSfitPeak_MCP1", TSfitPeak_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("TSfitPeak_MCP2", TSfitPeak_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("TScf_MCP1", TScf_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("TScf_MCP2", TScf_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("charge5nsS_MCP1", charge5nsS_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("charge5nsS_MCP2", charge5nsS_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("charge5nsB_MCP1", charge5nsB_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("charge5nsB_MCP2", charge5nsB_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("TS_toClock_FE_MCP1", TS_toClock_FE_MCP1_loaded[event_nr]); + rd_full->doubleUserRecords.add("TS_toClock_FE_MCP2", TS_toClock_FE_MCP2_loaded[event_nr]); + rd_full->doubleUserRecords.add("meanClockFE", meanClockFE_loaded[event_nr]); + rd_full->doubleUserRecords.add("rmsClockFE", rmsClockFE_loaded[event_nr]); } } (*dwcs)[0] = *dwc1; @@ -258,9 +269,29 @@ void HGCalTBBeamWireChamberProducer::loadRun(int loading_run) { reco3_x_loaded.clear(); reco3_y_loaded.clear(); z3_loaded.clear(); averageHitMultiplicty3_loaded.clear(); reco4_x_loaded.clear(); reco4_y_loaded.clear(); z4_loaded.clear(); averageHitMultiplicty4_loaded.clear(); XCET_021507_signal_loaded.clear(); XCET_021523_signal_loaded.clear(); scintillator_coincidences_loaded.clear(); scintillator_vetos_loaded.clear(); - valid_TS_MCP1_loaded.clear(); valid_TS_MCP2_loaded.clear(); - TS_15PercentRise_MCP1_loaded.clear(); TS_15PercentRise_MCP2_loaded.clear(); TS_30PercentRise_MCP2_loaded.clear(); TS_30PercentRise_MCP1_loaded.clear(); TS_45PercentRise_MCP1_loaded.clear(); TS_45PercentRise_MCP2_loaded.clear(); TS_60PercentRise_MCP1_loaded.clear(); TS_60PercentRise_MCP2_loaded.clear(); amp_MCP1_loaded.clear(); amp_MCP2_loaded.clear(); - TS_MCP1_loaded.clear(); TS_MCP2_loaded.clear(); TS_MCP1_to_last_falling_Edge_loaded.clear(); TS_MCP2_to_last_falling_Edge_loaded.clear(); + + valid_TS_MCP1_loaded.clear(); + valid_TS_MCP2_loaded.clear(); + noise_MCP1_loaded.clear(); + noise_MCP2_loaded.clear(); + TSpeak_MCP1_loaded.clear(); + TSpeak_MCP2_loaded.clear(); + amp_MCP1_loaded.clear(); + amp_MCP2_loaded.clear(); + ampFit_MCP1_loaded.clear(); + ampFit_MCP2_loaded.clear(); + TSfitPeak_MCP1_loaded.clear(); + TSfitPeak_MCP2_loaded.clear(); + TScf_MCP1_loaded.clear(); + TScf_MCP2_loaded.clear(); + charge5nsS_MCP1_loaded.clear(); + charge5nsS_MCP2_loaded.clear(); + charge5nsB_MCP1_loaded.clear(); + charge5nsB_MCP2_loaded.clear(); + TS_toClock_FE_MCP1_loaded.clear(); + TS_toClock_FE_MCP2_loaded.clear(); + meanClockFE_loaded.clear(); + rmsClockFE_loaded.clear(); if (tree == NULL) return; @@ -312,22 +343,28 @@ void HGCalTBBeamWireChamberProducer::loadRun(int loading_run) { scintillator_coincidences_loaded[eventId] = scintillator_coincidences; scintillator_vetos_loaded[eventId] = scintillator_vetos; - valid_TS_MCP1_loaded[eventId] = valid_TS_MCP1; - valid_TS_MCP2_loaded[eventId] = valid_TS_MCP2; - TS_MCP1_loaded[eventId] = TS_MCP1; - TS_MCP2_loaded[eventId] = TS_MCP2; - TS_15PercentRise_MCP1_loaded[eventId] = TS_15PercentRise_MCP1; - TS_15PercentRise_MCP2_loaded[eventId] = TS_15PercentRise_MCP2; - TS_30PercentRise_MCP2_loaded[eventId] = TS_30PercentRise_MCP2; - TS_30PercentRise_MCP1_loaded[eventId] = TS_30PercentRise_MCP1; - TS_45PercentRise_MCP1_loaded[eventId] = TS_45PercentRise_MCP1; - TS_45PercentRise_MCP2_loaded[eventId] = TS_45PercentRise_MCP2; - TS_60PercentRise_MCP1_loaded[eventId] = TS_60PercentRise_MCP1; - TS_60PercentRise_MCP2_loaded[eventId] = TS_60PercentRise_MCP2; - amp_MCP1_loaded[eventId] = amp_MCP1; - amp_MCP2_loaded[eventId] = amp_MCP2; - TS_MCP1_to_last_falling_Edge_loaded[eventId] = TS_MCP1_to_last_falling_Edge; - TS_MCP2_to_last_falling_Edge_loaded[eventId] = TS_MCP2_to_last_falling_Edge; + valid_TS_MCP1_loaded[eventId] = valid_TS_MCP1; + valid_TS_MCP2_loaded[eventId] = valid_TS_MCP2; + noise_MCP1_loaded[eventId] = noise_MCP1; + noise_MCP2_loaded[eventId] = noise_MCP2; + TSpeak_MCP1_loaded[eventId] = TSpeak_MCP1; + TSpeak_MCP2_loaded[eventId] = TSpeak_MCP2; + amp_MCP1_loaded[eventId] = amp_MCP1; + amp_MCP2_loaded[eventId] = amp_MCP2; + ampFit_MCP1_loaded[eventId] = ampFit_MCP1; + ampFit_MCP2_loaded[eventId] = ampFit_MCP2; + TSfitPeak_MCP1_loaded[eventId] = TSfitPeak_MCP1; + TSfitPeak_MCP2_loaded[eventId] = TSfitPeak_MCP2; + TScf_MCP1_loaded[eventId] = TScf_MCP1; + TScf_MCP2_loaded[eventId] = TScf_MCP2; + charge5nsS_MCP1_loaded[eventId] = charge5nsS_MCP1; + charge5nsS_MCP2_loaded[eventId] = charge5nsS_MCP2; + charge5nsB_MCP1_loaded[eventId] = charge5nsB_MCP1; + charge5nsB_MCP2_loaded[eventId] = charge5nsB_MCP2; + TS_toClock_FE_MCP1_loaded[eventId] = TS_toClock_FE_MCP1; + TS_toClock_FE_MCP2_loaded[eventId] = TS_toClock_FE_MCP2; + meanClockFE_loaded[eventId] = meanClockFE; + rmsClockFE_loaded[eventId] = rmsClockFE; } #ifdef DEBUG std::cout << "Loaded run " << loading_run << std::endl; diff --git a/RawToDigi/plugins/HGCalTBCAENProducer.h b/RawToDigi/plugins/HGCalTBCAENProducer.h index 5a30375..1b6245d 100644 --- a/RawToDigi/plugins/HGCalTBCAENProducer.h +++ b/RawToDigi/plugins/HGCalTBCAENProducer.h @@ -42,10 +42,30 @@ class HGCalTBBeamWireChamberProducer : public edm::EDProducer{ short XCET_021507_signal, XCET_021523_signal; short scintillator_coincidences; short scintillator_vetos; - short valid_TS_MCP1, valid_TS_MCP2; - float TS_15PercentRise_MCP1; float TS_15PercentRise_MCP2; float TS_30PercentRise_MCP2; float TS_30PercentRise_MCP1; float TS_45PercentRise_MCP1; float TS_45PercentRise_MCP2; float TS_60PercentRise_MCP1; float TS_60PercentRise_MCP2; float amp_MCP1; float amp_MCP2; - float TS_MCP1, TS_MCP2, TS_MCP1_to_last_falling_Edge, TS_MCP2_to_last_falling_Edge; - + + short valid_TS_MCP1; + short valid_TS_MCP2; + float noise_MCP1; + float noise_MCP2; + float TSpeak_MCP1; + float TSpeak_MCP2; + float amp_MCP1; + float amp_MCP2; + float ampFit_MCP1; + float ampFit_MCP2; + float TSfitPeak_MCP1; + float TSfitPeak_MCP2; + float TScf_MCP1; + float TScf_MCP2; + float charge5nsS_MCP1; + float charge5nsS_MCP2; + float charge5nsB_MCP1; + float charge5nsB_MCP2; + float TS_toClock_FE_MCP1; + float TS_toClock_FE_MCP2; + float meanClockFE; + float rmsClockFE; + TBranch *b_run; TBranch *b_event; TBranch *b_goodDWC_Measurement; @@ -79,23 +99,29 @@ class HGCalTBBeamWireChamberProducer : public edm::EDProducer{ TBranch *b_XCET_021523_signal; TBranch *b_scintillator_coincidences; TBranch *b_scintillator_vetos; + TBranch *b_valid_TS_MCP1; TBranch *b_valid_TS_MCP2; - TBranch *b_TS_MCP1; - TBranch *b_TS_MCP2; - TBranch *b_TS_15PercentRise_MCP1; - TBranch *b_TS_15PercentRise_MCP2; - TBranch *b_TS_30PercentRise_MCP2; - TBranch *b_TS_30PercentRise_MCP1; - TBranch *b_TS_45PercentRise_MCP1; - TBranch *b_TS_45PercentRise_MCP2; - TBranch *b_TS_60PercentRise_MCP1; - TBranch *b_TS_60PercentRise_MCP2; + TBranch *b_noise_MCP1; + TBranch *b_noise_MCP2; + TBranch *b_TSpeak_MCP1; + TBranch *b_TSpeak_MCP2; TBranch *b_amp_MCP1; - TBranch *b_amp_MCP2; - TBranch *b_TS_MCP1_to_last_falling_Edge; - TBranch *b_TS_MCP2_to_last_falling_Edge; - + TBranch *b_amp_MCP2; + TBranch *b_ampFit_MCP1; + TBranch *b_ampFit_MCP2; + TBranch *b_TSfitPeak_MCP1; + TBranch *b_TSfitPeak_MCP2; + TBranch *b_TScf_MCP1; + TBranch *b_TScf_MCP2; + TBranch *b_charge5nsS_MCP1; + TBranch *b_charge5nsS_MCP2; + TBranch *b_charge5nsB_MCP1; + TBranch *b_charge5nsB_MCP2; + TBranch *b_TS_toClock_FE_MCP1; + TBranch *b_TS_toClock_FE_MCP2; + TBranch *b_meanClockFE; + TBranch *b_rmsClockFE; int loaded_run; std::map goodDWC_Measurement_loaded; @@ -136,21 +162,26 @@ class HGCalTBBeamWireChamberProducer : public edm::EDProducer{ std::map valid_TS_MCP1_loaded; std::map valid_TS_MCP2_loaded; - std::map TS_MCP1_loaded; - std::map TS_MCP2_loaded; - std::map TS_15PercentRise_MCP1_loaded; - std::map TS_15PercentRise_MCP2_loaded; - std::map TS_30PercentRise_MCP2_loaded; - std::map TS_30PercentRise_MCP1_loaded; - std::map TS_45PercentRise_MCP1_loaded; - std::map TS_45PercentRise_MCP2_loaded; - std::map TS_60PercentRise_MCP1_loaded; - std::map TS_60PercentRise_MCP2_loaded; + std::map noise_MCP1_loaded; + std::map noise_MCP2_loaded; + std::map TSpeak_MCP1_loaded; + std::map TSpeak_MCP2_loaded; std::map amp_MCP1_loaded; std::map amp_MCP2_loaded; - - std::map TS_MCP1_to_last_falling_Edge_loaded; - std::map TS_MCP2_to_last_falling_Edge_loaded; + std::map ampFit_MCP1_loaded; + std::map ampFit_MCP2_loaded; + std::map TSfitPeak_MCP1_loaded; + std::map TSfitPeak_MCP2_loaded; + std::map TScf_MCP1_loaded; + std::map TScf_MCP2_loaded; + std::map charge5nsS_MCP1_loaded; + std::map charge5nsS_MCP2_loaded; + std::map charge5nsB_MCP1_loaded; + std::map charge5nsB_MCP2_loaded; + std::map TS_toClock_FE_MCP1_loaded; + std::map TS_toClock_FE_MCP2_loaded; + std::map meanClockFE_loaded; + std::map rmsClockFE_loaded; void loadRun(int loading_run); }; diff --git a/RawToDigi/plugins/HGCalTBCAENSource.cc b/RawToDigi/plugins/HGCalTBCAENSource.cc index dfdd45a..5731033 100644 --- a/RawToDigi/plugins/HGCalTBCAENSource.cc +++ b/RawToDigi/plugins/HGCalTBCAENSource.cc @@ -557,15 +557,16 @@ void HGCalTBWireChamberSource::produce(edm::Event & event) { pedestal_ch7[i] = digi_samples.at(7)->at(i); } + MCPwaveform mcpX; //determine and subtract the baseline for all samples - substractBaseline(N_digi_samples, MCP1_waveform, getBaseline(findAbsolutePeak(N_digi_samples, MCP1_waveform, "pos"), MCP1_waveform, 90, -1)); - substractBaseline(N_digi_samples, MCP2_waveform, getBaseline(findAbsolutePeak(N_digi_samples, MCP2_waveform, "pos"), MCP2_waveform, 90, -1)); - substractBaseline(N_digi_samples, pedestal_ch4, getBaseline(findAbsolutePeak(N_digi_samples, pedestal_ch4, "pos"), pedestal_ch4, 90, -1)); - substractBaseline(N_digi_samples, pedestal_ch5, getBaseline(findAbsolutePeak(N_digi_samples, pedestal_ch5, "pos"), pedestal_ch5, 90, -1)); - substractBaseline(N_digi_samples, pedestal_ch6, getBaseline(findAbsolutePeak(N_digi_samples, pedestal_ch6, "pos"), pedestal_ch6, 90, -1)); - substractBaseline(N_digi_samples, pedestal_ch7, getBaseline(findAbsolutePeak(N_digi_samples, pedestal_ch7, "pos"), pedestal_ch7, 90, -1)); - + //range chosen is 5 - 50 + mcpX.substractBaseline(N_digi_samples, MCP1_waveform); + mcpX.substractBaseline(N_digi_samples, MCP2_waveform); + mcpX.substractBaseline(N_digi_samples, pedestal_ch4); + mcpX.substractBaseline(N_digi_samples, pedestal_ch5); + mcpX.substractBaseline(N_digi_samples, pedestal_ch6); + mcpX.substractBaseline(N_digi_samples, pedestal_ch7); //step1: obtain common mode noise from inactive channels short* commonBaseline = new short[N_digi_samples]; @@ -583,62 +584,118 @@ void HGCalTBWireChamberSource::produce(edm::Event & event) { short* MCP2_waveform_cleared = new short[N_digi_samples]; for (int i = 0; i < N_digi_samples; i++) MCP1_waveform_cleared[i] = +commonBaseline[i] - MCP1_waveform[i]; for (int i = 0; i < N_digi_samples; i++) MCP2_waveform_cleared[i] = +commonBaseline[i] - MCP2_waveform[i]; + + //3. apply pol2 fits around the maximum (+/- 5 sample) + // and the linear fit to the rising edge (with 5 points = 1ns = about the rise time) + MCPwaveform* mcp1signl = new MCPwaveform(0, N_digi_samples, MCP1_waveform_cleared, "pos"); + MCPwaveform* mcp2signl = new MCPwaveform(0, N_digi_samples, MCP2_waveform_cleared, "pos"); +#ifdef DEBUG + std::cout << " mcp1 baseline = " << mcp1signl->getBaseline(90, -1) << std::endl; + std::cout << " mcp2 baseline = " << mcp2signl->getBaseline(90, -1) << std::endl; +#endif + mcp1signl->analysePeak(); + mcp2signl->analysePeak(); - //3. apply gaussian fits around the maximum (+/- 1 sample) and the linear fit to the rising edge - peakValues* MCPSignal1 = analysePeak(N_digi_samples, MCP1_waveform_cleared); - peakValues* MCPSignal2 = analysePeak(N_digi_samples, MCP2_waveform_cleared); #ifdef DEBUG - std::cout << "MCP1: " << MCPSignal1->fQuality << " " << MCPSignal1->peak << " " << MCPSignal1->amp << " " << MCPSignal1->amppeak << " " << MCPSignal1->tpeak << " " << MCPSignal1->base << " " << std::endl; - std::cout << "MCP2: " << MCPSignal2->fQuality << " " << MCPSignal2->peak << " " << MCPSignal2->amp << " " << MCPSignal2->amppeak << " " << MCPSignal2->tpeak << " " << MCPSignal2->base << " " << std::endl; + std::cout << "MCP1: " << mcp1signl->getQuality() << " fitAmp " << mcp1signl->getFitAmp() << " amp " << mcp1signl->getAmp() + << " peak = " << mcp1signl->getPeak() << " fit peak = " << mcp1signl->getFitPeak() + << " CF " << mcp1signl->getTimeCF(0.5) << " chargeS " << mcp1signl->getCharge5nsS() + << " chargeB " << mcp1signl->getCharge5nsB() << " " << std::endl; + + std::cout << "MCP2: " << mcp2signl->getQuality() << " fitAmp " << mcp2signl->getFitAmp() << " amp " << mcp2signl->getAmp() + << " peak = " << mcp2signl->getPeak() << " fit peak = " << mcp2signl->getFitPeak() + << " CF " << mcp2signl->getTimeCF(0.5) << " chargeS " << mcp2signl->getCharge5nsS() + << " chargeB " << mcp2signl->getCharge5nsB() << " " << std::endl; #endif - - //4. determine distance to the prior falling clock edge - float priorFallingClockEdge_MCP1 = 0; - for (int sample = 2; sample < N_digi_samples - 1; sample++) { - if (sample > MCPSignal1->tpeak) break; - if ((digi_clock[sample - 2] == 4095) && (digi_clock[sample - 1] > digi_clock[sample]) && (digi_clock[sample] > digi_clock[sample + 1])) { - float under3200 = sample + (3200. - digi_clock[sample]) / (digi_clock[sample + 1] - digi_clock[sample]); -#ifdef DEBUG - std::cout << digi_clock[sample - 1] << " " << sample << "," << digi_clock[sample] << " " << sample + 1 << "," << digi_clock[sample + 1] << ": " << under3200 << std::endl; -#endif - if (sample < MCPSignal1->tpeak) priorFallingClockEdge_MCP1 = under3200; - } - } - float priorFallingClockEdge_MCP2 = 0; - for (int sample = 2; sample < N_digi_samples - 1; sample++) { - if (sample > MCPSignal2->tpeak) break; - if ((digi_clock[sample - 2] == 4095) && (digi_clock[sample - 1] > digi_clock[sample]) && (digi_clock[sample] > digi_clock[sample + 1])) { - float under3200 = sample + (3200. - digi_clock[sample]) / (digi_clock[sample + 1] - digi_clock[sample]); -#ifdef DEBUG - std::cout << digi_clock[sample - 1] << " " << sample << "," << digi_clock[sample] << " " << sample + 1 << "," << digi_clock[sample + 1] << ": " << under3200 << std::endl; -#endif - if (sample < MCPSignal2->tpeak) priorFallingClockEdge_MCP2 = under3200; - } + //timeCF available only for goodQuality waveform (-1 otherwise) + //loose goodQuality criteria (amp > 3.*noise && amp > 100.) + //not to loose acceptance => need study and tuning for analysis + float MCP1sampleTimeCF = mcp1signl->getTimeCF(0.5); + float MCP2sampleTimeCF = mcp2signl->getTimeCF(0.5); + + + //4. extract time stamp from synch board clock + // take mean from 8 falling edges and save also rms + // assume y = A + B x + // B = (NSxy Sy - Sx Sy)/(NSxx - SxSx) + // A = Sy/N - B (Sx/N) = (Sxx Sy - Sx Sxy)/(NSxx - SxSx) + std::vector tsClock_FE_at3200; + for(int startSample = 2; startSample < 1000; startSample += 125){ + float xx = 0.; + float xy = 0.; + float Sx = 0.; + float Sy = 0.; + float Sxx = 0.; + float Sxy = 0.; + int usedSamples=0; + for (int sample = startSample; sample < N_digi_samples; ++sample){ + if(digi_clock[sample] == 4095 && digi_clock[sample+1] < 4095 && digi_clock[sample+2] < 4095){ + for(int iSample=sample+1; iSamplebooleanUserRecords.add("valid_TS_MCP1", MCPSignal1->fQuality); - rd->booleanUserRecords.add("valid_TS_MCP2", MCPSignal2->fQuality); - rd->doubleUserRecords.add("TS_MCP1", 0.2 * MCPSignal1->tpeak); - rd->doubleUserRecords.add("TS_MCP2", 0.2 * MCPSignal2->tpeak); - rd->doubleUserRecords.add("TS_15PercentRise_MCP1", 0.2 * MCPSignal1->tlinear15); - rd->doubleUserRecords.add("TS_15PercentRise_MCP2", 0.2 * MCPSignal2->tlinear15); - rd->doubleUserRecords.add("TS_30PercentRise_MCP1", 0.2 * MCPSignal1->tlinear30); - rd->doubleUserRecords.add("TS_30PercentRise_MCP2", 0.2 * MCPSignal2->tlinear30); - rd->doubleUserRecords.add("TS_45PercentRise_MCP1", 0.2 * MCPSignal1->tlinear45); - rd->doubleUserRecords.add("TS_45PercentRise_MCP2", 0.2 * MCPSignal2->tlinear45); - rd->doubleUserRecords.add("TS_60PercentRise_MCP1", 0.2 * MCPSignal1->tlinear60); - rd->doubleUserRecords.add("TS_60PercentRise_MCP2", 0.2 * MCPSignal2->tlinear60); - rd->doubleUserRecords.add("amp_MCP1", MCPSignal1->amppeak); - rd->doubleUserRecords.add("amp_MCP2", MCPSignal2->amppeak); - rd->doubleUserRecords.add("TS_MCP1_to_last_falling_Edge", 0.2 * (MCPSignal1->tpeak - priorFallingClockEdge_MCP1)); - rd->doubleUserRecords.add("TS_MCP2_to_last_falling_Edge", 0.2 * (MCPSignal2->tpeak - priorFallingClockEdge_MCP2)); + rd->intUserRecords.add("valid_TS_MCP1", mcp1signl->getQuality()); + rd->intUserRecords.add("valid_TS_MCP2", mcp2signl->getQuality()); + rd->doubleUserRecords.add("noise_MCP1", mcp1signl->getNoise()); + rd->doubleUserRecords.add("noise_MCP2", mcp2signl->getNoise()); + rd->doubleUserRecords.add("TSpeak_MCP1", mcp1signl->getPeak()*0.2); + rd->doubleUserRecords.add("TSpeak_MCP2", mcp2signl->getPeak()*0.2); + rd->doubleUserRecords.add("amp_MCP1", mcp1signl->getAmp()); + rd->doubleUserRecords.add("amp_MCP2", mcp2signl->getAmp()); + rd->doubleUserRecords.add("ampFit_MCP1", mcp1signl->getFitAmp()); + rd->doubleUserRecords.add("ampFit_MCP2", mcp2signl->getFitAmp()); + rd->doubleUserRecords.add("TSfitPeak_MCP1", mcp1signl->getFitPeak()*0.2); + rd->doubleUserRecords.add("TSfitPeak_MCP2", mcp2signl->getFitPeak()*0.2); + rd->doubleUserRecords.add("TScf_MCP1", mcp1signl->getTimeCF(0.5)*0.2); + rd->doubleUserRecords.add("TScf_MCP2", mcp2signl->getTimeCF(0.5)*0.2); + + rd->doubleUserRecords.add("charge5nsS_MCP1", mcp1signl->getCharge5nsS()); + rd->doubleUserRecords.add("charge5nsS_MCP2", mcp2signl->getCharge5nsS()); + rd->doubleUserRecords.add("charge5nsB_MCP1", mcp1signl->getCharge5nsB()); + rd->doubleUserRecords.add("charge5nsB_MCP2", mcp2signl->getCharge5nsB()); + + rd->doubleUserRecords.add("TS_toClock_FE_MCP1", TS_fromMeanClockFE_MCP1*0.2); + rd->doubleUserRecords.add("TS_toClock_FE_MCP2", TS_fromMeanClockFE_MCP2*0.2); + rd->doubleUserRecords.add("meanClockFE", meanClockFE_at3200); + rd->doubleUserRecords.add("rmsClockFE", rmsClockFE_at3200); + + delete mcp1signl; + delete mcp2signl; } } diff --git a/RawToDigi/plugins/HGCalTBCAENSource.h b/RawToDigi/plugins/HGCalTBCAENSource.h index cb2f2a9..14fa814 100644 --- a/RawToDigi/plugins/HGCalTBCAENSource.h +++ b/RawToDigi/plugins/HGCalTBCAENSource.h @@ -20,7 +20,7 @@ #include "TBranch.h" #include "TTree.h" -#include "HGCal/Reco/interface/MCP_waveform_reco.h" +#include "HGCal/Reco/interface/MCPwaveform.h" //DWC_to_TDC_MAP int DWC1_LEFT; diff --git a/Reco/interface/MCPwaveform.h b/Reco/interface/MCPwaveform.h new file mode 100644 index 0000000..24b05ff --- /dev/null +++ b/Reco/interface/MCPwaveform.h @@ -0,0 +1,95 @@ +#ifndef MCPwaveform_h +#define MCPwaveform_h + +#include +#include +#include +#include +#include +#include +#include + +#include "TTree.h" +#include "TBranch.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TFile.h" +#include "TTree.h" +#include "TBranch.h" +#include "TCanvas.h" +#include "TMath.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TF1.h" +#include "TLinearFitter.h" +#include "TGraph.h" +#include "TGraphErrors.h" + + + +class MCPwaveform{ + public: + + MCPwaveform(); + MCPwaveform(int nsamples, short* pulse, std::string); + MCPwaveform(int verbosity, int nsamples, short* pulse, std::string); + + void analysePeak(); + + int findAbsolutePeak(std::string polarity="pos"); + + void substractBaseline(int nSamples, short* pulse); + float getBaseline(int nbinsExcludedLeftOfPeak, int nbinsExcludedRightOfPeak); + + float getIntegral(int min, int max); + float getSignalIntegral(int riseWin, int fallWin); + float getModIntegral(int min, int max); + + float getTimeCF(float frac, int nFitSamples=0, int min=0, int max=0); + float linearInterpolation(const int& min, const int& max, const int& sampleToSkip=0); + float pol2Interpolation(const int& min, const int& max, const int& sampleToSkip=0); + float getInterpolatedAmpMax(int min, int max); + + + inline float getNoise(){ return noise_; }; + inline int getPeak(){ return peak_; }; + inline float getAmp(){ return amp_; }; + inline float getFitAmp(){ return fitAmp_; }; + inline float getFitPeak(){ return fitTimeMax_; }; + inline int getQuality(){ return fQuality_; }; + inline float getCharge5nsS(){ return charge5nsS_; }; + inline float getCharge5nsB(){ return charge5nsB_; }; + + + + private: + + int verbosity_; + std::vector wf_; + int nSamples_; + std::string polarity_; + int sigWinMin_; + int sigWinMax_; + + Int_t peak_; + Float_t amp_; + Float_t base_; + Float_t noise_; + + Float_t charge5nsS_; + Float_t charge5nsB_; + + UInt_t fQuality_; + Float_t Aparam_; + Float_t Bparam_; + Float_t fitAmp_; + Float_t fitTimeMax_; + Float_t fitTimeCF_; + + int nFitSamples_; + int minForAmpFit_; + int maxForAmpFit_; +}; + +#endif diff --git a/Reco/plugins/CAEN_NTupelizer.cc b/Reco/plugins/CAEN_NTupelizer.cc index e22213b..b3c0089 100644 --- a/Reco/plugins/CAEN_NTupelizer.cc +++ b/Reco/plugins/CAEN_NTupelizer.cc @@ -83,26 +83,37 @@ class DWC_NTupelizer : public edm::one::EDAnalyzer { std::vector digi_clock; std::vector digi_MCP1; std::vector digi_MCP2; + std::vector digi_MCP4; + std::vector digi_MCP5; + std::vector digi_MCP6; + std::vector digi_MCP7; + std::vector digi_MCP1_cleared; + std::vector digi_MCP2_cleared; std::vector digi_scintillator_big; std::vector digi_synchboard_trigger; - short valid_TS_MCP1; - short valid_TS_MCP2; - float TS_MCP1; - float TS_MCP2; - float TS_15PercentRise_MCP1; - float TS_15PercentRise_MCP2; - float TS_30PercentRise_MCP2; - float TS_30PercentRise_MCP1; - float TS_45PercentRise_MCP1; - float TS_45PercentRise_MCP2; - float TS_60PercentRise_MCP1; - float TS_60PercentRise_MCP2; - float amp_MCP1; - float amp_MCP2; - float TS_MCP1_to_last_falling_Edge; - float TS_MCP2_to_last_falling_Edge; - + short valid_TS_MCP1; + short valid_TS_MCP2; + float noise_MCP1; + float noise_MCP2; + float TSpeak_MCP1; + float TSpeak_MCP2; + float amp_MCP1; + float amp_MCP2; + float ampFit_MCP1; + float ampFit_MCP2; + float TSfitPeak_MCP1; + float TSfitPeak_MCP2; + float TScf_MCP1; + float TScf_MCP2; + float charge5nsS_MCP1; + float charge5nsS_MCP2; + float charge5nsB_MCP1; + float charge5nsB_MCP2; + float TS_toClock_FE_MCP1; + float TS_toClock_FE_MCP2; + float meanClockFE; + float rmsClockFE; std::map Sensors; ParticleTrack* TrackFull; @@ -261,28 +272,40 @@ void DWC_NTupelizer::analyze(const edm::Event& event, const edm::EventSetup& set if (rd->shortVectorUserRecords.has("digi_clock")) digi_clock = rd->shortVectorUserRecords.get("digi_clock"); if (rd->shortVectorUserRecords.has("digi_MCP1")) digi_MCP1 = rd->shortVectorUserRecords.get("digi_MCP1"); if (rd->shortVectorUserRecords.has("digi_MCP2")) digi_MCP2 = rd->shortVectorUserRecords.get("digi_MCP2"); + if (rd->shortVectorUserRecords.has("digi_MCP4")) digi_MCP4 = rd->shortVectorUserRecords.get("digi_MCP4"); + if (rd->shortVectorUserRecords.has("digi_MCP5")) digi_MCP5 = rd->shortVectorUserRecords.get("digi_MCP5"); + if (rd->shortVectorUserRecords.has("digi_MCP6")) digi_MCP6 = rd->shortVectorUserRecords.get("digi_MCP6"); + if (rd->shortVectorUserRecords.has("digi_MCP6")) digi_MCP7 = rd->shortVectorUserRecords.get("digi_MCP7"); + if (rd->shortVectorUserRecords.has("digi_MCP1_cleared")) digi_MCP1_cleared = rd->shortVectorUserRecords.get("digi_MCP1_cleared"); + if (rd->shortVectorUserRecords.has("digi_MCP2_cleared")) digi_MCP2_cleared = rd->shortVectorUserRecords.get("digi_MCP2_cleared"); if (rd->shortVectorUserRecords.has("digi_scintillator_4x4")) digi_scintillator_big = rd->shortVectorUserRecords.get("digi_scintillator_4x4"); if (rd->shortVectorUserRecords.has("digi_synchboard_trigger")) digi_synchboard_trigger = rd->shortVectorUserRecords.get("digi_synchboard_trigger"); - if(rd->booleanUserRecords.has("valid_TS_MCP1")) valid_TS_MCP1 = rd->booleanUserRecords.get("valid_TS_MCP1"); else valid_TS_MCP1 = -999; - if(rd->booleanUserRecords.has("valid_TS_MCP2")) valid_TS_MCP2 = rd->booleanUserRecords.get("valid_TS_MCP2"); else valid_TS_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_MCP1")) TS_MCP1 = rd->doubleUserRecords.get("TS_MCP1"); else TS_MCP1 = -999; - if(rd->doubleUserRecords.has("TS_MCP2")) TS_MCP2 = rd->doubleUserRecords.get("TS_MCP2"); else TS_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_15PercentRise_MCP1")) TS_15PercentRise_MCP1 = rd->doubleUserRecords.get("TS_15PercentRise_MCP1"); else TS_15PercentRise_MCP1 = -999; - if(rd->doubleUserRecords.has("TS_15PercentRise_MCP2")) TS_15PercentRise_MCP2 = rd->doubleUserRecords.get("TS_15PercentRise_MCP2"); else TS_15PercentRise_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_30PercentRise_MCP2")) TS_30PercentRise_MCP2 = rd->doubleUserRecords.get("TS_30PercentRise_MCP2"); else TS_30PercentRise_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_30PercentRise_MCP1")) TS_30PercentRise_MCP1 = rd->doubleUserRecords.get("TS_30PercentRise_MCP1"); else TS_30PercentRise_MCP1 = -999; - if(rd->doubleUserRecords.has("TS_45PercentRise_MCP1")) TS_45PercentRise_MCP1 = rd->doubleUserRecords.get("TS_45PercentRise_MCP1"); else TS_45PercentRise_MCP1 = -999; - if(rd->doubleUserRecords.has("TS_45PercentRise_MCP2")) TS_45PercentRise_MCP2 = rd->doubleUserRecords.get("TS_45PercentRise_MCP2"); else TS_45PercentRise_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_60PercentRise_MCP1")) TS_60PercentRise_MCP1 = rd->doubleUserRecords.get("TS_60PercentRise_MCP1"); else TS_60PercentRise_MCP1 = -999; - if(rd->doubleUserRecords.has("TS_60PercentRise_MCP2")) TS_60PercentRise_MCP2 = rd->doubleUserRecords.get("TS_60PercentRise_MCP2"); else TS_60PercentRise_MCP2 = -999; - if(rd->doubleUserRecords.has("amp_MCP1")) amp_MCP1 = rd->doubleUserRecords.get("amp_MCP1"); else amp_MCP1 = -999; - if(rd->doubleUserRecords.has("amp_MCP2")) amp_MCP2 = rd->doubleUserRecords.get("amp_MCP2"); else amp_MCP2 = -999; - if(rd->doubleUserRecords.has("TS_MCP1_to_last_falling_Edge")) TS_MCP1_to_last_falling_Edge = rd->doubleUserRecords.get("TS_MCP1_to_last_falling_Edge"); else TS_MCP1_to_last_falling_Edge = -999; - if(rd->doubleUserRecords.has("TS_MCP2_to_last_falling_Edge")) TS_MCP2_to_last_falling_Edge = rd->doubleUserRecords.get("TS_MCP2_to_last_falling_Edge"); else TS_MCP2_to_last_falling_Edge = -999; - tree->Fill(); + valid_TS_MCP1 = (rd->intUserRecords.has("valid_TS_MCP1")) ? rd->intUserRecords.get("valid_TS_MCP1") : -999; + valid_TS_MCP2 = (rd->intUserRecords.has("valid_TS_MCP2")) ? rd->intUserRecords.get("valid_TS_MCP2") : -999; + noise_MCP1 = (rd->doubleUserRecords.has("noise_MCP1")) ? rd->doubleUserRecords.get("noise_MCP1") : -999; + noise_MCP2 = (rd->doubleUserRecords.has("noise_MCP2")) ? rd->doubleUserRecords.get("noise_MCP2") : -999; + TSpeak_MCP1 = (rd->doubleUserRecords.has("TSpeak_MCP1")) ? rd->doubleUserRecords.get("TSpeak_MCP1") : -999; + TSpeak_MCP2 = (rd->doubleUserRecords.has("TSpeak_MCP2")) ? rd->doubleUserRecords.get("TSpeak_MCP2") : -999; + amp_MCP1 = (rd->doubleUserRecords.has("amp_MCP1")) ? rd->doubleUserRecords.get("amp_MCP1") : -999; + amp_MCP2 = (rd->doubleUserRecords.has("amp_MCP2")) ? rd->doubleUserRecords.get("amp_MCP2") : -999; + ampFit_MCP1 = (rd->doubleUserRecords.has("ampFit_MCP1")) ? rd->doubleUserRecords.get("ampFit_MCP1") : -999; + ampFit_MCP2 = (rd->doubleUserRecords.has("ampFit_MCP2")) ? rd->doubleUserRecords.get("ampFit_MCP2") : -999; + TSfitPeak_MCP1 = (rd->doubleUserRecords.has("TSfitPeak_MCP1")) ? rd->doubleUserRecords.get("TSfitPeak_MCP1") : -999; + TSfitPeak_MCP2 = (rd->doubleUserRecords.has("TSfitPeak_MCP2")) ? rd->doubleUserRecords.get("TSfitPeak_MCP2") : -999; + TScf_MCP1 = (rd->doubleUserRecords.has("TScf_MCP1")) ? rd->doubleUserRecords.get("TScf_MCP1") : -999; + TScf_MCP2 = (rd->doubleUserRecords.has("TScf_MCP2")) ? rd->doubleUserRecords.get("TScf_MCP2") : -999; + charge5nsS_MCP1 = (rd->doubleUserRecords.has("charge5nsS_MCP1")) ? rd->doubleUserRecords.get("charge5nsS_MCP1") : -999; + charge5nsS_MCP2 = (rd->doubleUserRecords.has("charge5nsS_MCP2")) ? rd->doubleUserRecords.get("charge5nsS_MCP2") : -999; + charge5nsB_MCP1 = (rd->doubleUserRecords.has("charge5nsB_MCP1")) ? rd->doubleUserRecords.get("charge5nsB_MCP1") : -999; + charge5nsB_MCP2 = (rd->doubleUserRecords.has("charge5nsB_MCP2")) ? rd->doubleUserRecords.get("charge5nsB_MCP2") : -999; + TS_toClock_FE_MCP1 = (rd->doubleUserRecords.has("TS_toClock_FE_MCP1")) ? rd->doubleUserRecords.get("TS_toClock_FE_MCP1") : -999; + TS_toClock_FE_MCP2 = (rd->doubleUserRecords.has("TS_toClock_FE_MCP2")) ? rd->doubleUserRecords.get("TS_toClock_FE_MCP2") : -999; + meanClockFE = (rd->doubleUserRecords.has("meanClockFE")) ? (rd->doubleUserRecords.get("meanClockFE")) : -999; + rmsClockFE = (rd->doubleUserRecords.has("rmsClockFE")) ? (rd->doubleUserRecords.get("rmsClockFE")) : -999; + tree->Fill(); }// analyze ends here @@ -392,27 +415,37 @@ void DWC_NTupelizer::beginJob() { tree->Branch("digi_clock", &digi_clock); tree->Branch("digi_MCP1", &digi_MCP1); tree->Branch("digi_MCP2", &digi_MCP2); + tree->Branch("digi_MCP4", &digi_MCP4); + tree->Branch("digi_MCP5", &digi_MCP5); + tree->Branch("digi_MCP6", &digi_MCP6); + tree->Branch("digi_MCP7", &digi_MCP7); + tree->Branch("digi_MCP1_cleared", &digi_MCP1_cleared); + tree->Branch("digi_MCP2_cleared", &digi_MCP2_cleared); tree->Branch("digi_scintillator_4x4", &digi_scintillator_big); tree->Branch("digi_synchboard_trigger", &digi_synchboard_trigger); - tree->Branch("valid_TS_MCP1", &valid_TS_MCP1); tree->Branch("valid_TS_MCP2", &valid_TS_MCP2); - tree->Branch("TS_MCP1", &TS_MCP1); - tree->Branch("TS_MCP2", &TS_MCP2); - tree->Branch("TS_15PercentRise_MCP1", &TS_15PercentRise_MCP1); - tree->Branch("TS_15PercentRise_MCP2", &TS_15PercentRise_MCP2); - tree->Branch("TS_30PercentRise_MCP2", &TS_30PercentRise_MCP2); - tree->Branch("TS_30PercentRise_MCP1", &TS_30PercentRise_MCP1); - tree->Branch("TS_45PercentRise_MCP1", &TS_45PercentRise_MCP1); - tree->Branch("TS_45PercentRise_MCP2", &TS_45PercentRise_MCP2); - tree->Branch("TS_60PercentRise_MCP1", &TS_60PercentRise_MCP1); - tree->Branch("TS_60PercentRise_MCP2", &TS_60PercentRise_MCP2); - tree->Branch("amp_MCP1", &_MCP1); - tree->Branch("amp_MCP2", &_MCP2); - tree->Branch("TS_MCP1_to_last_falling_Edge", &TS_MCP1_to_last_falling_Edge); - tree->Branch("TS_MCP2_to_last_falling_Edge", &TS_MCP2_to_last_falling_Edge); - + tree->Branch("noise_MCP1", &noise_MCP1); + tree->Branch("noise_MCP2", &noise_MCP2); + tree->Branch("TSpeak_MCP1", &TSpeak_MCP1); + tree->Branch("TSpeak_MCP2", &TSpeak_MCP2); + tree->Branch("amp_MCP1", &_MCP1); + tree->Branch("amp_MCP2", &_MCP2); + tree->Branch("ampFit_MCP1", &Fit_MCP1); + tree->Branch("ampFit_MCP2", &Fit_MCP2); + tree->Branch("TSfitPeak_MCP1", &TSfitPeak_MCP1); + tree->Branch("TSfitPeak_MCP2", &TSfitPeak_MCP2); + tree->Branch("TScf_MCP1", &TScf_MCP1); + tree->Branch("TScf_MCP2", &TScf_MCP2); + tree->Branch("charge5nsS_MCP1", &charge5nsS_MCP1); + tree->Branch("charge5nsS_MCP2", &charge5nsS_MCP2); + tree->Branch("charge5nsB_MCP1", &charge5nsB_MCP1); + tree->Branch("charge5nsB_MCP2", &charge5nsB_MCP2); + tree->Branch("TS_toClock_FE_MCP1", &TS_toClock_FE_MCP1); + tree->Branch("TS_toClock_FE_MCP2", &TS_toClock_FE_MCP2); + tree->Branch("meanClockFE", &meanClockFE); + tree->Branch("rmsClockFE", &rmsClockFE); } } @@ -427,4 +460,4 @@ void DWC_NTupelizer::fillDescriptions(edm::ConfigurationDescriptions& descriptio } //define this as a plug-in -DEFINE_FWK_MODULE(DWC_NTupelizer); \ No newline at end of file +DEFINE_FWK_MODULE(DWC_NTupelizer); diff --git a/Reco/plugins/MCPNtupler.cc b/Reco/plugins/MCPNtupler.cc index dfc7261..8d91b01 100644 --- a/Reco/plugins/MCPNtupler.cc +++ b/Reco/plugins/MCPNtupler.cc @@ -50,21 +50,26 @@ class MCPNtupler : public edm::one::EDAnalyzer short valid_TS_MCP1; short valid_TS_MCP2; - float TS_MCP1; - float TS_MCP2; - float TS_15PercentRise_MCP1; - float TS_15PercentRise_MCP2; - float TS_30PercentRise_MCP2; - float TS_30PercentRise_MCP1; - float TS_45PercentRise_MCP1; - float TS_45PercentRise_MCP2; - float TS_60PercentRise_MCP1; - float TS_60PercentRise_MCP2; + float noise_MCP1; + float noise_MCP2; + float TSpeak_MCP1; + float TSpeak_MCP2; float amp_MCP1; - float amp_MCP2; - float TS_MCP1_to_last_falling_Edge; - float TS_MCP2_to_last_falling_Edge; - + float amp_MCP2; + float ampFit_MCP1; + float ampFit_MCP2; + float TSfitPeak_MCP1; + float TSfitPeak_MCP2; + float TScf_MCP1; + float TScf_MCP2; + float charge5nsS_MCP1; + float charge5nsS_MCP2; + float charge5nsB_MCP1; + float charge5nsB_MCP2; + float TS_toClock_FE_MCP1; + float TS_toClock_FE_MCP2; + float meanClockFE; + float rmsClockFE; }; void MCPNtupler::clearVariables() { @@ -87,24 +92,29 @@ MCPNtupler::MCPNtupler(const edm::ParameterSet& iConfig) // event info tree_->Branch("event", &ev_event_); tree_->Branch("run", &ev_run_); + tree_->Branch("valid_TS_MCP1", &valid_TS_MCP1); tree_->Branch("valid_TS_MCP2", &valid_TS_MCP2); - tree_->Branch("TS_MCP1", &TS_MCP1); - tree_->Branch("TS_MCP2", &TS_MCP2); - tree_->Branch("TS_15PercentRise_MCP1", &TS_15PercentRise_MCP1); - tree_->Branch("TS_15PercentRise_MCP2", &TS_15PercentRise_MCP2); - tree_->Branch("TS_30PercentRise_MCP2", &TS_30PercentRise_MCP2); - tree_->Branch("TS_30PercentRise_MCP1", &TS_30PercentRise_MCP1); - tree_->Branch("TS_45PercentRise_MCP1", &TS_45PercentRise_MCP1); - tree_->Branch("TS_45PercentRise_MCP2", &TS_45PercentRise_MCP2); - tree_->Branch("TS_60PercentRise_MCP1", &TS_60PercentRise_MCP1); - tree_->Branch("TS_60PercentRise_MCP2", &TS_60PercentRise_MCP2); + tree_->Branch("noise_MCP1", &noise_MCP1); + tree_->Branch("noise_MCP2", &noise_MCP2); + tree_->Branch("TSpeak_MCP1", &TSpeak_MCP1); + tree_->Branch("TSpeak_MCP2", &TSpeak_MCP2); tree_->Branch("amp_MCP1", &_MCP1); - tree_->Branch("amp_MCP2", &_MCP2); - - tree_->Branch("TS_MCP1_to_last_falling_Edge", &TS_MCP1_to_last_falling_Edge); - tree_->Branch("TS_MCP2_to_last_falling_Edge", &TS_MCP2_to_last_falling_Edge); - + tree_->Branch("amp_MCP2", &_MCP2); + tree_->Branch("ampFit_MCP1", &Fit_MCP1); + tree_->Branch("ampFit_MCP2", &Fit_MCP2); + tree_->Branch("TSfitPeak_MCP1", &TSfitPeak_MCP1); + tree_->Branch("TSfitPeak_MCP2", &TSfitPeak_MCP2); + tree_->Branch("TScf_MCP1", &TScf_MCP1); + tree_->Branch("TScf_MCP2", &TScf_MCP2); + tree_->Branch("charge5nsS_MCP1", &charge5nsS_MCP1); + tree_->Branch("charge5nsS_MCP2", &charge5nsS_MCP2); + tree_->Branch("charge5nsB_MCP1", &charge5nsB_MCP1); + tree_->Branch("charge5nsB_MCP2", &charge5nsB_MCP2); + tree_->Branch("TS_toClock_FE_MCP1", &TS_toClock_FE_MCP1); + tree_->Branch("TS_toClock_FE_MCP2", &TS_toClock_FE_MCP2); + tree_->Branch("meanClockFE", &meanClockFE); + tree_->Branch("rmsClockFE", &rmsClockFE); } @@ -126,22 +136,29 @@ void MCPNtupler::analyze(const edm::Event& event, const edm::EventSetup& setup) ev_run_ = rd->run; ev_event_ = rd->event; - if (rd->intUserRecords.has("valid_TS_MCP1")) valid_TS_MCP1 = rd->intUserRecords.get("valid_TS_MCP1"); else valid_TS_MCP1 = -999; - if (rd->intUserRecords.has("valid_TS_MCP2")) valid_TS_MCP2 = rd->intUserRecords.get("valid_TS_MCP2"); else valid_TS_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_15PercentRise_MCP1")) TS_15PercentRise_MCP1 = rd->doubleUserRecords.get("TS_15PercentRise_MCP1"); else TS_15PercentRise_MCP1 = -999; - if (rd->doubleUserRecords.has("TS_15PercentRise_MCP2")) TS_15PercentRise_MCP2 = rd->doubleUserRecords.get("TS_15PercentRise_MCP2"); else TS_15PercentRise_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_30PercentRise_MCP2")) TS_30PercentRise_MCP2 = rd->doubleUserRecords.get("TS_30PercentRise_MCP2"); else TS_30PercentRise_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_30PercentRise_MCP1")) TS_30PercentRise_MCP1 = rd->doubleUserRecords.get("TS_30PercentRise_MCP1"); else TS_30PercentRise_MCP1 = -999; - if (rd->doubleUserRecords.has("TS_45PercentRise_MCP1")) TS_45PercentRise_MCP1 = rd->doubleUserRecords.get("TS_45PercentRise_MCP1"); else TS_45PercentRise_MCP1 = -999; - if (rd->doubleUserRecords.has("TS_45PercentRise_MCP2")) TS_45PercentRise_MCP2 = rd->doubleUserRecords.get("TS_45PercentRise_MCP2"); else TS_45PercentRise_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_60PercentRise_MCP1")) TS_60PercentRise_MCP1 = rd->doubleUserRecords.get("TS_60PercentRise_MCP1"); else TS_60PercentRise_MCP1 = -999; - if (rd->doubleUserRecords.has("TS_60PercentRise_MCP2")) TS_60PercentRise_MCP2 = rd->doubleUserRecords.get("TS_60PercentRise_MCP2"); else TS_60PercentRise_MCP2 = -999; - if (rd->doubleUserRecords.has("amp_MCP1")) amp_MCP1 = rd->doubleUserRecords.get("amp_MCP1"); else amp_MCP1 = -999; - if (rd->doubleUserRecords.has("amp_MCP2")) amp_MCP2 = rd->doubleUserRecords.get("amp_MCP2"); else amp_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_MCP1")) TS_MCP1 = rd->doubleUserRecords.get("TS_MCP1"); else TS_MCP1 = -999; - if (rd->doubleUserRecords.has("TS_MCP2")) TS_MCP2 = rd->doubleUserRecords.get("TS_MCP2"); else TS_MCP2 = -999; - if (rd->doubleUserRecords.has("TS_MCP1_to_last_falling_Edge")) TS_MCP1_to_last_falling_Edge = rd->doubleUserRecords.get("TS_MCP1_to_last_falling_Edge"); else TS_MCP1_to_last_falling_Edge = -999; - if (rd->doubleUserRecords.has("TS_MCP2_to_last_falling_Edge")) TS_MCP2_to_last_falling_Edge = rd->doubleUserRecords.get("TS_MCP2_to_last_falling_Edge"); else TS_MCP2_to_last_falling_Edge = -999; + + valid_TS_MCP1 = (rd->intUserRecords.has("valid_TS_MCP1")) ? rd->intUserRecords.get("valid_TS_MCP1") : -999; + valid_TS_MCP2 = (rd->intUserRecords.has("valid_TS_MCP2")) ? rd->intUserRecords.get("valid_TS_MCP2") : -999; + noise_MCP1 = (rd->doubleUserRecords.has("noise_MCP1")) ? rd->doubleUserRecords.get("noise_MCP1") : -999; + noise_MCP2 = (rd->doubleUserRecords.has("noise_MCP2")) ? rd->doubleUserRecords.get("noise_MCP2") : -999; + TSpeak_MCP1 = (rd->doubleUserRecords.has("TSpeak_MCP1")) ? rd->doubleUserRecords.get("TSpeak_MCP1") : -999; + TSpeak_MCP2 = (rd->doubleUserRecords.has("TSpeak_MCP2")) ? rd->doubleUserRecords.get("TSpeak_MCP2") : -999; + amp_MCP1 = (rd->doubleUserRecords.has("amp_MCP1")) ? rd->doubleUserRecords.get("amp_MCP1") : -999; + amp_MCP2 = (rd->doubleUserRecords.has("amp_MCP2")) ? rd->doubleUserRecords.get("amp_MCP2") : -999; + ampFit_MCP1 = (rd->doubleUserRecords.has("ampFit_MCP1")) ? rd->doubleUserRecords.get("ampFit_MCP1") : -999; + ampFit_MCP2 = (rd->doubleUserRecords.has("ampFit_MCP2")) ? rd->doubleUserRecords.get("ampFit_MCP2") : -999; + TSfitPeak_MCP1 = (rd->doubleUserRecords.has("TSfitPeak_MCP1")) ? rd->doubleUserRecords.get("TSfitPeak_MCP1") : -999; + TSfitPeak_MCP2 = (rd->doubleUserRecords.has("TSfitPeak_MCP2")) ? rd->doubleUserRecords.get("TSfitPeak_MCP2") : -999; + TScf_MCP1 = (rd->doubleUserRecords.has("TScf_MCP1")) ? rd->doubleUserRecords.get("TScf_MCP1") : -999; + TScf_MCP2 = (rd->doubleUserRecords.has("TScf_MCP2")) ? rd->doubleUserRecords.get("TScf_MCP2") : -999; + charge5nsS_MCP1 = (rd->doubleUserRecords.has("charge5nsS_MCP1")) ? rd->doubleUserRecords.get("charge5nsS_MCP1") : -999; + charge5nsS_MCP2 = (rd->doubleUserRecords.has("charge5nsS_MCP2")) ? rd->doubleUserRecords.get("charge5nsS_MCP2") : -999; + charge5nsB_MCP1 = (rd->doubleUserRecords.has("charge5nsB_MCP1")) ? rd->doubleUserRecords.get("charge5nsB_MCP1") : -999; + charge5nsB_MCP2 = (rd->doubleUserRecords.has("charge5nsB_MCP2")) ? rd->doubleUserRecords.get("charge5nsB_MCP2") : -999; + TS_toClock_FE_MCP1 = (rd->doubleUserRecords.has("TS_toClock_FE_MCP1")) ? rd->doubleUserRecords.get("TS_toClock_FE_MCP1") : -999; + TS_toClock_FE_MCP2 = (rd->doubleUserRecords.has("TS_toClock_FE_MCP2")) ? rd->doubleUserRecords.get("TS_toClock_FE_MCP2") : -999; + meanClockFE = (rd->doubleUserRecords.has("meanClockFE")) ? rd->doubleUserRecords.get("meanClockFE") : -999; + rmsClockFE = (rd->doubleUserRecords.has("rmsClockFE")) ? rd->doubleUserRecords.get("rmsClockFE") : -999; tree_->Fill(); } diff --git a/Reco/src/MCPwaveform.cc b/Reco/src/MCPwaveform.cc new file mode 100644 index 0000000..c960f3f --- /dev/null +++ b/Reco/src/MCPwaveform.cc @@ -0,0 +1,390 @@ +#include "HGCal/Reco/interface/MCPwaveform.h" + + +MCPwaveform::MCPwaveform(){ + verbosity_ = 0; + nSamples_ = 0; + polarity_ = "neg"; + sigWinMin_ = 0; + sigWinMax_ = 1000; + peak_ = -1; + amp_ = -1; + base_ = -1; + noise_ = -1; + charge5nsS_ = -1; + charge5nsB_ = -1; + fQuality_ = 0; + Aparam_ = 0.; + Bparam_ = 0.; + fitAmp_ = -1; + fitTimeMax_ = -1; + fitTimeCF_ = -1; + nFitSamples_ = 0; + minForAmpFit_ = 0; + maxForAmpFit_ = 0; +} + + + +MCPwaveform::MCPwaveform(int verbosity, int nsamples, short* pulse, std::string polarity): + verbosity_(verbosity), + nSamples_(nsamples), polarity_(polarity), + sigWinMin_(0), sigWinMax_(1000), + peak_(-1), amp_(-1), base_(-1), noise_(-1), + charge5nsS_(-1), charge5nsB_(-1), + fQuality_(0), Aparam_(-1), Bparam_(-1), + fitAmp_(-1), fitTimeMax_(-1), fitTimeCF_(-1), + nFitSamples_(5), minForAmpFit_(5), maxForAmpFit_(5) +{ + for(int ij=0; ij 3. * charge5nsB_ && amp_ > 200.) fQuality_ = 1; + if(amp_ > 3. * noise_ && amp_ > 100.) fQuality_ = 1; + + if(verbosity_ == 1) + std::cout << " charge5nsS_ = " << charge5nsS_ << " charge5nsB_ = " << charge5nsB_ << std::endl; +} + + +//Get the max/min amplitude wrt polarity +//maybe add criteria to identify double genuine signal (multiple particles) +int MCPwaveform::findAbsolutePeak(std::string polarity){ + if(verbosity_ == 1) + std::cout << " in findAbsolutePeak polarity_ = " << polarity_ << std::endl; + + //loop over all samples looking for min value + //taken by Florian and Thorben version: exclude some first samples + int startSample = 5; + amp_ = wf_[startSample-1]; + + for(int i = startSample; i < nSamples_; ++i){ + if((polarity_ == "neg" && wf_[i] < amp_) || + (polarity_ == "pos" && wf_[i] > amp_)){ + amp_ = wf_[i]; + peak_ = i; + } + else if(polarity_ != "pos" && polarity_ != "neg"){ + if(wf_[i] > amp_){ + amp_ = wf_[i]; + peak_ = i; + } + } + } + + if(verbosity_ == 1) + std::cout << " amp_ = " << amp_ << " peak_ = " << peak_ << std::endl; + + return peak_; +} + + +// calculate mean and variance of all samples outside of peak region +// keep original call with (90, -1) => consider range [10, peak_-90] +float MCPwaveform::getBaseline(int nbinsExcludedLeftOfPeak, int nbinsExcludedRightOfPeak) { + if(verbosity_ == 1) std::cout << " in getBaseline fQuality_ = " << fQuality_ << std::endl; + + float sum = 0; + float cnt = 0; + float var = 0; + + if(peak_ == -1) findAbsolutePeak(polarity_); + + int startBin = sigWinMin_; + int lastBin = sigWinMax_; + + //treat "excluded" as good range interval + if(nbinsExcludedRightOfPeak != -1 && nbinsExcludedLeftOfPeak != -1){ + startBin = nbinsExcludedLeftOfPeak; + lastBin = nbinsExcludedRightOfPeak; + } + else if(nbinsExcludedRightOfPeak != -1 && peak_ + nbinsExcludedRightOfPeak < 1000){ + startBin = peak_ + nbinsExcludedRightOfPeak; + lastBin = 1000; + } + else if(nbinsExcludedLeftOfPeak != -1 && peak_ - nbinsExcludedLeftOfPeak > 0){ + startBin = 10; + lastBin = peak_ - nbinsExcludedLeftOfPeak; + } + else { + startBin = nbinsExcludedLeftOfPeak; + lastBin = nbinsExcludedRightOfPeak; + } + + for(int i = startBin; i < lastBin; ++i){ + sum += wf_[i]; + var += pow(wf_[i], 2); + cnt += 1; + } + + // break if no cnts + if (cnt == 0) { + base_ = -1; + noise_ = -1; + return -1; + } + + // calculate variance + base_ = sum / cnt; + var = (cnt * var - pow(sum, 2)) / pow(cnt, 2); + noise_ = sqrt(var); + + return base_; +} + + +void MCPwaveform::substractBaseline(int nSamples, short* pulse){ + float sum = 0.; + float cnt = 0.; + + for(int i = 5; i < 50; i++) { + sum += pulse[i]; + cnt += 1; + } + + float baseline = sum/cnt; + if(verbosity_ == 1) + std::cout << " >>> baseline = " << baseline << std::endl; + + for(int i = 0; i < nSamples; i++) { + pulse[i] = pulse[i] - baseline; + } + return; +} + + +// Get the waveform integral in the given range +float MCPwaveform::getIntegral(int min, int max){ + if(verbosity_ == 1) std::cout << " in getIntegral fQuality_ = " << fQuality_ << std::endl; + + float integral = 0.; + for(int iSample=min; iSample= nSamples_) break; + integral += wf_.at(iSample); + } + + charge5nsB_ = integral; + return integral; +} + + +//Get the signal integral around the the max +float MCPwaveform::getSignalIntegral(int riseWin, int fallWin){ + if(verbosity_ == 1) std::cout << " in getSignalIntegral fQuality_ = " << fQuality_ << std::endl; + + //compute position of the max + if(peak_ == -1) findAbsolutePeak(); + if(peak_-riseWin < 0) { + riseWin = 0; + fallWin = riseWin+25; + } + + //compute integral + float integral = 0.; + for(int iSample=peak_-riseWin; iSample= nSamples_) break; + integral += wf_.at(iSample); + } + + charge5nsS_ = integral; + return integral; +} + + +//Get the integral of Abs(WF) over the given range +float MCPwaveform::getModIntegral(int min, int max){ + float integral = 0.; + for(int iSample=min; iSample= nSamples_) break; + + if(wf_.at(iSample) < 0) integral -= wf_.at(iSample); + else integral += wf_.at(iSample); + } + + charge5nsB_ = integral; + return integral; +} + + +// Get CF time for a given fraction and in a given range +// note rise time is about 1ns = 5samples => nFitSamples to adapt +float MCPwaveform::getTimeCF(float frac, int nFitSamples, int min, int max){ + if(verbosity_ == 1) + std::cout << " in getTimeCF fQuality_ = " << fQuality_ << " fitAmp_ = " << fitAmp_ << " frac = " << frac << std::endl; + + if(!fQuality_) return -1; + + if(frac == 1) return fitTimeMax_; + if(Aparam_ != -1 && Bparam_ != -1) return (fitAmp_ * frac - Aparam_) / Bparam_; + + if(nFitSamples != 0) nFitSamples_ = nFitSamples; + if(min != 0 && max != 0){ + minForAmpFit_ = min; + maxForAmpFit_ = max; + } + + //good with min, max == 3 + if(fitAmp_ == -1) getInterpolatedAmpMax(minForAmpFit_, maxForAmpFit_); + + int cfSample = (sigWinMin_ == -1) ? 0 : sigWinMin_; + //find first sample above Amax*frac + for(int iSample=peak_; iSample>min; --iSample){ + if(wf_.at(iSample) < fitAmp_*frac){ + cfSample = iSample; + break; + } + } + + if(verbosity_ == 1) + std::cout << " in getTimeCF cfSample = " << cfSample << " peak_ = " << peak_ << std::endl; + + //rescale nFitSamples in case bigger than useful interval + //maybe restrict to 5 samples anyway + if(std::abs(cfSample - peak_) < (nFitSamples_-1)/2){ + nFitSamples_ = 2*std::abs(cfSample - peak_)+1; + } + + //interpolate with A+Bx = frac * amp + float chi2cf = linearInterpolation(cfSample - (nFitSamples_-1)/2, cfSample + (nFitSamples_-1)/2); + fitTimeCF_ = (fitAmp_ * frac - Aparam_) / Bparam_; + + if(verbosity_ == 1) + std::cout << " fit time = " << fitTimeCF_ << " peak_ = " << peak_ << std::endl; + + return fitTimeCF_; +} + + +//Linear interpolation with A+Bx +float MCPwaveform::linearInterpolation(const int& min, const int& max, const int& sampleToSkip){ + if(verbosity_ == 1) + std::cout << " in linearInterpolation " << std::endl; + + // assume y = A + B x + // B = (NSxy Sy - Sx Sy)/(NSxx - SxSx) + // A = Sy/N - B (Sx/N) = (Sxx Sy - Sx Sxy)/(NSxx - SxSx) + float xx = 0.; + float xy = 0.; + float Sx = 0.; + float Sy = 0.; + float Sxx = 0.; + float Sxy = 0.; + + int usedSamples=0; + for(int iSample=min; iSample<=max; ++iSample){ + + if(iSample<0 || iSample >= int(wf_.size())) continue; + + xx = iSample * iSample; + xy = iSample * wf_[iSample]; + Sx += iSample; + Sy += wf_[iSample]; + Sxx += xx; + Sxy += xy; + ++usedSamples; + } + + float Delta = usedSamples * Sxx - Sx * Sx; + Aparam_ = (Sxx * Sy - Sx * Sxy) / Delta; + Bparam_ = (usedSamples * Sxy - Sx * Sy) / Delta; + + //---compute chi2--- + float chi2 = 0.; + float sigma2 = pow(noise_, 2); + for(int iSample=min; iSample<=max; ++iSample){ + if(iSample<0 || iSample >= int(wf_.size())) continue; + + chi2 += pow(wf_[iSample] - Aparam_ - Bparam_ * iSample, 2)/sigma2; + } + + return chi2/(usedSamples-2); +} + + +//Linear interpolation with A+Bx+Cxx +float MCPwaveform::pol2Interpolation(const int& min, const int& max, const int& nFitSamples){ + if(verbosity_ == 1) + std::cout << " in pol2Interpolation fQuality_ = " << fQuality_ + << " min = " << min << "max = " << max << std::endl; + + TH1F h_max("h_max", "", nFitSamples, min, max+1); + TF1 f_max("f_max", "pol2", min, max+1); + + int bin = 1; + for(int iSample=min; iSample<=max; ++iSample){ + h_max.SetBinContent(bin, wf_[iSample]); + h_max.SetBinError(bin, noise_); + ++bin; + if(verbosity_ == 1) + std::cout << "tg->SetPoint(" << bin-1 << ", " << iSample << ", " << wf_[iSample] << ");" << std::endl; + } + + auto fit_result = h_max.Fit(&f_max, "QRS"); + fitTimeMax_ = -1. * f_max.GetParameter(1) / (2.*f_max.GetParameter(2)); + fitAmp_ = f_max.Eval(fitTimeMax_); + + float chi2 = f_max.GetChisquare()/(nFitSamples-3); + + return chi2; +} + + +//Get the interpolated max/min amplitude +float MCPwaveform::getInterpolatedAmpMax(int min, int max){ + if(verbosity_ == 1) std::cout << " in getInterpolatedAmpMax fQuality_ = " << fQuality_ << std::endl; + + if(!fQuality_) return -1; + + //get maximum sample + if(peak_ == -1) + findAbsolutePeak(); + + if(peak_ - min < sigWinMin_) min = minForAmpFit_; + if(peak_ + max > sigWinMax_) max = maxForAmpFit_; + + pol2Interpolation(peak_- min, peak_+ max+1, max+min+1); + + return fitAmp_; +}