32 #include <JANA/JApplication.h>
33 #include <JANA/JFactory.h>
49 gPARMS->SetDefaultParameter(
"BCAL_SiPM_saturation:VERBOSE",
VERBOSE,
"Verbosity level");
78 attenuation_parameters.clear();
79 eventLoop->GetCalib(
"/BCAL/attenuation_parameters", attenuation_parameters);
105 vector<const DBCALShower*> BCALShowers;
106 loop->Get(BCALShowers);
113 vector<const DMCThrown*> MCThrowns;
114 loop->Get(MCThrowns);
115 unsigned int NumThrown=MCThrowns.size();
117 double thetathrown=0;
118 if (NumThrown == 1) {
119 const DMCThrown* locMCThrown = MCThrowns[0];
120 Ethrown = locMCThrown->
energy();
123 thetathrown = 180./3.14159*pgen.Theta();
126 if (NumThrown==0 || thetathrown > 12) {
129 vector<const DNeutralShower*> NeutralShowers;
130 loop->Get(NeutralShowers);
132 unsigned int NumShowers=NeutralShowers.size();
136 for (
unsigned int i = 0; i < NumShowers; i++){
140 if (locDetector !=
SYS_BCAL)
continue;
142 vector<const DBCALShower*> BCALShowers;
143 locNeutralShower->Get(BCALShowers);
147 Eshower = locBCALShower->
E;
150 if (NumThrown == 1 && Eshower < Ethrown/2.)
167 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Ethrown", Ethrown,
168 "BCAL SiPM Saturation; Thrown Energy (GeV)",4*nbins,0,10);
170 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Eshower", Eshower,
171 "BCAL SiPM Saturation; Shower Energy (GeV)",10*nbins,0,10);
172 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Thrown Theta", thetathrown,
173 "BCAL SiPM Saturation; Thrown Theta (degrees)",4*90,0,90);
176 Fill2DHistogram (
"BCAL_SiPM_saturation",
"Hists2D",
"Eshower_vs_Ethrown", Ethrown, Eshower,
177 "BCAL SiPM Saturation; Thrown Energy (GeV); Shower Energy (GeV)",
178 4*nbins,0,10,4*nbins,0,10);
179 Fill2DHistogram (
"BCAL_SiPM_saturation",
"Hists2D",
"EDiff_vs_Ethrown", Ethrown, Eshower - Ethrown,
180 "BCAL SiPM Saturation; Thrown Energy (GeV); (Shower - Thrown) Energy (GeV)",
181 4*nbins,0,10,nbins,-0.5,0.5);
182 Fill2DHistogram (
"BCAL_SiPM_saturation",
"Hists2D",
"EDiff/Ethrown_vs_Ethrown", Ethrown, Ethrown>0? (Eshower - Ethrown)/Ethrown : 0,
183 "BCAL SiPM Saturation; Thrown Energy (GeV); (EShower - EThrown)/Ethrown",
184 4*nbins,0,10,nbins,-0.2,0.2);
187 vector<const DBCALPoint*> Points;
188 locNeutralShower->Get(Points);
189 uint Ncell = Points.size();
192 "BCAL SiPM Saturation; Number of cells",nbins,0,100);
194 for (
unsigned int j = 0; j < Ncell; j++){
196 float t = locPoint->
t();
197 float z = locPoint->
z();
198 float Ept = locPoint->
E();
199 int module= locPoint->
module();
201 int sector = locPoint->
sector();
205 float z_bcal_center = 212;
207 float zLocal = z - z_bcal_center;
210 float dUp = 0.5 * fibLen + zLocal;
211 float dDown = 0.5 * fibLen - zLocal;
212 if (dUp>fibLen) dUp=fibLen;
214 if (dDown>fibLen) dDown=fibLen;
215 if (dDown<0) dDown=0;
216 int channel = (module-1)*16 + (layer-1)*4 + (sector-1);
219 float attenuation_length = attenuation_parameters[channel][0];
220 float attUp = exp( -dUp / attenuation_length );
221 float attDown = exp( -dDown / attenuation_length );
224 if (
VERBOSE>=3) cout <<
" VERBOSE >=3" <<
" t=" << t <<
" z=" << z << endl;
230 "BCAL SiPM Saturation; Layer Number",5,0,5);
233 "BCAL SiPM Saturation; Point Energy (GeV)",10*nbins,0,10);
241 vector<const DBCALHit*> Hits;
243 uint Nhits = Hits.size();
245 for (
unsigned int j = 0; j < Nhits; j++){
247 float Ehit = locHit->
E;
251 int end = locHit->
end;
262 downpeak = pulse_peak;
270 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit pulse_peak layer=1", pulse_peak,
271 "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
273 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit integral layer=1", Ehit,
274 "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
276 else if (layer == 2) {
277 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit pulse_peak layer=2", pulse_peak,
278 "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
280 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit integral layer=2", Ehit,
281 "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
283 else if (layer == 3) {
284 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit pulse_peak layer=3", pulse_peak,
285 "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
287 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit integral layer=3", Ehit,
288 "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
290 else if (layer == 4) {
291 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit pulse_peak layer=4", pulse_peak,
292 "BCAL SiPM Saturation; Hit Pulse_peak (counts)",nbins,-100,5000);
294 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Hit integral layer=4", Ehit,
295 "BCAL SiPM Saturation; Hit integral (GeV)",nbins,0, 10);
298 cout <<
" ***Illegal layer=" << layer << endl;
307 float E_US = ( upHit / attUp );
308 float E_DS = ( downHit / attDown );
309 float Ecalc = ( E_US + E_DS ) / 2;
316 "BCAL SiPM Saturation; Calc Energy (GeV)",4*nbins,0,10);
317 Fill1DHistogram (
"BCAL_SiPM_saturation",
"Hists1D",
"Ecalc-Ept", Ecalc-Ept,
318 "BCAL SiPM Saturation; Calc-Pt Energy (GeV)",nbins,-0.05,0.05);
double energy(void) const
DetectorSystem_t dDetectorSystem
JEventProcessor_BCAL_SiPM_saturation()
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
const DVector3 & momentum(void) const
jerror_t fini(void)
Called after last event of last event source has been processed.
jerror_t init(void)
Called once at program start.
~JEventProcessor_BCAL_SiPM_saturation()