#include <MuonBackGenerator.h>

Inheritance diagram for MuonBackGenerator:

Collaboration diagram for MuonBackGenerator:

Public Member Functions
	MuonBackGenerator ()

	~MuonBackGenerator () override

Bool_t	ReadEvent (FairPrimaryGenerator *) override

Bool_t	Init (const char *, int) override

Bool_t	Init (const char *) override

Bool_t	Init (const std::vector< std::string > &, int) override

Bool_t	Init (const std::vector< std::string > &) override

Int_t	GetNevents ()

void	CloseFile ()

void	FollowAllParticles ()

void	SetSmearBeam (Double_t sb)

void	SetPaintRadius (Double_t r)

void	SetSameSeed (Int_t s)

void	SetPhiRandomize (Bool_t phiRandomize)

Bool_t	checkDiMuon (Int_t muIndex)

void	SetDownScaleDiMuon ()

Public Member Functions inherited from SHiP::Generator
	Generator ()

virtual	~Generator ()

virtual Bool_t	Init (const char *, int)=0

virtual Bool_t	Init (const char *)=0

virtual Bool_t	Init (const std::vector< std::string > &inFiles, int startNumber)

virtual Bool_t	Init (const std::vector< std::string > &inFiles)

virtual void	UseExternalFile (std::string x, Int_t i)

virtual void	UseExternalFile (std::vector< std::string > &inFiles, Int_t i)

Protected Attributes
Float_t	id

Float_t	parentid

Float_t	pythiaid

Float_t	w

Float_t	px

Float_t	py

Float_t	pz

Float_t	vx

Float_t	vy

Float_t	vz

Float_t	ecut

TClonesArray *	MCTrack

TClonesArray *	vetoPoints

std::vector< ShipMCTrack > *	MCTrack_vec

std::vector< vetoPoint > *	vetoPoints_vec

Bool_t	fUseSTL

TFile *	fInputFile
	flag to indicate if using STL vectors

TChain *	fTree

int	fNevents

float	f_zOffset

int	fn

Bool_t	fPhiRandomize

Double_t	fPaintBeam

Bool_t	fdownScaleDiMuon

Bool_t	followMuons

Int_t	fSameSeed

Double_t	fsmearBeam

Protected Attributes inherited from SHiP::Generator
std::optional< std::string >	fextFile

Int_t	firstEvent = 0

Detailed Description

Definition at line 20 of file MuonBackGenerator.h.

Constructor & Destructor Documentation

◆ MuonBackGenerator()

MuonBackGenerator::MuonBackGenerator ( )

default constructor

Definition at line 33 of file MuonBackGenerator.cxx.

    : MCTrack_vec(nullptr), vetoPoints_vec(nullptr), fUseSTL(false) {
  followMuons = true;
}

◆ ~MuonBackGenerator()

MuonBackGenerator::~MuonBackGenerator ( )

override

destructor

Definition at line 148 of file MuonBackGenerator.cxx.

                                      {
  if (!fUseSTL) {
    delete MCTrack;
    delete vetoPoints;
  }
}

Member Function Documentation

◆ checkDiMuon()

Bool_t MuonBackGenerator::checkDiMuon ( Int_t muIndex )

Definition at line 155 of file MuonBackGenerator.cxx.

                                                   {
  ShipMCTrack* mu = fUseSTL ? &(*MCTrack_vec)[muIndex]
                            : dynamic_cast<ShipMCTrack*>(MCTrack->At(muIndex));
  TString pName = mu->GetProcName();
  if (strncmp("Hadronic inelastic", pName.Data(), 18) == 0 ||
      strncmp("Positron annihilation", pName.Data(), 21) == 0 ||
      strncmp("Lepton pair production", pName.Data(), 22) == 0) {
    return true;
  }
  Int_t motherId = mu->GetMotherId();
  if (motherId < 0) {
    return false;
  }
  ShipMCTrack* mother = fUseSTL
                            ? &(*MCTrack_vec)[motherId]
                            : dynamic_cast<ShipMCTrack*>(MCTrack->At(motherId));
  Int_t Pcode = TMath::Abs(mother->GetPdgCode());
  return (Pcode == 221 || Pcode == 223 || Pcode == 333 || Pcode == 113 ||
          Pcode == 331);
}

◆ CloseFile()

void MuonBackGenerator::CloseFile ( )

Definition at line 333 of file MuonBackGenerator.cxx.

                                  {
  fInputFile->Close();
  fInputFile->Delete();
  delete fInputFile;
}

