Skip to content
Snippets Groups Projects
Commit e9f33e22 authored by Vardan Tadevosyan's avatar Vardan Tadevosyan
Browse files

Added per plane energy deposition calculations

parent 7fab8e8f
Branches
Tags
No related merge requests found
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 );
......
......@@ -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
......
......@@ -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;
......
......@@ -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
......
......@@ -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);
......
......@@ -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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment