MuDISGenerator Class Reference

#include <MuDISGenerator.h>

Inheritance diagram for MuDISGenerator:

Collaboration diagram for MuDISGenerator:

Public Member Functions
	MuDISGenerator ()
	~MuDISGenerator () override
Bool_t	ReadEvent (FairPrimaryGenerator *) override
Bool_t	Init (const char *, int) override
Bool_t	Init (const char *) override
Int_t	GetNevents ()
void	SetPositions (Double_t z_start, Double_t z_end)
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
Double_t	startZ
Double_t	endZ
TClonesArray *	iMuon
TClonesArray *	dPart
TClonesArray *	dPartSoft
FairLogger *	fLogger
TFile *	fInputFile
	don't make it persistent, magic ROOT command
TTree *	fTree
int	fNevents
int	fn
bool	fFirst
Protected Attributes inherited from SHiP::Generator
std::optional< std::string >	fextFile
Int_t	firstEvent = 0

Detailed Description

Definition at line 19 of file MuDISGenerator.h.

Constructor & Destructor Documentation

MuDISGenerator()

MuDISGenerator::MuDISGenerator

(

)

default constructor

Definition at line 31 of file MuDISGenerator.cxx.

{}

~MuDISGenerator()

MuDISGenerator::~MuDISGenerator ( )

override

destructor

Definition at line 59 of file MuDISGenerator.cxx.

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

Member Function Documentation

GetNevents()

Int_t MuDISGenerator::GetNevents

(

)

Definition at line 227 of file MuDISGenerator.cxx.

{ return fNevents; }

Init() [1/2]

Bool_t MuDISGenerator::Init ( const char *  fileName )

overridevirtual

Implements SHiP::Generator.

Definition at line 34 of file MuDISGenerator.cxx.

{ return Init(fileName, 0); }

Init() [2/2]

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

overridevirtual

Implements SHiP::Generator.

Definition at line 35 of file MuDISGenerator.cxx.

                                                                      {
  LOGF(info, "Opening input file %s", fileName);
 
  iMuon = 0;
  dPart = 0;
  dPartSoft = 0;
  fInputFile = TFile::Open(fileName);
  if (!fInputFile) {
    LOG(fatal) << "Error opening input file";
    return kFALSE;
  }
  fTree = fInputFile->Get<TTree>("DIS");
  fNevents = fTree->GetEntries();
  fn = startEvent;
 
  fTree->SetBranchAddress("InMuon", &iMuon);  // incoming muon
  fTree->SetBranchAddress("DISParticles", &dPart);
  fTree->SetBranchAddress("SoftParticles",
                          &dPartSoft);  // Soft interaction particles
  LOG(info) << "MuDISGenerator: Initialization successful.";
  return kTRUE;
}

ReadEvent()

Bool_t MuDISGenerator::ReadEvent ( FairPrimaryGenerator *  cpg )

override

