~maus-scifi/maus/tracker_devel

« back to all changes in this revision

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

  • Committer: cheid001
  • Date: 2015-04-27 21:34:57 UTC
  • Revision ID: cheid001-20150427213457-cqvn7u73q1715sfm
style updates and test corrections, still doesn't pass all tests

Show diffs side-by-side

added added

removed removed

Lines of Context:
48
48
  _configJSON = Globals::GetConfigurationCards();
49
49
  _vlpc_sim_on = (*_configJSON)["SciFi_BG_VLPCOn"].asInt();
50
50
  _rf_sim_on   = (*_configJSON)["SciFi_BG_RFOn"].asInt();
51
 
  
 
51
 
52
52
  if (_vlpc_sim_on == 1) {
53
53
    _vlpc_rate   = (*_configJSON)["SciFi_BG_VLPCRate"].asDouble();
54
54
    vlpc_mean_npe = -log(1-_vlpc_rate);
55
55
  }
56
 
  
 
56
 
57
57
  if (_rf_sim_on == 1) {
58
58
    _rf_ave_act = (*_configJSON)["SciFi_BG_RFAct"].asDouble();
59
59
    _rf_mag     = (*_configJSON)["SciFi_BG_RFMag"].asDouble();
60
60
    _rf_sdev    = (*_configJSON)["SciFi_BG_RFSDev"].asDouble();
61
61
    /***** Total number of events we would expect at the source *****/
62
 
    flux = _rf_ave_act / (1/pow(1.95,2.0)+1/pow(2.15,2.0)+1/pow(2.4,2.0)
63
 
           +1/pow(2.7,2.0)+1/pow(3.05,2.0))*(5.0);
 
62
    flux = _rf_ave_act / (1/pow(1.95, 2.0)+1/pow(2.15, 2.0)+1/pow(2.4, 2.0)
 
63
           +1/pow(2.7, 2.0)+1/pow(3.05, 2.0))*(5.0);
64
64
  }
65
65
}
66
66
 
80
80
                          "MC event array not initialised, aborting noise",
81
81
                          "MapCppTrackerMCNoise::_process");
82
82
  }
83
 
  
 
83
 
84
84
  /***** Examine each particle event in a spill individually *****/
85
85
  for ( unsigned int event_i = 0; event_i < spill.GetMCEvents()->size(); event_i++ ) {
86
86
    int spill_n = spill.GetSpillNumber();
87
87
    SciFiNoiseHitArray* noise = new SciFiNoiseHitArray();
88
 
  
 
88
 
89
89
    // ===== Noise Implementation ===============================
90
90
 
91
91
    /***** Simulates noise from thermal activation of VLPCs ****/
126
126
  * Individual channels are then selected and assigned a randomized PE value.
127
127
  * Likewise, the NPE is determined then from a Poisson calculation.
128
128
  **************************************************************************************/
129
 
        
 
129
 
130
130
  /***** Examine each tracker plane individually *****/
131
131
  for ( unsigned int mod_i = 0; mod_i < SF_modules.size(); mod_i++ ) {
132
132
    int nChannels = static_cast <int>
133
133
                (2*((SF_modules[mod_i]->propertyDouble("CentralFibre"))+0.5));
134
134
    double mean_channels = static_cast<double>(nChannels*_vlpc_rate);
135
135
    int active_channels = CLHEP::RandPoisson::shoot(mean_channels);
136
 
          
 
136
 
137
137
    /***** Fills array with all channels *****/
138
138
    chanArray.resize(0);
139
139
    for (int ich = 0; ich < nChannels; ich++) {
150
150
      int chan_ran = std::rand() % (nChannels - fch);
151
151
      int chan_i = chanArray[chan_ran];
152
152
      chanArray.erase(chanArray.begin()+chan_ran);
153
 
         
 
153
 
154
154
      /***** Setting PE value *****/
155
155
      int D_NPE = 0;
156
156
      Poisson_Hits(D_NPE, vlpc_mean_npe);
157
 
                
 
157
 
158
158
      /***** Writing out result *****/
159
159
      int tracker = SF_modules[mod_i]->propertyInt("Tracker");
160
160
      int station = SF_modules[mod_i]->propertyInt("Station");
174
174
******************************************************************************/ 
175
175
void MapCppTrackerMCNoise::Poisson_Hits(int& hits, const double& mean) const {
176
176
  double p_one     = mean*exp(-mean);
177
 
  double p_two     = pow(mean,2.0)*exp(-mean)/2;
178
 
  double p_three   = pow(mean,3.0)*exp(-mean)/6;
179
 
  double p_four    = pow(mean,4.0)*exp(-mean)/24;
 
177
  double p_two     = pow(mean, 2.0)*exp(-mean)/2;
 
178
  double p_three   = pow(mean, 3.0)*exp(-mean)/6;
 
179
  double p_four    = pow(mean, 4.0)*exp(-mean)/24;
180
180
  double p_sum     = p_one + p_two + p_three + p_four;
181
181
  double rand_pick = static_cast<double>(std::rand()) / static_cast<double>(RAND_MAX) * p_sum;
182
182
 
187
187
    hits = 2;
188
188
  } else if (rand_pick < p_one+p_two+p_three) {
189
189
    hits = 3;
190
 
  } else hits = 4;
 
190
  } else {
 
191
    hits = 4;
 
192
    }
191
193
}
192
194
 
193
195
void MapCppTrackerMCNoise::RF_X_Ray(Spill& spill, SciFiNoiseHitArray* noise,
211
213
    /***** Mean channels activated determined by distance from RF *****/
212
214
    double mean_channels = 0.0;
213
215
          if (SF_modules[mod_i]->propertyInt("Station") == 1) {
214
 
           mean_channels = flux/pow(1.95,2.0);    
 
216
           mean_channels = flux/pow(1.95, 2.0);
215
217
          } else if (SF_modules[mod_i]->propertyInt("Station") == 2) {
216
 
           mean_channels = flux/pow(2.15,2.0);
 
218
           mean_channels = flux/pow(2.15, 2.0);
217
219
          } else if (SF_modules[mod_i]->propertyInt("Station") == 3) {
218
 
           mean_channels = flux/pow(2.4,2.0);
 
220
           mean_channels = flux/pow(2.4, 2.0);
219
221
          } else if (SF_modules[mod_i]->propertyInt("Station") == 4) {
220
 
           mean_channels = flux/pow(2.7,2.0);
 
222
           mean_channels = flux/pow(2.7, 2.0);
221
223
          } else if (SF_modules[mod_i]->propertyInt("Station") == 5) {
222
 
           mean_channels = flux/pow(3.05,2.0);
223
 
          } 
 
224
           mean_channels = flux/pow(3.05, 2.0);
 
225
          }
224
226
    int active_channels = CLHEP::RandPoisson::shoot(mean_channels);
225
 
          
 
227
 
226
228
    /***** Fills array with all channels *****/
227
229
    chanArray.resize(0);
228
230
    for (int ich = 0; ich < nChannels; ich++) {
239
241
      int chan_ran = std::rand() % (nChannels - fch);
240
242
      int chan_i = chanArray[chan_ran];
241
243
      chanArray.erase(chanArray.begin()+chan_ran);
242
 
                
 
244
 
243
245
      /***** Setting PE value, toy values! *****/
244
246
      double D_NPE = floor(CLHEP::RandGauss::shoot(_rf_mag, _rf_sdev));
245
 
                
 
247
 
246
248
      /***** Creating the noise hit for writing to spill *****/
247
249
      int tracker = SF_modules[mod_i]->propertyInt("Tracker");
248
250
      int station = SF_modules[mod_i]->propertyInt("Station");