← Back to branch summary

~mwinter4/maus/ckov_0_9_3

« back to all changes in this revision

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

Merged from trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerMCDigitization/MapCppTrackerMCDigitization.cc
18
18
#include "src/common_cpp/DataStructure/ReconEvent.hh"
19
19
#include "src/common_cpp/JsonCppProcessors/SpillProcessor.hh"
20
20
#include "src/map/MapCppTrackerMCDigitization/MapCppTrackerMCDigitization.hh"
 
21
#include "src/common_cpp/Recon/SciFi/SciFiLookup.hh"
21
22
 
22
23
namespace MAUS {
 
24
 
23
25
MapCppTrackerMCDigitization::MapCppTrackerMCDigitization()
24
26
    : _spill_json(NULL), _spill_cpp(NULL) {
25
27
}
44
46
    modules = mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
45
47
    Json::Value *json = Globals::GetConfigurationCards();
46
48
    // Load constants.
47
 
    _SciFiNPECut        = (*json)["SciFiDigitNPECut"].asDouble();
 
49
    _SciFiNPECut        = (*json)["SciFiNoiseNPECut"].asDouble();
48
50
    _SciFivlpcEnergyRes = (*json)["SciFivlpcEnergyRes"].asDouble();
49
51
    _SciFiadcFactor     = (*json)["SciFiadcFactor"].asDouble();
50
52
    _SciFitdcBits       = (*json)["SciFitdcBits"].asDouble();
64
66
                 _SciFivlpcQE;
65
67
    // ______________________________________________
66
68
    return true;
67
 
  } catch(Exception& exception) {
 
69
  } catch (Exception& exception) {
68
70
    MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(exception, _classname);
69
 
  } catch(std::exception& exc) {
 
71
  } catch (std::exception& exc) {
70
72
    MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(exc, _classname);
71
73
  }
72
74
  return false;
102
104
    spill.SetReconEvents(revts);
103
105
  }
104
106
  // Construct digits from hits for each MC event
105
 
  for ( unsigned int event_i = 0; event_i < spill.GetMCEventSize(); event_i++ ) {
 
107
  for ( size_t event_i = 0; event_i < spill.GetMCEventSize(); event_i++ ) {
106
108
    MCEvent *mc_evt = spill.GetMCEvents()->at(event_i);
107
109
    SciFiDigitPArray digits;
108
110
    if ( mc_evt->GetSciFiHits() ) {
109
 
      construct_digits(mc_evt->GetSciFiHits(), spill.GetSpillNumber(),
110
 
                       static_cast<int>(event_i), digits);
 
111
      construct_digits(mc_evt->GetSciFiHits(), spill.GetSpillNumber(), event_i, digits);
 
112
    }
 
113
    if ( mc_evt->GetSciFiNoiseHits() ) {
 
114
      add_noise(mc_evt->GetSciFiNoiseHits(), digits);
111
115
    }
112
116
 
113
117
    // Make a SciFiEvent to hold the digits produced
142
146
    json_root = JsonWrapper::StringToJson(json_data);
143
147
    SpillProcessor spill_proc;
144
148
    _spill_cpp = spill_proc.JsonToCpp(json_root);
145
 
  } catch(...) {
 
149
  } catch (...) {
146
150
    Squeak::mout(Squeak::error) << "Bad json document" << std::endl;
147
151
    _spill_cpp = new Spill();
148
152
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
155
159
 
156
160
void MapCppTrackerMCDigitization::construct_digits(SciFiHitArray *hits, int spill_num,
157
161
                                                   int event_num, SciFiDigitPArray &digits) {
158
 
 
 
162
  SciFiLookup lookup;
159
163
  for ( unsigned int hit_i = 0; hit_i < hits->size(); hit_i++ ) {
160
164
    if ( !hits->at(hit_i).GetChannelId()->GetUsed() ) {
 
165
 
161
166
      SciFiHit *a_hit = &hits->at(hit_i);
162
167
      a_hit->GetChannelId()->SetUsed(true);
163
168
 
171
176
      double time1   = a_hit->GetTime();
172
177
      // int tdcCounts = compute_tdc_counts(time1);
173
178
 
 
179
      // Pull tracker/station/plane information from geometry
 
180
      int tracker = a_hit->GetChannelId()->GetTrackerNumber();
 
181
      int station = a_hit->GetChannelId()->GetStationNumber();
 
182
      int plane   = a_hit->GetChannelId()->GetPlaneNumber();
 
183
 
 
184
      // Create Digit for use with MC lookup, will update NPE later
 
185
      SciFiDigit *a_digit = new SciFiDigit(spill_num, event_num,
 
186
                                           tracker, station, plane, chanNo, nPE, time1);
 
187
 
 
188
      a_hit->GetChannelId()->SetID(lookup.get_digit_id(a_digit));
 
189
 
174
190
      // Loop over all the other hits.
175
191
      for ( unsigned int hit_j = hit_i+1; hit_j < hits->size(); hit_j++ ) {
176
192
        if ( check_param(&(*hits)[hit_i], &(*hits)[hit_j]) ) {
178
194
          double edep_j = same_digit->GetEnergyDeposited();
179
195
          nPE += edep_j*_eV_to_phe;
180
196
          same_digit->GetChannelId()->SetUsed(true);
 
197
          same_digit->GetChannelId()->SetID(lookup.get_digit_id(a_digit));
181
198
        }
182
199
      } // ends loop over all the array
183
 
      int tracker = a_hit->GetChannelId()->GetTrackerNumber();
184
 
      int station = a_hit->GetChannelId()->GetStationNumber();
185
 
      int plane   = a_hit->GetChannelId()->GetPlaneNumber();
186
 
 
187
 
      SciFiDigit *a_digit = new SciFiDigit(spill_num, event_num,
188
 
                                           tracker, station, plane, chanNo, nPE, time1);
 
200
 
 
201
      a_digit->set_npe(nPE);
 
202
 
189
203
      // .start. TO BE REMOVED .start.//
190
204
      ThreeVector position = a_hit->GetPosition();
191
205
      ThreeVector momentum = a_hit->GetMomentum();
192
206
      a_digit->set_true_position(position);
193
207
      a_digit->set_true_momentum(momentum);
194
208
      // .end. TO BE REMOVED .end.//
 
209
 
195
210
      digits.push_back(a_digit);
196
211
    }
197
212
  } // ends 'for' loop over hits
198
213
}
199
214
 
 
215
void MapCppTrackerMCDigitization::add_noise(SciFiNoiseHitArray *noises,
 
216
                                            SciFiDigitPArray &digits) {
 
217
 
 
218
    /**************************************************************************
 
219
    *  Function checks which channel has noise against which channel has a
 
220
    *  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.
 
225
    **************************************************************************/
 
226
 
 
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());
 
270
      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
    }
 
274
  }
 
275
}
 
276
 
200
277
int MapCppTrackerMCDigitization::compute_tdc_counts(double time1) {
201
278
  double tmpcounts;
202
279