From 26ca39acb4a20d6f06d4d87ed634fc14b0f4a821 Mon Sep 17 00:00:00 2001
From: Shinsuke OTA <ota@cns.s.u-tokyo.ac.jp>
Date: Tue, 12 Jun 2018 22:12:11 +0900
Subject: [PATCH] [add] charge compensation with rise time in pulse finder

---
 src/cat/TPulseFinder.cc | 29 ++++++++++++++++++++++++++++-
 src/cat/TPulseFinder.h  |  7 +++++--
 2 files changed, 33 insertions(+), 3 deletions(-)

diff --git a/src/cat/TPulseFinder.cc b/src/cat/TPulseFinder.cc
index 98c80e8..eb15b64 100644
--- a/src/cat/TPulseFinder.cc
+++ b/src/cat/TPulseFinder.cc
@@ -2,7 +2,7 @@
 /**
  * @file   TPulseFinder.cc
  * @date   Created : Feb 02, 2014 16:02:27 JST
- *   Last Modified : 2018-06-06 22:11:58 JST (ota)
+ *   Last Modified : 2018-06-12 21:28:46 JST (ota)
 
  * @author Shinsuke OTA <ota@cns.s.u-tokyo.ac.jp>
  *
@@ -16,12 +16,14 @@
 #include <numeric>
 #include <functional>
 #include <TArtSystemOfUnit.h>
+#include <TF1.h>
 
 ClassImp(art::TPulseFinder);
 
 using art::TPulseFinder;
 
 TPulseFinder::TPulseFinder()
+   : fFunChargeRiseCompensation(NULL)
 {
    StringVec_t def;
    def.push_back("at");
@@ -48,6 +50,8 @@ TPulseFinder::TPulseFinder()
    RegisterProcessorParameter("Time2Drift",
                               "type of time to calculate drift length (0: offset, 1: cfd (0.5), 2: peak",
                               fTime2Drift,0);
+   RegisterProcessorParameter("ChargeRiseCompensation","3 parameters for charge compensation with rise time",
+                              fChargeRiseCompensation,DoubleVec_t(0));
    FloatVec_t range;
    range.push_back(-500.);
    range.push_back(500.);
@@ -63,6 +67,14 @@ TPulseFinder::~TPulseFinder()
 void TPulseFinder::Init(TEventCollection *col)
 {
    TEventObject *obj = NULL;
+   if (fChargeRiseCompensation.size() > 0 && fChargeRiseCompensation.size() != 3) {
+      SetStateError(TString::Format("3 parameters required for charge compensation with rise time : %d parameters are given",
+                                    fChargeRiseCompensation.size()));
+      return;
+   } 
+      
+                    
+   
    if (!(obj = col->Get(fInputColName.front()))) {
       Info("Init","No such input %s",fInputColName.front().Data());
       return;
@@ -72,6 +84,14 @@ void TPulseFinder::Init(TEventCollection *col)
    col->Add(fOutputColName,fOutput,fOutputIsTransparent);
 
    fDriftVelocity *= TArtSystemOfUnit::cm / TArtSystemOfUnit::us;
+
+   if (fChargeRiseCompensation.size() == 3) {
+      fFunChargeRiseCompensation = new TF1("funChargeRiseCompensation",
+                                           "[2]-1/(1+TMath::Power(x/[0],[1]))",
+                                           0.,100.);
+      fFunChargeRiseCompensation->SetParameters(&fChargeRiseCompensation[0]);
+   }
+
 }
 
 
@@ -187,6 +207,13 @@ void TPulseFinder::Process()
                                   
 
             out->SetCharge(std::accumulate(itfirst,itlast,0.));
+
+            if (fFunChargeRiseCompensation) {
+               const Double_t c = out->GetCharge();
+               const Double_t dt = out->GetMaxSampleOffset() - out->GetLeadingEdgeOffset();
+               out->SetCharge(c / fFunChargeRiseCompensation->Eval(dt));
+            }
+            
             out->ResetQualityBit(TCatPulseShape::kInvalid);
 
 
diff --git a/src/cat/TPulseFinder.h b/src/cat/TPulseFinder.h
index 0c99ce5..725c1bd 100644
--- a/src/cat/TPulseFinder.h
+++ b/src/cat/TPulseFinder.h
@@ -2,7 +2,7 @@
 /**
  * @file   TPulseFinder.h
  * @date   Created : Feb 02, 2014 13:02:42 JST
- *   Last Modified : 2018-06-06 21:35:39 JST (ota)
+ *   Last Modified : 2018-06-12 21:21:46 JST (ota)
  * @author Shinsuke OTA <ota@cns.s.u-tokyo.ac.jp>
  *  
  *  
@@ -16,7 +16,7 @@
 namespace art {
    class TPulseFinder;
 }
-
+class TF1;
 class TClonesArray;
 
 class art::TPulseFinder  : public TProcessor {
@@ -46,6 +46,9 @@ protected:
    Float_t fDriftVelocity;
    Double_t fClockPeriod;
    Int_t fTime2Drift;
+
+   TF1 *fFunChargeRiseCompensation; //!
+   DoubleVec_t fChargeRiseCompensation;
    
    ClassDef(TPulseFinder,1);
 };
-- 
GitLab