From e9f33e2294770e4386008615fa9247fde0e95bca Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Tue, 28 May 2013 10:19:52 -0400
Subject: [PATCH] Added per plane energy deposition calculations
---
examples/hodtest.C | 5 ++--
examples/output.def | 6 +++++
src/THcShower.cxx | 57 +++++++++++++++++++++---------------------
src/THcShower.h | 13 ++++++++--
src/THcShowerPlane.cxx | 33 +++++++++++++++++++++++-
src/THcShowerPlane.h | 5 ++++
6 files changed, 86 insertions(+), 33 deletions(-)
diff --git a/examples/hodtest.C b/examples/hodtest.C
index dc58902..8968b15 100644
--- a/examples/hodtest.C
+++ b/examples/hodtest.C
@@ -1,4 +1,4 @@
-
+void hodtest(Int_t Nevents)
{
//
@@ -66,7 +66,8 @@
// run->SetEventRange(1,2000);// Physics Event number, does not
// include scaler or control events
// run->SetEventRange(1,999999);
- run->SetEventRange(1,1);
+ // run->SetEventRange(1,1);
+ run->SetEventRange(1,Nevents);
// Define the analysis parameters
analyzer->SetEvent( event );
diff --git a/examples/output.def b/examples/output.def
index a8fb853..ccbf7cd 100644
--- a/examples/output.def
+++ b/examples/output.def
@@ -34,6 +34,12 @@ TH1F chnegadc2 'HMS Cal B- ADC hits' H.cal.2ta.negadchits 13 0.5 13.5
TH1F chposadc3 'HMS Cal C+ ADC hits' H.cal.3ta.posadchits 13 0.5 13.5
TH1F chposadc4 'HMS Cal D+ ADC hits' H.cal.4ta.posadchits 13 0.5 13.5
+# Deposited Energies per Calorimeter layer
+TH1F edep1 'HMS Cal A Edep' H.cal.1pr.eplane 100 -0.1 2.074
+TH1F edep2 'HMS Cal B Edep' H.cal.2ta.eplane 100 -0.1 2.074
+TH1F edep3 'HMS Cal C Edep' H.cal.3ta.eplane 100 -0.1 2.074
+TH1F edep4 'HMS Cal D Edep' H.cal.4ta.eplane 100 -0.1 2.074
+
#Calorimeter ADC channels
TH1F hcaladc_A1p 'HMS Cal ADC A1p - PED' H.cal.1pr.apos_p[0] 150 50 500
TH1F hcaladc_A2p 'HMS Cal ADC A2p - PED' H.cal.1pr.apos_p[1] 150 50 500
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 23dde2c..7525a5f 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -215,15 +215,15 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
fNtotBlocks=0; //total number of blocks
for (Int_t i=0; i<fNLayers; i++) fNtotBlocks += fNBlocks[i];
- cout << "Total number of blocks in he calorimeter: " << fNtotBlocks << endl;
+ cout << "Total number of blocks in the calorimeter: " << fNtotBlocks << endl;
//Pedestal limits from hcal.param.
fShPosPedLimit = new Int_t [fNtotBlocks];
fShNegPedLimit = new Int_t [fNtotBlocks];
//Calibration constants
- fPosCalConst = new Double_t [fNtotBlocks];
- fNegCalConst = new Double_t [fNtotBlocks];
+ fPosGain = new Double_t [fNtotBlocks];
+ fNegGain = new Double_t [fNtotBlocks];
//Read in parameters from hcal.param
Double_t hcal_pos_cal_const[fNtotBlocks];
@@ -341,22 +341,22 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
//Calibration constants in GeV per ADC channel.
for (Int_t i=0; i<fNtotBlocks; i++) {
- fPosCalConst[i] = hcal_pos_cal_const[i] * hcal_pos_gain_cor[i];
- fNegCalConst[i] = hcal_neg_cal_const[i] * hcal_neg_gain_cor[i];
+ fPosGain[i] = hcal_pos_cal_const[i] * hcal_pos_gain_cor[i];
+ fNegGain[i] = hcal_neg_cal_const[i] * hcal_neg_gain_cor[i];
}
- cout << "fPosCalConst:" << endl;
+ cout << "fPosGain:" << endl;
for (Int_t j=0; j<fNLayers; j++) {
for (Int_t i=0; i<fNBlocks[j]; i++) {
- cout << fPosCalConst[j*fNBlocks[j]+i] << " ";
+ cout << fPosGain[j*fNBlocks[j]+i] << " ";
};
cout << endl;
};
- cout << "fNegCalConst:" << endl;
+ cout << "fNegGain:" << endl;
for (Int_t j=0; j<fNLayers; j++) {
for (Int_t i=0; i<fNBlocks[j]; i++) {
- cout << fNegCalConst[j*fNBlocks[j]+i] << " ";
+ cout << fNegGain[j*fNBlocks[j]+i] << " ";
};
cout << endl;
};
@@ -448,38 +448,37 @@ Int_t THcShower::Decode( const THaEvData& evdata )
// Get the Hall C style hitlist (fRawHitList) for this event
Int_t nhits = THcHitList::DecodeToHitList(evdata);
-if(gHaCuts->Result("Pedestal_event")) {
+ if(gHaCuts->Result("Pedestal_event")) {
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
-//cout << "nexthit = " << nexthit << endl;
+ //cout << "nexthit = " << nexthit << endl;
}
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return(0);
}
- if(fAnalyzePedestals) {
- for(Int_t ip=0;ip<fNLayers;ip++) {
- fPlanes[ip]->CalculatePedestals();
- }
- fAnalyzePedestals = 0; // Don't analyze pedestals next event
- }
-
-
+ if(fAnalyzePedestals) {
+ for(Int_t ip=0;ip<fNLayers;ip++) {
+ fPlanes[ip]->CalculatePedestals();
+ }
+ fAnalyzePedestals = 0; // Don't analyze pedestals next event
+ }
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
}
-/*
-// fRawHitList is TClones array of THcShowerHit objects
- for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
- THcShowerHit* hit = (THcShowerHit *) fRawHitList->At(ihit);
- cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
- << hit->fADC_pos << " " << hit->fADC_neg << " " << endl;
- }
- cout << endl;
-*/
+
+ // fRawHitList is TClones array of THcShowerHit objects
+ // cout << "THcShower::Decode: Shower raw hit list:" << endl;
+ // for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+ // THcShowerHit* hit = (THcShowerHit *) fRawHitList->At(ihit);
+ // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+ // << hit->fADC_pos << " " << hit->fADC_neg << " " << endl;
+ // }
+ // cout << endl;
+
return nhits;
}
@@ -507,6 +506,8 @@ Int_t THcShower::CoarseProcess( TClonesArray& ) //tracks
//
// static const Double_t sqrt2 = TMath::Sqrt(2.);
+ cout << "THcShower::CoarseProcess called ---------------------------" <<endl;
+
ApplyCorrections();
return 0;
diff --git a/src/THcShower.h b/src/THcShower.h
index 9c641fb..be952c9 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -47,6 +47,15 @@ public:
return ( Side == 0 ? fShPosPedLimit[nelem] : fShNegPedLimit[nelem]);
}
+ Double_t fGetGain(Int_t NBlock, Int_t NLayer, Int_t Side) {
+ if (Side!=0&&Side!=1) {
+ cout << "*** Wrong Side in fGetGain:" << Side << " ***" << endl;
+ return -1;
+ }
+ Int_t nelem = NBlock+NLayer*fNBlocks[NLayer];
+ return ( Side == 0 ? fPosGain[nelem] : fNegGain[nelem]);
+ }
+
Int_t fGetMinPeds() {
return fShMinPeds;
}
@@ -64,8 +73,8 @@ protected:
Int_t fShMinPeds; //Min.number of events to analize/update pedestals.
// Calibration constants
- Double_t* fPosCalConst;
- Double_t* fNegCalConst;
+ Double_t* fPosGain;
+ Double_t* fNegGain;
// Per-event data
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index 07299db..35b5a31 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -57,6 +57,9 @@ THcShowerPlane::~THcShowerPlane()
delete [] fA_Pos_p;
delete [] fA_Neg_p;
+ delete [] fEpos;
+ delete [] fEneg;
+ delete [] fEmean;
}
THaAnalysisObject::EStatus THcShowerPlane::Init( const TDatime& date )
@@ -116,6 +119,10 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
fA_Pos_p = new Float_t[fNelem];
fA_Neg_p = new Float_t[fNelem];
+ fEpos = new Float_t[fNelem];
+ fEneg = new Float_t[fNelem];
+ fEmean= new Float_t[fNelem];
+
// fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer();
//
// parname[plen]='\0';
@@ -172,6 +179,10 @@ Int_t THcShowerPlane::DefineVariables( EMode mode )
{"aneg", "Raw Negative ADC Amplitudes", "fA_Neg"},
{"apos_p", "Ped-subtracted Positive ADC Amplitudes", "fA_Pos_p"},
{"aneg_p", "Ped-subtracted Negative ADC Amplitudes", "fA_Neg_p"},
+ {"epos", "Energy Depositions from Positive Side PMTs", "fEpos"},
+ {"eneg", "Energy Depositions from Negative Side PMTs", "fEneg"},
+ {"emean", "Mean Energy Depositions", "fEMean"},
+ {"eplane", "Energy Deposition per plane", "fEplane"},
{ 0 }
};
@@ -202,6 +213,8 @@ Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks )
// HitCount();
+ cout << "THcShowerPlane::CoarseProcess called ---------------------" << endl;
+
return 0;
}
@@ -227,8 +240,13 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fA_Neg[i] = 0;
fA_Pos_p[i] = 0;
fA_Neg_p[i] = 0;
+ fEpos[i] = 0;
+ fEneg[i] = 0;
+ fEmean[i] = 0;
}
+ fEplane = 0;
+
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = nexthit;
@@ -249,19 +267,32 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fA_Pos_p[hit->fCounter-1] = hit->fADC_pos - fPosPed[hit->fCounter -1];
fA_Neg_p[hit->fCounter-1] = hit->fADC_neg - fNegPed[hit->fCounter -1];
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+
double thresh_pos = fPosThresh[hit->fCounter -1];
if(hit->fADC_pos > thresh_pos) {
- THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
sighit->Set(hit->fCounter, hit->fADC_pos);
+
+ fEpos[hit->fCounter-1] += fA_Pos_p[hit->fCounter-1]*
+ fParent->fGetGain(hit->fCounter-1,fLayerNum,0);
}
double thresh_neg = fNegThresh[hit->fCounter -1];
if(hit->fADC_neg > thresh_neg) {
+
THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
sighit->Set(hit->fCounter, hit->fADC_neg);
+
+ fEneg[hit->fCounter-1] += fA_Neg_p[hit->fCounter-1]*
+ fParent->fGetGain(hit->fCounter-1,fLayerNum,1);
}
+ fEmean[hit->fCounter-1] += (fEpos[hit->fCounter-1] + fEneg[hit->fCounter-1]);
+ fEplane += fEmean[hit->fCounter-1];
+
ihit++;
}
return(ihit);
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index baf5953..2e04130 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -53,6 +53,11 @@ class THcShowerPlane : public THaSubDetector {
Float_t* fA_Pos_p; // [fNelem] Array of pedestal subtracted ADC amplitudes
Float_t* fA_Neg_p; // [fNelem] Array of pedestal subtracted ADC amplitudes
+ Float_t* fEpos; // [fNelem] Array of energy depositions seen by positive PMTs
+ Float_t* fEneg; // [fNelem] Array of energy depositions seen by negative PMTs
+ Float_t* fEmean; // [fNelem] Array of mean energy depositions (pos + neg)
+ Float_t fEplane; // Energy deposition per plane
+
TClonesArray* fPosADCHits;
TClonesArray* fNegADCHits;
--
GitLab