◆ FollowAllParticles()

void MuonBackGenerator::FollowAllParticles ( )

inline

Definition at line 37 of file MuonBackGenerator.h.

37{ followMuons = false; };

◆ GetNevents()

Int_t MuonBackGenerator::GetNevents ( )

Definition at line 332 of file MuonBackGenerator.cxx.

332{ return fNevents; }

MuonBackGenerator::fNevents

int fNevents

Definition: MuonBackGenerator.h:58

◆ Init() [1/4]

Bool_t MuonBackGenerator::Init ( const char * fileName )

overridevirtual

Implements SHiP::Generator.

Definition at line 39 of file MuonBackGenerator.cxx.

                                                   {
  return Init(fileName, 0);
}

◆ Init() [2/4]

Bool_t MuonBackGenerator::Init	(	const char *	fileName,
		int	startEvent
	)

overridevirtual

Implements SHiP::Generator.

Definition at line 143 of file MuonBackGenerator.cxx.

                                                                         {
  std::vector<std::string> fileNames = {fileName};
  return Init(fileNames, startEvent);
}

◆ Init() [3/4]

Bool_t MuonBackGenerator::Init ( const std::vector< std::string > & fileNames )

overridevirtual

Reimplemented from SHiP::Generator.

Definition at line 43 of file MuonBackGenerator.cxx.

                                                                    {
  return Init(fileNames, 0);
}

◆ Init() [4/4]

Bool_t MuonBackGenerator::Init	(	const std::vector< std::string > &	fileNames,
		int	startEvent
	)

overridevirtual

Reimplemented from SHiP::Generator.

Definition at line 47 of file MuonBackGenerator.cxx.

                                                     {
  LOG(info) << "Opening input file to find keys " << fileNames.at(0);
  TFile testFile(fileNames.at(0).c_str());
  auto testKeys = testFile.GetListOfKeys();
  if (testKeys == nullptr) {
    LOG(fatal) << "Error opening the Signal file: " << fileNames.at(0);
  }
  fn = startEvent;
  fPaintBeam = 5 * cm;  // default value for painting beam
  fSameSeed = 0;
  fPhiRandomize = false;      // default value for phi randomization
  fsmearBeam = 8 * mm;        // default value for smearing beam
  fdownScaleDiMuon = kFALSE;  // only needed for muflux simulation
  if (testKeys->FindObject("pythia8-Geant4")) {
    fTree = new TChain("pythia8-Geant4");
    for (auto& f : fileNames) {
      LOG(info) << "Opening input file " << f;
      fTree->Add(f.c_str());
    }
    fNevents = fTree->GetEntries();
    LOG(info) << "Reading " << fNevents << " entries";
    // count only events with muons
    fTree->SetBranchAddress("id", &id);  // particle id
    fTree->SetBranchAddress("parentid",
                            &parentid);  // parent id, could be different
    fTree->SetBranchAddress("pythiaid",
                            &pythiaid);      // pythiaid original particle
    fTree->SetBranchAddress("ecut", &ecut);  // energy cut used in simulation
    fTree->SetBranchAddress("w", &w);        // weight of event
    //  check if ntuple has information of momentum at origin
    if (fTree->GetListOfLeaves()->GetSize() < 17) {
      fTree->SetBranchAddress(
          "x", &vx);  // position with respect to startOfTarget at -89.27m
      fTree->SetBranchAddress("y", &vy);
      fTree->SetBranchAddress("z", &vz);
      fTree->SetBranchAddress("px", &px);  // momentum
      fTree->SetBranchAddress("py", &py);
      fTree->SetBranchAddress("pz", &pz);
    } else {
      fTree->SetBranchAddress(
          "ox", &vx);  // position with respect to startOfTarget at -50m
      fTree->SetBranchAddress("oy", &vy);
      fTree->SetBranchAddress("oz", &vz);
      fTree->SetBranchAddress("opx", &px);  // momentum
      fTree->SetBranchAddress("opy", &py);
      fTree->SetBranchAddress("opz", &pz);
    }
  } else {
    id = -1;
    fTree = new TChain("cbmsim");
    for (auto& f : fileNames) {
      LOG(info) << "Opening input file " << f;
      fTree->Add(f.c_str());
    }
    fNevents = fTree->GetEntries();
    LOG(info) << "Reading " << fNevents << " entries";
    // Detect format by checking branch name:
    // STL format uses PlaneHAPoint, TClonesArray uses vetoPoint
    TBranch* mcBranch = fTree->GetBranch("MCTrack");
    if (!mcBranch) {
      LOG(fatal) << "MCTrack branch not found in input file";
    }
 
    if (fTree->GetBranch("PlaneHAPoint")) {
      // New STL format
      fUseSTL = true;
      MCTrack_vec = nullptr;
      vetoPoints_vec = nullptr;
      auto mcStatus = fTree->SetBranchAddress("MCTrack", &MCTrack_vec);
      auto vetoStatus =
          fTree->SetBranchAddress("PlaneHAPoint", &vetoPoints_vec);
      if (mcStatus < 0 || vetoStatus < 0) {
        LOG(fatal) << "Failed to set branch addresses for STL vector format";
      }
      LOG(info) << "Using STL vector format (PlaneHAPoint)";
    } else if (fTree->GetBranch("vetoPoint")) {
      // Old TClonesArray format
      fUseSTL = false;
      MCTrack = new TClonesArray("ShipMCTrack");
      vetoPoints = new TClonesArray("vetoPoint");
      auto mcStatus = fTree->SetBranchAddress("MCTrack", &MCTrack);
      auto vetoStatus = fTree->SetBranchAddress("vetoPoint", &vetoPoints);
      if (mcStatus < 0 || vetoStatus < 0) {
        LOG(fatal) << "Failed to set branch addresses for TClonesArray format";
      }
      LOG(info) << "Using TClonesArray format (vetoPoint)";
    } else {
      LOG(fatal)
          << "Neither PlaneHAPoint nor vetoPoint branch found in input file";
    }
  }
  return kTRUE;
}

