CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HPDNoiseLibraryReader.cc
Go to the documentation of this file.
1 
2 // --------------------------------------------------------
3 // A class to read HPD noise from the library.
4 // The deliverable of the class is the collection of
5 // noisy HcalDetIds with associated noise in units of fC for
6 // 10 time samples. During the library production a higher
7 // theshold is used to find a noisy HPD. A lower threshold is
8 // used to eliminate adding unnecessary quite channels to HPD
9 // noise event collection. Therefore user may not see whole 18
10 // channels for noisy HPD.
11 //
12 // Project: HPD noise library reader
13 // Author: T.Yetkin University of Iowa, Feb. 7, 2008
14 // $Id: HPDNoiseLibraryReader.cc,v 1.4 2011/06/17 12:34:22 abdullin Exp $
15 // --------------------------------------------------------
16 
18 
19 using namespace edm;
20 using namespace std;
21 
23 :theDischargeNoiseRate(0),
24 theIonFeedbackFirstPeakRate(0),
25 theIonFeedbackSecondPeakRate(0),
26 theNoisyPhi(0),
27 theRandFlat(0),
28 theRandGaussQ(0) {
30 
31  ParameterSet pSet = iConfig.getParameter < edm::ParameterSet > ("HPDNoiseLibrary");
32  FileInPath filepath = pSet.getParameter < edm::FileInPath > ("FileName");
33 
34  theHPDName = pSet.getUntrackedParameter < string > ("HPDName", "HPD");
35  string pName = filepath.fullPath();
36 
37  if (pName.find(".") == 0)
38  pName.erase(0, 2);
39  theReader = new HPDNoiseReader(pName);
40  theNames = theReader->allNames(); // all 72x2 HPDs
41 
42  fillRates();
43 }
44 
46  if (theRandFlat)
47  delete theRandFlat;
48 
49  if (theRandGaussQ)
50  delete theRandGaussQ;
51 }
52 
56  }
57 
58  std::string config =
59  "process CorrNoise = {"
60  "service = RandomNumberGeneratorService" "{" "untracked uint32 sourceSeed = 123456789" "}" "}";
61 
62  // create the services
64 
65  // make the services available
67 }
70  if (!rng.isAvailable()) {
71  throw cms::Exception("Configuration") << "HcalHPDNoiseLibrary requires the RandomNumberGeneratorService\n"
72  "which is not present in the configuration file. You must add the service\n"
73  "in the configuration file or remove the modules that require it.";
74  }
75  setRandomEngine(rng->getEngine());
76 }
77 void HPDNoiseLibraryReader::setRandomEngine(CLHEP::HepRandomEngine & engine) {
78  if (theRandGaussQ)
79  delete theRandGaussQ;
80 
81  if (theRandFlat)
82  delete theRandFlat;
83 
84  theRandGaussQ = new CLHEP::RandGaussQ(engine);
85  theRandFlat = new CLHEP::RandFlat(engine);
86 }
87 
89  for (size_t i = 0; i < theNames.size(); ++i) {
91  if (theReader->valid(hpdObj)) {
92  theDischargeNoiseRate.push_back(theReader->dischargeRate(hpdObj));
95  } else {
96  std::cerr << "HPD Handle Object is not valid!" << endl;
97  }
98  }
99 }
100 
102 
103 
104  HPDNoiseData *data = 0; //data->size() is checked wherever actually used
105 
106  // make sure that iphi from HcalDetId is found noisy at first.
107  // In other words, be sure that iphi is in the collection of
108  // noisy Phis
109  if (!(IsNoiseApplicable(iphi)))
110  return data;
111 
112  int zside = 1;
113 
114  if (iphi > 72) {
115  iphi = iphi - 72;
116  zside = -1;
117  }
118  std::string name;
119  if (zside == 1) {
120  name = "ZPlus" + theHPDName + itos(iphi);
121  } else if (zside == -1) {
122  name = "ZMinus" + theHPDName + itos(iphi);
123  } else {
124  cerr << " ZSide Calculation Error." << endl;
125  }
127  if (theReader->valid(hpdObj)) {
128  // randomly select one entry from library for this HPD
129  unsigned long entry = theRandFlat->fireInt(theReader->totalEntries(hpdObj));
130 
131  theReader->getEntry(hpdObj, entry, &data);
132  } else {
133  std::cerr << " HPD Name in the library is not valid." << std::endl;
134  }
135  return data;
136 }
137 
138 
140 
141  clearPhi();
142  double rndm[144];
143 
144  theRandFlat->shootArray(144, rndm);
145 
146  for (int i = 0; i < 144; ++i) {
147  if (rndm[i] < theDischargeNoiseRate[i]) {
148  theNoisyPhi.push_back(i + 1);
149  }
150  }
151 }
152 
154 
155  clearPhi();
156  double rndm[144];
157 
158  theRandFlat->shootArray(144, rndm);
159  for (int i = 0; i < 144; ++i) {
160  if (rndm[i] < theDischargeNoiseRate[i]) {
161  theNoisyPhi.push_back(i + 1);
162  }
163  }
164  // make sure one HPD is always noisy
165  if (theNoisyPhi.size() == 0) {
166  int iPhi = (theRandFlat->fireInt(144)) + 1; // integer from interval [0-144[ + 1
167 
168  theNoisyPhi.push_back(iPhi);
169  }
170 }
171 
172 vector <pair<HcalDetId, const float *> >HPDNoiseLibraryReader::getNoisyHcalDetIds() {
173 
174  vector <pair< HcalDetId, const float *> >result;
175 
176  // decide which phi are noisy by using noise rates
177  // and random numbers (to be called for each event)
178  getNoisyPhis();
179  for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
180  int iphi = theNoisyPhi[i];
182 
183  data = getNoiseData(iphi);
184  for (unsigned int i = 0; i < data->size(); ++i) {
185  pair < HcalDetId, const float *>tmp_pair(data->getDataFrame(i).id(), data->getDataFrame(i).getFrame());
186 
187  result.push_back(tmp_pair);
188  }
189  }
190  return result;
191 }
192 
193 vector <pair<HcalDetId, const float *> >HPDNoiseLibraryReader::getNoisyHcalDetIds(int timeSliceId)
194 {
195  vector <pair< HcalDetId, const float *> >result;
196  // decide which phi are noisy by using noise rates
197  // and random numbers (to be called for each event)
198  getNoisyPhis();
199  for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
200  int iphi = theNoisyPhi[i];
202 
203  data = getNoiseData(iphi);
204  for (unsigned int i = 0; i < data->size(); ++i) {
205  float* data_ = const_cast<float*>(data->getDataFrame(i).getFrame());
206  shuffleData(timeSliceId, data_);
207  const float* _data_ =const_cast<const float*>(data_);
208  pair < HcalDetId, const float *>tmp_pair(data->getDataFrame(i).id(), _data_);
209  result.push_back(tmp_pair);
210  }
211  }
212  return result;
213 
214 }
215 vector < pair < HcalDetId, const float *> >HPDNoiseLibraryReader::getBiasedNoisyHcalDetIds(int timeSliceId) {
216 
217  vector < pair < HcalDetId, const float *> >result;
218 
219  // decide which phi are noisy by using noise rates
220  // and random numbers (to be called for each event)
221  // at least one Phi is always noisy.
223  for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
224  int iphi = theNoisyPhi[i];
226 
227  data = getNoiseData(iphi);
228  for (unsigned int i = 0; i < data->size(); ++i) {
229  float* data_ = const_cast<float*>(data->getDataFrame(i).getFrame());
230  shuffleData(timeSliceId, data_);
231  const float* _data_ =const_cast<const float*>(data_);
232  pair < HcalDetId, const float *>tmp_pair(data->getDataFrame(i).id(), _data_);
233  result.push_back(tmp_pair);
234  }
235  }
236  return result;
237 }
238 
239 vector < pair < HcalDetId, const float *> >HPDNoiseLibraryReader::getBiasedNoisyHcalDetIds() {
240 
241  vector < pair < HcalDetId, const float *> >result;
242 
243  // decide which phi are noisy by using noise rates
244  // and random numbers (to be called for each event)
245  // at least one Phi is always noisy.
247  for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
248  int iphi = theNoisyPhi[i];
250 
251  data = getNoiseData(iphi);
252  for (unsigned int i = 0; i < data->size(); ++i) {
253  pair < HcalDetId, const float *>tmp_pair(data->getDataFrame(i).id(), data->getDataFrame(i).getFrame());
254 
255  result.push_back(tmp_pair);
256  }
257  }
258  return result;
259 }
260 
262 
263  // constants for simulation/parameterization
264  double pe2Charge = 0.333333; // fC/p.e.
265  double GeVperfC = 0.177; // in unit GeV/fC and this will be updated when it start reading from DB.
266  double PedSigma = 0.8;
267  double noise = 0.; // fC
268 
269  int iphi = (id.ieta() > 0) ? (id.iphi()) : (id.iphi() + 72);
270  double rateInTail = theIonFeedbackFirstPeakRate[iphi - 1];
271  double rateInSecondTail = theIonFeedbackSecondPeakRate[iphi - 1];
272 
273  if (bias != 0.) {
274  rateInTail = rateInTail * bias;
275  rateInSecondTail = rateInSecondTail * bias;
276  } else {
277  edm::LogError("HPDNoise") << "HPDNoise: ion feedback error (biased or unbiased selection)." << bias << " failed";
278  throw cms::Exception("Unknown", "biase selection ")
279  << "Usage of " << bias << " fails\n";
280  }
281  double Charge = energy / GeVperfC;
282 
283  // three gauss fit is applied to data to get ion feedback distribution
284  // the script is at neutralino: /home/tyetkin/work/hpd_noise/PlotFromPelin.C
285  // a new fit woth double-sigmoids is under way.
286  // parameters (in fC)
287  // first gaussian
288  // double p0 = 9.53192e+05;
289  // double p1 = -3.13653e-01;
290  // double p2 = 2.78350e+00;
291 
292  // second gaussian
293  // double p3 = 2.41611e+03;
294  double p4 = 2.06117e+01;
295  double p5 = 1.09239e+01;
296 
297  // third gaussian
298  // double p6 = 3.42793e+01;
299  double p7 = 5.45548e+01;
300  double p8 = 1.59696e+01;
301 
302  if (Charge > 3 * PedSigma) { // 3 sigma away from pedestal mean
303  int npe = int (Charge / pe2Charge);
304  double a = 0.;
305  double b = 0.;
306 
307  for (int j = 0; j < npe; ++j) {
308  double probability = theRandFlat->shoot();
309 
310  if (probability < rateInTail) { // total tail
311  if (probability < rateInSecondTail) { // second tail
312  Rannor(a, b);
313  noise += b * p8 + p7;
314  } else {
315  Rannor(a, b);
316  noise += b * p5 + p4;
317  }
318  }
319  }
320  // add pedestal
321  if (noise > 0.)
322  noise += theRandGaussQ->fire(0, 2 * PedSigma);
323  }
324  return (noise * GeVperfC); // returns noise in GeV.
325 
326 }
327 
329 
330  bool isAccepted = false;
331  vector < int >::iterator phi_iter;
332 
333  phi_iter = find(theNoisyPhi.begin(), theNoisyPhi.end(), iphi);
334  if (phi_iter != theNoisyPhi.end()) {
335  isAccepted = true;
336  }
337  return isAccepted;
338 }
339 void HPDNoiseLibraryReader::shuffleData(int timeSliceId, float* &data)
340 {
341  if(timeSliceId == -1 || (timeSliceId>=10)) return;
342  //make a local copy of input data
343  float Data[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
344  for(int i=0;i<10;++i){
345  Data[i] = data[i];
346  }
347  int ts_max = -1;
348  float max = -999.;
349  for(int i=0;i<10;++i){
350  if(Data[i]>max){
351  max = data[i];
352  ts_max = i;
353  }
354  }
355  if((ts_max == -1)){//couldn't find ts_max, return the same value.
356  return;
357  }else{
358  // always shift the noise to the right by putting zeroes to the previous slices.
359  // the noise is pedestal subtracted. 0 value is acceptable.
360  int k = -1;
361  for(int i=0;i<10;++i){
362  data[i] = 0.;
363  int newIdx = timeSliceId+k;
364  float dd = 0.;
365  if(newIdx < 10){
366  data[newIdx] = Data[ts_max+k];
367  dd = Data[ts_max+k];
368  i = newIdx;
369  }
370  data[i] = dd;
371  ++k;
372  }
373 
374  }
375 }
376 
377 //I couldn't find Rannor in CLHEP/Random. For now, use it from ROOT (copy/paste) by little modification.
378 void HPDNoiseLibraryReader::Rannor(double &a, double &b) {
379  double r,
380  x,
381  y,
382  z;
383 
384  y = theRandFlat->shoot();
385  z = theRandFlat->shoot();
386  x = z * 6.28318530717958623;
387  r = TMath::Sqrt(-2 * TMath::Log(y));
388  a = r * TMath::Sin(x);
389  b = r * TMath::Cos(x);
390 }
392  stringstream s;
393 
394  s << i;
395  return s.str();
396 }
397 
399  theNoisyPhi.clear();
400 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
std::vector< float > theIonFeedbackFirstPeakRate
const HPDNoiseDataFrame & getDataFrame(size_t i) const
retrive frame for the given index
Definition: HPDNoiseData.cc:22
static PluginManager & configure(const Config &)
const float * getFrame() const
array of 10 charges corresponding to one channel
float dischargeRate(Handle fHandle) const
discharge rate for the instance
CLHEP::RandGaussQ * theRandGaussQ
HPDNoiseLibraryReader(const edm::ParameterSet &)
std::vector< std::pair< HcalDetId, const float * > > getBiasedNoisyHcalDetIds()
Handle getHandle(const std::string &fName)
get handle to access data for one HPD instance
double getIonFeedbackNoise(HcalDetId id, double energy, double bias)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
double pe2Charge
unsigned long totalEntries(Handle fHandle) const
total number of noise events for the HPD instance in the library
double double double z
std::vector< int > theNoisyPhi
float ionFeedbackSecondPeakRate(Handle fHandle) const
std::vector< float > theIonFeedbackSecondPeakRate
static ServiceToken createServicesFromConfig(std::string const &config)
bool valid(Handle fHandle) const
test if handle is valid
PluginManager::Config config()
Definition: standard.cc:22
const T & max(const T &a, const T &b)
double p4[4]
Definition: TauolaWrapper.h:92
std::vector< std::string > theNames
std::pair< std::string, MonitorElement * > entry
Definition: ME_MAP.h:8
tuple result
Definition: query.py:137
bool isAvailable() const
Definition: Service.h:47
int j
Definition: DBlmapReader.cc:9
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
HcalDetId id() const
detId for the frame
std::vector< std::pair< HcalDetId, const float * > > getNoisyHcalDetIds()
int k[5][pyjets_maxn]
std::vector< std::string > allNames() const
all HPD instances in the library
void Rannor(double &a, double &b)
double b
Definition: hdecay.h:120
std::vector< float > theDischargeNoiseRate
CLHEP::RandFlat * theRandFlat
float ionFeedbackFirstPeakRate(Handle fHandle) const
ionfeedback rate for the instance
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
double a
Definition: hdecay.h:121
void getEntry(Handle fHandle, HPDNoiseData **fData)
retrive one entry from the sequentially
Definition: DDAxes.h:10
HPDNoiseData * getNoiseData(int iphi)
unsigned size() const
number of noise channels in the event
Definition: HPDNoiseData.h:25
void shuffleData(int timeSliceId, float *&data)