eclComputePulseTemplates_Step1.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#include <TF1.h>
#include <ecl/digitization/OfflineFitFunction.h>
#include <TTree.h>
#include <TFile.h>
#include <iostream>
#include <assert.h>
 
//
/*
 See eclComputePulseTemplates_Step0.cc for README instructions.
 */
//
// cppcheck-suppress constParameter ; main() keeps its standard signature
int main(int argc, char* argv[])
{
  //
  TString OutputDirectory = "";
  if (OutputDirectory == ".") {
    std::cout << "Error set output directory" << std::endl;
    return -1;
  }
  //
  assert(argc == 3);
  int LowIDLimit = atoi(argv[1]);
  int HighIDLimit = atoi(argv[2]);
  //
  double PhotonWaveformPar[11];
  TFile* PhotonParFile = new TFile(OutputDirectory + "PhotonWaveformParameters.root");
  TTree* chain = static_cast<TTree*>(PhotonParFile->Get("ParTree"));
  chain->SetBranchAddress("PhotonPar", &PhotonWaveformPar);
  //
  TFile* f = new TFile(OutputDirectory + Form("PhotonShapes_Low%d_High%d.root", LowIDLimit, HighIDLimit), "RECREATE");
  f->cd();
  TTree* mtree = new TTree("mtree", "");
  std::vector<double> PhotonWaveformArray(100000);
  mtree->Branch("PhotonArray", PhotonWaveformArray.data(), "PhotonWaveformArray[100000]/D");
 
  //
  for (Long64_t jentry = LowIDLimit; jentry < HighIDLimit; jentry++) {
    chain->GetEntry(jentry);
    TF1 PhotonShapeFunc = TF1(Form("photonShape_%lld", jentry), Belle2::ECL::WaveFuncTwoComponent, 0, 20, 26);;
    PhotonShapeFunc.SetNpx(1000);
    std::cout << PhotonWaveformPar[0] << std::endl;
    PhotonShapeFunc.SetParameter(0, 0);
    PhotonShapeFunc.SetParameter(1, 0);
    PhotonShapeFunc.SetParameter(2, 1);
    PhotonShapeFunc.SetParameter(3, 0);
    for (int k = 0; k < 10; k++) {
      PhotonShapeFunc.SetParameter(4 + k, PhotonWaveformPar[k + 1]);
      PhotonShapeFunc.SetParameter(10 + 4 + k, PhotonWaveformPar[k + 1]);
    }
    PhotonShapeFunc.SetParameter(24, PhotonWaveformPar[0]);
    PhotonShapeFunc.SetParameter(25, 1);
    //
    if (PhotonWaveformPar[0] > 0 && jentry <= 8737) {
      for (unsigned int k = 0; k < PhotonWaveformArray.size();
           k++) PhotonWaveformArray[k] = PhotonShapeFunc.Eval(((double)k) * (1. / 1000.)) ;
    } else {
      for (unsigned int k = 0; k < PhotonWaveformArray.size(); k++) PhotonWaveformArray[k] = -999;
    }
    mtree->Fill();
    //
  }
  //
  f->cd();
  mtree->Write();
  f->Write();
  //
  return 0;
}