◆ ReadEvent()

Bool_t MuonBackGenerator::ReadEvent ( FairPrimaryGenerator * cpg )

override

public method ReadEvent

Definition at line 177 of file MuonBackGenerator.cxx.

                                                             {
  auto* pdgBase = TDatabasePDG::Instance();
  Double_t mass = 0., e = 0., tof = 0.;
  std::unordered_map<int, int> muList;
  std::unordered_map<int, std::vector<int>> moList;
 
  // Helper lambdas to abstract STL vs TClonesArray access
  auto getVetoSize = [this]() -> int {
    return fUseSTL ? vetoPoints_vec->size() : vetoPoints->GetEntries();
  };
  auto getVetoPoint = [this](int i) -> vetoPoint* {
    return fUseSTL ? &(*vetoPoints_vec)[i]
                   : dynamic_cast<vetoPoint*>(vetoPoints->At(i));
  };
  auto getMCTrack = [this](int i) -> ShipMCTrack* {
    return fUseSTL ? &(*MCTrack_vec)[i]
                   : dynamic_cast<ShipMCTrack*>(MCTrack->At(i));
  };
  auto getMCTrackSize = [this]() -> int {
    return fUseSTL ? MCTrack_vec->size() : MCTrack->GetEntries();
  };
 
  bool found_muon = false;
  while (fn < fNevents) {
    fTree->GetEntry(fn);
    muList.clear();
    moList.clear();
    fn++;
    if (fn % 100000 == 0) {
      LOGF(info, "Reading event %i", fn);
    }
    // test if we have a muon, don't look at neutrinos:
    if (TMath::Abs(static_cast<int>(id)) == 13) {
      mass = pdgBase->GetParticle(id)->Mass();
      e = TMath::Sqrt(px * px + py * py + pz * pz + mass * mass);
      tof = 0;
      found_muon = true;
      break;
    }
    if (id == -1) {  // use tree as input file
      Bool_t found = false;
      for (int i = 0; i < getVetoSize(); i++) {
        auto* v = getVetoPoint(i);
        Int_t abspid = TMath::Abs(v->PdgCode());
        if (abspid == 13 || (!followMuons && abspid != 12 && abspid != 14)) {
          found = true;
          Int_t muIndex = v->GetTrackID();
          if (!fdownScaleDiMuon) {
            muList.insert({muIndex, i});
          } else if (abspid == 13) {
            if (checkDiMuon(muIndex)) {
              moList[getMCTrack(muIndex)->GetMotherId()].push_back(i);
            } else {
              muList.insert({muIndex, i});
            }
          }
        }
      }
      // reject muon if comes from boosted channel
 
      for (auto it = moList.begin(); it != moList.end(); it++) {
        if (gRandom->Uniform(0., 1.) > 0.99) {
          std::vector<int> list = it->second;
          for (unsigned i = 0; i < list.size(); i++) {
            auto* v = getVetoPoint(list.at(i));
            Int_t muIndex = v->GetTrackID();
            muList.insert({muIndex, i});
          }
        }
      }
      if (!found) {
        LOGF(debug, "No muon found %i", fn - 1);
      }
      if (found) {
        found_muon = true;
        break;
      }
    }
  }
  if (fn >= fNevents && !found_muon) {
    LOGF(info, "End of tree reached %i", fNevents);
    gMC->StopRun();
    return kTRUE;
  }
  if (fSameSeed) {
    Int_t theSeed = fn + fSameSeed * fNevents;
    LOGF(debug, "Seed: %d", theSeed);
    gRandom->SetSeed(theSeed);
  }
  auto [dx, dy] = CalculateBeamOffset(fsmearBeam, fPaintBeam);
  if (id == -1) {
    for (int i = 0; i < getMCTrackSize(); i++) {
      auto* track = getMCTrack(i);
      Int_t abspid = TMath::Abs(track->GetPdgCode());
      px = track->GetPx();
      py = track->GetPy();
      pz = track->GetPz();
      if (fPhiRandomize) {
        double phi_random = gRandom->Uniform(0., 2 * TMath::Pi());
        Double_t pt = track->GetPt();
        px = pt * TMath::Cos(phi_random);
        py = pt * TMath::Sin(phi_random);
      }
      vx = track->GetStartX() + dx;
      vy = track->GetStartY() + dy;
      vz = track->GetStartZ();
      tof = track->GetStartT() / 1E9;  // convert back from ns to sec;
      e = track->GetEnergy();
      Bool_t wanttracking = false;  // only transport muons
      for (std::pair<int, int> element : muList) {
        if (element.first == i) {
          wanttracking = true;
          if (!followMuons) {
            auto* v = getVetoPoint(element.second);
            TVector3 lpv = v->LastPoint();
            TVector3 lmv = v->LastMom();
            if (abspid == 22) {
              e = lmv.Mag();
            } else {
              e = TMath::Sqrt(lmv.Mag2() +
                              (track->GetMass()) * (track->GetMass()));
            }
            px = lmv[0];
            py = lmv[1];
            pz = lmv[2];
            vx = lpv[0];
            vy = lpv[1];
            vz = lpv[2];
            tof = v->GetTime() / 1E9;  // convert back from ns to sec
          }
          break;
        }
      }
      cpg->AddTrack(track->GetPdgCode(), px, py, pz, vx, vy, vz,
                    track->GetMotherId(), wanttracking, e, tof,
                    track->GetWeight(), (TMCProcess)track->GetProcID());
    }
  } else {
    vx += dx / 100.;
    vy += dy / 100.;
    if (fPhiRandomize) {
      double phi_random = gRandom->Uniform(0., 2 * TMath::Pi());
      Double_t pt = TMath::Sqrt(px * px + py * py);
      px = pt * TMath::Cos(phi_random);
      py = pt * TMath::Sin(phi_random);
    }
    cpg->AddTrack(static_cast<int>(pythiaid), px, py, pz, vx * 100., vy * 100.,
                  vz * 100., -1., false, e, pythiaid, parentid);
    cpg->AddTrack(static_cast<int>(id), px, py, pz, vx * 100., vy * 100.,
                  vz * 100., -1., true, e, tof, w);
  }
  return kTRUE;
}

