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 bool CosmicSpray::InDetector(double x, double y, double z)
0031 {
0032   double gap = 5;
0033   if (x > _x_max)
0034   {
0035     return false;
0036   }
0037   if (x < -_x_max)
0038   {
0039     return false;
0040   }
0041   if (z > _z_max + gap)
0042   {
0043     return false;
0044   }
0045   if (z < -_z_max - gap)
0046   {
0047     return false;
0048   }
0049   if (y > _y_fix + gap)
0050   {
0051     return false;
0052   }
0053   if (y < -_y_fix - gap)
0054   {
0055     return false;
0056   }
0057   return true;
0058 }
0059
0060 CosmicSpray::CosmicSpray(const std::string &name, const double R)
0061   : SubsysReco(name)
0062 {
0063   _x_max = 264.71;
0064   _x_min = 183.3;
0065   _z_max = 304.91;
0066   _z_min = -304.91;
0067   _y_fix = _x_max;
0068
0069   unsigned int seed = PHRandomSeed();
0070   gen.SetSeed(seed);
0071   gen.SetUseHSphere();            // half-spherical surface generation
0072   gen.SetHSphereRadius(R / 100);  // half-sphere radius
0073   gen.SetHSphereCenterPosition({{0., 0., -_y_fix / 100}});
0074   gen.SetMinimumMomentum(0.5);
0075   _R = R;
0076   return;
0077 }
0078
0079 int CosmicSpray::InitRun(PHCompositeNode *topNode)
0080 {
0081   PHG4InEvent *inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0082   if (!inevent)
0083   {
0084     PHNodeIterator iter(topNode);
0085     PHCompositeNode *dstNode;
0086     dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0087
0088     inevent = new PHG4InEvent();
0089     PHDataNode<PHObject> *newNode = new PHDataNode<PHObject>(inevent, "PHG4INEVENT", "PHObject");
0090     dstNode->addNode(newNode);
0091   }
0092   return Fun4AllReturnCodes::EVENT_OK;
0093 }
0094
0095 int CosmicSpray::process_event(PHCompositeNode *topNode)
0096 {
0097   // set_vertex
0098   if (Verbosity() > 0)
0099   {
0100     std::cout << "Processing Event" << std::endl;
0101   }
0102   std::string pdgname = "mu-";
0103   int pdgcode = 13;
0104   int trackid = 0;
0105   double gun_t = 0.0;
0106   double gun_x = 0, gun_y = 0, gun_z = 0;
0107   double gun_px = 0, gun_py = 0, gun_pz = 0;
0108   bool GoodEvent = false;
0109   while (!GoodEvent)
0110   {
0111     gen.Generate();
0112     std::array<double, 3> muon_position = gen.GetGenerationPosition();
0113     double muon_p = gen.GetGenerationMomentum();
0114     _gun_e = sqrt(0.105658 * 0.105658 + muon_p * muon_p);
0115     double tr = gen.GetGenerationTheta();
0116     double pr = gen.GetGenerationPhi();
0117     double muon_charge = gen.GetCharge();
0118
0119     if (muon_charge > 0)
0120     {
0121       pdgcode = -13;
0122       pdgname = "mu+";
0123     }
0124
0125     gun_px = muon_p * sin(tr) * sin(pr);
0126     gun_py = -1 * fabs(muon_p * cos(tr));
0127     gun_pz = muon_p * sin(tr) * cos(pr);
0128
0129     gun_x = muon_position[1] * 100;
0130     gun_y = muon_position[2] * 100;
0131     gun_z = muon_position[0] * 100;
0132
0133     // unit vectors of muon momentum
0134     double upx = gun_px / muon_p;
0135     double upy = gun_py / muon_p;
0136     double upz = gun_pz / muon_p;
0137
0138     // check if the muon goes in the detector
0139
0140     double x1 = gun_x;
0141     double y1 = gun_y;
0142     double z1 = gun_z;
0143
0144     double L = 0;
0145
0146     while (y1 > -_y_fix && L < _R)
0147     {
0148       L += 10;
0149       x1 += upx * 10;
0150       y1 += upy * 10;
0151       z1 += upz * 10;
0152
0153       if (InDetector(x1, y1, z1))
0154       {
0155         GoodEvent = true;
0156         break;
0157       }
0158     }
0159   }
0160
0161   PHG4InEvent *inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0162   int vtxindex = inevent->AddVtx(gun_x, gun_y, gun_z, gun_t);
0163
0164   PHG4Particle *particle = new PHG4Particlev2();
0165   particle->set_track_id(trackid);
0166   particle->set_vtx_id(vtxindex);
0167   particle->set_parent_id(0);
0168   particle->set_name(pdgname);
0169   particle->set_pid(pdgcode);
0170   particle->set_px(gun_px);
0171   particle->set_py(gun_py);
0172   particle->set_pz(gun_pz);
0173   particle->set_e(_gun_e);
0174   inevent->AddParticle(vtxindex, particle);
0175   if (Verbosity() > 0)
0176   {
0177     std::cout << "Momentum: " << gun_px << " / " << gun_py << " / " << gun_pz << std::endl;
0178     std::cout << "total mom: " << _gun_e << std::endl;
0179     std::cout << "track_id: " << trackid << std::endl;
0180     std::cout << "vtxindex: " << vtxindex << std::endl;
0181     std::cout << "pdgname: " << pdgname << std::endl;
0182     std::cout << "pdgcode: " << pdgcode << std::endl;
0183   }
0184   return 0;
0185 }