HNLPythia8Generator Class Reference

#include <HNLPythia8Generator.h>

Inheritance diagram for HNLPythia8Generator:

Collaboration diagram for HNLPythia8Generator:

Public Member Functions
	HNLPythia8Generator ()
	~HNLPythia8Generator () override
Bool_t	ReadEvent (FairPrimaryGenerator *) override
void	SetParameters (char *)
void	Print ()
void	List (int id)
Bool_t	Init (const char *inFile) override
Bool_t	Init (const char *inFile, int startEvent) override
Bool_t	Init () override
void	SetMom (Double_t mom)
void	SetId (Double_t id)
void	SetHNLId (Int_t id)
void	SetLmin (Double_t z)
void	SetLmax (Double_t z)
void	SetSmearBeam (Double_t sb)
void	SetPaintRadius (Double_t r)
void	SetfFDs (Double_t z)
void	UseRandom1 ()
void	UseRandom3 ()
void	UseDeepCopy ()
Int_t	nrOfRetries ()
Pythia8::Pythia *	getPythiaInstance ()
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)

Public Attributes
Pythia8::Pythia *	fPythia

Protected Attributes
Double_t	fMom
Int_t	fHNL
Int_t	fId
Bool_t	fUseRandom1
Bool_t	fUseRandom3
Double_t	fLmin
Double_t	fLmax
Int_t	fnRetries
Double_t	fctau
Double_t	fFDs
Double_t	fsmearBeam
Double_t	fPaintBeam
TFile *	fInputFile
TTree *	fTree
	pointer to a file
Int_t	fNevents
Int_t	fn
Int_t	fShipEventNr
Float_t	hpx [1]
Float_t	hpy [1]
Float_t	hpz [1]
Float_t	hE [1]
Float_t	hM [1]
Float_t	mpx [1]
Float_t	mpy [1]
Float_t	mpz [1]
Float_t	mE [1]
Float_t	hid [1]
Float_t	mid [1]
Bool_t	fDeepCopy
FairLogger *	fLogger
Protected Attributes inherited from SHiP::Generator
std::optional< std::string >	fextFile
Int_t	firstEvent = 0

Private Attributes
std::shared_ptr< Pythia8::RndmEngine >	fRandomEngine

Detailed Description

Definition at line 42 of file HNLPythia8Generator.h.

Constructor & Destructor Documentation

HNLPythia8Generator()

HNLPythia8Generator::HNLPythia8Generator

(

)

default constructor

Definition at line 23 of file HNLPythia8Generator.cxx.

                                         {
  fId = 2212;           // proton
  fMom = 400;           // proton
  fHNL = 9900015;       // HNL  pdg code
  fLmin = 5000. * cm;   // mm minimum  decay position z  ROOT units !
  fLmax = 12000. * cm;  // mm maximum decay position z
  fFDs = 7.7 / 10.4;  // correction for Pythia6 to match measured Ds production
  fsmearBeam = 8 * mm;  // default value for smearing beam (8 mm)
  fPaintBeam = 5 * cm;  // default value for painting beam (5 cm)
  fInputFile = nullptr;
  fnRetries = 0;
  fShipEventNr = 0;
  fPythia = new Pythia8::Pythia();
}

~HNLPythia8Generator()

HNLPythia8Generator::~HNLPythia8Generator ( )

override

destructor

Definition at line 94 of file HNLPythia8Generator.cxx.

{}

Member Function Documentation

getPythiaInstance()

Pythia8::Pythia * HNLPythia8Generator::getPythiaInstance ( )

inline

Definition at line 90 of file HNLPythia8Generator.h.

{ return fPythia; };

Init() [1/3]

Bool_t HNLPythia8Generator::Init ( )

override