◆ SetDownScaleDiMuon()

void MuonBackGenerator::SetDownScaleDiMuon ( )

inline

Definition at line 46 of file MuonBackGenerator.h.

46{ fdownScaleDiMuon = kTRUE; };

◆ SetPaintRadius()

void MuonBackGenerator::SetPaintRadius ( Double_t r )

inline

Definition at line 39 of file MuonBackGenerator.h.

39{ fPaintBeam = r; };

◆ SetPhiRandomize()

void MuonBackGenerator::SetPhiRandomize ( Bool_t phiRandomize )

inline

Definition at line 44 of file MuonBackGenerator.h.

44{ fPhiRandomize = phiRandomize; };

◆ SetSameSeed()

void MuonBackGenerator::SetSameSeed ( Int_t s )

inline

Definition at line 40 of file MuonBackGenerator.h.

                            {
    LOGF(info, "Seed: %d", s);
    fSameSeed = s;
  };

◆ SetSmearBeam()

void MuonBackGenerator::SetSmearBeam ( Double_t sb )

inline

Definition at line 38 of file MuonBackGenerator.h.

38{ fsmearBeam = sb; };

Member Data Documentation

◆ ecut

Float_t MuonBackGenerator::ecut

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ f_zOffset

float MuonBackGenerator::f_zOffset

