← Back to branch summary

~chris-rogers/maus/emr_mc_digitization

« back to all changes in this revision

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

  • Committer: Chris Rogers
  • Date: 2014-04-16 11:48:45 UTC
  • mfrom: (707 merge)
  • mto: This revision was merged to the branch mainline in revision 711.
  • Revision ID: chris.rogers@stfc.ac.uk-20140416114845-h3u3q7pdcxkxvovs
Update to trunk

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/map/MapCppTrackerMCNoise/MapCppTrackerMCNoise.cc
1
 
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
2
 
 *
3
 
 * MAUS is free software: you can redistribute it and/or modify
4
 
 * it under the terms of the GNU General Public License as published by
5
 
 * the Free Software Foundation, either version 3 of the License, or
6
 
 * (at your option) any later version.
7
 
 *
8
 
 * MAUS is distributed in the hope that it will be useful,
9
 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10
 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11
 
 * GNU General Public License for more details.
12
 
 *
13
 
 * You should have received a copy of the GNU General Public License
14
 
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
15
 
 *
16
 
 */
17
 
 
18
 
#include "src/map/MapCppTrackerMCNoise/MapCppTrackerMCNoise.hh"
19
 
 
20
 
namespace MAUS {
21
 
MapCppTrackerMCNoise::MapCppTrackerMCNoise()
22
 
