← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

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

  • Committer: Durga Rajaram
  • Date: 2014-07-16 15:13:05 UTC
  • mfrom: (659.1.92 cand)
  • Revision ID: durga@fnal.gov-20140716151305-q27rv1y9p03v9lks
Tags: MAUS-v0.9, MAUS-v0.9.0
MAUS-v0.9.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerMCDigitization/MapCppTrackerMCDigitization.cc
14
14
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
15
15
 *
16
16
 */
 
17
#include <map>
17
18
#include "src/common_cpp/DataStructure/MCEvent.hh"
18
19
#include "src/common_cpp/DataStructure/ReconEvent.hh"
19
20
#include "src/common_cpp/JsonCppProcessors/SpillProcessor.hh"
20
21
#include "src/map/MapCppTrackerMCDigitization/MapCppTrackerMCDigitization.hh"
21
22
#include "src/common_cpp/Recon/SciFi/SciFiLookup.hh"
 
23
#include "src/common_cpp/API/PyWrapMapBase.hh"
22
24
 
23
25
namespace MAUS {
 
26
PyMODINIT_FUNC init_MapCppTrackerMCDigitization(void) {
 
27
  PyWrapMapBase<MAUS::MapCppTrackerMCDigitization>::PyWrapMapBaseModInit
 
28
                                ("MapCppTrackerMCDigitization", "", "", "", "");
 
29
}
24
30
 
25
31
MapCppTrackerMCDigitization::MapCppTrackerMCDigitization()
26
 
