BSEmcMCHitToTimebasedWaveforms.h

Go to the documentation of this file.

//----------------------------------------------------------------------
// Author List:
//      Phil Strother                   Original author
//  Dima Melnichuk          - adaption for PANDA
//  Philipp Mahlberg        - integration in timebased simulation concept
//----------------------------------------------------------------------

//#pragma once
#ifndef BSEMCMCHITTOTIMEBASEDWAVEFORMS_HH
#define BSEMCMCHITTOTIMEBASEDWAVEFORMS_HH

#include <PndPersistencyTask.h>
#include <cassert>
#include <string>
#include <vector>

#include "TClonesArray.h"
#include "TROOT.h"
#include "TRandom.h"
#include "TStopwatch.h"
#include "TString.h"

#include "FairLink.h"
#include "FairLogger.h"
#include "FairRootManager.h"
#include "FairRun.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "PndTCAInputContainer.h"

#include "BSEmcAbsWaveformSimulator.h"
#include "BSEmcDigiPar.h"
#include "BSEmcMCHit.h"
#include "BSEmcMCHitToTimebasedWaveforms.h"
#include "BSEmcWaveform.h"
#include "BSEmcWaveformBuffer.h"
#include "BSEmcWaveformData.h"

template <class ParSet>
class BSEmcMCHitToTimebasedWaveforms : public PndPersistencyTask {
 public:
  // Constructors
  BSEmcMCHitToTimebasedWaveforms(const std::string &t_detectorname, Bool_t t_storewaves = kTRUE)
    : PndPersistencyTask("BSEmcMCHitToTimebasedWaveforms"), fDetectorName(t_detectorname)
  {
    SetPersistency(t_storewaves);
  }

  // Destructor
  virtual ~BSEmcMCHitToTimebasedWaveforms(){};

  virtual InitStatus Init() /*override*/
  {

    fHitBranchName = BSEmcDataBranchNames::fgMCHitBranchName + fDetectorName;
    if (fWaveformBranchName == "") {
      fWaveformBranchName = BSEmcDataBranchNames::fgMultiWaveformBranchName + this->fDetectorName;
    }
    if (fWaveformDataBranchName == "") {
      fWaveformDataBranchName = BSEmcDataBranchNames::fgWaveformDataBranchName + this->fDetectorName;
    }
    fHitArray.Init(fHitBranchName);

    if (fSimulator == nullptr) {
      LOG(debug) << "BSEmcMCHitToTimebasedWaveforms for " << fHitBranchName << " has no Simulator set. Calling SetupSimulator!";
      SetupSimulator();
      fSimulator->SetHitBranchName(fHitBranchName);
    }
    if (fSimulator == nullptr) {
      LOG(error) << "BSEmcMCHitToTimebasedWaveforms for " << fHitBranchName << " has no Simulator set. Aborting!";
      exit(-1);
    }

    // Get RootManager
    FairRootManager *ioman = FairRootManager::Instance();
    if (ioman == nullptr) {
      LOG(error) << "BSEmcMCHitToTimebasedWaveforms<ParSet>::Init: "
                 << "RootManager not instantiated!";
      return kFATAL;
    }

    // Create and activiate output Buffer....choose between BSEmcWaveform and BSEmcMultiWaveform
    fWaveformBuffer = new BSEmcWaveformBuffer(fWaveformBranchName, "BSEmcMultiWaveform", "Emc", GetPersistency());
    LOG(debug) << "Using BSEmcMultiWaveform in WaveformBuffer.";

    fWaveformBuffer = (BSEmcWaveformBuffer *)ioman->RegisterWriteoutBuffer(fWaveformBranchName, fWaveformBuffer);
    fWaveformBuffer->ActivateBuffering(fActivateBuffering);

    if (fStoreDataClass) {
      fWaveformBuffer->StoreWaveformData(fWaveformDataBranchName, "Emc", kTRUE);
    }

    fUse_photon_statistic = fPhotonStatisticPar->GetUseDigiEffectiveSmearingMode() == 2;

    if (fUse_photon_statistic) {
      LOG(debug) << "BSEmcMCHitToTimebasedWaveforms: " << fDetectorName << " using photon statistic";
      fNPhotoElectronsPerMeV = fPhotonStatisticPar->GetDetectedPhotonsPerMeV() * fPhotonStatisticPar->GetSensitiveArea() / fPhotonStatisticPar->GetRearCrystalSurfaceArea() *
                               fPhotonStatisticPar->GetQuantumEfficiency();
      fExcessNoiseFactor = fPhotonStatisticPar->GetExcessNoiseFactor();
    } else {
      fNPhotoElectronsPerMeV = 0;
      fExcessNoiseFactor = 1;
    }
    fOverlapTime = fPhotonStatisticPar->GetOverlapTime();
    LOG(debug) << "BSEmcMCHitToTimebasedWaveforms::Init() for " << fHitBranchName << " done.";
    return kSUCCESS;
  }

