add subdirectory for time calibration (TCAL)

R3BRootGroup · Aug 27, 2015 · 4c00a38 · 4c00a38
1 parent 42f919c
commit 4c00a38
Show file tree

Hide file tree

Showing 11 changed files with 719 additions and 0 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -280,6 +280,7 @@ add_subdirectory (mfi)
 add_subdirectory (psp)
 add_subdirectory (timestitcher)
 add_subdirectory (compilehelper)
+add_subdirectory (tcal)
 ENDIF (NOT MODULE)
 
 Execute_Process(COMMAND ${ROOT_CONFIG_EXECUTABLE} --has-opengl

diff --git a/tcal/CMakeLists.txt b/tcal/CMakeLists.txt
@@ -0,0 +1,42 @@
+# Create a library called "libR3BTCal" which includes the source files given in
+# the array .
+# The extension is already found.  Any number of sources could be listed here.
+
+Set(SYSTEM_INCLUDE_DIRECTORIES
+${SYSTEM_INCLUDE_DIRECTORIES}
+${BASE_INCLUDE_DIRECTORIES}
+)
+
+set(INCLUDE_DIRECTORIES
+#put here all directories where header files are located
+${R3BROOT_SOURCE_DIR}/tcal
+)
+
+include_directories( ${INCLUDE_DIRECTORIES})
+include_directories(SYSTEM ${SYSTEM_INCLUDE_DIRECTORIES})
+
+set(LINK_DIRECTORIES
+${ROOT_LIBRARY_DIR}
+${FAIRROOT_LIBRARY_DIR}
+)
+
+link_directories( ${LINK_DIRECTORIES})
+
+set(SRCS
+R3BTCalModulePar.cxx
+R3BTCalPar.cxx
+R3BTCalContFact.cxx
+R3BTCalEngine.cxx
+)
+
+# fill list of header files from list of source files
+# by exchanging the file extension
+CHANGE_FILE_EXTENSION(*.cxx *.h HEADERS "${SRCS}")
+
+Set(LINKDEF TCalLinkDef.h)
+Set(LIBRARY_NAME R3BTCal)
+Set(DEPENDENCIES
+    Base ParBase)
+
+GENERATE_LIBRARY()
+
diff --git a/tcal/R3BTCalContFact.cxx b/tcal/R3BTCalContFact.cxx
@@ -0,0 +1,70 @@
+/////////////////////////////////////////////////////////////
+//
+//  R3BTCalContFact
+//
+//  Factory for the parameter containers
+//
+/////////////////////////////////////////////////////////////
+#include "R3BTCalContFact.h"
+
+#include "R3BTCalPar.h"    // for R3BLandGeometryPar
+#include "FairParSet.h"    // for FairParSet
+#include "FairRuntimeDb.h" // for FairRuntimeDb
+#include "FairLogger.h"
+
+#include "Riosfwd.h" // for ostream
+#include "TList.h"   // for TList
+#include "TString.h" // for TString
+
+#include <string.h> // for strcmp, NULL
+
+using namespace std;
+
+ClassImp(R3BTCalContFact);
+
+static R3BTCalContFact gR3BTCalContFact;
+
+R3BTCalContFact::R3BTCalContFact()
+{
+    // Constructor (called when the library is loaded)
+    fName = "R3BTCalContFact";
+    fTitle = "Tutorial factory for parameter containers";
+    setAllContainers();
+    FairRuntimeDb::instance()->addContFactory(this);
+}
+
+void R3BTCalContFact::setAllContainers()
+{
+    /** Creates the Container objects with all accepted contexts and adds them to
+     *  the list of containers.*/
+
+    FairContainer* p1 = new FairContainer("LandTCalPar", "NeuLAND TCAL Calibration Parameters", "TestDefaultContext");
+    p1->addContext("TestNonDefaultContext");
+    containers->Add(p1);
+
+    FairContainer* p2 = new FairContainer("LosTCalPar", "LOS TCAL Calibration Parameters", "TestDefaultContext");
+    p2->addContext("TestNonDefaultContext");
+    containers->Add(p2);
+}
+
+FairParSet* R3BTCalContFact::createContainer(FairContainer* c)
+{
+    /** Calls the constructor of the corresponding parameter container.
+     * For an actual context, which is not an empty string and not the default context
+     * of this container, the name is concatinated with the context. */
+
+    const char* name = c->GetName();
+    LOG(INFO) << "R3BTCalContFact::createContainer : " << name << FairLogger::endl;
+    FairParSet* p = NULL;
+
+    if (strcmp(name, "LandTCalPar") == 0)
+    {
+        p = new R3BTCalPar(c->getConcatName().Data(), c->GetTitle(), c->getContext());
+    }
+    else if (strcmp(name, "LosTCalPar") == 0)
+    {
+        p = new R3BTCalPar(c->getConcatName().Data(), c->GetTitle(), c->getContext());
+    }
+
+    return p;
+}
diff --git a/tcal/R3BTCalContFact.h b/tcal/R3BTCalContFact.h
@@ -0,0 +1,25 @@
+#ifndef R3BTCALCONTFACT_H
+#define R3BTCALCONTFACT_H
+
+#include "FairContFact.h"
+
+#include "Rtypes.h"
+
+class FairParSet;
+
+class R3BTCalContFact : public FairContFact
+{
+  public:
+    R3BTCalContFact();
+    ~R3BTCalContFact()
+    {
+    }
+    FairParSet* createContainer(FairContainer* c);
+
+  private:
+    void setAllContainers();
+
+    ClassDef(R3BTCalContFact, 1);
+};
+
+#endif /* !R3BTCALCONTFACT_H */
diff --git a/tcal/R3BTCalEngine.cxx b/tcal/R3BTCalEngine.cxx
@@ -0,0 +1,195 @@
+
+#include <string>
+
+#include "TH1F.h"
+#include "TF1.h"
+
+#include "FairLogger.h"
+
+#include "R3BTCalPar.h"
+#include "R3BTCalEngine.h"
+
+R3BTCalEngine::R3BTCalEngine(R3BTCalPar* param, Int_t nModules, Int_t minStats)
+    : fMinStats(minStats)
+    , fNModules(nModules)
+    , fhData(new TH1F*[nModules])
+    , fhData100(new TH1F*[nModules])
+    , fhTime(new TH1F*[nModules])
+    , f1(new TF1("f1", "[0]", 1500., 2500.))
+    , fCal_Par(param)
+    , iMin(0)
+    , iMax(0)
+    , fClockFreq(0.)
+{
+    char strName[255];
+    for (Int_t i = 0; i < fNModules; i++)
+    {
+        sprintf(strName, "%s_tcaldata_%d", fCal_Par->GetName(), i);
+        fhData[i] = new TH1F(strName, "", 4096, 0.5, 4096.5);
+        sprintf(strName, "%s_time_%d", fCal_Par->GetName(), i);
+        fhTime[i] = new TH1F(strName, "", 4096, 0.5, 4096.5);
+    }
+}
+
+R3BTCalEngine::~R3BTCalEngine()
+{
+    for (Int_t i = 0; i < fNModules; i++)
+    {
+        if(fhData[i])
+        {
+            delete fhData[i];
+        }
+        if(fhTime[i])
+        {
+            delete fhTime[i];
+        }
+    }
+    delete fhData;
+    delete fhTime;
+    delete fhData100;
+    delete f1;
+}
+
+void R3BTCalEngine::Fill(Int_t iModule, Int_t tdc)
+{
+    if (iModule < fNModules && iModule >= 0)
+    {
+        fhData[iModule]->Fill(tdc);
+    }
+}
+
+void R3BTCalEngine::CalculateParamTacquila()
+{
+    fClockFreq = 1. / TACQUILA_CLOCK_MHZ * 1000.;
+
+    for (Int_t iModule = 0; iModule < fNModules; iModule++)
+    {
+        if(fhData[iModule]->GetEntries() < fMinStats)
+        {
+            continue;
+        }
+
+        // Define range of channels
+        fhData100[iModule] = (TH1F*)fhData[iModule]->Clone();
+        fhData100[iModule]->Rebin(8);
+        fhData100[iModule]->Fit(f1, "QNR");
+        for (Int_t i = 256; i >= 0; i--)
+        {
+            if (fhData100[iModule]->GetBinContent(i + 1) < 0.1 * f1->GetParameter(0))
+            {
+                iMin = i - 1;
+                break;
+            }
+        }
+        for (Int_t j = 256; j < 512; j++)
+        {
+            if (fhData100[iModule]->GetBinContent(j + 1) < 0.1 * f1->GetParameter(0))
+            {
+                iMax = j + 1;
+                break;
+            }
+        }
+        if (iMax <= iMin)
+        {
+            LOG(ERROR) << "Error in definition of TCAL range" << FairLogger::endl;
+            return;
+        }
+        iMin = Int_t(((Double_t)iMin - 0.5) * 8);
+        iMax = Int_t(((Double_t)iMax + 0.5) * 8);
+        if (iMin < 0 || iMax > 4095)
+        {
+            return;
+        }
+        LOG(INFO) << "R3BTCalEngine::CalculateParamTacquila() : Range of channels: " << iMin << " - " << iMax
+                  << FairLogger::endl;
+
+        Int_t nch = 0;
+        Int_t ibin = iMin;
+        Int_t group;
+        Double_t prev_time = 0.;
+
+        R3BTCalModulePar* pTCal = NULL;
+        pTCal = new R3BTCalModulePar();
+        pTCal->SetModuleId(iModule);
+
+        Int_t incr = 0;
+        while (ibin <= iMax)
+        {
+            // Iteratively compute parameter
+            group = CalculateBinTacquila(iModule, prev_time, ibin, 1);
+            // Fill time calibration parameter
+            for (Int_t i1 = ibin; i1 < (ibin + group); i1++)
+            {
+                fhTime[iModule]->SetBinContent(i1 + 1, prev_time);
+            }
+
+            LOG(DEBUG) << " Module: " << iModule << " bin range: " << ibin << " : " << ibin + group
+                       << " dbin: " << group << " time set: " << prev_time << FairLogger::endl;
+
+            if (pTCal)
+            {
+                pTCal->SetBinLowAt(ibin, incr);
+                pTCal->SetBinUpAt(ibin + group - 1, incr);
+                pTCal->SetTimeAt(prev_time, incr);
+                pTCal->IncrementNofChannels();
+            }
+
+            // Next range of channels
+            ibin += group;
+            nch += 1;
+            incr++;
+        }
+
+        fCal_Par->AddModulePar(pTCal);
+
+        LOG(INFO) << "R3BTCalEngine::CalculateParamTacquila() : Number of parameters: " << nch << FairLogger::endl;
+
+        fhData[iModule]->Write();
+        fhTime[iModule]->Write();
+
+        LOG(INFO) << "R3BTCalEngine::CalculateParamTacquila() : Module: " << iModule << " is calibrated." << FairLogger::endl << FairLogger::endl;
+
+    }
+
+    fCal_Par->setChanged();
+}
+
+Int_t R3BTCalEngine::CalculateBinTacquila(Int_t iModule, Double_t& prev_time, Int_t ibin, Int_t ngroup)
+{
+    if ((ibin + ngroup) > iMax)
+    {
+        Double_t total = fhData[iModule]->Integral(1, 4096);
+        Double_t itot = fhData[iModule]->Integral(1, (ibin + 1) + ngroup);
+        if (itot > 0. && total > 0.)
+        {
+            Double_t time = fClockFreq * itot / total; // time of channel in [ns]
+            LOG(DEBUG) << "R3BTCalEngine::CalculateBin() : bin=" << ibin << "  time=" << time << "  ngroup=" << ngroup
+                       << FairLogger::endl;
+            prev_time = time;
+        }
+        return ngroup;
+    }
+    Double_t total = fhData[iModule]->Integral(1, 4096);
+    Double_t itot = fhData[iModule]->Integral(1, (ibin + 1) + ngroup);
+    if (itot > 0. && total > 0.)
+    {
+        Double_t time = fClockFreq * itot / total; // time of channel in [ns]
+        Double_t diff = time - prev_time;        // time difference to previous range
+        if (diff * 1e3 < 70.)                    // check if below resolution
+        {
+            // Next iteration
+            return CalculateBinTacquila(iModule, prev_time, ibin, ngroup + 1);
+        }
+        else
+        {
+            // Finalize
+            LOG(DEBUG) << "R3BTCalEngine::CalculateBin() : bin=" << ibin << "  time=" << time << "  ngroup=" << ngroup
+                       << FairLogger::endl;
+            prev_time = time;
+            return ngroup;
+        }
+    }
+    return ngroup;
+}
+
+ClassImp(R3BTCalEngine)
diff --git a/tcal/R3BTCalEngine.h b/tcal/R3BTCalEngine.h
@@ -0,0 +1,43 @@
+
+#ifndef _R3BTCAL_ENGINE_
+#define _R3BTCAL_ENGINE_
+
+#define TACQUILA_CLOCK_MHZ 40.002903
+#define VFTX_CLOCK_MHZ 200
+
+#include "TObject.h"
+
+class TH1F;
+class TF1;
+class R3BTCalPar;
+
+class R3BTCalEngine : public TObject
+{
+  public:
+    R3BTCalEngine(R3BTCalPar* param, Int_t nModules, Int_t minStats = 10000);
+    virtual ~R3BTCalEngine();
+
+    void Fill(Int_t iModule, Int_t tdc);
+
+    void CalculateParamTacquila();
+
+  protected:
+    Int_t CalculateBinTacquila(Int_t iModule, Double_t& prev_time, Int_t ibin, Int_t ngroup);
+
+  private:
+    Int_t fMinStats;
+    Int_t fNModules;
+    TH1F** fhData;
+    TH1F** fhData100;
+    TH1F** fhTime;
+    TF1* f1;
+    R3BTCalPar* fCal_Par;
+    Int_t iMin;
+    Int_t iMax;
+    Double_t fClockFreq;
+
+  public:
+    ClassDef(R3BTCalEngine, 1)
+};
+
+#endif