Hall C ROOT/C++ Analyzer (hcana)
THcShowerPlane.cxx
Go to the documentation of this file.
1 
8 #include "THcShowerPlane.h"
9 #include "THcHodoscope.h"
10 #include "TClonesArray.h"
11 #include "THcSignalHit.h"
12 #include "THcGlobals.h"
13 #include "THcParmList.h"
14 #include "THcHitList.h"
15 #include "THcShower.h"
16 #include "THcRawShowerHit.h"
17 #include "TClass.h"
18 #include "math.h"
19 #include "THaTrack.h"
20 #include "THaTrackProj.h"
21 #include "THcHallCSpectrometer.h"
22 
23 #include <cstring>
24 #include <cstdio>
25 #include <cstdlib>
26 #include <iostream>
27 
28 #include <fstream>
29 
30 using namespace std;
31 
32 ClassImp(THcShowerPlane)
33 
34 //______________________________________________________________________________
35 THcShowerPlane::THcShowerPlane( const char* name,
36  const char* description,
37  const Int_t layernum,
38  THaDetectorBase* parent )
39  : THaSubDetector(name,description,parent)
40 {
41  // Normal constructor with name and description
42  fPosADCHits = new TClonesArray("THcSignalHit",fNelem);
43  fNegADCHits = new TClonesArray("THcSignalHit",fNelem);
44 
45  frPosAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
46  frPosAdcPedRaw = new TClonesArray("THcSignalHit", 16);
47  frPosAdcThreshold = new TClonesArray("THcSignalHit", 16);
48  frPosAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
49  frPosAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
50  frPosAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
51 
52  frPosAdcPed = new TClonesArray("THcSignalHit", 16);
53  frPosAdcPulseInt = new TClonesArray("THcSignalHit", 16);
54  frPosAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
55  frPosAdcPulseTime = new TClonesArray("THcSignalHit", 16);
56 
57  frNegAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
58  frNegAdcPedRaw = new TClonesArray("THcSignalHit", 16);
59  frNegAdcThreshold = new TClonesArray("THcSignalHit", 16);
60  frNegAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
61  frNegAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
62  frNegAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
63 
64  frNegAdcPed = new TClonesArray("THcSignalHit", 16);
65  frNegAdcPulseInt = new TClonesArray("THcSignalHit", 16);
66  frNegAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
67  frNegAdcPulseTime = new TClonesArray("THcSignalHit", 16);
68 
69  //#if ROOT_VERSION_CODE < ROOT_VERSION(5,32,0)
70  // fPosADCHitsClass = fPosADCHits->GetClass();
71  // fNegADCHitsClass = fNegADCHits->GetClass();
72  //#endif
73 
74  fLayerNum = layernum;
75 }
76 
77 //______________________________________________________________________________
78 THcShowerPlane::~THcShowerPlane()
79 {
80  // Destructor
81  delete fPosADCHits; fPosADCHits = NULL;
82  delete fNegADCHits; fNegADCHits = NULL;
83 
84  delete frPosAdcErrorFlag; frPosAdcErrorFlag = NULL;
85  delete frPosAdcPedRaw; frPosAdcPedRaw = NULL;
86  delete frPosAdcThreshold; frPosAdcThreshold = NULL;
87  delete frPosAdcPulseIntRaw; frPosAdcPulseIntRaw = NULL;
88  delete frPosAdcPulseAmpRaw; frPosAdcPulseAmpRaw = NULL;
89  delete frPosAdcPulseTimeRaw; frPosAdcPulseTimeRaw = NULL;
90 
91  delete frPosAdcPed; frPosAdcPed = NULL;
92  delete frPosAdcPulseInt; frPosAdcPulseInt = NULL;
93  delete frPosAdcPulseAmp; frPosAdcPulseAmp = NULL;
94  delete frPosAdcPulseTime; frPosAdcPulseTime = NULL;
95 
96  delete frNegAdcErrorFlag; frNegAdcErrorFlag = NULL;
97  delete frNegAdcPedRaw; frNegAdcPedRaw = NULL;
98  delete frNegAdcThreshold; frNegAdcThreshold = NULL;
99  delete frNegAdcPulseIntRaw; frNegAdcPulseIntRaw = NULL;
100  delete frNegAdcPulseAmpRaw; frNegAdcPulseAmpRaw = NULL;
101  delete frNegAdcPulseTimeRaw; frNegAdcPulseTimeRaw = NULL;
102 
103  delete frNegAdcPed; frNegAdcPed = NULL;
104  delete frNegAdcPulseInt; frNegAdcPulseInt = NULL;
105  delete frNegAdcPulseAmp; frNegAdcPulseAmp = NULL;
106  delete frNegAdcPulseTime; frNegAdcPulseTime = NULL;
107 
108  delete [] fPosPedSum;
109  delete [] fPosPedSum2;
110  delete [] fPosPedLimit;
111  delete [] fPosPedCount;
112 
113  delete [] fNegPedSum;
114  delete [] fNegPedSum2;
115  delete [] fNegPedLimit;
116  delete [] fNegPedCount;
117 
118  delete [] fPosPed;
119  delete [] fPosSig;
120  delete [] fPosThresh;
121 
122  delete [] fNegPed;
123  delete [] fNegSig;
124  delete [] fNegThresh;
125 
126 // delete [] fEpos;
127  // delete [] fEneg;
128  // delete [] fEmean;
129 }
130 
131 //_____________________________________________________________________________
132 THaAnalysisObject::EStatus THcShowerPlane::Init( const TDatime& date )
133 {
134  // Extra initialization for shower layer: set up DataDest map
135 
136  if( IsZombie())
137  return fStatus = kInitError;
138 
139  // How to get information for parent
140  // if( GetParent() )
141  // fOrigin = GetParent()->GetOrigin();
142  fParent = GetParent();
143 
144  // Get pointer to Cherenkov object. "hgcer" if SHMS, "cer" if other
145  THcHallCSpectrometer *app=dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
146  if (fParent->GetPrefix()[0] == 'P') {
147  fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("hgcer"));
148  } else {
149  fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("cer"));
150  }
151  if (!fCherenkov) {
152  cout << "****** THcShowerPlane::Init Cherenkov not found! ******" << endl;
153  cout << "****** THcShowerPlane::Accumulate will be skipped ******" << endl;
154  }
155 
156  if( !app ||
157  !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector("hod"))) ) {
158  static const char* const here = "ReadDatabase()";
159  Warning(Here(here),"Hodoscope \"%s\" not found. ","hod");
160  }
161 
162  EStatus status;
163  if( (status=THaSubDetector::Init( date )) )
164  return fStatus = status;
165 
166  return fStatus = kOK;
167 
168 }
169 
170 //_____________________________________________________________________________
171 Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
172 {
173 
174  // Retrieve FADC parameters. In principle may want different dynamic
175  // pedestal and integration range for preshower and shower, but for now
176  // use same parameters
177  char prefix[2];
178  prefix[0]=tolower(fParent->GetPrefix()[0]);
179  prefix[1]='\0';
180  fPedSampLow=0;
181  fPedSampHigh=9;
182  fDataSampLow=23;
183  fDataSampHigh=49;
184  fAdcNegThreshold=0.;
185  fAdcPosThreshold=0.;
186  fStatCerMin=1.;
187  fStatSlop=2.;
188  fStatMaxChi2=10.;
189  DBRequest list[]={
190  {"cal_AdcNegThreshold", &fAdcNegThreshold, kDouble, 0, 1},
191  {"cal_AdcPosThreshold", &fAdcPosThreshold, kDouble, 0, 1},
192  {"cal_ped_sample_low", &fPedSampLow, kInt, 0, 1},
193  {"cal_ped_sample_high", &fPedSampHigh, kInt, 0, 1},
194  {"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
195  {"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1},
196  {"cal_debug_adc", &fDebugAdc, kInt, 0, 1},
197  {"stat_cermin", &fStatCerMin, kDouble, 0, 1},
198  {"stat_slop", &fStatSlop, kDouble, 0, 1},
199  {"stat_maxchisq", &fStatMaxChi2, kDouble, 0, 1},
200  {0}
201  };
202 
203  fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file
204 
205  gHcParms->LoadParmValues((DBRequest*)&list, prefix);
206 
207  // Retrieve more parameters we need from parent class
208  //
209  THcShower* parent = static_cast<THcShower*>(fParent);
210 
211  // Find the number of elements
212  fNelem = parent->GetNBlocks(fLayerNum-1);
213 
214  // Origin of the plane:
215  //
216  // X is average of top X coordinates of the top and bottom blocks
217  // shifted by a half of the block thickness;
218  // Y is average of left and right edges;
219  // Z is _front_ position of the plane along the beam.
220 
221  Double_t BlockThick = parent->GetBlockThick(fLayerNum-1);
222 
223  Double_t xOrig = (parent->GetXPos(fLayerNum-1,0) +
224  parent->GetXPos(fLayerNum-1,fNelem-1))/2 +
225  BlockThick/2;
226 
227  Double_t yOrig = (parent->GetYPos(fLayerNum-1,0) +
228  parent->GetYPos(fLayerNum-1,1))/2;
229 
230  Double_t zOrig = parent->GetZPos(fLayerNum-1);
231 
232  fOrigin.SetXYZ(xOrig, yOrig, zOrig);
233 
234  fAdcTdcOffset=parent->GetAdcTdcOffset();
235 
236  // Create arrays to hold results here
237  //
238 
239  // Pedestal limits per channel.
240 
241  fPosPedLimit = new Int_t [fNelem];
242  fNegPedLimit = new Int_t [fNelem];
243 
244  for(Int_t i=0;i<fNelem;i++) {
245  fPosPedLimit[i] = parent->GetPedLimit(i,fLayerNum-1,0);
246  fNegPedLimit[i] = parent->GetPedLimit(i,fLayerNum-1,1);
247  }
248 
249  fMinPeds = parent->GetMinPeds();
250 
251  InitializePedestals();
252 
253  fNumGoodPosAdcHits = vector<Int_t> (fNelem, 0.0);
254  fNumGoodNegAdcHits = vector<Int_t> (fNelem, 0.0);
255 
256  fGoodPosAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0);
257  fGoodNegAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0);
258 
259  fGoodPosAdcPed = vector<Double_t> (fNelem, 0.0);
260  fGoodPosAdcPulseInt = vector<Double_t> (fNelem, 0.0);
261  fGoodPosAdcPulseAmp = vector<Double_t> (fNelem, 0.0);
262  fGoodPosAdcPulseTime = vector<Double_t> (fNelem, 0.0);
263  fGoodPosAdcTdcDiffTime = vector<Double_t> (fNelem, 0.0);
264 
265  fGoodNegAdcPed = vector<Double_t> (fNelem, 0.0);
266  fGoodNegAdcPulseInt = vector<Double_t> (fNelem, 0.0);
267  fGoodNegAdcPulseAmp = vector<Double_t> (fNelem, 0.0);
268  fGoodNegAdcPulseTime = vector<Double_t> (fNelem, 0.0);
269  fGoodNegAdcTdcDiffTime = vector<Double_t> (fNelem, 0.0);
270 
271  fGoodPosAdcMult = vector<Double_t> (fNelem, 0.0);
272  fGoodNegAdcMult = vector<Double_t> (fNelem, 0.0);
273 
274  // Energy depositions per block (not corrected for track coordinate)
275 
276  fEpos = vector<Double_t> (fNelem, 0.0);
277  fEneg = vector<Double_t> (fNelem, 0.0);
278  fEmean = vector<Double_t> (fNelem, 0.0);
279  // fEpos = new Double_t[fNelem];
280  // fEneg = new Double_t[fNelem];
281  // fEmean= new Double_t[fNelem];
282 
283  // Numbers of tracks and hits , for efficiency calculations.
284 
285  fStatNumTrk = vector<Int_t> (fNelem, 0);
286  fStatNumHit = vector<Int_t> (fNelem, 0);
287  fTotStatNumTrk = 0;
288  fTotStatNumHit = 0;
289 
290  // Debug output.
291 
292  if (parent->fdbg_init_cal) {
293  cout << "---------------------------------------------------------------\n";
294  cout << "Debug output from THcShowerPlane::ReadDatabase for "
295  << fParent->GetPrefix() << ":" << endl;
296 
297  cout << " Layer #" << fLayerNum << ", number of elements " << dec << fNelem
298  << endl;
299 
300  cout << " Origin of Layer at X = " << fOrigin.X()
301  << " Y = " << fOrigin.Y() << " Z = " << fOrigin.Z() << endl;
302 
303  cout << " fPosPedLimit:";
304  for(Int_t i=0;i<fNelem;i++) cout << " " << fPosPedLimit[i];
305  cout << endl;
306  cout << " fNegPedLimit:";
307  for(Int_t i=0;i<fNelem;i++) cout << " " << fNegPedLimit[i];
308  cout << endl;
309 
310  cout << " fMinPeds = " << fMinPeds << endl;
311  cout << "---------------------------------------------------------------\n";
312  }
313 
314  return kOK;
315 }
316 
317 //_____________________________________________________________________________
318 Int_t THcShowerPlane::DefineVariables( EMode mode )
319 {
320  // Initialize global variables and lookup table for decoder
321 
322  if( mode == kDefine && fIsSetup ) return kOK;
323  fIsSetup = ( mode == kDefine );
324 
325  // Register variables in global list
326 
327  if (fDebugAdc) {
328  RVarDef vars[] = {
329  {"posAdcPedRaw", "List of positive raw ADC pedestals", "frPosAdcPedRaw.THcSignalHit.GetData()"},
330  {"posAdcPulseIntRaw", "List of positive raw ADC pulse integrals.", "frPosAdcPulseIntRaw.THcSignalHit.GetData()"},
331  {"posAdcPulseAmpRaw", "List of positive raw ADC pulse amplitudes.", "frPosAdcPulseAmpRaw.THcSignalHit.GetData()"},
332  {"posAdcPulseTimeRaw", "List of positive raw ADC pulse times.", "frPosAdcPulseTimeRaw.THcSignalHit.GetData()"},
333 
334  {"posAdcPed", "List of positive ADC pedestals", "frPosAdcPed.THcSignalHit.GetData()"},
335  {"posAdcPulseInt", "List of positive ADC pulse integrals.", "frPosAdcPulseInt.THcSignalHit.GetData()"},
336  {"posAdcPulseAmp", "List of positive ADC pulse amplitudes.", "frPosAdcPulseAmp.THcSignalHit.GetData()"},
337  {"posAdcPulseTime", "List of positive ADC pulse times.", "frPosAdcPulseTime.THcSignalHit.GetData()"},
338 
339  {"negAdcPedRaw", "List of negative raw ADC pedestals", "frNegAdcPedRaw.THcSignalHit.GetData()"},
340  {"negAdcPulseIntRaw", "List of negative raw ADC pulse integrals.", "frNegAdcPulseIntRaw.THcSignalHit.GetData()"},
341  {"negAdcPulseAmpRaw", "List of negative raw ADC pulse amplitudes.", "frNegAdcPulseAmpRaw.THcSignalHit.GetData()"},
342  {"negAdcPulseTimeRaw", "List of negative raw ADC pulse times.", "frNegAdcPulseTimeRaw.THcSignalHit.GetData()"},
343 
344  {"negAdcPed", "List of negative ADC pedestals", "frNegAdcPed.THcSignalHit.GetData()"},
345  {"negAdcPulseInt", "List of negative ADC pulse integrals.", "frNegAdcPulseInt.THcSignalHit.GetData()"},
346  {"negAdcPulseAmp", "List of negative ADC pulse amplitudes.", "frNegAdcPulseAmp.THcSignalHit.GetData()"},
347  {"negAdcPulseTime", "List of negative ADC pulse times.", "frNegAdcPulseTime.THcSignalHit.GetData()"},
348  { 0 }
349  };
350  DefineVarsFromList( vars, mode);
351  } //end debug statement
352 
353  // Register counters for efficiency calculations in gHcParms so that the
354  // variables can be used in end of run reports.
355 
356  gHcParms->Define(Form("%sstat_trksum%d", fParent->GetPrefix(), fLayerNum),
357  Form("Number of tracks in calo. layer %d",fLayerNum), fTotStatNumTrk);
358  gHcParms->Define(Form("%sstat_hitsum%d", fParent->GetPrefix(), fLayerNum),
359  Form("Number of hits in calo. layer %d", fLayerNum), fTotStatNumHit);
360 
361  cout << "THcShowerPlane::DefineVariables: registered counters "
362  << Form("%sstat_trksum%d",fParent->GetPrefix(),fLayerNum) << " and "
363  << Form("%sstat_hitsum%d",fParent->GetPrefix(),fLayerNum) << endl;
364  // getchar();
365 
366  RVarDef vars[] = {
367  {"posAdcErrorFlag", "List of positive raw ADC Error Flags", "frPosAdcErrorFlag.THcSignalHit.GetData()"},
368  {"negAdcErrorFlag", "List of negative raw ADC Error Flags ", "frNegAdcErrorFlag.THcSignalHit.GetData()"},
369 
370  {"posAdcCounter", "List of positive ADC counter numbers.", "frPosAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"}, //PreSh+ raw occupancy
371  {"negAdcCounter", "List of negative ADC counter numbers.", "frNegAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"}, //PreSh- raw occupancy
372 
373  {"totNumPosAdcHits", "Total Number of Positive ADC Hits", "fTotNumPosAdcHits"}, // PreSh+ raw multiplicity
374  {"totNumNegAdcHits", "Total Number of Negative ADC Hits", "fTotNumNegAdcHits"}, // PreSh+ raw multiplicity
375  {"totnumAdcHits", "Total Number of ADC Hits Per PMT", "fTotNumAdcHits"}, // PreSh raw multiplicity
376 
377  {"numGoodPosAdcHits", "Number of Good Positive ADC Hits Per PMT", "fNumGoodPosAdcHits"}, // PreSh occupancy
378  {"numGoodNegAdcHits", "Number of Good Negative ADC Hits Per PMT", "fNumGoodNegAdcHits"}, // PreSh occupancy
379  {"totNumGoodPosAdcHits", "Total Number of Good Positive ADC Hits", "fTotNumGoodPosAdcHits"}, // PreSh multiplicity
380  {"totNumGoodNegAdcHits", "Total Number of Good Negative ADC Hits", "fTotNumGoodNegAdcHits"}, // PreSh multiplicity
381  {"totnumGoodAdcHits", "TotalNumber of Good ADC Hits Per PMT", "fTotNumGoodAdcHits"}, // PreSh multiplicity
382 
383  {"goodPosAdcPulseIntRaw", "Good Positive Raw ADC Pulse Integrals", "fGoodPosAdcPulseIntRaw"},
384  {"goodNegAdcPulseIntRaw", "Good Negative Raw ADC Pulse Integrals", "fGoodNegAdcPulseIntRaw"},
385 
386  {"goodPosAdcPed", "Good Positive ADC pedestals", "fGoodPosAdcPed"},
387  {"goodPosAdcPulseInt", "Good Positive ADC integrals", "fGoodPosAdcPulseInt"},
388  {"goodPosAdcPulseAmp", "Good Positive ADC amplitudes", "fGoodPosAdcPulseAmp"},
389  {"goodPosAdcPulseTime","Good Positive ADC times", "fGoodPosAdcPulseTime"},
390  {"goodPosAdcTdcDiffTime","Good Positive Hodo Start time-ADC times", "fGoodPosAdcTdcDiffTime"},
391 
392  {"goodNegAdcPed", "Good Negative ADC pedestals", "fGoodNegAdcPed"},
393  {"goodNegAdcPulseInt", "Good Negative ADC integrals", "fGoodNegAdcPulseInt"},
394  {"goodNegAdcPulseAmp", "Good Negative ADC amplitudes", "fGoodNegAdcPulseAmp"},
395  {"goodNegAdcPulseTime","Good Negative ADC times", "fGoodNegAdcPulseTime"},
396  {"goodNegAdcTdcDiffTime","Good Negative Hodo Start time-ADC times", "fGoodNegAdcTdcDiffTime"},
397  {"goodPosAdcMult", "Good Positive ADC Multiplicity", "fGoodPosAdcMult"},
398  {"goodNegAdcMult", "Good Negative ADC Multiplicity", "fGoodNegAdcMult"},
399  {"epos", "Energy Depositions from Positive Side PMTs", "fEpos"},
400  {"eneg", "Energy Depositions from Negative Side PMTs", "fEneg"},
401  {"emean", "Mean Energy Depositions", "fEmean"},
402  {"eplane", "Energy Deposition per plane", "fEplane"},
403  {"eplane_pos", "Energy Deposition per plane from pos. PMTs", "fEplane_pos"},
404  {"eplane_neg", "Energy Deposition per plane from neg. PMTs", "fEplane_neg"},
405  { 0 }
406  };
407 
408  return DefineVarsFromList(vars, mode);
409 }
410 
411 //_____________________________________________________________________________
412 void THcShowerPlane::Clear( Option_t* )
413 {
414  // Clears the hit lists
415 
416  fPosADCHits->Clear();
417  fNegADCHits->Clear();
418 
419  frPosAdcErrorFlag->Clear();
420  frPosAdcPedRaw->Clear();
421  frPosAdcThreshold->Clear();
422  frPosAdcPulseIntRaw->Clear();
423  frPosAdcPulseAmpRaw->Clear();
424  frPosAdcPulseTimeRaw->Clear();
425 
426  frPosAdcPed->Clear();
427  frPosAdcPulseInt->Clear();
428  frPosAdcPulseAmp->Clear();
429  frPosAdcPulseTime->Clear();
430 
431  frNegAdcErrorFlag->Clear();
432  frNegAdcPedRaw->Clear();
433  frNegAdcThreshold->Clear();
434  frNegAdcPulseIntRaw->Clear();
435  frNegAdcPulseAmpRaw->Clear();
436  frNegAdcPulseTimeRaw->Clear();
437 
438  frNegAdcPed->Clear();
439  frNegAdcPulseInt->Clear();
440  frNegAdcPulseAmp->Clear();
441  frNegAdcPulseTime->Clear();
442 
443  for (UInt_t ielem = 0; ielem < fGoodPosAdcPed.size(); ielem++) {
444  fGoodPosAdcPed.at(ielem) = 0.0;
445  fGoodPosAdcPulseIntRaw.at(ielem) = 0.0;
446  fGoodPosAdcPulseInt.at(ielem) = 0.0;
447  fGoodPosAdcPulseAmp.at(ielem) = 0.0;
448  fGoodPosAdcPulseTime.at(ielem) = kBig;
449  fGoodPosAdcTdcDiffTime.at(ielem) = kBig;
450  fGoodPosAdcMult.at(ielem) = 0.0;
451  fEpos.at(ielem) = 0.0;
452  fNumGoodPosAdcHits.at(ielem) = 0.0;
453  }
454 
455  for (UInt_t ielem = 0; ielem < fGoodNegAdcPed.size(); ielem++) {
456  fGoodNegAdcPed.at(ielem) = 0.0;
457  fGoodNegAdcPulseIntRaw.at(ielem) = 0.0;
458  fGoodNegAdcPulseInt.at(ielem) = 0.0;
459  fGoodNegAdcPulseAmp.at(ielem) = 0.0;
460  fGoodNegAdcPulseTime.at(ielem) = kBig;
461  fGoodNegAdcTdcDiffTime.at(ielem) = kBig;
462  fGoodNegAdcMult.at(ielem) = 0.0;
463  fEneg.at(ielem) = 0.0;
464  fNumGoodNegAdcHits.at(ielem) = 0.0;
465  }
466 
467  for (UInt_t ielem = 0; ielem < fEmean.size(); ielem++) {
468  fEmean.at(ielem) = 0.0;
469  }
470 
471  fTotNumAdcHits = 0;
472  fTotNumPosAdcHits = 0;
473  fTotNumNegAdcHits = 0;
474 
475  fTotNumGoodAdcHits = 0;
476  fTotNumGoodPosAdcHits = 0;
477  fTotNumGoodNegAdcHits = 0;
478 
479 
480  // Debug output.
481  if ( static_cast<THcShower*>(GetParent())->fdbg_decoded_cal ) {
482  cout << "---------------------------------------------------------------\n";
483  cout << "Debug output from THcShowerPlane::Clear for "
484  << GetParent()->GetPrefix() << ":" << endl;
485 
486  cout << " Cleared ADC hits for plane " << GetName() << endl;
487  cout << "---------------------------------------------------------------\n";
488  }
489 }
490 
491 //_____________________________________________________________________________
492 Int_t THcShowerPlane::Decode( const THaEvData& evdata )
493 {
494  // Doesn't actually get called. Use Fill method instead
495 
496  //Debug output.
497  if ( static_cast<THcShower*>(fParent)->fdbg_decoded_cal ) {
498  cout << "---------------------------------------------------------------\n";
499  cout << "Debug output from THcShowerPlane::Decode for "
500  << fParent->GetPrefix() << ":" << endl;
501 
502  cout << " Called for plane " << GetName() << endl;
503  cout << "---------------------------------------------------------------\n";
504  }
505 
506  return 0;
507 }
508 
509 //_____________________________________________________________________________
510 Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks )
511 {
512 
513  // Nothing is done here. See ProcessHits method instead.
514  //
515 
516  return 0;
517 }
518 
519 //_____________________________________________________________________________
520 Int_t THcShowerPlane::FineProcess( TClonesArray& tracks )
521 {
522  return 0;
523 }
524 
525 //_____________________________________________________________________________
526 Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
527 {
528  // Extract the data for this layer from hit list
529  // Assumes that the hit list is sorted by layer, so we stop when the
530  // plane doesn't agree and return the index for the next hit.
531 
532  // Initialize variables.
533 
534  fPosADCHits->Clear();
535  fNegADCHits->Clear();
536 
537  frPosAdcErrorFlag->Clear();
538  frPosAdcPedRaw->Clear();
539  frPosAdcThreshold->Clear();
540  frPosAdcPulseIntRaw->Clear();
541  frPosAdcPulseAmpRaw->Clear();
542  frPosAdcPulseTimeRaw->Clear();
543 
544  frPosAdcPed->Clear();
545  frPosAdcPulseInt->Clear();
546  frPosAdcPulseAmp->Clear();
547  frPosAdcPulseTime->Clear();
548 
549  frNegAdcErrorFlag->Clear();
550  frNegAdcPedRaw->Clear();
551  frNegAdcThreshold->Clear();
552  frNegAdcPulseIntRaw->Clear();
553  frNegAdcPulseAmpRaw->Clear();
554  frNegAdcPulseTimeRaw->Clear();
555 
556  frNegAdcPed->Clear();
557  frNegAdcPulseInt->Clear();
558  frNegAdcPulseAmp->Clear();
559  frNegAdcPulseTime->Clear();
560 
561  /*
562  for(Int_t i=0;i<fNelem;i++) {
563 
564  fEpos[i] = 0;
565  fEneg[i] = 0;
566  fEmean[i] = 0;
567  }
568  */
569 
570  fEplane = 0;
571  fEplane_pos = 0;
572  fEplane_neg = 0;
573 
574  UInt_t nrPosAdcHits = 0;
575  UInt_t nrNegAdcHits = 0;
576 
577  // Process raw hits. Get ADC hits for the plane, assign variables for each
578  // channel.
579 
580  Int_t nrawhits = rawhits->GetLast()+1;
581 
582  Int_t ihit = nexthit;
583 
584  while(ihit < nrawhits) {
585  THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
586 
587  // This is OK as far as the hit list is sorted by layer.
588  //
589  if(hit->fPlane > fLayerNum) {
590  break;
591  }
592 
593  Int_t padnum = hit->fCounter;
594 
595  THcRawAdcHit& rawPosAdcHit = hit->GetRawAdcHitPos();
596  for (UInt_t thit=0; thit<rawPosAdcHit.GetNPulses(); ++thit) {
597  ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPedRaw());
598  ((THcSignalHit*) frPosAdcThreshold->ConstructedAt(nrPosAdcHits))->Set(padnum,rawPosAdcHit.GetPedRaw()*rawPosAdcHit.GetF250_PeakPedestalRatio()+fAdcPosThreshold);
599  ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPed());
600 
601  ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseIntRaw(thit));
602  ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseInt(thit));
603 
604  ((THcSignalHit*) frPosAdcPulseAmpRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseAmpRaw(thit));
605  ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseAmp(thit));
606 
607  ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseTimeRaw(thit));
608  ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseTime(thit)+fAdcTdcOffset);
609 
610  if (rawPosAdcHit.GetPulseAmp(thit)>0&&rawPosAdcHit.GetPulseIntRaw(thit)>0) {
611  ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,0);
612  } else {
613  ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,1);
614  }
615  if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) {
616  Double_t PeakPedRatio= rawPosAdcHit.GetF250_PeakPedestalRatio();
617  Int_t NPedSamples= rawPosAdcHit.GetF250_NPedestalSamples();
618  Double_t AdcToC = rawPosAdcHit.GetAdcTopC();
619  Double_t AdcToV = rawPosAdcHit.GetAdcTomV();
620  Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,0);
621  if (PedDefaultTemp !=0) {
622  Double_t tPulseInt = AdcToC*(rawPosAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
623  ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(padnum, tPulseInt);
624  ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, PedDefaultTemp);
625  ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
626 
627  }
628  ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, 0.);
629  }
630  ++nrPosAdcHits;
631  fTotNumAdcHits++;
632  fTotNumPosAdcHits++;
633 
634  }
635  THcRawAdcHit& rawNegAdcHit = hit->GetRawAdcHitNeg();
636  for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
637  ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPedRaw());
638  ((THcSignalHit*) frNegAdcThreshold->ConstructedAt(nrNegAdcHits))->Set(padnum,rawNegAdcHit.GetPedRaw()*rawNegAdcHit.GetF250_PeakPedestalRatio()+fAdcNegThreshold);
639  ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPed());
640 
641  ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseIntRaw(thit));
642  ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseInt(thit));
643 
644  ((THcSignalHit*) frNegAdcPulseAmpRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseAmpRaw(thit));
645  ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseAmp(thit));
646 
647  ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseTimeRaw(thit));
648  ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseTime(thit)+fAdcTdcOffset);
649 
650  if (rawNegAdcHit.GetPulseAmp(thit)>0&&rawNegAdcHit.GetPulseIntRaw(thit)>0) {
651  ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,0);
652  } else {
653  ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,1);
654  }
655  if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) {
656  Double_t PeakPedRatio= rawNegAdcHit.GetF250_PeakPedestalRatio();
657  Int_t NPedSamples= rawNegAdcHit.GetF250_NPedestalSamples();
658  Double_t AdcToC = rawNegAdcHit.GetAdcTopC();
659  Double_t AdcToV = rawNegAdcHit.GetAdcTomV();
660  Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,1);
661  if (PedDefaultTemp !=0) {
662  Double_t tPulseInt = AdcToC*(rawNegAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
663  ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(padnum, tPulseInt);
664  ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, PedDefaultTemp);
665  ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
666 
667  }
668  ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, 0.);
669  }
670  ++nrNegAdcHits;
671  fTotNumAdcHits++;
672  fTotNumNegAdcHits++;
673  }
674  ihit++;
675  }
676  return(ihit);
677 }
678 //_____________________________________________________________________________
679 Int_t THcShowerPlane::CoarseProcessHits()
680 {
681  Int_t ADCMode=static_cast<THcShower*>(fParent)->GetADCMode();
682  if(ADCMode == kADCDynamicPedestal) {
683  FillADC_DynamicPedestal();
684  } else if (ADCMode == kADCSampleIntegral) {
685  FillADC_SampleIntegral();
686  } else if (ADCMode == kADCSampIntDynPed) {
687  FillADC_SampIntDynPed();
688  } else {
689  FillADC_Standard();
690  }
691  //
692  if (static_cast<THcShower*>(fParent)->fdbg_decoded_cal) {
693 
694  cout << "---------------------------------------------------------------\n";
695  cout << "Debug output from THcShowerPlane::ProcessHits for "
696  << fParent->GetPrefix() << ":" << endl;
697 
698  cout << " Sparsified hits for HMS calorimeter plane #" << fLayerNum
699  << ", " << GetName() << ":" << endl;
700 
701  Int_t nspar = 0;
702 
703  for (Int_t i=0; i<fNelem; i++) {
704 
705  if (GetAposP(i) > 0 || GetAnegP(i) >0) { //hit
706  cout << " counter = " << i
707  << " Emean = " << fEmean[i]
708  << " Epos = " << fEpos[i]
709  << " Eneg = " << fEneg[i]
710  << endl;
711  nspar++;
712  }
713  }
714 
715  if (nspar == 0) cout << " No hits\n";
716 
717  cout << " Eplane = " << fEplane
718  << " Eplane_pos = " << fEplane_pos
719  << " Eplane_neg = " << fEplane_neg
720  << endl;
721  cout << "---------------------------------------------------------------\n";
722  }
723  //
724  return 1;
725  //
726 }
727 //_____________________________________________________________________________
728 void THcShowerPlane::FillADC_SampIntDynPed()
729 {
730  // adc_pos = hit->GetRawAdcHitPos().GetSampleInt();
731  // adc_neg = hit->GetRawAdcHitNeg().GetSampleInt();
732  // adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw();
733  // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw();
734  // Need to create
735 }
736 //_____________________________________________________________________________
737 void THcShowerPlane::FillADC_SampleIntegral()
738 {
740  // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw() - fNegPed[hit->fCounter -1];
741  // adc_pos = hit->GetRawAdcHitPos().GetSampleIntRaw();
742  // adc_neg = hit->GetRawAdcHitNeg().GetSampleIntRaw();
743  // need to create
744 }
745 //_____________________________________________________________________________
746 void THcShowerPlane::FillADC_Standard()
747 {
748  for (Int_t ielem=0;ielem<frNegAdcPulseIntRaw->GetEntries();ielem++) {
749  Int_t npad = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetPaddleNumber() - 1;
750  Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
751  fGoodNegAdcPulseIntRaw.at(npad) = pulseIntRaw;
752  if(fGoodNegAdcPulseIntRaw.at(npad) > fNegThresh[npad]) {
753  fGoodNegAdcPulseInt.at(npad) = pulseIntRaw-fNegPed[npad];
754  fEneg.at(npad) = fGoodNegAdcPulseInt.at(npad)*static_cast<THcShower*>(fParent)->GetGain(npad,fLayerNum-1,1);
755  fEmean.at(npad) += fEneg.at(npad);
756  fEplane_neg += fEneg.at(npad);
757  }
758  }
759  for (Int_t ielem=0;ielem<frPosAdcPulseIntRaw->GetEntries();ielem++) {
760  Int_t npad = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetPaddleNumber() - 1;
761  Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
762  fGoodPosAdcPulseIntRaw.at(npad) =pulseIntRaw;
763  if(fGoodPosAdcPulseIntRaw.at(npad) > fPosThresh[npad]) {
764  fGoodPosAdcPulseInt.at(npad) =pulseIntRaw-fPosPed[npad] ;
765  fEpos.at(npad) =fGoodPosAdcPulseInt.at(npad)*static_cast<THcShower*>(fParent)->GetGain(npad,fLayerNum-1,0);
766  fEmean.at(npad) += fEpos.at(npad);
767  fEplane_pos += fEpos.at(npad);
768  }
769  }
770  fEplane= fEplane_neg+fEplane_pos;
771 }
772 //_____________________________________________________________________________
773 void THcShowerPlane::FillADC_DynamicPedestal()
774 {
775  Double_t StartTime = 0.0;
776  if( fglHod ) StartTime = fglHod->GetStartTime();
777  Double_t OffsetTime = 0.0;
778  if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
779  for (Int_t ielem=0;ielem<frNegAdcPulseInt->GetEntries();ielem++) {
780  Int_t npad = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
781  Double_t pulseInt = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
782  Double_t pulsePed = ((THcSignalHit*) frNegAdcPed->ConstructedAt(ielem))->GetData();
783  Double_t pulseAmp = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
784  Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
785  Double_t pulseTime = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
786  Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
787  Double_t threshold = ((THcSignalHit*) frNegAdcThreshold->ConstructedAt(ielem))->GetData();
788  Bool_t pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,1)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,1) );
789  fGoodNegAdcMult.at(npad) += 1;
790  if (pulseTimeCut) {
791  fGoodNegAdcPulseIntRaw.at(npad) =pulseIntRaw;
792 
793  if(fGoodNegAdcPulseIntRaw.at(npad) > threshold && fGoodNegAdcPulseInt.at(npad)==0) {
794  fGoodNegAdcPulseInt.at(npad) =pulseInt ;
795  fEneg.at(npad) = fGoodNegAdcPulseInt.at(npad)*static_cast<THcShower*>(fParent)->GetGain(npad,fLayerNum-1,1);
796  fEmean.at(npad) += fEneg.at(npad);
797  fEplane_neg += fEneg.at(npad);
798 
799  fGoodNegAdcPed.at(npad) = pulsePed;
800  fGoodNegAdcPulseAmp.at(npad) = pulseAmp;
801  fGoodNegAdcPulseTime.at(npad) = pulseTime;
802  fGoodNegAdcTdcDiffTime.at(npad) = adctdcdiffTime;
803 
804  fTotNumGoodAdcHits++;
805  fTotNumGoodNegAdcHits++;
806  fNumGoodNegAdcHits.at(npad) = npad + 1;
807 
808  }
809 
810  }
811  }
812  //
813  for (Int_t ielem=0;ielem<frPosAdcPulseInt->GetEntries();ielem++) {
814  Int_t npad = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
815  Double_t pulsePed = ((THcSignalHit*) frPosAdcPed->ConstructedAt(ielem))->GetData();
816  Double_t threshold = ((THcSignalHit*) frPosAdcThreshold->ConstructedAt(ielem))->GetData();
817  Double_t pulseAmp = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(ielem))->GetData();
818  Double_t pulseInt = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData();
819  Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
820  Double_t pulseTime = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
821  Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
822  Bool_t pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,0)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,0) );
823  fGoodPosAdcMult.at(npad) += 1;
824  if (pulseTimeCut) {
825  fGoodPosAdcPulseIntRaw.at(npad) = pulseIntRaw;
826 
827  if(fGoodPosAdcPulseIntRaw.at(npad) > threshold && fGoodPosAdcPulseInt.at(npad)==0) {
828 
829  fGoodPosAdcPulseInt.at(npad) =pulseInt ;
830  fEpos.at(npad) = fGoodPosAdcPulseInt.at(npad)*static_cast<THcShower*>(fParent)->GetGain(npad,fLayerNum-1,0);
831  fEmean.at(npad) += fEpos[npad];
832  fEplane_pos += fEpos.at(npad);
833 
834  fGoodPosAdcPed.at(npad) = pulsePed;
835  fGoodPosAdcPulseAmp.at(npad) = pulseAmp;
836  fGoodPosAdcPulseTime.at(npad) = pulseTime;
837  fGoodPosAdcTdcDiffTime.at(npad) = adctdcdiffTime;
838 
839  fTotNumGoodAdcHits++;
840  fTotNumGoodPosAdcHits++;
841  fNumGoodPosAdcHits.at(npad) = npad + 1;
842 
843  }
844  }
845  }
846  //
847  fEplane= fEplane_neg+fEplane_pos;
848 
849 }
850 //_____________________________________________________________________________
851 Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
852 {
853  // Extract the data for this plane from hit list, accumulating into
854  // arrays for calculating pedestals.
855 
856  Int_t nrawhits = rawhits->GetLast()+1;
857 
858  Int_t ihit = nexthit;
859  while(ihit < nrawhits) {
860 
861  THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
862 
863  // OK for hit list sorted by layer.
864  if(hit->fPlane > fLayerNum) {
865  break;
866  }
867  Int_t element = hit->fCounter - 1; // Should check if in range
868  Int_t adcpos = hit->GetData(0);
869  Int_t adcneg = hit->GetData(1);
870 
871  if(adcpos <= fPosPedLimit[element]) {
872  fPosPedSum[element] += adcpos;
873  fPosPedSum2[element] += adcpos*adcpos;
874  fPosPedCount[element]++;
875  if(fPosPedCount[element] == fMinPeds/5) {
876  fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
877  }
878  }
879  if(adcneg <= fNegPedLimit[element]) {
880  fNegPedSum[element] += adcneg;
881  fNegPedSum2[element] += adcneg*adcneg;
882  fNegPedCount[element]++;
883  if(fNegPedCount[element] == fMinPeds/5) {
884  fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
885  }
886  }
887  ihit++;
888  }
889 
890  fNPedestalEvents++;
891 
892  // Debug output.
893 
894  if ( static_cast<THcShower*>(fParent)->fdbg_raw_cal ) {
895 
896  cout << "---------------------------------------------------------------\n";
897  cout << "Debug output from THcShowerPlane::AcculatePedestals for "
898  << fParent->GetPrefix() << ":" << endl;
899 
900  cout << "Processed hit list for plane " << GetName() << ":\n";
901 
902  for (Int_t ih=nexthit; ih<nrawhits; ih++) {
903 
904  THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ih);
905 
906  // OK for hit list sorted by layer.
907  if(hit->fPlane > fLayerNum) {
908  break;
909  }
910 
911  cout << " hit " << ih << ":"
912  << " plane = " << hit->fPlane
913  << " counter = " << hit->fCounter
914  << " ADCpos = " << hit->GetData(0)
915  << " ADCneg = " << hit->GetData(1)
916  << endl;
917  }
918 
919  cout << "---------------------------------------------------------------\n";
920 
921  }
922 
923  return(ihit);
924 }
925 
926 //_____________________________________________________________________________
927 void THcShowerPlane::CalculatePedestals( )
928 {
929  // Use the accumulated pedestal data to calculate pedestals
930  // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
931 
932  for(Int_t i=0; i<fNelem;i++) {
933 
934  // Positive tubes
935  fPosPed[i] = ((Float_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
936  fPosSig[i] = sqrt(((Float_t)fPosPedSum2[i])/TMath::Max(1, fPosPedCount[i])
937  - fPosPed[i]*fPosPed[i]);
938  fPosThresh[i] = fPosPed[i] + TMath::Min(50., TMath::Max(10., 3.*fPosSig[i]));
939 
940  // Negative tubes
941  fNegPed[i] = ((Float_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
942  fNegSig[i] = sqrt(((Float_t)fNegPedSum2[i])/TMath::Max(1, fNegPedCount[i])
943  - fNegPed[i]*fNegPed[i]);
944  fNegThresh[i] = fNegPed[i] + TMath::Min(50., TMath::Max(10., 3.*fNegSig[i]));
945 
946  }
947 
948  // Debug output.
949 
950  if ( static_cast<THcShower*>(fParent)->fdbg_raw_cal ) {
951 
952  cout << "---------------------------------------------------------------\n";
953  cout << "Debug output from THcShowerPlane::CalculatePedestals for "
954  << fParent->GetPrefix() << ":" << endl;
955 
956  cout << " ADC pedestals and thresholds for calorimeter plane "
957  << GetName() << endl;
958  for(Int_t i=0; i<fNelem;i++) {
959  cout << " element " << i << ": "
960  << " Pos. pedestal = " << fPosPed[i]
961  << " Pos. threshold = " << fPosThresh[i]
962  << " Neg. pedestal = " << fNegPed[i]
963  << " Neg. threshold = " << fNegThresh[i]
964  << endl;
965  }
966  cout << "---------------------------------------------------------------\n";
967 
968  }
969 
970 }
971 
972 //_____________________________________________________________________________
973 void THcShowerPlane::InitializePedestals( )
974 {
975  fNPedestalEvents = 0;
976  fPosPedSum = new Int_t [fNelem];
977  fPosPedSum2 = new Int_t [fNelem];
978  fPosPedCount = new Int_t [fNelem];
979  fNegPedSum = new Int_t [fNelem];
980  fNegPedSum2 = new Int_t [fNelem];
981  fNegPedCount = new Int_t [fNelem];
982 
983  fPosSig = new Float_t [fNelem];
984  fNegSig = new Float_t [fNelem];
985  fPosPed = new Float_t [fNelem];
986  fNegPed = new Float_t [fNelem];
987  fPosThresh = new Float_t [fNelem];
988  fNegThresh = new Float_t [fNelem];
989  for(Int_t i=0;i<fNelem;i++) {
990  fPosPedSum[i] = 0;
991  fPosPedSum2[i] = 0;
992  fPosPedCount[i] = 0;
993  fNegPedSum[i] = 0;
994  fNegPedSum2[i] = 0;
995  fNegPedCount[i] = 0;
996  }
997 }
998 
999 //_____________________________________________________________________________
1000 Int_t THcShowerPlane::AccumulateStat(TClonesArray& tracks )
1001 {
1002  // Accumumate statistics for efficiency calculations.
1003  //
1004  // Choose electron events in gas Cherenkov with good Chisq of the best track.
1005  // Project best track to the plane,
1006  // calculate row number for the track,
1007  // accrue number of tracks for the row,
1008  // accrue number of hits for the row, if row is hit.
1009  // Accrue total numbers of tracks and hits for plane.
1010 
1011  if(!fCherenkov) return 0;
1012 
1013  THaTrack* BestTrack = static_cast<THaTrack*>( tracks[0]);
1014  if (BestTrack->GetChi2()/BestTrack->GetNDoF() > fStatMaxChi2) return 0;
1015 
1016  if (fCherenkov->GetCerNPE() < fStatCerMin) return 0;
1017 
1018  Double_t XTrk = kBig;
1019  Double_t YTrk = kBig;
1020  Double_t pathl = kBig;
1021 
1022  // Track interception with plane. The coordinates are in the calorimeter's
1023  // local system.
1024 
1025  fOrigin = GetOrigin();
1026  static_cast<THcShower*>(fParent)->CalcTrackIntercept(BestTrack, pathl, XTrk, YTrk);
1027 
1028  // Transform coordiantes to the spectrometer's coordinate system.
1029  XTrk += GetOrigin().X();
1030  YTrk += GetOrigin().Y();
1031 
1032  for (Int_t i=0; i<fNelem; i++) {
1033 
1034  if (TMath::Abs(XTrk - static_cast<THcShower*>(fParent)->GetXPos(fLayerNum-1,i)) < fStatSlop &&
1035  YTrk > static_cast<THcShower*>(fParent)->GetYPos(fLayerNum-1,1) &&
1036  YTrk < static_cast<THcShower*>(fParent)->GetYPos(fLayerNum-1,0) ) {
1037 
1038  fStatNumTrk.at(i)++;
1039  fTotStatNumTrk++;
1040 
1041  if (fGoodPosAdcPulseInt.at(i) > 0. || fGoodNegAdcPulseInt.at(i) > 0.) {
1042  fStatNumHit.at(i)++;
1043  fTotStatNumHit++;
1044  }
1045 
1046  }
1047 
1048  }
1049 
1050  if ( static_cast<THcShower*>(fParent)->fdbg_tracks_cal ) {
1051  cout << "---------------------------------------------------------------\n";
1052  cout << "THcShowerPlane::AccumulateStat:" << endl;
1053  cout << " Chi2/NDF = " <<BestTrack->GetChi2()/BestTrack->GetNDoF() << endl;
1054  cout << " HGCER Npe = " << fCherenkov->GetCerNPE() << endl;
1055  cout << " XTrk, YTrk = " << XTrk << " " << YTrk << endl;
1056  for (Int_t i=0; i<fNelem; i++) {
1057  if (TMath::Abs(XTrk - static_cast<THcShower*>(fParent)->GetXPos(fLayerNum-1,i)) < fStatSlop) {
1058 
1059  cout << " Module " << i << ", X=" << static_cast<THcShower*>(fParent)->GetXPos(fLayerNum-1,i)
1060  << " matches track" << endl;
1061 
1062  if (fGoodPosAdcPulseInt.at(i) > 0. || fGoodNegAdcPulseInt.at(i) > 0.)
1063  cout << " PulseIntegrals = " << fGoodPosAdcPulseInt.at(i) << " "
1064  << fGoodNegAdcPulseInt.at(i) << endl;
1065  }
1066  }
1067  cout << "---------------------------------------------------------------\n";
1068  // getchar();
1069  }
1070 
1071  return 1;
1072 }
THcRawAdcHit::GetF250_PeakPedestalRatio
Double_t GetF250_PeakPedestalRatio()
Definition: THcRawAdcHit.h:24
GetName
std::string GetName(const std::string &scope_name)
THcShower::GetBlockThick
Double_t GetBlockThick(Int_t NLayer)
Definition: THcShower.h:50
THcShowerPlane::ReadDatabase
virtual Int_t ReadDatabase(const TDatime &date)
Definition: THcShowerPlane.cxx:171
Float_t
float Float_t
TObjArray::GetLast
Int_t GetLast() const
Option_t
const char Option_t
THcCherenkov
Class for gas Cherenkov detectors.
Definition: THcCherenkov.h:16
THcShower::GetYPos
Double_t GetYPos(Int_t NLayer, Int_t Side) const
Definition: THcShower.h:40
THcShowerPlane::DefineVariables
virtual Int_t DefineVariables(EMode mode=kDefine)
Definition: THcShowerPlane.cxx:318
THcShowerPlane::Init
virtual EStatus Init(const TDatime &run_time)
Definition: THcShowerPlane.cxx:132
THcRawAdcHit::GetPulseTimeRaw
Int_t GetPulseTimeRaw(UInt_t iPulse=0) const
Gets raw pulse time. In subsamples.
Definition: THcRawAdcHit.cxx:399
TMath::Min
Short_t Min(Short_t a, Short_t b)
Int_t
int Int_t
Bool_t
bool Bool_t
THcShowerPlane::FineProcess
virtual Int_t FineProcess(TClonesArray &tracks)
Definition: THcShowerPlane.cxx:520
std
STL namespace.
THcShowerPlane::ProcessHits
virtual Int_t ProcessHits(TClonesArray *rawhits, Int_t nexthit)
Definition: THcShowerPlane.cxx:526
THcShowerPlane::AccumulatePedestals
virtual Int_t AccumulatePedestals(TClonesArray *rawhits, Int_t nexthit)
Definition: THcShowerPlane.cxx:851
THcRawHit::fCounter
Int_t fCounter
Definition: THcRawHit.h:55
TMath::Abs
Short_t Abs(Short_t d)
threshold
UInt_t threshold[NUMSLOTS][NADCCHAN]
THcShowerPlane::FillADC_DynamicPedestal
virtual void FillADC_DynamicPedestal()
Definition: THcShowerPlane.cxx:773
THcShowerPlane::AccumulateStat
Int_t AccumulateStat(TClonesArray &tracks)
Definition: THcShowerPlane.cxx:1000
THcShower::GetZPos
Double_t GetZPos(Int_t NLayer) const
Definition: THcShower.h:48
THcRawShowerHit::GetData
virtual Int_t GetData(Int_t signal)
Definition: THcRawShowerHit.cxx:106
THcShower::fdbg_init_cal
Int_t fdbg_init_cal
Definition: THcShower.h:273
THcShower::GetAdcTdcOffset
Double_t GetAdcTdcOffset()
Definition: THcShower.h:122
app
tuple app
THcRawAdcHit::GetPed
Double_t GetPed() const
Gets sample pedestal. In channels.
Definition: THcRawAdcHit.cxx:428
THcShower::GetMinPeds
Int_t GetMinPeds()
Definition: THcShower.h:125
sqrt
VecExpr< UnaryOp< Sqrt< T >, SVector< T, D >, T >, T, D > sqrt(const SVector< T, D > &rhs)
ClassImp
ClassImp(THcShowerPlane) THcShowerPlane
Definition: THcShowerPlane.cxx:32
THcShowerPlane::FillADC_Standard
virtual void FillADC_Standard()
Definition: THcShowerPlane.cxx:746
THcRawShowerHit::GetRawAdcHitPos
THcRawAdcHit & GetRawAdcHitPos()
Definition: THcRawShowerHit.cxx:147
THcRawAdcHit::GetPulseAmpRaw
Int_t GetPulseAmpRaw(UInt_t iPulse=0) const
Gets raw pulse amplitude. In channels.
Definition: THcRawAdcHit.cxx:383
THcRawAdcHit::GetPulseAmp
Double_t GetPulseAmp(UInt_t iPulse=0) const
Gets pedestal subtracted pulse amplitude. In channels.
Definition: THcRawAdcHit.cxx:436
THcRawAdcHit::GetPulseInt
Double_t GetPulseInt(UInt_t iPulse=0) const
Gets pedestal subtracted pulse integral. In channels.
Definition: THcRawAdcHit.cxx:432
UInt_t
unsigned int UInt_t
Form
char * Form(const char *fmt,...)
THcShowerPlane::InitializePedestals
virtual void InitializePedestals()
Definition: THcShowerPlane.cxx:973
Warning
void Warning(const char *location, const char *msgfmt,...)
THcShowerPlane::Decode
virtual Int_t Decode(const THaEvData &)
Definition: THcShowerPlane.cxx:492
SConscript.list
tuple list
Definition: SConscript.py:9
THcRawAdcHit::GetPulseIntRaw
Int_t GetPulseIntRaw(UInt_t iPulse=0) const
Gets raw pulse integral. In channels.
Definition: THcRawAdcHit.cxx:367
THcShowerPlane::FillADC_SampIntDynPed
virtual void FillADC_SampIntDynPed()
Definition: THcShowerPlane.cxx:728
THcShowerPlane::CalculatePedestals
virtual void CalculatePedestals()
Definition: THcShowerPlane.cxx:927
THcParmList::LoadParmValues
Int_t LoadParmValues(const DBRequest *list, const char *prefix="")
Retrieve parameter values from the parameter cache.
Definition: THcParmList.cxx:522
Double_t
double Double_t
THcRawAdcHit::GetNPulses
UInt_t GetNPulses() const
Gets number of set pulses.
Definition: THcRawAdcHit.cxx:351
THcShowerPlane::FillADC_SampleIntegral
virtual void FillADC_SampleIntegral()
Definition: THcShowerPlane.cxx:737
THcShowerPlane::CoarseProcess
virtual Int_t CoarseProcess(TClonesArray &tracks)
Definition: THcShowerPlane.cxx:510
THcRawHit::fPlane
Int_t fPlane
Definition: THcRawHit.h:54
THcRawShowerHit
Class representing a single raw hit for a shower paddle.
Definition: THcRawShowerHit.h:8
THcShower
Generic segmented shower detector.
Definition: THcShower.h:18
THcShower::GetPedLimit
Int_t GetPedLimit(Int_t NBlock, Int_t NLayer, Int_t Side)
Definition: THcShower.h:52
THcShowerPlane
One plane of shower blocks with side readout.
Definition: THcShowerPlane.h:29
THcShowerPlane::~THcShowerPlane
virtual ~THcShowerPlane()
Definition: THcShowerPlane.cxx:78
THcRawAdcHit::GetPedRaw
Int_t GetPedRaw() const
Gets raw signal pedestal. In channels.
Definition: THcRawAdcHit.cxx:363
TMath::Max
Short_t Max(Short_t a, Short_t b)
THcShowerPlane::CoarseProcessHits
virtual Int_t CoarseProcessHits()
Definition: THcShowerPlane.cxx:679
THcRawAdcHit::GetPulseTime
Double_t GetPulseTime(UInt_t iPulse=0) const
Definition: THcRawAdcHit.cxx:440
gHcParms
R__EXTERN class THcParmList * gHcParms
Definition: THcGlobals.h:11
THaEvData
THcRawShowerHit::GetRawAdcHitNeg
THcRawAdcHit & GetRawAdcHitNeg()
Definition: THcRawShowerHit.cxx:152
TObjArray::At
TObject * At(Int_t idx) const
THcRawAdcHit::GetF250_NPedestalSamples
Int_t GetF250_NPedestalSamples()
Definition: THcRawAdcHit.h:25
GetData
void GetData(std::string s, double *x, double *y, double *ey)
kBig
static const Double_t kBig
Definition: THcFormula.cxx:31
THcShower::GetNBlocks
Int_t GetNBlocks(Int_t NLayer) const
Definition: THcShower.h:34
THcRawAdcHit::GetAdcTopC
Double_t GetAdcTopC() const
Definition: THcRawAdcHit.cxx:483
THcRawAdcHit
Class representing a single raw ADC hit.
Definition: THcRawAdcHit.h:7
THcRawAdcHit::GetAdcTomV
Double_t GetAdcTomV() const
Definition: THcRawAdcHit.cxx:477
THcShower::GetXPos
Double_t GetXPos(Int_t NLayer, Int_t NRow) const
Definition: THcShower.h:36
THcShowerPlane::Clear
virtual void Clear(Option_t *opt="")
Definition: THcShowerPlane.cxx:412
THcHallCSpectrometer
A standard Hall C spectrometer apparatus.
Definition: THcHallCSpectrometer.h:36