CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalSiPMHitResponse.cc
Go to the documentation of this file.
14 
15 #include "CLHEP/Random/RandPoissonQ.h"
16 
17 #include <math.h>
18 #include <list>
19 
21  const CaloShapes * shapes, bool PreMix1) :
22  CaloHitResponse(parameterMap, shapes), theSiPM(), PreMixDigis(PreMix1),
23  nbins(BUNCHSPACE*HcalPulseShapes::invDeltaTSiPM_), dt(HcalPulseShapes::deltaTSiPM_), invdt(HcalPulseShapes::invDeltaTSiPM_) {}
24 
26 
28  precisionTimedPhotons.clear();
29 }
30 
31 void HcalSiPMHitResponse::finalizeHits(CLHEP::HepRandomEngine* engine) {
32  //do not add PE noise for initial premix
33  if(!PreMixDigis) addPEnoise(engine);
34 
35  photonTimeMap::iterator channelPhotons;
36  for (channelPhotons = precisionTimedPhotons.begin();
37  channelPhotons != precisionTimedPhotons.end();
38  ++channelPhotons) {
39  CaloSamples signal(makeSiPMSignal(channelPhotons->first,
40  channelPhotons->second,
41  engine));
42  bool keep( keepBlank() );
43  if (!keep) {
44  const unsigned int size ( signal.size() ) ;
45  if( 0 != size ) {
46  for( unsigned int i ( 0 ) ; i != size ; ++i ) {
47  keep = keep || signal[i] > 1.e-7 ;
48  }
49  }
50  }
51 
52  LogDebug("HcalSiPMHitResponse") << HcalDetId(signal.id()) << ' ' << signal;
53 
54  if (keep) add(signal);
55  }
56 }
57 
59  DetId id(signal.id());
60  CaloSamples * oldSignal = findSignal(id);
61  if (oldSignal == 0) {
62  theAnalogSignalMap[id] = signal;
63  } else {
64  (*oldSignal) += signal;
65  }
66 }
67 
68 void HcalSiPMHitResponse::add(const PCaloHit& hit, CLHEP::HepRandomEngine* engine) {
69  if (!edm::isNotFinite(hit.time()) &&
70  ((theHitFilter == 0) || (theHitFilter->accepts(hit)))) {
71  HcalDetId id(hit.id());
72  const HcalSimParameters& pars = dynamic_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
73  //divide out mean of crosstalk distribution 1/(1-lambda) = multiply by (1-lambda)
74  double signal(analogSignalAmplitude(id, hit.energy(), pars, engine)*(1-pars.sipmCrossTalk(id)));
75  unsigned int photons(signal + 0.5);
76  double tof( timeOfFlight(id) );
77  double time( hit.time() );
78  if(ignoreTime) time = tof;
79 
80  if (photons > 0)
81  if (precisionTimedPhotons.find(id)==precisionTimedPhotons.end()) {
82  precisionTimedPhotons.insert(
83  std::pair<DetId, photonTimeHist >(id,
84  photonTimeHist(nbins * pars.readoutFrameSize(), 0)
85  )
86  );
87  }
88 
89  LogDebug("HcalSiPMHitResponse") << id;
90  LogDebug("HcalSiPMHitResponse") << " fCtoGeV: " << pars.fCtoGeV(id)
91  << " samplingFactor: " << pars.samplingFactor(id)
92  << " photoelectronsToAnalog: " << pars.photoelectronsToAnalog(id)
93  << " simHitToPhotoelectrons: " << pars.simHitToPhotoelectrons(id);
94  LogDebug("HcalSiPMHitResponse") << " energy: " << hit.energy()
95  << " photons: " << photons
96  << " time: " << time;
97  LogDebug("HcalSiPMHitResponse") << " timePhase: " << pars.timePhase()
98  << " tof: " << tof
99  << " binOfMaximum: " << pars.binOfMaximum()
100  << " phaseShift: " << thePhaseShift_;
101  double tzero(0.0 + pars.timePhase() -
102  (time - tof) -
103  BUNCHSPACE*( pars.binOfMaximum() - thePhaseShift_));
104  LogDebug("HcalSiPMHitResponse") << " tzero: " << tzero;
105  double tzero_bin(-tzero*invdt);
106  LogDebug("HcalSiPMHitResponse") << " corrected tzero: " << tzero_bin << '\n';
107  double t_pe(0.);
108  int t_bin(0);
109  for (unsigned int pe(0); pe<photons; ++pe) {
111  t_bin = int(t_pe*invdt + tzero_bin + 0.5);
112  LogDebug("HcalSiPMHitResponse") << "t_pe: " << t_pe << " t_pe + tzero: " << (t_pe+tzero_bin*dt)
113  << " t_bin: " << t_bin << '\n';
114  if ((t_bin >= 0) &&
115  (static_cast<unsigned int>(t_bin) < precisionTimedPhotons[id].size()))
116  precisionTimedPhotons[id][t_bin] += 1;
117  }
118  }
119 }
120 
121 void HcalSiPMHitResponse::addPEnoise(CLHEP::HepRandomEngine* engine)
122 {
123  // Add SiPM dark current noise to all cells
124  for(std::vector<DetId>::const_iterator idItr = theDetIds->begin();
125  idItr != theDetIds->end(); ++idItr) {
126  HcalDetId id(*idItr);
127  const HcalSimParameters& pars =
128  static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
129 
130  // uA * ns / (fC/pe) = pe!
131  double dc_pe_avg =
132  pars.sipmDarkCurrentuA(id) * dt /
133  pars.photoelectronsToAnalog(id);
134 
135  if (dc_pe_avg <= 0.) continue;
136 
137  int nPreciseBins = nbins * pars.readoutFrameSize();
138 
139  unsigned int sumnoisePE(0);
140  double elapsedTime(0.);
141  for (int tprecise(0); tprecise < nPreciseBins; ++tprecise) {
142  int noisepe = CLHEP::RandPoissonQ::shoot(engine, dc_pe_avg); // add dark current noise
143 
144  if (noisepe > 0) {
145  if (precisionTimedPhotons.find(id)==precisionTimedPhotons.end()) {
146  photonTimeHist photons(nPreciseBins, 0);
147  photons[tprecise] = noisepe;
148  precisionTimedPhotons.insert
149  (std::pair<DetId, photonTimeHist >(id, photons ) );
150  } else {
151  precisionTimedPhotons[id][tprecise] += noisepe;
152  }
153 
154  sumnoisePE += noisepe;
155  }
156  elapsedTime += dt;
157 
158  } // precise time loop
159 
160  LogDebug("HcalSiPMHitResponse") << id;
161  LogDebug("HcalSiPMHitResponse") << " total noise (PEs): " << sumnoisePE;
162 
163  } // detId loop
164 } // HcalSiPMHitResponse::addPEnoise()
165 
168  int preciseSize(parameters.readoutFrameSize() * nbins);
169  CaloSamples result(detId, parameters.readoutFrameSize(), preciseSize);
170  result.setPresamples(parameters.binOfMaximum()-1);
171  result.setPrecise(result.presamples() * nbins, dt);
172  return result;
173 }
174 
176  photonTimeHist const& photonTimeBins,
177  CLHEP::HepRandomEngine* engine) {
178  const HcalSimParameters& pars = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(id));
179  theSiPM.setNCells(pars.pixels(id));
180  theSiPM.setTau(pars.sipmTau());
183 
184  //use to make signal
185  CaloSamples signal( makeBlankSignal(id) );
186  int sampleBin(0), preciseBin(0);
187  signal.resetPrecise();
188  unsigned int pe(0);
189  double hitPixels(0.), elapsedTime(0.);
190  unsigned int sumPE(0);
191  double sumHits(0.);
192 
193  HcalSiPMShape sipmPulseShape(pars.signalShape(id));
194 
195  std::list< std::pair<double, double> > pulses;
196  std::list< std::pair<double, double> >::iterator pulse;
197  double timeDiff, pulseBit;
198  LogDebug("HcalSiPMHitResponse") << "makeSiPMSignal for " << HcalDetId(id);
199 
200  for (unsigned int tbin(0); tbin < photonTimeBins.size(); ++tbin) {
201  pe = photonTimeBins[tbin];
202  sumPE += pe;
203  preciseBin = tbin;
204  sampleBin = preciseBin/nbins;
205  if (pe > 0) {
206  hitPixels = theSiPM.hitCells(engine, pe, 0., elapsedTime);
207  sumHits += hitPixels;
208  LogDebug("HcalSiPMHitResponse") << " elapsedTime: " << elapsedTime
209  << " sampleBin: " << sampleBin
210  << " preciseBin: " << preciseBin
211  << " pe: " << pe
212  << " hitPixels: " << hitPixels ;
213  if (pars.doSiPMSmearing()) {
214  pulses.push_back( std::pair<double, double>(elapsedTime, hitPixels) );
215  } else {
216  signal[sampleBin] += hitPixels;
217  hitPixels *= invdt;
218  signal.preciseAtMod(preciseBin) += 0.6*hitPixels;
219  if (preciseBin > 0)
220  signal.preciseAtMod(preciseBin-1) += 0.2*hitPixels;
221  if (preciseBin < signal.preciseSize() -1)
222  signal.preciseAtMod(preciseBin+1) += 0.2*hitPixels;
223  }
224  }
225 
226  if (pars.doSiPMSmearing()) {
227  pulse = pulses.begin();
228  while (pulse != pulses.end()) {
229  timeDiff = elapsedTime - pulse->first;
230  pulseBit = sipmPulseShape(timeDiff)*pulse->second;
231  LogDebug("HcalSiPMHitResponse") << " pulse t: " << pulse->first
232  << " pulse A: " << pulse->second
233  << " timeDiff: " << timeDiff
234  << " pulseBit: " << pulseBit;
235  signal[sampleBin] += pulseBit;
236  signal.preciseAtMod(preciseBin) += pulseBit*invdt;
237 
238  if (timeDiff > 1 && sipmPulseShape(timeDiff) < 1e-7)
239  pulse = pulses.erase(pulse);
240  else
241  ++pulse;
242  }
243  }
244  elapsedTime += dt;
245  }
246 
247  return signal;
248 }
249 
250 void HcalSiPMHitResponse::setDetIds(const std::vector<DetId> & detIds) {
251  theDetIds = &detIds;
252 }
#define LogDebug(id)
std::vector< float > sipmNonlinearity(const DetId &detId) const
double sipmTau() const
float dt
Definition: AMPTWrapper.h:126
int i
Definition: DBlmapReader.cc:9
double time() const
Definition: PCaloHit.h:36
const std::vector< DetId > * theDetIds
virtual void finalizeHits(CLHEP::HepRandomEngine *) override
Finalize hits.
void setSaturationPars(const std::vector< float > &pars)
Definition: HcalSiPM.cc:189
virtual double photoelectronsToAnalog(const DetId &detId) const
double energy() const
Definition: PCaloHit.h:29
virtual bool keepBlank() const
void resetPrecise()
Definition: CaloSamples.cc:27
int pixels(const DetId &detId) const
int preciseSize() const
get the size
Definition: CaloSamples.h:62
Main class for Parameters in different subdetectors.
virtual void addPEnoise(CLHEP::HepRandomEngine *engine)
tuple result
Definition: mps_fire.py:84
const int keep
bool isNotFinite(T x)
Definition: isFinite.h:10
Creates electronics signals from hits.
virtual double hitCells(CLHEP::HepRandomEngine *, unsigned int pes, double tempDiff=0., double photonTime=0.)
Definition: HcalSiPM.cc:100
unsigned int signalShape(const DetId &detId) const
void setNCells(int nCells)
Definition: HcalSiPM.cc:144
void setTau(double tau)
Definition: HcalSiPM.cc:151
bool doSiPMSmearing() const
virtual const CaloSimParameters & simParameters(const DetId &id) const =0
virtual void add(const PCaloHit &hit, CLHEP::HepRandomEngine *) override
process a single SimHit
void setCrossTalk(double xtalk)
Definition: HcalSiPM.cc:157
virtual void setDetIds(const std::vector< DetId > &detIds)
virtual bool accepts(const PCaloHit &hit) const =0
unsigned int id() const
Definition: PCaloHit.h:43
double timeOfFlight(const DetId &detId) const
static double generatePhotonTime(CLHEP::HepRandomEngine *engine)
int readoutFrameSize() const
for now, the LinearFrames and trhe digis will be one-to-one.
AnalogSignalMap theAnalogSignalMap
Definition: DetId.h:18
double analogSignalAmplitude(const DetId &id, float energy, const CaloSimParameters &parameters, CLHEP::HepRandomEngine *) const
std::vector< unsigned int > photonTimeHist
virtual CaloSamples makeBlankSignal(const DetId &detId) const
double sipmDarkCurrentuA(const DetId &detId) const
const CaloVSimParameterMap * theParameterMap
int size() const
get the size
Definition: CaloSamples.h:24
virtual void initializeHits() override
Initialize hits.
virtual CaloSamples makeSiPMSignal(DetId const &id, photonTimeHist const &photons, CLHEP::HepRandomEngine *)
CaloSamples * findSignal(const DetId &detId)
users can look for the signal for a given cell
float & preciseAtMod(int i)
mutable function to access precise samples
Definition: CaloSamples.h:31
double sipmCrossTalk(const DetId &detId) const
static const double tzero[3]
const CaloVHitFilter * theHitFilter
DetId id() const
get the (generic) id
Definition: CaloSamples.h:21
photonTimeMap precisionTimedPhotons
int binOfMaximum() const
tuple size
Write out results.
HcalSiPMHitResponse(const CaloVSimParameterMap *parameterMap, const CaloShapes *shapes, bool PreMix1=false)