    : _spill_json(NULL), _spill_cpp(NULL) {
23
 
}
24
 
 
25
 
MapCppTrackerMCNoise::~MapCppTrackerMCNoise() {
26
 
    if (_spill_json != NULL) {
27
 
        delete _spill_json;
28
 
    }
29
 
    if (_spill_cpp != NULL) {
30
 
        delete _spill_cpp;
31
 
    }
32
 
}
33
 
 
34
 
bool MapCppTrackerMCNoise::birth(std::string argJsonConfigDocument) {
35
 
  _classname = "MapCppTrackerMCNoise";
36
 
 
37
 
  try {
38
 
    if (!Globals::HasInstance()) {
39
 
      GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
40
 
    }
41
 
 
42
 
    static MiceModule* _mice_modules = Globals::GetMonteCarloMiceModules();
43
 
    SF_modules = _mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
44
 
    _configJSON = Globals::GetConfigurationCards();
45
 
    poisson_mean = -log(1.0-(*_configJSON)["SciFiDarkCountProababilty"].asDouble());
46
 
    return true;
47
 
  } catch(Squeal& squee) {
48
 
    MAUS::CppErrorHandler::getInstance()->HandleSquealNoJson(squee, _classname);
49
 
  } catch(std::exception& exc) {
50
 
    MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(exc, _classname);
51
 
  }
52
 
  return false;
53
 
}
54
 
 
55
 
bool MapCppTrackerMCNoise::death() {
56
 
  return true;
57
 
}
58
 
 
59
 
std::string MapCppTrackerMCNoise::process(std::string document) {
60
 
  Json::FastWriter writer;
61
 
 
62
 
  // Set up a spill object, then continue only if MC event array is initialised
63
 
  read_in_json(document);
64
 
  Spill& spill = *_spill_cpp;
65
 
 
66
 
  if ( spill.GetReconEvents() ) {
67
 
  } else {
68
 
    std::cerr << "Recon event array not initialised, aborting noise for this spill\n";
69
 
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
70
 
    std::stringstream ss;
71
 
    ss << _classname << " says:" << "Recon event array not initialised, aborting noise";
72
 
    errors["missing_branch"] = ss.str();
73
 
    save_to_json(spill);
74
 
    return writer.write(*_spill_json);
75
 
  }
76
 
 
77
 
  // ================= Noise ==================================
78
 
 
79
 
  // Adds Effects of Noise from Electrons to MC
80
 
  dark_count(spill);
81
 
  // ==========================================================
82
 
 
83
 
  save_to_json(spill);
84
 
  return writer.write(*_spill_json);
85
 
}
86
 
 
87
 
void MapCppTrackerMCNoise::read_in_json(std::string json_data) {
88
 
  Json::FastWriter writer;
89
 
  Json::Value json_root;
90
 
  if (_spill_cpp != NULL) {
91
 
    delete _spill_cpp;
92
 
    _spill_cpp = NULL;
93
 
  }
94
 
 
95
 
  try {
96
 
    json_root = JsonWrapper::StringToJson(json_data);
97
 
    SpillProcessor spill_proc;
98
 
    _spill_cpp = spill_proc.JsonToCpp(json_root);
99
 
  } catch(...) {
100
 
    Squeak::mout(Squeak::error) << "Bad json document" << std::endl;
101
 
    _spill_cpp = new Spill();
102
 
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
103
 
    std::stringstream ss;
104
 
    ss << _classname << " says:" << reader.getFormatedErrorMessages();
105
 
    errors["bad_json_document"] = ss.str();
106
 
    _spill_cpp->GetErrors();
107
 
  }
108
 
}
109
 
 
110
 
void MapCppTrackerMCNoise::save_to_json(Spill &spill) {
111
 
    SpillProcessor spill_proc;
112
 
    if (_spill_json != NULL) {
113
 
        delete _spill_json;
114
 
        _spill_json = NULL;
115
 
    }
116
 
    _spill_json = spill_proc.CppToJson(spill, "");
117
 
}
118
 
 
119
 
void MapCppTrackerMCNoise::dark_count(Spill &spill) {
120
 
 
121
 
  int spill_n = spill.GetSpillNumber();
122
 
  int time1 = 0;  // Fix: Need to assign a value to this!
123
 
  int exist_flag, D_NPE;
124
 
 
125
 
  /*************************************************************************************
126
 
  * Description:
127
 
  * This block of code examines each SciFi channel and determines how many dark count
128
 
  * PE are present as described by a Poisson distribution (as determined by a CLHEP
129
 
  * library).  The Poisson mean is determined when the map is first called by the by a
130
 
  * a simple calculation in the birth function
131
 
  **************************************************************************************/
132
 
 
133
 
  for ( unsigned int event_i = 0; event_i < spill.GetReconEvents()->size(); event_i++ ) {
134
 
    SciFiDigitPArray digits = spill.GetReconEvents()->at(event_i)->GetSciFiEvent()->digits();
135
 
    SciFiDigitPArray temp_digits;
136
 
    for ( unsigned int mod_i = 0; mod_i < SF_modules.size(); mod_i++ ) {
137
 
      int nChannels = static_cast <int>
138
 
                      (2*((SF_modules[mod_i]->propertyDouble("CentralFibre"))+0.5));
139
 
      for ( int chan_i = 0; chan_i < nChannels; chan_i++ ) {
140
 
        exist_flag = 0;
141
 
        D_NPE = CLHEP::RandPoisson::shoot(poisson_mean);
142
 
 
143
 
        /********************************************************************************
144
 
        * Description:
145
 
        * If there is a dark count a new digit is created.  A check is made to determine
146
 
        * if a digit already exist here, if so the npe is just added, otherwise the
147
 
        * noise digit is added to a temp digit array.
148
 
        ********************************************************************************/
149
 
 
150
 
        if ( D_NPE ) {
151
 
          int tracker = SF_modules[mod_i]->propertyInt("Tracker");
152
 
          int station = SF_modules[mod_i]->propertyInt("Station");
153
 
          int plane   = SF_modules[mod_i]->propertyInt("Plane");
154
 
          SciFiDigit *a_digit = new SciFiDigit(spill_n, event_i,
155
 
                                             tracker, station, plane,
156
 
                                             chan_i, D_NPE, time1);
157
 
          for ( unsigned int dig_i = 0; dig_i < digits.size(); dig_i++ ) {
158
 
            if ( digits.at(dig_i)->get_tracker() == a_digit->get_tracker() &&
159
 
                 digits.at(dig_i)->get_station() == a_digit->get_station() &&
160
 
                 digits.at(dig_i)->get_plane()   == a_digit->get_plane()   &&
161
 
                 digits.at(dig_i)->get_channel() == a_digit->get_channel() ) {
162
 
              digits.at(dig_i)->set_npe(digits.at(dig_i)->get_npe()+a_digit->get_npe());
163
 
              exist_flag = 1;
164
 
              continue;
165
 
            }
166
 
          }
167
 
          if ( !exist_flag ) {
168
 
            temp_digits.push_back(a_digit);
169
 
          }
170
 
        }
171
 
      }
172
 
    }
173
 
 
174
 
    /************************************************************************************
175
 
    * Description:
176
 
    * All the digits that do not overlap with a digit from an SciFiHit are now added
177
 
    * into the SciFiDigit array and merged back into the spill.
178
 
    ************************************************************************************/
179
 
 
180
 
    for ( unsigned int temp_i = 0; temp_i < temp_digits.size(); temp_i++ ) {
181
 
      digits.push_back(temp_digits.at(temp_i));
182
 
    }
183
 
    spill.GetReconEvents()->at(event_i)->GetSciFiEvent()->set_digits(digits);
184
 
  } // end of event_i of spill_n, start of event_i+1.
185
 
}
186
 
}
 
