Capacitor pedestal calibration

Linear regression sums are collected inside PadmeReco and are exported to a json file. This behaveour can be enabled with `GenADCCellCalib` keyword in PadmeReconstruction.cfg under `[RUNNINGMODE]` group. The produced sums need to be interpreted by utils/v1742-offline-cell-calib-linear-regression.py -- poly5 method is slowest but provides best results. On a second iteration, the calibration data can be provided to the reconstruction via `ADCCellCalibData` keyword again under `[RUNNINGMODE]` group Fixes #238
PADME-Experiment · Dec 17, 2021 · 123b533 · 123b533
1 parent 9c0e826
commit 123b533
Show file tree

Hide file tree

Showing 16 changed files with 578 additions and 4 deletions.
diff --git a/PadmeReco/GNUmakefile b/PadmeReco/GNUmakefile
@@ -25,6 +25,7 @@ LIBS += -LRecoBase/lib -lRecoBase$(STATIC)
 LIBS += -L$(PADMEROOT)/lib -lPadmeRoot
 #LIBS += -L/usr/lib64/mysql -lmysqlclient
 LIBS += -Llib -lPadmeReco$(STATIC)
+LIBS += -ljsoncpp
 
 
 #EXTRALIBS += -LEventDisplay/lib -lEventDisplay$(STATIC) -lEve -lGeom -lRGL -lFTGL -lGui -lGed -lTreePlayer -lEG 

diff --git a/PadmeReco/RecoBase/include/PadmeVReconstruction.hh b/PadmeReco/RecoBase/include/PadmeVReconstruction.hh
@@ -22,6 +22,8 @@
 #include "PadmeVRecoConfig.hh"
 #include "RecoVChannelID.hh"
 #include "PadmeVCalibration.hh"
+#include"ADCCellCalibData.hh"
+#include"ADCCellCalibHistograms.hh"
 
 #include "TH1.h"
 #include <vector>
@@ -142,8 +144,12 @@ protected:
   PadmeVTrigger *fTriggerProcessor;
   PadmeVGeometry *fGeometry;
 
+  static std::shared_ptr<ADCCellCalibData> fADCCellCalibData;
+  static std::shared_ptr<ADCCellCalibHistograms> fADCCellCalibHist;
+
   Bool_t fWriteHits;
   Bool_t fWriteClusters;
