FixedTargetGenerator.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 "FixedTargetGenerator.h"
 
#include <TGeoManager.h>
 
#include <array>
#include <cmath>
 
#include "BeamSmearingUtils.h"
#include "FairMCEventHeader.h"
#include "FairPrimaryGenerator.h"
#include "HNLPythia8Generator.h"
#include "MeanMaterialBudget.h"
#include "Pythia8Plugins/EvtGen.h"
#include "ShipUnit.h"
#include "TGeoBBox.h"
#include "TGeoNode.h"
#include "TGeoVolume.h"
#include "TMCProcess.h"
#include "TMath.h"
#include "TROOT.h"
 
using ShipUnit::cm;
using ShipUnit::mm;
const Double_t c_light = 2.99792458e+10;  // speed of light in cm/sec (c_light
                                          // = 2.99792458e+8 * m/s)
const Double_t mbarn = 1E-3 * 1E-24 * TMath::Na();  // cm^2 * Avogadro
 
// -----   Default constructor   -------------------------------------------
FixedTargetGenerator::FixedTargetGenerator() {
  fUseRandom1 = kFALSE;
  fUseRandom3 = kTRUE;
  fSeed = 0;
  fMom = 400;  // proton
  fLogger = FairLogger::GetLogger();
  targetName = "cave_1/target_vacuum_box_1/TargetArea_1/HeVolume_1";
  xOff = 0;
  yOff = 0;
  zOff = 0;
  fsmearBeam = 8 * mm;  // default value for smearing beam (8 mm)
  fPaintBeam = 5 * cm;  // default value for painting beam (5 cm)
  tauOnly = false;
  JpsiMainly = false;
  DrellYan = false;
  PhotonCollision = false;
  G4only = false;
  OnlyMuons = false;
  maxCrossSection = 0.;
  EMax = 0;
  fBoost = 1.;
  withEvtGen = kFALSE;
  withNtuple = kFALSE;
  chicc = 1.7e-3;      // prob to produce primary ccbar pair/pot
  chibb = 1.6e-7;      // prob to produce primary bbbar pair/pot
  nrpotspill = 5.E13;  // nr of protons / spill
  setByHand = kFALSE;
  Debug = kFALSE;
  targetFromGeometry = kFALSE;
  Option = "Primary";
  wspill = 1.;  // event weight == 1 for primary events
  heartbeat = 1000;
}
Bool_t FixedTargetGenerator::InitForCharmOrBeauty(TString fInName, Int_t nev,
                                                  Double_t npot, Int_t nStart) {
  Option = "charm";
  nEntry = nStart;
  // open input file with charm or beauty
  fin = TFile::Open(fInName);
  if (!fin) {
    LOG(fatal) << "Could not open input file: " << fInName.Data();
    return kFALSE;
  }
  nTree = dynamic_cast<TNtuple*>(
      fin->FindObjectAny("pythia6"));  // old format, simple ntuple
  nEvents = nTree->GetEntries();
  nTree->SetBranchAddress("id", &n_id);
  nTree->SetBranchAddress("px", &n_px);
  nTree->SetBranchAddress("py", &n_py);
  nTree->SetBranchAddress("pz", &n_pz);
  nTree->SetBranchAddress("M", &n_M);
  nTree->SetBranchAddress("E", &n_E);
  nTree->SetBranchAddress("mid", &n_mid);
  nTree->SetBranchAddress("mpx", &n_mpx);
  nTree->SetBranchAddress("mpy", &n_mpy);
  nTree->SetBranchAddress("mpz", &n_mpz);
  nTree->SetBranchAddress("mE", &n_mE);
  if (nTree->GetBranch("k")) {
    LOG(info) << "+++has branch+++";
    nTree->SetBranchAddress("k", &ck);
  }
  // check if we deal with charm or beauty:
  nTree->GetEvent(0);
  if (!setByHand && n_M > 5) {
    chicc = chibb;
    LOG(info) << "automatic detection of beauty, configured for beauty";
    LOG(info) << "bb cross section / mbias " << chicc;
  } else {
    LOG(info) << "cc cross section / mbias " << chicc;
  }
  // convert pot to weight corresponding to one spill of 5e13 pot
  // get histogram with number of pot to normalise
  // pot are counted double, i.e. for each signal, i.e. pot/2.
  auto* potHist = dynamic_cast<TH1F*>(fin->Get("2"));
  if (!potHist) {
    LOG(fatal) << "Histogram '2' not found in input file";
    return kFALSE;
  }
  Int_t nrcpot =
      potHist->GetBinContent(1) / 2.;  // number of primary interactions
  wspill = nrpotspill * chicc / nrcpot * nEvents / nev;
  LOG(info) << "Input file: " << fInName.Data() << " with " << nEvents
            << " entries, corresponding to nr-pot=" << (nrcpot / chicc);
  LOG(info) << "weight " << wspill << " corresponding to " << nrpotspill
            << " p.o.t. per spill for " << nev << " events to process";
 
  pot = 0.;
  // Determine fDs on this file for primaries
  nDsprim = 0;
  ntotprim = 0;
 
  return kTRUE;
}
// -----   Default Init   -------------------------------------------
Bool_t FixedTargetGenerator::Init() {
  fPythiaP =
      new Pythia8::Pythia();  // pythia needed also for G4only, to copy 2mu BRs
  if (Option == "Primary" && !G4only) {
    fPythiaN = new Pythia8::Pythia();
  } else if (Option != "charm" && Option != "beauty" && !G4only) {
    LOG(error) << "Option not known " << Option.Data() << ", abort";
  }
  if (fUseRandom1) fRandomEngine = std::make_shared<PyTr1Rng>();
  if (fUseRandom3) fRandomEngine = std::make_shared<PyTr3Rng>();
  std::vector<int> r = {221, 221, 223, 223, 113, 331, 333};
  std::vector<int> c = {6, 7, 5, 7, 5, 6, 9};  // decay channel mumu mumuX
 
  if (Option == "Primary" && !G4only) {
    fPythiaP->settings.mode("Beams:idB", 2212);
    fPythiaN->settings.mode("Beams:idB", 2112);
  } else {
    fPythiaP->readString("ProcessLevel:all = off");
  }
  std::array<Pythia8::Pythia*, 2> plist = {fPythiaP, fPythiaN};
  Int_t pcount = 0;
  for (const auto& fPythia : plist) {
    if (pcount > 0 && (Option != "Primary" || G4only)) {
      continue;
    }
    pcount += 1;
    fPythia->setRndmEnginePtr(fRandomEngine);
    fPythia->settings.mode("Random:seed", fSeed);
    fPythia->settings.mode("Next:numberCount", heartbeat);
    if (Option == "Primary") {
      fPythia->settings.mode("Beams:idA", 2212);
      fPythia->settings.mode("Beams:frameType", 2);
      fPythia->settings.parm("Beams:eA", fMom);  // codespell:ignore parm
      fPythia->settings.parm("Beams:eB", 0.);    // codespell:ignore parm
    }
    if (JpsiMainly) {
      // use this for all onia productions
      fPythia->readString("Onia:all(3S1) = on");
      fPythia->readString(
          "443:new  J/psi  J/psi  3   0   0    3.09692    0.00009    3.09602   "
          " 3.09782  0.   1   1   0   1   0");
      fPythia->readString("443:addChannel = 1   1.    0      -13       13");
      fPythia->readString(
          "553:new  Upsilon  Upsilon  3   0   0    9.46030    0.00005    "
          "9.45980    9.46080  0.00000e+00   0   1   0   1   0");
      fPythia->readString("553:addChannel = 1   1.    0      -13       13");
    }
    if (DrellYan) {
      fPythia->readString("WeakSingleBoson:ffbar2gmZ = on");
      //   Z0/gamma  -> mu+ mu-
      fPythia->readString("23:onMode = off");
      fPythia->readString("23:onIfAll = 13 13");
      fPythia->readString("23:mMin = 0.5");
      fPythia->readString("PhaseSpace:mHatMin = 0.5");
      fPythia->readString("PartonLevel:FSR = on");
      fPythia->readString("PDF:pSet = 4");
    }
    if (PhotonCollision) {
      fPythia->readString("PhotonCollision:gmgm2mumu = on");
    }
    if (!DrellYan && !PhotonCollision) {
      fPythia->readString("SoftQCD:inelastic = on");
      fPythia->readString("PhotonCollision:gmgm2mumu = on");
      fPythia->readString("PromptPhoton:all = on");
      fPythia->readString("WeakBosonExchange:all = on");
      fPythia->readString("WeakSingleBoson:all = on");
    }
    if (tauOnly) {
      fPythia->readString(
          "431:new  D_s+  D_s-  1   3   0    1.96849    0.00000    0.00000    "
          "0.00000  1.49900e-01   0   1   0   1   0");
      fPythia->readString(
          "431:addChannel = 1   0.0640000    0      -15       16");
    }
 
    // find all long lived particles in pythia
    Int_t n = 1;
    while (n != 0) {
      n = fPythia->particleData.nextId(n);
      std::shared_ptr<Pythia8::ParticleDataEntry> p =
          fPythia->particleData.particleDataEntryPtr(n);
      if (p->tau0() > 1) {
        std::string particle = std::to_string(n) + ":mayDecay = false";
        fPythia->readString(particle);
        LOG(info) << "Made " << p->name().c_str()
                  << " stable for Pythia, should decay in Geant4";
      }
    }
    // boost branching fraction of rare di-muon decays
    //                       eta  omega rho0  eta' phi
    if (fBoost != 1.) {
      LOG(info) << "Rescale BRs of dimuon decays in Pythia: " << fBoost;
      for (unsigned int i = 0; i < r.size(); ++i) {
        std::shared_ptr<Pythia8::ParticleDataEntry> V =
            fPythia->particleData.particleDataEntryPtr(r[i]);
        Pythia8::DecayChannel ch = V->channel(c[i]);
        if (TMath::Abs(ch.product(0)) != 13 ||
            TMath::Abs(ch.product(1)) != 13) {
          LOG(info) << "this is not the right decay channel: " << r[i] << " "
                    << c[i];
        } else {
          TString tmp = "";
          tmp += r[i];
          tmp += ":";
          tmp += c[i];
          tmp += ":bRatio =";
          tmp += fBoost * ch.bRatio();
          fPythia->readString(tmp.Data());
        }
      }
    }
    fPythia->init();
  }
  // Initialize EvtGen.
  if (withEvtGen) {
    // Use EVTGENDATA environment variable to find EvtGen data files
    const char* evtgendata = getenv("EVTGENDATA");
    if (!evtgendata) {
      LOG(fatal) << "EVTGENDATA environment variable not set";
    }
    TString DecayFile = TString(evtgendata) + "/DECAY.DEC";
 
    TString ParticleFile = TString(evtgendata) + "/evt.pdl";
    std::cout << "Using $EVTGENDATA " << evtgendata << std::endl;
    EvtExternalGenList* extPtr = new EvtExternalGenList();
    std::list<EvtDecayBase*> models = extPtr->getListOfModels();
    TString UdecayFile = getenv("FAIRSHIP");
    UdecayFile += "/gconfig/USERDECAY.DEC";
    evtgenP = new Pythia8::EvtGenDecays(fPythiaP, DecayFile.Data(),
                                        ParticleFile.Data(), extPtr);
    evtgenP->readDecayFile(
        UdecayFile.Data());  // will make update of EvtGen with user decay file
    // use one instance of EvtGen, requires patch to Pythia8Plugins/EvtGen.h
    if (Option == "Primary") {
      evtgenN =
          new Pythia8::EvtGenDecays(fPythiaN, DecayFile.Data(), "", extPtr);
    }
  }
  if (targetFromGeometry) {
    // Use geometry-based coordinates passed from run_simScript.py
    fMaterialInvestigator = new GenieGenerator();
    start[0] = xOff;
    start[1] = yOff;
    start[2] = startZ + zOff;
    end[0] = xOff;
    end[1] = yOff;
    end[2] = endZ;
    LOG(info) << "FixedTargetGenerator: Using geometry-based target "
                 "coordinates startZ="
              << startZ << " endZ=" << endZ;
    // find maximum interaction length
    bparam = shipgen::MeanMaterialBudget(start, end, mparam);
    maxCrossSection = mparam[9];
  } else if (targetName != "") {
    // Fallback to fragile TGeo navigation for backward compatibility
    fMaterialInvestigator = new GenieGenerator();
    TGeoNavigator* nav = gGeoManager->GetCurrentNavigator();
    if (nav->CheckPath(targetName)) {
      nav->cd(targetName);
    } else {
      LOG(fatal) << "Invalid target volume specified";
    }
    TGeoNode* target = nav->GetCurrentNode();
    TObjArray* nodes = target->GetVolume()->GetNodes();
    // Get the first and last node of the target to calculate the material seen
    TGeoNode* first = static_cast<TGeoNode*>(nodes->At(0));
    TGeoNode* last = static_cast<TGeoNode*>(nodes->At(nodes->GetSize() - 1));
    nav->cd(targetName + "/" + first->GetName());
    TGeoBBox* sha = static_cast<TGeoBBox*>(first->GetVolume()->GetShape());
    Double_t dz = sha->GetDZ();
    Double_t origin[3] = {0, 0, -dz};
    Double_t master[3] = {0, 0, 0};
    nav->LocalToMaster(origin, master);
    startZ = master[2];
    nav->cd(targetName + "/" + last->GetName());
    sha = static_cast<TGeoBBox*>(first->GetVolume()->GetShape());
    dz = sha->GetDZ();
    origin[2] = +dz;
    nav->LocalToMaster(origin, master);
    endZ = master[2];
    start[0] = xOff;
    start[1] = yOff;
    start[2] = startZ + zOff;
    end[0] = xOff;
    end[1] = yOff;
    end[2] = endZ;
    // find maximum interaction length
    bparam = shipgen::MeanMaterialBudget(start, end, mparam);
    maxCrossSection = mparam[9];
  } else {
    LOG(fatal) << "No target set.";
  }
 
  return kTRUE;
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
FixedTargetGenerator::~FixedTargetGenerator() {}
// -------------------------------------------------------------------------
 
// -----   Passing the event   ---------------------------------------------
Bool_t FixedTargetGenerator::ReadEvent(FairPrimaryGenerator* cpg) {
  // Calculate beam smearing and painting
  auto [dx, dy] = CalculateBeamOffset(fsmearBeam, fPaintBeam);
 
  Double_t zinter = 0;
  Double_t ZoverA = 1.;
  if (!targetName.IsNull()) {
    // calculate primary proton interaction point:
    // loop over trajectory between start and end to pick an interaction point,
    // copied from GenieGenerator and adapted to hadrons
    Double_t prob2int = -1.;
    Double_t rndm = 0.;
    Double_t sigma;
    Double_t zinterStart = start[2];
    if (Option == "charm" || Option == "beauty") {
      // simulate more downstream interaction points for interactions down in
      // the cascade
      if (!(nTree->GetBranch("k"))) {
        ck = 1;
      }
    } else {
      ck = 1;
    }
    while (ck > 0.5) {
      while (prob2int < rndm) {
        // place x,y,z uniform along path
        zinter = gRandom->Uniform(zinterStart, end[2]);
        Double_t point[3] = {xOff, yOff, zinter};
        bparam = shipgen::MeanMaterialBudget(start, point, mparam);
        Double_t interLength = mparam[8];
        TGeoNode* node = gGeoManager->FindNode(point[0], point[1], point[2]);
        TGeoMaterial* mat = nullptr;
        if (node && !gGeoManager->IsOutside()) {
          mat = node->GetVolume()->GetMaterial();
          Double_t n = mat->GetDensity() / mat->GetA();
          ZoverA = mat->GetZ() / mat->GetA();
          sigma = 1. / (n * mat->GetIntLen()) / mbarn;
          prob2int = TMath::Exp(-interLength) * sigma / maxCrossSection;
        } else {
          prob2int = 0.;
        }
        rndm = gRandom->Uniform(0., 1.);
      }
      zinterStart = zinter;
      ck -= 1;
    }
    zinter = zinter * cm;
  }
  Pythia8::Pythia* fPythia;
  if (G4only) {
    cpg->AddTrack(2212, 0., 0., fMom, (xOff + dx) / cm, (yOff + dy) / cm,
                  start[2], -1, kTRUE, -1., 0., 1.);
    return kTRUE;
  } else if (Option == "Primary") {
    if (gRandom->Uniform(0., 1.) < ZoverA) {
      fPythiaP->next();
      if (withEvtGen) {
        evtgenP->decay();
      }
      fPythia = fPythiaP;
    } else {
      fPythiaN->next();
      if (withEvtGen) {
        evtgenN->decay();
      }
      fPythia = fPythiaN;
    }
  } else {
    if (nEntry == nEvents) {
      LOG(info) << "Rewind input file: " << nEntry;
      nEntry = 0;
    }
    nTree->GetEvent(nEntry);
    nEntry += 1;
    // sanity check, count number of p.o.t. on input file.
    Double_t pt = TMath::Sqrt((n_mpx * n_mpx) + (n_mpy * n_mpy));
    // every event appears twice, i.e.
    if (pt < 1.e-5 && n_mid == 2212) {
      pot += +0.5;
      ntotprim += 1;
    }
    Int_t idabs = static_cast<int>(TMath::Abs(n_id));
    if (idabs == 431) {
      nDsprim += 1;
    }
    fPythiaP->event.reset();
    Int_t nID1 = n_id;
    fPythiaP->event.append(static_cast<int>(n_id), 1, 0, 0, n_px, n_py, n_pz,
                           n_E, n_M, 0., 9.);
    TMCProcess procID = kPTransportation;
    if (n_mid == 2212 && (n_mpx * n_mpx + n_mpy * n_mpy) < 1E-5) {
      procID = kPPrimary;
    }  // probably primary and not from cascade
    cpg->AddTrack(static_cast<int>(n_mid), n_mpx, n_mpy, n_mpz,
                  (xOff + dx) * cm, (yOff + dy) * cm, zinter * cm, -1, kFALSE,
                  n_mE, 0., wspill, procID);
    // second charm hadron in the event
    nTree->GetEvent(nEntry);
    if (nID1 * n_id > 0) {
      LOG(info) << "same sign charm: " << nEntry << ", " << nID1 << ", "
                << n_id;
    }
    nEntry += 1;
    fPythiaP->event.append(static_cast<int>(n_id), 1, 0, 0, n_px, n_py, n_pz,
                           n_E, n_M, 0., 9.);
    fPythiaP->next();
    fPythia = fPythiaP;
  }
  if (withEvtGen) {
    fPythia->moreDecays();
  }  // let the very short lived resonances decay via Pythia8
  if (Debug && !G4only) {
    fPythia->event.list();
  }
  TMCProcess procID;
  for (Int_t ii = 1; ii < fPythia->event.size(); ii++) {
    Double_t e = fPythia->event[ii].e();
    Double_t m = fPythia->event[ii].m();
    Double_t pz = fPythia->event[ii].pz();
    Int_t id = fPythia->event[ii].id();
    Bool_t wanttracking = kTRUE;
    if (e - m < EMax || !fPythia->event[ii].isFinal() || pz < 0) {
      wanttracking = kFALSE;
    }
    if (DrellYan || PhotonCollision || OnlyMuons) {
      // don't track underlying event
      if (fabs(id) != 13) {
        wanttracking = kFALSE;
      }
    }
    Double_t z = fPythia->event[ii].zProd() * mm + zinter * cm;
    Double_t x = fPythia->event[ii].xProd() * mm + xOff * cm + dx * cm;
    Double_t y = fPythia->event[ii].yProd() * mm + yOff * cm + dy * cm;
    Double_t tof =
        fPythia->event[ii].tProd() / (10 * c_light);  // to go from mm to s
    Double_t px = fPythia->event[ii].px();
    Double_t py = fPythia->event[ii].py();
    Int_t im = fPythia->event[ii].mother1() - 1;
    procID = kPPrimary;
    if (Option != "Primary") {
      procID = kPDecay;
      im += 1;
      if (ii == 1) {
        procID = kPHadronic;
      }
    } else {
      if (ii < 3) {
        im = -1;
      }
    }
    cpg->AddTrack(id, px, py, pz, x * cm, y * cm, z * cm, im, wanttracking, e,
                  tof, wspill, procID);
    if (withNtuple) {
      int idabs = TMath::Abs(id);
      if (idabs < 10) {
        continue;
      }
      if (idabs < 18 || idabs == 22 || idabs == 111 || idabs == 221 ||
          idabs == 223 || idabs == 331 || idabs == 211 || idabs == 113 ||
          idabs == 333 || idabs == 321 || idabs == 2212) {
        Int_t moID = fPythia->event[im + 1].id();
        fNtuple->Fill(0, id, px, py, pz, e, x * cm, y * cm, z * cm, moID);
      }
    }
  }
  return kTRUE;
}
// -------------------------------------------------------------------------