Definition at line 65 of file MuDISGenerator.cxx.

                                                          {
  if (fn == fNevents) {
    LOG(warning) << "End of input file. Rewind.";
  }
  fTree->GetEntry(fn % fNevents);
  fn++;
  if (fn % 10 == 0) {
    LOG(info) << "Info MuDISGenerator: MuDIS event-nr " << fn;
  }
 
  int nf = dPart->GetEntries();
  LOG(debug) << "*********************************************************";
  LOG(debug) << "muon DIS Generator debug " << iMuon->GetEntries() << " "
             << iMuon->AddrAt(0) << " nf " << nf << " fn=" << fn;
 
  // some start/end positions in z (f.i. emulsion to Tracker 1)
  Double_t start[3] = {0., 0., startZ};
  Double_t end[3] = {0., 0., endZ};
 
  // incoming muon  array('d',[pid,px,py,pz,E,x,y,z,w,t])
  TVectorD* mu = dynamic_cast<TVectorD*>(iMuon->AddrAt(0));
  LOG(debug) << "muon DIS Generator in muon " << int((*mu)[0]);
  Double_t x = (*mu)[5] * 100.;  // come in m -> cm
  Double_t y = (*mu)[6] * 100.;  // come in m -> cm
  Double_t z = (*mu)[7] * 100.;  // come in m -> cm
  Double_t w =
      (*mu)[8];  // weight of the original muon ( normalised to a spill)
  Double_t cross_sec = (*mu)[10];             // in mbarns
  Double_t t_muon = (*mu)[11];                // in ns
  Double_t DIS_multiplicity = 1 / (*mu)[12];  // 1/nDIS
 
  // calculate start/end positions along this muon, and amount of material in
  // between
 
  Double_t txmu = (*mu)[1] / (*mu)[3];
  Double_t tymu = (*mu)[2] / (*mu)[3];
  start[0] = x - (z - start[2]) * txmu;
  start[1] = y - (z - start[2]) * tymu;
  end[0] = x - (z - end[2]) * txmu;
  end[1] = y - (z - end[2]) * tymu;
  LOG(debug) << "MuDIS: mu xyz position " << x << ", " << y << ", " << z;
  LOG(debug) << "MuDIS: mu pxyz position " << (*mu)[1] << ", " << (*mu)[2]
             << ", " << (*mu)[3];
  LOG(debug) << "MuDIS: mu weight*DISmultiplicity  " << w;
  LOG(debug) << "MuDIS: mu DIS cross section " << cross_sec;
  LOG(debug) << "MuDIS: start position " << start[0] << ", " << start[1] << ", "
             << start[2];
  LOG(debug) << "MuDIS: end position " << end[0] << ", " << end[1] << ", "
             << end[2];
 
  Double_t bparam;
  Double_t mparam[10];
  bparam = shipgen::MeanMaterialBudget(start, end, mparam);
  // loop over trajectory between start and end to pick an interaction point
  Double_t prob2int = 0.;
  Double_t xmu = 0.;
  Double_t ymu = 0.;
  Double_t zmu = 0.;
  Int_t count = 0;
  LOG(debug) << "Info MuDISGenerator Start prob2int while loop, bparam= "
             << bparam << ", " << bparam * 1.e8;
  LOG(debug) << "Info MuDISGenerator What was maximum density, mparam[7]= "
             << mparam[7] << ", " << mparam[7] * 1.e8;
 
  while (prob2int < gRandom->Uniform(0., 1.)) {
    zmu = gRandom->Uniform(start[2], end[2]);
    xmu = x - (z - zmu) * txmu;
    ymu = y - (z - zmu) * tymu;
    // get local material at this point
    TGeoNode* node = gGeoManager->FindNode(xmu, ymu, zmu);
    TGeoMaterial* mat = nullptr;
    if (node && !gGeoManager->IsOutside())
      mat = node->GetVolume()->GetMaterial();
    LOG(debug) << "Info MuDISGenerator: mat " << count << ", " << mat->GetName()
               << ", " << mat->GetDensity();
    // density relative to Prob largest density along this trajectory, i.e. use
    // rho(Pt)
    prob2int = mat->GetDensity() / mparam[7];
    if (prob2int > 1.) {
      LOG(warning) << "MuDISGenerator: prob2int > Maximum density????";
      LOG(warning) << "prob2int: " << prob2int;
      LOG(warning) << "maxrho: " << mparam[7];
      LOG(warning) << "material: " << mat->GetName();
    }
    count += 1;
  }
 
  LOG(debug) << "Info MuDISGenerator: prob2int " << prob2int << ", " << count;
  LOG(debug) << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu;
 
  Double_t total_mom =
      TMath::Sqrt(TMath::Power((*mu)[1], 2) + TMath::Power((*mu)[2], 2) +
                  TMath::Power((*mu)[3], 2));  // in GeV
 
  Double_t distance =
      TMath::Sqrt(TMath::Power(x - xmu, 2) + TMath::Power(y - ymu, 2) +
                  TMath::Power(z - zmu, 2));  // in cm
 
  Double_t v =
      c_light * total_mom /
      TMath::Sqrt(TMath::Power(total_mom, 2) + TMath::Power(muon_mass, 2));
  Double_t t_rmu = distance / v;
 
  // Adjust time based on the relative positions
  if (zmu < z) {
    t_rmu = -t_rmu;
  }
 
  Double_t t_DIS =
      (t_muon + t_rmu) / 1e9;  // time taken in seconds to reach [xmu,ymu,zmu]
 
  cpg->AddTrack(static_cast<int>((*mu)[0]),  // incoming muon track ()
                (*mu)[1], (*mu)[2], (*mu)[3], xmu, ymu, zmu, -1,
                false,  // tracking disabled
                (*mu)[4],
                t_DIS,  // shift time of the incoming muon track wrt t_muon from
                        // the input file.
                w * DIS_multiplicity);  // muon weight associated with a spill*
                                        // DISmultiplicity
 
  // outgoing DIS particles, [did,dpx,dpy,dpz,E], put density along trajectory
  // as weight, g/cm^2
 
  w = mparam[0] * mparam[4];  // modify weight, by multiplying with average
                              // density * track length
  int index = 0;
  for (auto&& particle : *dPart) {
    TVectorD* Part = dynamic_cast<TVectorD*>(particle);
    LOG(debug) << "muon DIS Generator out part " << int((*Part)[0]);
    LOG(debug) << "muon DIS Generator out part mom " << (*Part)[1] << " "
               << (*Part)[2] << " " << (*Part)[3] << " " << (*Part)[4];
    LOG(debug) << "muon DIS Generator out part pos " << xmu << " " << ymu << ""
               << zmu;
    LOG(debug) << "muon DIS Generator out part w " << w;
 
    if (index == 0) {
      cpg->AddTrack(static_cast<int>((*Part)[0]), (*Part)[1], (*Part)[2],
                    (*Part)[3], xmu, ymu, zmu, 0, true, (*Part)[4], t_DIS,
                    cross_sec);  // save DIS cross section in MCTrack[1]
    } else {
      cpg->AddTrack(static_cast<int>((*Part)[0]), (*Part)[1], (*Part)[2],
                    (*Part)[3], xmu, ymu, zmu, 0, true, (*Part)[4], t_DIS, w);
    }
    index += 1;
  }
 
  // Soft interaction tracks
  for (auto&& softParticle : *dPartSoft) {
    TVectorD* SoftPart = dynamic_cast<TVectorD*>(softParticle);
    if ((*SoftPart)[7] > zmu) {
      continue;
    }  // Soft interactions after the DIS point are not saved
    Double_t t_soft = (*SoftPart)[8] / 1e9;  // Time in seconds
    cpg->AddTrack(static_cast<int>((*SoftPart)[0]), (*SoftPart)[1],
                  (*SoftPart)[2], (*SoftPart)[3], (*SoftPart)[5],
                  (*SoftPart)[6], (*SoftPart)[7], 0, true, (*SoftPart)[4],
                  t_soft, w);
  }
 
  return kTRUE;
}

