g4simulation/g4bbc/MbdDigitization.cc

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