    : _spill_json(NULL), _spill_cpp(NULL) {
 
32
    : MapBase<Data>("MapCppTrackerMCDigitization") {
27
33
}
28
34
 
29
35
MapCppTrackerMCDigitization::~MapCppTrackerMCDigitization() {
30
 
  if (_spill_json != NULL) {
31
 
    delete _spill_json;
32
 
  }
33
 
  if (_spill_cpp != NULL) {
34
 
    delete _spill_cpp;
35
 
  }
36
 
}
37
 
 
38
 
bool MapCppTrackerMCDigitization::birth(std::string argJsonConfigDocument) {
39
 
  _classname = "MapCppTrackerMCDigitization";
40
 
 
41
 
  try {
42
 
    if (!Globals::HasInstance()) {
43
 
      GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
44
 
    }
45
 
    static MiceModule* mice_modules = Globals::GetMonteCarloMiceModules();
46
 
    modules = mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
47
 
    Json::Value *json = Globals::GetConfigurationCards();
48
 
    // Load constants.
49
 
    _SciFiNPECut        = (*json)["SciFiNoiseNPECut"].asDouble();
50
 
    _SciFivlpcEnergyRes = (*json)["SciFivlpcEnergyRes"].asDouble();
51
 
    _SciFiadcFactor     = (*json)["SciFiadcFactor"].asDouble();
52
 
    _SciFitdcBits       = (*json)["SciFitdcBits"].asDouble();
53
 
    _SciFivlpcTimeRes   = (*json)["SciFivlpcTimeRes"].asDouble();
54
 
    _SciFitdcFactor     = (*json)["SciFitdcFactor"].asDouble();
55
 
    _SciFiFiberConvFactor  = (*json)["SciFiFiberConvFactor"].asDouble();
56
 
    _SciFiFiberTrappingEff = (*json)["SciFiFiberTrappingEff"].asDouble();
57
 
    _SciFiFiberMirrorEff   = (*json)["SciFiFiberMirrorEff"].asDouble();
58
 
    _SciFivlpcQE           = (*json)["SciFivlpcQE"].asDouble();
59
 
    _SciFiFiberTransmissionEff = (*json)["SciFiFiberTransmissionEff"].asDouble();
60
 
    _SciFiMUXTransmissionEff   = (*json)["SciFiMUXTransmissionEff"].asDouble();
61
 
    _eV_to_phe = _SciFiFiberConvFactor *
62
 
                 _SciFiFiberTrappingEff *
63
 
                 ( 1.0 + _SciFiFiberMirrorEff ) *
64
 
                 _SciFiFiberTransmissionEff *
65
 
                 _SciFiMUXTransmissionEff *
66
 
                 _SciFivlpcQE;
67
 
    // ______________________________________________
68
 
    return true;
69
 
  } catch (Exception& exception) {
70
 
    MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(exception, _classname);
71
 
  } catch (std::exception& exc) {
72
 
    MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(exc, _classname);
73
 
  }
74
 
  return false;
75
 
}
76
 
 
77
 
bool MapCppTrackerMCDigitization::death() {
78
 
  return true;
79
 
}
80
 
 
81
 
std::string MapCppTrackerMCDigitization::process(std::string document) {
82
 
  Json::FastWriter writer;
83
 
 
84
 
  // Set up a spill object, then continue only if MC event array is initialised
85
 
  read_in_json(document);
86
 
  Spill& spill = *_spill_cpp;
87
 
 
88
 
  if ( spill.GetMCEvents() ) {
89
 
  } else {
90
 
    std::cerr << "MC event array not initialised, aborting digitisation for this spill\n";
91
 
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
92
 
    std::stringstream ss;
93
 
    ss << _classname << " says:" << "MC event array not initialised, aborting digitisation";
94
 
    errors["missing_branch"] = ss.str();
95
 
    save_to_json(spill);
96
 
    return writer.write(*_spill_json);
 
36
}
 
37
 
 
38
void MapCppTrackerMCDigitization::_birth(const std::string& argJsonConfigDocument) {
 
39
  if (!Globals::HasInstance()) {
 
40
    GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
 
41
  }
 
42
  static MiceModule* mice_modules = Globals::GetMonteCarloMiceModules();
 
43
  modules = mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
 
44
  Json::Value *json = Globals::GetConfigurationCards();
 
45
  // Load constants.
 
46
  _SciFiNPECut        = (*json)["SciFiNoiseNPECut"].asDouble();
 
47
  _SciFivlpcEnergyRes = (*json)["SciFivlpcEnergyRes"].asDouble();
 
48
  _SciFiadcFactor     = (*json)["SciFiadcFactor"].asDouble();
 
49
  _SciFitdcBits       = (*json)["SciFitdcBits"].asDouble();
 
50
  _SciFivlpcTimeRes   = (*json)["SciFivlpcTimeRes"].asDouble();
 
51
  _SciFitdcFactor     = (*json)["SciFitdcFactor"].asDouble();
 
52
  _SciFiFiberConvFactor  = (*json)["SciFiFiberConvFactor"].asDouble();
 
53
  _SciFiFiberTrappingEff = (*json)["SciFiFiberTrappingEff"].asDouble();
 
54
  _SciFiFiberMirrorEff   = (*json)["SciFiFiberMirrorEff"].asDouble();
 
55
  _SciFivlpcQE           = (*json)["SciFivlpcQE"].asDouble();
 
56
  _SciFiFiberTransmissionEff = (*json)["SciFiFiberTransmissionEff"].asDouble();
 
57
  _SciFiMUXTransmissionEff   = (*json)["SciFiMUXTransmissionEff"].asDouble();
 
58
  _eV_to_phe = _SciFiFiberConvFactor *
 
59
               _SciFiFiberTrappingEff *
 
60
               ( 1.0 + _SciFiFiberMirrorEff ) *
 
61
               _SciFiFiberTransmissionEff *
 
62
               _SciFiMUXTransmissionEff *
 
63
               _SciFivlpcQE;
 
64
}
 
65
 
 
66
void MapCppTrackerMCDigitization::_death() {
 
67
}
 
68
 
 
69
void MapCppTrackerMCDigitization::_process(MAUS::Data* data) const {
 
70
  Spill& spill = *(data->GetSpill());
 
71
 
 
72
  if (!spill.GetMCEvents()) {
 
73
    throw MAUS::Exception(Exception::recoverable,
 
74
            "MC event array not initialised.",
 
75
            "MapCppTrackerMCDigitization::process");
97
76
  }
98
77
 
99
78
  // ================= Reconstruction =========================
110
89
    if ( mc_evt->GetSciFiHits() ) {
111
90
      construct_digits(mc_evt->GetSciFiHits(), spill.GetSpillNumber(), event_i, digits);
112
91
    }
 
92
    // Adding in effects from noise in VLPC
113
93
    if ( mc_evt->GetSciFiNoiseHits() ) {
114
94
      add_noise(mc_evt->GetSciFiNoiseHits(), digits);
115
95
    }
116
 
 
 
96
    // Smearing NPE results from ADC resolution
 
97
    // for (size_t digit_j = 0; digit_j < digits.size(); digit_j++ ) {
 
98
    //  digits.at(digit_j)->set_npe(compute_adc_counts(digits.at(digit_j)->get_npe()));
 
99
    // }
117
100
    // Make a SciFiEvent to hold the digits produced
118
101
    SciFiEvent *sf_evt = new SciFiEvent();
119
102
    sf_evt->set_digits(digits);
129
112
      spill.GetReconEvents()->push_back(revt);
130
113
    }
131
114
  }
132
 
  // ==========================================================
133
 
  save_to_json(spill);
134
 
  return writer.write(*_spill_json);
135
 
}
136
 
 
137
 
void MapCppTrackerMCDigitization::read_in_json(std::string json_data) {
138
 
  Json::FastWriter writer;
139
 
  Json::Value json_root;
140
 
  if (_spill_cpp != NULL) {
141
 
    delete _spill_cpp;
142
 
    _spill_cpp = NULL;
143
 
  }
144
 
 
145
 
  try {
146
 
    json_root = JsonWrapper::StringToJson(json_data);
147
 
    SpillProcessor spill_proc;
148
 
    _spill_cpp = spill_proc.JsonToCpp(json_root);
149
 
  } catch (...) {
150
 
    Squeak::mout(Squeak::error) << "Bad json document" << std::endl;
151
 
    _spill_cpp = new Spill();
152
 
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
153
 
    std::stringstream ss;
154
 
    ss << _classname << " says:" << reader.getFormatedErrorMessages();
155
 
    errors["bad_json_document"] = ss.str();
156
 
    _spill_cpp->GetErrors();
157
 
  }
158
115
}
159
116
 
160
117
void MapCppTrackerMCDigitization::construct_digits(SciFiHitArray *hits, int spill_num,
161
 
                                                   int event_num, SciFiDigitPArray &digits) {
 
118
                                                   int event_num, SciFiDigitPArray &digits) const {
162
119
  SciFiLookup lookup;
163
120
  for ( unsigned int hit_i = 0; hit_i < hits->size(); hit_i++ ) {
164
121
    if ( !hits->at(hit_i).GetChannelId()->GetUsed() ) {
199
156
      } // ends loop over all the array
200
157
 
201
158
      a_digit->set_npe(nPE);
202
 
 
203
 
      // .start. TO BE REMOVED .start.//
204
 
      ThreeVector position = a_hit->GetPosition();
205
 
      ThreeVector momentum = a_hit->GetMomentum();
206
 
      a_digit->set_true_position(position);
207
 
      a_digit->set_true_momentum(momentum);
208
 
      // .end. TO BE REMOVED .end.//
209
 
 
210
159
      digits.push_back(a_digit);
211
160
    }
212
161
  } // ends 'for' loop over hits
213
162
}
214
163
 
