benchmarks/zdc_lyso/gen_particles.cxx

0001 #include "HepMC3/GenEvent.h"
0002 #include "HepMC3/ReaderAscii.h"
0003 #include "HepMC3/WriterAscii.h"
0004 #include "HepMC3/Print.h"
0005
0006 #include "TRandom3.h"
0007 #include "TVector3.h"
0008
0009 #include <iostream>
0010 #include <random>
0011 #include <cmath>
0012 #include <math.h>
0013 #include <TMath.h>
0014 #include <TDatabasePDG.h>
0015 #include <TParticlePDG.h>
0016
0017 using namespace HepMC3;
0018
0019 // Generate single electron as input to the Insert simulation.
0020 // --
0021 // We generate events with a constant polar angle with respect to
0022 // the proton direction and then rotate so that the events are given
0023 // in normal lab coordinate system
0024 // --
0025 void gen_particles(
0026                     int n_events = 10,
0027                     const char* out_fname = "gen_particles.hepmc",
0028                     TString particle_name = "e-",
0029                     double th_min = 3., // Minimum polar angle, in degrees
0030             double th_max = 3., // Maximum polar angle, in degrees
0031             double phi_min = 0., // Minimum azimuthal angle, in degrees
0032                     double phi_max = 360., // Maximum azimuthal angle, in degrees
0033                     double p = 10.,  // Momentum in GeV/c
0034             int dist = 0  //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian
0035                   )
0036 {
0037   WriterAscii hepmc_output(out_fname);
0038   int events_parsed = 0;
0039   GenEvent evt(Units::GEV, Units::MM);
0040
0041   // Random number generator
0042   TRandom3 *r1 = new TRandom3(0); //Use time as random seed
0043
0044   // Getting generated particle information
0045   TDatabasePDG *pdg = new TDatabasePDG();
0046   TParticlePDG *particle = pdg->GetParticle(particle_name);
0047   const double mass = particle->Mass();
0048   const int pdgID = particle->PdgCode();
0049
0050   for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
0051
0052     //Set the event number
0053     evt.set_event_number(events_parsed);
0054
0055     // FourVector(px,py,pz,e,pdgid,status)
0056     // type 4 is beam
0057     // pdgid 11 - electron
0058     // pdgid 111 - pi0
0059     // pdgid 2212 - proton
0060     GenParticlePtr p1 =
0061         std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
0062     GenParticlePtr p2 = std::make_shared<GenParticle>(
0063         FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
0064
0065     // Define momentum with respect to proton direction
0066     double phi   = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad());
0067     double th    = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad());
0068
0069     //Total momentum distribution
0070     double pevent = -1;
0071     if(dist==0){ //fixed
0072     pevent = p;
0073     }
0074     else if(dist==1){ //Uniform: +-50% variation
0075     pevent = p*(1. + r1->Uniform(-0.5,0.5) );
0076     }
0077     else if(dist==2){  //Gaussian: Sigma = 0.1*mean
0078     while(pevent<0) //Avoid negative values
0079         pevent = r1->Gaus(p,0.1*p);
0080     }
0081
0082     double px    = pevent * std::cos(phi) * std::sin(th);
0083     double py    = pevent * std::sin(phi) * std::sin(th);
0084     double pz    = pevent * std::cos(th);
0085     TVector3 pvec(px,py,pz);
0086
0087     //Rotate to lab coordinate system
0088     double cross_angle = -25./1000.; //in Rad
0089     TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction
0090     pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X
0091     // type 1 is final state
0092     // pdgid 11 - electron 0.510 MeV/c^2
0093     GenParticlePtr p3 = std::make_shared<GenParticle>(
0094         FourVector(
0095             pvec.X(), pvec.Y(), pvec.Z(),
0096             sqrt(pevent*pevent + (mass * mass))),
0097         pdgID, 1);
0098
0099     //If wanted, set non-zero vertex
0100     double vx = 0.;
0101     double vy = 0.;
0102     double vz = 0.;
0103     double vt = 0.;
0104
0105     GenVertexPtr v1 = std::make_shared<GenVertex>();
0106     evt.shift_position_by(FourVector(vx, vy, vz, vt));
0107     v1->add_particle_in(p1);
0108     v1->add_particle_in(p2);
0109
0110     v1->add_particle_out(p3);
0111     evt.add_vertex(v1);
0112
0113     if (events_parsed == 0) {
0114       std::cout << "First event: " << std::endl;
0115       Print::listing(evt);
0116     }
0117
0118     hepmc_output.write_event(evt);
0119     if (events_parsed % 100 == 0) {
0120       std::cout << "Event: " << events_parsed << std::endl;
0121     }
0122     evt.clear();
0123   }
0124   hepmc_output.close();
0125   std::cout << "Events parsed and written: " << events_parsed << std::endl;
0126 }