From 6ad9abcf99e42a751ef054b4ef98d161b159f8d0 Mon Sep 17 00:00:00 2001
From: hallc-online <hallc-online@jlab.org>
Date: Mon, 27 Mar 2017 22:33:55 -0400
Subject: [PATCH] Modify THcShower, THcShowerPlane and THcShowerArray classes
Main purposes:
1) Add the possibility of multiple ADC hits per channel and
select best ADC by a time window.
2) Change ProcessHits method in ShowerPlane and ShowerArray to only
fill "raw" TClonesArrays.
3) Add CoarseProcessHits method to ShowerPlane and ShowerArray
to fill the "good" data arrays from call in THcShower:CoarseProcess
4) Moved code to fill Track object with energy from FineProcess
to CoarseProcess so that the Track->Get(Energy) can be used by other
detectors in their FineProcess.
5) FineProcess loops over Track object add gets Track energy
6) Add AdcErrorFlag TClonesArray to be filled to check for
problems with FADC
THcShower.h
1) Add methods GetADCMode(),GetAdcTimeWindowMin(),GetAdcTimeWindowMax()
so that the ShowerPlane can access these parameters
2) Add parameters ADCMode,fAdcTimeWindowMin,fAdcTimeWindowMax
3) Add variables
fNclustTrack : NUmber of cluster that matches best track
fXclustTrack : X pos of cluster that matches best track
fXTrack : X pos of best track that matches cluster
fEtrack : Cluster energy associated with best track
fEtrackNorm : Cluster energy/Track momentum associated with best track
THcShower.cxx
1) Add calls to THcShowerPlane and THcShowerArray CoarseProcessHits
in CoarseProcess
2) Moved filling of Energy in Track object to CoarseProcess
3) In FineProcess get the "best" Track object and match
Track energy and trajectory to cluster.
THcShowerArray.h and cxx
1)add parameters fADCMode, fAdcTimeWindowMin and fAdcTimeWindowMax
2) Add method CoarseProcessHits
3) In ProcessHits fill fPosThresh and fNegThresh with a fix
value of 250 integrated channels above pedestal.
Need to make this a parameter.
THcShowerPlane.h and cxx
1)Gets fADCMode, fAdcTimeWindowMin and fAdcTimeWindowMax from THcShower
2) Add method CoarseProcessHits
3) In ProcessHits fill fPosThresh and fNegThresh with a fix
value of 250 integrated channels above pedestal.
Need to make this a parameter.
---
src/THcShower.cxx | 119 +++++++++++++---------
src/THcShower.h | 28 ++++-
src/THcShowerArray.cxx | 199 +++++++++++++++++++++---------------
src/THcShowerArray.h | 13 ++-
src/THcShowerHit.cxx | 2 +-
src/THcShowerPlane.cxx | 226 +++++++++++++++++++++++++----------------
src/THcShowerPlane.h | 8 +-
7 files changed, 376 insertions(+), 219 deletions(-)
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 6c7d120..60fc502 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -209,9 +209,12 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
{
DBRequest list[]={
{"cal_num_neg_columns", &fNegCols, kInt, 0, 1},
- {"cal_slop", &fSlop, kDouble},
+ {"cal_slop", &fSlop, kDouble,0,1},
{"cal_fv_test", &fvTest, kInt,0,1},
- {"cal_fv_delta", &fvDelta, kDouble},
+ {"cal_fv_delta", &fvDelta, kDouble,0,1},
+ {"cal_ADCMode", &fADCMode, kInt, 0, 1},
+ {"cal_AdcTimeWindowMin", &fAdcTimeWindowMin, kDouble, 0, 1},
+ {"cal_AdcTimeWindowMax", &fAdcTimeWindowMax, kDouble, 0, 1},
{"dbg_raw_cal", &fdbg_raw_cal, kInt, 0, 1},
{"dbg_decoded_cal", &fdbg_decoded_cal, kInt, 0, 1},
{"dbg_sparsified_cal", &fdbg_sparsified_cal, kInt, 0, 1},
@@ -227,7 +230,9 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
fdbg_clusters_cal = 0;
fdbg_tracks_cal = 0;
fdbg_init_cal = 0;
-
+ fAdcTimeWindowMin=0;
+ fAdcTimeWindowMax=10000;
+ fADCMode=kADCDynamicPedestal;
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
}
@@ -343,6 +348,10 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
//Pedestal limits from hcal.param.
fShPosPedLimit = new Int_t [fNTotBlocks];
fShNegPedLimit = new Int_t [fNTotBlocks];
+ for (UInt_t i;i<fNTotBlocks;i++) {
+ fShPosPedLimit[i]=0.;
+ fShNegPedLimit[i]=0.;
+ }
//Calibration constants
fPosGain = new Double_t [fNTotBlocks];
@@ -350,31 +359,22 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
//Read in parameters from hcal.param
Double_t hcal_pos_cal_const[fNTotBlocks];
- // Double_t hcal_pos_gain_ini[fNTotBlocks]; not used
- // Double_t hcal_pos_gain_cur[fNTotBlocks]; not used
- // Int_t hcal_pos_ped_limit[fNTotBlocks]; not used
Double_t hcal_pos_gain_cor[fNTotBlocks];
Double_t hcal_neg_cal_const[fNTotBlocks];
- // Double_t hcal_neg_gain_ini[fNTotBlocks]; not used
- // Double_t hcal_neg_gain_cur[fNTotBlocks]; not used
- // Int_t hcal_neg_ped_limit[fNTotBlocks]; not used
Double_t hcal_neg_gain_cor[fNTotBlocks];
DBRequest list[]={
{"cal_pos_cal_const", hcal_pos_cal_const, kDouble, fNTotBlocks},
- // {"cal_pos_gain_ini", hcal_pos_gain_ini, kDouble, fNTotBlocks},
- // {"cal_pos_gain_cur", hcal_pos_gain_cur, kDouble, fNTotBlocks},
- {"cal_pos_ped_limit", fShPosPedLimit, kInt, fNTotBlocks},
+ {"cal_pos_ped_limit", fShPosPedLimit, kInt, fNTotBlocks,1},
{"cal_pos_gain_cor", hcal_pos_gain_cor, kDouble, fNTotBlocks},
{"cal_neg_cal_const", hcal_neg_cal_const, kDouble, fNTotBlocks},
- // {"cal_neg_gain_ini", hcal_neg_gain_ini, kDouble, fNTotBlocks},
- // {"cal_neg_gain_cur", hcal_neg_gain_cur, kDouble, fNTotBlocks},
- {"cal_neg_ped_limit", fShNegPedLimit, kInt, fNTotBlocks},
+ {"cal_neg_ped_limit", fShNegPedLimit, kInt, fNTotBlocks,1},
{"cal_neg_gain_cor", hcal_neg_gain_cor, kDouble, fNTotBlocks},
- {"cal_min_peds", &fShMinPeds, kInt},
+ {"cal_min_peds", &fShMinPeds, kInt,0,1},
{0}
};
+ fShMinPeds=0.;
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
// Debug output.
@@ -507,9 +507,13 @@ Int_t THcShower::DefineVariables( EMode mode )
RVarDef vars[] = {
{ "nhits", "Number of hits", "fNhits" },
{ "nclust", "Number of clusters", "fNclust" },
+ { "nclusttrack", "Number of cluster which matched best track","fNclustTrack" },
+ { "xclusttrack", "X pos of cluster which matched best track","fXclustTrack" },
+ { "xtrack", "X pos of track which matched cluster", "fXTrack" },
{ "etot", "Total energy", "fEtot" },
{ "etotnorm", "Total energy divided by Central Momentum", "fEtotNorm" },
{ "etrack", "Track energy", "fEtrack" },
+ { "etracknorm", "Total energy divided by track momentum", "fEtrackNorm" },
{ "ntracks", "Number of shower tracks", "fNtracks" },
{ 0 }
};
@@ -560,10 +564,14 @@ void THcShower::Clear(Option_t* opt)
fNhits = 0;
fNclust = 0;
+ fNclustTrack = -2;
+ fXclustTrack = -1000;
+ fXTrack = -1000;
fNtracks = 0;
fEtot = 0.;
fEtotNorm = 0.;
fEtrack = 0.;
+ fEtrackNorm = 0.;
// Purge cluster list
@@ -584,7 +592,7 @@ Int_t THcShower::Decode( const THaEvData& evdata )
// Get the Hall C style hitlist (fRawHitList) for this event
Int_t nhits = DecodeToHitList(evdata);
-
+
fEvent = evdata.GetEvNum();
if(gHaCuts->Result("Pedestal_event")) {
@@ -612,14 +620,10 @@ Int_t THcShower::Decode( const THaEvData& evdata )
Int_t nexthit = 0;
for(UInt_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
- fEtot += fPlanes[ip]->GetEplane();
}
if(fHasArray) {
nexthit = fArray->ProcessHits(fRawHitList, nexthit);
- fEtot += fArray->GetEarray();
}
- THcHallCSpectrometer *app = static_cast<THcHallCSpectrometer*>(GetApparatus());
- fEtotNorm=fEtot/(app->GetPcentral());
return nhits;
}
@@ -637,27 +641,25 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
//
// Fill set of unclustered hits.
-
+ for(UInt_t ip=0;ip<fNLayers;ip++) {
+ fPlanes[ip]->CoarseProcessHits();
+ fEtot += fPlanes[ip]->GetEplane();
+ }
+ if(fHasArray) {
+ fArray->CoarseProcessHits();
+ fEtot += fArray->GetEarray();
+ }
+ THcHallCSpectrometer *app = static_cast<THcHallCSpectrometer*>(GetApparatus());
+ fEtotNorm=fEtot/(app->GetPcentral());
+ //
THcShowerHitSet HitSet;
for(UInt_t j=0; j < fNLayers; j++) {
for (UInt_t i=0; i<fNBlocks[j]; i++) {
- //May be should be done this way.
- //
- // Double_t Edep = fPlanes[j]->GetEmean(i);
- // if (Edep > 0.) { //hit
- // Double_t x = fYPos[j][i] + BlockThick[j]/2.; //top + thick/2
- // Double_t z = fLayerZPos[j] + BlockThick[j]/2.;//front + thick/2
- // THcShowerHit* hit = new THcShowerHit(i,j,x,z,Edep);
-
- //ENGINE way.
//
- if (fPlanes[j]->GetApos(i) - fPlanes[j]->GetPosPed(i) >
- fPlanes[j]->GetPosThr(i) - fPlanes[j]->GetPosPed(i) ||
- fPlanes[j]->GetAneg(i) - fPlanes[j]->GetNegPed(i) >
- fPlanes[j]->GetNegThr(i) - fPlanes[j]->GetNegPed(i)) { //hit
+ if (fPlanes[j]->GetAposP(i) > 0 || fPlanes[j]->GetAnegP(i) >0) { //hit
Double_t Edep = fPlanes[j]->GetEmean(i);
Double_t Epos = fPlanes[j]->GetEpos(i);
@@ -682,6 +684,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
cout << "Debug output from THcShower::CoarseProcess for "
<< GetApparatus()->GetName() << endl;
+ cout << " event = " << fEvent << endl;
cout << " List of unclustered hits. Total hits: " << fNhits << endl;
THcShowerHitIt it = HitSet.begin(); //<set> version
for (Int_t i=0; i!=fNhits; i++) {
@@ -700,7 +703,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
if (fdbg_clusters_cal) {
- cout << " Clustered hits. Number of clusters: " << fNclust << endl;
+ cout << " event = " << fEvent << " Clustered hits. Number of clusters: " << fNclust << endl;
UInt_t i = 0;
for (THcShowerClusterListIt ppcl = (*fClusterList).begin();
@@ -729,7 +732,18 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
// Do same for Array.
if(fHasArray) fArray->CoarseProcess(tracks);
+ //
+ Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+ Double_t save_energy=0;
+ for (Int_t itrk=0; itrk<Ntracks; itrk++) {
+ THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
+ save_energy = GetShEnergy(theTrack);
+ if (fHasArray) save_energy += fArray->GetShEnergy(theTrack);
+ theTrack->SetEnergy(save_energy);
+ } //over tracks
+
+ //
return 0;
}
@@ -761,7 +775,6 @@ void THcShower::ClusterHits(THcShowerHitSet& HitSet,
k++) {
if ((**i).isNeighbour(*k)) {
-
(*cluster).insert(*i); //If the hit #i is neighbouring a hit
HitSet.erase(i); //in the cluster, then move it
//into the cluster.
@@ -945,7 +958,7 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
Double_t dx = TMath::Abs( clX(cluster) - XTrFront );
if (dx <= (0.5*BlockThick[0] + fSlop)) {
- fNtracks++; // number of shower tracks (Consistent with engine)
+ fNtracks++; // lumber of shower tracks (Consistent with engine)
if (dx < deltaX) {
mclust = i;
deltaX = dx;
@@ -960,6 +973,7 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShower::MatchCluster for "
<< GetApparatus()->GetName() << endl;
+ cout << " event = " << fEvent << endl;
cout << " Track at DC:"
<< " X = " << Track->GetX()
@@ -1026,6 +1040,7 @@ Float_t THcShower::GetShEnergy(THaTrack* Track) {
corneg = 0.;
}
+ // cout << ip << " clust energy pos = " << clEplane(cluster,ip,0)<< " clust energy pos = " << clEplane(cluster,ip,1) << endl;
Etrk += clEplane(cluster,ip,0) * corpos;
Etrk += clEplane(cluster,ip,1) * corneg;
@@ -1039,6 +1054,7 @@ Float_t THcShower::GetShEnergy(THaTrack* Track) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShower::GetShEnergy for "
<< GetApparatus()->GetName() << endl;
+ cout << " event = " << fEvent << endl;
cout << " Track at the calorimeter: X = " << Xtr << " Y = " << Ytr;
if (mclust >= 0)
@@ -1061,15 +1077,23 @@ Int_t THcShower::FineProcess( TClonesArray& tracks )
//
Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
-
for (Int_t itrk=0; itrk<Ntracks; itrk++) {
-
THaTrack* theTrack = static_cast<THaTrack*>( tracks[itrk] );
-
- fEtrack = GetShEnergy(theTrack);
- if (fHasArray) fEtrack += fArray->GetShEnergy(theTrack);
- theTrack->SetEnergy(fEtrack);
-
+ if (theTrack->GetIndex()==0) {
+ fEtrack=theTrack->GetEnergy();
+ fEtrackNorm=fEtrack/theTrack->GetP();
+
+ Double_t Xtr = -100.;
+ Double_t Ytr = -100.;
+ fNclustTrack = MatchCluster(theTrack, Xtr, Ytr);
+ fXTrack=Xtr;
+ if (fNclustTrack>=0) {
+ THcShowerCluster* cluster = *(fClusterList->begin()+fNclustTrack);
+ Double_t dx = TMath::Abs( clX(cluster) - Xtr );
+ fXclustTrack=clX(cluster);
+ //cout << fNclustTrack << " " << Xtr << " " << clX(cluster) << " " << dx << " " << Ytr << " " << fEtrack<< endl;
+ }
+ }
} //over tracks
//Debug output.
@@ -1078,7 +1102,7 @@ Int_t THcShower::FineProcess( TClonesArray& tracks )
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShower::FineProcess for "
<< GetApparatus()->GetName() << endl;
-
+ cout << " event = " << fEvent << endl;
cout << " Number of tracks = " << Ntracks << endl;
for (Int_t itrk=0; itrk<Ntracks; itrk++) {
@@ -1088,7 +1112,8 @@ Int_t THcShower::FineProcess( TClonesArray& tracks )
<< " Y = " << theTrack->GetY()
<< " Theta = " << theTrack->GetTheta()
<< " Phi = " << theTrack->GetPhi()
- << " Energy = " << theTrack->GetEnergy() << endl;
+ << " Energy = " << theTrack->GetEnergy()
+ << " Energy/Ptrack = " << fEtrackNorm << endl;
}
cout << "---------------------------------------------------------------\n";
diff --git a/src/THcShower.h b/src/THcShower.h
index b4ad6aa..5770cb4 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -71,6 +71,15 @@ public:
return ( Side == 0 ? fPosGain[nelem] : fNegGain[nelem]);
}
+ Int_t GetADCMode() {
+ return fADCMode;
+ }
+ Double_t GetAdcTimeWindowMin() {
+ return fAdcTimeWindowMin;
+ }
+ Double_t GetAdcTimeWindowMax() {
+ return fAdcTimeWindowMax;
+ }
Int_t GetMinPeds() {
return fShMinPeds;
}
@@ -78,6 +87,9 @@ public:
Int_t GetNLayers() {
return fNLayers;
}
+ Int_t GetNBlocks(Int_t layer) {
+ return fNBlocks[layer];
+ }
//Coordinate correction for single PMT modules.
//PMT attached at right (positive) side.
@@ -109,6 +121,16 @@ public:
protected:
Int_t fEvent;
+ Int_t fADCMode; // != 0 if using FADC
+ // 1 == Use the pulse int - pulse ped
+ // 2 == Use the sample integral - known ped
+ // 3 == Use the sample integral - sample ped
+ static const Int_t kADCStandard=0;
+ static const Int_t kADCDynamicPedestal=1;
+ static const Int_t kADCSampleIntegral=2;
+ static const Int_t kADCSampIntDynPed=3;
+ Double_t fAdcTimeWindowMin;
+ Double_t fAdcTimeWindowMax;
Int_t fAnalyzePedestals; // Flag for pedestal analysis.
@@ -124,11 +146,15 @@ protected:
Int_t fNhits; // Total number of hits
Int_t fNclust; // Number of clusters
+ Int_t fNclustTrack; // NUmber of cluster that match best track
+ Double_t fXclustTrack; // X pos of cluster that match best track
+ Double_t fXTrack; // X pos of best track that match cluster
Int_t fNtracks; // Number of shower tracks, i.e. number of
// cluster-to-track association
Double_t fEtot; // Total energy
Double_t fEtotNorm; // Total energy divided by spec central momentum
- Double_t fEtrack; // Cluster energy associated to the last track
+ Double_t fEtrack; // Cluster energy associated to the best track
+ Double_t fEtrackNorm; // Cluster energy associated to the best track
THcShowerClusterList* fClusterList; // List of hit clusters
diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
index da42572..63e5357 100644
--- a/src/THcShowerArray.cxx
+++ b/src/THcShowerArray.cxx
@@ -39,7 +39,8 @@ THcShowerArray::THcShowerArray( const char* name,
fADCHits = new TClonesArray("THcSignalHit",100);
fLayerNum = layernum;
- frAdcPedRaw = new TClonesArray("THcSignalHit", 16);
+ frAdcPedRaw = new TClonesArray("THcSignalHit", 16);
+ frAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
frAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
frAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
frAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
@@ -61,6 +62,7 @@ THcShowerArray::~THcShowerArray()
delete fADCHits;
delete frAdcPedRaw; frAdcPedRaw = NULL;
+ delete frAdcErrorFlag; frAdcErrorFlag = NULL;
delete frAdcPulseIntRaw; frAdcPulseIntRaw = NULL;
delete frAdcPulseAmpRaw; frAdcPulseAmpRaw = NULL;
delete frAdcPulseTimeRaw; frAdcPulseTimeRaw = NULL;
@@ -123,6 +125,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
{"cal_arr_ystep", &fYStep, kDouble},
{"cal_arr_zsize", &fZSize, kDouble},
{"cal_using_fadc", &fUsingFADC, kInt, 0, 1},
+ {"cal_arr_ADCMode", &fADCMode, kInt, 0, 1},
+ {"cal_arr_AdcTimeWindowMin", &fAdcTimeWindowMin, kDouble, 0, 1},
+ {"cal_arr_AdcTimeWindowMax", &fAdcTimeWindowMax, kDouble, 0, 1},
{"cal_ped_sample_low", &fPedSampLow, kInt, 0, 1},
{"cal_ped_sample_high", &fPedSampHigh, kInt, 0, 1},
{"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
@@ -130,7 +135,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
-
+ fADCMode=kADCDynamicPedestal;
+ fAdcTimeWindowMin=0;
+ fAdcTimeWindowMax=10000;
fNelem = fNRows*fNColumns;
fXPos = new Double_t* [fNRows];
@@ -216,10 +223,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
Double_t cal_arr_gain_cor[fNelem];
DBRequest list1[]={
- {"cal_arr_ped_limit", fPedLimit, kInt, static_cast<UInt_t>(fNelem)},
+ {"cal_arr_ped_limit", fPedLimit, kInt, static_cast<UInt_t>(fNelem),1},
{"cal_arr_cal_const", cal_arr_cal_const, kDouble, static_cast<UInt_t>(fNelem)},
{"cal_arr_gain_cor", cal_arr_gain_cor, kDouble, static_cast<UInt_t>(fNelem)},
- // {"cal_min_peds", &fShMinPeds, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list1, prefix);
@@ -261,7 +267,7 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
// Calibration constants (GeV / ADC channel).
fGain = new Double_t [fNelem];
- for (UInt_t i=0; i<fNelem; i++) {
+ for (Int_t i=0; i<fNelem; i++) {
fGain[i] = cal_arr_cal_const[i] * cal_arr_gain_cor[i];
}
@@ -328,6 +334,7 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
{"p", "Dynamic ADC Pedestal", "fP"},
{"a_p", "Sparsified, ped-subtracted ADC Amplitudes", "fA_p"},
{ "nhits", "Number of hits", "fNhits" },
+ { "nghits", "Number of good hits ( pass threshold on raw ADC)", "fNgoodhits" },
{ "nclust", "Number of clusters", "fNclust" },
{"e", "Energy Depositions per block", "fE"},
{"earray", "Energy Deposition in array", "fEarray"},
@@ -336,6 +343,7 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
{"adcCounter", "List of ADC counter numbers.", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"adcPedRaw", "List of raw ADC pedestals", "frAdcPedRaw.THcSignalHit.GetData()"},
+ {"adcErrorFlag", "List of raw ADC pedestals", "frAdcErrorFlag.THcSignalHit.GetData()"},
{"adcPulseIntRaw", "List of raw ADC pulse integrals.", "frAdcPulseIntRaw.THcSignalHit.GetData()"},
{"adcPulseAmpRaw", "List of raw ADC pulse amplitudes.", "frAdcPulseAmpRaw.THcSignalHit.GetData()"},
{"adcPulseTimeRaw", "List of raw ADC pulse times.", "frAdcPulseTimeRaw.THcSignalHit.GetData()"},
@@ -357,6 +365,7 @@ void THcShowerArray::Clear( Option_t* )
fADCHits->Clear();
fNhits = 0;
+ fNgoodhits = 0;
fNclust = 0;
fNtracks = 0;
@@ -368,6 +377,7 @@ void THcShowerArray::Clear( Option_t* )
fClusterList->clear();
frAdcPedRaw->Clear();
+ frAdcErrorFlag->Clear();
frAdcPulseIntRaw->Clear();
frAdcPulseAmpRaw->Clear();
frAdcPulseTimeRaw->Clear();
@@ -413,9 +423,6 @@ Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks )
k++;
}
}
-
- fNhits = HitSet.size();
-
//Debug output, print out hits before clustering.
THcShower* fParent = (THcShower*) GetParent();
@@ -427,7 +434,7 @@ Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks )
cout << " List of unclustered hits. Total hits: " << fNhits << endl;
THcShowerHitIt it = HitSet.begin(); //<set> version
- for (Int_t i=0; i!=fNhits; i++) {
+ for (Int_t i=0; i!=fNgoodhits; i++) {
cout << " hit " << i << ": ";
(*(it++))->show();
}
@@ -643,6 +650,96 @@ Int_t THcShowerArray::FineProcess( TClonesArray& tracks )
return 0;
}
+//_____________________________________________________________________________
+Int_t THcShowerArray::CoarseProcessHits()
+{
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+ Int_t ADCMode=fParent->GetADCMode();
+ if(ADCMode == kADCDynamicPedestal) {
+ FillADC_DynamicPedestal();
+ } else if (ADCMode == kADCSampleIntegral) {
+ FillADC_SampleIntegral();
+ } else if (ADCMode == kADCSampIntDynPed) {
+ FillADC_SampIntDynPed();
+ } else {
+ FillADC_Standard();
+ }
+ //
+ if (fParent->fdbg_decoded_cal) {
+
+ cout << "---------------------------------------------------------------\n";
+ cout << "Debug output from THcShowerArray::ProcessHits for "
+ << fParent->GetPrefix() << ":" << endl;
+
+ cout << " Sparsified hits for shower array, plane #" << fLayerNum
+ << ", " << GetName() << ":" << endl;
+
+ Int_t nspar = 0;
+ Int_t k=0;
+ for(UInt_t j=0; j < fNColumns; j++) {
+ for (UInt_t i=0; i<fNRows; i++) {
+ if(fA[k] > fThresh[k]) {
+ cout << " counter = " << k
+ << " E = " << fE[k]
+ << endl;
+ nspar++;
+ }
+ k++;
+ }
+ }
+
+ if (nspar == 0) cout << " No hits\n";
+
+ cout << " Earray = " << fEarray << endl;
+ cout << "---------------------------------------------------------------\n";
+ }
+ //
+ return 1;
+}
+//_____________________________________________________________________________
+void THcShowerArray::FillADC_SampIntDynPed()
+{
+ // adc_pos = hit->GetRawAdcHitPos().GetSampleInt();
+ // adc_neg = hit->GetRawAdcHitNeg().GetSampleInt();
+ // adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw();
+ // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw();
+ // Need to create
+}
+//_____________________________________________________________________________
+void THcShowerArray::FillADC_SampleIntegral()
+{
+ /// adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw() - fPosPed[hit->fCounter -1];
+ // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw() - fNegPed[hit->fCounter -1];
+ // adc_pos = hit->GetRawAdcHitPos().GetSampleIntRaw();
+ // adc_neg = hit->GetRawAdcHitNeg().GetSampleIntRaw();
+ // need to create
+}
+//_____________________________________________________________________________
+void THcShowerArray::FillADC_Standard()
+{
+}
+//_____________________________________________________________________________
+void THcShowerArray::FillADC_DynamicPedestal()
+{
+ for (Int_t ielem=0;ielem<frAdcPulseInt->GetEntries();ielem++) {
+ Int_t npad = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+ Double_t pulseInt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData();
+ Double_t pulseIntRaw = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+ Double_t pulseTime = ((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+ Double_t errorflag = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData();
+ Bool_t pulseTimeCut = (pulseTime > fAdcTimeWindowMin) && (pulseTime < fAdcTimeWindowMax);
+ if (errorflag==0 && pulseTimeCut) {
+ fA[npad] =pulseIntRaw;
+ if(fA[npad] > fThresh[npad]) {
+ fA_p[npad] =pulseInt ;
+ fE[npad] = fA_p[npad]*fGain[npad];
+ fEarray += fE[npad];
+ }
+ }
+ }
+ //
+}
//_____________________________________________________________________________
Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
{
@@ -656,6 +753,7 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fADCHits->Clear();
frAdcPedRaw->Clear();
+ frAdcErrorFlag->Clear();
frAdcPulseIntRaw->Clear();
frAdcPulseAmpRaw->Clear();
frAdcPulseTimeRaw->Clear();
@@ -679,9 +777,6 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
Int_t ihit = nexthit;
- Int_t ngood = 0;
- Int_t threshold = 100;
-
UInt_t nrAdcHits = 0;
while(ihit < nrawhits) {
@@ -692,11 +787,12 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
Int_t padnum = hit->fCounter;
-
THcRawAdcHit& rawAdcHit = hit->GetRawAdcHitPos();
+ //
for (UInt_t thit=0; thit<rawAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPedRaw());
- ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPed());
+ fThresh[padnum]=rawAdcHit.GetPedRaw()*rawAdcHit.GetF250_PeakPedestalRatio()+250.;
+ ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPed());
((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseIntRaw(thit));
((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseInt(thit));
@@ -705,48 +801,18 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseAmp(thit));
((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseTimeRaw(thit));
-
+ if (rawAdcHit.GetPulseAmp(thit)>0&&rawAdcHit.GetPulseIntRaw(thit)>0) {
+ ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(padnum,0);
+ } else {
+ ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(padnum,1);
+ }
++nrAdcHits;
}
- fADCMode=kADCDynamicPedestal;
- Double_t adc;
- Double_t adc_pedsub;
- if(fADCMode == kADCDynamicPedestal) {
- adc_pedsub = hit->GetRawAdcHitPos().GetPulseInt();
- adc= hit->GetRawAdcHitPos().GetPulseIntRaw();
- } else if (fADCMode == kADCSampleIntegral) {
- adc_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw() - fPed[hit->fCounter -1];
- adc = hit->GetRawAdcHitPos().GetSampleIntRaw();
- } else if (fADCMode == kADCSampIntDynPed) {
- adc = hit->GetRawAdcHitPos().GetSampleInt();
- adc_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw();
- } else {
- adc_pedsub = hit->GetRawAdcHitPos().GetPulseIntRaw()-fPed[hit->fCounter -1];
- adc = hit->GetRawAdcHitPos().GetPulseIntRaw();
- }
-
- fA[hit->fCounter-1] = adc;
- threshold=hit->GetRawAdcHitPos().GetPedRaw()+100;
- if(fA[hit->fCounter-1] > threshold) {
- ngood++;
- }
-
- // Sparsify hits, fill the hit list, compute the energy depostion.
- fThresh[hit->fCounter -1] = hit->GetRawAdcHitPos().GetPedRaw()+100;
- if(fA[hit->fCounter-1] > fThresh[hit->fCounter -1]) {
- fA_p[hit->fCounter-1] = adc_pedsub;
- fE[hit->fCounter-1] += fA_p[hit->fCounter-1] * fGain[hit->fCounter-1];
- }
-
- // Accumulate energies in the plane.
-
- fEarray += fE[hit->fCounter-1];
-
ihit++;
}
#if 0
- if(ngood > 0) {
+ if(fNgoodhits > 0) {
cout << "+";
for(Int_t column=0;column<fNColumns;column++) {
cout << "-";
@@ -767,7 +833,7 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
#endif
#ifdef HITPIC
- if(ngood > 0) {
+ if(fNgoodhits > 0) {
for(Int_t row=0;row<fNRows;row++) {
if(piccolumn==0) {
hitpic[row][0] = '|';
@@ -803,37 +869,6 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
//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 << " Earray = " << fEarray << endl;
- cout << "---------------------------------------------------------------\n";
- }
return(ihit);
}
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
index ad01670..4b178cb 100644
--- a/src/THcShowerArray.h
+++ b/src/THcShowerArray.h
@@ -45,10 +45,15 @@ public:
virtual Bool_t IsPid() { return kFALSE; }
virtual Int_t ProcessHits(TClonesArray* rawhits, Int_t nexthit);
+ virtual Int_t CoarseProcessHits();
virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit);
virtual void CalculatePedestals( );
virtual void InitializePedestals( );
-
+ virtual void FillADC_DynamicPedestal( );
+ virtual void FillADC_SampleIntegral( );
+ virtual void FillADC_SampIntDynPed( );
+ virtual void FillADC_Standard( );
+
// Cluster to track association method.
Int_t MatchCluster(THaTrack*, Double_t&, Double_t&);
@@ -106,7 +111,9 @@ protected:
static const Int_t kADCDynamicPedestal=1;
static const Int_t kADCSampleIntegral=2;
static const Int_t kADCSampIntDynPed=3;
- Int_t fPedSampLow; // Sample range for
+ Double_t fAdcTimeWindowMin ;
+ Double_t fAdcTimeWindowMax ;
+Int_t fPedSampLow; // Sample range for
Int_t fPedSampHigh; // dynamic pedestal
Int_t fDataSampLow; // Sample range for
Int_t fDataSampHigh; // sample integration
@@ -137,6 +144,7 @@ protected:
Double_t fEarray; // Total Energy deposition in the array.
Int_t fNhits; // Total number of hits
+ Int_t fNgoodhits; // Total number of good hits (pass threshold)
Int_t fNclust; // Number of hit clusters
Int_t fNtracks; // Number of shower tracks, i.e. number of
// cluster-to-track associations
@@ -144,6 +152,7 @@ protected:
THcShowerClusterList* fClusterList; // List of hit clusters
TClonesArray* frAdcPedRaw;
+ TClonesArray* frAdcErrorFlag;
TClonesArray* frAdcPulseIntRaw;
TClonesArray* frAdcPulseAmpRaw;
TClonesArray* frAdcPulseTimeRaw;
diff --git a/src/THcShowerHit.cxx b/src/THcShowerHit.cxx
index 574d056..e7bc7bf 100644
--- a/src/THcShowerHit.cxx
+++ b/src/THcShowerHit.cxx
@@ -44,7 +44,7 @@ bool THcShowerHit::isNeighbour(THcShowerHit* hit1) {
//Print out hit information
//
void THcShowerHit::show() {
- cout << "row=" << fRow << " column=" << fCol
+ cout << "row=" << fRow+1 << " column=" << fCol+1
<< " x=" << fX << " z=" << fZ
<< " E=" << fE << " Epos=" << fEpos << " Eneg=" << fEneg << endl;
}
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index 44eba1a..0a9452d 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -39,6 +39,7 @@ THcShowerPlane::THcShowerPlane( const char* name,
fPosADCHits = new TClonesArray("THcSignalHit",fNelem);
fNegADCHits = new TClonesArray("THcSignalHit",fNelem);
+ frPosAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
frPosAdcPedRaw = new TClonesArray("THcSignalHit", 16);
frPosAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
frPosAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
@@ -48,6 +49,7 @@ THcShowerPlane::THcShowerPlane( const char* name,
frPosAdcPulseInt = new TClonesArray("THcSignalHit", 16);
frPosAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
+ frNegAdcErrorFlag = new TClonesArray("THcSignalHit", 16);
frNegAdcPedRaw = new TClonesArray("THcSignalHit", 16);
frNegAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
frNegAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
@@ -72,6 +74,7 @@ THcShowerPlane::~THcShowerPlane()
delete fPosADCHits;
delete fNegADCHits;
+ frPosAdcErrorFlag = NULL;
frPosAdcPedRaw = NULL;
frPosAdcPulseIntRaw = NULL;
frPosAdcPulseAmpRaw = NULL;
@@ -81,6 +84,7 @@ THcShowerPlane::~THcShowerPlane()
frPosAdcPulseInt = NULL;
frPosAdcPulseAmp = NULL;
+ frNegAdcErrorFlag = NULL;
frNegAdcPedRaw = NULL;
frNegAdcPulseIntRaw = NULL;
frNegAdcPulseAmpRaw = NULL;
@@ -256,6 +260,7 @@ Int_t THcShowerPlane::DefineVariables( EMode mode )
{"posAdcCounter", "List of positive ADC counter numbers.", "frPosAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"negAdcCounter", "List of negative ADC counter numbers.", "frNegAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
+ {"posAdcErrorFlag", "List of positive raw ADC Error Flags", "frPosAdcErrorFlag.THcSignalHit.GetData()"},
{"posAdcPedRaw", "List of positive raw ADC pedestals", "frPosAdcPedRaw.THcSignalHit.GetData()"},
{"posAdcPulseIntRaw", "List of positive raw ADC pulse integrals.", "frPosAdcPulseIntRaw.THcSignalHit.GetData()"},
{"posAdcPulseAmpRaw", "List of positive raw ADC pulse amplitudes.", "frPosAdcPulseAmpRaw.THcSignalHit.GetData()"},
@@ -265,6 +270,7 @@ Int_t THcShowerPlane::DefineVariables( EMode mode )
{"posAdcPulseInt", "List of positive ADC pulse integrals.", "frPosAdcPulseInt.THcSignalHit.GetData()"},
{"posAdcPulseAmp", "List of positive ADC pulse amplitudes.", "frPosAdcPulseAmp.THcSignalHit.GetData()"},
+ {"negAdcErrorFlag", "List of negative raw ADC Error Flag ", "frNegAdcErrorFlag.THcSignalHit.GetData()"},
{"negAdcPedRaw", "List of negative raw ADC pedestals", "frNegAdcPedRaw.THcSignalHit.GetData()"},
{"negAdcPulseIntRaw", "List of negative raw ADC pulse integrals.", "frNegAdcPulseIntRaw.THcSignalHit.GetData()"},
{"negAdcPulseAmpRaw", "List of negative raw ADC pulse amplitudes.", "frNegAdcPulseAmpRaw.THcSignalHit.GetData()"},
@@ -287,6 +293,7 @@ void THcShowerPlane::Clear( Option_t* )
fPosADCHits->Clear();
fNegADCHits->Clear();
+ frPosAdcErrorFlag->Clear();
frPosAdcPedRaw->Clear();
frPosAdcPulseIntRaw->Clear();
frPosAdcPulseAmpRaw->Clear();
@@ -296,6 +303,7 @@ void THcShowerPlane::Clear( Option_t* )
frPosAdcPulseInt->Clear();
frPosAdcPulseAmp->Clear();
+ frNegAdcErrorFlag->Clear();
frNegAdcPedRaw->Clear();
frNegAdcPulseIntRaw->Clear();
frNegAdcPulseAmpRaw->Clear();
@@ -365,6 +373,7 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fPosADCHits->Clear();
fNegADCHits->Clear();
+ frPosAdcErrorFlag->Clear();
frPosAdcPedRaw->Clear();
frPosAdcPulseIntRaw->Clear();
frPosAdcPulseAmpRaw->Clear();
@@ -374,6 +383,7 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
frPosAdcPulseInt->Clear();
frPosAdcPulseAmp->Clear();
+ frNegAdcErrorFlag->Clear();
frNegAdcPedRaw->Clear();
frNegAdcPulseIntRaw->Clear();
frNegAdcPulseAmpRaw->Clear();
@@ -421,6 +431,7 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
THcRawAdcHit& rawPosAdcHit = hit->GetRawAdcHitPos();
for (UInt_t thit=0; thit<rawPosAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPedRaw());
+ fPosThresh[hit->fCounter -1]=rawPosAdcHit.GetPedRaw()*rawPosAdcHit.GetF250_PeakPedestalRatio()+250.;
((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPed());
((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseIntRaw(thit));
@@ -430,12 +441,17 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseAmp(thit));
((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseTimeRaw(thit));
-
+ if (rawPosAdcHit.GetPulseAmp(thit)>0&&rawPosAdcHit.GetPulseIntRaw(thit)>0) {
+ ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,0);
+ } else {
+ ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,1);
+ }
++nrPosAdcHits;
}
THcRawAdcHit& rawNegAdcHit = hit->GetRawAdcHitNeg();
for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPedRaw());
+ fNegThresh[hit->fCounter -1]=rawNegAdcHit.GetPedRaw()*rawNegAdcHit.GetF250_PeakPedestalRatio()+250.;
((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPed());
((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseIntRaw(thit));
@@ -445,107 +461,54 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseAmp(thit));
((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseTimeRaw(thit));
-
+ if (rawNegAdcHit.GetPulseAmp(thit)>0&&rawNegAdcHit.GetPulseIntRaw(thit)>0) {
+ ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,0);
+ } else {
+ ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,1);
+ }
++nrNegAdcHits;
}
- fADCMode=kADCDynamicPedestal;
- Double_t adc_pos;
- Double_t adc_neg;
- Double_t adc_pos_pedsub;
- Double_t adc_neg_pedsub;
- if(fADCMode == kADCDynamicPedestal) {
- adc_pos_pedsub = hit->GetRawAdcHitPos().GetPulseInt();
- adc_neg_pedsub = hit->GetRawAdcHitNeg().GetPulseInt();
- adc_pos = hit->GetRawAdcHitPos().GetPulseIntRaw();
- adc_neg = hit->GetRawAdcHitNeg().GetPulseIntRaw();
- } else if (fADCMode == kADCSampleIntegral) {
- adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw() - fPosPed[hit->fCounter -1];
- adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw() - fNegPed[hit->fCounter -1];
- adc_pos = hit->GetRawAdcHitPos().GetSampleIntRaw();
- adc_neg = hit->GetRawAdcHitNeg().GetSampleIntRaw();
- } else if (fADCMode == kADCSampIntDynPed) {
- adc_pos = hit->GetRawAdcHitPos().GetSampleInt();
- adc_neg = hit->GetRawAdcHitNeg().GetSampleInt();
- adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw();
- adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw();
- } else {
- adc_pos_pedsub = hit->GetRawAdcHitPos().GetPulseIntRaw()-fPosPed[hit->fCounter -1];
- adc_neg_pedsub = hit->GetRawAdcHitNeg().GetPulseIntRaw()-fNegPed[hit->fCounter -1];
- adc_pos = hit->GetRawAdcHitPos().GetPulseIntRaw();
- adc_neg = hit->GetRawAdcHitNeg().GetPulseIntRaw();
- }
- //
- // Sparsify positive side hits, fill the hit list, compute the
- // energy depostion from positive side for the counter.
- fA_Pos[hit->fCounter-1] =adc_pos;
- fA_Neg[hit->fCounter-1] =adc_neg;
-
- Double_t thresh_pos = fPosThresh[hit->fCounter -1];
- thresh_pos =0;
- if(fA_Pos[hit->fCounter-1] > thresh_pos) {
-
- fA_Pos_p[hit->fCounter-1] =adc_pos_pedsub ;
-
- // cout << " pos " << hit->fCounter << " " << fA_Pos_p[hit->fCounter-1] << " " << fA_Pos[hit->fCounter-1] << " " << fPosPed[hit->fCounter -1] << "" << frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits-1))->Get(padnum)<< endl;
- fEpos[hit->fCounter-1] += fA_Pos_p[hit->fCounter-1]*
- fParent->GetGain(hit->fCounter-1,fLayerNum-1,0);
- }
-
- // Sparsify negative side hits, fill the hit list, compute the
- // energy depostion from negative side for the counter.
-
- Double_t thresh_neg = fNegThresh[hit->fCounter -1];
- thresh_neg=0;
- if(fA_Neg[hit->fCounter-1] > thresh_neg) {
-
-
- fA_Neg_p[hit->fCounter-1] =adc_neg_pedsub ;
-
- fEneg[hit->fCounter-1] += fA_Neg_p[hit->fCounter-1]*
- fParent->GetGain(hit->fCounter-1,fLayerNum-1,1);
- }
-
- // Mean energy in the counter.
-
- fEmean[hit->fCounter-1] += (fEpos[hit->fCounter-1] + fEneg[hit->fCounter-1]);
-
- // Accumulate energies in the plane.
-
- fEplane += fEmean[hit->fCounter-1];
- fEplane_pos += fEpos[hit->fCounter-1];
- fEplane_neg += fEneg[hit->fCounter-1];
-
ihit++;
}
- //Debug output.
+ return(ihit);
+}
+//_____________________________________________________________________________
+Int_t THcShowerPlane::CoarseProcessHits()
+{
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+ Int_t ADCMode=fParent->GetADCMode();
+ if(ADCMode == kADCDynamicPedestal) {
+ FillADC_DynamicPedestal();
+ } else if (ADCMode == kADCSampleIntegral) {
+ FillADC_SampleIntegral();
+ } else if (ADCMode == kADCSampIntDynPed) {
+ FillADC_SampIntDynPed();
+ } else {
+ FillADC_Standard();
+ }
+ //
if (fParent->fdbg_decoded_cal) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerPlane::ProcessHits for "
<< fParent->GetPrefix() << ":" << endl;
- cout << " nrawhits = " << nrawhits << " nexthit = " << nexthit << endl;
cout << " Sparsified hits for HMS calorimeter 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;
- }
+ for (UInt_t i=0; i<fParent->GetNBlocks(fLayerNum-1); i++) {
- if(fA_Pos[hit->fCounter-1] > fPosThresh[hit->fCounter -1] ||
- fA_Neg[hit->fCounter-1] > fNegThresh[hit->fCounter -1]) {
- cout << " plane = " << hit->fPlane
- << " counter = " << hit->fCounter
- << " Emean = " << fEmean[hit->fCounter-1]
- << " Epos = " << fEpos[hit->fCounter-1]
- << " Eneg = " << fEneg[hit->fCounter-1]
+ if (GetAposP(i) > 0 || GetAnegP(i) >0) { //hit
+ cout << " counter = " << i
+ << " Emean = " << fEmean[i]
+ << " Epos = " << fEpos[i]
+ << " Eneg = " << fEneg[i]
<< endl;
nspar++;
}
@@ -559,10 +522,103 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
<< endl;
cout << "---------------------------------------------------------------\n";
}
-
- return(ihit);
+ //
+ return 1;
+ //
+}
+//_____________________________________________________________________________
+void THcShowerPlane::FillADC_SampIntDynPed()
+{
+ // adc_pos = hit->GetRawAdcHitPos().GetSampleInt();
+ // adc_neg = hit->GetRawAdcHitNeg().GetSampleInt();
+ // adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw();
+ // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw();
+ // Need to create
}
+//_____________________________________________________________________________
+void THcShowerPlane::FillADC_SampleIntegral()
+{
+ /// adc_pos_pedsub = hit->GetRawAdcHitPos().GetSampleIntRaw() - fPosPed[hit->fCounter -1];
+ // adc_neg_pedsub = hit->GetRawAdcHitNeg().GetSampleIntRaw() - fNegPed[hit->fCounter -1];
+ // adc_pos = hit->GetRawAdcHitPos().GetSampleIntRaw();
+ // adc_neg = hit->GetRawAdcHitNeg().GetSampleIntRaw();
+ // need to create
+}
+//_____________________________________________________________________________
+void THcShowerPlane::FillADC_Standard()
+{
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+ for (Int_t ielem=0;ielem<frNegAdcPulseIntRaw->GetEntries();ielem++) {
+ Int_t npad = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+ Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+ fA_Neg[npad] =pulseIntRaw;
+ if(fA_Neg[npad] > fNegThresh[npad]) {
+ fA_Neg_p[npad] = pulseIntRaw-fNegPed[npad];
+ fEneg[npad] = fA_Neg_p[npad]*fParent->GetGain(npad,fLayerNum-1,1);
+ fEmean[npad] += fEneg[npad];
+ fEplane_neg += fEneg[npad];
+ }
+ }
+ for (Int_t ielem=0;ielem<frPosAdcPulseIntRaw->GetEntries();ielem++) {
+ Int_t npad = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+ Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+ fA_Pos[npad] =pulseIntRaw;
+ if(fA_Pos[npad] > fPosThresh[npad]) {
+ fA_Pos_p[npad] =pulseIntRaw-fPosPed[npad] ;
+ fEpos[npad] = fA_Pos_p[npad]*fParent->GetGain(npad,fLayerNum-1,0);
+ fEmean[npad] += fEpos[npad];
+ fEplane_pos += fEpos[npad];
+ }
+ }
+ fEplane= fEplane_neg+fEplane_pos;
+}
+//_____________________________________________________________________________
+void THcShowerPlane::FillADC_DynamicPedestal()
+{
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+ Double_t AdcTimeWindowMin=fParent->GetAdcTimeWindowMin();
+ Double_t AdcTimeWindowMax=fParent->GetAdcTimeWindowMax();
+ for (Int_t ielem=0;ielem<frNegAdcPulseInt->GetEntries();ielem++) {
+ Int_t npad = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+ Double_t pulseInt = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
+ Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+ Double_t pulseTime = ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+ Double_t errorflag = ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
+ Bool_t pulseTimeCut = (pulseTime > AdcTimeWindowMin) && (pulseTime < AdcTimeWindowMax);
+ if (errorflag==0 && pulseTimeCut) {
+ fA_Neg[npad] =pulseIntRaw;
+ if(fA_Neg[npad] > fNegThresh[npad]) {
+ fA_Neg_p[npad] =pulseInt ;
+ fEneg[npad] = fA_Neg_p[npad]*fParent->GetGain(npad,fLayerNum-1,1);
+ fEmean[npad] += fEneg[npad];
+ fEplane_neg += fEneg[npad];
+ }
+ }
+ }
+ //
+ for (Int_t ielem=0;ielem<frPosAdcPulseInt->GetEntries();ielem++) {
+ Int_t npad = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+ Double_t pulseInt = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData();
+ Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+ Double_t pulseTime = ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+ Double_t errorflag = ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
+ Bool_t pulseTimeCut = pulseTime > AdcTimeWindowMin && pulseTime < AdcTimeWindowMax;
+ if (errorflag==0 && pulseTimeCut) {
+ fA_Pos[npad] =pulseIntRaw;
+ if(fA_Pos[npad] > fPosThresh[npad]) {
+ fA_Pos_p[npad] =pulseInt ;
+ fEpos[npad] = fA_Pos_p[npad]*fParent->GetGain(npad,fLayerNum-1,0);
+ fEmean[npad] += fEpos[npad];
+ fEplane_pos += fEpos[npad];
+ }
+ }
+ }
+ //
+ fEplane= fEplane_neg+fEplane_pos;
+}
//_____________________________________________________________________________
Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
{
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 8cb4346..fe00766 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -39,6 +39,7 @@ public:
virtual Bool_t IsPid() { return kFALSE; }
virtual Int_t ProcessHits(TClonesArray* rawhits, Int_t nexthit);
+ virtual Int_t CoarseProcessHits();
virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit);
virtual void CalculatePedestals( );
@@ -110,7 +111,6 @@ protected:
// Flash ADC parameters
Int_t fUsingFADC; // != 0 if using FADC in sample mode
- Int_t fADCMode; //
// 1 == Use the pulse int - pulse ped
// 2 == Use the sample integral - known ped
// 3 == Use the sample integral - sample ped
@@ -159,6 +159,7 @@ protected:
Float_t *fNegSig;
Float_t *fNegThresh;
+ TClonesArray* frPosAdcErrorFlag;
TClonesArray* frPosAdcPedRaw;
TClonesArray* frPosAdcPulseIntRaw;
TClonesArray* frPosAdcPulseAmpRaw;
@@ -168,6 +169,7 @@ protected:
TClonesArray* frPosAdcPulseInt;
TClonesArray* frPosAdcPulseAmp;
+ TClonesArray* frNegAdcErrorFlag;
TClonesArray* frNegAdcPedRaw;
TClonesArray* frNegAdcPulseIntRaw;
TClonesArray* frNegAdcPulseAmpRaw;
@@ -180,6 +182,10 @@ protected:
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
virtual void InitializePedestals( );
+ virtual void FillADC_DynamicPedestal( );
+ virtual void FillADC_SampleIntegral( );
+ virtual void FillADC_SampIntDynPed( );
+ virtual void FillADC_Standard( );
ClassDef(THcShowerPlane,0); // Calorimeter bars in a plane
};
#endif
--
GitLab