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