NtupleGenerator.cxx

Go to the documentation of this file.

// SPDX-License-Identifier: LGPL-3.0-or-later
// SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP
// Collaboration
 
#include "NtupleGenerator.h"
 
#include <cmath>
 
#include "FairPrimaryGenerator.h"
#include "TDatabasePDG.h"  // for TDatabasePDG
#include "TFile.h"
#include "TMath.h"  // for Sqrt
#include "TROOT.h"
 
using std::cout;
using std::endl;
// read events from ntuples produced
 
// -----   Default constructor   -------------------------------------------
NtupleGenerator::NtupleGenerator() {}
// -------------------------------------------------------------------------
// -----   Default constructor   -------------------------------------------
Bool_t NtupleGenerator::Init(const char* fileName) { return Init(fileName, 0); }
// -----   Default constructor   -------------------------------------------
Bool_t NtupleGenerator::Init(const char* fileName, const int startEvent) {
  cout << "Info NtupleGenerator: Opening input file " << fileName << endl;
  fInputFile = new TFile(fileName);
  if (fInputFile->IsZombie()) {
    cout << "-E NtupleGenerator: Error opening the Signal file" << fileName
         << endl;
  }
  fTree = dynamic_cast<TTree*>(fInputFile->Get("ntuple"));
  fNevents = fTree->GetEntries();
  fn = startEvent;
  fTree->SetBranchAddress("id", &id);  // particle id
  if (fTree->FindBranch("parentid")) {
    fTree->SetBranchAddress("parentid", &parentid);
  }  // parent id
  if (fTree->FindBranch("tof")) {
    fTree->SetBranchAddress("tof", &tof);
  }  // time of flight
  fTree->SetBranchAddress("Nmeas", &Nmeas);  // number of Geant4 points
  fTree->SetBranchAddress("Ezero", &Ezero);  // incoming muon energy
  fTree->SetBranchAddress("w", &w);          // weight of event
  fTree->SetBranchAddress("x", &vx);         // position
  fTree->SetBranchAddress("y", &vy);
  fTree->SetBranchAddress("z", &vz);
  fTree->SetBranchAddress("px", &px);  // momentum
  fTree->SetBranchAddress("py", &py);
  fTree->SetBranchAddress("pz", &pz);
  fTree->SetBranchAddress("volid", &volid);    // which volume
  fTree->SetBranchAddress("procid", &procid);  // which process
  return kTRUE;
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
NtupleGenerator::~NtupleGenerator() {
  // cout << "destroy Ntuple" << endl;
  fInputFile->Close();
  fInputFile->Delete();
  delete fInputFile;
}
// -------------------------------------------------------------------------
 
// -----   Passing the event   ---------------------------------------------
Bool_t NtupleGenerator::ReadEvent(FairPrimaryGenerator* cpg) {
  while (fn < fNevents) {
    fTree->GetEntry(fn);
    fn++;
    if (fn % 10000 == 0) {
      cout << "reading event " << fn << endl;
    }
    // test if muon survives:
    Int_t i = Nmeas - 3;
    Float_t r2 = (vx[i] * vx[i] + vy[i] * vy[i]);
    if (procid[Nmeas - 1] == 2 && r2 < 9) {
      break;
    }
  }
  if (fn == fNevents) {
    cout << "No more input events" << endl;
    return kFALSE;
  }
  TDatabasePDG* pdgBase = TDatabasePDG::Instance();
  Double_t mass = pdgBase->GetParticle(id)->Mass();
  Double_t e =
      TMath::Sqrt(px[0] * px[0] + py[0] * py[0] + pz[0] * pz[0] + mass * mass);
  tof = 0;
  // first, original muon
  cpg->AddTrack(id, px[0], py[0], pz[0], vx[0] * 100., vy[0] * 100.,
                vz[0] * 100., -1., false, e, tof, w);
  Int_t i = Nmeas - 1;
  // second, surviving muon, extrapolate back to end of muon shield, z=20m
  Double_t zscor = 20.;
  Double_t lam = (zscor - vz[i]) / pz[i];
  Double_t xscor = vx[i] + lam * px[i];
  Double_t yscor = vy[i] + lam * py[i];
  e = TMath::Sqrt(px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] + mass * mass);
  cpg->AddTrack(id, px[i], py[i], pz[i], xscor * 100., yscor * 100.,
                zscor * 100., 0, true, e, tof, w);
  return kTRUE;
}
 
// -------------------------------------------------------------------------
Int_t NtupleGenerator::GetNevents() { return fNevents; }