  virtual void Exec(Option_t */*t_opt*/) /*override*/
  {
    LOG(debug) << "BSEmcMCHitToTimebasedWaveforms<ParSet> " << fDetectorName << " Exec() ";
    FairRootManager *ioman = FairRootManager::Instance();

    TStopwatch timer;
    timer.Start();

    const BSEmcMCHit *theHit = nullptr;

    Int_t nHits = fHitArray.GetSize();
    LOG(debug) << "BSEmcMCHitToTimebasedWaveforms<ParSet> " << fDetectorName << " Hit array contains " << nHits << " hits";

    for (Int_t iHit = 0; iHit < nHits; iHit++) {

      theHit = fHitArray.GetConstElementPtr(iHit);
      Double_t energy = theHit->GetEnergy();
      if (energy == 0) {
        continue;
      }
      LOG(DEBUG3) << "hit energy : " << energy;
      if (fUse_photon_statistic) {
        Double_t crystalPhotonsMeV = 1.0e3 * energy * fNPhotoElectronsPerMeV;
        Double_t photonStatFactor = gRandom->Gaus(1, sqrt(fExcessNoiseFactor / crystalPhotonsMeV));
        LOG(DEBUG3) << "photonStatFactor: " << photonStatFactor;
        energy *= photonStatFactor;
        LOG(DEBUG3) << "energy now: " << energy;
        // energy *= 1e3;
        LOG(DEBUG3) << "energy now in MeV: " << energy;
      }
      // construct corresponding waveform data Object
      BSEmcWaveformData wfData(theHit->GetDetectorID(), fSimulator);
      wfData.SetOverlapTime(fOverlapTime);
      // register hit...timebased framework uses ns, whereas emc deals with seconds as time unit
      FairLink linkToHit(-1, ioman->GetEntryNr(), fHitBranchName, iHit, 1.0);
      if (fActivateBuffering) {
        wfData.AddHit(linkToHit, ioman->GetEventTime() + theHit->GetTime() * 1.0e9, energy);
      } else {
        wfData.AddHit(linkToHit, theHit->GetTime() * 1.0e9, energy);
      }
      LOG(debug) << "BSEmcMCHitToTimebasedWaveforms " << fDetectorName << " adding hit with event time " << ioman->GetEventTime() << " and hit time: " << theHit->GetTime() * 1.0e9
                 << " and energy: " << energy;
      fWaveformBuffer->FillNewData(&wfData);
    }
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    LOG(debug) << "BSEmcMCHitToTimebasedWaveforms, Real time " << rtime << " s, CPU time " << ctime << " s";
  }

  void RunTimebased(Bool_t t_timebased = kTRUE) { fActivateBuffering = t_timebased; };
  void SetSimulator(BSEmcAbsWaveformSimulator *t_simulator) { fSimulator = t_simulator; }
  void SetStorageOfData(Bool_t t_storeWaves = kTRUE) { SetPersistency(t_storeWaves); }; 
  void StoreDataClass(Bool_t t_storeData = kTRUE) { fStoreDataClass = t_storeData; };

  void SetWaveformBranchName(const TString &t_branchName) { fWaveformBranchName = t_branchName; }
  void SetWaveformDataBranchName(const TString &t_branchName) { fWaveformDataBranchName = t_branchName; }
  void SetPhotonStatisticParName(const std::string &t_photonStatParName) { fPhotonStatisticParName = t_photonStatParName; }
  void SetWaveformGenParName(const std::string &t_waveformGenParName) { fWaveformGenParName = t_waveformGenParName; }

 protected:
  virtual void SetParContainers() /*override*/
  {
    // Get run and runtime database
    FairRun *run = FairRun::Instance();
    if (run == nullptr) {
      Fatal("SetParContainers", "No analysis run");
    }
    FairRuntimeDb *db = run->GetRuntimeDb();
    if (db == nullptr) {
      Fatal("SetParContainers", "No runtime database");
    }
    // Get Emc digitisation parameter container
    if (fWaveformGenParName == "") {
      fWaveformGenParName = ParSet::fgParameterName;
    }
    fDigiPar = dynamic_cast<ParSet *>(db->getContainer(fWaveformGenParName.c_str()));

    if (fPhotonStatisticParName == "") {
      fPhotonStatisticParName = BSEmcDigiPar::fgParameterName + fDetectorName;
    }
    fPhotonStatisticPar = dynamic_cast<BSEmcDigiPar *>(db->getContainer(fPhotonStatisticParName.c_str()));
  }
  virtual void SetupSimulator() = 0;

  TString fDetectorName{""};
  TString fHitBranchName{""};
  TString fWaveformBranchName{""};
  TString fWaveformDataBranchName{""};
  PndTCAInputContainer<BSEmcMCHit> fHitArray{};
  BSEmcWaveformBuffer *fWaveformBuffer{nullptr};

  Bool_t fStoreDataClass{kFALSE};
  Bool_t fActivateBuffering{kFALSE};
  std::string fPhotonStatisticParName{""};
  BSEmcDigiPar *fPhotonStatisticPar{nullptr};

  std::string fWaveformGenParName{""};
  ParSet *fDigiPar{nullptr}; 
  Bool_t fUse_photon_statistic{kFALSE};
  Double_t fNPhotoElectronsPerMeV{0};
  Double_t fExcessNoiseFactor{0};
  Double_t fOverlapTime{0};
  BSEmcAbsWaveformSimulator *fSimulator{nullptr};

  ClassDef(BSEmcMCHitToTimebasedWaveforms, 2);
};

#endif /*BSEMCMCHITTOTIMEBASEDWAVEFORMS_HH*/