Definition at line 40 of file HNLPythia8Generator.cxx.

                                 {
  if (debug) {
    List(9900015);
  }
  if (fUseRandom1) fRandomEngine = std::make_shared<PyTr1Rng>();
  if (fUseRandom3) fRandomEngine = std::make_shared<PyTr3Rng>();
  fPythia->setRndmEnginePtr(fRandomEngine);
  fn = 0;
  if (fextFile) {
    fInputFile = TFile::Open(fextFile->c_str());
    LOG(info) << "Open external file with charm or beauty hadrons: "
              << *fextFile;
    if (!fInputFile) {
      LOG(fatal) << "Error opening input file.";
      return kFALSE;
    }
 
    fTree = fInputFile->Get<TTree>("pythia6");
    fNevents = fTree->GetEntries();
    fn = firstEvent;
    fTree->SetBranchAddress("id", &hid);  // particle id
    fTree->SetBranchAddress("px", &hpx);  // momentum
    fTree->SetBranchAddress("py", &hpy);
    fTree->SetBranchAddress("pz", &hpz);
    fTree->SetBranchAddress("E", &hE);
    fTree->SetBranchAddress("M", &hM);
    fTree->SetBranchAddress("mid", &mid);  // mother
    fTree->SetBranchAddress("mpx", &mpx);  // momentum
    fTree->SetBranchAddress("mpy", &mpy);
    fTree->SetBranchAddress("mpz", &mpz);
    fTree->SetBranchAddress("mE", &mE);
  } else {
    LOG(debug) << "Beam Momentum " << fMom;
    fPythia->settings.mode("Beams:idA", fId);
    fPythia->settings.mode("Beams:idB", 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
  }
  TDatabasePDG* pdgBase = TDatabasePDG::Instance();
  Double_t root_ctau = pdgBase->GetParticle(fHNL)->Lifetime();
  fctau = fPythia->particleData.tau0(fHNL);  //* 3.3333e-12
  LOG(debug) << "tau root " << root_ctau
             << "[s] ctau root = " << root_ctau * 3e10 << "[cm]";
  LOG(debug) << "ctau pythia " << fctau << "[mm]";
  if (debug) {
    List(9900015);
  }
  fPythia->init();
  return kTRUE;
}

Init() [2/3]

Bool_t HNLPythia8Generator::Init ( const char *  inFile )

inlineoverridevirtual

Implements SHiP::Generator.

Definition at line 58 of file HNLPythia8Generator.h.

{ return Init(inFile, 0); };

Init() [3/3]

Bool_t HNLPythia8Generator::Init ( const char *  inFile, int  startEvent  )

inlineoverridevirtual

Implements SHiP::Generator.

Definition at line 60 of file HNLPythia8Generator.h.

                                                           {
    LOG(warning) << "Init with files not implemented for HNLPythia8Generator. "
                    "Using default Init() instead";
    return Init();
  };

List()

void HNLPythia8Generator::List ( int  id )

inline

Definition at line 54 of file HNLPythia8Generator.h.

{ fPythia->particleData.list(id); };

nrOfRetries()

Int_t HNLPythia8Generator::nrOfRetries ( )

inline

Definition at line 89 of file HNLPythia8Generator.h.

{ return fnRetries; };

Print()

void HNLPythia8Generator::Print ( )

inline

Definition at line 53 of file HNLPythia8Generator.h.

{ fPythia->settings.listAll(); };

ReadEvent()

Bool_t HNLPythia8Generator::ReadEvent ( FairPrimaryGenerator *  cpg )

override

public method ReadEvent

Definition at line 98 of file HNLPythia8Generator.cxx.

                                                               {
  Double_t tp, tS, zp, xp, yp, zS, xS, yS, pz, px, py, e, w;
  Double_t tm, zm, xm, ym, pmz, pmx, pmy, em;
  Int_t im;
  // take HNL decay of Pythia, move it to the SHiP decay region
  // recalculate decay time
  // weight event with exp(-t_ship/tau_HNL) / exp(-t_pythia/tau_HNL)
 
  int iHNL = 0;  // index of the chosen HNL (the 1st one), also ensures that at
                 // least 1 HNL is produced
  std::vector<int>
      dec_chain;  // pythia indices of the particles to be stored on the stack
  std::vector<int> hnls;  // pythia indices of HNL particles
  do {
    if (fextFile) {
      // take charm or beauty hadron from external file
      // correct for too much Ds produced by pythia6
      bool x = true;
      while (x) {
        if (fn == fNevents) {
          LOG(warning) << "End of input file. Rewind.";
        }
        fTree->GetEntry(fn % fNevents);
        fn++;
        if (static_cast<int>(fabs(hid[0])) != 431) {
          x = false;
        } else {
          Double_t rnr = gRandom->Uniform(0, 1);
          if (rnr < fFDs) {
            x = false;
          };
          // cout<<"what is x "<<x<<" id "<<int(fabs(hid[0]))<<" rnr " << rnr
          // <<" "<< fFDs <<std::endl ;
        }
      }
      fPythia->event.reset();
      fPythia->event.append((Int_t)hid[0], 1, 0, 0, hpx[0], hpy[0], hpz[0],
                            hE[0], hM[0], 0., 9.);
    }
    fPythia->next();
    for (int i = 0; i < fPythia->event.size(); i++) {
      // find first HNL
      if (abs(fPythia->event[i].id()) == fHNL) {
        hnls.push_back(i);
      }
    }
    iHNL = hnls.size();
    if (iHNL == 0) {
      // fLogger->Info(MESSAGE_ORIGIN,"Event without HNL. Retry.");
      // fPythia->event.list();
      fnRetries += 1;  // can happen if phasespace does not allow charm hadron
                       // to decay to HNL
    } else {
      int r = static_cast<int>(gRandom->Uniform(0, iHNL));
      // cout << " ----> debug 2 " << r  <<  endl;
      int i = hnls[r];
      // production vertex
      zp = fPythia->event[i].zProd();
      xp = fPythia->event[i].xProd();
      yp = fPythia->event[i].yProd();
      tp = fPythia->event[i].tProd();
      // momentum
      pz = fPythia->event[i].pz();
      px = fPythia->event[i].px();
      py = fPythia->event[i].py();
      e = fPythia->event[i].e();
      // new decay vertex
      Double_t LS = gRandom->Uniform(fLmin, fLmax);  // mm, G4 and Pythia8 units
      Double_t p = TMath::Sqrt(px * px + py * py + pz * pz);
      Double_t lam = LS / p;
      xS = xp + lam * px;
      yS = yp + lam * py;
      zS = zp + lam * pz;
      Double_t gam = e / TMath::Sqrt(e * e - p * p);
      Double_t beta = p / e;
      tS = tp + LS / beta;  // units ? [mm/c] + [mm/beta] (beta is dimensionless
                            // speed, and c=1 here) if one would use [s], then
                            // tS = tp/(cm*c_light) + (LS/cm)/(beta*c_light) =
                            // tS/(cm*c_light) i.e. units look consistent
      w = TMath::Exp(-LS / (beta * gam * fctau)) *
          ((fLmax - fLmin) / (beta * gam * fctau));
      im = (Int_t)fPythia->event[i].mother1();
      zm = fPythia->event[im].zProd();
      xm = fPythia->event[im].xProd();
      ym = fPythia->event[im].yProd();
      pmz = fPythia->event[im].pz();
      pmx = fPythia->event[im].px();
      pmy = fPythia->event[im].py();
      em = fPythia->event[im].e();
      tm = fPythia->event[im].tProd();
      // foresee finite beam size
      auto [dx, dy] = CalculateBeamOffset(fsmearBeam, fPaintBeam);
      if (fextFile) {
        // take grand mother particle from input file, to know if primary or
        // secondary production
        cpg->AddTrack((Int_t)mid[0], mpx[0], mpy[0], mpz[0], xm / cm + dx / cm,
                      ym / cm + dy / cm, zm / cm, -1, false, mE[0], 0., 1.);
        cpg->AddTrack(
            (Int_t)fPythia->event[im].id(), pmx, pmy, pmz, xm / cm + dx / cm,
            ym / cm + dy / cm, zm / cm, 0, false, em, tm / cm / c_light,
            w);  // convert pythia's (x,y,z[mm], t[mm/c]) to ([cm], [s])
        cpg->AddTrack(fHNL, px, py, pz, xp / cm + dx / cm, yp / cm + dy / cm,
                      zp / cm, 1, false, e, tp / cm / c_light, w);
      } else {
        cpg->AddTrack(
            (Int_t)fPythia->event[im].id(), pmx, pmy, pmz, xm / cm + dx / cm,
            ym / cm + dy / cm, zm / cm, -1, false, em, tm / cm / c_light,
            w);  // convert pythia's (x,y,z[mm], t[mm/c]) to ([cm], [s])
        cpg->AddTrack(fHNL, px, py, pz, xp / cm + dx / cm, yp / cm + dy / cm,
                      zp / cm, 0, false, e, tp / cm / c_light, w);
      }
      // bookkeep the indices of stored particles
      dec_chain.push_back(im);
      dec_chain.push_back(i);
      // cout << endl << " insert mother pdg=" <<fPythia->event[im].id() << "
      // pmz = " << pmz << " [GeV],  zm = " << zm << " [mm] tm = " << tm << "
      // [mm/c]" << endl; cout << " ----> insert HNL =" << fHNL << " pz = " <<
      // pz << " [GeV] zp = " << zp << " [mm] tp = " << tp << " [mm/c]" << endl;
      iHNL = i;
    }
  } while (iHNL ==
           0);  // ----------- avoid rare empty events w/o any HNL's produced
 
  if (fShipEventNr % 100 == 0) {
    LOGF(info, "ship event %i / pythia event-nr %i", fShipEventNr, fn);
  }
  fShipEventNr += 1;
  // fill a container with pythia indices of the HNL decay chain
  for (int k = 0; k < fPythia->event.size(); k++) {
    // if daughter of HNL, copy
    im = fPythia->event[k].mother1();
    while (im > 0) {
      if (im == iHNL) {
        break;
      }  // pick the decay products of only 1 chosen HNL
      // if ( abs(fPythia->event[im].id())==fHNL && im == iHNL ){break;}
      else {
        im = fPythia->event[im].mother1();
      }
    }
    if (im < 1) {
      continue;
    }
    dec_chain.push_back(k);
  }
 
  // go over daughters and store them on the stack, starting from 2 to account
  // for HNL and its mother
  for (std::vector<int>::iterator it = dec_chain.begin() + 2;
       it != dec_chain.end(); ++it) {
    // pythia index of the particle to store
    int k = *it;
    // find mother position on the output stack: impy -> im
    int impy = fPythia->event[k].mother1();
    std::vector<int>::iterator itm =
        std::find(dec_chain.begin(), dec_chain.end(), impy);
    im = -1;  // safety
    if (itm != dec_chain.end())
      im = itm - dec_chain.begin();  // convert iterator into sequence number
 
    Bool_t wanttracking = false;
    if (fPythia->event[k].isFinal()) {
      wanttracking = true;
    }
    pz = fPythia->event[k].pz();
    px = fPythia->event[k].px();
    py = fPythia->event[k].py();
    e = fPythia->event[k].e();
    if (fextFile) {
      im += 1;
    }
    cpg->AddTrack((Int_t)fPythia->event[k].id(), px, py, pz, xS / cm, yS / cm,
                  zS / cm, im, wanttracking, e, tS / cm / c_light, w);
    // std::cout <<k<< " insert pdg =" <<fPythia->event[k].id() << " pz = " <<
    // pz << " [GeV] zS = " << zS << " [mm] tS = " << tS << "[mm/c]" <<  endl;
  }
  return kTRUE;
}

SetfFDs()

void HNLPythia8Generator::SetfFDs ( Double_t  z )

inline

Definition at line 78 of file HNLPythia8Generator.h.

{ fFDs = z; };

SetHNLId()

void HNLPythia8Generator::SetHNLId ( Int_t  id )

inline

Definition at line 69 of file HNLPythia8Generator.h.

{ fHNL = id; };

SetId()

void HNLPythia8Generator::SetId ( Double_t  id )

inline

Definition at line 68 of file HNLPythia8Generator.h.

{ fId = id; };

SetLmax()

void HNLPythia8Generator::SetLmax ( Double_t  z )

inline

Definition at line 71 of file HNLPythia8Generator.h.

{ fLmax = z * 10; };  // Conversion from cm to mm

SetLmin()

void HNLPythia8Generator::SetLmin ( Double_t  z )

inline

Definition at line 70 of file HNLPythia8Generator.h.

{ fLmin = z * 10; };  // Conversion from cm to mm

SetMom()

void HNLPythia8Generator::SetMom ( Double_t  mom )

inline

Definition at line 67 of file HNLPythia8Generator.h.

{ fMom = mom; };

SetPaintRadius()

void HNLPythia8Generator::SetPaintRadius ( Double_t  r )

inline

Definition at line 75 of file HNLPythia8Generator.h.

                                  {
    fPaintBeam = r * 10;
  };  // Conversion from cm to mm

SetParameters()

void HNLPythia8Generator::SetParameters

(

char *

par

)

Definition at line 277 of file HNLPythia8Generator.cxx.

                                                 {
  // Set Parameters
  fPythia->readString(par);
  LOG(debug) << "fPythia->readString(\"" << par << "\")";
}

SetSmearBeam()

void HNLPythia8Generator::SetSmearBeam ( Double_t  sb )

inline

Definition at line 72 of file HNLPythia8Generator.h.

                                 {
    fsmearBeam = sb * 10;
  };  // Conversion from cm to mm

UseDeepCopy()

void HNLPythia8Generator::UseDeepCopy ( )

inline

Definition at line 88 of file HNLPythia8Generator.h.

{ fDeepCopy = kTRUE; };

UseRandom1()

void HNLPythia8Generator::UseRandom1 ( )

inline

Definition at line 79 of file HNLPythia8Generator.h.

                    {
    fUseRandom1 = kTRUE;
    fUseRandom3 = kFALSE;
  };

UseRandom3()

void HNLPythia8Generator::UseRandom3 ( )

inline

Definition at line 83 of file HNLPythia8Generator.h.

                    {
    fUseRandom1 = kFALSE;
    fUseRandom3 = kTRUE;
  };

Member Data Documentation

fctau

Double_t HNLPythia8Generator::fctau

protected

Definition at line 104 of file HNLPythia8Generator.h.

fDeepCopy

Bool_t HNLPythia8Generator::fDeepCopy

protected

Definition at line 113 of file HNLPythia8Generator.h.

fFDs

Double_t HNLPythia8Generator::fFDs

protected

Definition at line 105 of file HNLPythia8Generator.h.

fHNL

Int_t HNLPythia8Generator::fHNL

protected

Definition at line 97 of file HNLPythia8Generator.h.

fId

Int_t HNLPythia8Generator::fId

protected

Definition at line 98 of file HNLPythia8Generator.h.

fInputFile

TFile* HNLPythia8Generator::fInputFile

protected

Definition at line 108 of file HNLPythia8Generator.h.

fLmax

Double_t HNLPythia8Generator::fLmax

protected

Definition at line 102 of file HNLPythia8Generator.h.

fLmin

Double_t HNLPythia8Generator::fLmin

protected

Definition at line 101 of file HNLPythia8Generator.h.

fLogger

FairLogger* HNLPythia8Generator::fLogger

protected

Definition at line 114 of file HNLPythia8Generator.h.

fMom

Double_t HNLPythia8Generator::fMom

protected

Definition at line 96 of file HNLPythia8Generator.h.

fn

Int_t HNLPythia8Generator::fn

protected

Definition at line 110 of file HNLPythia8Generator.h.

fNevents

Int_t HNLPythia8Generator::fNevents

protected

Definition at line 110 of file HNLPythia8Generator.h.

fnRetries

Int_t HNLPythia8Generator::fnRetries

protected

Definition at line 103 of file HNLPythia8Generator.h.

fPaintBeam

Double_t HNLPythia8Generator::fPaintBeam

protected

Definition at line 107 of file HNLPythia8Generator.h.

fPythia

Pythia8::Pythia* HNLPythia8Generator::fPythia

Definition at line 91 of file HNLPythia8Generator.h.

fRandomEngine

std::shared_ptr<Pythia8::RndmEngine> HNLPythia8Generator::fRandomEngine

private

Definition at line 93 of file HNLPythia8Generator.h.

fShipEventNr

Int_t HNLPythia8Generator::fShipEventNr

protected

Definition at line 110 of file HNLPythia8Generator.h.

fsmearBeam

Double_t HNLPythia8Generator::fsmearBeam

protected

Definition at line 106 of file HNLPythia8Generator.h.

fTree

TTree* HNLPythia8Generator::fTree

protected

pointer to a file

Definition at line 109 of file HNLPythia8Generator.h.

fUseRandom1

Bool_t HNLPythia8Generator::fUseRandom1

protected

Definition at line 99 of file HNLPythia8Generator.h.

fUseRandom3

Bool_t HNLPythia8Generator::fUseRandom3

protected

Definition at line 100 of file HNLPythia8Generator.h.

hE

Float_t HNLPythia8Generator::hE[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

hid

Float_t HNLPythia8Generator::hid[1]

protected

Definition at line 112 of file HNLPythia8Generator.h.

hM

Float_t HNLPythia8Generator::hM[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

hpx

Float_t HNLPythia8Generator::hpx[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

hpy

Float_t HNLPythia8Generator::hpy[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

hpz

Float_t HNLPythia8Generator::hpz[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

mE

Float_t HNLPythia8Generator::mE[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

mid

Float_t HNLPythia8Generator::mid[1]

protected

Definition at line 112 of file HNLPythia8Generator.h.

mpx

Float_t HNLPythia8Generator::mpx[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

mpy

Float_t HNLPythia8Generator::mpy[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

mpz

Float_t HNLPythia8Generator::mpz[1]

protected

Definition at line 111 of file HNLPythia8Generator.h.

---

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

• shipgen/HNLPythia8Generator.h
• shipgen/HNLPythia8Generator.cxx