protected

Definition at line 59 of file MuonBackGenerator.h.

◆ fdownScaleDiMuon

Bool_t MuonBackGenerator::fdownScaleDiMuon

protected

Definition at line 63 of file MuonBackGenerator.h.

◆ fInputFile

TFile* MuonBackGenerator::fInputFile

protected

flag to indicate if using STL vectors

Definition at line 56 of file MuonBackGenerator.h.

◆ fn

int MuonBackGenerator::fn

protected

Definition at line 60 of file MuonBackGenerator.h.

◆ fNevents

int MuonBackGenerator::fNevents

protected

Definition at line 58 of file MuonBackGenerator.h.

◆ followMuons

Bool_t MuonBackGenerator::followMuons

protected

Definition at line 64 of file MuonBackGenerator.h.

◆ fPaintBeam

Double_t MuonBackGenerator::fPaintBeam

protected

Definition at line 62 of file MuonBackGenerator.h.

◆ fPhiRandomize

Bool_t MuonBackGenerator::fPhiRandomize

protected

Definition at line 61 of file MuonBackGenerator.h.

◆ fSameSeed

Int_t MuonBackGenerator::fSameSeed

protected

Definition at line 65 of file MuonBackGenerator.h.

◆ fsmearBeam

Double_t MuonBackGenerator::fsmearBeam

protected

Definition at line 66 of file MuonBackGenerator.h.

◆ fTree

TChain* MuonBackGenerator::fTree

protected

Definition at line 57 of file MuonBackGenerator.h.

◆ fUseSTL

Bool_t MuonBackGenerator::fUseSTL

protected

Definition at line 55 of file MuonBackGenerator.h.

◆ id

Float_t MuonBackGenerator::id

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ MCTrack

TClonesArray* MuonBackGenerator::MCTrack

protected

Definition at line 51 of file MuonBackGenerator.h.

◆ MCTrack_vec

std::vector<ShipMCTrack>* MuonBackGenerator::MCTrack_vec

protected

Definition at line 53 of file MuonBackGenerator.h.

◆ parentid

Float_t MuonBackGenerator::parentid

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ px

Float_t MuonBackGenerator::px

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ py

Float_t MuonBackGenerator::py

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ pythiaid

Float_t MuonBackGenerator::pythiaid

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ pz

Float_t MuonBackGenerator::pz

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ vetoPoints

TClonesArray* MuonBackGenerator::vetoPoints

protected

Definition at line 52 of file MuonBackGenerator.h.

◆ vetoPoints_vec

std::vector<vetoPoint>* MuonBackGenerator::vetoPoints_vec

protected

Definition at line 54 of file MuonBackGenerator.h.

◆ vx

Float_t MuonBackGenerator::vx

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ vy

Float_t MuonBackGenerator::vy

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ vz

Float_t MuonBackGenerator::vz

protected

Definition at line 50 of file MuonBackGenerator.h.

◆ w

Float_t MuonBackGenerator::w

protected

Definition at line 50 of file MuonBackGenerator.h.

The documentation for this class was generated from the following files:

shipgen/MuonBackGenerator.h
shipgen/MuonBackGenerator.cxx

Public Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ MuonBackGenerator()

◆ ~MuonBackGenerator()

Member Function Documentation

◆ checkDiMuon()

◆ CloseFile()

◆ FollowAllParticles()

◆ GetNevents()

◆ Init() [1/4]

◆ Init() [2/4]

◆ Init() [3/4]

◆ Init() [4/4]

◆ ReadEvent()

◆ SetDownScaleDiMuon()

◆ SetPaintRadius()

◆ SetPhiRandomize()

◆ SetSameSeed()

◆ SetSmearBeam()

Member Data Documentation

◆ ecut

◆ f_zOffset

◆ fdownScaleDiMuon

◆ fInputFile

◆ fn

◆ fNevents

◆ followMuons

◆ fPaintBeam

◆ fPhiRandomize

◆ fSameSeed

◆ fsmearBeam

◆ fTree

◆ fUseSTL

◆ id

◆ MCTrack

◆ MCTrack_vec

◆ parentid

◆ px

◆ py

◆ pythiaid

◆ pz

◆ vetoPoints

◆ vetoPoints_vec

◆ vx

◆ vy

◆ vz

◆ w