diff --git a/.gitmodules b/.gitmodules index 3e8cb207..bd057b97 100755 --- a/.gitmodules +++ b/.gitmodules @@ -1,6 +1,6 @@ [submodule "DynamicTTree"] path = DynamicTTree - url = https://github.com/simonepigazzini/DynamicTTree.git + url = https://github.com/pmeridian/DynamicTTree.git [submodule "CfgManager"] path = CfgManager - url = https://github.com/simonepigazzini/CfgManager.git + url = https://github.com/pmeridian/CfgManager.git diff --git a/DynamicTTree b/DynamicTTree index d08f6279..d9fce079 160000 --- a/DynamicTTree +++ b/DynamicTTree @@ -1 +1 @@ -Subproject commit d08f6279fcc89435d6750d6a51d5450635d616a8 +Subproject commit d9fce079bdbc887d3f15c36b1c6781c634bb8911 diff --git a/cfg/ECAL_H4_Oct2021/ECAL_H4_Phase1_base.cfg b/cfg/ECAL_H4_Oct2021/ECAL_H4_Phase1_base.cfg index 5bb677a1..c2c8e9f6 100644 --- a/cfg/ECAL_H4_Oct2021/ECAL_H4_Phase1_base.cfg +++ b/cfg/ECAL_H4_Oct2021/ECAL_H4_Phase1_base.cfg @@ -1,10 +1,12 @@ -outNameSuffix /eos/cms/store/group/dpg_ecal/comm_ecal/upgrade/testbeam/ECALTB_H4_Oct2021/Reco/ +outNameSuffix /tmp/meridian/ path2data /eos/cms/store/group/dpg_ecal/comm_ecal/upgrade/testbeam/ECALTB_H4_Oct2021/DataTree/ pluginList H4Hodo DigiReco WFReco run 10467 maxEvents 1000 maxFiles -1 +dataType H4Tree +preProcessorType H4DAQPreProcessor VFEs B1 B2 B3 B4 B5 \ diff --git a/cfg/MTDTB_FNAL_Mar2023/SCOPE_reco_base.cfg b/cfg/MTDTB_FNAL_Mar2023/SCOPE_reco_base.cfg new file mode 100644 index 00000000..5a075fff --- /dev/null +++ b/cfg/MTDTB_FNAL_Mar2023/SCOPE_reco_base.cfg @@ -0,0 +1,105 @@ + +outNameSuffix ntuples/ +path2data /eos/cms/store/group/dpg_mtd/comm_mtd/TB/MTDTB_FNAL_Mar2023/KeySightScope/TimingDAQRECO/RecoWithTracks/v1/ +run 66353 +maxEvents 100 +maxFiles -1 +dataType scopeFNALTree +preProcessorType ScopeFNALPreProcessor +pluginList DigiReco WFReco + + +SCOPE TOFCLK_P TOFCLK_P_1 TOFCLK_P_2 TOFCLK_M TOFCLK_M_1 TOFCLK_M_2 MCP + +#---VFE config + +pluginType DigitizerReco +channelsNames= SCOPE + + +#---WF reconstruction for VFEs + +pluginType WFAnalyzer +srcInstanceName DigiReco +channelsNames= DigiReco.channelsNames +timeRecoTypes LED CFD CLK +fillWFtree 1 +WFtreePrescale 10 + + +#---Scope channels + +digiBoard 0 +digiGroup 0 +digiChannel 0 +type Clock +polarity 1 +nSamples 800 +tUnit 0.05 +#CFD +CLK -0.8 0.8 +LED 0 2 2 100 400 +baselineFit 1 +baselineWin 150 600 +baselineInt 150 600 +signalWin 100 700 5 pol2 +signalInt 100 700 + +file data/tmpl_CLK_FNAL_Mar2023_large.root tmpl_TOFCLK_large +fitWin 0. 15 45 + + + + +copyChannel TOFCLK_P 0 400 +nSamples 400 +CLK -0.8 0.8 50 350 +LED 0 2 2 10 350 +baselineWin 50 350 +baselineInt 50 350 +signalWin 50 350 5 pol2 +signalInt 50 350 + + + +copyChannel TOFCLK_P 400 800 + + + +polarity -1 + +file data/tmpl_CLK_M_FNAL_Mar2023_large.root tmpl_TOFCLK_M_large +fitWin 0. 15 45 + + + + +copyChannel TOFCLK_M 0 400 + +file data/tmpl_CLK_M_FNAL_Mar2023_large.root tmpl_TOFCLK_M_large +fitWin 0. 15 45 + + + + +copyChannel TOFCLK_M 400 800 + + + +digiBoard 0 +digiGroup 0 +digiChannel 3 +polarity -1 +nSamples 800 +tUnit 0.05 +baselineWin 1 100 +baselineInt 10 90 +signalWin 50 700 5 gaus +signalInt 40 40 +CFD 0.5 5 +LED 20 2 2 + +file data/tmpl_MCP_FNAL_Mar2023.root tmpl_MCP +fitWin 0. 10 21 + + diff --git a/interface/DataLoader.h b/interface/DataLoader.h index 9a455afe..28c88fda 100644 --- a/interface/DataLoader.h +++ b/interface/DataLoader.h @@ -7,8 +7,9 @@ #include #include "CfgManager/interface/CfgManager.h" +#include "DynamicTTree/interface/DynamicTTreeBase.h" -#include "interface/H4Tree.h" +using namespace std; class DataLoader { @@ -21,7 +22,7 @@ class DataLoader ~DataLoader() {}; //----getters--- - inline H4Tree& GetTree() {return *inTree_;}; + inline DynamicTTreeBase* GetTree() {return inTree_;}; inline int GetNFiles() {return fileList_.size();}; inline int GetNFilesProcessed() {return iFile_;}; inline long int GetCurrentEntry() {return inTree_->GetCurrentEntry();}; @@ -39,9 +40,9 @@ class DataLoader vector fileList_; int iFile_; TFile* currentFile_; - H4Tree* inTree_; + DynamicTTreeBase* inTree_; bool firstEventInSpill_; - + string dataType_; }; #endif diff --git a/interface/H4Tree.h b/interface/H4Tree.h index d1a518fd..3f87e300 100644 --- a/interface/H4Tree.h +++ b/interface/H4Tree.h @@ -12,10 +12,19 @@ using namespace std; #define MAX_ADC_CHANNELS 500000 +#define MAX_DIGI_SAMPLES 100000 #define MAX_TDC_CHANNELS 200 +#define MAX_SCALER_WORDS 16 +#define MAX_PATTERNS 16 +#define MAX_PATTERNS_SHODO 16 +#define SMALL_HODO_X_NFIBERS 8 +#define SMALL_HODO_Y_NFIBERS 8 +#define MAX_TRIG_WORDS 32 +#define MAX_RO 100 + typedef unsigned long int uint32; -typedef unsigned long long int uint64; +typedef unsigned long long uint64; //**************************************************************************************** //----------Helper functions-------------------------------------------------------------- @@ -63,24 +72,25 @@ typedef std::unordered_map bgc_map_t; DATA(unsigned int, nTriggerWords) #define DATA_VECT_TABLE \ - DATA(unsigned int, evtTimeBoard, nEvtTimes) \ - DATA(uint64, evtTime, nEvtTimes) \ - DATA(unsigned int, adcBoard, MAX_ADC_CHANNELS) \ - DATA(unsigned int, adcChannel, MAX_ADC_CHANNELS) \ - DATA(unsigned int, adcData, MAX_ADC_CHANNELS) \ - DATA(unsigned int, tdcChannel, MAX_TDC_CHANNELS) \ - DATA(unsigned int, tdcData, MAX_TDC_CHANNELS) \ - DATA(unsigned int, pattern, nPatterns) \ - DATA(unsigned int, patternBoard, nPatterns) \ - DATA(unsigned int, patternChannel, nPatterns) \ - DATA(unsigned int, triggerWords, nTriggerWords) \ - DATA(unsigned int, triggerWordsBoard, nTriggerWords) \ - DATA(int, digiBoard, nDigiSamples) \ - DATA(unsigned int, digiGroup, nDigiSamples) \ - DATA(unsigned int, digiChannel, nDigiSamples) \ - DATA(unsigned int, digiStartIndexCell, nDigiSamples) \ - DATA(float, digiSampleValue, nDigiSamples) \ - DATA(float, digiSampleGain, nDigiSamples) + DATA(unsigned int, evtTimeBoard, nEvtTimes, MAX_RO) \ + DATA(unsigned long long, evtTime, nEvtTimes, MAX_RO) \ + DATA(unsigned int, adcBoard, nAdcChannels, MAX_ADC_CHANNELS) \ + DATA(unsigned int, adcChannel, nAdcChannels, MAX_ADC_CHANNELS) \ + DATA(unsigned int, adcData, nAdcChannels, MAX_ADC_CHANNELS) \ + DATA(unsigned int, tdcChannel, nTdcChannels, MAX_TDC_CHANNELS) \ + DATA(unsigned int, tdcData, nTdcChannels, MAX_TDC_CHANNELS) \ + DATA(unsigned int, pattern, nPatterns, MAX_PATTERNS) \ + DATA(unsigned int, patternBoard, nPatterns, MAX_PATTERNS) \ + DATA(unsigned int, patternChannel, nPatterns, MAX_PATTERNS) \ + DATA(unsigned int, triggerWords, nTriggerWords, MAX_TRIG_WORDS) \ + DATA(unsigned int, triggerWordsBoard, nTriggerWords, MAX_TRIG_WORDS) \ + DATA(int, digiBoard, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(unsigned int, digiGroup, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(unsigned int, digiChannel, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(unsigned int, digiStartIndexCell, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(float, digiSampleValue, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(float, digiSampleTime, nDigiSamples, MAX_DIGI_SAMPLES) \ + DATA(float, digiSampleGain, nDigiSamples, MAX_DIGI_SAMPLES) #include "DynamicTTree/interface/DynamicTTreeInterface.h" @@ -88,10 +98,16 @@ typedef std::unordered_map bgc_map_t; #undef DATA_TABLE #undef DATA_VECT_TABLE + class H4Tree : public H4TreeBase { public: //---ctors + H4Tree(const char* name="", const char* title=""): + H4TreeBase(name,title) + { + } + H4Tree(TChain* t): H4TreeBase(t) { diff --git a/interface/PluginLoader.h b/interface/PluginLoader.h index 8de0d962..a1c29041 100644 --- a/interface/PluginLoader.h +++ b/interface/PluginLoader.h @@ -49,7 +49,10 @@ template inline PluginLoader

