← Back to branch summary

~durga/maus/rel709

« back to all changes in this revision

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

  • Committer: Durga Rajaram
  • Date: 2013-08-27 04:36:50 UTC
  • mfrom: (659.1.73 rc)
  • Revision ID: durga@fnal.gov-20130827043650-me0hgsbzlzikdoik
Tags: MAUS-v0.7.0
MAUS-v0.7.0

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
}