SetPositions()

void MuDISGenerator::SetPositions ( Double_t  z_start, Double_t  z_end  )

inline

Definition at line 34 of file MuDISGenerator.h.

                                                      {
    startZ = z_start;
    endZ = z_end;
  }

Member Data Documentation

dPart

TClonesArray* MuDISGenerator::dPart

protected

Definition at line 42 of file MuDISGenerator.h.

dPartSoft

TClonesArray* MuDISGenerator::dPartSoft

protected

Definition at line 43 of file MuDISGenerator.h.

endZ

Double_t MuDISGenerator::endZ

protected

Definition at line 40 of file MuDISGenerator.h.

fFirst

bool MuDISGenerator::fFirst

protected

Definition at line 49 of file MuDISGenerator.h.

fInputFile

TFile* MuDISGenerator::fInputFile

protected

don't make it persistent, magic ROOT command

Definition at line 45 of file MuDISGenerator.h.

fLogger

FairLogger* MuDISGenerator::fLogger

protected

Definition at line 44 of file MuDISGenerator.h.

fn

int MuDISGenerator::fn

protected

Definition at line 48 of file MuDISGenerator.h.

fNevents

int MuDISGenerator::fNevents

protected

Definition at line 47 of file MuDISGenerator.h.

fTree

TTree* MuDISGenerator::fTree

protected

Definition at line 46 of file MuDISGenerator.h.

iMuon

TClonesArray* MuDISGenerator::iMuon

protected

Definition at line 41 of file MuDISGenerator.h.

startZ

Double_t MuDISGenerator::startZ

protected

Definition at line 40 of file MuDISGenerator.h.

---

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

• shipgen/MuDISGenerator.h
• shipgen/MuDISGenerator.cxx