~mwinter4/maus/ckov-update

« back to all changes in this revision

Viewing changes to src/map/MapCppTOFMCDigitizer/MapCppTOFMCDigitizer.cc

  • Committer: Chris Rogers
  • Date: 2013-02-01 13:08:06 UTC
  • mfrom: (659.1.57 release-candidate)
  • Revision ID: chris.rogers@stfc.ac.uk-20130201130806-kkxs626x4loudqsm
Tags: MAUS-v0.4, MAUS-v0.4.3
MAUS-v0.4.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
82
82
  Json::Value tof_evt;
83
83
 
84
84
  // loop over events
 
85
  std::cout << "mc numevts = " << mc.size() << std::endl;
85
86
  if (fDebug) std::cout << "mc numevts = " << mc.size() << std::endl;
86
87
  for ( unsigned int i = 0; i < mc.size(); i++ ) {
87
88
    Json::Value particle = mc[i];
130
131
 
131
132
  for ( unsigned int j = 0; j < hits.size(); j++ ) {  //  j-th hit
132
133
      Json::Value hit = hits[j];
133
 
      if (fDebug) std::cout << "hit# " << j << hit << std::endl;
 
134
      // if (fDebug) std::cout << "hit# " << j << hit << std::endl;
 
135
      if (fDebug) std::cout << "=================== hit# " << j << std::endl;
134
136
 
135
137
      // make sure we can get the station/slab info
136
138
      if (!hit.isMember("channel_id"))
156
158
      if (fDebug) {
157
159
         std::cout << "tofhit: " << hit["channel_id"] << " "
158
160
                   << hit["position"] << " " << hit["momentum"]
159
 
                   << " " << hit["time"] << std::endl;
 
161
                   << " " << hit["time"] << " " << hit["energy_deposited"] << std::endl;
160
162
      }
161
163
 
162
164
      int stn = hit["channel_id"]["station_number"].asInt();
221
223
      // convert edep to photoelectrons for this slab/pmt
222
224
      // can't convert to adc yet since we need to add up ph.el's
223
225
      //   from other hits if any
 
226
      if (fDebug) std::cout << "edep= " << edep << std::endl;
224
227
      double npe1 = get_npe(dist1, edep);
225
228
      double npe2 = get_npe(dist2, edep);
 
229
      if (fDebug) printf("npe# %3.15f %3.4f %3.4f\n", edep, npe1, npe2);
 
230
      if (fDebug) printf("npe# %3.15f %3.4f %3.4f\n", edep, npe1, npe2);
226
231
 
227
232
      // get the hit time
228
233
      double csp = _configJSON["TOFscintLightSpeed"].asDouble();
334
339
}
335
340
 
336
341
//////////////////////////////////////////////////////////////////////
337
 
double MapCppTOFMCDigitizer::get_npe(double edep, double dist) {
338
 
      double peRes = 1e-4;
 
342
double MapCppTOFMCDigitizer::get_npe(double dist, double edep) {
339
343
      double nphot = 0;
 
344
      double nptmp = 0.;
 
345
      if (fDebug) std::cout << "get_npe::edep= " << edep << std::endl;
340
346
 
341
347
      if (!_configJSON.isMember("TOFattenuationLength"))
342
348
          throw(Squeal(Squeal::recoverable,
350
356
          throw(Squeal(Squeal::recoverable,
351
357
                       "Could not find TOFconversionFactor in config",
352
358
                       "MapCppTOFMCDigitizer::birth"));
 
359
      if (fDebug)
 
360
          printf("nphot0: %3.12f %3.12f\n", edep, _configJSON["TOFconversionFactor"].asDouble());
 
361
      // convert energy deposited to number of photoelectrons
353
362
      nphot = edep / (_configJSON["TOFconversionFactor"].asDouble());
354
363
 
 
364
      TRandom* rnd = new TRandom();
 
365
      rnd = new TRandom();
 
366
      nptmp = rnd->Poisson(nphot);
 
367
      nphot = nptmp;
 
368
      // attenuate the yield
355
369
      nphot *= exp(-dist / (_configJSON["TOFattenuationLength"].asDouble()));
 
370
      rnd = new TRandom();
 
371
      nptmp = rnd->Poisson(nphot);
 
372
      nphot = nptmp;
 
373
      // correct for phototube quantum efficiency
356
374
      nphot *= (_configJSON["TOFpmtQuantumEfficiency"].asDouble());
357
 
      nphot = CLHEP::RandGauss::shoot(nphot, peRes);
358
 
      // cannot convert to adc yet -- need to add up yields from multi hits in bar
359
 
      // take care of adc conversion after weeding out duplicate/multi hits
360
 
      // int adc = (int)( nPE / (_configJSON.["TOFadcConversionFactor"].asDouble()) );
 
375
      rnd = new TRandom();
 
376
      nptmp = rnd->Poisson(nphot);
 
377
      nphot = nptmp;
 
378
 
361
379
      return nphot;
362
380
}
363
381
 
421
439
 
422
440
      // convert light yield to adc & set the charge
423
441
      int adc = static_cast<int>(npe / (_configJSON["TOFadcConversionFactor"].asDouble()));
 
442
      if (fDebug) std::cout << "npe-adc: " << npe << " " << adc << std::endl;
424
443
      // ROGERS = changed from "charge" to "charge_pm" for data integrity
425
444
      digit["charge_pm"] = adc;
426
445
      // NOTE: needs tweaking/verifying -- DR 3/15