diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
index a171e940910a2fed19ee0c4fd5e24e9c09a80ad1..6581c5abbf6cadb2a09757911956c45d3148d3b3 100644
--- a/src/THcShowerArray.cxx
+++ b/src/THcShowerArray.cxx
@@ -48,6 +48,9 @@ THcShowerArray::~THcShowerArray()
delete [] fA;
delete [] fP;
+ delete [] fA_p;
+
+ delete [] fE;
}
//_____________________________________________________________________________
@@ -98,22 +101,103 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
fNelem = fNRows*fNColumns;
+ // Debug output.
+
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+
+ if (fParent->fdbg_init_cal) {
+ cout << "---------------------------------------------------------------\n";
+ cout << "Debug output from THcShowerArray::ReadDatabase for "
+ << GetParent()->GetPrefix() << ":" << endl;
+
+ cout << " Layer #" << fLayerNum << ", number of elements " << dec << fNelem
+ << endl;
+ cout << " Columns " << fNColumns << ", Rows " << fNRows << endl;
+
+ cout << " Using FADC " << fUsingFADC << endl;
+ if (fUsingFADC) {
+ cout << " FADC pedestal sample low = " << fPedSampLow << ", high = "
+ << fPedSampHigh << endl;
+ cout << " FADC data sample low = " << fDataSampLow << ", high = "
+ << fDataSampHigh << endl;
+ }
+
+ }
+
// Here read the 2-D arrays of pedestals, gains, etc.
// Pedestal limits per channel.
fPedLimit = new Int_t [fNelem];
+ Double_t cal_arr_cal_const[fNelem];
+ Double_t cal_arr_gain_cor[fNelem];
+
DBRequest list1[]={
- // {"cal_arr_cal_const", hcal_pos_cal_const, kDouble, fNelem},
{"cal_arr_ped_limit", fPedLimit, kInt, fNelem},
- // {"cal_arr_gain_cor", hcal_pos_gain_cor, kDouble, fNelem},
+ {"cal_arr_cal_const", cal_arr_cal_const, kDouble, fNelem},
+ {"cal_arr_gain_cor", cal_arr_gain_cor, kDouble, fNelem},
// {"cal_min_peds", &fShMinPeds, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list1, prefix);
- THcShower* fParent;
- fParent = (THcShower*) GetParent();
+ // Debug output.
+ if (fParent->fdbg_init_cal) {
+
+ cout << " fPedLimit:" << endl;
+ Int_t el=0;
+ for (UInt_t j=0; j<fNColumns; j++) {
+ cout << " ";
+ for (UInt_t i=0; i<fNRows; i++) {
+ cout << fPedLimit[el++] << " ";
+ };
+ cout << endl;
+ };
+
+ cout << " cal_arr_cal_const:" << endl;
+ el=0;
+ for (UInt_t j=0; j<fNColumns; j++) {
+ cout << " ";
+ for (UInt_t i=0; i<fNRows; i++) {
+ cout << cal_arr_cal_const[el++] << " ";
+ };
+ cout << endl;
+ };
+
+ cout << " cal_arr_gain_cor:" << endl;
+ el=0;
+ for (UInt_t j=0; j<fNColumns; j++) {
+ cout << " ";
+ for (UInt_t i=0; i<fNRows; i++) {
+ cout << cal_arr_gain_cor[el++] << " ";
+ };
+ cout << endl;
+ };
+
+ } // end of debug output
+
+ // Calibration constants (GeV / ADC channel).
+ fGain = new Double_t [fNelem];
+ for (UInt_t i=0; i<fNelem; i++) {
+ fGain[i] = cal_arr_cal_const[i] * cal_arr_gain_cor[i];
+ }
+
+ // Debug output.
+ if (fParent->fdbg_init_cal) {
+
+ cout << " fGain:" << endl;
+ Int_t el=0;
+ for (UInt_t j=0; j<fNColumns; j++) {
+ cout << " ";
+ for (UInt_t i=0; i<fNRows; i++) {
+ cout << fGain[el++] << " ";
+ };
+ cout << endl;
+ };
+
+ }
+
fMinPeds = fParent->GetMinPeds();
InitializePedestals();
@@ -121,6 +205,11 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
// Event by event amplitude and pedestal
fA = new Double_t[fNelem];
fP = new Double_t[fNelem];
+ fA_p = new Double_t[fNelem];
+
+ // Energy depositions per block.
+
+ fE = new Double_t[fNelem];
#ifdef HITPIC
hitpic = new char*[fNRows];
@@ -133,30 +222,12 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
// Debug output.
if (fParent->fdbg_init_cal) {
- cout << "---------------------------------------------------------------\n";
- cout << "Debug output from THcShowerArray::ReadDatabase for "
- << GetParent()->GetPrefix() << ":" << endl;
- cout << " Layer #" << fLayerNum << ", number of elements " << dec << fNelem
- << endl;
- cout << " Columns " << fNColumns << ", Rows " << fNRows << endl;
-
- cout << " Using FADC " << fUsingFADC << endl;
- if (fUsingFADC) {
- cout << " FADC pedestal sample low = " << fPedSampLow << ", high = "
- << fPedSampHigh << endl;
- cout << " FADC data sample low = " << fDataSampLow << ", high = "
- << fDataSampHigh << endl;
- }
+ cout << " fMinPeds = " << fMinPeds << endl;
// cout << " Origin of Layer at X = " << fOrigin.X()
// << " Y = " << fOrigin.Y() << " Z = " << fOrigin.Z() << endl;
- cout << " fPedLimit:";
- for(Int_t i=0;i<fNelem;i++) cout << " " << fPedLimit[i];
- cout << endl;
-
- cout << " fMinPeds = " << fMinPeds << endl;
cout << "---------------------------------------------------------------\n";
}
@@ -176,6 +247,8 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
{"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"},
{"a", "Raw ADC Amplitude", "fA"},
{"p", "Dynamic ADC Pedestal", "fP"},
+ {"a_p", "Sparsified, ped-subtracted ADC Amplitudes", "fA_p"},
+ {"e", "Energy Depositions per block", "fE"},
{ 0 }
};
@@ -228,8 +301,12 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
for(Int_t i=0;i<fNelem;i++) {
fA[i] = 0;
+ fA_p[i] = 0;
+ fE[i] = 0;
}
+ fETot = 0;
+
// Process raw hits. Get ADC hits for the plane, assign variables for each
// channel.
@@ -242,6 +319,10 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
while(ihit < nrawhits) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
+ if(hit->fPlane != fLayerNum) {
+ break;
+ }
+
// Should probably check that counter # is in range
if(fUsingFADC) {
fA[hit->fCounter-1] = hit->GetData(0,fPedSampLow,fPedSampHigh,
@@ -250,11 +331,30 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else {
fA[hit->fCounter-1] = hit->GetData(0);
}
+
if(fA[hit->fCounter-1] > threshold) {
ngood++;
}
- // Do other stuff like comparison to thresholds, signal hits, energy sums
+ // Sparsify hits, fill the hit list, compute the energy depostion.
+
+ Double_t thresh = fThresh[hit->fCounter -1];
+ if(fA[hit->fCounter-1] > thresh) {
+
+ // THcSignalHit *sighit =
+ // (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ // sighit->Set(hit->fCounter, fA_Pos[hit->fCounter-1]);
+
+ fUsingFADC ?
+ fA_p[hit->fCounter-1] = fA[hit->fCounter-1] :
+ fA_p[hit->fCounter-1] = fA[hit->fCounter-1] - fP[hit->fCounter -1];
+
+ fE[hit->fCounter-1] += fA_p[hit->fCounter-1] * fGain[hit->fCounter-1];
+ }
+
+ // Accumulate energies in the plane.
+
+ fETot += fE[hit->fCounter-1];
ihit++;
}
@@ -315,6 +415,40 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
#endif
+ //Debug output.
+
+ if (fParent->fdbg_decoded_cal) {
+
+ cout << "---------------------------------------------------------------\n";
+ cout << "Debug output from THcShowerArray::ProcessHits for "
+ << fParent->GetPrefix() << ":" << endl;
+
+ cout << " nrawhits = " << nrawhits << " nexthit = " << nexthit << endl;
+ cout << " Sparsified hits for shower array, plane #" << fLayerNum
+ << ", " << GetName() << ":" << endl;
+
+ Int_t nspar = 0;
+ for (Int_t jhit = nexthit; jhit < nrawhits; jhit++) {
+
+ THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(jhit);
+ if(hit->fPlane != fLayerNum) {
+ break;
+ }
+
+ if(fA[hit->fCounter-1] > fThresh[hit->fCounter -1]) {
+ cout << " counter = " << hit->fCounter
+ << " E = " << fE[hit->fCounter-1]
+ << endl;
+ nspar++;
+ }
+ }
+
+ if (nspar == 0) cout << " No hits\n";
+
+ cout << " E total = " << fETot << endl;
+ cout << "---------------------------------------------------------------\n";
+ }
+
return(ihit);
}
//_____________________________________________________________________________
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
index d3ca918d87a9ef58a9153bf177178848c0f8b380..48e3bb1a2821254a0db0004ca0e093996f70c363 100644
--- a/src/THcShowerArray.h
+++ b/src/THcShowerArray.h
@@ -63,8 +63,10 @@ protected:
char **hitpic;
Int_t piccolumn;
#endif
- Double_t* fA; // [fNelem] ADC amplitude of blocks
- Double_t* fP; // [fNelem] Event by event pedestals
+ Double_t* fA; // [fNelem] ADC amplitudes of blocks
+ Double_t* fP; // [fNelem] Event by event (FADC) pedestals
+ Double_t* fA_p; // [fNelem] sparsified, pedestal subtracted
+ // (FASTBUS) ADC amplitudes
TClonesArray* fADCHits; // List of ADC hits
@@ -95,6 +97,13 @@ protected:
Float_t *fSig; // [fNelem] pedestal rms-s
Float_t *fThresh; // [fNelem] ADC thresholds
+ Double_t* fGain; // [fNelem] Gain constants from calibration
+
+ //Energy depositions.
+
+ Double_t* fE; // [fNelem] energy depositions in the blocks.
+ Double_t fETot; // Total Energy deposition in the array.
+
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
ClassDef(THcShowerArray,0); // Fly;s Eye calorimeter array