g4simulation/g4main/CosmicSpray.cc

0001 // CosmicSpray class
0002 //  Author: Daniel Lis
0003 //  Brief: Particel generator Class that sources a muon with a vertex and momentum that should mimic real life
0004 //  modified by Shuhang on 03/2022: now this class serves as a wrapper class that drives "EcoMug"
0005
0006 // For the muon rate calculation, if using the default setting: time(second) = nevents * 1.434e-4, then scale the histogram by 1/time to get the rate.
0007
0008 #include "CosmicSpray.h"
0009 #include "EcoMug.h"
0010
0011 #include "PHG4InEvent.h"
0012 #include "PHG4Particle.h"
0013 #include "PHG4Particlev2.h"
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <fun4all/SubsysReco.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHDataNode.h>      // for PHDataNode
0020 #include <phool/PHNode.h>          // for PHNode
0021 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0022 #include <phool/PHObject.h>        // for PHObject
0023 #include <phool/PHRandomSeed.h>
0024 #include <phool/getClass.h>
0025
0026 #include <array>  // for array
0027 #include <cmath>
0028 #include <iostream>  // for operator<<, endl, basic_ostream
0029
0030 namespace
0031 {
0032   constexpr double muon_mass = 0.1056583745;  // pdg value 2014
0033 }
0034
0035 bool CosmicSpray::InDetector(double x, double y, double z) const
0036 {
0037   double gap = 5;
0038   if (x > _x_max)
0039   {
0040     return false;
0041   }
0042   if (x < -_x_max)
0043   {
0044     return false;
0045   }
0046   if (z > _z_max + gap)
0047   {
0048     return false;
0049   }
0050   if (z < -_z_max - gap)
0051   {
0052     return false;
0053   }
0054   if (y > _y_fix + gap)
0055   {
0056     return false;
0057   }
0058   if (y < -_y_fix - gap)
0059   {
0060     return false;
0061   }
0062   return true;
0063 }
0064
0065 CosmicSpray::CosmicSpray(const std::string &name, const double R)
0066   : SubsysReco(name)
0067   //  , _x_min(183.3)
0068   , _x_max(264.71)
0069   //  , _z_min(-304.91)
0070   , _z_max(304.91)
0071   , _y_fix(_x_max)
0072   , _R(R)
0073 {
0074   unsigned int seed = PHRandomSeed();
0075   gen.SetSeed(seed);
0076   gen.SetUseHSphere();            // half-spherical surface generation
0077   gen.SetHSphereRadius(R / 100);  // half-sphere radius
0078   gen.SetHSphereCenterPosition({{0., 0., -_y_fix / 100}});
0079   gen.SetMinimumMomentum(0.5);
0080
0081   return;
0082 }
0083
0084 int CosmicSpray::InitRun(PHCompositeNode *topNode)
0085 {
0086   PHG4InEvent *inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0087   if (!inevent)
0088   {
0089     PHNodeIterator iter(topNode);
0090     PHCompositeNode *dstNode;
0091     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0092
0093     inevent = new PHG4InEvent();
0094     PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(inevent, "PHG4INEVENT", "PHObject");
0095     dstNode->addNode(newNode);
0096   }
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int CosmicSpray::process_event(PHCompositeNode *topNode)
0101 {
0102   // set_vertex
0103   if (Verbosity() > 0)
0104   {
0105     std::cout << "Processing Event" << std::endl;
0106   }
0107   std::string pdgname = "mu-";
0108   int pdgcode = 13;
0109   int trackid = 0;
0110   double gun_t = 0.0;
0111   double gun_x = 0;
0112   double gun_y = 0;
0113   double gun_z = 0;
0114   double gun_px = 0;
0115   double gun_py = 0;
0116   double gun_pz = 0;
0117   bool GoodEvent = false;
0118   while (!GoodEvent)
0119   {
0120     gen.Generate();
0121     std::array<double, 3> muon_position = gen.GetGenerationPosition();
0122     double muon_p = gen.GetGenerationMomentum();
0123     _gun_e = sqrt(muon_mass * muon_mass + muon_p * muon_p);
0124     double tr = gen.GetGenerationTheta();
0125     double pr = gen.GetGenerationPhi();
0126     double muon_charge = gen.GetCharge();
0127
0128     if (muon_charge > 0)
0129     {
0130       pdgcode = -13;
0131       pdgname = "mu+";
0132     }
0133
0134     gun_px = muon_p * sin(tr) * sin(pr);
0135     gun_py = -1 * fabs(muon_p * cos(tr));
0136     gun_pz = muon_p * sin(tr) * cos(pr);
0137
0138     gun_x = muon_position[1] * 100;
0139     gun_y = muon_position[2] * 100;
0140     gun_z = muon_position[0] * 100;
0141
0142     // unit vectors of muon momentum
0143     double upx = gun_px / muon_p;
0144     double upy = gun_py / muon_p;
0145     double upz = gun_pz / muon_p;
0146
0147     // check if the muon goes in the detector
0148
0149     double x1 = gun_x;
0150     double y1 = gun_y;
0151     double z1 = gun_z;
0152
0153     double L = 0;
0154
0155     while (y1 > -_y_fix && L < _R)
0156     {
0157       L += 10;
0158       x1 += upx * 10;
0159       y1 += upy * 10;
0160       z1 += upz * 10;
0161
0162       if (InDetector(x1, y1, z1))
0163       {
0164         GoodEvent = true;
0165         break;
0166       }
0167     }
0168   }
0169
0170   PHG4InEvent *inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0171   int vtxindex = inevent->AddVtx(gun_x, gun_y, gun_z, gun_t);
0172
0173   PHG4Particle *particle = new PHG4Particlev2();
0174   particle->set_track_id(trackid);
0175   particle->set_vtx_id(vtxindex);
0176   particle->set_parent_id(0);
0177   particle->set_name(pdgname);
0178   particle->set_pid(pdgcode);
0179   particle->set_px(gun_px);
0180   particle->set_py(gun_py);
0181   particle->set_pz(gun_pz);
0182   particle->set_e(_gun_e);
0183   inevent->AddParticle(vtxindex, particle);
0184   if (Verbosity() > 0)
0185   {
0186     std::cout << "Momentum: " << gun_px << " / " << gun_py << " / " << gun_pz << std::endl;
0187     std::cout << "total mom: " << _gun_e << std::endl;
0188     std::cout << "track_id: " << trackid << std::endl;
0189     std::cout << "vtxindex: " << vtxindex << std::endl;
0190     std::cout << "pdgname: " << pdgname << std::endl;
0191     std::cout << "pdgcode: " << pdgcode << std::endl;
0192   }
0193   return 0;
0194 }