1
 
 
2
/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 
3
 *
 
4
 * MAUS is free software: you can redistribute it and/or modify
 
5
 * it under the terms of the GNU General Public License as published by
 
6
 * the Free Software Foundation, either version 3 of the License, or
 
7
 * (at your option) any later version.
 
8
 *
 
9
 * MAUS is distributed in the hope that it will be useful,
 
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
 * GNU General Public License for more details.
 
13
 *
 
14
 * You should have received a copy of the GNU General Public License
 
15
 * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
 
16
 *
 
17
 */
 
18
 
 
19
#include "src/map/MapCppTrackerMCNoise/MapCppTrackerMCNoise.hh"
 
20
 
 
21
namespace MAUS {
 
22
MapCppTrackerMCNoise::MapCppTrackerMCNoise()
 
23
    : _spill_json(NULL), _spill_cpp(NULL) {
 
24
}
 
25
 
 
26
MapCppTrackerMCNoise::~MapCppTrackerMCNoise() {
 
27
    if (_spill_json != NULL) {
 
28
        delete _spill_json;
 
29
    }
 
30
    if (_spill_cpp != NULL) {
 
31
        delete _spill_cpp;
 
32
    }
 
33
}
 
34
 
 
35
bool MapCppTrackerMCNoise::birth(std::string argJsonConfigDocument) {
 
36
  _classname = "MapCppTrackerMCNoise";
 
37
 
 
38
  try {
 
39
    if (!Globals::HasInstance()) {
 
40
      GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
 
41
    }
 
42
 
 
43
    static MiceModule* _mice_modules = Globals::GetMonteCarloMiceModules();
 
44
    SF_modules = _mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi");
 
45
    _configJSON = Globals::GetConfigurationCards();
 
46
    poisson_mean = -log(1.0-(*_configJSON)["SciFiDarkCountProababilty"].asDouble());
 
47
    return true;
 
48
  } catch (Exception& exception) {
 
49
    MAUS::CppErrorHandler::getInstance()->HandleExceptionNoJson(exception, _classname);
 
50
  } catch (std::exception& exc) {
 
51
    MAUS::CppErrorHandler::getInstance()->HandleStdExcNoJson(exc, _classname);
 
52
  }
 
53
  return false;
 
54
}
 
55
 
 
56
bool MapCppTrackerMCNoise::death() {
 
57
  return true;
 
58
}
 
59
 
 
60
std::string MapCppTrackerMCNoise::process(std::string document) {
 
61
  Json::FastWriter writer;
 
62
 
 
63
  // Set up a spill object, then continue only if MC event array is initialised
 
64
  read_in_json(document);
 
65
  Spill& spill = *_spill_cpp;
 
66
 
 
67
  if ( spill.GetMCEvents() ) {
 
68
  } else {
 
69
    std::cerr << "MC event array not initialised, aborting noise for this spill\n";
 
70
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
 
71
    std::stringstream ss;
 
72
    ss << _classname << " says:" << "MC event array not initialised, aborting noise";
 
73
    errors["missing_branch"] = ss.str();
 
74
    save_to_json(spill);
 
75
    return writer.write(*_spill_json);
 
76
  }
 
77
 
 
78
  // ================= Noise ==================================
 
79
 
 
80
  // Adds Effects of Noise from Electrons to MC
 
81
  dark_count(spill);
 
82
  // ==========================================================
 
83
 
 
84
  save_to_json(spill);
 
85
  return writer.write(*_spill_json);
 
86
}
 
87
 
 
88
void MapCppTrackerMCNoise::read_in_json(std::string json_data) {
 
89
  Json::Reader reader;
 
90
  Json::Value json_root;
 
91
  Json::FastWriter writer;
 
92
 
 
93
  if (_spill_cpp != NULL) {
 
94
    delete _spill_cpp;
 
95
    _spill_cpp = NULL;
 
96
  }
 
97
 
 
98
  try {
 
99
    json_root = JsonWrapper::StringToJson(json_data);
 
100
    SpillProcessor spill_proc;
 
101
    _spill_cpp = spill_proc.JsonToCpp(json_root);
 
102
  } catch (...) {
 
103
    std::cerr << "Bad json document" << std::endl;
 
104
    _spill_cpp = new Spill();
 
105
    MAUS::ErrorsMap errors = _spill_cpp->GetErrors();
 
106
    std::stringstream ss;
 
107
    ss << _classname << " says:" << reader.getFormatedErrorMessages();
 
108
    errors["bad_json_document"] = ss.str();
 
109
    _spill_cpp->GetErrors();
 
110
  }
 
111
}
 
112
 
 
113
void MapCppTrackerMCNoise::save_to_json(Spill &spill) {
 
114
    SpillProcessor spill_proc;
 
115
    if (_spill_json != NULL) {
 
116
        delete _spill_json;
 
117
        _spill_json = NULL;
 
118
    }
 
119
    _spill_json = spill_proc.CppToJson(spill, "");
 
120
}
 
121
 
 
122
void MapCppTrackerMCNoise::dark_count(Spill &spill) {
 
123
  int spill_n = spill.GetSpillNumber();
 
124
  double time1 = 0.;  // Fix: Need to assign a value to this!
 
125
  double D_NPE;
 
126
 
 
127
  /*************************************************************************************
 
128
  * Description:
 
129
  * This block of code examines each SciFi channel and determines how many dark count
 
130
  * PE are present as described by a Poisson distribution (as determined by a CLHEP
 
131
  * library).  The Poisson mean is determined when the map is first called by the by a
 
132
  * a simple calculation in the birth function
 
133
  **************************************************************************************/
 
134
 
 
135
  for ( unsigned int event_i = 0; event_i < spill.GetMCEvents()->size(); event_i++ ) {
 
136
    SciFiNoiseHitArray* noise_hits = new SciFiNoiseHitArray();
 
137
    for ( unsigned int mod_i = 0; mod_i < SF_modules.size(); mod_i++ ) {
 
138
      int nChannels = static_cast <int>
 
139
                      (2*((SF_modules[mod_i]->propertyDouble("CentralFibre"))+0.5));
 
140
      for ( int chan_i = 0; chan_i < nChannels; chan_i++ ) {
 
141
        D_NPE = CLHEP::RandPoisson::shoot(poisson_mean);
 
142
 
 
143
        /********************************************************************************
 
144
        * Description:
 
145
        * If there is a dark count a new digit is created.  A check is made to determine
 
146
        * if a digit already exist here, if so the npe is just added, otherwise the
 
147
        * noise digit is added to a temp digit array.
 
148
        ********************************************************************************/
 
149
 
 
150
        if ( D_NPE ) {
 
151
          int tracker = SF_modules[mod_i]->propertyInt("Tracker");
 
152
          int station = SF_modules[mod_i]->propertyInt("Station");
 
153
          int plane   = SF_modules[mod_i]->propertyInt("Plane");
 
154
          SciFiNoiseHit a_noise_hit(spill_n, event_i,
 
155
                                    tracker, station, plane,
 
156
                                    chan_i, D_NPE, time1);
 
157
                  noise_hits->push_back(a_noise_hit);
 
158
        }
 
159
      }
 
160
    }
 
161
 
 
162
    /************************************************************************************
 
163
    * Description:
 
164
    * All the digits that do not overlap with a digit from an SciFiHit are now added
 
165
    * into the SciFiDigit array and merged back into the spill.
 
166
    ************************************************************************************/
 
167
 
 
168
    spill.GetMCEvents()->at(event_i)->SetSciFiNoiseHits(noise_hits);
 
169
  } // end of event_i of spill_n, start of event_i+1.
 
170
}
 
171
 
 
172
} // ~namespace MAUS
 
173