215
164
void MapCppTrackerMCDigitization::add_noise(SciFiNoiseHitArray *noises,
216
 
                                            SciFiDigitPArray &digits) {
 
165
                                            SciFiDigitPArray &digits) const {
217
166
 
218
167
    /**************************************************************************
219
168
    *  Function checks which channel has noise against which channel has a
220
169
    *  digit.  If noise is in the same channel as a digit then the noise is 
221
 
    *  added to the digit.  If the noise is removed from the digit by one 
222
 
    *  channel, then a new digit is created.
223
 
    *  Regardless, if noise is over the NPE cut off, 2 NPE, then a new digit
224
 
    *  is created.
 
170
    *  added to the digit, else noise is added as a new digit.
225
171
    **************************************************************************/
226
172
 
227
 
  SciFiLookup lookup;
228
 
 
229
 
  for (unsigned int j = 0; j < noises->size(); j++) {
230
 
    for (unsigned int i = 0; i < digits.size(); i++) {
231
 
 
232
 
      // Checks if noise is in the same channel as a digit
233
 
      if (digits[i]->get_tracker() == noises->at(j).GetTracker() &&
234
 
          digits[i]->get_station() == noises->at(j).GetStation() &&
235
 
          digits[i]->get_plane()   == noises->at(j).GetPlane()   &&
236
 
          digits[i]->get_channel() == noises->at(j).GetChannel() &&
237
 
          noises->at(j).GetUsed() == false) {
238
 
        double digit_npe = digits[i]->get_npe();
239
 
        double noise_npe = noises->at(j).GetNPE();
240
 
        digits[i]->set_npe(digit_npe + noise_npe);
241
 
        noises->at(j).SetUsed(true);
242
 
        noises->at(j).SetID(static_cast<double>(lookup.get_digit_id(digits[i])));
243
 
        continue;
244
 
 
245
 
        // Checks if noise is one channel removed from a digit
246
 
      } else if (digits[i]->get_tracker() == noises->at(j).GetTracker() &&
247
 
                 digits[i]->get_station() == noises->at(j).GetStation() &&
248
 
                 digits[i]->get_plane()   == noises->at(j).GetPlane()   &&
249
 
                 abs(digits[i]->get_channel() - noises->at(j).GetChannel()) <= 1.0 &&
250
 
                 noises->at(j).GetUsed() == false) {
251
 
        SciFiNoiseHit a_noise = noises->at(j);
252
 
        SciFiDigit* a_digit = new SciFiDigit(a_noise.GetSpill(), a_noise.GetEvent(),
253
 
                           a_noise.GetTracker(), a_noise.GetStation(), a_noise.GetPlane(),
254
 
                           a_noise.GetChannel(), a_noise.GetNPE(), a_noise.GetTime());
255
 
        digits.push_back(a_digit);
256
 
        noises->at(j).SetUsed(true);
257
 
        noises->at(j).SetID(static_cast<double>(lookup.get_digit_id(a_digit)));
258
 
        continue;
259
 
      }
260
 
    }
261
 
 
262
 
    // Checks if noise is above NPE cutoff
263
 
    if (noises->at(j).GetNPE() >= _SciFiNPECut &&
264
 
        noises->at(j).GetUsed() == false) {
265
 
      SciFiNoiseHit a_noise = noises->at(j);
266
 
      SciFiDigit* a_digit = new SciFiDigit(a_noise.GetSpill(), a_noise.GetEvent(),
267
 
                                           a_noise.GetTracker(), a_noise.GetStation(),
268
 
                                           a_noise.GetPlane(), a_noise.GetChannel(),
269
 
                                           a_noise.GetNPE(), a_noise.GetTime());
 
173
  std::map<int, int> digit_map;
 
174
  for (size_t digit_i = 0; digit_i < digits.size(); digit_i++) {
 
175
    int track_id = digits.at(digit_i)->get_tracker();
 
176
    int stat_id  = digits.at(digit_i)->get_station();
 
177
    int plane_id = digits.at(digit_i)->get_plane();
 
178
    int chan_id  = digits.at(digit_i)->get_channel();
 
179
    int key = chan_id + 1000*plane_id + 10000*stat_id + 100000*track_id;
 
180
    digit_map[key] = digit_i;
 
181
  }
 
182
  for (size_t noise_j = 0; noise_j < noises->size(); noise_j++) {
 
183
    int track_id = noises->at(noise_j).GetTracker();
 
184
    int stat_id  = noises->at(noise_j).GetStation();
 
185
    int plane_id = noises->at(noise_j).GetPlane();
 
186
    int chan_id  = noises->at(noise_j).GetChannel();
 
187
    int key = chan_id + 1000*plane_id + 10000*stat_id + 100000*track_id;
 
188
 
 
189
    if (digit_map.find(key) != digit_map.end()) {
 
190
      double digit_npe = digits.at(digit_map[key])->get_npe();
 
191
      double noise_npe = noises->at(noise_j).GetNPE();
 
192
      digits.at(digit_map[key])->set_npe(digit_npe + noise_npe);
 
193
    } else {
 
194
      SciFiDigit* a_digit = new SciFiDigit(noises->at(noise_j).GetSpill(),
 
195
                                           noises->at(noise_j).GetEvent(),
 
196
                                           noises->at(noise_j).GetTracker(),
 
197
                                           noises->at(noise_j).GetStation(),
 
198
                                           noises->at(noise_j).GetPlane(),
 
199
                                           noises->at(noise_j).GetChannel(),
 
200
                                           noises->at(noise_j).GetNPE(),
 
201
                                           noises->at(noise_j).GetTime());
270
202
      digits.push_back(a_digit);
271
 
      noises->at(j).SetUsed(true);
272
 
      noises->at(j).SetID(static_cast<double>(lookup.get_digit_id(a_digit)));
273
203
    }
274
204
  }
275
205
}
276
206
 
277
 
int MapCppTrackerMCDigitization::compute_tdc_counts(double time1) {
 
207
int MapCppTrackerMCDigitization::compute_tdc_counts(double time1) const {
278
208
  double tmpcounts;
279
209
 
280
210
  tmpcounts = CLHEP::RandGauss::shoot(time1, _SciFivlpcTimeRes)*_SciFitdcFactor;
288
218
  return tdcCounts;
289
219
}
290
220
 
291
 
int MapCppTrackerMCDigitization::compute_chan_no(MAUS::SciFiHit *ahit) {
 
221
int MapCppTrackerMCDigitization::compute_chan_no(MAUS::SciFiHit *ahit) const {
292
222
  // This is the channel number computed from the fibre number.
293
223
  int fiberNumber = ahit->GetChannelId()->GetFibreNumber();
294
224
  int chanNo      = static_cast<int> (floor(fiberNumber/7.0));
296
226
  return chanNo;
297
227
}
298
228
 
299
 
int MapCppTrackerMCDigitization::compute_adc_counts(double numb_pe) {
 
229
int MapCppTrackerMCDigitization::compute_adc_counts(double numb_pe) const {
300
230
  double tmpcounts;
301
231
  if ( numb_pe == 0 ) return 0;
302
232
 
314
244
  return adcCounts;
315
245
}
316
246
 
317
 
bool MapCppTrackerMCDigitization::check_param(MAUS::SciFiHit *hit1, MAUS::SciFiHit *hit2) {
 
247
bool MapCppTrackerMCDigitization::check_param(MAUS::SciFiHit *hit1, MAUS::SciFiHit *hit2) const {
318
248
  if ( hit2->GetChannelId()->GetUsed() ) {
319
249
    return false;
320
250
  } else {
337
267
    }
338
268
  }
339
269
}
340
 
 
341
 
void MapCppTrackerMCDigitization::save_to_json(Spill &spill) {
342
 
    SpillProcessor spill_proc;
343
 
    if (_spill_json != NULL) {
344
 
        delete _spill_json;
345
 
        _spill_json = NULL;
346
 
    }
347
 
    _spill_json = spill_proc.CppToJson(spill, "");
348
 
}
349
 
 
350
270
} // ~namespace MAUS
 
271