g4simulation/g4tpc/PHG4TpcElectronDrift.cc

0001 // this is the new containers version
0002 // it uses the same MapToPadPlane as the old containers version
0003
0004 #include "PHG4TpcElectronDrift.h"
0005 #include "PHG4TpcDistortion.h"
0006 #include "PHG4TpcPadPlane.h"  // for PHG4TpcPadPlane
0007 #include "TpcClusterBuilder.h"
0008
0009 #include <trackbase/ClusHitsVerbosev1.h>
0010 #include <trackbase/TpcDefs.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrClusterContainerv4.h>
0014 #include <trackbase/TrkrDefs.h>
0015 #include <trackbase/TrkrHit.h>  // for TrkrHit
0016 #include <trackbase/TrkrHitSet.h>
0017 #include <trackbase/TrkrHitSetContainerv1.h>
0018 #include <trackbase/TrkrHitTruthAssoc.h>  // for TrkrHitTruthA...
0019 #include <trackbase/TrkrHitTruthAssocv1.h>
0020 #include <trackbase/TrkrHitv2.h>
0021
0022 #include <g4tracking/TrkrTruthTrackContainerv1.h>
0023 #include <g4tracking/TrkrTruthTrackv1.h>
0024
0025 #include <g4detectors/PHG4TpcGeom.h>
0026 #include <g4detectors/PHG4TpcGeomContainer.h>
0027
0028 #include <g4main/PHG4Hit.h>
0029 #include <g4main/PHG4HitContainer.h>
0030 #include <g4main/PHG4Particlev3.h>
0031 #include <g4main/PHG4TruthInfoContainer.h>
0032
0033 #include <phparameter/PHParameterInterface.h>  // for PHParameterIn...
0034 #include <phparameter/PHParameters.h>
0035 #include <phparameter/PHParametersContainer.h>
0036
0037 #include <pdbcalbase/PdbParameterMapContainer.h>
0038
0039 #include <fun4all/Fun4AllReturnCodes.h>
0040 #include <fun4all/Fun4AllServer.h>
0041 #include <fun4all/SubsysReco.h>  // for SubsysReco
0042
0043 #include <phool/PHCompositeNode.h>
0044 #include <phool/PHDataNode.h>  // for PHDataNode
0045 #include <phool/PHIODataNode.h>
0046 #include <phool/PHNode.h>  // for PHNode
0047 #include <phool/PHNodeIterator.h>
0048 #include <phool/PHObject.h>  // for PHObject
0049 #include <phool/PHRandomSeed.h>
0050 #include <phool/getClass.h>
0051 #include <phool/phool.h>  // for PHWHERE
0052
0053 #include <TFile.h>
0054 #include <TH1.h>
0055 #include <TH2.h>
0056 #include <TNtuple.h>
0057 #include <TSystem.h>
0058
0059 #include <gsl/gsl_randist.h>
0060 #include <gsl/gsl_rng.h>  // for gsl_rng_alloc
0061
0062 #include <array>
0063 #include <cassert>
0064 #include <cmath>    // for sqrt, abs, NAN
0065 #include <cstdlib>  // for exit
0066 #include <format>
0067 #include <iostream>
0068 #include <map>      // for _Rb_tree_cons...
0069 #include <utility>  // for pair
0070
0071 namespace
0072 {
0073   template <class T>
0074   constexpr T square(const T &x)
0075   {
0076     return x * x;
0077   }
0078 }  // namespace
0079
0080 PHG4TpcElectronDrift::PHG4TpcElectronDrift(const std::string &name)
0081   : SubsysReco(name)
0082   , PHParameterInterface(name)
0083   , temp_hitsetcontainer(new TrkrHitSetContainerv1)
0084   , single_hitsetcontainer(new TrkrHitSetContainerv1)
0085 {
0086   InitializeParameters();
0087   RandomGenerator.reset(gsl_rng_alloc(gsl_rng_mt19937));
0088   set_seed(PHRandomSeed());
0089 }
0090
0091 //_____________________________________________________________
0092 int PHG4TpcElectronDrift::Init(PHCompositeNode *topNode)
0093 {
0094   padplane->Init(topNode);
0095   event_num = 0;
0096   return Fun4AllReturnCodes::EVENT_OK;
0097 }
0098
0099 //_____________________________________________________________
0100 int PHG4TpcElectronDrift::InitRun(PHCompositeNode *topNode)
0101 {
0102   PHNodeIterator iter(topNode);
0103
0104   // Looking for the DST node
0105   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0106   if (!dstNode)
0107   {
0108     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0109     exit(1);
0110   }
0111   auto *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0112   auto *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0113   const std::string paramnodename = "G4CELLPARAM_" + detector;
0114   const std::string geonodename = "G4CELLPAR_" + detector;
0115   const std::string tpcgeonodename = "G4GEO_" + detector;
0116   hitnodename = "G4HIT_" + detector;
0117   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0118   if (!g4hit)
0119   {
0120     std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
0121     topNode->print();
0122     gSystem->Exit(1);
0123     exit(1);
0124   }
0125   // new containers
0126   hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0127   if (!hitsetcontainer)
0128   {
0129     PHNodeIterator dstiter(dstNode);
0130     auto *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0131     if (!DetNode)
0132     {
0133       DetNode = new PHCompositeNode("TRKR");
0134       dstNode->addNode(DetNode);
0135     }
0136
0137     hitsetcontainer = new TrkrHitSetContainerv1;
0138     auto *newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
0139     DetNode->addNode(newNode);
0140   }
0141
0142   hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0143   if (!hittruthassoc)
0144   {
0145     PHNodeIterator dstiter(dstNode);
0146     auto *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0147     if (!DetNode)
0148     {
0149       DetNode = new PHCompositeNode("TRKR");
0150       dstNode->addNode(DetNode);
0151     }
0152
0153     hittruthassoc = new TrkrHitTruthAssocv1;
0154     auto *newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
0155     DetNode->addNode(newNode);
0156   }
0157
0158   truthtracks = findNode::getClass<TrkrTruthTrackContainer>(topNode, "TRKR_TRUTHTRACKCONTAINER");
0159   if (!truthtracks)
0160   {
0161     PHNodeIterator dstiter(dstNode);
0162     auto *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0163     if (!DetNode)
0164     {
0165       DetNode = new PHCompositeNode("TRKR");
0166       dstNode->addNode(DetNode);
0167     }
0168
0169     truthtracks = new TrkrTruthTrackContainerv1();
0170     auto *newNode = new PHIODataNode<PHObject>(truthtracks, "TRKR_TRUTHTRACKCONTAINER", "PHObject");
0171     DetNode->addNode(newNode);
0172   }
0173
0174   truthclustercontainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_TRUTHCLUSTERCONTAINER");
0175   if (!truthclustercontainer)
0176   {
0177     PHNodeIterator dstiter(dstNode);
0178     auto *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0179     if (!DetNode)
0180     {
0181       DetNode = new PHCompositeNode("TRKR");
0182       dstNode->addNode(DetNode);
0183     }
0184
0185     truthclustercontainer = new TrkrClusterContainerv4;
0186     auto *newNode = new PHIODataNode<PHObject>(truthclustercontainer, "TRKR_TRUTHCLUSTERCONTAINER", "PHObject");
0187     DetNode->addNode(newNode);
0188   }
0189
0190   seggeonodename = "TPCGEOMCONTAINER";  // + detector;
0191   seggeo = findNode::getClass<PHG4TpcGeomContainer>(topNode, seggeonodename);
0192   assert(seggeo);
0193
0194     // the z geometry is the same for all layers
0195   PHG4TpcGeom *layergeom = seggeo->GetLayerCellGeom(20);
0196   // from top of GEM stack to top of GEM stack
0197   tpc_length = 2 * (layergeom->get_max_driftlength() + layergeom->get_CM_halfwidth());
0198
0199   UpdateParametersWithMacro();
0200   PHNodeIterator runIter(runNode);
0201   auto *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
0202   if (!RunDetNode)
0203   {
0204     RunDetNode = new PHCompositeNode(detector);
0205     runNode->addNode(RunDetNode);
0206   }
0207   SaveToNodeTree(RunDetNode, paramnodename);
0208
0209   // save this to the parNode for use
0210   PHNodeIterator parIter(parNode);
0211   auto *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
0212   if (!ParDetNode)
0213   {
0214     ParDetNode = new PHCompositeNode(detector);
0215     parNode->addNode(ParDetNode);
0216   }
0217   PutOnParNode(ParDetNode, geonodename);
0218
0219   // find Tpc Geo
0220   PHNodeIterator tpcpariter(ParDetNode);
0221   auto *tpcparams = findNode::getClass<PHParametersContainer>(ParDetNode, tpcgeonodename);
0222   if (!tpcparams)
0223   {
0224     const std::string runparamname = "G4GEOPARAM_" + detector;
0225     auto *tpcpdbparams = findNode::getClass<PdbParameterMapContainer>(RunDetNode, runparamname);
0226     if (tpcpdbparams)
0227     {
0228       tpcparams = new PHParametersContainer(detector);
0229       if (Verbosity())
0230       {
0231         tpcpdbparams->print();
0232       }
0233       tpcparams->CreateAndFillFrom(tpcpdbparams, detector);
0234       ParDetNode->addNode(new PHDataNode<PHParametersContainer>(tpcparams, tpcgeonodename));
0235     }
0236     else
0237     {
0238       std::cout << "PHG4TpcElectronDrift::InitRun - failed to find " << runparamname << " in order to initialize " << tpcgeonodename << ". Aborting run ..." << std::endl;
0239       return Fun4AllReturnCodes::ABORTRUN;
0240     }
0241   }
0242   assert(tpcparams);
0243
0244   if (Verbosity())
0245   {
0246     tpcparams->Print();
0247   }
0248   const PHParameters *tpcparam = tpcparams->GetParameters(0);
0249   assert(tpcparam);
0250
0251   diffusion_long = get_double_param("diffusion_long");
0252   added_smear_sigma_long = get_double_param("added_smear_long");
0253   diffusion_trans = get_double_param("diffusion_trans");
0254   if (zero_bfield)
0255   {
0256     diffusion_trans *= zero_bfield_diffusion_factor;
0257   }
0258   added_smear_sigma_trans = get_double_param("added_smear_trans");
0259
0260   // Data on gasses @20 C and 760 Torr from the following source:
0261   // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
0262   // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
0263   // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
0264
0265   double Ne_dEdx = 1.56;  // keV/cm
0266   double Ne_NTotal = 43;  // Number/cm
0267   double Ne_frac = tpcparam->get_double_param("Ne_frac");
0268
0269   double Ar_dEdx = 2.44;  // keV/cm
0270   double Ar_NTotal = 94;  // Number/cm
0271   double Ar_frac = tpcparam->get_double_param("Ar_frac");
0272
0273   double CF4_dEdx = 7;      // keV/cm
0274   double CF4_NTotal = 100;  // Number/cm
0275   double CF4_frac = tpcparam->get_double_param("CF4_frac");
0276
0277   double N2_dEdx = 2.127;  // keV/cm https://pdg.lbl.gov/2024/AtomicNuclearProperties/HTML/nitrogen_gas.html
0278   double N2_NTotal = 25;   // Number/cm (probably not right but has a very small impact)
0279   double N2_frac = tpcparam->get_double_param("N2_frac");
0280
0281   double isobutane_dEdx = 5.93;   // keV/cm
0282   double isobutane_NTotal = 195;  // Number/cm
0283   double isobutane_frac = tpcparam->get_double_param("isobutane_frac");
0284
0285   if (m_use_PDG_gas_params)
0286   {
0287     Ne_dEdx = 1.446;
0288     Ne_NTotal = 40;
0289
0290     Ar_dEdx = 2.525;
0291     Ar_NTotal = 97;
0292
0293     CF4_dEdx = 6.382;
0294     CF4_NTotal = 120;
0295   }
0296
0297   double Tpc_NTot = (Ne_NTotal * Ne_frac) + (Ar_NTotal * Ar_frac) + (CF4_NTotal * CF4_frac) + (N2_NTotal * N2_frac) + (isobutane_NTotal * isobutane_frac);
0298
0299   double Tpc_dEdx = (Ne_dEdx * Ne_frac) + (Ar_dEdx * Ar_frac) + (CF4_dEdx * CF4_frac) + (N2_dEdx * N2_frac) + (isobutane_dEdx * isobutane_frac);
0300
0301   electrons_per_gev = (Tpc_NTot / Tpc_dEdx) * 1e6;
0302
0303   // min_time to max_time is the time window for accepting drifted electrons after the trigger
0304   min_time = 0.0;
0305   max_time = layergeom->get_max_driftlength() / layergeom->get_drift_velocity_sim() + layergeom->get_extended_readout_time();
0306   min_active_radius = get_double_param("min_active_radius");
0307   max_active_radius = get_double_param("max_active_radius");
0308
0309   if (Verbosity() > 0)
0310   {
0311     std::cout << PHWHERE << " drift velocity " << layergeom->get_drift_velocity_sim() << " extended_readout_time " << layergeom->get_extended_readout_time() << " max time cutoff " << max_time << std::endl;
0312   }
0313
0314   auto *se = Fun4AllServer::instance();
0315   dlong = new TH1F("difflong", "longitudinal diffusion", 100, diffusion_long - diffusion_long / 2., diffusion_long + diffusion_long / 2.);
0316   se->registerHisto(dlong);
0317   dtrans = new TH1F("difftrans", "transversal diffusion", 100, diffusion_trans - diffusion_trans / 2., diffusion_trans + diffusion_trans / 2.);
0318   se->registerHisto(dtrans);
0319
0320   do_ElectronDriftQAHistos = false;  // Whether or not to produce an ElectronDriftQA.root file with useful info
0321   if (do_ElectronDriftQAHistos)
0322   {
0323     hitmapstart = new TH2F("hitmapstart", "g4hit starting X-Y locations", 1560, -78, 78, 1560, -78, 78);
0324     hitmapend = new TH2F("hitmapend", "g4hit final X-Y locations", 1560, -78, 78, 1560, -78, 78);
0325     hitmapstart_z = new TH2F("hitmapstart_z", "g4hit starting Z-R locations", 2000, -100, 100, 780, 0, 78);
0326     hitmapend_z = new TH2F("hitmapend_z", "g4hit final Z-R locations", 2000, -100, 100, 780, 0, 78);
0327     z_startmap = new TH2F("z_startmap", "g4hit starting Z vs. R locations", 2000, -100, 100, 780, 0, 78);
0328     deltaphi = new TH2F("deltaphi", "Total delta phi; phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
0329     deltaRphinodiff = new TH2F("deltaRphinodiff", "Total delta R*phi, no diffusion; r (cm);#Delta R*phi (cm)", 600, 20, 80, 1000, -3, 5);
0330     deltaphivsRnodiff = new TH2F("deltaphivsRnodiff", "Total delta phi vs. R; phi (rad);#Delta phi (rad)", 600, 20, 80, 1000, -.2, .2);
0331     deltaz = new TH2F("deltaz", "Total delta z; z (cm);#Delta z (cm)", 1000, 0, 100, 1000, -.5, 5);
0332     deltaphinodiff = new TH2F("deltaphinodiff", "Total delta phi (no diffusion, only SC distortion); phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
0333     deltaphinodist = new TH2F("deltaphinodist", "Total delta phi (no SC distortion, only diffusion); phi (rad);#Delta phi (rad)", 600, -M_PI, M_PI, 1000, -.2, .2);
0334     deltar = new TH2F("deltar", "Total Delta r; r (cm);#Delta r (cm)", 580, 20, 78, 1000, -3, 5);
0335     deltarnodiff = new TH2F("deltarnodiff", "Delta r (no diffusion, only SC distortion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5);
0336     deltarnodist = new TH2F("deltarnodist", "Delta r (no SC distortion, only diffusion); r (cm);#Delta r (cm)", 580, 20, 78, 1000, -2, 5);
0337     ratioElectronsRR = new TH1F("ratioElectronsRR", "Ratio of electrons reach readout vs all in acceptance", 1561, -0.0325, 1.0465);
0338   }
0339
0340   if (Verbosity())
0341   {
0342     // eval tree only when verbosity is on
0343     m_outf.reset(new TFile("nt_out.root", "recreate"));
0344     nt = new TNtuple("nt", "electron drift stuff", "hit:ts:tb:tsig:rad:zstart:zfinal");
0345     nthit = new TNtuple("nthit", "TrkrHit collecting", "layer:phipad:zbin:neffelectrons");
0346     ntfinalhit = new TNtuple("ntfinalhit", "TrkrHit collecting", "layer:phipad:zbin:neffelectrons");
0347     ntpad = new TNtuple("ntpad", "electron by electron pad centroid", "layer:phigem:phiclus:zgem:zclus");
0348     se->registerHisto(nt);
0349     se->registerHisto(nthit);
0350     se->registerHisto(ntpad);
0351   }
0352
0353   padplane->InitRun(topNode);
0354
0355   // print all layers radii
0356   if (Verbosity())
0357   {
0358     const auto range = seggeo->get_begin_end();
0359     std::cout << "PHG4TpcElectronDrift::InitRun - layers: " << std::distance(range.first, range.second) << std::endl;
0360     int counter = 0;
0361     for (auto layeriter = range.first; layeriter != range.second; ++layeriter)
0362     {
0363       const auto radius = layeriter->second->get_radius();
0364       std::cout << std::format("{:.3f} ", radius);
0365       if (++counter == 8)
0366       {
0367         counter = 0;
0368         std::cout << std::endl;
0369       }
0370     }
0371     std::cout << std::endl;
0372   }
0373
0374   if (record_ClusHitsVerbose)
0375   {
0376     // get the node
0377     mClusHitsVerbose = findNode::getClass<ClusHitsVerbosev1>(topNode, "Trkr_TruthClusHitsVerbose");
0378     if (!mClusHitsVerbose)
0379     {
0380       PHNodeIterator dstiter(dstNode);
0381       auto *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0382       if (!DetNode)
0383       {
0384         DetNode = new PHCompositeNode("TRKR");
0385         dstNode->addNode(DetNode);
0386       }
0387       mClusHitsVerbose = new ClusHitsVerbosev1();
0388       auto *newNode = new PHIODataNode<PHObject>(mClusHitsVerbose, "Trkr_TruthClusHitsVerbose", "PHObject");
0389       DetNode->addNode(newNode);
0390     }
0391   }
0392
0393   return Fun4AllReturnCodes::EVENT_OK;
0394 }
0395
0396 int PHG4TpcElectronDrift::process_event(PHCompositeNode *topNode)
0397 {
0398   truth_track = nullptr;  // track to which truth clusters are built
0399
0400   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0401   if (!m_tGeometry)
0402   {
0403     std::cout << PHWHERE << "ActsGeometry not found on node tree. Exiting" << std::endl;
0404     return Fun4AllReturnCodes::ABORTRUN;
0405   }
0406
0407   PHG4TpcGeom *layergeom = seggeo->GetLayerCellGeom(20);
0408
0409   if (truth_clusterer.needs_input_nodes())
0410   {
0411     truth_clusterer.set_input_nodes(truthclustercontainer, m_tGeometry,
0412                                     seggeo, mClusHitsVerbose);
0413   }
0414
0415   static constexpr unsigned int print_layer = 18;
0416
0417   // tells m_distortionMap which event to look at
0418   if (m_distortionMap)
0419   {
0420     m_distortionMap->load_event(event_num);
0421   }
0422
0423   // g4hits
0424   auto *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0425   if (!g4hit)
0426   {
0427     std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
0428     gSystem->Exit(1);
0429   }
0430   PHG4TruthInfoContainer *truthinfo =
0431       findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0432
0433   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0434   unsigned int count_g4hits = 0;
0435   //  int count_electrons = 0;
0436
0437   //  double ecollectedhits = 0.0;
0438   //  int ncollectedhits = 0;
0439   double ihit = 0;
0440   unsigned int dump_interval = 5000;  // dump temp_hitsetcontainer to the node tree after this many g4hits
0441   unsigned int dump_counter = 0;
0442
0443   int trkid = -1;
0444
0445   PHG4Hit *prior_g4hit = nullptr;  // used to check for jumps in g4hits;
0446   // if there is a big jump (such as crossing into the INTT area or out of the TPC)
0447   // then cluster the truth clusters before adding a new hit. This prevents
0448   // clustering loopers in the same HitSetKey surfaces in multiple passes
0449   for (auto hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0450   {
0451     count_g4hits++;
0452     dump_counter++;
0453
0454     const double t0 = std::fmax(hiter->second->get_t(0), hiter->second->get_t(1));
0455     if (t0 > max_time)
0456     {
0457       continue;
0458     }
0459
0460     int trkid_new = hiter->second->get_trkid();
0461     if (trkid != trkid_new)
0462     {  // starting a new track
0463       prior_g4hit = nullptr;
0464       if (truth_track)
0465       {
0466         truth_clusterer.cluster_hits(truth_track);
0467       }
0468       trkid = trkid_new;
0469
0470       if (Verbosity() > 1000)
0471       {
0472         std::cout << " New track : " << trkid << " is embed? : ";
0473       }
0474
0475       if (truthinfo->isEmbeded(trkid))
0476       {
0477         truth_track = truthtracks->getTruthTrack(trkid, truthinfo);
0478         truth_clusterer.b_collect_hits = true;
0479         if (Verbosity() > 1000)
0480         {
0481           std::cout << " YES embedded" << std::endl;
0482         }
0483       }
0484       else
0485       {
0486         truth_track = nullptr;
0487         truth_clusterer.b_collect_hits = false;
0488         if (Verbosity() > 1000)
0489         {
0490           std::cout << " NOT embedded" << std::endl;
0491         }
0492       }
0493     }
0494
0495     // see if there is a jump in x or y relative to previous PHG4Hit
0496     if (truth_clusterer.b_collect_hits)
0497     {
0498       if (prior_g4hit)
0499       {
0500         // if the g4hits jump in x or y by > max_g4hit_jump, cluster the truth tracks
0501         if (std::abs(prior_g4hit->get_x(0) - hiter->second->get_x(0)) > max_g4hitstep || std::abs(prior_g4hit->get_y(0) - hiter->second->get_y(0)) > max_g4hitstep)
0502         {
0503           if (truth_track)
0504           {
0505             truth_clusterer.cluster_hits(truth_track);
0506           }
0507         }
0508       }
0509       prior_g4hit = hiter->second;
0510     }
0511
0512     // for very high occupancy events, accessing the TrkrHitsets on the node tree
0513     // for every drifted electron seems to be very slow
0514     // Instead, use a temporary map to accumulate the charge from all
0515     // drifted electrons, then copy to the node tree later
0516
0517     double eion = hiter->second->get_eion();
0518     unsigned int n_electrons = gsl_ran_poisson(RandomGenerator.get(), eion * electrons_per_gev);
0519     //    count_electrons += n_electrons;
0520
0521     if (Verbosity() > 100)
0522     {
0523       std::cout << "  new hit with t0, " << t0 << " g4hitid " << hiter->first
0524                 << " eion " << eion << " n_electrons " << n_electrons
0525                 << " entry z " << hiter->second->get_z(0) << " exit z "
0526                 << hiter->second->get_z(1) << " avg z"
0527                 << (hiter->second->get_z(0) + hiter->second->get_z(1)) / 2.0
0528                 << std::endl;
0529     }
0530
0531     if (n_electrons == 0)
0532     {
0533       continue;
0534     }
0535
0536     if (Verbosity() > 100)
0537     {
0538       std::cout << std::endl
0539                 << "electron drift: g4hit " << hiter->first << " created electrons: "
0540                 << n_electrons << " from " << eion * 1000000 << " keV" << std::endl;
0541       std::cout << " entry x,y,z = " << hiter->second->get_x(0) << "  "
0542                 << hiter->second->get_y(0) << "  " << hiter->second->get_z(0)
0543                 << " radius " << sqrt(pow(hiter->second->get_x(0), 2) + pow(hiter->second->get_y(0), 2)) << std::endl;
0544       std::cout << " exit x,y,z = " << hiter->second->get_x(1) << "  "
0545                 << hiter->second->get_y(1) << "  " << hiter->second->get_z(1)
0546                 << " radius " << sqrt(pow(hiter->second->get_x(1), 2) + pow(hiter->second->get_y(1), 2)) << std::endl;
0547     }
0548
0549     int notReachingReadout = 0;
0550     //    int notInAcceptance = 0;
0551     for (unsigned int i = 0; i < n_electrons; i++)
0552     {
0553       // We choose the electron starting position at random from a flat
0554       // distribution along the path length the parameter t is the fraction of
0555       // the distance along the path betwen entry and exit points, it has
0556       // values between 0 and 1
0557       const double f = gsl_ran_flat(RandomGenerator.get(), 0.0, 1.0);
0558
0559       const double x_start = hiter->second->get_x(0) + f * (hiter->second->get_x(1) - hiter->second->get_x(0));
0560       const double y_start = hiter->second->get_y(0) + f * (hiter->second->get_y(1) - hiter->second->get_y(0));
0561       const double z_start = hiter->second->get_z(0) + f * (hiter->second->get_z(1) - hiter->second->get_z(0));
0562       const double t_start = hiter->second->get_t(0) + f * (hiter->second->get_t(1) - hiter->second->get_t(0));
0563
0564       unsigned int side = 0;
0565       if (z_start > 0)
0566       {
0567         side = 1;
0568       }
0569
0570       const double r_sigma = diffusion_trans * sqrt(tpc_length / 2. - std::abs(z_start));
0571       const double rantrans =
0572           gsl_ran_gaussian(RandomGenerator.get(), r_sigma) +
0573           gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_trans);
0574
0575       const double t_path = (tpc_length / 2. - std::abs(z_start)) / layergeom->get_drift_velocity_sim();
0576       const double t_sigma = diffusion_long * sqrt(tpc_length / 2. - std::abs(z_start)) / layergeom->get_drift_velocity_sim();
0577       const double rantime =
0578           gsl_ran_gaussian(RandomGenerator.get(), t_sigma) +
0579     gsl_ran_gaussian(RandomGenerator.get(), added_smear_sigma_long) / layergeom->get_drift_velocity_sim();
0580       double t_final = t_start + t_path + rantime;
0581
0582       if (t_final < min_time || t_final > max_time)
0583       {
0584         continue;
0585       }
0586
0587       double z_final;
0588       if (z_start < 0)
0589       {
0590         z_final = -tpc_length / 2. + t_final * layergeom->get_drift_velocity_sim();
0591       }
0592       else
0593       {
0594         z_final = tpc_length / 2. - t_final * layergeom->get_drift_velocity_sim();
0595       }
0596
0597       const double radstart = std::sqrt(square(x_start) + square(y_start));
0598       const double phistart = std::atan2(y_start, x_start);
0599       const double ranphi = gsl_ran_flat(RandomGenerator.get(), -M_PI, M_PI);
0600
0601       double x_final = x_start + rantrans * std::cos(ranphi);  // Initialize these to be only diffused first, will be overwritten if doing SC distortion
0602       double y_final = y_start + rantrans * std::sin(ranphi);
0603
0604       double rad_final = sqrt(square(x_final) + square(y_final));
0605       double phi_final = atan2(y_final, x_final);
0606
0607       if (do_ElectronDriftQAHistos)
0608       {
0609         z_startmap->Fill(z_start, radstart);                   // map of starting location in Z vs. R
0610         deltaphinodist->Fill(phistart, rantrans / rad_final);  // delta phi no distortion, just diffusion+smear
0611         deltarnodist->Fill(radstart, rantrans);                // delta r no distortion, just diffusion+smear
0612       }
0613
0614       if (m_distortionMap)
0615       {
0616         // zhangcanyu
0617         const double reaches = m_distortionMap->get_reaches_readout(radstart, phistart, z_start);
0618         if (reaches < thresholdforreachesreadout)
0619         {
0620           notReachingReadout++;
0621           continue;
0622         }
0623
0624         const double r_distortion = m_distortionMap->get_r_distortion(radstart, phistart, z_start);
0625         const double phi_distortion = m_distortionMap->get_rphi_distortion(radstart, phistart, z_start) / radstart;
0626         const double z_distortion = m_distortionMap->get_z_distortion(radstart, phistart, z_start);
0627
0628         rad_final += r_distortion;
0629         phi_final += phi_distortion;
0630         z_final += z_distortion;
0631         if (z_start < 0)
0632         {
0633           t_final = (z_final + tpc_length / 2.0) / layergeom->get_drift_velocity_sim();
0634         }
0635         else
0636         {
0637           t_final = (tpc_length / 2.0 - z_final) / layergeom->get_drift_velocity_sim();
0638         }
0639
0640         x_final = rad_final * std::cos(phi_final);
0641         y_final = rad_final * std::sin(phi_final);
0642
0643         //  if(i < 1)
0644         //{std::cout << " electron " << i << " r_distortion " << r_distortion << " phi_distortion " << phi_distortion << " rad_final " << rad_final << " phi_final " << phi_final << " r*dphi distortion " << rad_final * phi_distortion << " z_distortion " << z_distortion << std::endl;}
0645
0646         if (do_ElectronDriftQAHistos)
0647         {
0648           const double phi_final_nodiff = phistart + phi_distortion;
0649           const double rad_final_nodiff = radstart + r_distortion;
0650           deltarnodiff->Fill(radstart, rad_final_nodiff - radstart);    // delta r no diffusion, just distortion
0651           deltaphinodiff->Fill(phistart, phi_final_nodiff - phistart);  // delta phi no diffusion, just distortion
0652           deltaphivsRnodiff->Fill(radstart, phi_final_nodiff - phistart);
0653           deltaRphinodiff->Fill(radstart, rad_final_nodiff * phi_final_nodiff - radstart * phistart);
0654
0655           // Fill Diagnostic plots, written into ElectronDriftQA.root
0656           hitmapstart->Fill(x_start, y_start);  // G4Hit starting positions
0657           hitmapend->Fill(x_final, y_final);    // INcludes diffusion and distortion
0658           hitmapstart_z->Fill(z_start, radstart);
0659           hitmapend_z->Fill(z_final, rad_final);
0660           deltar->Fill(radstart, rad_final - radstart);    // total delta r
0661           deltaphi->Fill(phistart, phi_final - phistart);  // total delta phi
0662           deltaz->Fill(z_start, z_distortion);             // map of distortion in Z (time)
0663         }
0664       }
0665
0666       // remove electrons outside of our acceptance. Careful though, electrons from just inside 30 cm can contribute in the 1st active layer readout, so leave a little margin
0667       if (rad_final < min_active_radius - 2.0 || rad_final > max_active_radius + 1.0)
0668       {
0669         //        notInAcceptance++;
0670         continue;
0671       }
0672
0673       if (Verbosity() > 1000)
0674       //      if(i < 1)
0675       {
0676         std::cout << "electron " << i << " g4hitid " << hiter->first << " f " << f << std::endl;
0677         std::cout << "radstart " << radstart << " x_start: " << x_start
0678                   << ", y_start: " << y_start
0679                   << ",z_start: " << z_start
0680                   << " t_start " << t_start
0681                   << " t_path " << t_path
0682                   << " t_sigma " << t_sigma
0683                   << " rantime " << rantime
0684                   << std::endl;
0685
0686         std::cout << "       rad_final " << rad_final << " x_final " << x_final
0687                   << " y_final " << y_final
0688                   << " z_final " << z_final << " t_final " << t_final
0689                   << " zdiff " << z_final - z_start << std::endl;
0690       }
0691
0692       if (Verbosity() > 0)
0693       {
0694         assert(nt);
0695         nt->Fill(ihit, t_start, t_final, t_sigma, rad_final, z_start, z_final);
0696       }
0697       padplane->MapToPadPlane(truth_clusterer, single_hitsetcontainer.get(),
0698                               temp_hitsetcontainer.get(), hittruthassoc, x_final, y_final, t_final,
0699                               side, hiter, ntpad, nthit);
0700     }  // end loop over electrons for this g4hit
0701
0702     if (do_ElectronDriftQAHistos)
0703     {
0704       ratioElectronsRR->Fill((double) (n_electrons - notReachingReadout) / n_electrons);
0705     }
0706
0707     TrkrHitSetContainer::ConstRange single_hitset_range = single_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
0708     for (TrkrHitSetContainer::ConstIterator single_hitset_iter = single_hitset_range.first;
0709          single_hitset_iter != single_hitset_range.second;
0710          ++single_hitset_iter)
0711     {
0712       // we have an itrator to one TrkrHitSet for the Tpc from the single_hitsetcontainer
0713       TrkrDefs::hitsetkey node_hitsetkey = single_hitset_iter->first;
0714       const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey);
0715       const int sector = TpcDefs::getSectorId(node_hitsetkey);
0716       const int side = TpcDefs::getSide(node_hitsetkey);
0717
0718       if (Verbosity() > 8)
0719       {
0720         std::cout << " hitsetkey " << node_hitsetkey << " layer " << layer << " sector " << sector << " side " << side << std::endl;
0721       }
0722       // get all of the hits from the single hitset
0723       TrkrHitSet::ConstRange single_hit_range = single_hitset_iter->second->getHits();
0724       for (TrkrHitSet::ConstIterator single_hit_iter = single_hit_range.first;
0725            single_hit_iter != single_hit_range.second;
0726            ++single_hit_iter)
0727       {
0728         TrkrDefs::hitkey single_hitkey = single_hit_iter->first;
0729
0730         // Add the hit-g4hit association
0731         // no need to check for duplicates, since the hit is new
0732         hittruthassoc->addAssoc(node_hitsetkey, single_hitkey, hiter->first);
0733         if (Verbosity() > 100)
0734         {
0735           std::cout << "        adding assoc for node_hitsetkey " << node_hitsetkey << " single_hitkey " << single_hitkey << " g4hitkey " << hiter->first << std::endl;
0736         }
0737       }
0738     }
0739
0740     // Dump the temp_hitsetcontainer to the node tree and reset it
0741     //    - after every "dump_interval" g4hits
0742     //    - if this is the last g4hit
0743     if (dump_counter >= dump_interval || count_g4hits == g4hit->size())
0744     {
0745       // std::cout << " dump_counter " << dump_counter << " count_g4hits " << count_g4hits << std::endl;
0746
0747       double eg4hit = 0.0;
0748       TrkrHitSetContainer::ConstRange temp_hitset_range = temp_hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
0749       for (TrkrHitSetContainer::ConstIterator temp_hitset_iter = temp_hitset_range.first;
0750            temp_hitset_iter != temp_hitset_range.second;
0751            ++temp_hitset_iter)
0752       {
0753         // we have an itrator to one TrkrHitSet for the Tpc from the temp_hitsetcontainer
0754         TrkrDefs::hitsetkey node_hitsetkey = temp_hitset_iter->first;
0755         const unsigned int layer = TrkrDefs::getLayer(node_hitsetkey);
0756         const int sector = TpcDefs::getSectorId(node_hitsetkey);
0757         const int side = TpcDefs::getSide(node_hitsetkey);
0758         if (Verbosity() > 100)
0759         {
0760           std::cout << "PHG4TpcElectronDrift: temp_hitset with key: " << node_hitsetkey << " in layer " << layer
0761                     << " with sector " << sector << " side " << side << std::endl;
0762         }
0763
0764         // find or add this hitset on the node tree
0765         TrkrHitSetContainer::Iterator node_hitsetit = hitsetcontainer->findOrAddHitSet(node_hitsetkey);
0766
0767         // get all of the hits from the temporary hitset
0768         TrkrHitSet::ConstRange temp_hit_range = temp_hitset_iter->second->getHits();
0769         for (TrkrHitSet::ConstIterator temp_hit_iter = temp_hit_range.first;
0770              temp_hit_iter != temp_hit_range.second;
0771              ++temp_hit_iter)
0772         {
0773           TrkrDefs::hitkey temp_hitkey = temp_hit_iter->first;
0774           TrkrHit *temp_tpchit = temp_hit_iter->second;
0775           if (Verbosity() > 10 && layer == print_layer)
0776           {
0777             std::cout << "      temp_hitkey " << temp_hitkey << " layer " << layer << " pad " << TpcDefs::getPad(temp_hitkey)
0778                       << " z bin " << TpcDefs::getTBin(temp_hitkey)
0779                       << "  energy " << temp_tpchit->getEnergy() << " eg4hit " << eg4hit << std::endl;
0780
0781             eg4hit += temp_tpchit->getEnergy();
0782             //            ecollectedhits += temp_tpchit->getEnergy();
0783             //            ncollectedhits++;
0784           }
0785
0786           // find or add this hit to the node tree
0787           TrkrHit *node_hit = node_hitsetit->second->getHit(temp_hitkey);
0788           if (!node_hit)
0789           {
0790             // Otherwise, create a new one
0791             node_hit = new TrkrHitv2();
0792             node_hitsetit->second->addHitSpecificKey(temp_hitkey, node_hit);
0793           }
0794
0795           // Either way, add the energy to it
0796           node_hit->addEnergy(temp_tpchit->getEnergy());
0797
0798         }  // end loop over temp hits
0799
0800         if (Verbosity() > 100 && layer == print_layer)
0801         {
0802           std::cout << "  ihit " << ihit << " collected energy = " << eg4hit << std::endl;
0803         }
0804
0805       }  // end loop over temp hitsets
0806
0807       // erase all entries in the temp hitsetcontainer
0808       temp_hitsetcontainer->Reset();
0809
0810       // reset the dump counter
0811       dump_counter = 0;
0812     }  // end copy of temp hitsetcontainer to node tree hitsetcontainer
0813
0814     ++ihit;
0815
0816     single_hitsetcontainer->Reset();
0817
0818   }  // end loop over g4hits
0819
0820   if (truth_track)
0821   {
0822     truth_clusterer.cluster_hits(truth_track);
0823   }
0824   truth_clusterer.clear_hitsetkey_cnt();
0825
0826   if (Verbosity() > 20)
0827   {
0828     std::cout << "From PHG4TpcElectronDrift: hitsetcontainer printout at end:" << std::endl;
0829     // We want all hitsets for the Tpc
0830     TrkrHitSetContainer::ConstRange hitset_range = hitsetcontainer->getHitSets(TrkrDefs::TrkrId::tpcId);
0831     for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
0832          hitset_iter != hitset_range.second;
0833          ++hitset_iter)
0834     {
0835       // we have an itrator to one TrkrHitSet for the Tpc from the trkrHitSetContainer
0836       TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
0837       const unsigned int layer = TrkrDefs::getLayer(hitsetkey);
0838       if (layer != print_layer)
0839       {
0840         continue;
0841       }
0842       const int sector = TpcDefs::getSectorId(hitsetkey);
0843       const int side = TpcDefs::getSide(hitsetkey);
0844
0845       std::cout << "PHG4TpcElectronDrift: hitset with key: " << hitsetkey << " in layer " << layer << " with sector " << sector << " side " << side << std::endl;
0846
0847       // get all of the hits from this hitset
0848       TrkrHitSet *hitset = hitset_iter->second;
0849       TrkrHitSet::ConstRange hit_range = hitset->getHits();
0850       for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
0851            hit_iter != hit_range.second;
0852            ++hit_iter)
0853       {
0854         TrkrDefs::hitkey hitkey = hit_iter->first;
0855         TrkrHit *tpchit = hit_iter->second;
0856         std::cout << "      hitkey " << hitkey << " pad " << TpcDefs::getPad(hitkey) << " z bin " << TpcDefs::getTBin(hitkey)
0857                   << "  energy " << tpchit->getEnergy() << std::endl;
0858       }
0859     }
0860   }
0861
0862   if (Verbosity() > 1000)
0863   {
0864     std::cout << "From PHG4TpcElectronDrift: hittruthassoc dump:" << std::endl;
0865     hittruthassoc->identify();
0866
0867     hittruthassoc->identify();
0868   }
0869
0870   ++event_num;  // if doing more than one event, event_num will be incremented.
0871
0872   if (Verbosity() > 500)
0873   {
0874     std::cout << " TruthTrackContainer results at end of event in PHG4TpcElectronDrift::process_event " << std::endl;
0875     truthtracks->identify();
0876   }
0877
0878   if (Verbosity() > 800)
0879   {
0880     truth_clusterer.print(truthtracks);
0881     truth_clusterer.print_file(truthtracks, "drift_clusters.txt");
0882   }
0883
0884   return Fun4AllReturnCodes::EVENT_OK;
0885 }
0886
0887 int PHG4TpcElectronDrift::End(PHCompositeNode * /*topNode*/)
0888 {
0889   if (Verbosity() > 0)
0890   {
0891     assert(m_outf);
0892     assert(nt);
0893     assert(ntpad);
0894     assert(nthit);
0895     assert(ntfinalhit);
0896
0897     m_outf->cd();
0898     nt->Write();
0899     ntpad->Write();
0900     nthit->Write();
0901     ntfinalhit->Write();
0902     m_outf->Close();
0903   }
0904   if (do_ElectronDriftQAHistos)
0905   {
0906     EDrift_outf.reset(new TFile("ElectronDriftQA.root", "recreate"));
0907     EDrift_outf->cd();
0908     deltar->Write();
0909     deltaphi->Write();
0910     deltaz->Write();
0911     deltarnodist->Write();
0912     deltaphinodist->Write();
0913     deltarnodiff->Write();
0914     deltaphinodiff->Write();
0915     deltaRphinodiff->Write();
0916     deltaphivsRnodiff->Write();
0917     hitmapstart->Write();
0918     hitmapend->Write();
0919     hitmapstart_z->Write();
0920     hitmapend_z->Write();
0921     z_startmap->Write();
0922     ratioElectronsRR->Write();
0923     EDrift_outf->Close();
0924   }
0925   return Fun4AllReturnCodes::EVENT_OK;
0926 }
0927
0928 void PHG4TpcElectronDrift::set_seed(const unsigned int seed)
0929 {
0930   gsl_rng_set(RandomGenerator.get(), seed);
0931 }
0932
0933 void PHG4TpcElectronDrift::SetDefaultParameters()
0934 {
0935   // longitudinal diffusion for 50:50 Ne:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008
0936   // longitudinal diffusion for 60:40 Ar:CF4 is 0.012, transverse is 0.004, drift velocity is 0.008 (chosen to be the same at 50/50 Ne:CF4)
0937   // longitudinal diffusion for 65:25:10 Ar:CF4:N2 is 0.013613, transverse is 0.005487, drift velocity is 0.006965
0938   // longitudinal diffusion for 75:20:05 Ar:CF4:i-C4H10 is 0.014596, transverse is 0.005313, drift velocity is 0.007550
0939
0940   set_default_double_param("diffusion_long", 0.014596);   // cm/SQRT(cm)
0941   set_default_double_param("diffusion_trans", 0.005313);  // cm/SQRT(cm)
0942   set_default_double_param("Ne_frac", 0.00);
0943   set_default_double_param("Ar_frac", 0.75);
0944   set_default_double_param("CF4_frac", 0.20);
0945   set_default_double_param("N2_frac", 0.00);
0946   set_default_double_param("isobutane_frac", 0.05);
0947   set_default_double_param("min_active_radius", 30.);        // cm
0948   set_default_double_param("max_active_radius", 78.);        // cm
0949
0950   // These are purely fudge factors, used to increase the resolution to 150 microns and 500 microns, respectively
0951   // override them from the macro to get a different resolution
0952   set_default_double_param("added_smear_trans", 0.0);  // cm (used to be 0.085 before sims got better)
0953   set_default_double_param("added_smear_long", 0.0);   // cm (used to be 0.105 before sims got better)
0954
0955   return;
0956 }
0957
0958 void PHG4TpcElectronDrift::setTpcDistortion(PHG4TpcDistortion *distortionMap)
0959 {
0960   m_distortionMap.reset(distortionMap);
0961 }
0962
0963 void PHG4TpcElectronDrift::registerPadPlane(PHG4TpcPadPlane *inpadplane)
0964 {
0965   if (Verbosity())
0966   {
0967     std::cout << "Registering padplane " << std::endl;
0968   }
0969   padplane.reset(inpadplane);
0970   padplane->Detector(Detector());
0971   padplane->UpdateInternalParameters();
0972   if (Verbosity())
0973   {
0974     std::cout << "padplane registered and parameters updated" << std::endl;
0975   }
0976
0977   return;
0978 }
0979
0980 void PHG4TpcElectronDrift::set_flag_threshold_distortion(bool setflag, float setthreshold)
0981 {
0982   std::cout << std::format("The logical status of threshold is now {}! and the value is set to {}", setflag, setthreshold)
0983             << std::endl
0984             << std::endl
0985             << std::endl;
0986   do_getReachReadout = setflag;
0987   thresholdforreachesreadout = setthreshold;
0988 }