/** \class THcCherenkov
\ingroup Detectors
\brief Class for gas Cherenkov detectors
*/
#include "THcCherenkov.h"
#include "TClonesArray.h"
#include "THaApparatus.h"
#include "THaCutList.h"
#include "THaDetMap.h"
#include "THaEvData.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcHodoscope.h"
#include "THcParmList.h"
#include "THcSignalHit.h"
#include "TMath.h"
#include "VarDef.h"
#include "VarType.h"
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <string>
using namespace std;
using std::cin;
using std::cout;
using std::endl;
using std::setprecision;
using std::setw;
//_____________________________________________________________________________
THcCherenkov::THcCherenkov(const char* name, const char* description, THaApparatus* apparatus)
: THaNonTrackingDetector(name, description, apparatus) {
// Normal constructor with name and description
frAdcPedRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseIntRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseAmpRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseTimeRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPed = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseInt = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseAmp = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
frAdcPulseTime = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
fAdcErrorFlag = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
fNumAdcHits = vector<Int_t>(MaxNumCerPmt, 0.0);
fNumGoodAdcHits = vector<Int_t>(MaxNumCerPmt, 0.0);
fNumTracksMatched = vector<Int_t>(MaxNumCerPmt, 0.0);
fNumTracksFired = vector<Int_t>(MaxNumCerPmt, 0.0);
fNpe = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcPed = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcMult = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcHitUsed = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcPulseInt = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcPulseIntRaw = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcPulseAmp = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcPulseTime = vector<Double_t>(MaxNumCerPmt, 0.0);
fGoodAdcTdcDiffTime = vector<Double_t>(MaxNumCerPmt, 0.0);
InitArrays();
}
//_____________________________________________________________________________
THcCherenkov::THcCherenkov() {
// Constructor
frAdcPedRaw = NULL;
frAdcPulseIntRaw = NULL;
frAdcPulseAmpRaw = NULL;
frAdcPulseTimeRaw = NULL;
frAdcPed = NULL;
frAdcPulseInt = NULL;
frAdcPulseAmp = NULL;
frAdcPulseTime = NULL;
fAdcErrorFlag = NULL;
InitArrays();
}
//_____________________________________________________________________________
THcCherenkov::~THcCherenkov() {
// Destructor
delete frAdcPedRaw;
frAdcPedRaw = NULL;
delete frAdcPulseIntRaw;
frAdcPulseIntRaw = NULL;
delete frAdcPulseAmpRaw;
frAdcPulseAmpRaw = NULL;
delete frAdcPulseTimeRaw;
frAdcPulseTimeRaw = NULL;
delete frAdcPed;
frAdcPed = NULL;
delete frAdcPulseInt;
frAdcPulseInt = NULL;
delete frAdcPulseAmp;
frAdcPulseAmp = NULL;
delete frAdcPulseTime;
frAdcPulseTime = NULL;
delete fAdcErrorFlag;
fAdcErrorFlag = NULL;
DeleteArrays();
}
//_____________________________________________________________________________
void THcCherenkov::InitArrays() {
fGain = NULL;
fPedSum = NULL;
fPedSum2 = NULL;
fPedLimit = NULL;
fPedMean = NULL;
fPedCount = NULL;
fPed = NULL;
fThresh = NULL;
}
//_____________________________________________________________________________
void THcCherenkov::DeleteArrays() {
delete[] fGain;
fGain = NULL;
// 6 Gev variables
delete[] fPedSum;
fPedSum = NULL;
delete[] fPedSum2;
fPedSum2 = NULL;
delete[] fPedLimit;
fPedLimit = NULL;
delete[] fPedMean;
fPedMean = NULL;
delete[] fPedCount;
fPedCount = NULL;
delete[] fPed;
fPed = NULL;
delete[] fThresh;
fThresh = NULL;
delete[] fAdcTimeWindowMin;
fAdcTimeWindowMin = 0;
delete[] fAdcTimeWindowMax;
fAdcTimeWindowMax = 0;
delete[] fRegionValue;
fRegionValue = 0;
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcCherenkov::Init(const TDatime& date) {
// cout << "THcCherenkov::Init for: " << GetName() << endl;
string EngineDID = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
std::transform(EngineDID.begin(), EngineDID.end(), EngineDID.begin(), ::toupper);
if (gHcDetectorMap->FillMap(fDetMap, EngineDID.c_str()) < 0) {
static const char* const here = "Init()";
Error(Here(here), "Error filling detectormap for %s.", EngineDID.c_str());
return kInitError;
}
EStatus status;
if ((status = THaNonTrackingDetector::Init(date)))
return fStatus = status;
// Should probably put this in ReadDatabase as we will know the
// maximum number of hits after setting up the detector map
InitHitList(fDetMap, "THcCherenkovHit", fDetMap->GetTotNumChan() + 1, 0, fADC_RefTimeCut);
THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
if (!app || !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector("hod")))) {
static const char* const here = "ReadDatabase()";
Warning(Here(here), "Hodoscope \"%s\" not found. ", "hod");
}
fPresentP = 0;
THaVar* vpresent = gHaVars->Find(Form("%s.present", GetApparatus()->GetName()));
if (vpresent) {
fPresentP = (Bool_t*)vpresent->GetValuePointer();
}
return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcCherenkov::ReadDatabase(const TDatime& date) {
// This function is called by THaDetectorBase::Init() once at the beginning
// of the analysis.
// cout << "THcCherenkov::ReadDatabase for: " << GetName() << endl; // Ahmed
string prefix = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
std::transform(prefix.begin(), prefix.end(), prefix.begin(), ::tolower);
CreateMissReportParms(prefix.c_str());
fNRegions = 4; // Defualt if not set in paramter file
DBRequest list_1[] = {{"_tot_pmts", &fNelem, kInt}, {"_num_regions", &fNRegions, kInt}, {0}};
gHcParms->LoadParmValues(list_1, prefix.c_str());
Bool_t optional = true;
_det_logger->info("Created Cherenkov detector {}.{} with {} PMTs", GetApparatus()->GetName(),
GetName(), fNelem);
// cout << "Created Cherenkov detector " << GetApparatus()->GetName() << "."
// << GetName() << " with " << fNelem << " PMTs" << endl;
// 6 GeV pedestal paramters
fPedLimit = new Int_t[fNelem];
fGain = new Double_t[fNelem];
fPedMean = new Double_t[fNelem];
fAdcTimeWindowMin = new Double_t[fNelem];
fAdcTimeWindowMax = new Double_t[fNelem];
// Region parameters
fRegionsValueMax = fNRegions * 8;
fRegionValue = new Double_t[fRegionsValueMax];
fAdcGoodElem = new Int_t[fNelem];
fAdcPulseAmpTest = new Double_t[fNelem];
DBRequest list[] = {{"_ped_limit", fPedLimit, kInt, (UInt_t)fNelem, optional},
{"_adc_to_npe", fGain, kDouble, (UInt_t)fNelem},
{"_red_chi2_min", &fRedChi2Min, kDouble},
{"_red_chi2_max", &fRedChi2Max, kDouble},
{"_beta_min", &fBetaMin, kDouble},
{"_beta_max", &fBetaMax, kDouble},
{"_enorm_min", &fENormMin, kDouble},
{"_enorm_max", &fENormMax, kDouble},
{"_dp_min", &fDpMin, kDouble},
{"_dp_max", &fDpMax, kDouble},
{"_mirror_zpos", &fMirrorZPos, kDouble},
{"_npe_thresh", &fNpeThresh, kDouble},
{"_debug_adc", &fDebugAdc, kInt, 0, 1},
{"_adcTimeWindowMin", fAdcTimeWindowMin, kDouble, (UInt_t)fNelem, 1},
{"_adcTimeWindowMax", fAdcTimeWindowMax, kDouble, (UInt_t)fNelem, 1},
{"_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
{"_region", &fRegionValue[0], kDouble, (UInt_t)fRegionsValueMax},
{"_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
{0}};
for (Int_t i = 0; i < fNelem; i++) {
fAdcTimeWindowMin[i] = -1000.;
fAdcTimeWindowMax[i] = 1000.;
}
fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file
fAdcTdcOffset = 0.0;
fADC_RefTimeCut = 0;
gHcParms->LoadParmValues((DBRequest*)&list, prefix.c_str());
// if (fDebugAdc) cout << "Cherenkov ADC Debug Flag Set To TRUE" << endl;
fIsInit = true;
// cout << "Track Matching Parameters for: " << GetName() << endl;
// for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
// cout << "Region = " << iregion + 1 << endl;
// for (Int_t ivalue = 0; ivalue < 8; ivalue++)
// cout << fRegionValue[GetIndex(iregion, ivalue)] << " ";
// cout << endl;
// }
// Create arrays to hold pedestal results
InitializePedestals();
return kOK;
}
//_____________________________________________________________________________
Int_t THcCherenkov::DefineVariables(EMode mode) {
// Initialize global variables for histogramming and tree
// cout << "THcCherenkov::DefineVariables called for: " << GetName() << endl;
if (mode == kDefine && fIsSetup)
return kOK;
fIsSetup = (mode == kDefine);
// Register variables in global list
if (fDebugAdc) {
RVarDef vars[] = {
{"numAdcHits", "Number of ADC Hits Per PMT", "fNumAdcHits"}, // Cherenkov occupancy
{"totNumAdcHits", "Total Number of ADC Hits", "fTotNumAdcHits"}, // Cherenkov multiplicity
{"adcPedRaw", "Raw ADC pedestals", "frAdcPedRaw.THcSignalHit.GetData()"},
{"adcPulseIntRaw", "Raw ADC pulse integrals", "frAdcPulseIntRaw.THcSignalHit.GetData()"},
{"adcPulseAmpRaw", "Raw ADC pulse amplitudes", "frAdcPulseAmpRaw.THcSignalHit.GetData()"},
{"adcPulseTimeRaw", "Raw ADC pulse times", "frAdcPulseTimeRaw.THcSignalHit.GetData()"},
{"adcPed", "ADC pedestals", "frAdcPed.THcSignalHit.GetData()"},
{"adcPulseInt", "ADC pulse integrals", "frAdcPulseInt.THcSignalHit.GetData()"},
{"adcPulseAmp", "ADC pulse amplitudes", "frAdcPulseAmp.THcSignalHit.GetData()"},
{"adcPulseTime", "ADC pulse times", "frAdcPulseTime.THcSignalHit.GetData()"},
{0}};
DefineVarsFromList(vars, mode);
} // end debug statement
RVarDef vars[] = {
{"adcCounter", "ADC counter numbers", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"},
{"numGoodAdcHits", "Number of Good ADC Hits Per PMT",
"fNumGoodAdcHits"}, // Cherenkov occupancy
{"totNumGoodAdcHits", "Total Number of Good ADC Hits",
"fTotNumGoodAdcHits"}, // Cherenkov multiplicity
{"numTracksMatched", "Number of Tracks Matched Per Region", "fNumTracksMatched"},
{"numTracksFired", "Number of Tracks that Fired Per Region", "fNumTracksFired"},
{"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"},
{"totNumTracksFired", "Total Number of Tracks that Fired", "fTotNumTracksFired"},
{"xAtCer", "Track X at Cherenkov mirror", "fXAtCer"},
{"yAtCer", "Track Y at Cherenkov mirror", "fYAtCer"},
{"npe", "Number of PEs", "fNpe"},
{"npeSum", "Total Number of PEs", "fNpeSum"},
{"goodAdcPed", "Good ADC pedestals", "fGoodAdcPed"},
{"goodAdcMult", "Good ADC Multiplicity", "fGoodAdcMult"},
{"goodAdcHitUsed", "Good ADC Hit Used", "fGoodAdcHitUsed"},
{"goodAdcPulseInt", "Good ADC pulse integrals", "fGoodAdcPulseInt"},
{"goodAdcPulseIntRaw", "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"},
{"goodAdcPulseAmp", "Good ADC pulse amplitudes", "fGoodAdcPulseAmp"},
{"goodAdcPulseTime", "Good ADC pulse times", "fGoodAdcPulseTime"},
{"goodAdcTdcDiffTime", "Good Hodo Start - ADC pulse times", "fGoodAdcTdcDiffTime"},
{0}};
return DefineVarsFromList(vars, mode);
}
//_____________________________________________________________________________
Int_t THcCherenkov::ManualInitTree(TTree* t) {
// The most direct path to the output tree!!!
std::string app_name = GetApparatus()->GetName();
std::string det_name = GetName();
// if (t) {
// std::string branch_name = (app_name + "_" + det_name + "_data");
// _det_logger->info("THcHodoscope::ManualInitTree : Adding branch, {}, to output tree",
// branch_name); t->Branch(branch_name.c_str(), &_basic_data, 32000, 99);
//}
if (t) {
std::string branch_name = (app_name + "_" + det_name + "_waveforms");
_det_logger->info("THcCherenkov::ManualInitTree : Adding branch, {}, to output tree",
branch_name);
t->Branch(branch_name.c_str(), &_waveforms, 32000, 99);
}
return 0;
}
//_____________________________________________________________________________
void THcCherenkov::Clear(Option_t* opt) {
// Clear the hit lists
fNhits = 0;
fTotNumAdcHits = 0;
fTotNumGoodAdcHits = 0;
fTotNumTracksMatched = 0;
fTotNumTracksFired = 0;
fXAtCer = 0.0;
fYAtCer = 0.0;
fNpeSum = 0.0;
frAdcPedRaw->Clear();
frAdcPulseIntRaw->Clear();
frAdcPulseAmpRaw->Clear();
frAdcPulseTimeRaw->Clear();
frAdcPed->Clear();
frAdcPulseInt->Clear();
frAdcPulseAmp->Clear();
frAdcPulseTime->Clear();
fAdcErrorFlag->Clear();
for (UInt_t ielem = 0; ielem < fNumAdcHits.size(); ielem++)
fNumAdcHits.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumGoodAdcHits.size(); ielem++)
fNumGoodAdcHits.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumTracksMatched.size(); ielem++)
fNumTracksMatched.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumTracksFired.size(); ielem++)
fNumTracksFired.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fGoodAdcPed.size(); ielem++) {
fGoodAdcPed.at(ielem) = 0.0;
fGoodAdcMult.at(ielem) = 0.0;
fGoodAdcHitUsed.at(ielem) = 0.0;
fGoodAdcPulseInt.at(ielem) = 0.0;
fGoodAdcPulseIntRaw.at(ielem) = 0.0;
fGoodAdcPulseAmp.at(ielem) = 0.0;
fGoodAdcPulseTime.at(ielem) = kBig;
fGoodAdcTdcDiffTime.at(ielem) = kBig;
fNpe.at(ielem) = 0.0;
}
for(auto& wf : _waveforms) {
wf.ZeroBuffer();
}
}
//_____________________________________________________________________________
Int_t THcCherenkov::Decode(const THaEvData& evdata) {
// Get the Hall C style hitlist (fRawHitList) for this event
Bool_t present = kTRUE; // Suppress reference time warnings
if (fPresentP) { // if this spectrometer not part of trigger
present = *fPresentP;
}
fNhits = DecodeToHitList(evdata, !present);
if (gHaCuts->Result("Pedestal_event")) {
AccumulatePedestals(fRawHitList);
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return (0);
}
if (fAnalyzePedestals) {
CalculatePedestals();
fAnalyzePedestals = 0; // Don't analyze pedestals next event
}
Int_t ihit = 0;
UInt_t nrAdcHits = 0;
_waveforms.clear();
while (ihit < fNhits) {
THcCherenkovHit* hit = (THcCherenkovHit*)fRawHitList->At(ihit);
Int_t npmt = hit->fCounter;
THcRawAdcHit& rawAdcHit = hit->GetRawAdcHitPos();
_waveforms.push_back({rawAdcHit.GetSampleBuffer()});
for (UInt_t thit = 0; thit < rawAdcHit.GetNPulses(); thit++) {
((THcSignalHit*)frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPedRaw());
((THcSignalHit*)frAdcPed->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPed());
((THcSignalHit*)frAdcPulseIntRaw->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseIntRaw(thit));
((THcSignalHit*)frAdcPulseInt->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseInt(thit));
((THcSignalHit*)frAdcPulseAmpRaw->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseAmpRaw(thit));
((THcSignalHit*)frAdcPulseAmp->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseAmp(thit));
((THcSignalHit*)frAdcPulseTimeRaw->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseTimeRaw(thit));
((THcSignalHit*)frAdcPulseTime->ConstructedAt(nrAdcHits))
->Set(npmt, rawAdcHit.GetPulseTime(thit) + fAdcTdcOffset);
if (rawAdcHit.GetPulseAmpRaw(thit) > 0)
((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 0);
if (rawAdcHit.GetPulseAmpRaw(thit) <= 0)
((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 1);
++nrAdcHits;
fTotNumAdcHits++;
fNumAdcHits.at(npmt - 1) = npmt;
}
ihit++;
}
return ihit;
}
//_____________________________________________________________________________
Int_t THcCherenkov::ApplyCorrections(void) { return (0); }
//_____________________________________________________________________________
Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
for(Int_t ipmt = 0; ipmt < fNelem; ipmt++) {
fAdcPulseAmpTest[ipmt] = -1000.;
fAdcGoodElem[ipmt]=-1;
}
//
for(Int_t ielem = 0; ielem < frAdcPulseInt->GetEntries(); ielem++) {
Int_t npmt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
Double_t pulseTime = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin[npmt] && adctdcdiffTime < fAdcTimeWindowMax[npmt];
if (!errorFlag)
{
fGoodAdcMult.at(npmt) += 1;
}
if (!errorFlag && pulseTimeCut && pulseAmp > fAdcPulseAmpTest[npmt]) {
fAdcGoodElem[npmt]=ielem;
fAdcPulseAmpTest[npmt] = pulseAmp;
}
}
// Loop over the npmt
for(Int_t npmt = 0; npmt < fNelem; npmt++) {
Int_t ielem = fAdcGoodElem[npmt];
if (ielem != -1) {
Double_t pulsePed = ((THcSignalHit*) frAdcPed->ConstructedAt(ielem))->GetData();
Double_t pulseInt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData();
Double_t pulseIntRaw = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
// By default, the last hit within the timing cut will be considered "good"
fGoodAdcPed.at(npmt) = pulsePed;
fGoodAdcHitUsed.at(npmt) = ielem + 1;
fGoodAdcPulseInt.at(npmt) = pulseInt;
fGoodAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodAdcPulseAmp.at(npmt) = pulseAmp;
fGoodAdcPulseTime.at(npmt) = pulseTime;
fGoodAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fNpe.at(npmt) = fGain[npmt] * fGoodAdcPulseInt.at(npmt);
fNpeSum += fNpe.at(npmt);
fTotNumGoodAdcHits++;
fNumGoodAdcHits.at(npmt) = npmt + 1;
}
}
return 0;
}
//_____________________________________________________________________________
Int_t THcCherenkov::FineProcess(TClonesArray& tracks) {
Int_t nTracks = tracks.GetLast() + 1;
for (Int_t itrack = 0; itrack < nTracks; itrack++) {
THaTrack* track = dynamic_cast<THaTrack*>(tracks[itrack]);
if (track->GetIndex() != 0)
continue; // Select the best track
Double_t trackChi2 = track->GetChi2();
Int_t trackNDoF = track->GetNDoF();
Double_t trackRedChi2 = trackChi2 / trackNDoF;
Double_t trackBeta = track->GetBeta();
Double_t trackEnergy = track->GetEnergy();
Double_t trackMom = track->GetP();
Double_t trackDp = track->GetDp();
Double_t trackENorm = trackEnergy / trackMom;
Double_t trackXfp = track->GetX();
Double_t trackYfp = track->GetY();
Double_t trackTheta = track->GetTheta();
Double_t trackPhi = track->GetPhi();
Bool_t trackRedChi2Cut = trackRedChi2 > fRedChi2Min && trackRedChi2 < fRedChi2Max;
Bool_t trackBetaCut = trackBeta > fBetaMin && trackBeta < fBetaMax;
Bool_t trackENormCut = trackENorm > fENormMin && trackENorm < fENormMax;
Bool_t trackDpCut = trackDp > fDpMin && trackDp < fDpMax;
fXAtCer = trackXfp + trackTheta * fMirrorZPos;
fYAtCer = trackYfp + trackPhi * fMirrorZPos;
if (trackRedChi2Cut && trackBetaCut && trackENormCut && trackDpCut) {
// Project the track to the Cherenkov mirror planes
// cout << "Cherenkov Detector: " << GetName() << " has fNRegions = " << fNRegions << endl;
// cout << "nTracks = " << nTracks << "\t" << "trackChi2 = " << trackChi2
// << "\t" << "trackNDof = " << trackNDoF << "\t" << "trackRedChi2 = " <<
// trackRedChi2 << endl; cout << "trackBeta = " << trackBeta << "\t" << "trackEnergy = " <<
// trackEnergy << "\t"
// << "trackMom = " << trackMom << "\t" << "trackENorm = " << trackENorm << endl;
// cout << "trackXfp = " << trackXfp << "\t" << "trackYfp = " << trackYfp << "\t"
// << "trackTheta = " << trackTheta << "\t" << "trackPhi = " << trackPhi << endl;
// cout << "fMirrorZPos = " << fMirrorZPos << "\t" << "fXAtCer = " << fXAtCer << "\t" <<
// "fYAtCer = " << fYAtCer << endl; cout <<
// "=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:" << endl;
for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
if ((TMath::Abs(fRegionValue[GetIndex(iregion, 0)] - fXAtCer) <
fRegionValue[GetIndex(iregion, 4)]) &&
(TMath::Abs(fRegionValue[GetIndex(iregion, 1)] - fYAtCer) <
fRegionValue[GetIndex(iregion, 5)]) &&
(TMath::Abs(fRegionValue[GetIndex(iregion, 2)] - trackTheta) <
fRegionValue[GetIndex(iregion, 6)]) &&
(TMath::Abs(fRegionValue[GetIndex(iregion, 3)] - trackPhi) <
fRegionValue[GetIndex(iregion, 7)])) {
fTotNumTracksMatched++;
fNumTracksMatched.at(iregion) = iregion + 1;
if (fNpeSum > fNpeThresh) {
fTotNumTracksFired++;
fNumTracksFired.at(iregion) = iregion + 1;
} // NPE threshold cut
} // Regional cuts
} // Loop over regions
} // Tracking cuts
} // Track loop
return 0;
}
//_____________________________________________________________________________
void THcCherenkov::InitializePedestals() {
fNPedestalEvents = 0;
fMinPeds = 500; // In engine, this is set in parameter file
fPedSum = new Int_t[fNelem];
fPedSum2 = new Int_t[fNelem];
fPedCount = new Int_t[fNelem];
fPed = new Double_t[fNelem];
fThresh = new Double_t[fNelem];
for (Int_t i = 0; i < fNelem; i++) {
fPedSum[i] = 0;
fPedSum2[i] = 0;
fPedCount[i] = 0;
}
}
//_____________________________________________________________________________
void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits) {
// Extract data from the hit list, accumulating into arrays for
// calculating pedestals
Int_t nrawhits = rawhits->GetLast() + 1;
Int_t ihit = 0;
while (ihit < nrawhits) {
THcCherenkovHit* hit = (THcCherenkovHit*)rawhits->At(ihit);
Int_t element = hit->fCounter - 1;
Int_t nadc = hit->GetRawAdcHitPos().GetPulseIntRaw();
if (nadc <= fPedLimit[element]) {
fPedSum[element] += nadc;
fPedSum2[element] += nadc * nadc;
fPedCount[element]++;
if (fPedCount[element] == fMinPeds / 5) {
fPedLimit[element] = 100 + fPedSum[element] / fPedCount[element];
}
}
ihit++;
}
fNPedestalEvents++;
return;
}
//_____________________________________________________________________________
void THcCherenkov::CalculatePedestals() {
// Use the accumulated pedestal data to calculate pedestals
// Later add check to see if pedestals have drifted ("Danger Will Robinson!")
// cout << "Plane: " << fPlaneNum << endl;
for (Int_t i = 0; i < fNelem; i++) {
// PMT tubes
fPed[i] = ((Double_t)fPedSum[i]) / TMath::Max(1, fPedCount[i]);
fThresh[i] = fPed[i] + 15;
// Just a copy for now, but allow the possibility that fXXXPedMean is set
// in a parameter file and only overwritten if there is a sufficient number of
// pedestal events. (So that pedestals are sensible even if the pedestal events were
// not acquired.)
if (fMinPeds > 0) {
if (fPedCount[i] > fMinPeds) {
fPedMean[i] = fPed[i];
}
}
}
// cout << " " << endl;
}
//_____________________________________________________________________________
Int_t THcCherenkov::GetIndex(Int_t nRegion, Int_t nValue) { return fNRegions * nValue + nRegion; }
//_____________________________________________________________________________
void THcCherenkov::Print(const Option_t* opt) const {
THaNonTrackingDetector::Print(opt);
// Print out the pedestals
cout << endl;
cout << "Cherenkov Pedestals" << endl;
// Ahmed
cout << "No. ADC" << endl;
for (Int_t i = 0; i < fNelem; i++)
cout << " " << i << " " << fPed[i] << endl;
cout << endl;
}
//_____________________________________________________________________________
Double_t THcCherenkov::GetCerNPE() { return fNpeSum; }
//_____________________________________________________________________________
Int_t THcCherenkov::End(THaRunBase* run) {
MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
return 0;
}
ClassImp(THcCherenkov)
////////////////////////////////////////////////////////////////////////////////