-
Vardan Tadevosyan authored
Add time windows for pos and neg ADC selection. Added fGoodPosAdcMult and fGoodNegAdcMult variables Remove cuts for calculation of fXAtAero and fYAtAero
Vardan Tadevosyan authoredAdd time windows for pos and neg ADC selection. Added fGoodPosAdcMult and fGoodNegAdcMult variables Remove cuts for calculation of fXAtAero and fYAtAero
THcAerogel.cxx 43.84 KiB
/** \class THcAerogel
\ingroup Detectors
\brief Class for an Aerogel detector consisting of pairs of PMT's
attached to a diffuser box
*/
#include "THcAerogel.h"
#include "THcHodoscope.h"
#include "TClonesArray.h"
#include "THcSignalHit.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THaCutList.h"
#include "THcParmList.h"
#include "THcHitList.h"
#include "THaApparatus.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"
#include "THaTrackProj.h"
#include "THcRawAdcHit.h"
#include "THcHallCSpectrometer.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
using namespace std;
//_____________________________________________________________________________
THcAerogel::THcAerogel( const char* name, const char* description,
THaApparatus* apparatus ) :
THaNonTrackingDetector(name,description,apparatus)
{
InitArrays();
}
//_____________________________________________________________________________
THcAerogel::THcAerogel( ) :
THaNonTrackingDetector()
{
// Constructor
frPosAdcPedRaw = NULL;
frPosAdcPulseIntRaw = NULL;
frPosAdcPulseAmpRaw = NULL;
frPosAdcPulseTimeRaw = NULL;
frPosAdcPed = NULL;
frPosAdcPulseInt = NULL;
frPosAdcPulseAmp = NULL;
frPosAdcPulseTime = NULL;
frNegAdcPedRaw = NULL;
frNegAdcPulseIntRaw = NULL;
frNegAdcPulseAmpRaw = NULL;
frNegAdcPulseTimeRaw = NULL;
frNegAdcPed = NULL;
frNegAdcPulseInt = NULL;
frNegAdcPulseAmp = NULL;
frNegAdcPulseTime = NULL;
fPosAdcErrorFlag = NULL;
fNegAdcErrorFlag = NULL;
// 6 GeV variables fPosTDCHits = NULL;
fNegTDCHits = NULL;
fPosADCHits = NULL;
fNegADCHits = NULL;
InitArrays();
}
//_____________________________________________________________________________
THcAerogel::~THcAerogel()
{
// Destructor
delete frPosAdcPedRaw; frPosAdcPedRaw = NULL;
delete frPosAdcPulseIntRaw; frPosAdcPulseIntRaw = NULL;
delete frPosAdcPulseAmpRaw; frPosAdcPulseAmpRaw = NULL;
delete frPosAdcPulseTimeRaw; frPosAdcPulseTimeRaw = NULL;
delete frPosAdcPed; frPosAdcPed = NULL;
delete frPosAdcPulseInt; frPosAdcPulseInt = NULL;
delete frPosAdcPulseAmp; frPosAdcPulseAmp = NULL;
delete frPosAdcPulseTime; frPosAdcPulseTime = NULL;
delete frNegAdcPedRaw; frNegAdcPedRaw = NULL;
delete frNegAdcPulseIntRaw; frNegAdcPulseIntRaw = NULL;
delete frNegAdcPulseAmpRaw; frNegAdcPulseAmpRaw = NULL;
delete frNegAdcPulseTimeRaw; frNegAdcPulseTimeRaw = NULL;
delete frNegAdcPed; frNegAdcPed = NULL;
delete frNegAdcPulseInt; frNegAdcPulseInt = NULL;
delete frNegAdcPulseAmp; frNegAdcPulseAmp = NULL;
delete frNegAdcPulseTime; frNegAdcPulseTime = NULL;
delete fPosAdcErrorFlag; fPosAdcErrorFlag = NULL;
delete fNegAdcErrorFlag; fNegAdcErrorFlag = NULL;
// 6 GeV variables
delete fPosTDCHits; fPosTDCHits = NULL;
delete fNegTDCHits; fNegTDCHits = NULL;
delete fPosADCHits; fPosADCHits = NULL;
delete fNegADCHits; fNegADCHits = NULL;
DeleteArrays();
}
//_____________________________________________________________________________
void THcAerogel::InitArrays()
{
fPosGain = NULL;
fNegGain = NULL;
// 6 GeV variables
fA_Pos = NULL;
fA_Neg = NULL;
fA_Pos_p = NULL;
fA_Neg_p = NULL;
fT_Pos = NULL;
fT_Neg = NULL;
fPosPedLimit = NULL;
fNegPedLimit = NULL;
fPosPedMean = NULL;
fNegPedMean = NULL;
fPosPedSum = NULL;
fPosPedSum2 = NULL;
fPosPedCount = NULL;
fNegPedSum = NULL;
fNegPedSum2 = NULL;
fNegPedCount = NULL;
fPosPed = NULL;
fPosSig = NULL;
fPosThresh = NULL;
fNegPed = NULL;
fNegSig = NULL; fNegThresh = NULL;
}
//_____________________________________________________________________________
void THcAerogel::DeleteArrays()
{
delete [] fPosGain; fPosGain = NULL;
delete [] fNegGain; fNegGain = NULL;
// 6 GeV variables
delete [] fA_Pos; fA_Pos = NULL;
delete [] fA_Neg; fA_Neg = NULL;
delete [] fA_Pos_p; fA_Pos_p = NULL;
delete [] fA_Neg_p; fA_Neg_p = NULL;
delete [] fT_Pos; fT_Pos = NULL;
delete [] fT_Neg; fT_Neg = NULL;
delete [] fPosPedLimit; fPosPedLimit = NULL;
delete [] fNegPedLimit; fNegPedLimit = NULL;
delete [] fPosPedMean; fPosPedMean = NULL;
delete [] fNegPedMean; fNegPedMean = NULL;
delete [] fPosPedSum; fPosPedSum = NULL;
delete [] fPosPedSum2; fPosPedSum2 = NULL;
delete [] fPosPedCount; fPosPedCount = NULL;
delete [] fNegPedSum; fNegPedSum = NULL;
delete [] fNegPedSum2; fNegPedSum2 = NULL;
delete [] fNegPedCount; fNegPedCount = NULL;
delete [] fPosPed; fPosPed = NULL;
delete [] fPosSig; fPosSig = NULL;
delete [] fPosThresh; fPosThresh = NULL;
delete [] fNegPed; fNegPed = NULL;
delete [] fNegSig; fNegSig = NULL;
delete [] fNegThresh; fNegThresh = NULL;
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcAerogel::Init( const TDatime& date )
{
// cout << "THcAerogel::Init for: " << GetName() << endl;
char EngineDID[] = "xAERO";
EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
static const char* const here = "Init()";
Error( Here(here), "Error filling detectormap for %s.", EngineDID );
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, "THcAerogelHit", 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 THcAerogel::ReadDatabase( const TDatime& date )
{
// This function is called by THaDetectorBase::Init() once at the beginning
// of the analysis.
// cout << "THcAerogel::ReadDatabase for: " << GetName() << endl;
char prefix[2];
prefix[0]=tolower(GetApparatus()->GetName()[0]);
prefix[1]='\0';
CreateMissReportParms(Form("%saero",prefix));
fNRegions = 1; // Default if not set in parameter file
DBRequest listextra[]={
{"aero_num_pairs", &fNelem, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&listextra, prefix);
Bool_t optional = true ;
cout << "Created aerogel detector " << GetApparatus()->GetName() << "."
<< GetName() << " with " << fNelem << " PMT pairs" << endl;
fPosGain = new Double_t[fNelem];
fNegGain = new Double_t[fNelem];
// 6 GeV variables
fTdcOffset = 0; // Offset to make reference time subtracted times positve
fPosPedLimit = new Int_t[fNelem];
fNegPedLimit = new Int_t[fNelem];
fA_Pos = new Float_t[fNelem];
fA_Neg = new Float_t[fNelem];
fA_Pos_p = new Float_t[fNelem];
fA_Neg_p = new Float_t[fNelem];
fT_Pos = new Float_t[fNelem];
fT_Neg = new Float_t[fNelem];
fPosPedMean = new Double_t[fNelem];
fNegPedMean = new Double_t[fNelem];
// Normal constructor with name and description
frPosAdcPedRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseIntRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseAmpRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseTimeRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPed = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseInt = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseAmp = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frPosAdcPulseTime = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPedRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseIntRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseAmpRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseTimeRaw = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPed = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseInt = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseAmp = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
frNegAdcPulseTime = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
fPosAdcErrorFlag = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
fNegAdcErrorFlag = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
fNumPosAdcHits = vector<Int_t> (fNelem, 0.0);
fNumGoodPosAdcHits = vector<Int_t> (fNelem, 0.0);
fNumNegAdcHits = vector<Int_t> (fNelem, 0.0);
fNumGoodNegAdcHits = vector<Int_t> (fNelem, 0.0);
fNumTracksMatched = vector<Int_t> (fNelem, 0.0); fNumTracksFired = vector<Int_t> (fNelem, 0.0);
fPosNpe = vector<Double_t> (fNelem, 0.0);
fNegNpe = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcPed = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcMult = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcPulseInt = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcPulseAmp = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcPulseTime = vector<Double_t> (fNelem, 0.0);
fGoodPosAdcTdcDiffTime = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcPed = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcMult = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcPulseInt = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcPulseIntRaw = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcPulseAmp = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcPulseTime = vector<Double_t> (fNelem, 0.0);
fGoodNegAdcTdcDiffTime = vector<Double_t> (fNelem, 0.0);
// 6 GeV variables
fPosTDCHits = new TClonesArray("THcSignalHit", fNelem*16);
fNegTDCHits = new TClonesArray("THcSignalHit", fNelem*16);
fPosADCHits = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
fNegADCHits = new TClonesArray("THcSignalHit", fNelem*MaxNumAdcPulse);
fPosNpeSixGev = vector<Double_t> (fNelem, 0.0);
fNegNpeSixGev = vector<Double_t> (fNelem, 0.0);
// Create arrays to hold pedestal results
if (fSixGevData) InitializePedestals();
// Region parameters
fRegionsValueMax = fNRegions * 8;
fRegionValue = new Double_t[fRegionsValueMax];
DBRequest list[]={
{"aero_num_regions", &fNRegions, kInt},
{"aero_red_chi2_min", &fRedChi2Min, kDouble},
{"aero_red_chi2_max", &fRedChi2Max, kDouble},
{"aero_beta_min", &fBetaMin, kDouble},
{"aero_beta_max", &fBetaMax, kDouble},
{"aero_enorm_min", &fENormMin, kDouble},
{"aero_enorm_max", &fENormMax, kDouble},
{"aero_dp_min", &fDpMin, kDouble},
{"aero_dp_max", &fDpMax, kDouble},
{"aero_diff_box_zpos", &fDiffBoxZPos, kDouble},
{"aero_npe_thresh", &fNpeThresh, kDouble},
//// {"aero_adcTimeWindowMin", &fAdcTimeWindowMin, kDouble},
//// {"aero_adcTimeWindowMax", &fAdcTimeWindowMax, kDouble},
{"aero_adcPosTimeWindowMin", &fAdcPosTimeWindowMin, kDouble, 0, 1},
{"aero_adcPosTimeWindowMax", &fAdcPosTimeWindowMax, kDouble, 0, 1},
{"aero_adcNegTimeWindowMin", &fAdcNegTimeWindowMin, kDouble, 0, 1},
{"aero_adcNegTimeWindowMax", &fAdcNegTimeWindowMax, kDouble, 0, 1},
{"aero_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
{"aero_debug_adc", &fDebugAdc, kInt, 0, 1},
{"aero_six_gev_data", &fSixGevData, kInt, 0, 1},
{"aero_pos_gain", fPosGain, kDouble, (UInt_t) fNelem},
{"aero_neg_gain", fNegGain, kDouble, (UInt_t) fNelem},
{"aero_pos_ped_limit", fPosPedLimit, kInt, (UInt_t) fNelem, optional},
{"aero_neg_ped_limit", fNegPedLimit, kInt, (UInt_t) fNelem, optional},
{"aero_pos_ped_mean", fPosPedMean, kDouble, (UInt_t) fNelem, optional},
{"aero_neg_ped_mean", fNegPedMean, kDouble, (UInt_t) fNelem, optional},
{"aero_tdc_offset", &fTdcOffset, kInt, 0, optional},
{"aero_min_peds", &fMinPeds, kInt, 0, optional},
{"aero_region", &fRegionValue[0], kDouble, (UInt_t) fRegionsValueMax},
{"aero_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
{0}
};
fAdcPosTimeWindowMin=-1000.;
fAdcNegTimeWindowMin=-1000.;
fAdcPosTimeWindowMax=1000.; fAdcNegTimeWindowMax=1000.;
fSixGevData = 0; // Set 6 GeV data parameter to false unless set in parameter file
fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file
fAdcTdcOffset = 0.0;
fADC_RefTimeCut = 0;
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
if (fSixGevData) cout << "6 GeV Data Analysis Flag Set To TRUE" << endl;
fIsInit = true;
// cout << "Track Matching Parameters for: " << GetApparatus()->GetName()
// << "." << 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;
// }
return kOK;
}
//_____________________________________________________________________________
Int_t THcAerogel::DefineVariables( EMode mode )
{
// Initialize global variables for histogramming and tree
// cout << "THcAerogel::DefineVariables called for: " << GetName() << endl;
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
// Register variables in global list
// Do we need to put the number of pos/neg TDC/ADC hits into the variables?
// No. They show up in tree as Ndata.H.aero.postdchits for example
if (fDebugAdc) {
RVarDef vars[] = {
{"posGain", "Positive PMT gains", "fPosGain"},
{"negGain", "Negative PMT gains", "fNegGain"},
{"numPosAdcHits", "Number of Positive ADC Hits Per PMT", "fNumPosAdcHits"}, // Aerogel occupancy
{"totNumPosAdcHits", "Total Number of Positive ADC Hits", "fTotNumPosAdcHits"}, // Aerogel multiplicity
{"numNegAdcHits", "Number of Negative ADC Hits Per PMT", "fNumNegAdcHits"}, // Aerogel occupancy
{"totNumNegAdcHits", "Total Number of Negative ADC Hits", "fTotNumNegAdcHits"}, // Aerogel multiplicity
{"totnumAdcHits", "Total Number of ADC Hits Per PMT", "fTotNumAdcHits"}, // Aerogel multiplicity
{"posAdcPedRaw", "Positive Raw ADC pedestals", "frPosAdcPedRaw.THcSignalHit.GetData()"},
{"posAdcPulseIntRaw", "Positive Raw ADC pulse integrals", "frPosAdcPulseIntRaw.THcSignalHit.GetData()"},
{"posAdcPulseAmpRaw", "Positive Raw ADC pulse amplitudes", "frPosAdcPulseAmpRaw.THcSignalHit.GetData()"},
{"posAdcPulseTimeRaw", "Positive Raw ADC pulse times", "frPosAdcPulseTimeRaw.THcSignalHit.GetData()"},
{"posAdcPed", "Positive ADC pedestals", "frPosAdcPed.THcSignalHit.GetData()"},
{"posAdcPulseInt", "Positive ADC pulse integrals", "frPosAdcPulseInt.THcSignalHit.GetData()"},
{"posAdcPulseAmp", "Positive ADC pulse amplitudes", "frPosAdcPulseAmp.THcSignalHit.GetData()"},
{"posAdcPulseTime", "Positive ADC pulse times", "frPosAdcPulseTime.THcSignalHit.GetData()"},
{"negAdcPedRaw", "Negative Raw ADC pedestals", "frNegAdcPedRaw.THcSignalHit.GetData()"},
{"negAdcPulseIntRaw", "Negative Raw ADC pulse integrals", "frNegAdcPulseIntRaw.THcSignalHit.GetData()"},
{"negAdcPulseAmpRaw", "Negative Raw ADC pulse amplitudes", "frNegAdcPulseAmpRaw.THcSignalHit.GetData()"},
{"negAdcPulseTimeRaw", "Negative Raw ADC pulse times", "frNegAdcPulseTimeRaw.THcSignalHit.GetData()"},
{"negAdcPed", "Negative ADC pedestals", "frNegAdcPed.THcSignalHit.GetData()"},
{"negAdcPulseInt", "Negative ADC pulse integrals", "frNegAdcPulseInt.THcSignalHit.GetData()"},
{"negAdcPulseAmp", "Negative ADC pulse amplitudes", "frNegAdcPulseAmp.THcSignalHit.GetData()"},
{"negAdcPulseTime", "Negative ADC pulse times", "frNegAdcPulseTime.THcSignalHit.GetData()"},
{ 0 }
};
DefineVarsFromList( vars, mode); } //end debug statement
if (fSixGevData) {
RVarDef vars[] = {
{"apos", "Positive Raw ADC Amplitudes", "fA_Pos"},
{"aneg", "Negative Raw ADC Amplitudes", "fA_Neg"},
{"apos_p", "Positive Ped-subtracted ADC Amplitudes", "fA_Pos_p"},
{"aneg_p", "Negative Ped-subtracted ADC Amplitudes", "fA_Neg_p"},
{"tpos", "Positive Raw TDC", "fT_Pos"},
{"tneg", "Negative Raw TDC", "fT_Neg"},
{"ntdc_pos_hits", "Number of Positive Tube Hits", "fNTDCPosHits"},
{"ntdc_neg_hits", "Number of Negative Tube Hits", "fNTDCNegHits"},
{"posadchits", "Positive ADC hits", "fPosADCHits.THcSignalHit.GetPaddleNumber()"},
{"negadchits", "Negative ADC hits", "fNegADCHits.THcSignalHit.GetPaddleNumber()"},
{"postdchits", "Positive TDC hits", "fPosTDCHits.THcSignalHit.GetPaddleNumber()"},
{"negtdchits", "Negative TDC hits", "fNegTDCHits.THcSignalHit.GetPaddleNumber()"},
{"nGoodHits", "Total number of good hits", "fNGoodHits"},
{"posNpeSixGev", "Number of Positive PEs", "fPosNpeSixGev"},
{"negNpeSixGev", "Number of Negative PEs", "fNegNpeSixGev"},
{"posNpeSumSixGev", "Total Number of Positive PEs", "fPosNpeSumSixGev"},
{"negNpeSumSixGev", "Total Number of Negative PEs", "fNegNpeSumSixGev"},
{"npeSumSixGev", "Total Number of PEs", "fNpeSumSixGev"},
{ 0 }
};
DefineVarsFromList( vars, mode);
} //end fSixGevData statement
RVarDef vars[] = {
{"posAdcCounter", "Positive ADC counter numbers", "frPosAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"negAdcCounter", "Negative ADC counter numbers", "frNegAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
{"posAdcErrorFlag", "Error Flag for When FPGA Fails", "fPosAdcErrorFlag.THcSignalHit.GetData()"},
{"negAdcErrorFlag", "Error Flag for When FPGA Fails", "fNegAdcErrorFlag.THcSignalHit.GetData()"},
{"numGoodPosAdcHits", "Number of Good Positive ADC Hits Per PMT", "fNumGoodPosAdcHits"}, // Aerogel occupancy
{"numGoodNegAdcHits", "Number of Good Negative ADC Hits Per PMT", "fNumGoodNegAdcHits"}, // Aerogel occupancy
{"totNumGoodPosAdcHits", "Total Number of Good Positive ADC Hits", "fTotNumGoodPosAdcHits"}, // Aerogel multiplicity
{"totNumGoodNegAdcHits", "Total Number of Good Negative ADC Hits", "fTotNumGoodNegAdcHits"}, // Aerogel multiplicity
{"totnumGoodAdcHits", "TotalNumber of Good ADC Hits Per PMT", "fTotNumGoodAdcHits"}, // Aerogel 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"},
{"xAtAero", "Track X at Aero diffusion box", "fXAtAero"},
{"yAtAero", "Track Y at Aero diffusion box", "fYAtAero"},
{"posNpe", "Number of Positive PEs", "fPosNpe"},
{"negNpe", "Number of Negative PEs", "fNegNpe"},
{"posNpeSum", "Total Number of Positive PEs", "fPosNpeSum"},
{"negNpeSum", "Total Number of Negative PEs", "fNegNpeSum"},
{"npeSum", "Total Number of PEs", "fNpeSum"},
{"goodPosAdcPed", "Good Negative ADC pedestals", "fGoodPosAdcPed"},
{"goodPosAdcMult", "Good Positive ADC mult", "fGoodPosAdcMult"},
{"goodPosAdcPulseInt", "Good Negative ADC pulse integrals", "fGoodPosAdcPulseInt"},
{"goodPosAdcPulseIntRaw", "Good Negative ADC raw pulse integrals", "fGoodPosAdcPulseIntRaw"},
{"goodPosAdcPulseAmp", "Good Negative ADC pulse amplitudes", "fGoodPosAdcPulseAmp"},
{"goodPosAdcPulseTime", "Good Negative ADC pulse times", "fGoodPosAdcPulseTime"},
{"goodPosAdcTdcDiffTime", "Good Positive hodo Start - ADC pulse times", "fGoodPosAdcTdcDiffTime"},
{"goodNegAdcPed", "Good Negative ADC pedestals", "fGoodNegAdcPed"},
{"goodNegAdcMult", "Good Negative ADC Mult", "fGoodNegAdcMult"},
{"goodNegAdcPulseInt", "Good Negative ADC pulse integrals", "fGoodNegAdcPulseInt"},
{"goodNegAdcPulseIntRaw", "Good Negative ADC raw pulse integrals", "fGoodNegAdcPulseIntRaw"},
{"goodNegAdcPulseAmp", "Good Negative ADC pulse amplitudes", "fGoodNegAdcPulseAmp"},
{"goodNegAdcPulseTime", "Good Negative ADC pulse times", "fGoodNegAdcPulseTime"},
{"goodNegAdcTdcDiffTime", "Good Negative hodo Start - ADC pulse times", "fGoodNegAdcTdcDiffTime"},
{ 0 }
};
return DefineVarsFromList(vars, mode);
}
//_____________________________________________________________________________
inline
void THcAerogel::Clear(Option_t* opt)
{
// Clear the hit lists
fNhits = 0;
fTotNumAdcHits = 0;
fTotNumGoodAdcHits = 0;
fTotNumPosAdcHits = 0;
fTotNumNegAdcHits = 0;
fTotNumGoodPosAdcHits = 0;
fTotNumGoodNegAdcHits = 0;
fTotNumTracksMatched = 0;
fTotNumTracksFired = 0;
fXAtAero = 0.0;
fYAtAero = 0.0;
fNpeSum = 0.0;
fPosNpeSum = 0.0;
fNegNpeSum = 0.0;
frPosAdcPedRaw->Clear();
frPosAdcPulseIntRaw->Clear();
frPosAdcPulseAmpRaw->Clear();
frPosAdcPulseTimeRaw->Clear();
frPosAdcPed->Clear();
frPosAdcPulseInt->Clear();
frPosAdcPulseAmp->Clear();
frPosAdcPulseTime->Clear();
frNegAdcPedRaw->Clear();
frNegAdcPulseIntRaw->Clear();
frNegAdcPulseAmpRaw->Clear();
frNegAdcPulseTimeRaw->Clear();
frNegAdcPed->Clear();
frNegAdcPulseInt->Clear();
frNegAdcPulseAmp->Clear();
frNegAdcPulseTime->Clear();
fPosAdcErrorFlag->Clear();
fNegAdcErrorFlag->Clear();
for (UInt_t ielem = 0; ielem < fNumPosAdcHits.size(); ielem++)
fNumPosAdcHits.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumNegAdcHits.size(); ielem++)
fNumNegAdcHits.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumGoodPosAdcHits.size(); ielem++)
fNumGoodPosAdcHits.at(ielem) = 0;
for (UInt_t ielem = 0; ielem < fNumGoodNegAdcHits.size(); ielem++)
fNumGoodNegAdcHits.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 < fGoodPosAdcPed.size(); ielem++) {
fGoodPosAdcPed.at(ielem) = 0.0;
fGoodPosAdcMult.at(ielem) = 0.0;
fGoodPosAdcPulseInt.at(ielem) = 0.0;
fGoodPosAdcPulseIntRaw.at(ielem) = 0.0;
fGoodPosAdcPulseAmp.at(ielem) = 0.0;
fGoodPosAdcPulseTime.at(ielem) = kBig;
fGoodPosAdcTdcDiffTime.at(ielem) = kBig;
fPosNpe.at(ielem) = 0.0;
}
for (UInt_t ielem = 0; ielem < fGoodNegAdcPed.size(); ielem++) { fGoodNegAdcPed.at(ielem) = 0.0;
fGoodNegAdcMult.at(ielem) = 0.0;
fGoodNegAdcPulseInt.at(ielem) = 0.0;
fGoodNegAdcPulseIntRaw.at(ielem) = 0.0;
fGoodNegAdcPulseAmp.at(ielem) = 0.0;
fGoodNegAdcPulseTime.at(ielem) = kBig;
fGoodNegAdcTdcDiffTime.at(ielem) = kBig;
fNegNpe.at(ielem) = 0.0;
}
// 6 GeV variables
fNGoodHits = 0;
fNADCPosHits = 0;
fNADCNegHits = 0;
fNTDCPosHits = 0;
fNTDCNegHits = 0;
fNpeSumSixGev = 0.0;
fPosNpeSumSixGev = 0.0;
fNegNpeSumSixGev = 0.0;
fPosTDCHits->Clear();
fNegTDCHits->Clear();
fPosADCHits->Clear();
fNegADCHits->Clear();
for (UInt_t ielem = 0; ielem < fPosNpeSixGev.size(); ielem++)
fPosNpeSixGev.at(ielem) = 0.0;
for (UInt_t ielem = 0; ielem < fNegNpeSixGev.size(); ielem++)
fNegNpeSixGev.at(ielem) = 0.0;
for(Int_t itube = 0;itube < fNelem;itube++) {
fA_Pos[itube] = 0;
fA_Neg[itube] = 0;
fA_Pos_p[itube] = 0;
fA_Neg_p[itube] = 0;
fT_Pos[itube] = 0;
fT_Neg[itube] = 0;
}
}
//_____________________________________________________________________________
Int_t THcAerogel::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 (fSixGevData) {
if(gHaCuts->Result("Pedestal_event")) {
AccumulatePedestals(fRawHitList);
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return(0);
}
if(fAnalyzePedestals) {
CalculatePedestals();
Print("");
fAnalyzePedestals = 0; // Don't analyze pedestals next event
}
}
Int_t ihit = 0;
UInt_t nrPosAdcHits = 0;
UInt_t nrNegAdcHits = 0;
while(ihit < fNhits) {
THcAerogelHit* hit = (THcAerogelHit*) fRawHitList->At(ihit);
Int_t npmt = hit->fCounter; THcRawAdcHit& rawPosAdcHit = hit->GetRawAdcHitPos();
THcRawAdcHit& rawNegAdcHit = hit->GetRawAdcHitNeg();
for (UInt_t thit=0; thit<rawPosAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPedRaw());
((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPed());
((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseIntRaw(thit));
((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseInt(thit));
((THcSignalHit*) frPosAdcPulseAmpRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseAmpRaw(thit));
((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseAmp(thit));
((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseTimeRaw(thit));
((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseTime(thit)+fAdcTdcOffset);
if (rawPosAdcHit.GetPulseAmpRaw(thit) > 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 0);
if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 1);
++nrPosAdcHits;
fTotNumAdcHits++;
fTotNumPosAdcHits++;
fNumPosAdcHits.at(npmt-1) = npmt;
}
for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPedRaw());
((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPed());
((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseIntRaw(thit));
((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseInt(thit));
((THcSignalHit*) frNegAdcPulseAmpRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseAmpRaw(thit));
((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseAmp(thit));
((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseTimeRaw(thit));
((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseTime(thit));
if (rawNegAdcHit.GetPulseAmpRaw(thit) > 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 0);
if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 1);
++nrNegAdcHits;
fTotNumAdcHits++;
fTotNumNegAdcHits++;
fNumNegAdcHits.at(npmt-1) = npmt;
}
ihit++;
}
return ihit;
}
//_____________________________________________________________________________
Int_t THcAerogel::ApplyCorrections( void )
{
return(0);
}
//_____________________________________________________________________________
Int_t THcAerogel::CoarseProcess( TClonesArray& ) //tracks
{
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
//cout << " starttime = " << StartTime << endl ;
// Loop over the elements in the TClonesArray
for(Int_t ielem = 0; ielem < frPosAdcPulseInt->GetEntries(); ielem++) {
Int_t npmt = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
Double_t pulsePed = ((THcSignalHit*) frPosAdcPed->ConstructedAt(ielem))->GetData();
Double_t pulseInt = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData(); Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
//// Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcPosTimeWindowMin && adctdcdiffTime < fAdcPosTimeWindowMax;
// By default, the last hit within the timing cut will be considered "good"
if (!errorFlag && pulseTimeCut) {
fGoodPosAdcPed.at(npmt) = pulsePed;
fGoodPosAdcMult.at(npmt) = frPosAdcPulseInt->GetEntries();
// cout << " out = " << npmt << " " << frPosAdcPulseInt->GetEntries() << " " <<fGoodPosAdcMult.at(npmt);
fGoodPosAdcPulseInt.at(npmt) = pulseInt;
fGoodPosAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodPosAdcPulseAmp.at(npmt) = pulseAmp;
fGoodPosAdcPulseTime.at(npmt) = pulseTime;
fGoodPosAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fPosNpe.at(npmt) = fPosGain[npmt]*fGoodPosAdcPulseInt.at(npmt);
fPosNpeSum += fPosNpe.at(npmt);
fTotNumGoodAdcHits++;
fTotNumGoodPosAdcHits++;
fNumGoodPosAdcHits.at(npmt) = npmt + 1;
}
}
// Loop over the elements in the TClonesArray
for(Int_t ielem = 0; ielem < frNegAdcPulseInt->GetEntries(); ielem++) {
Int_t npmt = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
Double_t pulsePed = ((THcSignalHit*) frNegAdcPed->ConstructedAt(ielem))->GetData();
Double_t pulseInt = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorFlag = ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
//// Bool_t pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
Bool_t pulseTimeCut = adctdcdiffTime > fAdcNegTimeWindowMin && adctdcdiffTime < fAdcNegTimeWindowMax;
// By default, the last hit within the timing cut will be considered "good"
if (!errorFlag && pulseTimeCut) {
fGoodNegAdcPed.at(npmt) = pulsePed;
fGoodNegAdcMult.at(npmt) = frNegAdcPulseInt->GetEntries();
fGoodNegAdcPulseInt.at(npmt) = pulseInt;
fGoodNegAdcPulseIntRaw.at(npmt) = pulseIntRaw;
fGoodNegAdcPulseAmp.at(npmt) = pulseAmp;
fGoodNegAdcPulseTime.at(npmt) = pulseTime;
fGoodNegAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
fNegNpe.at(npmt) = fNegGain[npmt]*fGoodNegAdcPulseInt.at(npmt);
fNegNpeSum += fNegNpe.at(npmt);
fTotNumGoodAdcHits++;
fTotNumGoodNegAdcHits++;
fNumGoodNegAdcHits.at(npmt) = npmt + 1;
}
}
fNpeSum = fNegNpeSum + fPosNpeSum;
for(Int_t ihit=0; ihit < fNhits; ihit++) {
Int_t nPosTDCHits = 0;
Int_t nNegTDCHits = 0;
Int_t nPosADCHits = 0;
Int_t nNegADCHits = 0;
// 6 GeV calculations
Int_t adc_pos;
Int_t adc_neg;
Int_t tdc_pos = -1;
Int_t tdc_neg = -1;
if (fSixGevData) {
THcAerogelHit* hit = (THcAerogelHit*) fRawHitList->At(ihit);
Int_t npmt = hit->fCounter - 1;
// Sum positive and negative hits to fill tot_good_hits
if(fPosNpe.at(npmt) > 0.3) {fNADCPosHits++; fNGoodHits++;}
if(fNegNpe.at(npmt) > 0.3) {fNADCNegHits++; fNGoodHits++;}
// ADC positive hit
if((adc_pos = hit->GetRawAdcHitPos().GetPulseInt()) > 0) {
THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
sighit->Set(hit->fCounter, adc_pos);
}
// ADC negative hit
if((adc_neg = hit->GetRawAdcHitNeg().GetPulseInt()) > 0) {
THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
sighit->Set(hit->fCounter, adc_neg);
}
// TDC positive hit
if(hit->GetRawTdcHitPos().GetNHits() > 0) {
THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
tdc_pos = hit->GetRawTdcHitPos().GetTime()+fTdcOffset;
sighit->Set(hit->fCounter, tdc_pos);
}
// TDC negative hit
if(hit->GetRawTdcHitNeg().GetNHits() > 0) {
THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
tdc_neg = hit->GetRawTdcHitNeg().GetTime()+fTdcOffset;
sighit->Set(hit->fCounter, tdc_neg);
}
// For each TDC, identify the first hit that is positive.
tdc_pos = -1;
tdc_neg = -1;
for(UInt_t thit=0; thit<hit->GetRawTdcHitPos().GetNHits(); thit++) {
Int_t tdc = hit->GetRawTdcHitPos().GetTime(thit);
if(tdc >=0 ) {
tdc_pos = tdc;
break;
}
}
for(UInt_t thit=0; thit<hit->GetRawTdcHitNeg().GetNHits(); thit++) {
Int_t tdc = hit->GetRawTdcHitNeg().GetTime(thit);
if(tdc >= 0) {
tdc_neg = tdc;
break;
}
}
fA_Pos[npmt] = adc_pos;
fA_Neg[npmt] = adc_neg;
fA_Pos_p[npmt] = fA_Pos[npmt] - fPosPedMean[npmt];
fA_Neg_p[npmt] = fA_Neg[npmt] - fNegPedMean[npmt];
fT_Pos[npmt] = tdc_pos;
fT_Neg[npmt] = tdc_neg;
if(fA_Pos[npmt] < 8000) fPosNpeSixGev[npmt] = fPosGain[npmt]*fA_Pos_p[npmt];
else fPosNpeSixGev[npmt] = 100.0;
if(fA_Neg[npmt] < 8000) fNegNpeSixGev[npmt] = fNegGain[npmt]*fA_Neg_p[npmt];
else fNegNpeSixGev[npmt] = 100.0;
fPosNpeSumSixGev += fPosNpeSixGev[npmt];
fNegNpeSumSixGev += fNegNpeSixGev[npmt];
// Sum positive and negative hits to fill tot_good_hits if(fPosNpeSixGev[npmt] > 0.3) {fNADCPosHits++; fNGoodHits++;}
if(fNegNpeSixGev[npmt] > 0.3) {fNADCNegHits++; fNGoodHits++;}
if(fT_Pos[npmt] > 0 && fT_Pos[npmt] < 8000) fNTDCPosHits++;
if(fT_Neg[npmt] > 0 && fT_Neg[npmt] < 8000) fNTDCNegHits++;
}
if (fPosNpeSumSixGev > 0.5 || fNegNpeSumSixGev > 0.5)
fNpeSumSixGev = fPosNpeSumSixGev + fNegNpeSumSixGev;
else fNpeSumSixGev = 0.0;
// If total hits are 0, then give a noticable ridiculous NPE
if (fNhits < 1) fNpeSumSixGev = 0.0;
}
return 0;
}
//_____________________________________________________________________________
Int_t THcAerogel::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 trackENorm = trackEnergy/trackMom;
Double_t trackDp = track->GetDp();
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;
fXAtAero = trackXfp + trackTheta * fDiffBoxZPos;
fYAtAero = trackYfp + trackPhi * fDiffBoxZPos;
if (trackRedChi2Cut && trackBetaCut && trackENormCut && trackDpCut) {
// Project the track to the Aerogel diffuser box plane
// cout << "Aerogel Detector: " << GetName() << 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 << "fDiffBoxZPos = " << fDiffBoxZPos << "\t" << "fXAtAero = " << fXAtAero << "\t" << "fYAtAero = " << fYAtAero << endl;
// cout << "=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:" << endl;
for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
if ((TMath::Abs(fRegionValue[GetIndex(iregion, 0)] - fXAtAero) < fRegionValue[GetIndex(iregion, 4)]) &&
(TMath::Abs(fRegionValue[GetIndex(iregion, 1)] - fYAtAero) < 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;
}
//_____________________________________________________________________________
// Method for initializing pedestals in the 6 GeV era
void THcAerogel::InitializePedestals()
{
fNPedestalEvents = 0;
fMinPeds = 0; // Do not calculate pedestals by default
fPosPedSum = new Int_t [fNelem];
fPosPedSum2 = new Int_t [fNelem];
fPosPedLimit = new Int_t [fNelem];
fPosPedCount = new Int_t [fNelem];
fNegPedSum = new Int_t [fNelem];
fNegPedSum2 = new Int_t [fNelem];
fNegPedLimit = new Int_t [fNelem];
fNegPedCount = new Int_t [fNelem];
fPosPed = new Double_t [fNelem];
fNegPed = new Double_t [fNelem];
fPosThresh = new Double_t [fNelem];
fNegThresh = new Double_t [fNelem];
for(Int_t i = 0;i < fNelem; i++) {
fPosPedSum[i] = 0;
fPosPedSum2[i] = 0;
fPosPedLimit[i] = 1000; // In engine, this are set in parameter file
fPosPedCount[i] = 0;
fNegPedSum[i] = 0;
fNegPedSum2[i] = 0;
fNegPedLimit[i] = 1000; // In engine, this are set in parameter file
fNegPedCount[i] = 0;
fPosPedMean[i] = 0; // Default pedestal values
fNegPedMean[i] = 0; // Default pedestal values
}
}
//_____________________________________________________________________________
// Method for accumulating pedestals in the 6 GeV era
void THcAerogel::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) {
THcAerogelHit* hit = (THcAerogelHit *) rawhits->At(ihit);
Int_t element = hit->fCounter - 1;
Int_t adcpos = hit->GetRawAdcHitPos().GetPulseInt();
Int_t adcneg = hit->GetRawAdcHitNeg().GetPulseInt();
if(adcpos <= fPosPedLimit[element]) {
fPosPedSum[element] += adcpos;
fPosPedSum2[element] += adcpos*adcpos;
fPosPedCount[element]++; if(fPosPedCount[element] == fMinPeds/5)
fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
}
if(adcneg <= fNegPedLimit[element]) {
fNegPedSum[element] += adcneg;
fNegPedSum2[element] += adcneg*adcneg;
fNegPedCount[element]++;
if(fNegPedCount[element] == fMinPeds/5)
fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
}
ihit++;
}
fNPedestalEvents++;
return;
}
//_____________________________________________________________________________
// Method for calculating pedestals in the 6 GeV era
void THcAerogel::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++) {
// Positive tubes
fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
fPosThresh[i] = fPosPed[i] + 15;
// Negative tubes
fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
fNegThresh[i] = fNegPed[i] + 15;
// cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
// 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(fPosPedCount[i] > fMinPeds)
fPosPedMean[i] = fPosPed[i];
if(fNegPedCount[i] > fMinPeds)
fNegPedMean[i] = fNegPed[i];
}
}
}
//_____________________________________________________________________________
Int_t THcAerogel::GetIndex(Int_t nRegion, Int_t nValue)
{
return fNRegions * nValue + nRegion;
}
//_____________________________________________________________________________
void THcAerogel::Print(const Option_t* opt) const
{
THaNonTrackingDetector::Print(opt);
// Print out the pedestals
cout << endl;
cout << "Aerogel Pedestals" << endl;
cout << "No. Neg Pos" << endl;
for(Int_t i=0; i<fNelem; i++)
cout << " " << i << "\t" << fNegPedMean[i] << "\t" << fPosPedMean[i] << endl;
cout << endl;
cout << " fMinPeds = " << fMinPeds << endl;
cout << endl;
}
//_____________________________________________________________________________
Int_t THcAerogel::End(THaRunBase* run)
{ MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
return 0;
}
ClassImp(THcAerogel)
////////////////////////////////////////////////////////////////////////////////