+  bool fMakeADCCellCalib;
 
   GlobalRecoConfigOptions* fGlobalRecoConfigOptions;
   /*

diff --git a/PadmeReco/RecoBase/src/PadmeVReconstruction.cc b/PadmeReco/RecoBase/src/PadmeVReconstruction.cc
@@ -12,6 +12,8 @@
 //#include "Debug.h" 
 
 //PadmeVReconstruction::PadmeVReconstruction(TFile* HistoFile, TString Name, TString ConfigFileName) : PadmeVNamedModule(Name), fHistoFile(HistoFile), fRecoEvent(0)
+std::shared_ptr<ADCCellCalibData> PadmeVReconstruction::fADCCellCalibData;
+std::shared_ptr<ADCCellCalibHistograms> PadmeVReconstruction::fADCCellCalibHist;
 PadmeVReconstruction::PadmeVReconstruction(TFile* HistoFile, TString Name, TString ConfigFileName) : PadmeVNamedModule(Name), fHistoFile(HistoFile)
 {
   /*
@@ -40,7 +42,8 @@ PadmeVReconstruction::PadmeVReconstruction(TFile* HistoFile, TString Name, TStri
   fWriteClusters = true;
   fWriteHits     = (Bool_t)fConfig->GetParOrDefault("Output", "Hits"          , 1 );
   fWriteClusters = (Bool_t)fConfig->GetParOrDefault("Output", "Clusters"      , 1 );
-
+  fMakeADCCellCalib = fConfig->GetParOrDefault("RECO", "ADCCellCalib", 1);
+
   //  InitChannelID(fConfig);
 
   //----------- Parse config file for common parameters ----------//
@@ -375,6 +378,15 @@ void PadmeVReconstruction::BuildHits(TRawEvent* rawEv){
   for(Int_t iBoard = 0; iBoard < nBoards; iBoard++) {
     ADC = rawEv->ADCBoard(iBoard);
     if(GetConfig()->BoardIsMine( ADC->GetBoardId())) {
+      if(fADCCellCalibHist.use_count()){
+        fADCCellCalibHist->FillHistograms(ADC, &fName[0], false);
+      }
+      if(fMakeADCCellCalib && fADCCellCalibData.use_count()){
+        fADCCellCalibData->CalibrateADC(ADC);
+      }
+      if(fADCCellCalibHist.use_count()){
+        fADCCellCalibHist->FillHistograms(ADC, &fName[0], true);
+      }
       //Loop over the channels and perform reco
       for(unsigned ich = 0; ich < ADC->GetNADCChannels();ich++) {
 	TADCChannel* chn = ADC->ADCChannel(ich);

diff --git a/PadmeReco/config.template/EVeto.cfg b/PadmeReco/config.template/EVeto.cfg
@@ -17,6 +17,7 @@ InputImpedance	50
 VoltageBin	0.000244
 
 [RECO]
+ADCCellCalib                     1
 SignalWindow       		150
 SignalPreSamples		30
 SignalPostSamples		120

diff --git a/PadmeReco/config.template/HEPVeto.cfg b/PadmeReco/config.template/HEPVeto.cfg
@@ -18,6 +18,7 @@ InputImpedance	50
 VoltageBin	0.000244
 
 [RECO]
+ADCCellCalib                     1
 SignalWindow       		150
 SignalPreSamples		30
 SignalPostSamples		150

diff --git a/PadmeReco/config.template/PVeto.cfg b/PadmeReco/config.template/PVeto.cfg
@@ -17,6 +17,7 @@ InputImpedance	50
 VoltageBin	0.000244
 
 [RECO]
+ADCCellCalib                     1
 SignalWindow                   151
 SignalPreSamples               30
 SignalPostSamples              120

diff --git a/PadmeReco/config.template/PadmeReconstruction.cfg b/PadmeReco/config.template/PadmeReconstruction.cfg
@@ -24,6 +24,9 @@ Monitor     0
 #------------------
 Cosmics     0
 Debug       0
+#GenADCCellCalib adc-raw-calib.json
+#ADCCellCalibData poly5.json
+#FillADCCellCalibHists 1
 
 [RECOCONFIG]
 ECal    config/ECal.cfg

diff --git a/PadmeReco/include/ADCCellCalibData.hh b/PadmeReco/include/ADCCellCalibData.hh
@@ -0,0 +1,18 @@
+#ifndef __PADME_RECO__ADC_CELL_CALIB_DATA_H__
+#define __PADME_RECO__ADC_CELL_CALIB_DATA_H__ 0
+
+#include"TADCBoard.hh"
+#include<map>
+#include<string>
+
+class ADCCellCalibData{
+  public:
+    ADCCellCalibData(const std::string &caldatjson);
+    void CalibrateADC(TADCBoard* raw_board);
+
+  private:
+    ADCCellCalibData();
+    std::map<unsigned int, double[32][1024][6]> fData;
+};
+
+#endif
diff --git a/PadmeReco/include/ADCCellCalibHistograms.hh b/PadmeReco/include/ADCCellCalibHistograms.hh
@@ -0,0 +1,33 @@
+#ifndef __PADME_RECO__ADC_CELL_CALIB_HISTOGRAMS_H__
+#define __PADME_RECO__ADC_CELL_CALIB_HISTOGRAMS_H__ 0
+
+#include"TADCBoard.hh"
+
+#include<TH2S.h>
+#include<TFile.h>
+
+#include<string>
+
+class ADCCellCalibHistograms{
+  public:
+    struct hists_t{
+      TH2S* nocal_cell;
+      TH2S* calib_cell;
+      TH2S* nocal_time;
+      TH2S* calib_time;
+    };
+    void FillHistograms(TADCBoard* raw_board, std::string det_name, bool is_calibrated);
+    void WriteHistograms(TFile*);
+    ADCCellCalibHistograms(const std::string& dirname, unsigned ylow=3300, unsigned yhigh=4100):
+      fYLow(ylow),
+      fYHigh(yhigh),
+      fDirName(dirname)
+  {}
+  private:
+    ADCCellCalibHistograms();
+    unsigned fYLow, fYHigh;
+    std::string fDirName;
+    std::map<int, hists_t[32]> fHists;
+    void CreateHisos(TADCBoard* raw_board, const std::string &det_name);
+};
+#endif
diff --git a/PadmeReco/include/ADCCellCalibSums.hh b/PadmeReco/include/ADCCellCalibSums.hh
@@ -0,0 +1,67 @@
+#ifndef __PADME_RECO__ADC_CELL_CALIB_SUMS_H__
+#define __PADME_RECO__ADC_CELL_CALIB_SUMS_H__ 0
+
+#include "TRawEvent.hh"
+
+#include<stdint.h>
+#include<string>
+#include<map>
+#include<sstream>
+
+class ADCCellCalibSums{
+  private:
+    typedef unsigned __int128 uint128_t;
+    struct capacitor_sum_t{
+      capacitor_sum_t():
+        n     (0),
+        t     (0),
+        t2    (0),
+        t3    (0),
+        t4    (0),
+        t5    (0),
+        t6    (0),
+        t7    (0),
+        t8    (0),
+        t9    (0),
+        t10   (0),
+        Ct    (0),
+        tCt   (0),
+        t2Ct  (0),
+        t3Ct  (0),
+        t4Ct  (0),
+        t5Ct  (0)
+      {}
+      // An overflow can occur when the number of
+      // processed events is bigger than 65535
+      uint16_t     n   ; // 1<<16 65535
+      uint32_t     t   ; // 1<<32/1024**1 4194304
+      uint64_t     t2  ;
+      uint64_t     t3  ;
+      uint64_t     t4  ; // 1<<64/1024**4 16777216
+      uint128_t    t5  ;
+      uint128_t    t6  ;
+      uint128_t    t7  ;
+      uint128_t    t8  ;
+      uint128_t    t9  ;
+      uint128_t    t10 ; // 1<<128/1024**10 268435456
+      uint64_t     Ct  ;
+      uint64_t     tCt ;
+      uint64_t     t2Ct;
+      uint64_t     t3Ct; // 1<<64/(4096*1024**3) 4194304
+      uint128_t    t4Ct;
+      uint128_t    t5Ct;
+    };
+  public:
+    ADCCellCalibSums(const std::string& sums_out_file_name):
+      fNEvtsProcessed(0) {
+        fSumsOutFileName = sums_out_file_name;
+      }
+    void Process(TRawEvent* ev);
+    void WriteFile();
+  private:
+    ADCCellCalibSums();
+    std::string fSumsOutFileName;
+    std::map<int, capacitor_sum_t[32][1024]> fAllBoardSums;
+    unsigned fNEvtsProcessed;
+};
+#endif
diff --git a/PadmeReco/include/PadmeReconstruction.hh b/PadmeReco/include/PadmeReconstruction.hh
@@ -1,6 +1,9 @@
 #ifndef PadmeReconstruction_H
 #define PadmeReconstruction_H
 
+
+#include"ADCCellCalibSums.hh"
+
 #include "TFile.h"
 #include "TChain.h"
 #include "TTree.h"
@@ -107,7 +110,8 @@ private:
   TECalRecoEvent*    fECalRecoEvent;
   TSACRecoEvent*     fSACRecoEvent;
   TTPixRecoEvent*    fTPixRecoEvent;
-
+
+  ADCCellCalibSums* fADCCellCalibSums;
   typedef std::vector<PadmeVReconstruction*> RecoVector;
   RecoVector fRecoLibrary;
 

diff --git a/PadmeReco/src/ADCCellCalibData.cc b/PadmeReco/src/ADCCellCalibData.cc
@@ -0,0 +1,48 @@
+#include"ADCCellCalibData.hh"
+#include<json/json.h>
+#include<fstream>
+
+ADCCellCalibData::ADCCellCalibData(const std::string &caldatjson){
+  if(caldatjson.size()==0){
+    return;
+  }
+  Json::Value caldat;
+  std::ifstream caldat_file(caldatjson.c_str(), std::ifstream::binary);
+  caldat_file >> caldat;
+  for(auto it=caldat.begin();it!=caldat.end();it++){
+    for(int ch=0;ch<32;++ch){
+      for(int cell=0;cell<1024;++cell){
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][0]=
+          caldat[it.key().asString()][ch][cell][0].asDouble();
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][1]=
+          caldat[it.key().asString()][ch][cell][1].asDouble();
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][2]=
+          caldat[it.key().asString()][ch][cell][2].asDouble();
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][3]=
+          caldat[it.key().asString()][ch][cell][3].asDouble();
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][4]=
+          caldat[it.key().asString()][ch][cell][4].asDouble();
+        fData[std::atoi(it.key().asString().c_str())][ch][cell][5]=
+          caldat[it.key().asString()][ch][cell][5].asDouble();
+      }
+    }
+  }
+}
+
+void ADCCellCalibData::CalibrateADC(TADCBoard* raw_board){
+  int n_chans = raw_board->GetNADCChannels();
+  for(int chan_i=0; chan_i<n_chans; chan_i++){
+    TADCChannel* raw_chan = raw_board->ADCChannel(chan_i);
+    int start_cell = raw_board->ADCTrigger(chan_i/8)->GetStartIndexCell();
+    double (*chan_caldat)[6] = fData[raw_board->GetBoardSN()][chan_i];
+    auto raw_samp = raw_chan->GetSamplesArray();
+    int n_samp = raw_chan->GetNSamples();
+    for(int t=0; t<n_samp; t++){
+      long long int tPower = 1;
+      for(int power=0; power<6; power++){
+        raw_samp[t] -= chan_caldat[(t+start_cell)%1024][power]*tPower;
+        tPower*=t;
+      }
+    }
+  }
+}
diff --git a/PadmeReco/src/ADCCellCalibHistograms.cc b/PadmeReco/src/ADCCellCalibHistograms.cc
@@ -0,0 +1,51 @@
+#include"ADCCellCalibHistograms.hh"
+
+#include<sstream>
+
+void ADCCellCalibHistograms::FillHistograms(TADCBoard* raw_board, std::string det_name, bool is_calibrated){
+  int board_sn = raw_board->GetBoardSN();
+  if(fHists.count(board_sn)==0){
+    CreateHisos(raw_board, det_name);
+  }
+  for(unsigned chan_i=0; chan_i<raw_board->GetNADCChannels(); chan_i++){
+    TADCChannel *raw_chan = raw_board->ADCChannel(chan_i);
+    hists_t& hist = fHists[board_sn][chan_i];
+    int start_index_cell = raw_board->ADCTrigger(chan_i/8)->GetStartIndexCell();
+    for(unsigned t=0; t < 1000; t++){ // ignore last 24 samples
+      const int p  = (t+start_index_cell)%raw_chan->GetNSamples();
+      const int Ct = raw_chan->GetSample(t);
+      if(is_calibrated){
+        hist.calib_cell->Fill(p, Ct);
+        hist.calib_time->Fill(t, Ct);
+      }else{
+        hist.nocal_cell->Fill(p, Ct);
+        hist.nocal_time->Fill(t, Ct);
+      }
+    }
+  }
+}
+
+void ADCCellCalibHistograms::CreateHisos(TADCBoard* raw_board, const std::string &det_name){
+  int board_sn = raw_board->GetBoardSN();
+  for(int chan_i=0; chan_i<32; chan_i++){
+    std::stringstream name;
+    name<<det_name<<"_sn"<<board_sn<<"_b"<<int(raw_board->GetBoardId())<<"_ch"<<chan_i<<"_";
+    fHists[board_sn][chan_i].nocal_cell = new TH2S((name.str()+"nocal_cell").c_str(),(name.str()+"nocal_cell;CellID; "       ).c_str(),1024,0-.5,1024-.5, fYHigh-fYLow, fYLow-.5,fYHigh-.5);
+    fHists[board_sn][chan_i].calib_cell = new TH2S((name.str()+"calib_cell").c_str(),(name.str()+"calib_cell;CellID; "       ).c_str(),1024,0-.5,1024-.5, fYHigh-fYLow, fYLow-.5,fYHigh-.5);
+    fHists[board_sn][chan_i].nocal_time = new TH2S((name.str()+"nocal_time").c_str(),(name.str()+"nocal_time;Time[samples]; ").c_str(),1024,0-.5,1024-.5, fYHigh-fYLow, fYLow-.5,fYHigh-.5);
+    fHists[board_sn][chan_i].calib_time = new TH2S((name.str()+"calib_time").c_str(),(name.str()+"calib_time;Time[samples]; ").c_str(),1024,0-.5,1024-.5, fYHigh-fYLow, fYLow-.5,fYHigh-.5);
+  }
+}
+
+void ADCCellCalibHistograms::WriteHistograms(TFile*rootfile){
+  rootfile->cd("/");
+  rootfile->mkdir(fDirName.c_str())->cd();
+  for(auto hist_it=fHists.begin(); hist_it!=fHists.end(); hist_it++){
+    for(int chan_i=0; chan_i<32; chan_i++){
+      hist_it->second[chan_i].nocal_cell->Write();
+      hist_it->second[chan_i].calib_cell->Write();
+      hist_it->second[chan_i].nocal_time->Write();
+      hist_it->second[chan_i].calib_time->Write();
+    }
+  }
+}