::PluginLoader(string plugin_name) cout << ">>> PluginLoader: " << dlerror() << endl; exit(1); } + pluginCreator = (P* (*)(void))dlsym(pluginHandle, "create"); + std::cout << "Loaded " << plugin_name << ": " << pluginCreator << std::endl; + if((error = dlerror()) != NULL) { cout << ">>> PluginLoader: " << error << endl; diff --git a/interface/PreProcessorBase.h b/interface/PreProcessorBase.h new file mode 100644 index 00000000..cf088e64 --- /dev/null +++ b/interface/PreProcessorBase.h @@ -0,0 +1,71 @@ +#ifndef __PREPROCESSOR_BASE__ +#define __PREPROCESSOR_BASE__ + +#include +#include + +#include "CfgManager/interface/CfgManager.h" +#include "CfgManager/interface/CfgManagerT.h" +#include "interface/H4Tree.h" + +//**********Helper macros***************************************************************** +//---Helper macro to keep track of the preProcessor running: must be inserted at the +// beginning of each method +#define CHECKPOINT() \ + SetCurrentMethod(__FUNCTION__); + + +//**********PREPROCESSOR BASE CLASS************************************************************* +class PreProcessorBase +{ +public: + + //**********LOGGER LEVELS***************************************************************** + enum LoggerLevel { + INFO = 0, + WARN = 1, + ERR = 2 + }; + + //---ctors--- + PreProcessorBase(){}; + + //---dtor--- + virtual ~PreProcessorBase(){}; + + //---setters--- + void SetPluginType(const std::string& preProcessor) { preProcessorType_ = preProcessor; }; + void SetInstanceName(const std::string& instance) { instanceName_ = instance; }; + void SetCurrentMethod(std::string method) { currentMethod_ = method; }; + + //---getters--- + std::string GetPluginType() { return preProcessorType_; }; + std::string GetInstanceName() { return instanceName_; }; + std::string GetCurrentMethod() { return currentMethod_; }; + + //---utils--- + virtual bool Begin(CfgManager& opts) + { CHECKPOINT(); return true; }; + virtual H4Tree* ProcessEvent(DynamicTTreeBase* event, CfgManager& opts) + { CHECKPOINT(); return 0; }; + + void Log(std::string message, LoggerLevel lv=INFO); + +protected: + //---keep track of the preProcessor type as defined in the cfg + std::string preProcessorType_; + //---keep track of the preProcessor name as defined in the cfg + std::string instanceName_; + //---keep track of the current running preProcessor method. Must be called manually in each function; + std::string currentMethod_; +}; + +#undef CHECKPOINT + +//---To be inserted at the end of each preProcessor definition +#define DEFINE_PREPROCESSOR(NAME) \ + extern "C" PreProcessorBase* create() { return new NAME; } \ + extern "C" void destroy(PreProcessorBase* preProcessor) { delete preProcessor; } + +#endif + diff --git a/interface/WFClass.h b/interface/WFClass.h index ab4f4ec2..54752954 100644 --- a/interface/WFClass.h +++ b/interface/WFClass.h @@ -95,6 +95,9 @@ class WFClass : public TObject //---dtor--- ~WFClass() {}; + //crop method + void CropWF(WFClass& WF, int firstSample=-1, int lastSample=-1); + //---getters--- inline const vector* GetSamples() {return &samples_;}; inline const vector* GetTimes() {return ×_;}; @@ -151,15 +154,17 @@ class WFClass : public TObject void SetSignalIntegralWindow(int min, int max); void SetBaselineWindow(int min, int max); void SetBaselineIntegralWindow(int min, int max); - void SetTemplate(TH1* templateWF=NULL); + virtual void SetTemplate(TH1* templateWF=NULL); //---utils--- virtual void Reset(); bool ApplyCalibration(); virtual void AddSample(float sample); - virtual void AddSample(float sample, float gain) {AddSample(sample);}; + virtual void AddSampleTime(float sample, float time); + virtual void AddSampleGain(float sample, float gain) {AddSample(sample);}; WFBaseline SubtractBaseline(int min=-1, int max=-1); WFBaseline SubtractBaseline(float baseline); + virtual WFBaseline SubtractBaselineFit(int min=-1, int max=-1) {}; virtual WFFitResults TemplateFit(float amp_threshold=0., float offset=0., int lW=0, int hW=0); @@ -178,7 +183,7 @@ class WFClass : public TObject //---utils--- float BaselineRMS(); float LinearInterpolation(float& A, float& B, const int& min, const int& max, const int& skipSample=-1); - double TemplateChi2(const double* par=NULL); + virtual double TemplateChi2(const double* par=NULL); double AnalyticChi2(const double* par=NULL); protected: diff --git a/interface/WFClassClock.h b/interface/WFClassClock.h index d1c02c71..71a18aad 100644 --- a/interface/WFClassClock.h +++ b/interface/WFClassClock.h @@ -11,16 +11,22 @@ class WFClassClock : public WFClass public: //---ctors--- WFClassClock() {}; - WFClassClock(float tUnit); + WFClassClock(int polarity, float tUnit); + + WFBaseline SubtractBaselineFit(int min=-1, int max=-1) override; + void SetTemplate(TH1* templateWF=NULL) override; //---getters--- void AddSample(float sample) override; WFFitResults GetTime(string method, vector& params) override; WFFitResults GetTimeLE(float thr = 0, int nmFitSamples=2, int npFitSamples=2, int min=-1, int max=-1) override; - WFFitResults GetTimeCLK(float wleft=-1.3, float wright=1.3, int min=130, int max=900); + WFFitResults GetTimeCLK(float wleft=-1.3, float wright=1.3, int min=100, int max=700); WFFitResults TemplateFit(float ampl_threshold=0, float offset=0., int lW=0, int hW=0) override; float GetPeriod() override { return clkPeriod_; }; - float GetTemplateFitPeriod() override { return clkPeriod_; }; + float GetTemplateFitPeriod() override { return tmplFitPeriod_; }; + + protected: + double TemplateChi2(const double* par=NULL) override; private: float clkPeriod_; diff --git a/interface/WFClassLiTeDTU.h b/interface/WFClassLiTeDTU.h index 02cab056..f17f3199 100644 --- a/interface/WFClassLiTeDTU.h +++ b/interface/WFClassLiTeDTU.h @@ -25,7 +25,7 @@ class WFClassLiTeDTU : public WFClass //---utils--- void Reset() override; double TemplatesChi2(const double* par=NULL); - void AddSample(float sample, float gain) override; + void AddSampleGain(float sample, float gain) override; WFFitResults TemplateFit(float ampl_threshold=0, float offset=0., int lW=0, int hW=0) override; WFFitResultsScintPlusSpike TemplateFitScintPlusSpike(float amp_threshold=0., float offset=0., int lW=0, int hW=0) override; diff --git a/interface/scopeFNALTree.h b/interface/scopeFNALTree.h new file mode 100644 index 00000000..c3b5c6d4 --- /dev/null +++ b/interface/scopeFNALTree.h @@ -0,0 +1,58 @@ +#ifndef __SCOPE_FNAL_TREE__ +#define __SCOPE_FNAL_TREE__ + +#include +#include +#include +#include +#include + +#include "DynamicTTree/interface/DynamicTTreeBase.h" + +#define NCHANNELS 4 +#define NDIGIS 800 + +typedef unsigned long int uint32; +typedef unsigned long long int uint64; + +//**************************************************************************************** +//----------Tree reader class------------------------------------------------------------- +#define DYNAMIC_TREE_NAME scopeFNALTreeBase + +#define DATA_TABLE \ + DATA(unsigned int, i_evt) \ + +#define DATA_VECT_TABLE \ + DATA(float,channel, NCHANNELS*NDIGIS, NCHANNELS*NDIGIS ) \ + DATA(float,time, NDIGIS, NDIGIS) + +#include "DynamicTTree/interface/DynamicTTreeInterface.h" + +#undef DYNAMIC_TREE_NAME +#undef DATA_TABLE +//#undef DATA_VECT_TABLE + +class scopeFNALTree : public scopeFNALTreeBase +{ +public: + //---ctors + scopeFNALTree(TChain* t): + scopeFNALTreeBase(t) + { + Init(); + } + scopeFNALTree(TTree* t): + scopeFNALTreeBase(t) + { + Init(); + } + //---dtor + ~scopeFNALTree(); + + //---members + void Init(); + uint64 GetEntries(){ return tree_->GetEntriesFast(); }; + +}; + +#endif diff --git a/main/H4Reco.cpp b/main/H4Reco.cpp index 3dc60d2e..c9c8d4a9 100644 --- a/main/H4Reco.cpp +++ b/main/H4Reco.cpp @@ -19,6 +19,7 @@ #include "interface/PluginLoader.h" #include "interface/PluginBase.h" +#include "interface/PreProcessorBase.h" #include "interface/DataLoader.h" #include "interface/RecoTree.h" @@ -77,6 +78,30 @@ void TrackProcess(float* cpu, float* mem, float* vsz, float* rss) << "time lasted: " << time << endl; } +//----------Exception handler------------------------------------------------------------- +//---This fuction is meant to provide debug information by printing +//---the last called Plugin::Method while re-throwing the original exception +void HandleException(std::exception_ptr eptr, PreProcessorBase* plugin) +{ + try + { + if (eptr) + { + std::rethrow_exception(eptr); + } + } + catch(const std::exception& e) + { + auto plugin_type = std::regex_replace(typeid(plugin).name(), std::regex(".*[0-9]+"), ""); + std::cout << "\033[1;31m" << ">>>>> H4Reco ERROR! <<<<<" << "\033[0m" << std::endl + << "Error in: " << "\033[1;33m" << plugin_type << "::" << plugin->GetCurrentMethod() << "\033[0m" << std::endl + << "Plugin type: " << "\033[1;33m" << plugin->GetPluginType() << "\033[0m" << std::endl + << "Instance name: " << "\033[1;33m" << plugin->GetInstanceName() << "\033[0m" << std::endl + << "Caught exception: " << e.what() << std::endl; + exit(-1); + } +} + //----------Exception handler------------------------------------------------------------- //---This fuction is meant to provide debug information by printing //---the last called Plugin::Method while re-throwing the original exception @@ -144,6 +169,7 @@ int main(int argc, char* argv[]) string run = opts.GetOpt("h4reco.run"); int totLoops= opts.OptExist("h4reco.totLoops") ? opts.GetOpt("h4reco.totLoops") : 1; + //-----Load raw data----- vector spillOpt(1, to_string(spill)); opts.SetOpt("h4reco.firstSpill", spillOpt); @@ -164,6 +190,32 @@ int main(int argc, char* argv[]) RecoTree mainTree(&index); + //--- PreProcessor --- + string preProcessorType = opts.OptExist("h4reco.preProcessorType") ? opts.GetOpt("h4reco.preProcessorType") : "H4DAQPreProcessor" ; + PluginLoader* preProcessLoader = new PluginLoader(preProcessorType); + preProcessLoader->Create(); + PreProcessorBase* preProcessor = preProcessLoader->CreateInstance(preProcessorType); + if(!preProcessor) + { + cout << ">>> ERROR: preprocessor type " << preProcessorType << " is not defined." << endl; + return 0; + } + try + { + preProcessor->PreProcessorBase::Begin(opts); + if (!preProcessor->Begin(opts)) + { + cout << ">>> ERROR: preProcessor returned bad flag from Begin() call: " << preProcessor->GetInstanceName() << endl; + exit(-1); + } + } + catch(...) + { + eptr = std::current_exception(); + } + HandleException(eptr, preProcessor); + + //---Get plugin sequence--- PluginLoader* loader; vector* > pluginLoaders; @@ -215,7 +267,6 @@ int main(int argc, char* argv[]) } HandleException(eptr, plugin); } - int iLoop=0; int maxEvents = opts.OptExist("h4reco.maxEvents") ? opts.GetOpt("h4reco.maxEvents") : -1; @@ -270,18 +321,39 @@ int main(int argc, char* argv[]) } + + H4Tree* event=0; //---events loop while((dataLoader.NextEvent() && (nEvents < maxEvents || maxEvents == -1)) || (isSim && (nEvents < maxEvents))) { + if(dataLoader.FirstEventInSpill()) { cout << "\033[1;36m" << ">>> Processed spills: " << dataLoader.GetNFilesProcessed() << "/" << dataLoader.GetNFiles() << endl; cout << ">>> Processed events: " << nEvents << "\033[0m" << endl; TrackProcess(cpu, mem, vsz, rss); } - + + + try + { + preProcessor->PreProcessorBase::ProcessEvent(dataLoader.GetTree(),opts); + event=preProcessor->ProcessEvent(dataLoader.GetTree(),opts); + + if(!event) + { + cout << ">>> ERROR: preprocess error call: " << preProcessor->GetInstanceName() << endl; + exit(-1); + } + } + catch(...) + { + eptr = std::current_exception(); + } + HandleException(eptr, preProcessor); + //---set index value run*1e10+spill*1e4+event - index = dataLoader.GetTree().runNumber*1e9 + dataLoader.GetTree().spillNumber*1e5 + dataLoader.GetTree().evtNumber; + index = (*event).runNumber*1e9 + (*event).spillNumber*1e5 + (*event).evtNumber; //---call ProcessEvent for each plugin and check the return status bool status=true; @@ -290,10 +362,10 @@ int main(int argc, char* argv[]) try { //---fake call to base class for debug porpouses - plugin->PluginBase::ProcessEvent(dataLoader.GetTree(), pluginMap, opts); + plugin->PluginBase::ProcessEvent((*event), pluginMap, opts); //---real call + check for filters if(status) - status &= plugin->ProcessEvent(dataLoader.GetTree(), pluginMap, opts); + status &= plugin->ProcessEvent((*event), pluginMap, opts); } catch(...) { @@ -304,12 +376,12 @@ int main(int argc, char* argv[]) //---Fill the main tree with info variables and increase event counter mainTree.time_stamps.clear(); - for(int iT=0; iT& plugins, CfgManager& opts, u { nSamples_[channel] = opts.GetOpt(channel+".nSamples"); auto tUnit = opts.GetOpt(channel+".tUnit"); + auto polarity = opts.OptExist(channel+".polarity") ? opts.GetOpt(channel+".polarity") : 1; if(opts.OptExist(channel+".type")) { if(opts.GetOpt(channel+".type") == "NINO") - WFs_[channel] = new WFClassNINO(opts.GetOpt(channel+".polarity"), tUnit); + WFs_[channel] = new WFClassNINO(polarity, tUnit); else if(opts.GetOpt(channel+".type") == "Clock") - WFs_[channel] = new WFClassClock(tUnit); + WFs_[channel] = new WFClassClock(polarity,tUnit); else if(opts.GetOpt(channel+".type") == "LiTeDTU") - WFs_[channel] = new WFClassLiTeDTU(opts.GetOpt(channel+".polarity"), tUnit); + WFs_[channel] = new WFClassLiTeDTU(polarity, tUnit); } else - WFs_[channel] = new WFClass(opts.GetOpt(channel+".polarity"), tUnit); + WFs_[channel] = new WFClass(polarity, tUnit); //---set channel calibration if available @@ -90,34 +91,49 @@ bool DigitizerReco::ProcessEvent(H4Tree& event, map& plugin { //---reset and read new WFs_ WFs_[channel]->Reset(); + auto digiBd = opts.GetOpt(channel+".digiBoard"); auto digiGr = opts.GetOpt(channel+".digiGroup"); auto digiCh = opts.GetOpt(channel+".digiChannel"); auto offset = event.digiMap.at(make_tuple(digiBd, digiGr, digiCh)); auto max_sample = offset+std::min(nSamples_[channel], event.digiNSamplesMap[make_tuple(digiBd, digiGr, digiCh)]); auto iSample = offset; - - while(iSample < max_sample && event.digiBoard[iSample] != -1) - { - //Set the start index cell - if (iSample==offset) - WFs_[channel]->SetStartIndexCell(event.digiStartIndexCell[iSample]); - //---H4DAQ bug: sometimes ADC value is out of bound. - //---skip everything if one channel is bad - if(event.digiSampleValue[iSample] > 1e6) - { - evtStatus = false; - WFs_[channel]->AddSample(4095); - } - else - WFs_[channel]->AddSample(event.digiSampleValue[iSample], event.digiSampleGain[iSample]); - iSample++; - } - if(opts.OptExist(channel+".useTrigRef") && opts.GetOpt(channel+".useTrigRef")) - WFs_[channel]->SetTrigRef(trigRef); - if(WFs_[channel]->GetCalibration()) - WFs_[channel]->ApplyCalibration(); + if (!opts.OptExist(channel+".copyChannel")) + { + while( iSample < max_sample && event.digiBoard[iSample] != -1) + { + //Set the start index cell + if (iSample==offset) + WFs_[channel]->SetStartIndexCell(event.digiStartIndexCell[iSample]); + + //---H4DAQ bug: sometimes ADC value is out of bound. + //---skip everything if one channel is bad + if(event.digiSampleValue[iSample] > 1e6) + { + evtStatus = false; + WFs_[channel]->AddSample(4095); + } + else + { + auto dataType = opts.OptExist("h4reco.dataType") ? opts.GetOpt("h4reco.dataType") : "H4Tree"; + if (dataType=="scopeFNALTree") + WFs_[channel]->AddSampleTime(event.digiSampleValue[iSample], event.digiSampleTime[iSample]); + else + WFs_[channel]->AddSampleGain(event.digiSampleValue[iSample], event.digiSampleGain[iSample]); + } + iSample++; + } + if(opts.OptExist(channel+".useTrigRef") && opts.GetOpt(channel+".useTrigRef")) + WFs_[channel]->SetTrigRef(trigRef); + if(WFs_[channel]->GetCalibration()) + WFs_[channel]->ApplyCalibration(); + } + else + { + //Crop from another WF (need to be set first in the list) + WFs_[channel]->CropWF(*(WFs_[opts.GetOpt(channel+".copyChannel",0)]),opts.GetOpt(channel+".copyChannel",1),opts.GetOpt(channel+".copyChannel",2)); + } } if(!evtStatus) diff --git a/plugins/H4DAQPreProcessor.cc b/plugins/H4DAQPreProcessor.cc new file mode 100644 index 00000000..b68f4daf --- /dev/null +++ b/plugins/H4DAQPreProcessor.cc @@ -0,0 +1,12 @@ +#include "H4DAQPreProcessor.h" + +bool H4DAQPreProcessor::Begin(CfgManager& opts) +{ + return true; +} + +H4Tree* H4DAQPreProcessor::ProcessEvent(DynamicTTreeBase* event, CfgManager& opts) +{ + H4Tree* thisEvent=dynamic_cast(event); + return thisEvent; +} diff --git a/plugins/H4DAQPreProcessor.h b/plugins/H4DAQPreProcessor.h new file mode 100644 index 00000000..56752ca1 --- /dev/null +++ b/plugins/H4DAQPreProcessor.h @@ -0,0 +1,28 @@ +#ifndef __H4DAQ_PREPROCESSOR__ +#define __H4DAQ_PREPROCESSOR__ + +#include +#include + +#include "interface/PreProcessorBase.h" +#include "interface/H4Tree.h" + +class H4DAQPreProcessor: public PreProcessorBase +{ +public: + //---ctors--- + H4DAQPreProcessor(){}; + + //---dtor--- + ~H4DAQPreProcessor(){}; + + //---utils--- + bool Begin(CfgManager& opts); + H4Tree* ProcessEvent(DynamicTTreeBase* event, CfgManager& opts); + +private: +}; + +DEFINE_PREPROCESSOR(H4DAQPreProcessor); + +#endif diff --git a/plugins/ScopeFNALPreProcessor.cc b/plugins/ScopeFNALPreProcessor.cc new file mode 100644 index 00000000..15ae0910 --- /dev/null +++ b/plugins/ScopeFNALPreProcessor.cc @@ -0,0 +1,38 @@ +#include "ScopeFNALPreProcessor.h" + +bool ScopeFNALPreProcessor::Begin(CfgManager& opts) +{ + runNumber_=std::stoi(opts.GetOpt("h4reco.run")); + + for (int ich=0;ich(event); + h4Tree_.runNumber=runNumber_; + h4Tree_.spillNumber=1; + h4Tree_.evtNumber=thisEvent->i_evt; + + //filling digis + h4Tree_.nDigiSamples=NCHANNELS*NDIGIS; + for (int ich=0;ich +#include + +#include "interface/PreProcessorBase.h" +#include "interface/H4Tree.h" +#include "interface/scopeFNALTree.h" + +class ScopeFNALPreProcessor: public PreProcessorBase +{ +public: + //---ctors--- + ScopeFNALPreProcessor(): h4Tree_("H4Tree","H4TreeBase"){}; + + //---dtor--- + ~ScopeFNALPreProcessor(){}; + + //---utils--- + bool Begin(CfgManager& opts); + H4Tree* ProcessEvent(DynamicTTreeBase* event, CfgManager& opts); + +private: + //---internal data + int runNumber_; + H4Tree h4Tree_; +}; + +DEFINE_PREPROCESSOR(ScopeFNALPreProcessor); + +#endif diff --git a/plugins/WFAnalyzer.cc b/plugins/WFAnalyzer.cc index f12f2c2f..9c7f7553 100644 --- a/plugins/WFAnalyzer.cc +++ b/plugins/WFAnalyzer.cc @@ -6,7 +6,7 @@ bool WFAnalyzer::Begin(map& plugins, CfgManager& opts, uint trg_ = ""; //---inputs--- - for(auto& src : {"srcInstanceName", "trgInstanceName"}) + for(auto& src : {"srcInstanceName"}) { if(!opts.OptExist(instanceName_+"."+src)) { @@ -14,8 +14,12 @@ bool WFAnalyzer::Begin(map& plugins, CfgManager& opts, uint return false; } } + srcInstance_ = opts.GetOpt(instanceName_+".srcInstanceName"); - trgInstance_ = opts.GetOpt(instanceName_+".trgInstanceName"); + if(opts.OptExist(instanceName_+".trgInstanceName")) + trgInstance_ = opts.GetOpt(instanceName_+".trgInstanceName"); + else + trgInstance_=""; channelsNames_ = opts.GetOpt >(instanceName_+".channelsNames"); timeRecoTypes_ = opts.GetOpt >(instanceName_+".timeRecoTypes"); @@ -153,17 +157,21 @@ bool WFAnalyzer::ProcessEvent(H4Tree& event, map& plugins, fillWFtree = *digiTree_.index % opts.GetOpt(instanceName_+".WFtreePrescale") == 0; //---Check if trigger bit has changed from previous event - auto ctrg = ((TObjString*)(plugins[trgInstance_]->GetSharedData(trgInstance_+"_trg_bit", "", false))[0].obj)->GetString().Data(); + auto ctrg = "PHYS"; + + if (trgInstance_!="") + ctrg=((TObjString*)(plugins[trgInstance_]->GetSharedData(trgInstance_+"_trg_bit", "", false))[0].obj)->GetString().Data(); + if(ctrg != trg_ && templates_.find(ctrg) != templates_.end()) - { - for(auto& channel : channelsNames_) - { - if(templates_[ctrg].find(channel) != templates_[ctrg].end()) - WFs_[channel]->SetTemplate(templates_[ctrg][channel]); - } - trg_ = ctrg; - } + { + for(auto& channel : channelsNames_) + { + if(templates_[ctrg].find(channel) != templates_[ctrg].end()) + WFs_[channel]->SetTemplate(templates_[ctrg][channel]); + } + trg_ = ctrg; + } //---compute reco variables for(auto& channel : channelsNames_) @@ -198,10 +206,12 @@ bool WFAnalyzer::ProcessEvent(H4Tree& event, map& plugins, opts.GetOpt(channel+".signalInt", 1)); WFBaseline baselineInfo; if (opts.OptExist(channel+".baseline")) - baselineInfo = WFs_[channel]->SubtractBaseline(opts.GetOpt(channel+".baseline")); + baselineInfo = WFs_[channel]->SubtractBaseline(opts.GetOpt(channel+".baseline")); + else if (opts.OptExist(channel+".baselineFit")) + baselineInfo = WFs_[channel]->SubtractBaselineFit(); else - baselineInfo = WFs_[channel]->SubtractBaseline(); - + baselineInfo = WFs_[channel]->SubtractBaseline(); + //FIXME MAREMMA MAIALA WFFitResults interpolAmpMax; if(opts.OptExist(channel+".signalWin", 4)) @@ -259,8 +269,8 @@ bool WFAnalyzer::ProcessEvent(H4Tree& event, map& plugins, digiTree_.time_slope[outCh+iT*channelsNames_.size()] = -99; } } - digiTree_.period[outCh] = WFs_[channel]->GetPeriod(); + digiTree_.period[outCh] = WFs_[channel]->GetPeriod(); //---template fit (only specified channels) if(opts.OptExist(channel+".templateFit.file")) { diff --git a/scripts/setup.sh b/scripts/setup.sh index e86a7db6..a4781c10 100644 --- a/scripts/setup.sh +++ b/scripts/setup.sh @@ -1,7 +1,8 @@ #!/bin/sh -eval "$(conda shell.bash hook)" -conda activate /eos/project/c/cms-ecal-calibration/ecal-conda/h4analysis +#eval "$(conda shell.bash hook)" +#conda activate /eos/project/c/cms-ecal-calibration/ecal-conda/h4analysis +source /cvmfs/sft.cern.ch/lcg/releases/LCG_100/ROOT/v6.24.00/x86_64-centos7-gcc8-opt/ROOT-env.sh export LD_LIBRARY_PATH=./lib:DynamicTTree/lib/:CfgManager/lib/:$LD_LIBRARY_PATH export ROOT_INCLUDE_PATH=./interface:DynamicTTree/interface/:CfgManager/interface/:$ROOT_INCLUDE_PATH diff --git a/scripts/submitBatch.py b/scripts/submitBatch.py index 370c0b57..ce87ae17 100644 --- a/scripts/submitBatch.py +++ b/scripts/submitBatch.py @@ -260,4 +260,4 @@ def resubmitFailed(runs, outdir): firstspill += nfiles jobctr += 1 if not args.dryrun: - print('submitted {j} jobs to {b}'.format(j=jobctr, b=args.batch)) \ No newline at end of file + print('submitted {j} jobs to {b}'.format(j=jobctr, b=args.batch)) diff --git a/src/DataLoader.cc b/src/DataLoader.cc index 20c062c1..3f541a34 100644 --- a/src/DataLoader.cc +++ b/src/DataLoader.cc @@ -1,4 +1,6 @@ #include "interface/DataLoader.h" +#include "interface/H4Tree.h" +#include "interface/scopeFNALTree.h" //**********Contructor******************************************************************** DataLoader::DataLoader(CfgManager& opts): @@ -13,6 +15,7 @@ DataLoader::DataLoader(CfgManager& opts): //---Get list of files from run folder and store the list bool DataLoader::ReadInputFiles() { + dataType_= opts_.OptExist("h4reco.dataType") ? opts_.GetOpt("h4reco.dataType") : "H4Tree"; int firstSpill=opts_.GetOpt("h4reco.firstSpill"); string ls_command; string file; @@ -25,6 +28,9 @@ bool DataLoader::ReadInputFiles() // if ( getMachineDomain() != "cern.ch" ) // ls_command = string("gfal-ls root://eoscms/"+path+run+" | grep 'root' > /tmp/"+run+".list"); // else + if (dataType_=="scopeFNALTree") + ls_command = string("ls "+path+" | grep "+run+" > /tmp/"+run+".list"); + else ls_command = string("ls "+path+run+" | grep 'root' > /tmp/"+run+".list"); } else if(path.find("srm://") != string::npos) @@ -38,13 +44,26 @@ bool DataLoader::ReadInputFiles() while(waveList >> file && (opts_.GetOpt("h4reco.maxFiles")<0 || fileList_.size()("h4reco.maxFiles")) ) { //---skip files before specified spill - auto currentSpill = std::stoi(file.substr(0, file.size()-4)); + int currentSpill; + if (dataType_=="scopeFNALTree") + currentSpill = 1; + else + currentSpill=std::stoi(file.substr(0, file.size()-4)); + if(firstSpill == -1 || currentSpill >= firstSpill) { if(path.find("/eos/cms") != string::npos) { - std::cout << "+++ Adding file " << (path+run+"/"+file).c_str() << std::endl; - fileList_.push_back((path+run+"/"+file).c_str()); + if (dataType_=="scopeFNALTree") + { + std::cout << "+++ Adding file " << (path+"/"+file).c_str() << std::endl; + fileList_.push_back((path+"/"+file).c_str()); + } + else + { + std::cout << "+++ Adding file " << (path+run+"/"+file).c_str() << std::endl; + fileList_.push_back((path+run+"/"+file).c_str()); + } } else if(path.find("srm://") != string::npos) { @@ -107,10 +126,12 @@ bool DataLoader::LoadNextFile() currentFile_ = TFile::Open(fileList_[iFile_].c_str(), "READ"); if(currentFile_) { - inTree_ = new H4Tree((TTree*)currentFile_->Get("H4tree")); - ++iFile_; - - return true; + if (dataType_ == "H4Tree") + inTree_ = new H4Tree((TTree*)currentFile_->Get("H4tree")); + else if (dataType_ == "scopeFNALTree") + inTree_ = new scopeFNALTree((TTree*)currentFile_->Get("pulse")); + ++iFile_; + return true; } else return false; diff --git a/src/H4Tree.cc b/src/H4Tree.cc index b904ae6b..f1f2b6b9 100644 --- a/src/H4Tree.cc +++ b/src/H4Tree.cc @@ -15,7 +15,7 @@ void H4Tree::Init() { currentDigiChannel = digiChannel[iSample]; currentDigiGroup = digiGroup[iSample]; - currentDigiBoard = digiBoard[iSample]; + currentDigiBoard = digiBoard[iSample]; digiMap[make_tuple(currentDigiBoard, currentDigiGroup, currentDigiChannel)] = iSample; digiNSamplesMap[make_tuple(currentDigiBoard, currentDigiGroup, currentDigiChannel)] = 1; } diff --git a/src/PreProcessorBase.cc b/src/PreProcessorBase.cc new file mode 100644 index 00000000..0fb55d42 --- /dev/null +++ b/src/PreProcessorBase.cc @@ -0,0 +1,13 @@ +#include "interface/PreProcessorBase.h" + +void PreProcessorBase::Log(std::string message, LoggerLevel lv) +{ + //---Format message with preProcessor type and instance name + std::map colors = { {INFO, "\033[1;36m"}, + {WARN, "\033[1;33m"}, + {ERR, "\033[1;31m"} }; + std::cout << colors[lv] << "["+preProcessorType_+"::"+instanceName_+"]: " + << "\033[0m" << message << std::endl; + + return; +} diff --git a/src/WFClass.cc b/src/WFClass.cc index 3b61561d..bc1061f0 100644 --- a/src/WFClass.cc +++ b/src/WFClass.cc @@ -460,8 +460,6 @@ void WFClass::Reset() tmplFitTimeErr_=-1; tmplFitAmp_=-1; tmplFitAmpShift_=0; - interpolatorMin_=-1; - interpolatorMax_=-1; uncalibSamples_.clear(); calibSamples_.clear(); gain_.clear(); @@ -509,6 +507,14 @@ void WFClass::AddSample(float sample) }; +void WFClass::AddSampleTime(float sample,float time) +{ + uncalibSamples_.push_back(polarity_*sample); + times_.push_back( time ); + samples_ = uncalibSamples_; +}; + + //---------estimate the baseline in a given range and then subtract it from the signal---- WFBaseline WFClass::SubtractBaseline(int min, int max) @@ -520,6 +526,7 @@ WFBaseline WFClass::SubtractBaseline(int min, int max) } //---compute baseline float baseline_=0; + int nSamples=0; for(int iSample=bWinMin_; iSample= samples_.size()) break; baseline_ += samples_.at(iSample); + ++nSamples; } - baseline_ = baseline_/((float)(bWinMax_-bWinMin_)); + + baseline_ = baseline_/((float)nSamples); //---subtract baseline for(unsigned int iSample=0; iSample(origin.samples_.begin()+firstSample,origin.samples_.begin()+lastSample); + times_ = std::vector(origin.times_.begin()+firstSample,origin.times_.begin()+lastSample); + startIndexCell_=origin.startIndexCell_; + trigRef_=origin.trigRef_; +} + //**********operators********************************************************************* //----------assignment-------------------------------------------------------------------- WFClass& WFClass::operator=(const WFClass& origin) diff --git a/src/WFClassClock.cc b/src/WFClassClock.cc index 714b97f8..52d7105f 100644 --- a/src/WFClassClock.cc +++ b/src/WFClassClock.cc @@ -1,11 +1,86 @@ #include "interface/WFClassClock.h" +#include "TError.h" //**********Constructors****************************************************************** -WFClassClock::WFClassClock(float tUnit): - WFClass(1., tUnit), clkPeriod_(-1), tmplFitPeriod_(-1) +WFClassClock::WFClassClock(int polarity,float tUnit): + WFClass(polarity, tUnit), clkPeriod_(-1), tmplFitPeriod_(-1) {} +//---------estimate the baseline from a template fit of the clock signal +WFBaseline WFClassClock::SubtractBaselineFit(int min, int max) +{ + if(min!=-1 && max==-1) + { + bWinMin_=min; + bWinMax_=max; + } + //---compute baseline + float baseline_=0; + bRMS_=10; //custom value as starting point... + + //Remove annoying Minuit2 warnings + gErrorIgnoreLevel = kError; + + //---setup minimization + fWinMin_=bWinMin_; + fWinMax_=bWinMax_; + int maxSample_=bWinMin_; + int minSample_=bWinMin_; + for(int iSample=bWinMin_; iSample= samples_.size()) + break; + if(samples_.at(iSample) > samples_.at(maxSample_)) + maxSample_ = iSample; + if(samples_.at(iSample) < samples_.at(minSample_)) + minSample_ = iSample; + } + + float baseline_mean=(samples_.at(maxSample_)+samples_.at(minSample_))/2.; + std::default_random_engine generator; + auto t0 = (times_.at(bWinMin_)+times_.at(bWinMax_))/2.; + std::uniform_real_distribution rnd_t0(-2,2); + ROOT::Math::Functor chi2(this, &WFClassClock::TemplateChi2, 2); + ROOT::Math::Minimizer* minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad"); + minimizer->SetMaxFunctionCalls(100000); + minimizer->SetMaxIterations(1000); + minimizer->SetTolerance(1e-4); + minimizer->SetPrintLevel(-1); + minimizer->SetFunction(chi2); + minimizer->SetLimitedVariable(0, "deltaV", baseline_mean, 1e-2, baseline_mean*0.7,baseline_mean*1.3); + minimizer->SetLimitedVariable(1, "deltaT", t0, 1e-3, t0-10, t0+10); + //---fit + minimizer->Minimize(); + int tries=0; + while(minimizer->Status() != 0 && tries < 1000) + { + minimizer->SetVariableValue(1, (times_[bWinMin_]+times_[bWinMax_])/2.+rnd_t0(generator)); + minimizer->Minimize(); + ++tries; + } + baseline_ = minimizer->X()[0]; + tmplFitAmp_ = minimizer->X()[0]; + tmplFitTime_ = minimizer->X()[1]; + delete minimizer; + + //---subtract baseline + for(unsigned int iSample=0; iSample& params) return GetTimeCLK(params[0], params[1]); else if(params.size()==3) return GetTimeCLK(params[0], params[1], params[2]); + else if(params.size()==4) + return GetTimeCLK(params[0], params[1], params[2],params[3]); else cout << ">>>ERROR: wrong number of arguments passed for CLK time computation: " << params.size() << ". usage: CLK [min max] [+/- 1]" << endl; @@ -70,8 +147,8 @@ WFFitResults WFClassClock::GetTimeLE(float thr, int nmFitSamples, int npFitSampl chi2le_ = LinearInterpolation(A, B, leSample_-nmFitSamples, leSample_+npFitSamples); leTime_ = (leThr_ - A) / B; } - return WFFitResults{leThr_, leTime_, chi2le_, B}; + } //----------Get time of a clock pulse----------------------------------------------------- @@ -127,7 +204,6 @@ WFFitResults WFClassClock::GetTimeCLK(float wleft, float wright, int min, int ma } //auto phase_err = h_phase_err.Fit("gaus", "QRSO", "goff")->Parameter(2); auto phase_err = h_phase_err.GetRMS(); - // return WFFitResults{0, phase, std::sqrt(phase_err/(N-1)), 0}; return WFFitResults{0, phase, phase_err, -1, 0}; } @@ -136,9 +212,12 @@ WFFitResults WFClassClock::GetTimeCLK(float wleft, float wright, int min, int ma //---Template fit of each single clock cycle WFFitResults WFClassClock::TemplateFit(float ampl_threshold, float offset, int lW, int hW) { + //Remove annoying Minuit2 warnings + gErrorIgnoreLevel = kError; + if(tmplFitAmp_ == -1) { - BaselineRMS(); + // BaselineRMS(); GetAmpMax(); bRMS_ = 10.; int N=0; @@ -163,10 +242,10 @@ WFFitResults WFClassClock::TemplateFit(float ampl_threshold, float offset, int l minimizer->SetMaxFunctionCalls(100000); minimizer->SetMaxIterations(1000); minimizer->SetTolerance(1e-2); - minimizer->SetPrintLevel(0); + minimizer->SetPrintLevel(-1); minimizer->SetFunction(chi2); - minimizer->SetLimitedVariable(0, "amplitude", GetAmpMax(), 1e-2, GetAmpMax()*0.8, GetAmpMax()*1.2); - minimizer->SetLimitedVariable(1, "deltaT", t0, 1e-3, times_[fWinMin_], times_[fWinMax_]); + minimizer->SetLimitedVariable(0, "deltaV", 0, 1e-2, -0.1*GetAmpMax(), 0.1*GetAmpMax() ); + minimizer->SetLimitedVariable(1, "deltaT", t0, 1e-3, times_[fWinMin_]-1., times_[fWinMax_]+1); //---fit minimizer->Minimize(); int tries=0; @@ -206,6 +285,7 @@ WFFitResults WFClassClock::TemplateFit(float ampl_threshold, float offset, int l //tmplFitTimeErr_ = h_phase_err.Fit("gaus", "QRSO", "goff")->Parameter(2); } + gErrorIgnoreLevel = kWarning; return WFFitResults{tmplFitAmp_, tmplFitTime_, tmplFitTimeErr_, TemplateChi2()/(hW-lW), 0}; } @@ -218,3 +298,71 @@ void WFClassClock::AddSample(float sample) times_.push_back( (samples_.size()-1.)*tUnit_ ); samples_ = uncalibSamples_; }; + +//----------Set the fit template---------------------------------------------------------- +void WFClassClock::SetTemplate(TH1* templateWF) +{ + //---check input + if(!templateWF) + { + cout << ">>>ERROR: template passed as input does not exist" << endl; + return; + } + + //---reset template fit variables + if(interpolator_) + delete interpolator_; + + interpolator_ = new ROOT::Math::Interpolator(0, ROOT::Math::Interpolation::kCSPLINE); + tmplFitTime_ = 0; + tmplFitAmp_ = -1; + + //---fill interpolator data + vector x, y; + for(int iBin=1; iBin<=templateWF->GetNbinsX(); ++iBin) + { + x.push_back(templateWF->GetBinCenter(iBin)-tmplFitTime_); + y.push_back(templateWF->GetBinContent(iBin)); + } + interpolator_->SetData(x, y); + interpolatorMin_ = templateWF->GetBinCenter(1)-tmplFitTime_; + interpolatorMax_ = templateWF->GetBinCenter(templateWF->GetNbinsX())-tmplFitTime_; + return; +} + +double WFClassClock::TemplateChi2(const double* par) +{ + double chi2 = 0; + double delta2 = 0; +#ifdef PARALLEL +#pragma omp parallel for reduction(+:chi2) +#endif + for(int iSample=fWinMin_; iSample<=fWinMax_; ++iSample) + { + if(iSample < 0 || iSample >= int(samples_.size()) || + (par && (times_.at(iSample)-par[1] < interpolatorMin_ || times_.at(iSample)-par[1] > interpolatorMax_) ) ) + { + cout << ">>>WARNING: template fit out of samples rage (chi2 set to 9999)" << endl; + chi2 += 9999; + } + else + { + //---fit: par[0]*ref_shape(t-par[1]) par[0]=amplitude, par[1]=DeltaT + //---if not fitting return chi2 value of best fit + if(par) + { + //auto deriv = par[0]*interpolator_->Deriv(times_[iSample]-par[1]); + auto err2 = bRMS_*bRMS_;// + pow(tUnit_/sqrt(12)*deriv/2, 2); + delta2 = pow((samples_.at(iSample) - par[0]+interpolator_->Eval(times_[iSample]-par[1])), 2)/err2; + } + else + { + //auto deriv = tmplFitAmp_*interpolator_->Deriv(times_[iSample]-tmplFitTime_); + auto err2 = bRMS_*bRMS_;// + pow(tUnit_/sqrt(12)*deriv/2, 2); + delta2 = pow((samples_.at(iSample) - tmplFitAmp_+interpolator_->Eval(times_[iSample]-tmplFitTime_)), 2)/err2; + } + chi2 += delta2; + } + } + return chi2; +} diff --git a/src/WFClassLiTeDTU.cc b/src/WFClassLiTeDTU.cc index 6184fc46..c861dcb4 100644 --- a/src/WFClassLiTeDTU.cc +++ b/src/WFClassLiTeDTU.cc @@ -13,7 +13,7 @@ WFClassLiTeDTU::WFClassLiTeDTU(int polarity, float tUnit, DigiChannelCalibration //---sample is inserted at the end of uncalibSamples_ //---the times vector is filled with the uncalibrated sample time computed from the time unit //---a gain for each sample can be added. This is stored in a separate vector as well multiplied to the sample value. -void WFClassLiTeDTU::AddSample(float sample, float gain) +void WFClassLiTeDTU::AddSampleGain(float sample, float gain) { uncalibSamples_.push_back(polarity_*sample*gain); gain_.push_back(gain); diff --git a/src/scopeFNALTree.cc b/src/scopeFNALTree.cc new file mode 100644 index 00000000..02713ba6 --- /dev/null +++ b/src/scopeFNALTree.cc @@ -0,0 +1,11 @@ +#include "interface/scopeFNALTree.h" + +void scopeFNALTree::Init() +{ +} + +scopeFNALTree::~scopeFNALTree() +{ + delete[] channel; + delete[] time; +}