g4simulation/g4bbc/MbdDigitization.cc

0001 #include "MbdDigitization.h"
0002
0003 #include <mbd/MbdDefs.h>
0004 #include <mbd/MbdPmtHit.h>
0005 #include <mbd/MbdPmtSimContainerV1.h>
0006
0007 #include <g4main/PHG4Hit.h>
0008 #include <g4main/PHG4HitContainer.h>
0009 #include <g4main/PHG4Particle.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011 #include <g4main/PHG4VtxPoint.h>
0012
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/PHTFileServer.h>
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>
0018 #include <phool/PHNode.h>
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHObject.h>
0021 #include <phool/PHRandomSeed.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>
0024 #include <phool/sphenix_constants.h>
0025
0026 #include <TDatabasePDG.h>
0027 #include <TF1.h>
0028 #include <TFile.h>
0029 #include <TH1F.h>
0030 #include <TH2F.h>
0031 #include <TLorentzVector.h>
0032 #include <TString.h>
0033 #include <TSystem.h>
0034 #include <TTree.h>
0035
0036 #include <gsl/gsl_randist.h>
0037 #include <gsl/gsl_rng.h>
0038
0039 #include <cmath>
0040 #include <iostream>
0041
0042 //____________________________________
0043 MbdDigitization::MbdDigitization(const std::string &name)
0044   : SubsysReco(name)
0045   , m_RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0046 {
0047   f_pmtq.fill(0.);
0048   f_pmtnpe.fill(0.);
0049   f_pmtt0.fill(std::numeric_limits<decltype(f_pmtt0)::value_type>::quiet_NaN());
0050   f_pmtt1.fill(std::numeric_limits<decltype(f_pmtt1)::value_type>::quiet_NaN());
0051   m_Seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0052   gsl_rng_set(m_RandomGenerator, m_Seed);
0053 }
0054
0055 MbdDigitization::~MbdDigitization()
0056 {
0057   gsl_rng_free(m_RandomGenerator);
0058   return;
0059 }
0060
0061 //___________________________________
0062 int MbdDigitization::Init(PHCompositeNode *topNode)
0063 {
0064   if (Verbosity())
0065   {
0066     std::cout << PHWHERE << std::endl;
0067   }
0068   CreateNodes(topNode);
0069
0070   _pdg = new TDatabasePDG();  // database of PDG info on particles
0071
0072   gaussian = new TF1("gaussian", "gaus", 0, 20);
0073   gaussian->FixParameter(2, _tres);  // set sigma to timing resolution
0074
0075   return 0;
0076 }
0077
0078 //___________________________________
0079 int MbdDigitization::InitRun(PHCompositeNode *topNode)
0080 {
0081   GetNodes(topNode);
0082
0083   _gains = {135.137, 135.035, 134.861, 134.543, 135.761, 134.942, 134.053, 134.398,
0084             135.568, 134.874, 134.677, 135.614, 135.366, 134.564, 135.143, 135.412, 135.308,
0085             135.752, 134.392, 135.372, 135.025, 135.247, 135.337, 135.396, 133.553, 134.93,
0086             135.917, 135.209, 135.886, 134.914, 134.766, 135.734, 135.343, 134.815, 134.68,
0087             135.214, 135.251, 134.617, 135.182, 134.659, 135.119, 135.022, 134.648, 134.814,
0088             135.1, 134.99, 134.82, 134.844, 135.817, 135.788, 135.469, 135.178, 135.17, 134.938,
0089             135.544, 135.445, 134.872, 134.4, 134.875, 134.852, 134.23, 133.922, 135.46, 134.44,
0090             135.475, 135.716, 135.084, 136.074, 136.198, 136.226, 136.214, 136.1, 136.727,
0091             135.422, 136.258, 136.394, 135.754, 136.072, 136.897, 135.775, 135.902, 135.959,
0092             136.306, 135.846, 135.986, 135.487, 136.404, 134.988, 135.858, 136.511, 136.128,
0093             135.535, 135.984, 136.162, 135.496, 135.104, 135.762, 135.653, 135.379, 135.868,
0094             136.747, 135.606, 136.074, 135.833, 136.382, 135.796, 135.364, 135.985, 135.46,
0095             135.869, 136.068, 135.876, 136.339, 135.48, 136.435, 135.997, 136.01, 136.675,
0096             136.105, 136.423, 135.739, 136.112, 136.523, 135.883, 136.972, 136.741, 135.419, 135.407};
0097
0098   return 0;
0099 }
0100
0101 //__________________________________
0102 // Call user instructions for every event
0103 int MbdDigitization::process_event(PHCompositeNode * /*topNode*/)
0104 {
0105   if (Verbosity())
0106   {
0107     std::cout << PHWHERE << std::endl;
0108   }
0109   //**** Initialize Variables
0110
0111   // PMT data
0112   std::array<float, MbdDefs::MBD_N_PMT> len{0.};
0113   std::array<float, MbdDefs::MBD_N_PMT> edep{0.};
0114   std::array<float, MbdDefs::MBD_N_PMT> first_time{};
0115   first_time.fill(1e12);
0116   f_pmtt0.fill(1e12);
0117   f_pmtt1.fill(1e12);
0118   f_pmtq.fill(0.);
0119   f_pmtnpe.fill(0.);
0120   // Get True Vertex
0121   // NB: Currently PrimaryVertexIndex is always 1, need to figure out how to handle collision pile-up
0122   PHG4VtxPoint *vtxp = _truth_container->GetPrimaryVtx(_truth_container->GetPrimaryVertexIndex());
0123   if (vtxp != nullptr)
0124   {
0125     f_vx = vtxp->get_x();
0126     f_vy = vtxp->get_y();
0127     f_vz = vtxp->get_z();
0128     f_vt = vtxp->get_t();
0129
0130     if (Verbosity())
0131     {
0132       std::cout << "VTXP "
0133                 << "\t" << f_vx << "\t" << f_vy << "\t" << f_vz << "\t" << f_vt << std::endl;
0134     }
0135   }
0136
0137   // Go through BBC G4 hits
0138   if (Verbosity() > 10)
0139   {
0140     std::cout << "Processing BBC G4 Hits" << std::endl;
0141   }
0142
0143   TLorentzVector v4;
0144   unsigned int nhits = 0;
0145
0146   PHG4HitContainer::ConstRange bbc_hit_range = _bbchits->getHits();
0147   for (auto hit_iter = bbc_hit_range.first; hit_iter != bbc_hit_range.second; hit_iter++)
0148   {
0149     PHG4Hit *this_hit = hit_iter->second;
0150
0151     unsigned int ch = this_hit->get_layer();  // pmt channel
0152     // int arm = ch/64;;                         // south=0, north=1
0153
0154     int trkid = this_hit->get_trkid();
0155     // if ( trkid>0 && f_evt<20 ) std::cout << "TRKID " << trkid << std::endl;
0156
0157     PHG4Particle *part = _truth_container->GetParticle(trkid);
0158     v4.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e());
0159
0160     int pid = part->get_pid();
0161     TParticlePDG *partinfo = _pdg->GetParticle(pid);
0162     Double_t charge = -9999.;
0163     if (partinfo)
0164     {
0165       charge = partinfo->Charge() / 3;  // PDG gives charge in 1/3 e
0166     }
0167     else if (part->isIon())
0168     {
0169       charge = part->get_IonCharge();
0170     }
0171
0172     // get the first time
0173     if (this_hit->get_t(1) < first_time[ch])
0174     {
0175       if (std::abs(this_hit->get_t(1)) < sphenix_constants::time_between_crossings)
0176       {
0177         first_time[ch] = this_hit->get_t(1) - vtxp->get_t();
0178         Float_t dt = gsl_ran_gaussian(m_RandomGenerator, _tres);  // get fluctuation in time
0179         first_time[ch] += dt;
0180       }
0181       else
0182       {
0183         if (Verbosity())
0184         {
0185           std::cout << "BAD " << ch << "\t" << this_hit->get_t(1) << std::endl;
0186         }
0187       }
0188     }
0189
0190     edep[ch] += this_hit->get_edep();
0191
0192     // get summed path length for particles that can create CKOV light
0193     // n.p.e. is determined from path length
0194     Double_t beta = v4.Beta();
0195     if (beta > MbdDefs::v_ckov && charge != 0.)
0196     {
0197       len[ch] += this_hit->get_path_length();
0198
0199       _pids[pid] += 1;
0200     }
0201
0202     nhits++;
0203   }
0204
0205   if (Verbosity() > 10)
0206   {
0207     std::cout << "Found " << nhits << " MBD hits" << std::endl;
0208     std::cout << "Calculating response and storing in MbdPmtHits" << std::endl;
0209   }
0210
0211   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0212   {
0213     // Fill charge and time info
0214     if (len[ipmt] > 0.)
0215     {
0216       if (Verbosity() > 0)
0217       {
0218         std::cout << "ipmt " << ipmt << "\t" << len[ipmt] << "\t" << edep[ipmt] << std::endl;
0219       }
0220
0221       // Get charge in BBC tube
0222       f_pmtnpe[ipmt] = len[ipmt] * (120 / 3.0);                                     // MBD/BBC gets 120 p.e. per 3 cm
0223       float dnpe = gsl_ran_gaussian(m_RandomGenerator, std::sqrt(f_pmtnpe[ipmt]));  // get fluctuation in npe
0224       // add detector resolution
0225       float dres = gsl_ran_gaussian(m_RandomGenerator, f_pmtnpe[ipmt] * 0.127);  // to match real data
0226
0227       f_pmtnpe[ipmt] += (dnpe + dres);  // apply the fluctuations
0228       f_pmtq[ipmt] = f_pmtnpe[ipmt] / _gains[ipmt];
0229
0230       // Now time
0231       if (first_time[ipmt] < 9999.)
0232       {
0233         f_pmtt0[ipmt] = first_time[ipmt] - 8.346;
0234         f_pmtt1[ipmt] = first_time[ipmt] - 8.346;
0235       }
0236       else  // should never happen
0237       {
0238         if (Verbosity() != 0)
0239         {
0240           std::cout << "ERROR, have hit but no time" << std::endl;
0241         }
0242       }
0243
0244       _bbcpmts->get_pmt(ipmt)->set_pmt(ipmt, f_pmtq[ipmt], f_pmtt0[ipmt], f_pmtt1[ipmt]);
0245       _bbcpmts->get_pmt(ipmt)->set_simpmt(f_pmtnpe[ipmt]);
0246
0247       if (Verbosity() > 0)
0248       {
0249         std::cout << "Adding " << ipmt << ", " << f_pmtq[ipmt] << ", " << f_pmtt0[ipmt] << " , " << f_pmtt1[ipmt] << std::endl;
0250       }
0251     }
0252     else
0253     {
0254       // empty channel
0255       _bbcpmts->get_pmt(ipmt)->set_pmt(ipmt, f_pmtq[ipmt], f_pmtt0[ipmt], f_pmtt1[ipmt]);
0256       _bbcpmts->get_pmt(ipmt)->set_simpmt(f_pmtnpe[ipmt]);
0257     }
0258   }
0259
0260   return 0;
0261 }
0262
0263 void MbdDigitization::CreateNodes(PHCompositeNode *topNode)
0264 {
0265   PHNodeIterator iter(topNode);
0266   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0267   if (!dstNode)
0268   {
0269     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0270     gSystem->Exit(1);
0271     exit(1);
0272   }
0273
0274   PHCompositeNode *bbcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "MBD"));
0275   if (!bbcNode)
0276   {
0277     bbcNode = new PHCompositeNode("MBD");
0278     dstNode->addNode(bbcNode);
0279   }
0280
0281   //-* contains info for each pmt (nmips, time, etc)
0282   _bbcpmts = findNode::getClass<MbdPmtContainer>(bbcNode, "MbdPmtContainer");
0283   if (!_bbcpmts)
0284   {
0285     _bbcpmts = new MbdPmtSimContainerV1();
0286     PHIODataNode<PHObject> *BbcPmtNode = new PHIODataNode<PHObject>(_bbcpmts, "MbdPmtContainer", "PHObject");
0287     bbcNode->addNode(BbcPmtNode);
0288   }
0289 }
0290
0291 //___________________________________
0292 void MbdDigitization::GetNodes(PHCompositeNode *topNode)
0293 {
0294   // Get the DST objects
0295
0296   // Truth container
0297   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0298   if (!_truth_container)
0299   {
0300     std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << std::endl;
0301     gSystem->Exit(1);
0302   }
0303
0304   // BBC hit container
0305   _bbchits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BBC");
0306   if (!_bbchits)
0307   {
0308     std::cout << PHWHERE << " G4HIT_BBC node not found on node tree" << std::endl;
0309     gSystem->Exit(1);
0310   }
0311
0312   /** DST Objects **/
0313
0314   // MbdPmtContainer
0315   _bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0316   if (!_bbcpmts)
0317   {
0318     std::cout << PHWHERE << " MbdPmtContainer node not found on node tree" << std::endl;
0319     gSystem->Exit(1);
0320   }
0321 }