DV Running LUND for Moller Distribution
Writing LUND files
Reading the GEANT4 simulation output in,
void MollerLUND()
{
struct evt_t
{
Long_t event;
Float_t IntKE, IntMom[3],IntPos[3],FnlKE,FnlMom[3],FnlPos[3],MolKE,MolMom[3],MolPos[3];
};
ifstream in;
in.open("MollerScattering_NH3_4e8incident.dat");
evt_t evt;
Int_t nlines=0;
TFile *f = new TFile("MollerScattering_NH3_4e7incident.root","RECREATE");
TTree *tree0 = new TTree("MollerScattering_NH3_4e7","Moller data from ascii file");
FILE *f0;
f0=fopen("MollerScattering_NH3_4e7.LUND","w");
tree0->Branch("evt",&evt.event,"event/I:IntKE/F:IntPx:IntPy:IntPz:IntPosx:IntPosy:IntPosz:FnlKE:FnlPx:FnlPy:FnlPz:FnlPosx:FnlPosy:FnlPosz:
MolKE:MolPx:MolPy:MolPz:MolPosx:MolPosy:MolPosz");
while(in.good())
{
evt.event=nlines;
in >> evt.IntKE >> evt.IntMom[0] >> evt.IntMom[1] >> evt.IntMom[2] >> evt.IntPos[0]
>> evt.IntPos[1] >> evt.IntPos[2] >> evt.FnlKE >> evt.FnlMom[0] >> evt.FnlMom[1] >> evt.FnlMom[2]
>> evt.FnlPos[0] >> evt.FnlPos[1] >> evt.FnlPos[2] >> evt.MolKE >> evt.MolMom[0] >> evt.MolMom[1]
>> evt.MolMom[2] >> evt.MolPos[0] >> evt.MolPos[1] >> evt.MolPos[2]; double Nu;
double x_bj; // Bjorken-x
double y;
double W ; // missing mass
// Momentum and Energy are measured in MeV!
evt.FnlKE=sqrt(evt.FnlMom[0]*evt.FnlMom[0]+evt.FnlMom[1]*evt.FnlMom[1]+evt.FnlMom[2]*evt.FnlMom[2]+0.511*0.511);
evt.MolKE=sqrt(evt.MolMom[0]*evt.MolMom[0]+evt.MolMom[1]*evt.MolMom[1]+evt.MolMom[2]*evt.MolMom[2]+0.511*0.511);
Nu=evt.IntKE-evt.FnlKE;
Qsqrd=4*evt.IntKE*evt.FnlKE*(1-evt.FnlMom[2]/evt.FnlKE)/2; /* should be final momentum and not final energy*/
x_bj=Qsqrd/(2*0.938*Nu);
y=Nu/evt.IntKE;
W=0.938*0.938+2*0.938*Nu-Qsqrd;
if(W>0)
W=sqrt(W);
//
// MeV MUST be converted to GeV for LUND format!
double px,py,pz;
double Px,Py,Pz;
double KE,ke;
Px=evt.FnlMom[0]/1000;
Py=evt.FnlMom[1]/1000;
Pz=evt.FnlMom[2]/1000;
px=evt.MolMom[0]/1000;
py=evt.MolMom[1]/1000;
pz=evt.MolMom[2]/1000;
KE=evt.FnlKE/1000;
ke=evt.MolKE/1000;
fprintf(f0,"%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",2,1.,1.,1.,1.,x_bj,y,W,Qsqrd,Nu);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%f\t%f\t%f\t%f\t%g\t%g\t%g\t%g\n",
1,-1.,1.,11.,0.,0.,Px,Py,Pz,KE,0.000511, 0.,0.,0.);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%f\t%f\t%f\t%f\t%g\t%g\t%g\t%g\n",
2,-1.,1.,11.,0.,0.,px,py,pz,ke,0.000511, 0.,0.,0.);
nlines++;
tree0->Fill();
}
printf("NumberLines=%d\n",nlines);
tree0->Print();
tree0->Write();
delete tree0;
}
The file is written as seen below.
Corrections
The prefered coordinates for the reconstruction within GEMC would have the vertex location at (0, 0, 0). The Moller event data from GEANT4 contains the location within that simulation where the scattering occured. These variables will need to be renormalized within the LUND file to zero. Only the 14th column within the 2nd and 3rd rows out of every 3 rows needs to be changed. Using a script:
awk '{if('NR%3==2'||'NR%3==0') {$14="0"}; print $0}' MollerScattering_NH3_4e7.LUND>MollerScattering.LUND
The LUND format has a limit of 75000 lines, so to account for this by creating 21 new files with the suffix of aa-au are created:
split -l 75000 MollerScattering.LUND MollerScattering_NH3_4e7_
We will use 5 of the 21 sets of 75000 entries in GEMC trials involving different solenoid field strengths. We can use the simple command:
mv MollerScattering_NH3_4e7_aa MollerScattering_NH3_4e7_aa.LUND
Repeating this for the other 4 files.
References
[1] M. Ungaro, et al., Jefferson Lab:The LUND format, https://gemc.jlab.org/gemc/html/documentation/generator/lund.html (2016)
[2] F. Halzen, and A. Martin, Quarks and Leptons:An Introductory Course in Modern Particle Physics, John Wiley & Sons (1984)
[3] W. Tung, Bjorken scaling, http://www.scholarpedia.org/article/Bjorken_scaling (2009)
[4] L. Garren, et al., MONTE CARLO PARTICLE NUMBERING SCHEME, http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf (2007)