Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jlab/hallc/analyzer_software/hcana
  • whit/hcana
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2129 additions and 1373 deletions
src/THcHitList.h
View file @ e974befb
......@@ -73,6 +73,7 @@ protected:
Int_t slot;
Int_t channel;
Int_t reftime;
Int_t refdifftime;
};
std::vector<RefIndexMap> fRefIndexMaps;
// Should this be a sparse list instead in case user
......
src/THcHodoHit.h
View file @ e974befb
......@@ -34,6 +34,9 @@ public:
Double_t GetNegADCpeak() const { return fNegADC_Peak; }
Double_t GetPosADCtime() const { return fPosADC_Time; }
Double_t GetNegADCtime() const { return fNegADC_Time; }
Double_t GetPosADCCorrtime() const { return fPosADC_CorrTime; }
Double_t GetNegADCCorrtime() const { return fNegADC_CorrTime; }
Double_t GetCalcPosition() const { return fCalcPosition; }
Int_t GetPosTDC() const { return fPosTDC; }
Int_t GetNegTDC() const { return fNegTDC; }
Double_t GetPosCorrectedTime() const { return fPosCorrectedTime;}
......@@ -46,8 +49,9 @@ public:
Int_t GetPaddleNumber() const { return fPaddleNumber; }
Double_t GetPaddleCenter() const { return fPaddleCenter; }
void SetCorrectedTimes(Double_t pos, Double_t neg, Double_t) {
void SetCorrectedTimes(Double_t pos, Double_t neg) {
fPosCorrectedTime = pos; fNegCorrectedTime = neg;
fHasCorrectedTimes = kFALSE;
}
void SetCorrectedTimes(Double_t pos, Double_t neg,
Double_t postof, Double_t negtof,
......@@ -75,6 +79,15 @@ public:
void SetNegADCtime( Double_t ptime) {
fNegADC_Time =ptime;
}
void SetPosADCCorrtime( Double_t ptime) {
fPosADC_CorrTime =ptime;
}
void SetNegADCCorrtime( Double_t ptime) {
fNegADC_CorrTime =ptime;
}
void SetCalcPosition( Double_t calcpos) {
fCalcPosition =calcpos;
}
protected:
static const Double_t kBig; //!
......@@ -86,6 +99,9 @@ protected:
Double_t fNegADC_Peak; // ADC peak amplitude
Double_t fPosADC_Time; // ADC time
Double_t fNegADC_Time; // ADC time
Double_t fPosADC_CorrTime; // ADC time
Double_t fNegADC_CorrTime; // ADC time
Double_t fCalcPosition; // Position along paddle calculated by time diff
Int_t fPaddleNumber;
Double_t fPosCorrectedTime; // Pulse height corrected time
......
src/THcHodoscope.cxx
View file @ e974befb
......@@ -242,9 +242,10 @@ Int_t THcHodoscope::ReadDatabase(const TDatime& date) {
fBetaNoTrk = 0.;
fBetaNoTrkChiSq = 0.;
fNPaddle = new UInt_t[fNPlanes];
fFPTime = new Double_t[fNPlanes];
fPlaneCenter = new Double_t[fNPlanes];
fNPaddle = new UInt_t [fNPlanes];
fFPTime = new Double_t [fNPlanes];
fPlaneCenter = new Double_t[fNPlanes];
fPlaneSpacing = new Double_t[fNPlanes];
prefix[0] = tolower(GetApparatus()->GetName()[0]);
......@@ -260,9 +261,13 @@ Int_t THcHodoscope::ReadDatabase(const TDatime& date) {
// GN added
// reading variables from *hodo.param
fMaxScinPerPlane = fNPaddle[0];
for (Int_t i = 1; i < fNPlanes; i++) {
fMaxScinPerPlane = (fMaxScinPerPlane > fNPaddle[i]) ? fMaxScinPerPlane : fNPaddle[i];
fMaxScinPerPlane=fNPaddle[0];
fTotHodScin=fNPaddle[0];
for (Int_t i=1;i<fNPlanes;i++) {
fMaxScinPerPlane=(fMaxScinPerPlane > fNPaddle[i])? fMaxScinPerPlane : fNPaddle[i];
fTotHodScin+=(fNPaddle[i]);
}
// need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
fMaxHodoScin = fMaxScinPerPlane * fNPlanes;
......@@ -316,52 +321,56 @@ Int_t THcHodoscope::ReadDatabase(const TDatime& date) {
fAdcTdcOffset[ip] = 0.0;
}
DBRequest list[] = {
{"cosmicflag", &fCosmicFlag, kInt, 0, 1},
{"NumPlanesBetaCalc", &fNumPlanesBetaCalc, kInt, 0, 1},
{"start_time_center", &fStartTimeCenter, kDouble},
{"start_time_slop", &fStartTimeSlop, kDouble},
{"scin_tdc_to_time", &fScinTdcToTime, kDouble},
{"scin_tdc_min", &fScinTdcMin, kDouble},
{"scin_tdc_max", &fScinTdcMax, kDouble},
{"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
{"pathlength_central", &fPathLengthCentral, kDouble},
{"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1},
{"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1},
{"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1},
{"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1},
{"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
{"xloscin", &fxLoScin[0], kInt, (UInt_t)fNHodoscopes},
{"xhiscin", &fxHiScin[0], kInt, (UInt_t)fNHodoscopes},
{"yloscin", &fyLoScin[0], kInt, (UInt_t)fNHodoscopes},
{"yhiscin", &fyHiScin[0], kInt, (UInt_t)fNHodoscopes},
{"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
{"trackeff_scint_ydiff_max", &trackeff_scint_ydiff_max, kDouble, 0, 1},
{"trackeff_scint_xdiff_max", &trackeff_scint_xdiff_max, kDouble, 0, 1},
{"cer_npe", &fNCerNPE, kDouble, 0, 1},
{"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
{"hodo_slop", fHodoSlop, kDouble, (UInt_t)fNPlanes},
{"debugprintscinraw", &fdebugprintscinraw, kInt, 0, 1},
{"hodo_tdc_offset", fTdcOffset, kInt, (UInt_t)fNPlanes, 1},
{"hodo_adc_tdc_offset", fAdcTdcOffset, kDouble, (UInt_t)fNPlanes, 1},
{"hodo_PosAdcTimeWindowMin", fHodoPosAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
{"hodo_PosAdcTimeWindowMax", fHodoPosAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
{"hodo_NegAdcTimeWindowMin", fHodoNegAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
{"hodo_NegAdcTimeWindowMax", fHodoNegAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
{"dumptof", &fDumpTOF, kInt, 0, 1},
{"TOFCalib_shtrk_lo", &fTOFCalib_shtrk_lo, kDouble, 0, 1},
{"TOFCalib_shtrk_hi", &fTOFCalib_shtrk_hi, kDouble, 0, 1},
{"TOFCalib_cer_lo", &fTOFCalib_cer_lo, kDouble, 0, 1},
{"TOFCalib_beta_lo", &fTOFCalib_beta_lo, kDouble, 0, 1},
{"TOFCalib_beta_hi", &fTOFCalib_beta_hi, kDouble, 0, 1},
{"dumptof_filename", &fTOFDumpFile, kString, 0, 1},
{0}};
DBRequest list[]={
{"cosmicflag", &fCosmicFlag, kInt, 0, 1},
{"NumPlanesBetaCalc", &fNumPlanesBetaCalc, kInt, 0, 1},
{"start_time_center", &fStartTimeCenter, kDouble},
{"start_time_slop", &fStartTimeSlop, kDouble},
{"scin_tdc_to_time", &fScinTdcToTime, kDouble},
{"scin_tdc_min", &fScinTdcMin, kDouble},
{"scin_tdc_max", &fScinTdcMax, kDouble},
{"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
{"pathlength_central", &fPathLengthCentral, kDouble},
{"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1},
{"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1},
{"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1},
{"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1},
{"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
{"xloscin", &fxLoScin[0], kInt, (UInt_t) fNHodoscopes},
{"xhiscin", &fxHiScin[0], kInt, (UInt_t) fNHodoscopes},
{"yloscin", &fyLoScin[0], kInt, (UInt_t) fNHodoscopes},
{"yhiscin", &fyHiScin[0], kInt, (UInt_t) fNHodoscopes},
{"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
{"trackeff_scint_ydiff_max", &trackeff_scint_ydiff_max, kDouble, 0, 1},
{"trackeff_scint_xdiff_max", &trackeff_scint_xdiff_max, kDouble, 0, 1},
{"cer_npe", &fNCerNPE, kDouble, 0, 1},
{"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
{"hodo_slop", fHodoSlop, kDouble, (UInt_t) fNPlanes},
{"debugprintscinraw", &fdebugprintscinraw, kInt, 0,1},
{"hodo_tdc_offset", fTdcOffset, kInt, (UInt_t) fNPlanes, 1},
{"hodo_adc_tdc_offset", fAdcTdcOffset, kDouble, (UInt_t) fNPlanes, 1},
{"hodo_PosAdcTimeWindowMin", fHodoPosAdcTimeWindowMin, kDouble, (UInt_t) fMaxHodoScin, 1},
{"hodo_PosAdcTimeWindowMax", fHodoPosAdcTimeWindowMax, kDouble, (UInt_t) fMaxHodoScin, 1},
{"hodo_NegAdcTimeWindowMin", fHodoNegAdcTimeWindowMin, kDouble, (UInt_t) fMaxHodoScin, 1},
{"hodo_NegAdcTimeWindowMax", fHodoNegAdcTimeWindowMax, kDouble, (UInt_t) fMaxHodoScin, 1},
{"dumptof", &fDumpTOF, kInt, 0, 1},
{"TOFCalib_shtrk_lo", &fTOFCalib_shtrk_lo, kDouble, 0, 1},
{"TOFCalib_shtrk_hi", &fTOFCalib_shtrk_hi, kDouble, 0, 1},
{"TOFCalib_cer_lo", &fTOFCalib_cer_lo, kDouble, 0, 1},
{"TOFCalib_beta_lo", &fTOFCalib_beta_lo, kDouble, 0, 1},
{"TOFCalib_beta_hi", &fTOFCalib_beta_hi, kDouble, 0, 1},
{"dumptof_filename", &fTOFDumpFile, kString, 0, 1},
{"TrackBetaIncludeSinglePmtHits", &fTrackBetaIncludeSinglePmtHits, kInt, 0, 1},
{0}
};
// Defaults if not defined in parameter file
fTrackBetaIncludeSinglePmtHits=0; // do not use paddles with only one hit in the TRack Beta calculation set ==1 to include
trackeff_scint_ydiff_max=20.;
trackeff_scint_xdiff_max=20.;
for(UInt_t ip=0;ip<fMaxHodoScin;ip++) {
trackeff_scint_ydiff_max = 20.;
trackeff_scint_xdiff_max = 20.;
for (UInt_t ip = 0; ip < fMaxHodoScin; ip++) {
fHodoPosAdcTimeWindowMin[ip] = -1000.;
fHodoPosAdcTimeWindowMax[ip] = 1000.;
fHodoNegAdcTimeWindowMin[ip] = -1000.;
......@@ -549,20 +558,26 @@ Int_t THcHodoscope::DefineVariables(EMode mode) {
// Register variables in global list
RVarDef vars[] = {
// Move these into THcHallCSpectrometer using track fTracks
{"beta", "Beta including track info", "fBeta"},
{"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
{"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
{"fpHitsTime", "Time at focal plane from all hits", "fFPTimeAll"},
{"starttime", "Hodoscope Start Time", "fStartTime"},
{"goodstarttime", "Hodoscope Good Start Time (logical flag)", "fGoodStartTime"},
{"goodscinhit", "Hit in fid area", "fGoodScinHits"},
{"TimeHist_Sigma", "", "fTimeHist_Sigma"},
{"TimeHist_Peak", "", "fTimeHist_Peak"},
{"TimeHist_Hits", "", "fTimeHist_Hits"},
{0}};
return DefineVarsFromList(vars, mode);
// Move these into THcHallCSpectrometer using track fTracks
{"beta", "Beta including track info", "fBeta"},
{"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
{"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
{"fpHitsTime", "Time at focal plane from all hits", "fFPTimeAll"},
{"starttime", "Hodoscope Start Time", "fStartTime"},
{"goodstarttime", "Hodoscope Good Start Time (logical flag)", "fGoodStartTime"},
{"goodscinhit", "Hit in fid area", "fGoodScinHits"},
{"adctdc_offset"," ","fOffsetTime"},
{"TimeHist_StartTime_Sigma", "", "fTimeHist_StartTime_Sigma"},
{"TimeHist_StartTime_Peak", "", "fTimeHist_StartTime_Peak"},
{"TimeHist_StartTime_NumPeaks", "", "fTimeHist_StartTime_NumPeaks"},
{"TimeHist_StartTime_Hits", "", "fTimeHist_StartTime_Hits"},
{"TimeHist_FpTime_Sigma", "", "fTimeHist_FpTime_Sigma"},
{"TimeHist_FpTime_Peak", "", "fTimeHist_FpTime_Peak"},
{"TimeHist_FpTime_NumPeaks", "", "fTimeHist_FpTime_NumPeaks"},
{"TimeHist_FpTime_Hits", "", "fTimeHist_FpTime_Hits"},
{ 0 }
};
return DefineVarsFromList( vars, mode );
// return kOK;
}
//_____________________________________________________________________________
......@@ -705,20 +720,31 @@ void THcHodoscope::Clear(Option_t* opt) {
* Called by THcHodoscope::Decode
*
*/
fTimeHist_Sigma = kBig;
fTimeHist_Peak = kBig;
fTimeHist_Hits = kBig;
fTimeHist_StartTime_Sigma= kBig;
fTimeHist_StartTime_Peak= kBig;
fTimeHist_StartTime_NumPeaks= 0;
fTimeHist_StartTime_Hits= kBig;
fTimeHist_FpTime_Sigma= kBig;
fTimeHist_FpTime_Peak= kBig;
fTimeHist_FpTime_NumPeaks= 0;
fTimeHist_FpTime_Hits= kBig;
fBeta = 0.0;
fBetaNoTrk = 0.0;
fBetaNoTrkChiSq = 0.0;
fStartTime = -1000.;
fFPTimeAll = -1000.;
fGoodStartTime = kFALSE;
fGoodScinHits = 0;
if (*opt != 'I') {
for (Int_t ip = 0; ip < fNPlanes; ip++) {
fStartTime = -1000.;
fADCStartTime = -1000.;
fOffsetTime = kBig;
fFPTimeAll= -1000.;
fGoodStartTime = kFALSE;
fGoodScinHits = 0;
if( *opt != 'I' ) {
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlanes[ip]->Clear();
fFPTime[ip] = 0.;
fPlaneCenter[ip] = 0.;
......@@ -733,6 +759,10 @@ void THcHodoscope::Clear(Option_t* opt) {
fNClust.clear();
fClustSize.clear();
fClustPos.clear();
fNCluster.clear();
fClusterSize.clear();
fClusterXPos.clear();
fClusterYPos.clear();
fThreeScin.clear();
fGoodScinHitsX.clear();
fGoodFlags.clear();
......@@ -760,11 +790,11 @@ Int_t THcHodoscope::Decode(const THaEvData& evdata) {
present = *fPresentP;
}
fNHits = DecodeToHitList(evdata, !present);
fEventNum = evdata.GetEvNum();
//
// GN: print event number so we can cross-check with engine
// if (evdata.GetEvNum()>1000)
// cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;
// cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;
fCheckEvent = evdata.GetEvNum();
fEventType = evdata.GetEvType();
......@@ -789,9 +819,12 @@ Int_t THcHodoscope::Decode(const THaEvData& evdata) {
// Let each plane get its hits
Int_t nexthit = 0;
fNfptimes = 0;
//THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
// cout << " event number = " << fEventNum << " Evtyp = " << fEventType<< " spec = " << app->GetName() << endl;
fNfptimes=0;
Int_t thits = 0;
for (Int_t ip = 0; ip < fNPlanes; ip++) {
for(Int_t ip=0;ip<fNPlanes;ip++) {
fPlaneCenter[ip] = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
fPlaneSpacing[ip] = fPlanes[ip]->GetSpacing();
......@@ -802,9 +835,12 @@ Int_t THcHodoscope::Decode(const THaEvData& evdata) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
thits += fPlanes[ip]->GetNScinHits();
}
fStartTime = -1000;
if (thits > 0)
EstimateFocalPlaneTime();
//
//
fStartTime=-1000;
if (thits>0 ) EstimateFocalPlaneTime();
if (fdebugprintscinraw == 1) {
for (UInt_t ihit = 0; ihit < fNRawHits; ihit++) {
......@@ -818,6 +854,129 @@ Int_t THcHodoscope::Decode(const THaEvData& evdata) {
return fNHits;
}
//_____________________________________________________________________________
Double_t THcHodoscope::DetermineTimePeak(Int_t FillFlag)
{
Double_t time_peak=-1000;
Int_t NBinsX=hTime->GetNbinsX();
Int_t hTimeScanRange = 10.; // Integrate over HtimeScanRange
vector<Double_t> hpeakCent;
vector<Double_t> hpeakNum;
vector<Double_t> hpeakRMS;
vector<Double_t> hpeakFlag;
vector<Int_t> hpeakBin;
Double_t MinimumNum=2.;
Double_t test_peakmax=0.;
Bool_t scanning_for_local_peak=kFALSE;
Double_t save_mean=0,save_rms=0,save_num=0;
Int_t save_bin;
Bool_t new_peak=kFALSE;
Bool_t replace_peak=kFALSE;
UInt_t best_peak_index=0;
Int_t best_peak_num=-1;
Double_t best_peak_diff=1000;
Int_t nfound=0;
for (Int_t nb=1;nb<NBinsX-hTimeScanRange;nb++) {
hTime->GetXaxis()->SetRange(nb,nb+hTimeScanRange);
Double_t test_int = hTime->Integral();
if (scanning_for_local_peak) {
if ( test_int <= test_peakmax) {
Int_t ps=hpeakCent.size();
replace_peak=kFALSE;
new_peak=kFALSE;
if (ps==0) new_peak=kTRUE;
if (ps!=0 && nb==hpeakBin[ps]+1 && save_mean!=hpeakCent[ps-1]) new_peak=kFALSE;
if (ps!=0 && save_num > hpeakNum[ps-1] && abs(save_mean-hpeakCent[ps-1])<5) replace_peak=kTRUE;
if (ps!=0 && nb!=hpeakBin[ps]+1 && save_num > MinimumNum && abs(save_mean-hpeakCent[ps-1])>=5) new_peak=kTRUE;
if (new_peak) {
hpeakCent.push_back(save_mean);
hpeakRMS.push_back(save_rms);
hpeakNum.push_back(save_num);
hpeakBin.push_back(save_bin);
hpeakFlag.push_back(1);
}
if (replace_peak) {
hpeakCent[ps-1]=save_mean;
hpeakRMS[ps-1]=save_rms;
hpeakNum[ps-1]=save_num;
hpeakBin[ps-1]=save_bin;
hpeakFlag[ps-1]=1;
}
scanning_for_local_peak = kFALSE;
test_peakmax = 0;
best_peak_index=-1;
best_peak_num=5;
nfound=0;
for (UInt_t np=0;np<hpeakNum.size();np++) {
hpeakFlag[np]=-1;
if ( hpeakNum[np] > 5 && (hpeakNum[np]>= best_peak_num || abs(hpeakNum[np] - best_peak_num)<= 4) ) {
if (nfound==0 || (hpeakNum[np]== best_peak_num || abs(hpeakNum[np] - best_peak_num)<= 4) ) {
hpeakFlag[np]=1;
if (nfound==0 ) best_peak_num =hpeakNum[np] ;
if (nfound==0 ) best_peak_index= np;
nfound++;
} else {
for (UInt_t nt=0;nt<np;nt++) {hpeakFlag[nt]=-1;}
nfound=1;
best_peak_num =hpeakNum[np] ;
best_peak_index= np;
}
}
}
if (nfound>1) {
best_peak_diff=1000;
for (UInt_t np=0;np<hpeakNum.size();np++) {
if (hpeakFlag[np]==1 && abs(hpeakCent[np]-fStartTimeCenter)<best_peak_diff) {
best_peak_diff = abs(hpeakCent[np]-fStartTimeCenter);
best_peak_index= np;
}
}}
if (nfound==0) {
best_peak_diff=1000;
for (UInt_t np=0;np<hpeakNum.size();np++) {
if (abs(hpeakCent[np]-fStartTimeCenter)<best_peak_diff) {
best_peak_diff = abs(hpeakCent[np]-fStartTimeCenter);
best_peak_index= np;
}
}}
} else {
test_peakmax = test_int;
save_mean=hTime->GetMean();
save_rms=hTime->GetRMS();
save_num=hTime->Integral();
save_bin=nb;
}
} else {
if ( test_int > MinimumNum) {
test_peakmax = test_int;
scanning_for_local_peak = kTRUE;
save_mean=hTime->GetMean();
save_rms=hTime->GetRMS();
save_num=hTime->Integral();
save_bin=nb;
}
}
}
//
if (hpeakNum.size() >0 && best_peak_index<hpeakNum.size() ) {
time_peak= hpeakCent[best_peak_index];
if (FillFlag==1) {
fTimeHist_StartTime_NumPeaks=hpeakNum.size() ;
fTimeHist_StartTime_Peak= time_peak;
fTimeHist_StartTime_Sigma= hpeakRMS[best_peak_index] ;
fTimeHist_StartTime_Hits= hpeakNum[best_peak_index];
}
if (FillFlag==2) {
fTimeHist_FpTime_NumPeaks=hpeakNum.size() ;
fTimeHist_FpTime_Peak= time_peak;
fTimeHist_FpTime_Sigma= hpeakRMS[best_peak_index] ;
fTimeHist_FpTime_Hits= hpeakNum[best_peak_index];
}
}
//
return time_peak;
}
//_____________________________________________________________________________
void THcHodoscope::EstimateFocalPlaneTime() {
......@@ -835,93 +994,160 @@ void THcHodoscope::EstimateFocalPlaneTime() {
Int_t nscinhits = 0; // Total # hits with at least one good tdc
hTime->Reset();
//
for (Int_t ip = 0; ip < fNPlanes; ip++) {
Int_t nphits = fPlanes[ip]->GetNScinHits();
//
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
nscinhits += nphits;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for (Int_t i = 0; i < nphits; i++) {
THcHodoHit* hit = (THcHodoHit*)hodoHits->At(i);
if (hit->GetHasCorrectedTimes()) {
Double_t postime = hit->GetPosTOFCorrectedTime();
Double_t negtime = hit->GetNegTOFCorrectedTime();
hTime->Fill(postime);
hTime->Fill(negtime);
for(Int_t i=0;i<nphits;i++) {
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
if(hit->GetHasCorrectedTimes()) {
Double_t postime=hit->GetPosTOFCorrectedTime();
Double_t negtime=hit->GetNegTOFCorrectedTime();
hTime->Fill(postime);
hTime->Fill(negtime);
}
}
}
}
//
Double_t TimePeak=DetermineTimePeak(1);
hTime->Reset();
//
ihit = 0;
Double_t fpTimeSum = 0.0;
fNfptimes = 0;
Int_t Ngood_hits_plane = 0;
Double_t Plane_fptime_sum = 0.0;
Bool_t goodplanetime[fNPlanes];
Bool_t twogoodtimes[nscinhits];
Double_t tmin = 0.5 * hTime->GetMaximumBin();
fTimeHist_Peak = tmin;
fTimeHist_Sigma = hTime->GetRMS();
fTimeHist_Hits = hTime->Integral();
_basic_data.fTimeHist_Peak = fTimeHist_Peak;
_basic_data.fTimeHist_Sigma = fTimeHist_Sigma;
_basic_data.fTimeHist_Hits = fTimeHist_Hits;
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
goodplanetime[ip] = kFALSE;
Int_t nphits = fPlanes[ip]->GetNScinHits();
Double_t AdcTdcDiffTimeSum=0;
Double_t NAdcTdcDiffTimeSum=0;
//
for(Int_t ip=0;ip<fNPlanes;ip++) {
Int_t nphits=fPlanes[ip]->GetNScinHits();
nscinhits += nphits;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
for(Int_t i=0;i<nphits;i++) {
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(i);
if(hit->GetHasCorrectedTimes()) {
NAdcTdcDiffTimeSum++;
AdcTdcDiffTimeSum+=(hit->GetPosADCtime()-hit->GetPosTDC()*fScinTdcToTime);
NAdcTdcDiffTimeSum++;
AdcTdcDiffTimeSum+=(hit->GetNegADCtime()-hit->GetNegTDC()*fScinTdcToTime);
Double_t postime=hit->GetPosADCCorrtime();
Double_t negtime=hit->GetNegADCCorrtime();
hTime->Fill(postime);
hTime->Fill(negtime);
}
}
}
if (NAdcTdcDiffTimeSum>0) AdcTdcDiffTimeSum=AdcTdcDiffTimeSum/NAdcTdcDiffTimeSum;
//
Double_t AdcTimePeak=DetermineTimePeak(3);
//
ihit = 0;
Double_t fpTimeSum = 0.0;
Double_t adcfpTimeSum = 0.0;
Double_t adcNfptimes=0;
fNfptimes=0;
Int_t Ngood_hits_plane=0;
Int_t Ngood_adchits_plane=0;
Double_t Plane_fptime_sum=0.0;
Double_t Plane_adcfptime_sum=0.0;
Bool_t goodplanetime[fNPlanes];
Bool_t twogoodtimes[nscinhits];
Int_t NumPlanesGoodHit=0;
Int_t NumPlanesGoodAdcHit=0;
if (TimePeak>0) {
for(Int_t ip=0;ip<fNumPlanesBetaCalc;ip++) {
goodplanetime[ip] = kFALSE;
Int_t nphits=fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Ngood_hits_plane = 0;
Plane_fptime_sum = 0.0;
for (Int_t i = 0; i < nphits; i++) {
THcHodoHit* hit = (THcHodoHit*)hodoHits->At(i);
twogoodtimes[ihit] = kFALSE;
if (hit->GetHasCorrectedTimes()) {
Double_t postime = hit->GetPosTOFCorrectedTime();
Double_t negtime = hit->GetNegTOFCorrectedTime();
if ((postime > (tmin - fTofTolerance)) && (postime < (tmin + fTofTolerance)) &&
(negtime > (tmin - fTofTolerance)) && (negtime < (tmin + fTofTolerance))) {
hit->SetTwoGoodTimes(kTRUE);
twogoodtimes[ihit] = kTRUE; // Both tubes fired
Int_t index = hit->GetPaddleNumber() - 1; //
Double_t fptime;
if (fCosmicFlag == 1) {
fptime = hit->GetScinCorrectedTime() +
(fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
(29.979 * fBetaNominal);
} else {
fptime = hit->GetScinCorrectedTime() -
(fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
(29.979 * fBetaNominal);
}
if(hit->GetHasCorrectedTimes()) {
Double_t postime=hit->GetPosTOFCorrectedTime();
Double_t negtime=hit->GetNegTOFCorrectedTime();
Double_t adcpostime=hit->GetPosADCCorrtime();
Double_t adcnegtime=hit->GetNegADCCorrtime();
if ((postime>(TimePeak-fTofTolerance)) && (postime<(TimePeak+fTofTolerance)) &&
(negtime>(TimePeak-fTofTolerance)) && (negtime<(TimePeak+fTofTolerance)) ) {
hit->SetTwoGoodTimes(kTRUE);
twogoodtimes[ihit] = kTRUE; // Both tubes fired
Int_t index=hit->GetPaddleNumber()-1; //
Double_t fptime;
if(fCosmicFlag==1) {
fptime = hit->GetScinCorrectedTime()
+ (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
/ (29.979 * fBetaNominal);
}else{
fptime = hit->GetScinCorrectedTime()
- (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
/ (29.979 * fBetaNominal);
}
Ngood_hits_plane++;
Plane_fptime_sum += fptime;
fpTimeSum += fptime;
fNfptimes++;
goodplanetime[ip] = kTRUE;
} else {
hit->SetTwoGoodTimes(kFALSE);
}
Plane_fptime_sum+=fptime;
fpTimeSum += fptime;
fNfptimes++;
goodplanetime[ip] = kTRUE;
} else {
hit->SetTwoGoodTimes(kFALSE);
}
//
if ((adcpostime>(AdcTimePeak-fTofTolerance)) && (adcpostime<(AdcTimePeak+fTofTolerance)) &&
(adcnegtime>(AdcTimePeak-fTofTolerance)) && (adcnegtime<(AdcTimePeak+fTofTolerance)) ) {
Int_t index=hit->GetPaddleNumber()-1; //
Double_t fptime;
if(fCosmicFlag==1) {
fptime = hit->GetScinCorrectedTime()
+ (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
/ (29.979 * fBetaNominal);
}else{
fptime = hit->GetScinCorrectedTime()
- (fPlanes[ip]->GetZpos()+(index%2)*fPlanes[ip]->GetDzpos())
/ (29.979 * fBetaNominal);
}
Ngood_adchits_plane++;
Plane_adcfptime_sum+=fptime;
adcfpTimeSum += fptime;
adcNfptimes++;
}
//
}
ihit++;
}
if (Ngood_hits_plane)
fPlanes[ip]->SetFpTime(Plane_fptime_sum / float(Ngood_hits_plane));
if (Ngood_hits_plane>0) NumPlanesGoodHit++;
if (Ngood_adchits_plane>0) NumPlanesGoodAdcHit++;
if (Ngood_hits_plane>0) fPlanes[ip]->SetFpTime(Plane_fptime_sum/float(Ngood_hits_plane));
fPlanes[ip]->SetNGoodHits(Ngood_hits_plane);
}
if (fNfptimes > 0) {
fStartTime = fpTimeSum / fNfptimes;
fGoodStartTime = kTRUE;
fFPTimeAll = fStartTime;
} // if TimePeak>0
//
if(NumPlanesGoodHit>=3) {
fStartTime = fpTimeSum/fNfptimes;
fGoodStartTime=kTRUE;
fFPTimeAll = fStartTime ;
fOffsetTime=0;
if(NumPlanesGoodAdcHit>=3) {
fADCStartTime = adcfpTimeSum/adcNfptimes-fStartTime;
}
fOffsetTime =AdcTdcDiffTimeSum;
} else {
fStartTime = fStartTimeCenter;
fGoodStartTime = kFALSE;
fFPTimeAll = fStartTime;
fStartTime = fStartTimeCenter;
fADCStartTime = fStartTimeCenter;
fGoodStartTime=kFALSE;
fFPTimeAll = fStartTime ;
fOffsetTime=AdcTdcDiffTimeSum;
}
//
//
hTime->Reset();
//
if ((goodplanetime[0] || goodplanetime[1]) && (goodplanetime[2] || goodplanetime[3])) {
if(fGoodStartTime && (goodplanetime[0]||goodplanetime[1]) &&(goodplanetime[2]||goodplanetime[3])) {
Double_t sumW = 0.;
Double_t sumT = 0.;
......@@ -1026,6 +1252,8 @@ void THcHodoscope::EstimateFocalPlaneTime() {
fGoodEventTOFCalib = kTRUE;
//
//
} else {
fBetaNoTrkChiSq = -10.; // Flag if does not try to find beta
}
}
......@@ -1080,7 +1308,9 @@ Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
Double_t sumFPTime = 0.; // Line 138
fNScinHit.push_back(0);
//! Calculate all corrected hit times and histogram
//! Calculate all corrected hit times and histogram
//! This uses a copy of code below. Results are save in time_pos,neg
//! including the z-pos. correction assuming nominal value of betap
//! Code is currently hard-wired to look for a peak in the
......@@ -1093,12 +1323,14 @@ Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
//! Default value in case user hasnt defined something reasonable
// Loop over scintillator planes.
// In ENGINE, its loop over good scintillator hits.
hTime->Reset();
fTOFCalc.clear(); // SAW - Can we
fTOFPInfo.clear(); // SAW - combine these two?
Int_t ihhit = 0; // Hit # overall
// In ENGINE, its loop over good scintillator hits. (comments from old version)
hTime->Reset();
fTOFCalc.clear(); // SAW - Can we
fTOFPInfo.clear(); // SAW - combine these two?
Int_t ihhit = 0; // Hit # overall
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
std::vector<GoodFlags> goodflagstmp2;
......@@ -1108,146 +1340,159 @@ Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
#else
fGoodFlags[itrack].push_back(goodflagstmp2);
#endif
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Double_t zPos = fPlanes[ip]->GetZpos();
Double_t dzPos = fPlanes[ip]->GetDzpos();
// first loop over hits with in a single plane
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++) {
// iphit is hit # within a plane
THcHodoHit* hit = (THcHodoHit*)hodoHits->At(iphit);
fTOFPInfo.push_back(TOFPInfo());
// Can remove these as we will initialize in the constructor
// fTOFPInfo[ihhit].time_pos = -99.0;
// fTOFPInfo[ihhit].time_neg = -99.0;
// fTOFPInfo[ihhit].keep_pos = kFALSE;
// fTOFPInfo[ihhit].keep_neg = kFALSE;
fTOFPInfo[ihhit].scin_pos_time = 0.0;
fTOFPInfo[ihhit].scin_neg_time = 0.0;
fTOFPInfo[ihhit].hit = hit;
fTOFPInfo[ihhit].planeIndex = ip;
fTOFPInfo[ihhit].hitNumInPlane = iphit;
fTOFPInfo[ihhit].onTrack = kFALSE;
Int_t paddle = hit->GetPaddleNumber() - 1;
Double_t zposition = zPos + (paddle % 2) * dzPos;
Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() * (zposition); // Line 183
Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() * (zposition); // Line 184
Double_t scinTrnsCoord, scinLongCoord;
if ((ip == 0) || (ip == 2)) { // !x plane. Line 185
scinTrnsCoord = xHitCoord;
scinLongCoord = yHitCoord;
} else if ((ip == 1) || (ip == 3)) { // !y plane. Line 188
scinTrnsCoord = yHitCoord;
scinLongCoord = xHitCoord;
} else {
return -1;
} // Line 195
fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;
Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// Index to access the 2d arrays of paddle/scintillator properties
Int_t fPIndex = GetScinIndex(ip, paddle);
Double_t betatrack = theTrack->GetP() / TMath::Sqrt(theTrack->GetP() * theTrack->GetP() +
fPartMass * fPartMass);
if (TMath::Abs(scinCenter - scinTrnsCoord) <
(fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop())) { // Line 293
fTOFPInfo[ihhit].onTrack = kTRUE;
Double_t zcor = zposition / (29.979 * betatrack) *
TMath::Sqrt(1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi());
fTOFPInfo[ihhit].zcor = zcor;
if (fCosmicFlag) {
Double_t zcor = -zposition / (29.979 * 1.0) *
TMath::Sqrt(1. + theTrack->GetTheta() * theTrack->GetTheta() +
theTrack->GetPhi() * theTrack->GetPhi());
fTOFPInfo[ihhit].zcor = zcor;
}
Double_t tdc_pos = hit->GetPosTDC();
if (tdc_pos >= fScinTdcMin && tdc_pos <= fScinTdcMax) {
Double_t adc_pos = hit->GetPosADC();
Double_t adcamp_pos = hit->GetPosADCpeak();
Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
fTOFPInfo[ihhit].pathp = pathp;
Double_t timep = tdc_pos * fScinTdcToTime;
if (fTofUsingInvAdc) {
timep -= fHodoPosInvAdcOffset[fPIndex] + pathp / fHodoPosInvAdcLinear[fPIndex] +
fHodoPosInvAdcAdc[fPIndex] / TMath::Sqrt(TMath::Max(20.0 * .020, adc_pos));
} else {
// Double_t tw_corr_pos =
// fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) -
// fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
Double_t tw_corr_pos = 0.;
pathp = scinLongCoord;
if (adcamp_pos > 0)
tw_corr_pos = 1. / pow(adcamp_pos / fTdc_Thrs, fHodoPos_c2[fPIndex]) -
1. / pow(200. / fTdc_Thrs, fHodoPos_c2[fPIndex]);
timep += -tw_corr_pos + fHodo_LCoeff[fPIndex] + pathp / fHodoVelFit[fPIndex];
}
fTOFPInfo[ihhit].scin_pos_time = timep;
timep -= zcor;
fTOFPInfo[ihhit].time_pos = timep;
fNScinHits[ip] = fPlanes[ip]->GetNScinHits();
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Double_t zPos = fPlanes[ip]->GetZpos();
Double_t dzPos = fPlanes[ip]->GetDzpos();
// first loop over hits with in a single plane
for (Int_t iphit = 0; iphit < fNScinHits[ip]; iphit++ ){
// iphit is hit # within a plane
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
fTOFPInfo.push_back(TOFPInfo());
// Can remove these as we will initialize in the constructor
// fTOFPInfo[ihhit].time_pos = -99.0;
// fTOFPInfo[ihhit].time_neg = -99.0;
// fTOFPInfo[ihhit].keep_pos = kFALSE;
// fTOFPInfo[ihhit].keep_neg = kFALSE;
fTOFPInfo[ihhit].scin_pos_time = 0.0;
fTOFPInfo[ihhit].scin_neg_time = 0.0;
fTOFPInfo[ihhit].hit = hit;
fTOFPInfo[ihhit].planeIndex = ip;
fTOFPInfo[ihhit].hitNumInPlane = iphit;
fTOFPInfo[ihhit].onTrack = kFALSE;
Int_t paddle = hit->GetPaddleNumber()-1;
Double_t zposition = zPos + (paddle%2)*dzPos;
Double_t xHitCoord = theTrack->GetX() + theTrack->GetTheta() *
( zposition ); // Line 183
Double_t yHitCoord = theTrack->GetY() + theTrack->GetPhi() *
( zposition ); // Line 184
Double_t scinTrnsCoord, scinLongCoord;
if ( ( ip == 0 ) || ( ip == 2 ) ){ // !x plane. Line 185
scinTrnsCoord = xHitCoord;
scinLongCoord = yHitCoord;
} else if ( ( ip == 1 ) || ( ip == 3 ) ){ // !y plane. Line 188
scinTrnsCoord = yHitCoord;
scinLongCoord = xHitCoord;
} else { return -1; } // Line 195
fTOFPInfo[ihhit].scinTrnsCoord = scinTrnsCoord;
fTOFPInfo[ihhit].scinLongCoord = scinLongCoord;
Double_t scinCenter = fPlanes[ip]->GetPosCenter(paddle) + fPlanes[ip]->GetPosOffset();
// Index to access the 2d arrays of paddle/scintillator properties
Int_t fPIndex = GetScinIndex(ip,paddle);
Double_t betatrack = theTrack->GetP()/TMath::Sqrt(theTrack->GetP()*theTrack->GetP()+fPartMass*fPartMass);
if ( TMath::Abs( scinCenter - scinTrnsCoord ) <
( fPlanes[ip]->GetSize() * 0.5 + fPlanes[ip]->GetHodoSlop() ) ){ // Line 293
fTOFPInfo[ihhit].onTrack = kTRUE;
Double_t zcor = zposition/(29.979*betatrack)*
TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
+ theTrack->GetPhi()*theTrack->GetPhi());
fTOFPInfo[ihhit].zcor = zcor;
if (fCosmicFlag) {
Double_t zcor = -zposition/(29.979*1.0)*
TMath::Sqrt(1. + theTrack->GetTheta()*theTrack->GetTheta()
+ theTrack->GetPhi()*theTrack->GetPhi());
fTOFPInfo[ihhit].zcor = zcor;
}
Double_t tdc_pos = hit->GetPosTDC();
Double_t tdc_neg = hit->GetNegTDC();
//
if( (tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax) &&
(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax )) {
fTOFPInfo[ihhit].scin_pos_time = hit->GetPosCorrectedTime();
Double_t timep = hit->GetPosCorrectedTime()-zcor;
fTOFPInfo[ihhit].time_pos = timep;
hTime->Fill(timep);
}
Double_t tdc_neg = hit->GetNegTDC();
if (tdc_neg >= fScinTdcMin && tdc_neg <= fScinTdcMax) {
Double_t adc_neg = hit->GetNegADC();
Double_t adcamp_neg = hit->GetNegADCpeak();
Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
fTOFPInfo[ihhit].pathn = pathn;
Double_t timen = tdc_neg * fScinTdcToTime;
if (fTofUsingInvAdc) {
timen -= fHodoNegInvAdcOffset[fPIndex] + pathn / fHodoNegInvAdcLinear[fPIndex] +
fHodoNegInvAdcAdc[fPIndex] / TMath::Sqrt(TMath::Max(20.0 * .020, adc_neg));
} else {
pathn = scinLongCoord;
Double_t tw_corr_neg = 0;
if (adcamp_neg > 0)
tw_corr_neg = 1. / pow(adcamp_neg / fTdc_Thrs, fHodoNeg_c2[fPIndex]) -
1. / pow(200. / fTdc_Thrs, fHodoNeg_c2[fPIndex]);
timen += -tw_corr_neg - 2 * fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex] -
pathn / fHodoVelFit[fPIndex];
}
fTOFPInfo[ihhit].scin_neg_time = timen;
timen -= zcor;
fTOFPInfo[ihhit].time_neg = timen;
fTOFPInfo[ihhit].scin_neg_time = hit->GetNegCorrectedTime();
Double_t timen = hit->GetNegCorrectedTime()-zcor;
fTOFPInfo[ihhit].time_neg = timen;
hTime->Fill(timen);
}
} // condition for cenetr on a paddle
ihhit++;
} // First loop over hits in a plane <---------
} else {
//
if (fTrackBetaIncludeSinglePmtHits==1) {
if(tdc_pos >=fScinTdcMin && tdc_pos <= fScinTdcMax ) {
Double_t adc_pos = hit->GetPosADC();
Double_t adcamp_pos = hit->GetPosADCpeak();
Double_t pathp = fPlanes[ip]->GetPosLeft() - scinLongCoord;
fTOFPInfo[ihhit].pathp = pathp;
Double_t timep = tdc_pos*fScinTdcToTime;
if(fTofUsingInvAdc) {
timep -= fHodoPosInvAdcOffset[fPIndex]
+ pathp/fHodoPosInvAdcLinear[fPIndex]
+ fHodoPosInvAdcAdc[fPIndex]
/TMath::Sqrt(TMath::Max(20.0*.020,adc_pos));
} else {
//Double_t tw_corr_pos = fHodoPos_c1[fPIndex]/pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - fHodoPos_c1[fPIndex]/pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
Double_t tw_corr_pos=0.;
pathp=scinLongCoord;
if (adcamp_pos>0) tw_corr_pos = 1./pow(adcamp_pos/fTdc_Thrs,fHodoPos_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoPos_c2[fPIndex]);
timep += -tw_corr_pos + fHodo_LCoeff[fPIndex]+ pathp/fHodoVelFit[fPIndex];
}
fTOFPInfo[ihhit].scin_pos_time = timep;
timep -= zcor;
fTOFPInfo[ihhit].time_pos = timep;
//-----------------------------------------------------------------------------------------------
//------------- First large loop over scintillator hits ends here --------------------
//-----------------------------------------------------------------------------------------------
hTime->Fill(timep);
}
if(tdc_neg >=fScinTdcMin && tdc_neg <= fScinTdcMax ) {
Double_t adc_neg = hit->GetNegADC();
Double_t adcamp_neg = hit->GetNegADCpeak();
Double_t pathn = scinLongCoord - fPlanes[ip]->GetPosRight();
fTOFPInfo[ihhit].pathn = pathn;
Double_t timen = tdc_neg*fScinTdcToTime;
if(fTofUsingInvAdc) {
timen -= fHodoNegInvAdcOffset[fPIndex]
+ pathn/fHodoNegInvAdcLinear[fPIndex]
+ fHodoNegInvAdcAdc[fPIndex]
/TMath::Sqrt(TMath::Max(20.0*.020,adc_neg));
} else {
pathn=scinLongCoord ;
Double_t tw_corr_neg =0 ;
if (adcamp_neg >0) tw_corr_neg= 1./pow(adcamp_neg/fTdc_Thrs,fHodoNeg_c2[fPIndex]) - 1./pow(200./fTdc_Thrs, fHodoNeg_c2[fPIndex]);
timen += -tw_corr_neg- 2*fHodoCableFit[fPIndex] + fHodo_LCoeff[fPIndex]- pathn/fHodoVelFit[fPIndex];
}
fTOFPInfo[ihhit].scin_neg_time = timen;
timen -= zcor;
fTOFPInfo[ihhit].time_neg = timen;
hTime->Fill(timen);
}
} // new fTrackBetaIncludeSinglePmtHits
} // matches else
} // condition for cenetr on a paddle
ihhit++;
} // First loop over hits in a plane <---------
//-----------------------------------------------------------------------------------------------
//------------- First large loop over scintillator hits ends here --------------------
//-----------------------------------------------------------------------------------------------
}
Int_t nhits = ihhit;
Int_t nhits=ihhit;
Double_t TimePeak = DetermineTimePeak(2);
if(TimePeak> 0) {
for(Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
if ( ( fTOFPInfo[ih].time_pos > (TimePeak-fTofTolerance) ) && ( fTOFPInfo[ih].time_pos < ( TimePeak + fTofTolerance ) ) ) {
fTOFPInfo[ih].keep_pos=kTRUE;
}
if ( ( fTOFPInfo[ih].time_neg > (TimePeak-fTofTolerance) ) && ( fTOFPInfo[ih].time_neg < ( TimePeak + fTofTolerance ) ) ){
fTOFPInfo[ih].keep_neg=kTRUE;
}
}
if (0.5 * hTime->GetMaximumBin() > 0) {
Double_t tmin = 0.5 * hTime->GetMaximumBin();
for (Int_t ih = 0; ih < nhits; ih++) { // loop over all scintillator hits
if ((fTOFPInfo[ih].time_pos > (tmin - fTofTolerance)) &&
(fTOFPInfo[ih].time_pos < (tmin + fTofTolerance))) {
fTOFPInfo[ih].keep_pos = kTRUE;
}
if ((fTOFPInfo[ih].time_neg > (tmin - fTofTolerance)) &&
(fTOFPInfo[ih].time_neg < (tmin + fTofTolerance))) {
fTOFPInfo[ih].keep_neg = kTRUE;
}
}
}
//---------------------------------------------------------------------------------------------
......@@ -1491,34 +1736,38 @@ Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
//
// ---------------------------------------------------------------------------
Double_t FPTimeSum = 0.0;
Int_t nFPTimeSum = 0;
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
if (fNPlaneTime[ip] != 0) {
fFPTime[ip] = (fSumPlaneTime[ip] / fNPlaneTime[ip]);
FPTimeSum += fSumPlaneTime[ip];
nFPTimeSum += fNPlaneTime[ip];
} else {
fFPTime[ip] = 1000. * (ip + 1);
}
Double_t FPTimeSum=0.0;
Int_t nFPTimeSum=0;
Int_t nGoodPlanesHit=0;
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
if ( fNPlaneTime[ip] != 0 ){
nGoodPlanesHit++;
fFPTime[ip] = ( fSumPlaneTime[ip] / fNPlaneTime[ip] );
FPTimeSum += fSumPlaneTime[ip];
nFPTimeSum += fNPlaneTime[ip];
} else {
fFPTime[ip] = 1000. * ( ip + 1 );
}
}
Double_t fptime = -1000;
if (nFPTimeSum > 0)
fptime = FPTimeSum / nFPTimeSum;
fFPTimeAll = fptime;
Double_t dedx = 0.0;
for (UInt_t ih = 0; ih < fTOFCalc.size(); ih++) {
if (fTOFCalc[ih].good_scin_time) {
dedx = fTOFCalc[ih].dedx;
break;
}
Double_t fptime=-2000;
fptime=fStartTime;
if (nGoodPlanesHit>=3) fptime = FPTimeSum/nFPTimeSum;
fFPTimeAll = fptime;
Double_t dedx=0.0;
for(UInt_t ih=0;ih<fTOFCalc.size();ih++) {
if(fTOFCalc[ih].good_scin_time) {
dedx = fTOFCalc[ih].dedx;
break;
}
}
theTrack->SetDedx(dedx);
theTrack->SetFPTime(fptime);
theTrack->SetBeta(beta);
theTrack->SetBetaChi2(betaChiSq);
theTrack->SetNPMT(nPmtHit[itrack]);
theTrack->SetFPTime(timeAtFP[itrack]);
} // Main loop over tracks ends here.
......@@ -1527,10 +1776,156 @@ Int_t THcHodoscope::CoarseProcess(TClonesArray& tracks) {
// OriginalTrackEffTest();
TrackEffTest();
//
CalcCluster();
return 0;
}
void THcHodoscope::TrackEffTest(void) {
//
void THcHodoscope::CalcCluster(void)
{
// THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
// cout << app->GetName() << endl;
const Int_t MaxNCluster=5;
std::vector<Int_t > iw(MaxNCluster,0);
std::vector<Double_t > dw(MaxNCluster,0);
for(Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ) {
fNCluster.push_back(0);
fClusterSize.push_back(iw);
fClusterXPos.push_back(dw);
fClusterYPos.push_back(dw);
}
for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++ ){
Double_t pl_xypos=0;
Double_t pl_calcpos=0;
Double_t pl_zpos=0;
Int_t num_good_pad=0;
Double_t pl_x=0,pl_y=0;
TClonesArray* hodoHits = fPlanes[ip]->GetHits();
Int_t prev_padnum=-100;
for (Int_t iphit = 0; iphit < fPlanes[ip]->GetNScinHits(); iphit++ ){
THcHodoHit *hit = (THcHodoHit*)hodoHits->At(iphit);
if ( hit->GetTwoGoodTimes() ) {
Int_t padind = hit->GetPaddleNumber()-1;
Int_t padnum = padind+1;
if (ip==0 || ip==2) pl_x = fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
if (ip==0 || ip==2) pl_y = ((THcHodoHit*)hodoHits->At(iphit))->GetCalcPosition();
if (ip==1 || ip==3) pl_y = fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
if (ip==1 || ip==3) pl_x = ((THcHodoHit*)hodoHits->At(iphit))->GetCalcPosition();
pl_xypos+=fPlanes[ip]->GetPosCenter(padind)+ fPlanes[ip]->GetPosOffset();
pl_calcpos=((THcHodoHit*)hodoHits->At(iphit))->GetCalcPosition();
pl_zpos+=fPlanes[ip]->GetZpos()+ (padind%2)*fPlanes[ip]->GetDzpos();
num_good_pad++;
if ( fNCluster[ip]>0 && abs(padnum-prev_padnum)==1 && fClusterSize[ip][fNCluster[ip]-1]==1) {
fClusterSize[ip][fNCluster[ip]-1]=fClusterSize[ip][fNCluster[ip]-1]+1;
fClusterXPos[ip][fNCluster[ip]-1]+=pl_x;
fClusterYPos[ip][fNCluster[ip]-1]+=pl_y;
// cout << "Add to cluster pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus = " << fNCluster[ip] << " cl size = " << fClusterSize[ip][fNCluster[ip]-1] << " Xpos " << pl_x << " Ypos = " << pl_y << " postime = " << hit->GetPosTOFCorrectedTime() << " negtime = " << hit->GetNegTOFCorrectedTime() << endl;
} else {
if (fNCluster[ip]<MaxNCluster) fNCluster[ip]++;
fClusterSize[ip][fNCluster[ip]-1]=1;
fClusterXPos[ip][fNCluster[ip]-1]=pl_x;
fClusterYPos[ip][fNCluster[ip]-1]=pl_y;
// cout << " New clus pl = " << ip+1 << " hit = " << iphit << " pad = " << padnum << " clus = " << fNCluster[ip] << " cl size = " << fClusterSize[ip][fNCluster[ip]-1] << " Xpos = " << pl_x << " Ypos = " << pl_y << " postime = " << hit->GetPosTOFCorrectedTime() << " negtime = " << hit->GetNegTOFCorrectedTime() << endl;
}
prev_padnum=padnum;
}
}
//
for (Int_t ic = 0; ic < fNCluster[ip]; ic++ ){
fClusterXPos[ip][ic]/=fClusterSize[ip][ic];
fClusterYPos[ip][ic]/=fClusterSize[ip][ic];
// cout << " Cluster = " << ic+1 << " Xpos = " << fClusterXPos[ip][ic] << " Ypos = " << fClusterYPos[ip][ic] << endl;
}
//
if (num_good_pad !=0 ) {
pl_xypos=pl_xypos/num_good_pad;
pl_calcpos=pl_calcpos/num_good_pad;
pl_zpos=pl_zpos/num_good_pad;
} else {
pl_xypos = kBig;
pl_calcpos = kBig;
pl_zpos = kBig;
}
// fPlanes[ip]->SetTranversePos(pl_xypos);
//fPlanes[ip]->SetLongPos(pl_calcpos);
}
//
// analyze clusters
Int_t best_cluster[4]={-1,-1,-1,-1};
Double_t diffx_test;
Double_t diffx;
Double_t diffy_test;
Double_t diffy;
Int_t pl1,pl2;
for (Int_t nch = 0; nch < 2; nch++ ){
if (nch==0) pl1=0;
if (nch==1) pl1=2;
pl2=pl1+1;
diffx_test=1000;
diffy_test=1000;
if ( fNCluster[pl1]>=1 && fNCluster[pl2]>=1 ) {
for (Int_t ic1 = 0; ic1 < fNCluster[pl1]; ic1++ ){
for (Int_t ic2 = 0; ic2 < fNCluster[pl2]; ic2++ ){
diffx= abs(fClusterXPos[pl1][ic1]-fClusterXPos[pl2][ic2]);
diffy= abs(fClusterYPos[pl1][ic1]-fClusterYPos[pl2][ic2]);
if ( (ic1==0 && ic2==0) || (diffx <=diffx_test && diffy <=diffy_test)) {
diffx_test=diffx;
diffy_test=diffy;
best_cluster[pl1]=ic1;
best_cluster[pl2]=ic2;
}
}}
} else {
if (fNCluster[pl1]==1) best_cluster[pl1]=0;
if (fNCluster[pl2]==1) best_cluster[pl2]=0;
}
}
//
Int_t pl_test1[4]={0,1,2,3};
Int_t pl_test2[4]={2,3,0,1};
for (Int_t npl = 0; npl < 4; npl++ ){
pl1=pl_test1[npl];
pl2=pl_test2[npl];
if (fNCluster[pl1]>0 && best_cluster[pl1]==-1 && fNCluster[pl2]>0 && best_cluster[pl2]>-1) {
diffx_test=1000;
diffy_test=1000;
for (Int_t ic1 = 0; ic1 < fNCluster[pl1]; ic1++ ){
diffx= abs(fClusterXPos[pl1][ic1]-fClusterXPos[pl2][best_cluster[pl2]]);
diffy= abs(fClusterYPos[pl1][ic1]-fClusterYPos[pl2][best_cluster[pl2]]);
if ( (diffx <=diffx_test && diffy <=diffy_test)) {
diffx_test=diffx;
diffy_test=diffy;
best_cluster[pl1]=ic1;
}
}
}
}
//
//
//
for (Int_t npl = 0; npl < 4; npl++ ){
/*
if (best_cluster[npl]==-1) {
cout << " PLane = " << npl+1 << " no best cluster " << endl;
} else {
cout << " plane = " << npl+1 << " xpos = " << fClusterXPos[npl][best_cluster[npl]] << " ypos = " << fClusterYPos[npl][best_cluster[npl]] << endl;
}
*/
if (best_cluster[npl]!=-1) fPlanes[npl]->SetScinYPos( fClusterYPos[npl][best_cluster[npl]] );
if (best_cluster[npl]!=-1) fPlanes[npl]->SetScinXPos( fClusterXPos[npl][best_cluster[npl]] );
}
//
}
//
void THcHodoscope::TrackEffTest(void)
{
//Double_t PadLow[4];
//Double_t PadHigh[4];
// assume X planes are 0,2 and Y planes are 1,3
std::array<int, 4> PadLow = {fxLoScin[0], fyLoScin[0], fxLoScin[1], fyLoScin[1]};
std::array<int, 4> PadHigh = {fxHiScin[0], fyHiScin[0], fxHiScin[1], fyHiScin[1]};
......
src/THcHodoscope.h
View file @ e974befb
......@@ -76,11 +76,14 @@ public:
virtual Int_t FineProcess( TClonesArray& tracks );
virtual Int_t End(THaRunBase* run=0);
Double_t DetermineTimePeak(Int_t FillFlag);
void EstimateFocalPlaneTime(void);
void OriginalTrackEffTest(void);
void TrackEffTest(void);
void CalcCluster(void);
virtual Int_t ApplyCorrections( void );
Double_t GetStartTime() const { return fStartTime; }
Double_t GetOffsetTime() const { return fOffsetTime; }
Bool_t IsStartTimeGood() const {return fGoodStartTime;};
Int_t GetNfptimes() const {return fNfptimes;};
Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle);
......@@ -192,12 +195,19 @@ protected:
Bool_t fSHMS;
Bool_t fGoodStartTime;
Double_t fStartTime;
Double_t fADCStartTime;
Double_t fOffsetTime;
Double_t fFPTimeAll;
Int_t fNfptimes;
Bool_t* fPresentP;
Double_t fTimeHist_Peak;
Double_t fTimeHist_Sigma;
Double_t fTimeHist_Hits;
Double_t fTimeHist_StartTime_NumPeaks;
Double_t fTimeHist_StartTime_Peak;
Double_t fTimeHist_StartTime_Sigma;
Double_t fTimeHist_StartTime_Hits;
Double_t fTimeHist_FpTime_NumPeaks;
Double_t fTimeHist_FpTime_Peak;
Double_t fTimeHist_FpTime_Sigma;
Double_t fTimeHist_FpTime_Hits;
Double_t fBeta;
......@@ -207,7 +217,7 @@ protected:
// Potential Hall C parameters. Mostly here for demonstration
Int_t fNPlanes; // Number of planes
UInt_t fMaxScinPerPlane,fMaxHodoScin; // max number of scin/plane; product of the first two
UInt_t fMaxScinPerPlane,fMaxHodoScin,fTotHodScin; // max number of scin/plane; product of the first two
Double_t fStartTimeCenter, fStartTimeSlop, fScinTdcToTime;
Double_t fTofTolerance;
Int_t fCosmicFlag; //
......@@ -216,6 +226,7 @@ protected:
Double_t fScinTdcMin, fScinTdcMax; // min and max TDC values
char** fPlaneNames;
UInt_t* fNPaddle; // Number of paddles per plane
Int_t fTrackBetaIncludeSinglePmtHits;
Double_t *fHodoNegAdcTimeWindowMin;
Double_t *fHodoNegAdcTimeWindowMax;
......@@ -265,6 +276,7 @@ protected:
Int_t fCheckEvent;
Int_t fEventType;
Int_t fEventNum;
Int_t fGoodTrack;
Double_t fScin2XZpos;
......@@ -407,6 +419,10 @@ protected:
std::vector<Int_t > fNClust; // # scins clusters for the plane
std::vector<std::vector<Int_t> > fClustSize; // # scin cluster size
std::vector<std::vector<Double_t> > fClustPos; // # scin cluster position
std::vector<Int_t > fNCluster; // # scins clusters for the plane
std::vector<std::vector<Int_t> > fClusterSize; // # scin cluster size
std::vector<std::vector<Double_t> > fClusterXPos; // # scin cluster position
std::vector<std::vector<Double_t> > fClusterYPos; // # scin cluster position
std::vector<Int_t > fThreeScin; // # scins three clusters for the plane
std::vector<Int_t > fGoodScinHitsX; // # hits in fid x range
// Could combine the above into a structure
......
src/THcParmList.cxx
View file @ e974befb
......@@ -43,50 +43,48 @@ An instance of THaTextvars is created to hold the string parameters.
#include <iomanip>
#include "TBufferJSON.h"
#include "nlohmann/json.hpp"
#include "TObjArray.h"
#include "TObjString.h"
#include "TSystem.h"
#include "nlohmann/json.hpp"
#include "THcParmList.h"
#include "THaVar.h"
#include "THaFormula.h"
#include "THaVar.h"
#include "THcParmList.h"
#include "TMath.h"
/* #incluce <algorithm> include <fstream> include <cstring> */
#include <iostream>
#include <fstream>
#include <cassert>
#include <cstdlib>
#include <stdexcept>
#include <fstream>
#include <iostream>
#include <memory>
#include <stdexcept>
using namespace std;
Int_t fDebug = 1; // Keep this at one while we're working on the code
Int_t fDebug = 1; // Keep this at one while we're working on the code
#if __cplusplus < 201103L
# define SMART_PTR auto_ptr
#define SMART_PTR auto_ptr
#else
# define SMART_PTR unique_ptr
#define SMART_PTR unique_ptr
#endif
ClassImp(THcParmList)
/// Create empty numerical and string parameter lists
THcParmList::THcParmList() : hcana::ConfigLogging<THaVarList>()
{
/// Create empty numerical and string parameter lists
THcParmList::THcParmList()
: hcana::ConfigLogging<THaVarList>() {
TextList = new THaTextvars;
}
inline static bool IsComment( const string& s, string::size_type pos )
{
return ( pos != string::npos && pos < s.length() &&
(s[pos] == '#' || s[pos] == ';' || s.substr(pos,2) == "//") );
inline static bool IsComment(const string& s, string::size_type pos) {
return (pos != string::npos && pos < s.length() &&
(s[pos] == '#' || s[pos] == ';' || s.substr(pos, 2) == "//"));
}
void THcParmList::Load( const char* fname, Int_t RunNumber )
{
void THcParmList::Load(const char* fname, Int_t RunNumber) {
/**
Most parameter files used in the ENGINE should work.
......@@ -127,138 +125,136 @@ The ENGINE CTP support parameter "blocks" which were marked with
static const char* const whtspc = " \t";
ifstream ifiles[100]; // Should use stack instead
ifstream ifiles[100]; // Should use stack instead
Int_t nfiles=0;
Int_t nfiles = 0;
ifiles[nfiles].open(fname);
if(ifiles[nfiles].is_open()) {
//cout << "Opening parameter file: [" << nfiles << "] " << fname << endl;
if (ifiles[nfiles].is_open()) {
// cout << "Opening parameter file: [" << nfiles << "] " << fname << endl;
nfiles++;
}
if(!nfiles) {
static const char* const here = "THcParmList::LoadFromFile";
//Error (here, "error opening parameter file %s",fname);
return; // Need a success argument returned
if (!nfiles) {
static const char* const here = "THcParmList::LoadFromFile";
// Error (here, "error opening parameter file %s",fname);
return; // Need a success argument returned
}
string line;
char varname[100];
Int_t InRunRange;
Int_t currentindex = 0;
Int_t linecount=0; // Count of non comment/blank lines
char varname[100];
Int_t InRunRange;
Int_t currentindex = 0;
Int_t linecount = 0; // Count of non comment/blank lines
varname[0] = '\0';
if(RunNumber > 0) {
InRunRange = 0; // Wait until run number range matching RunNumber is found
//cout << "Reading Parameters for run " << RunNumber << endl;
if (RunNumber > 0) {
InRunRange = 0; // Wait until run number range matching RunNumber is found
// cout << "Reading Parameters for run " << RunNumber << endl;
} else {
InRunRange = 1; // Interpret all lines
InRunRange = 1; // Interpret all lines
}
while(nfiles) {
while (nfiles) {
string current_comment("");
// EJB_Note: existing_comment is never used.
// string existing_comment("");
string::size_type start, pos = 0;
if(!getline(ifiles[nfiles-1],line)) {
ifiles[nfiles-1].close();
if (!getline(ifiles[nfiles - 1], line)) {
ifiles[nfiles - 1].close();
nfiles--;
// cout << nfiles << ": " << "Closed" << endl;
continue;
}
// Look for include statement
if(line.compare(0,strlen(INCLUDESTR),INCLUDESTR)==0) {
line.erase(0,strlen(INCLUDESTR));
if (line.compare(0, strlen(INCLUDESTR), INCLUDESTR) == 0) {
line.erase(0, strlen(INCLUDESTR));
pos = line.find_first_not_of(whtspc);
// Strip leading white space
if(pos != string::npos && pos > 0 && pos < line.length()) {
line.erase(0,pos);
if (pos != string::npos && pos > 0 && pos < line.length()) {
line.erase(0, pos);
}
char quotechar=line[0];
if(quotechar == '"' || quotechar == '\'') {
line.erase(0,1);
line.erase(line.find_first_of(quotechar));
char quotechar = line[0];
if (quotechar == '"' || quotechar == '\'') {
line.erase(0, 1);
line.erase(line.find_first_of(quotechar));
} else {
line.erase(line.find_first_of(whtspc));
line.erase(line.find_first_of(whtspc));
}
// cout << line << endl;
ifiles[nfiles].open(line.c_str());
if(ifiles[nfiles].is_open()) {
if (ifiles[nfiles].is_open()) {
_logger->info("Opening parameter file: [{}] {} ", nfiles, line);
//cout << "Opening parameter file: [" << nfiles << "] " << line << endl;
nfiles++;
// cout << "Opening parameter file: [" << nfiles << "] " << line << endl;
nfiles++;
}
continue;
}
// Blank line or comment?
if( line.empty()
|| (start = line.find_first_not_of( whtspc )) == string::npos
|| IsComment(line, start) )
if (line.empty() || (start = line.find_first_not_of(whtspc)) == string::npos ||
IsComment(line, start))
continue;
// Get rid of trailing comments and leading and trailing whitespace
// Need to save the comment and put it in the thVar
while( (pos = line.find_first_of("#;/", pos+1)) != string::npos ) {
if( IsComment(line, pos) ) {
current_comment.assign(line,pos+1,line.length());
line.erase(pos); // Strip off comment
// Strip leading white space from comment
//cout << "CommentA: " << current_comment << endl;
pos = current_comment.find_first_not_of(whtspc);
if(pos!=string::npos && pos > 0 && pos < current_comment.length()) {
current_comment.erase(0,pos);
}
//cout << "CommentB: " << current_comment << endl;
break;
while ((pos = line.find_first_of("#;/", pos + 1)) != string::npos) {
if (IsComment(line, pos)) {
current_comment.assign(line, pos + 1, line.length());
line.erase(pos); // Strip off comment
// Strip leading white space from comment
// cout << "CommentA: " << current_comment << endl;
pos = current_comment.find_first_not_of(whtspc);
if (pos != string::npos && pos > 0 && pos < current_comment.length()) {
current_comment.erase(0, pos);
}
// cout << "CommentB: " << current_comment << endl;
break;
}
}
pos = line.find_last_not_of( whtspc );
assert( pos != string::npos );
if( pos != string::npos && ++pos < line.length() )
pos = line.find_last_not_of(whtspc);
assert(pos != string::npos);
if (pos != string::npos && ++pos < line.length())
line.erase(pos);
pos = line.find_first_not_of(whtspc);
// Strip leading white space
if(pos != string::npos && pos > 0 && pos < line.length()) {
line.erase(0,pos);
if (pos != string::npos && pos > 0 && pos < line.length()) {
line.erase(0, pos);
}
// Ignore begin and end statements
if(line.compare(0,5,"begin")==0 ||
line.compare(0,3,"end")==0) {
if (line.compare(0, 5, "begin") == 0 || line.compare(0, 3, "end") == 0) {
cout << "Skipping: " << line << endl;
continue;
}
// Get rid of all white space not in quotes
// Step through one char at a time
pos = 0;
int inquote=0;
char quotechar=' ';
pos = 0;
int inquote = 0;
char quotechar = ' ';
// cout << "Unstripped line: |" << line << "|" << endl;
while(pos<line.length()) {
if(inquote) {
if(line[pos++] == quotechar) { // Possibly end of quoted string
if(line[pos] == quotechar) { // Protected quote
pos++; // Skip the protected quote
} else { // End of quoted string
inquote = 0;
quotechar = ' ';
pos++;
}
}
while (pos < line.length()) {
if (inquote) {
if (line[pos++] == quotechar) { // Possibly end of quoted string
if (line[pos] == quotechar) { // Protected quote
pos++; // Skip the protected quote
} else { // End of quoted string
inquote = 0;
quotechar = ' ';
pos++;
}
}
} else {
if(line[pos] == ' ' || line[pos] == '\t') {
line.erase(pos,1);
} else if(line[pos] == '"' || line[pos] == '\'') {
quotechar = line[pos++];\
inquote = 1;
} else {
pos++;
}
if (line[pos] == ' ' || line[pos] == '\t') {
line.erase(pos, 1);
} else if (line[pos] == '"' || line[pos] == '\'') {
quotechar = line[pos++];
inquote = 1;
} else {
pos++;
}
}
}
// cout << "Stripped line: |" << line << "|" << endl;
......@@ -269,96 +265,97 @@ The ENGINE CTP support parameter "blocks" which were marked with
linecount++;
// If RunNumber>0 and first line we encounter is not a run range, need to
// print an error
if(RunNumber>0 && nfiles==1) {
if(line.find_first_not_of("0123456789-,")==string::npos) { // Interpret as runnum range
// Interpret line as a list of comma separated run numbers or ranges
TString runnums(line.c_str());
SMART_PTR<TObjArray> runnumarr( runnums.Tokenize(",") );
Int_t nranges=runnumarr->GetLast()+1;
InRunRange = 0;
Int_t ind;
for(Int_t i=0;i<nranges;i++) {
TString runstr = ((TObjString *)runnumarr->At(i))->GetString();
if(runstr.IsDec()) { // A single run number
if(RunNumber == runstr.Atoi()) {
InRunRange = 1;
break;
}
} else if ((ind=runstr.First('-'))>=0) { // A run range
TString start=runstr(0,ind);
TString end=runstr(ind+1,runstr.Length());
if(start.IsDec() && end.IsDec()) {
if((RunNumber >= start.Atoi()) && (RunNumber <= end.Atoi())) {
InRunRange = 1;
break;
}
}
}
}
continue; // Skip to next line
if (RunNumber > 0 && nfiles == 1) {
if (line.find_first_not_of("0123456789-,") == string::npos) { // Interpret as runnum range
// Interpret line as a list of comma separated run numbers or ranges
TString runnums(line.c_str());
SMART_PTR<TObjArray> runnumarr(runnums.Tokenize(","));
Int_t nranges = runnumarr->GetLast() + 1;
InRunRange = 0;
Int_t ind;
for (Int_t i = 0; i < nranges; i++) {
TString runstr = ((TObjString*)runnumarr->At(i))->GetString();
if (runstr.IsDec()) { // A single run number
if (RunNumber == runstr.Atoi()) {
InRunRange = 1;
break;
}
} else if ((ind = runstr.First('-')) >= 0) { // A run range
TString start = runstr(0, ind);
TString end = runstr(ind + 1, runstr.Length());
if (start.IsDec() && end.IsDec()) {
if ((RunNumber >= start.Atoi()) && (RunNumber <= end.Atoi())) {
InRunRange = 1;
break;
}
}
}
}
continue; // Skip to next line
} else {
if(linecount==1) {
_logger->warn("THcParmList::Load in database mode but first line is not\n"
" a run number or run number range. Parameter definitions\n"
" will be ignored until a run number or range is specified.\n");
}
if (linecount == 1) {
_logger->warn("THcParmList::Load in database mode but first line is not\n"
" a run number or run number range. Parameter definitions\n"
" will be ignored until a run number or range is specified.\n");
}
}
}
if(!InRunRange) continue;
if (!InRunRange)
continue;
// Interpret left of = as var name
Int_t valuestartpos=0; // Stays zero if no = found
Int_t ttype = 0; // Are any of the values floating point?
if((pos=line.find_first_of("="))!=string::npos) {
strcpy(varname, (line.substr(0,pos)).c_str());
valuestartpos = pos+1;
currentindex = 0;
Int_t valuestartpos = 0; // Stays zero if no = found
Int_t ttype = 0; // Are any of the values floating point?
if ((pos = line.find_first_of("=")) != string::npos) {
strcpy(varname, (line.substr(0, pos)).c_str());
valuestartpos = pos + 1;
currentindex = 0;
}
// If first char after = is a quote, then this is a string assignment
if(line[valuestartpos] == '"' || line[valuestartpos] == '\'') {
if (line[valuestartpos] == '"' || line[valuestartpos] == '\'') {
quotechar = line[valuestartpos++];
// Scan until end of line or terminating quote
// valuestartpos++;
pos = valuestartpos;
while(pos<line.length()) {
if(line[pos++] == quotechar) { // Possibly end of quoted string
if(line[pos] == quotechar) { // Protected quote
pos++;
} else {
pos--;
break;
}
}
while (pos < line.length()) {
if (line[pos++] == quotechar) { // Possibly end of quoted string
if (line[pos] == quotechar) { // Protected quote
pos++;
} else {
pos--;
break;
}
}
}
if(TextList) {
// Should check that a numerical assignment doesn't exist, but for
// now, the same variable name can be used for strings and numbers
string varnames(varname);
AddString(varnames, line.substr(valuestartpos,pos-valuestartpos));
if (TextList) {
// Should check that a numerical assignment doesn't exist, but for
// now, the same variable name can be used for strings and numbers
string varnames(varname);
AddString(varnames, line.substr(valuestartpos, pos - valuestartpos));
}
continue;
}
TString values((line.substr(valuestartpos)).c_str());
TObjArray *vararr = values.Tokenize(",");
Int_t nvals = vararr->GetLast()+1;
TString values((line.substr(valuestartpos)).c_str());
TObjArray* vararr = values.Tokenize(",");
Int_t nvals = vararr->GetLast() + 1;
Int_t* ip=0;
Double_t* fp=0;
Int_t* ip = 0;
Double_t* fp = 0;
// or expressions
for(Int_t i=0;(ttype==0&&i<nvals);i++) {
TString valstr = ((TObjString *)vararr->At(i))->GetString();
if(valstr.IsFloat()) { // Is a number
if(valstr.Contains(".") || valstr.Contains("e",TString::kIgnoreCase)) {
ttype = 1;
break;
}
for (Int_t i = 0; (ttype == 0 && i < nvals); i++) {
TString valstr = ((TObjString*)vararr->At(i))->GetString();
if (valstr.IsFloat()) { // Is a number
if (valstr.Contains(".") || valstr.Contains("e", TString::kIgnoreCase)) {
ttype = 1;
break;
}
} else {
ttype = 2; // Force float if expression or Var
break;
ttype = 2; // Force float if expression or Var
break;
}
}
......@@ -382,150 +379,145 @@ The ENGINE CTP support parameter "blocks" which were marked with
//
// There is some code duplication here. Refactor later
Int_t newlength = currentindex + nvals;
THaVar* existingvar=Find(varname);
if(existingvar) {
Int_t newlength = currentindex + nvals;
THaVar* existingvar = Find(varname);
if (existingvar) {
string existingcomment;
existingcomment.assign(existingvar->GetTitle());
if(!existingcomment.empty()) {
current_comment.assign(existingcomment);
if (!existingcomment.empty()) {
current_comment.assign(existingcomment);
}
Int_t existingtype=existingvar->GetType();
Int_t existinglength=existingvar->GetLen();
if(newlength > existinglength ||
(existingtype == kInt && ttype > 0)) { // Length or type change needed
if(newlength < existinglength) newlength = existinglength;
Int_t newtype=-1;
if(ttype>0 || existingtype == kDouble) {
newtype = kDouble;
fp = new Double_t[newlength];
if(existingtype == kDouble) {
Double_t* existingp= (Double_t*) existingvar->GetValuePointer();
for(Int_t i=0;i<existinglength;i++) {
fp[i] = existingp[i];
}
} else if(existingtype == kInt) {
Int_t* existingp= (Int_t*) existingvar->GetValuePointer();
for(Int_t i=0;i<existinglength;i++) {
fp[i] = existingp[i];
}
} else {
cout << "Whoops!" << endl;
}
} else if(existingtype == kInt) { // Stays int
newtype = kInt;
ip = new Int_t[newlength];
Int_t* existingp= (Int_t*) existingvar->GetValuePointer();
for(Int_t i=0;i<existinglength;i++) {
ip[i] = existingp[i];
}
} else {
cout << "Whoops!" << endl;
}
// Now copy new values in
for(Int_t i=0;i<nvals;i++) {
TString valstr = ((TObjString *)vararr->At(i))->GetString();
if(newtype == kInt) {
ip[currentindex+i] = valstr.Atoi();
} else {
if(valstr.IsFloat()) {
fp[currentindex+i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula
("temp",valstr.Data(), (Bool_t) 0, this, 0);
fp[currentindex+i] = formula->Eval();
delete formula;
}
}
}
currentindex += nvals;
// Remove old variable and recreate
if(existingtype == kDouble) {
delete [] (Double_t*) existingvar->GetValuePointer();
} else if (existingtype == kInt) {
delete [] (Int_t*) existingvar->GetValuePointer();
}
RemoveName(varname);
char *arrayname=new char [strlen(varname)+20];
sprintf(arrayname,"%s[%d]",varname,newlength);
if(newtype == kInt) {
Define(arrayname, current_comment.c_str(), *ip);
} else {
Define(arrayname, current_comment.c_str(), *fp);
}
delete[] arrayname;
Int_t existingtype = existingvar->GetType();
Int_t existinglength = existingvar->GetLen();
if (newlength > existinglength ||
(existingtype == kInt && ttype > 0)) { // Length or type change needed
if (newlength < existinglength)
newlength = existinglength;
Int_t newtype = -1;
if (ttype > 0 || existingtype == kDouble) {
newtype = kDouble;
fp = new Double_t[newlength];
if (existingtype == kDouble) {
Double_t* existingp = (Double_t*)existingvar->GetValuePointer();
for (Int_t i = 0; i < existinglength; i++) {
fp[i] = existingp[i];
}
} else if (existingtype == kInt) {
Int_t* existingp = (Int_t*)existingvar->GetValuePointer();
for (Int_t i = 0; i < existinglength; i++) {
fp[i] = existingp[i];
}
} else {
cout << "Whoops!" << endl;
}
} else if (existingtype == kInt) { // Stays int
newtype = kInt;
ip = new Int_t[newlength];
Int_t* existingp = (Int_t*)existingvar->GetValuePointer();
for (Int_t i = 0; i < existinglength; i++) {
ip[i] = existingp[i];
}
} else {
cout << "Whoops!" << endl;
}
// Now copy new values in
for (Int_t i = 0; i < nvals; i++) {
TString valstr = ((TObjString*)vararr->At(i))->GetString();
if (newtype == kInt) {
ip[currentindex + i] = valstr.Atoi();
} else {
if (valstr.IsFloat()) {
fp[currentindex + i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula("temp", valstr.Data(), (Bool_t)0, this, 0);
fp[currentindex + i] = formula->Eval();
delete formula;
}
}
}
currentindex += nvals;
// Remove old variable and recreate
if (existingtype == kDouble) {
delete[](Double_t*) existingvar->GetValuePointer();
} else if (existingtype == kInt) {
delete[](Int_t*) existingvar->GetValuePointer();
}
RemoveName(varname);
char* arrayname = new char[strlen(varname) + 20];
sprintf(arrayname, "%s[%d]", varname, newlength);
if (newtype == kInt) {
Define(arrayname, current_comment.c_str(), *ip);
} else {
Define(arrayname, current_comment.c_str(), *fp);
}
delete[] arrayname;
} else {
// Existing array long enough and of right type, just copy to it.
if(ttype == 0 && existingtype == kInt) {
Int_t* existingp= (Int_t*) existingvar->GetValuePointer();
for(Int_t i=0;i<nvals;i++) {
TString valstr = ((TObjString *)vararr->At(i))->GetString();
existingp[currentindex+i] = valstr.Atoi();
}
} else {
Double_t* existingp= (Double_t*) existingvar->GetValuePointer();
for(Int_t i=0;i<nvals;i++) {
TString valstr = ((TObjString *)vararr->At(i))->GetString();
if(valstr.IsFloat()) {
existingp[currentindex+i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula
("temp",valstr.Data(), (Bool_t) 0, this, 0);
existingp[currentindex+i] = formula->Eval();
delete formula;
}
}
}
currentindex += nvals;
// Existing array long enough and of right type, just copy to it.
if (ttype == 0 && existingtype == kInt) {
Int_t* existingp = (Int_t*)existingvar->GetValuePointer();
for (Int_t i = 0; i < nvals; i++) {
TString valstr = ((TObjString*)vararr->At(i))->GetString();
existingp[currentindex + i] = valstr.Atoi();
}
} else {
Double_t* existingp = (Double_t*)existingvar->GetValuePointer();
for (Int_t i = 0; i < nvals; i++) {
TString valstr = ((TObjString*)vararr->At(i))->GetString();
if (valstr.IsFloat()) {
existingp[currentindex + i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula("temp", valstr.Data(), (Bool_t)0, this, 0);
existingp[currentindex + i] = formula->Eval();
delete formula;
}
}
}
currentindex += nvals;
}
} else {
if(currentindex !=0) {
cout << "currentindex=" << currentindex << " shouldn't be!" << endl;
if (currentindex != 0) {
cout << "currentindex=" << currentindex << " shouldn't be!" << endl;
}
if(ttype==0) {
ip = new Int_t[nvals];
if (ttype == 0) {
ip = new Int_t[nvals];
} else {
fp = new Double_t[nvals];
fp = new Double_t[nvals];
}
for(Int_t i=0;i<nvals;i++) {
TString valstr = ((TObjString *)vararr->At(i))->GetString();
if(ttype==0) {
ip[i] = valstr.Atoi();
} else {
if(valstr.IsFloat()) {
fp[i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula
("temp",valstr.Data(), (Bool_t) 0, this, 0);
fp[i] = formula->Eval();
delete formula;
}
}
for (Int_t i = 0; i < nvals; i++) {
TString valstr = ((TObjString*)vararr->At(i))->GetString();
if (ttype == 0) {
ip[i] = valstr.Atoi();
} else {
if (valstr.IsFloat()) {
fp[i] = valstr.Atof();
} else {
THaFormula* formula = new THaFormula("temp", valstr.Data(), (Bool_t)0, this, 0);
fp[i] = formula->Eval();
delete formula;
}
}
}
currentindex = nvals;
char *arrayname=new char [strlen(varname)+20];
sprintf(arrayname,"%s[%d]",varname,nvals);
if(ttype==0) {
Define(arrayname, current_comment.c_str(), *ip);
char* arrayname = new char[strlen(varname) + 20];
sprintf(arrayname, "%s[%d]", varname, nvals);
if (ttype == 0) {
Define(arrayname, current_comment.c_str(), *ip);
} else {
Define(arrayname, current_comment.c_str(), *fp);
Define(arrayname, current_comment.c_str(), *fp);
}
delete[] arrayname;
}
delete vararr; // Discard result of Tokenize
delete vararr; // Discard result of Tokenize
// cout << line << endl;
}
return;
}
//_____________________________________________________________________________
Int_t THcParmList::LoadParmValues(const DBRequest* list, const char* prefix)
{
Int_t THcParmList::LoadParmValues(const DBRequest* list, const char* prefix) {
/**
The following example loads several parameters held in the `gHcParms`
parameter cache into scalar variables or arrays.
......@@ -550,77 +542,80 @@ zero), then there will be no error if the parameter is missing.
*/
const DBRequest *ti = list;
Int_t cnt=0;
Int_t this_cnt=0;
const DBRequest* ti = list;
Int_t cnt = 0;
Int_t this_cnt = 0;
if( !prefix ) prefix = "";
if (!prefix)
prefix = "";
while ( ti && ti->name ) {
string keystr(prefix); keystr.append(ti->name);
while (ti && ti->name) {
string keystr(prefix);
keystr.append(ti->name);
const char* key = keystr.c_str();
// cout <<"Now at "<<ti->name<<endl;
this_cnt = 0;
if(this->Find(key)) {
if (this->Find(key)) {
VarType ty = this->Find(key)->GetType();
if (ti->nelem>1) {
// it is an array, use the appropriateinterface
switch (ti->type) {
case (kDouble) :
this_cnt = GetArray(key,static_cast<Double_t*>(ti->var),ti->nelem);
break;
case (kInt) :
this_cnt = GetArray(key,static_cast<Int_t*>(ti->var),ti->nelem);
break;
default:
Error("THcParmList","Invalid type to read %s",key);
break;
}
if (ti->nelem > 1) {
// it is an array, use the appropriateinterface
switch (ti->type) {
case (kDouble):
this_cnt = GetArray(key, static_cast<Double_t*>(ti->var), ti->nelem);
break;
case (kInt):
this_cnt = GetArray(key, static_cast<Int_t*>(ti->var), ti->nelem);
break;
default:
Error("THcParmList", "Invalid type to read %s", key);
break;
}
} else {
switch (ti->type) {
case (kDouble) :
if(ty == kInt) {
*static_cast<Double_t*>(ti->var)=*(Int_t *)this->Find(key)->GetValuePointer();
} else if (ty == kDouble) {
*static_cast<Double_t*>(ti->var)=*(Double_t *)this->Find(key)->GetValuePointer();
} else {
cout << "*** ERROR!!! Type Mismatch " << key << endl;
}
this_cnt=1;
break;
case (kInt) :
if(ty == kInt) {
*static_cast<Int_t*>(ti->var)=*(Int_t *)this->Find(key)->GetValuePointer();
} else if (ty == kDouble) {
*static_cast<Int_t*>(ti->var)=TMath::Nint(*(Double_t *)this->Find(key)->GetValuePointer());
_logger->warn("Rounded {} to nearest integer ",key);
} else {
_logger->error("Type Mismatch {} ",key);
}
this_cnt=1;
break;
default:
_logger->error("THcParmList: Invalid type to read {} ",key);
switch (ti->type) {
case (kDouble):
if (ty == kInt) {
*static_cast<Double_t*>(ti->var) = *(Int_t*)this->Find(key)->GetValuePointer();
} else if (ty == kDouble) {
*static_cast<Double_t*>(ti->var) = *(Double_t*)this->Find(key)->GetValuePointer();
} else {
cout << "*** ERROR!!! Type Mismatch " << key << endl;
}
this_cnt = 1;
break;
}
case (kInt):
if (ty == kInt) {
*static_cast<Int_t*>(ti->var) = *(Int_t*)this->Find(key)->GetValuePointer();
} else if (ty == kDouble) {
*static_cast<Int_t*>(ti->var) =
TMath::Nint(*(Double_t*)this->Find(key)->GetValuePointer());
_logger->warn("Rounded {} to nearest integer ", key);
} else {
_logger->error("Type Mismatch {} ", key);
}
this_cnt = 1;
break;
default:
_logger->error("THcParmList: Invalid type to read {} ", key);
break;
}
}
} else { // See if it is a text variable
} else { // See if it is a text variable
const char* value = GetString(key);
if(value) {
this_cnt = 1;
if(ti->type == kString) {
*((string*)ti->var) = string(value);
} else if (ti->type == kTString) {
*((TString*)ti->var) = (TString) value;
} else {
Error("THcParmList","No conversion for strings: %s",key);
}
if (value) {
this_cnt = 1;
if (ti->type == kString) {
*((string*)ti->var) = string(value);
} else if (ti->type == kTString) {
*((TString*)ti->var) = (TString)value;
} else {
Error("THcParmList", "No conversion for strings: %s", key);
}
}
}
if (this_cnt<=0) {
if ( !ti->optional ) {
if (this_cnt <= 0) {
if (!ti->optional) {
string msg = string("Could not find `") + key + "` in database!";
throw std::runtime_error("<THcParmList::LoadParmValues>: " + msg);
}
......@@ -633,67 +628,68 @@ zero), then there will be no error if the parameter is missing.
// READING AN ARRAY INTO A C-style ARRAY
//_____________________________________________________________________________
Int_t THcParmList::GetArray(const char* attr, Int_t* array, Int_t size)
{
return ReadArray(attr,array,size);
Int_t THcParmList::GetArray(const char* attr, Int_t* array, Int_t size) {
return ReadArray(attr, array, size);
}
//_____________________________________________________________________________
Int_t THcParmList::GetArray(const char* attr, Double_t* array, Int_t size)
{
return ReadArray(attr,array,size);
Int_t THcParmList::GetArray(const char* attr, Double_t* array, Int_t size) {
return ReadArray(attr, array, size);
}
//_____________________________________________________________________________
template<class T>
Int_t THcParmList::ReadArray(const char* attrC, T* array, Int_t size)
{
template <class T>
Int_t THcParmList::ReadArray(const char* attrC, T* array, Int_t size) {
/**
No resizing is done, so only 'size' elements may be stored.
*/
Int_t cnt=0;
Int_t cnt = 0;
THaVar *var = Find(attrC);
if(!var) return(cnt);
THaVar* var = Find(attrC);
if (!var)
return (cnt);
VarType ty = var->GetType();
if( ty != kInt && ty != kDouble) {
if (ty != kInt && ty != kDouble) {
cout << "*** ERROR: " << attrC << " has unsupported type " << ty << endl;
return(cnt);
return (cnt);
}
Int_t sz = var->GetLen();
const void *vp = var->GetValuePointer();
if(size != sz) {
Int_t sz = var->GetLen();
const void* vp = var->GetValuePointer();
if (size > sz) {
_logger->error("requested {} elements of {} which has only {} elements", size, attrC, sz);
exit(EXIT_FAILURE);
} else if (size < sz) {
_logger->warn("requested {} elements of {} which has length {}", size, attrC, sz);
}
if(size<sz) sz = size;
if (size < sz)
sz = size;
Int_t donint = 0;
if(ty == kDouble && typeid(array[0]) == typeid(Int_t)) {
donint = 1; // Use nint when putting doubles in nint
if (ty == kDouble && typeid(array[0]) == typeid(Int_t)) {
donint = 1; // Use nint when putting doubles in nint
cout << "*** WARNING!!! Rounded " << attrC << " elements to nearest integer " << endl;
}
for(cnt=0;cnt<sz;cnt++) {
if(ty == kInt) {
for (cnt = 0; cnt < sz; cnt++) {
if (ty == kInt) {
array[cnt] = ((Int_t*)vp)[cnt];
} else
if(donint) {
array[cnt] = TMath::Nint(((Double_t*)vp)[cnt]);
} else {
array[cnt] = ((Double_t*)vp)[cnt];
}
} else if (donint) {
array[cnt] = TMath::Nint(((Double_t*)vp)[cnt]);
} else {
array[cnt] = ((Double_t*)vp)[cnt];
}
}
return(cnt);
return (cnt);
}
//_____________________________________________________________________________
std::string THcParmList::PrintJSON(int run_number ) const {
TIter next(this);
std::string THcParmList::PrintJSON(int run_number) const {
TIter next(this);
nlohmann::json j;
while( THaVar* obj = (THaVar*) next() ) {
while (THaVar* obj = (THaVar*)next()) {
if( obj->IsBasic() ) {
if( obj->IsArray() ) {
if( obj->GetLen() == 1 ) {
if (obj->IsBasic()) {
if (obj->IsArray()) {
if (obj->GetLen() == 1) {
j[obj->GetName()] = obj->GetValue();
} else {
j[obj->GetName()] = obj->GetValues();
......@@ -701,133 +697,127 @@ std::string THcParmList::PrintJSON(int run_number ) const {
} else {
obj->Print();
}
} else {
} else {
obj->Print();
}
}
nlohmann::json jrun;
jrun[std::to_string(run_number)] = j;
//std::cout << j.dump(2) << "\n";
// std::cout << j.dump(2) << "\n";
// write prettified JSON to another file
//std::ofstream o("pretty.json");
//o << std::setw(4) << jrun << std::endl;
// std::ofstream o("pretty.json");
// o << std::setw(4) << jrun << std::endl;
return jrun.dump(2);
}
void THcParmList::PrintFull( Option_t* option ) const
{
void THcParmList::PrintFull(Option_t* option) const {
THaVarList::PrintFull(option);
TextList->Print();
}
#ifdef WITH_CCDB
//_____________________________________________________________________________
Int_t THcParmList::OpenCCDB(Int_t runnum)
{
Int_t THcParmList::OpenCCDB(Int_t runnum) {
// Use the Environment variable "CCDB_CONNECTION" as the
// connection string
const char* connection_string = gSystem->Getenv("CCDB_CONNECTION");
return(OpenCCDB(runnum,connection_string));
return (OpenCCDB(runnum, connection_string));
}
Int_t THcParmList::OpenCCDB(Int_t runnum, const char* connection_string)
{
Int_t THcParmList::OpenCCDB(Int_t runnum, const char* connection_string) {
// Connect to a CCDB database pointed to by connection_string
std::string s (connection_string);
CCDB_obj = new SQLiteCalibration(runnum);
std::string s(connection_string);
CCDB_obj = new SQLiteCalibration(runnum);
Int_t result = CCDB_obj->Connect(s);
if(!result) return -1; // Need some error codes
if (!result)
return -1; // Need some error codes
cout << "Opened " << s << " for run " << runnum << endl;
return 0;
}
Int_t THcParmList::CloseCCDB()
{
Int_t THcParmList::CloseCCDB() {
delete CCDB_obj;
return(0);
return (0);
}
Int_t THcParmList::LoadCCDBDirectory(const char* directory,
const char* prefix)
{
Int_t THcParmList::LoadCCDBDirectory(const char* directory, const char* prefix) {
// Load all parameters in directory
// Prepend prefix onto the name of each
std::string dirname (directory);
std::string dirname(directory);
if(dirname[dirname.length()-1]!='/') {
if (dirname[dirname.length() - 1] != '/') {
dirname.append("/");
}
Int_t dirlen=dirname.length();
Int_t dirlen = dirname.length();
vector<string> namepaths;
CCDB_obj->GetListOfNamepaths(namepaths);
for(UInt_t iname=0;iname<namepaths.size();iname++) {
std::string varname (namepaths[iname]);
if(varname.compare(0,dirlen,dirname) == 0) {
varname.replace(0,dirlen,prefix);
for (UInt_t iname = 0; iname < namepaths.size(); iname++) {
std::string varname(namepaths[iname]);
if (varname.compare(0, dirlen, dirname) == 0) {
varname.replace(0, dirlen, prefix);
// cout << namepaths[iname] << " -> " << varname << endl;
// To what extent is there duplication here with Load() method?
// Retrieve assignment
Assignment* assignment = CCDB_obj->GetAssignment(namepaths[iname], true);
ConstantsTypeColumn::ColumnTypes ccdbtype=assignment->GetValueType(0);
Int_t ccdbncolumns=assignment->GetColumnsCount();
Int_t ccdbnrows=assignment->GetRowsCount();
std::string title = assignment->GetTypeTable()->GetComment();
Assignment* assignment = CCDB_obj->GetAssignment(namepaths[iname], true);
ConstantsTypeColumn::ColumnTypes ccdbtype = assignment->GetValueType(0);
Int_t ccdbncolumns = assignment->GetColumnsCount();
Int_t ccdbnrows = assignment->GetRowsCount();
std::string title = assignment->GetTypeTable()->GetComment();
// Only load single column tables
if(ccdbncolumns == 1) {
THaVar* existingvar=Find(varname.c_str());
// Need to append [size] to end of varname
char sizestring[20];
sprintf(sizestring,"[%d]",ccdbnrows);
std::string size_str (sizestring);
std::string varnamearray (varname);
varnamearray.append(size_str);
// Select data type
if(ccdbtype==ConstantsTypeColumn::cIntColumn) {
vector<vector<int> > data;
CCDB_obj->GetCalib(data, namepaths[iname]);
if(existingvar) {
RemoveName(varname.c_str());
}
Int_t* ip = new Int_t[data.size()];
for(UInt_t row=0;row<data.size(); row++) {
ip[row] = data[row][0];
}
Define(varnamearray.c_str(), title.c_str(), *ip);
} else if (ccdbtype==ConstantsTypeColumn::cDoubleColumn) {
vector<vector<double> > data;
CCDB_obj->GetCalib(data, namepaths[iname]);
if(existingvar) {
RemoveName(varname.c_str());
}
Double_t* fp = new Double_t[data.size()];
for(UInt_t row=0;row<data.size(); row++) {
fp[row] = data[row][0];
}
Define(varnamearray.c_str(), title.c_str(), *fp);
} else if (ccdbtype==ConstantsTypeColumn::cStringColumn) {
if(ccdbnrows > 1) {
cout << namepaths[iname] << ": Only first element of CCDB string array loaded." << endl;
}
vector<vector<string> > data;
CCDB_obj->GetCalib(data, namepaths[iname]);
AddString(varname, data[0][0]);
} else {
cout << namepaths[iname] << ": Unsupported CCDB data type: " << ccdbtype << endl;
}
if (ccdbncolumns == 1) {
THaVar* existingvar = Find(varname.c_str());
// Need to append [size] to end of varname
char sizestring[20];
sprintf(sizestring, "[%d]", ccdbnrows);
std::string size_str(sizestring);
std::string varnamearray(varname);
varnamearray.append(size_str);
// Select data type
if (ccdbtype == ConstantsTypeColumn::cIntColumn) {
vector<vector<int>> data;
CCDB_obj->GetCalib(data, namepaths[iname]);
if (existingvar) {
RemoveName(varname.c_str());
}
Int_t* ip = new Int_t[data.size()];
for (UInt_t row = 0; row < data.size(); row++) {
ip[row] = data[row][0];
}
Define(varnamearray.c_str(), title.c_str(), *ip);
} else if (ccdbtype == ConstantsTypeColumn::cDoubleColumn) {
vector<vector<double>> data;
CCDB_obj->GetCalib(data, namepaths[iname]);
if (existingvar) {
RemoveName(varname.c_str());
}
Double_t* fp = new Double_t[data.size()];
for (UInt_t row = 0; row < data.size(); row++) {
fp[row] = data[row][0];
}
Define(varnamearray.c_str(), title.c_str(), *fp);
} else if (ccdbtype == ConstantsTypeColumn::cStringColumn) {
if (ccdbnrows > 1) {
cout << namepaths[iname] << ": Only first element of CCDB string array loaded." << endl;
}
vector<vector<string>> data;
CCDB_obj->GetCalib(data, namepaths[iname]);
AddString(varname, data[0][0]);
} else {
cout << namepaths[iname] << ": Unsupported CCDB data type: " << ccdbtype << endl;
}
} else {
cout << namepaths[iname] << ": Multicolumn CCDB variables not supported" << endl;
cout << namepaths[iname] << ": Multicolumn CCDB variables not supported" << endl;
}
}
}
......
src/THcRawAdcHit.cxx
View file @ e974befb
......@@ -258,6 +258,10 @@ void THcRawAdcHit::SetRefTime(Int_t refTime) {
fHasRefTime = kTRUE;
}
void THcRawAdcHit::SetRefDiffTime(Int_t refDiffTime) {
fRefDiffTime = refDiffTime;
}
void THcRawAdcHit::SetSample(Int_t data) {
if (fNSamples >= fMaxNSamples) {
// throw std::out_of_range("`THcRawAdcHit::SetSample`: too many samples!");
......@@ -481,6 +485,22 @@ Int_t THcRawAdcHit::GetRefTime() const {
}
}
Bool_t THcRawAdcHit::HasRefTime() const { return fHasRefTime; }
Int_t THcRawAdcHit::GetRefDiffTime() const {
if (fHasRefTime) {
return fRefDiffTime;
}
else {
TString msg = TString::Format(
"`THcRawAdcHit::GetRefTime`: Reference time not available!"
);
throw std::runtime_error(msg.Data());
}
}
Bool_t THcRawAdcHit::HasRefTime() const {
return fHasRefTime;
}
ClassImp(THcRawAdcHit)
src/THcRawAdcHit.h
View file @ e974befb
......@@ -16,12 +16,14 @@ class THcRawAdcHit : public podd2::HitLogging<TObject> {
void SetData(Int_t data);
void SetSample(Int_t data);
void SetRefTime(Int_t refTime);
void SetRefDiffTime(Int_t refDiffTime);
void SetDataTimePedestalPeak(
Int_t data, Int_t time, Int_t pedestal, Int_t peak
);
Int_t GetRawData(UInt_t iPulse=0) const;
Double_t GetF250_PeakPedestalRatio() {return fPeakPedestalRatio;};
Int_t GetF250_NPedestalSamples() {return fNPedestalSamples;};
Double_t GetAverage(UInt_t iSampleLow, UInt_t iSampleHigh) const;
Int_t GetIntegral(UInt_t iSampleLow, UInt_t iSampleHigh) const;
......@@ -41,6 +43,7 @@ class THcRawAdcHit : public podd2::HitLogging<TObject> {
Int_t GetPulseTimeRaw(UInt_t iPulse=0) const;
Int_t GetSampleRaw(UInt_t iSample=0) const;
Int_t GetRefTime() const;
Int_t GetRefDiffTime() const;
Double_t GetPed() const;
Double_t GetPulseInt(UInt_t iPulse=0) const;
......@@ -82,6 +85,7 @@ class THcRawAdcHit : public podd2::HitLogging<TObject> {
Int_t fPulseTime[fMaxNPulses];
Int_t fSample[fMaxNSamples]; // the big buffer
Int_t fRefTime;
Int_t fRefDiffTime;
Bool_t fHasMulti;
Bool_t fHasRefTime;
......
src/THcRawDCHit.cxx
View file @ e974befb
......@@ -59,6 +59,17 @@ void THcRawDCHit::SetReference(Int_t signal, Int_t reference) {
}
}
void THcRawDCHit::SetReferenceDiff(Int_t signal, Int_t reference) {
if (signal == 0) {
fTdcHit.SetRefDiffTime(reference);
}
else {
throw std::out_of_range(
"`THcRawDCHit::SetReference`: only signal `0` available!"
);
}
}
Int_t THcRawDCHit::GetData(Int_t signal) {
if (signal == 0) {
......@@ -95,6 +106,17 @@ Int_t THcRawDCHit::GetReference(Int_t signal) {
}
}
Int_t THcRawDCHit::GetReferenceDiff(Int_t signal) {
if (signal == 0) {
return fTdcHit.GetRefDiffTime();
}
else {
throw std::out_of_range(
"`THcRawDCHit::GetReference`: only signal `0` available!"
);
}
}
THcRawHit::ESignalType THcRawDCHit::GetSignalType(Int_t signal) {
if (signal == 0) {
......
src/THcRawDCHit.h
View file @ e974befb
......@@ -18,10 +18,12 @@ class THcRawDCHit : public THcRawHit {
virtual void SetData(Int_t signal, Int_t data);
virtual void SetReference(Int_t signal, Int_t reference);
virtual void SetReferenceDiff(Int_t signal, Int_t reference);
virtual Int_t GetData(Int_t signal);
virtual Int_t GetRawData(Int_t signal);
virtual Int_t GetReference(Int_t signal);
virtual Int_t GetReferenceDiff(Int_t signal);
virtual ESignalType GetSignalType(Int_t signal);
virtual Int_t GetNSignals();
......
src/THcRawHit.h
View file @ e974befb
......@@ -29,20 +29,21 @@ public:
// virtual Bool_t operator==( const THcRawHit& ) = 0;
// virtual Bool_t operator!=( const THcRawHit& ) = 0;
virtual void SetData(Int_t signal, Int_t data) {};
virtual void SetSample(Int_t signal, Int_t data) {};
virtual void SetDataTimePedestalPeak(Int_t signal, Int_t data,
Int_t time, Int_t pedestal, Int_t peak) {};
virtual Int_t GetData(Int_t signal) {return 0;}; /* Ref time subtracted */
virtual Int_t GetRawData(Int_t signal) {return 0;} /* Ref time not subtracted */
virtual ESignalType GetSignalType(Int_t signal) {return kUndefined;}
virtual void SetData(Int_t /* signal */, Int_t /* data */) {};
virtual void SetSample(Int_t /* signal */, Int_t /* data */) {};
virtual void SetDataTimePedestalPeak(Int_t /* signal */, Int_t /* data */,
Int_t /* time */, Int_t /* pedestal */, Int_t /* peak */) {};
virtual Int_t GetData(Int_t /* signal */) {return 0;}; /* Ref time subtracted */
virtual Int_t GetRawData(Int_t /* signal */) {return 0;} /* Ref time not subtracted */
virtual ESignalType GetSignalType(Int_t /* signal */) {return kUndefined;}
virtual Int_t GetNSignals() { return 1;}
virtual void SetReference(Int_t signal, Int_t reference) {};
virtual void SetReferenceDiff(Int_t signal, Int_t reference) {};
virtual Bool_t HasReference(Int_t signal) {return kFALSE;};
virtual Int_t GetReference(Int_t signal) {return 0;};
virtual void SetF250Params(Int_t NSA, Int_t NSB, Int_t NPED) {};
virtual void SetF250Params(Int_t /* NSA */, Int_t /* NSB */, Int_t /* NPED */) {};
// Derived objects must be sortable and supply Compare method
// virtual Bool_t IsSortable () const {return kFALSE; }
......
src/THcRawHodoHit.cxx
View file @ e974befb
......@@ -105,6 +105,18 @@ void THcRawHodoHit::SetReference(Int_t signal, Int_t reference) {
}
}
void THcRawHodoHit::SetReferenceDiff(Int_t signal, Int_t referenceDiff) {
if (signal < fNAdcSignals) {
fAdcHits[signal].SetRefDiffTime(referenceDiff);
} else if (signal < fNAdcSignals+fNTdcSignals) {
fTdcHits[signal-fNAdcSignals].SetRefDiffTime(referenceDiff);
} else {
throw std::out_of_range(
"`THcRawHodoHit::SetReference`: only signals `2` and `3` available!"
);
}
}
Int_t THcRawHodoHit::GetData(Int_t signal) {
if (0 <= signal && signal < fNAdcSignals) {
......@@ -147,6 +159,17 @@ Int_t THcRawHodoHit::GetReference(Int_t signal) {
}
}
Int_t THcRawHodoHit::GetReferenceDiff(Int_t signal) {
if (fNAdcSignals <= signal && signal < fNAdcSignals+fNTdcSignals) {
return fTdcHits[signal-fNAdcSignals].GetRefDiffTime();
}
else {
throw std::out_of_range(
"`THcRawHodoHit::GetReference`: only signals `2` and `3` available!"
);
}
}
THcRawHit::ESignalType THcRawHodoHit::GetSignalType(Int_t signal) {
if (0 <= signal && signal < fNAdcSignals) {
......
src/THcRawHodoHit.h
View file @ e974befb
......@@ -25,10 +25,12 @@ class THcRawHodoHit : public THcRawHit {
Int_t signal, Int_t data, Int_t time, Int_t pedestal, Int_t peak
);
virtual void SetReference(Int_t signal, Int_t reference);
virtual void SetReferenceDiff(Int_t signal, Int_t referenceDiff);
virtual Int_t GetData(Int_t signal);
virtual Int_t GetRawData(Int_t signal);
virtual Int_t GetReference(Int_t signal);
virtual Int_t GetReferenceDiff(Int_t signal);
virtual ESignalType GetSignalType(Int_t signal);
virtual Int_t GetNSignals();
......
src/THcRawShowerHit.cxx
View file @ e974befb
......@@ -92,6 +92,16 @@ void THcRawShowerHit::SetReference(Int_t signal, Int_t reference) {
}
}
void THcRawShowerHit::SetReferenceDiff(Int_t signal, Int_t reference) {
if (signal < fNAdcSignals) {
fAdcHits[signal].SetRefDiffTime(reference);
} else {
throw std::out_of_range(
"`THcRawHodoHit::SetReference`: only signals `2` and `3` available!"
);
}
}
Int_t THcRawShowerHit::GetData(Int_t signal) {
if (0 <= signal && signal < fNAdcSignals) {
......
src/THcRawShowerHit.h
View file @ e974befb
......@@ -22,6 +22,7 @@ class THcRawShowerHit : public THcRawHit {
Int_t signal, Int_t data, Int_t time, Int_t pedestal, Int_t peak
);
virtual void SetReference(Int_t signal, Int_t reference);
virtual void SetReferenceDiff(Int_t signal, Int_t reference);
virtual Int_t GetData(Int_t signal);
virtual Int_t GetRawData(Int_t signal);
......
src/THcRawTdcHit.cxx
View file @ e974befb
......@@ -95,6 +95,7 @@ THcRawTdcHit& THcRawTdcHit::operator=(const THcRawTdcHit& right) {
fTime[iHit] = right.fTime[iHit];
}
fRefTime = right.fRefTime;
fRefDiffTime = right.fRefDiffTime;
fHasRefTime = right.fHasRefTime;
fNHits = right.fNHits;
}
......@@ -113,6 +114,7 @@ void THcRawTdcHit::Clear(Option_t* opt) {
fTime[iHit] = 0;
}
fRefTime = 0;
fRefDiffTime = 0;
fHasRefTime = kFALSE;
fNHits = 0;
}
......@@ -138,6 +140,10 @@ void THcRawTdcHit::SetRefTime(Int_t refTime) {
fHasRefTime = kTRUE;
}
void THcRawTdcHit::SetRefDiffTime(Int_t refDiffTime) {
fRefDiffTime = refDiffTime;
}
Int_t THcRawTdcHit::GetTimeRaw(UInt_t iHit) const {
if (iHit < fNHits) {
......@@ -177,6 +183,18 @@ Int_t THcRawTdcHit::GetRefTime() const {
}
}
Int_t THcRawTdcHit::GetRefDiffTime() const {
if (fHasRefTime) {
return fRefDiffTime;
}
else {
TString msg = TString::Format(
"`THcRawTdcHit::GetRefDiffTime`: Reference time not available!"
);
throw std::runtime_error(msg.Data());
}
}
Bool_t THcRawTdcHit::HasRefTime() const {
return fHasRefTime;
......
src/THcRawTdcHit.h
View file @ e974befb
......@@ -14,10 +14,12 @@ class THcRawTdcHit : public TObject {
void SetTime(Int_t time);
void SetRefTime(Int_t refTime);
void SetRefDiffTime(Int_t refDiffTime);
Int_t GetTimeRaw(UInt_t iHit=0) const;
Int_t GetTime(UInt_t iHit=0) const;
Int_t GetRefTime() const;
Int_t GetRefDiffTime() const;
Bool_t HasRefTime() const;
......@@ -30,6 +32,7 @@ class THcRawTdcHit : public TObject {
Int_t fTime[fMaxNHits];
Int_t fRefTime;
Int_t fRefDiffTime;
Bool_t fHasRefTime;
UInt_t fNHits;
......
src/THcScalerEvtHandler.cxx
View file @ e974befb
......@@ -32,55 +32,71 @@ To enable debugging you may try this in the setup script
~~~
\author E. Brash based on THaScalerEvtHandler by R. Michaels
*/
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <iterator>
#include <map>
#include <sstream>
#include "THaEvtTypeHandler.h"
#include "THcScalerEvtHandler.h"
#include "GenScaler.h"
#include "Helper.h"
#include "Scaler1151.h"
#include "Scaler3800.h"
#include "Scaler3801.h"
#include "Scaler1151.h"
#include "Scaler560.h"
#include "Scaler9001.h"
#include "Scaler9250.h"
#include "THaCodaData.h"
#include "THaEvData.h"
#include "THcParmList.h"
#include "THcGlobals.h"
#include "THaEvtTypeHandler.h"
#include "THaGlobals.h"
#include "TNamed.h"
#include "THaVarList.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THcScalerEvtHandler.h"
#include "TMath.h"
#include "TString.h"
#include "TNamed.h"
#include "TROOT.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <map>
#include <iterator>
#include "THaVarList.h"
#include "TString.h"
#include "VarDef.h"
#include "Helper.h"
using namespace std;
using namespace Decoder;
static const UInt_t ICOUNT = 1;
static const UInt_t IRATE = 2;
static const UInt_t ICURRENT = 3;
static const UInt_t ICURRENT = 3;
static const UInt_t ICHARGE = 4;
static const UInt_t ITIME = 5;
static const UInt_t ICUT = 6;
static const UInt_t ITIME = 5;
static const UInt_t ICUT = 6;
static const UInt_t MAXCHAN = 32;
static const UInt_t defaultDT = 4;
THcScalerEvtHandler::THcScalerEvtHandler(const char *name, const char* description)
: hcana::ConfigLogging<THaEvtTypeHandler>(name,description), evcount(0), evcountR(0.0), ifound(0), fNormIdx(-1),
fNormSlot(-1),
dvars(0),dvars_prev_read(0), dvarsFirst(0), fScalerTree(0), fUseFirstEvent(kTRUE),
fOnlySyncEvents(kFALSE), fOnlyBanks(kFALSE), fDelayedType(-1),
fClockChan(-1), fLastClock(0), fClockOverflows(0)
{
THcScalerEvtHandler::THcScalerEvtHandler(const char* name, const char* description)
: hcana::ConfigLogging<THaEvtTypeHandler>(name, description),
fBCM_Gain(0),
fBCM_Offset(0),
fBCM_SatOffset(0),
fBCM_SatQuadratic(0),
fBCM_delta_charge(0),
evcount(0),
evcountR(0.0),
ifound(0),
fNormIdx(-1),
fNormSlot(-1),
dvars(0),
dvars_prev_read(0),
dvarsFirst(0),
fScalerTree(0),
fUseFirstEvent(kTRUE),
fOnlySyncEvents(kFALSE),
fOnlyBanks(kFALSE),
fDelayedType(-1),
fClockChan(-1),
fLastClock(0),
fClockOverflows(0) {
fRocSet.clear();
fModuleSet.clear();
scal_prev_read.clear();
......@@ -88,91 +104,94 @@ THcScalerEvtHandler::THcScalerEvtHandler(const char *name, const char* descripti
scal_overflows.clear();
}
THcScalerEvtHandler::~THcScalerEvtHandler()
{
THcScalerEvtHandler::~THcScalerEvtHandler() {
// The tree object is owned by ROOT since it gets associated wth the output
// file, so DO NOT delete it here.
// file, so DO NOT delete it here.
if (!TROOT::Initialized()) {
delete fScalerTree;
}
Podd::DeleteContainer(scalers);
Podd::DeleteContainer(scalerloc);
delete [] dvars_prev_read;
delete [] dvars;
delete [] dvarsFirst;
delete [] fBCM_Gain;
delete [] fBCM_Offset;
delete [] fBCM_delta_charge;
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
delete[] dvars_prev_read;
delete[] dvars;
delete[] dvarsFirst;
delete[] fBCM_Gain;
delete[] fBCM_Offset;
delete[] fBCM_SatOffset;
delete[] fBCM_SatQuadratic;
delete[] fBCM_delta_charge;
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
}
Int_t THcScalerEvtHandler::End( THaRunBase* )
{
Int_t THcScalerEvtHandler::End(THaRunBase*) {
// Process any delayed events in order received
_param_logger->info("THcScalerEvtHandler::End Analyzing {} delayed scaler events", fDelayedEvents.size());
_param_logger->info("THcScalerEvtHandler::End Analyzing {} delayed scaler events",
fDelayedEvents.size());
for(std::vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it) {
for (std::vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end();
++it) {
UInt_t* rdata = *it;
AnalyzeBuffer(rdata,kFALSE);
AnalyzeBuffer(rdata, kFALSE);
}
if (fDebugFile) *fDebugFile << "scaler tree ptr "<<fScalerTree<<endl;
evNumberR = -1;
if (fScalerTree) fScalerTree->Fill();
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
if (fDebugFile)
*fDebugFile << "scaler tree ptr " << fScalerTree << endl;
evNumberR = -1;
if (fScalerTree)
fScalerTree->Fill();
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
if (fScalerTree) fScalerTree->Write();
if (fScalerTree)
fScalerTree->Write();
return 0;
}
Int_t THcScalerEvtHandler::ReadDatabase(const TDatime& date )
{
Int_t THcScalerEvtHandler::ReadDatabase(const TDatime& date) {
char prefix[2];
prefix[0]='g';
prefix[1]='\0';
fNumBCMs = 0;
DBRequest list[]={
{"NumBCMs",&fNumBCMs, kInt, 0, 1},
{0}
};
prefix[0] = 'g';
prefix[1] = '\0';
fNumBCMs = 0;
DBRequest list[] = {{"NumBCMs", &fNumBCMs, kInt, 0, 1}, {0}};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
//cout << " NUmber of BCMs = " << fNumBCMs << endl;
// cout << " NUmber of BCMs = " << fNumBCMs << endl;
//
if(fNumBCMs > 0) {
fBCM_Gain = new Double_t[fNumBCMs];
fBCM_Offset = new Double_t[fNumBCMs];
fBCM_delta_charge= new Double_t[fNumBCMs];
string bcm_namelist;
DBRequest list2[]={
{"BCM_Gain", fBCM_Gain, kDouble, (UInt_t) fNumBCMs},
{"BCM_Offset", fBCM_Offset, kDouble,(UInt_t) fNumBCMs},
{"BCM_Names", &bcm_namelist, kString},
{"BCM_Current_threshold", &fbcm_Current_Threshold, kDouble,0, 1},
{"BCM_Current_threshold_index", &fbcm_Current_Threshold_Index, kInt,0,1},
{0}
};
fbcm_Current_Threshold = 0.0;
if (fNumBCMs > 0) {
fBCM_Gain = new Double_t[fNumBCMs];
fBCM_Offset = new Double_t[fNumBCMs];
fBCM_SatOffset = new Double_t[fNumBCMs];
fBCM_SatQuadratic = new Double_t[fNumBCMs];
fBCM_delta_charge = new Double_t[fNumBCMs];
string bcm_namelist;
DBRequest list2[] = {{"BCM_Gain", fBCM_Gain, kDouble, (UInt_t)fNumBCMs},
{"BCM_Offset", fBCM_Offset, kDouble, (UInt_t)fNumBCMs},
{"BCM_SatQuadratic", fBCM_SatQuadratic, kDouble, (UInt_t)fNumBCMs, 1},
{"BCM_SatOffset", fBCM_SatOffset, kDouble, (UInt_t)fNumBCMs, 1},
{"BCM_Names", &bcm_namelist, kString},
{"BCM_Current_threshold", &fbcm_Current_Threshold, kDouble, 0, 1},
{"BCM_Current_threshold_index", &fbcm_Current_Threshold_Index, kInt, 0, 1},
{0}};
fbcm_Current_Threshold = 0.0;
fbcm_Current_Threshold_Index = 0;
for (Int_t i = 0; i < fNumBCMs; i++) {
fBCM_SatOffset[i] = 0.;
fBCM_SatQuadratic[i] = 0.;
}
gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
vector<string> bcm_names = vsplit(bcm_namelist);
for(Int_t i=0;i<fNumBCMs;i++) {
fBCM_Name.push_back(bcm_names[i]+".scal");
fBCM_delta_charge[i]=0.;
for (Int_t i = 0; i < fNumBCMs; i++) {
fBCM_Name.push_back(bcm_names[i] + ".scal");
fBCM_delta_charge[i] = 0.;
}
}
fTotalTime=0.;
fPrevTotalTime=0.;
fDeltaTime=-1.;
fTotalTime = 0.;
fPrevTotalTime = 0.;
fDeltaTime = -1.;
//
//
return kOK;
......@@ -188,23 +207,22 @@ void THcScalerEvtHandler::SetDelayedType(int evtype) {
*/
fDelayedType = evtype;
}
Int_t THcScalerEvtHandler::Analyze(THaEvData *evdata)
{
Int_t lfirst=1;
if ( !IsMyEvent(evdata->GetEvType()) ) return -1;
Int_t THcScalerEvtHandler::Analyze(THaEvData* evdata) {
Int_t lfirst = 1;
if (!IsMyEvent(evdata->GetEvType()))
return -1;
if (fDebugFile) {
*fDebugFile << endl << "---------------------------------- "<<endl<<endl;
*fDebugFile << "\nEnter THcScalerEvtHandler for fName = "<<fName<<endl;
*fDebugFile << endl << "---------------------------------- " << endl << endl;
*fDebugFile << "\nEnter THcScalerEvtHandler for fName = " << fName << endl;
EvDump(evdata);
}
/// \todo : Put this first event stuff in separte function
if (lfirst && !fScalerTree) {
lfirst = 0; // Can't do this in Init for some reason
TString sname1 = "TS";
......@@ -212,71 +230,73 @@ Int_t THcScalerEvtHandler::Analyze(THaEvData *evdata)
TString sname3 = fName + " Scaler Data";
if (fDebugFile) {
*fDebugFile << "\nAnalyze 1st time for fName = "<<fName<<endl;
*fDebugFile << sname2 << " " <<sname3<<endl;
*fDebugFile << "\nAnalyze 1st time for fName = " << fName << endl;
*fDebugFile << sname2 << " " << sname3 << endl;
}
fScalerTree = new TTree(sname2.Data(),sname3.Data());
fScalerTree = new TTree(sname2.Data(), sname3.Data());
fScalerTree->SetAutoSave(200000000);
TString name, tinfo;
name = "evcount";
name = "evcount";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &evcountR, tinfo.Data(), 4000);
name = "evNumber";
name = "evNumber";
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &evNumberR, tinfo.Data(), 4000);
for (size_t i = 0; i < scalerloc.size(); i++) {
name = scalerloc[i]->name;
name = scalerloc[i]->name;
tinfo = name + "/D";
fScalerTree->Branch(name.Data(), &dvars[i], tinfo.Data(), 4000);
}
} // if (lfirst && !fScalerTree)
} // if (lfirst && !fScalerTree)
UInt_t *rdata = (UInt_t*) evdata->GetRawDataBuffer();
UInt_t* rdata = (UInt_t*)evdata->GetRawDataBuffer();
if( evdata->GetEvType() == fDelayedType) { // Save this event for processing later
if (evdata->GetEvType() == fDelayedType) { // Save this event for processing later
Int_t evlen = evdata->GetEvLength();
UInt_t *datacopy = new UInt_t[evlen];
UInt_t* datacopy = new UInt_t[evlen];
fDelayedEvents.push_back(datacopy);
memcpy(datacopy,rdata,evlen*sizeof(UInt_t));
memcpy(datacopy, rdata, evlen * sizeof(UInt_t));
return 1;
} else { // A normal event
if (fDebugFile) *fDebugFile<<"\n\nTHcScalerEvtHandler :: Debugging event type "<<dec<<evdata->GetEvType()<< " event num = " << evdata->GetEvNum() << endl<<endl;
evNumber=evdata->GetEvNum();
} else { // A normal event
if (fDebugFile)
*fDebugFile << "\n\nTHcScalerEvtHandler :: Debugging event type " << dec
<< evdata->GetEvType() << " event num = " << evdata->GetEvNum() << endl
<< endl;
evNumber = evdata->GetEvNum();
evNumberR = evNumber;
Int_t ret;
if((ret=AnalyzeBuffer(rdata,fOnlySyncEvents))) {
if (fDebugFile) *fDebugFile << "scaler tree ptr "<<fScalerTree<<endl;
if (fScalerTree) fScalerTree->Fill();
if ((ret = AnalyzeBuffer(rdata, fOnlySyncEvents))) {
if (fDebugFile)
*fDebugFile << "scaler tree ptr " << fScalerTree << endl;
if (fScalerTree)
fScalerTree->Fill();
fScalerTree->AutoSave("SaveSelf");
}
return ret;
}
}
Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
{
Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync) {
// Parse the data, load local data arrays.
UInt_t *p = (UInt_t*) rdata;
UInt_t* p = (UInt_t*)rdata;
UInt_t *plast = p+*p; // Index to last word in the bank
UInt_t* plast = p + *p; // Index to last word in the bank
ifound=0;
while(p<plast) {
p++; // point to header
ifound = 0;
while (p < plast) {
p++; // point to header
if (fDebugFile) {
*fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p-1) << endl;
*fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p - 1) << endl;
}
if((*p & 0xff00) == 0x1000) { // Bank Containing banks
p++; // Now pointing to a bank in the bank
if ((*p & 0xff00) == 0x1000) { // Bank Containing banks
p++; // Now pointing to a bank in the bank
} else if (((*p & 0xff00) == 0x100) && (*p != 0xC0000100)) {
// Bank containing integers. Look for scalers
// This is either ROC bank containing integers or
......@@ -285,470 +305,539 @@ Int_t THcScalerEvtHandler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync)
// Assume that very first word is a scaler header
// At any point in the bank where the word is not a matching
// header, we stop.
UInt_t tag = (*p>>16) & 0xffff;
UInt_t num = (*p) & 0xff;
UInt_t *pnext = p+*(p-1); // Next bank
p++; // First data word
UInt_t tag = (*p >> 16) & 0xffff;
UInt_t num = (*p) & 0xff;
UInt_t* pnext = p + *(p - 1); // Next bank
p++; // First data word
// Skip over banks that can't contain scalers
// If SetOnlyBanks(kTRUE) called, fRocSet will be empty
// so only bank tags matching module types will be considered.
if(fModuleSet.find(tag)!=fModuleSet.end()) {
if(onlysync && num==0) {
ifound = 0;
return 0;
}
} else if (fRocSet.find(tag)==fRocSet.end()) {
p = pnext; // Fall through to end of the above else if
if (fModuleSet.find(tag) != fModuleSet.end()) {
if (onlysync && num == 0) {
ifound = 0;
return 0;
}
} else if (fRocSet.find(tag) == fRocSet.end()) {
p = pnext; // Fall through to end of the above else if
}
// Look for normalization scaler module first.
if(fNormIdx >= 0) {
UInt_t *psave = p;
while(p < pnext) {
if(scalers[fNormIdx]->IsSlot(*p)) {
scalers[fNormIdx]->Decode(p);
ifound = 1;
break;
}
p += scalers[fNormIdx]->GetNumChan() + 1;
}
p = psave;
if (fNormIdx >= 0) {
UInt_t* psave = p;
while (p < pnext) {
if (scalers[fNormIdx]->IsSlot(*p)) {
scalers[fNormIdx]->Decode(p);
ifound = 1;
break;
}
p += scalers[fNormIdx]->GetNumChan() + 1;
}
p = psave;
}
while(p < pnext) {
Int_t nskip = 0;
if(fDebugFile) {
*fDebugFile << "Scaler Header: " << hex << *p << dec;
}
for(size_t j=0; j<scalers.size(); j++) {
if(scalers[j]->IsSlot(*p)) {
nskip = scalers[j]->GetNumChan() + 1;
if((Int_t) j != fNormIdx) {
if(fDebugFile) {
*fDebugFile << " found (" << j << ") skip " << nskip << endl;
}
scalers[j]->Decode(p);
ifound = 1;
}
break;
}
}
if(nskip == 0) {
if(fDebugFile) {
*fDebugFile << endl;
}
break; // Didn't find a matching header
}
p = p + nskip;
while (p < pnext) {
Int_t nskip = 0;
if (fDebugFile) {
*fDebugFile << "Scaler Header: " << hex << *p << dec;
}
for (size_t j = 0; j < scalers.size(); j++) {
if (scalers[j]->IsSlot(*p)) {
nskip = scalers[j]->GetNumChan() + 1;
if ((Int_t)j != fNormIdx) {
if (fDebugFile) {
*fDebugFile << " found (" << j << ") skip " << nskip << endl;
}
scalers[j]->Decode(p);
ifound = 1;
}
break;
}
}
if (nskip == 0) {
if (fDebugFile) {
*fDebugFile << endl;
}
break; // Didn't find a matching header
}
p = p + nskip;
}
p = pnext;
} else {
p = p+*(p-1); // Skip to next bank
p = p + *(p - 1); // Skip to next bank
}
}
if (fDebugFile) {
*fDebugFile << "Finished with decoding. "<<endl;
*fDebugFile << " Found flag = "<<ifound<<endl;
*fDebugFile << "Finished with decoding. " << endl;
*fDebugFile << " Found flag = " << ifound << endl;
}
// HMS has headers which are different from SOS, but both are
// event type 0 and come here. If you found no headers, return.
if (!ifound) return 0;
if (!ifound)
return 0;
// The correspondance between dvars and the scaler and the channel
// will be driven by a scaler.map file -- later
Double_t scal_current=0;
UInt_t thisClock = scalers[fNormIdx]->GetData(fClockChan);
if(thisClock < fLastClock) { // Count clock scaler wrap arounds
Double_t scal_current = 0;
UInt_t thisClock = scalers[fNormIdx]->GetData(fClockChan);
if (thisClock < fLastClock) { // Count clock scaler wrap arounds
fClockOverflows++;
}
fTotalTime = (thisClock+(((Double_t) fClockOverflows)*kMaxUInt+fClockOverflows))/fClockFreq;
fTotalTime =
(thisClock + (((Double_t)fClockOverflows) * kMaxUInt + fClockOverflows)) / fClockFreq;
fLastClock = thisClock;
fDeltaTime= fTotalTime - fPrevTotalTime;
if (fDeltaTime==0) {
fDeltaTime = fTotalTime - fPrevTotalTime;
if (fDeltaTime == 0) {
cout << " ******************* Severe Warning ****************************" << endl;
cout << " In THcScalerEvtHandler have found fDeltaTime is zero !! " << endl;
cout << " ******************* Alert DAQ experts ****************************" << endl;
cout << " ******************* Alert DAQ experts ****************************" << endl;
}
fPrevTotalTime=fTotalTime;
Int_t nscal=0;
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t idx = scalerloc[i]->index;
fPrevTotalTime = fTotalTime;
Int_t nscal = 0;
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t idx = scalerloc[i]->index;
size_t ichan = scalerloc[i]->ichan;
if (evcount==0) {
if (fDebugFile) *fDebugFile << "Debug dvarsFirst "<<i<<" "<<ivar<<" "<<idx<<" "<<ichan<<endl;
if ((ivar < scalerloc.size()) &&
(idx < scalers.size()) &&
(ichan < MAXCHAN)){
if(fUseFirstEvent){
if (scalerloc[ivar]->ikind == ICOUNT){
UInt_t scaldata = scalers[idx]->GetData(ichan);
dvars[ivar] = scaldata;
scal_present_read.push_back(scaldata);
scal_prev_read.push_back(0);
scal_overflows.push_back(0);
dvarsFirst[ivar] = 0.0;
}
if (scalerloc[ivar]->ikind == ITIME){
dvars[ivar] =fTotalTime;
dvarsFirst[ivar] = 0;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvars[ivar] = (scalers[idx]->GetData(ichan))/fDeltaTime;
dvarsFirst[ivar] = dvars[ivar];
//printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan)); //checks
}
if(scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE){
Int_t bcm_ind=-1;
for(Int_t itemp =0; itemp<fNumBCMs;itemp++)
{
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match!=string::npos)
{
bcm_ind=itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvars[ivar]=0.;
if (bcm_ind != -1) dvars[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvars[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
fBCM_delta_charge[bcm_ind]=fDeltaTime*((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
dvars[ivar]+=fBCM_delta_charge[bcm_ind];
}
}
// printf("1st event %i index %i fBCMname %s scalerloc %s offset %f gain %f computed %f\n",evcount, bcm_ind, fBCM_Name[bcm_ind],scalerloc[ivar]->name.Data(),fBCM_Offset[bcm_ind],fBCM_Gain[bcm_ind],dvars[ivar]);
}
if (fDebugFile) *fDebugFile << " dvarsFirst "<<scalerloc[ivar]->ikind<<" "<<dvarsFirst[ivar]<<endl;
} else { //ifnotusefirstevent
if (scalerloc[ivar]->ikind == ICOUNT) {
dvarsFirst[ivar] = scalers[idx]->GetData(ichan) ;
scal_present_read.push_back(dvarsFirst[ivar]);
scal_prev_read.push_back(0);
}
if (scalerloc[ivar]->ikind == ITIME){
dvarsFirst[ivar] = fTotalTime;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvarsFirst[ivar] = (scalers[idx]->GetData(ichan))/fDeltaTime;
//printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan)); //checks
}
if(scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE)
{
Int_t bcm_ind=-1;
for(Int_t itemp =0; itemp<fNumBCMs;itemp++)
{
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match!=string::npos)
{
bcm_ind=itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvarsFirst[ivar]=0.0;
if (bcm_ind != -1) dvarsFirst[ivar]=((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvarsFirst[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
fBCM_delta_charge[bcm_ind]=fDeltaTime*((scalers[idx]->GetData(ichan))/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
dvarsFirst[ivar]+=fBCM_delta_charge[bcm_ind];
}
}
}
if (fDebugFile) *fDebugFile << " dvarsFirst "<<scalerloc[ivar]->ikind<<" "<<dvarsFirst[ivar]<<endl;
}
}
else {
cout << "THcScalerEvtHandler:: ERROR:: incorrect index "<<ivar<<" "<<idx<<" "<<ichan<<endl;
if (evcount == 0) {
if (fDebugFile)
*fDebugFile << "Debug dvarsFirst " << i << " " << ivar << " " << idx << " " << ichan
<< endl;
if ((ivar < scalerloc.size()) && (idx < scalers.size()) && (ichan < MAXCHAN)) {
if (fUseFirstEvent) {
if (scalerloc[ivar]->ikind == ICOUNT) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
dvars[ivar] = scaldata;
scal_present_read.push_back(scaldata);
scal_prev_read.push_back(0);
scal_overflows.push_back(0);
dvarsFirst[ivar] = 0.0;
}
if (scalerloc[ivar]->ikind == ITIME) {
dvars[ivar] = fTotalTime;
dvarsFirst[ivar] = 0;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvars[ivar] = (scalers[idx]->GetData(ichan)) / fDeltaTime;
dvarsFirst[ivar] = dvars[ivar];
// printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan)); //checks
}
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvars[ivar] = 0.;
if (bcm_ind != -1) {
dvars[ivar] = ((scalers[idx]->GetData(ichan)) / fDeltaTime - fBCM_Offset[bcm_ind]) /
fBCM_Gain[bcm_ind];
dvars[ivar] =
dvars[ivar] +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(dvars[ivar] - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
}
if (bcm_ind == fbcm_Current_Threshold_Index)
scal_current = dvars[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
Double_t cur_temp =
((scalers[idx]->GetData(ichan)) / fDeltaTime - fBCM_Offset[bcm_ind]) /
fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.0);
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
dvars[ivar] += fBCM_delta_charge[bcm_ind];
}
}
// printf("1st event %i index %i fBCMname %s scalerloc %s offset %f gain %f
// computed %f\n",evcount, bcm_ind,
// fBCM_Name[bcm_ind],scalerloc[ivar]->name.Data(),fBCM_Offset[bcm_ind],fBCM_Gain[bcm_ind],dvars[ivar]);
}
if (fDebugFile)
*fDebugFile << " dvarsFirst " << scalerloc[ivar]->ikind << " " << dvarsFirst[ivar]
<< endl;
} else { // ifnotusefirstevent
if (scalerloc[ivar]->ikind == ICOUNT) {
dvarsFirst[ivar] = scalers[idx]->GetData(ichan);
scal_present_read.push_back(dvarsFirst[ivar]);
scal_prev_read.push_back(0);
}
if (scalerloc[ivar]->ikind == ITIME) {
dvarsFirst[ivar] = fTotalTime;
}
if (scalerloc[ivar]->ikind == IRATE) {
dvarsFirst[ivar] = (scalers[idx]->GetData(ichan)) / fDeltaTime;
// printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan)); //checks
}
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvarsFirst[ivar] = 0.0;
if (bcm_ind != -1) {
dvarsFirst[ivar] =
((scalers[idx]->GetData(ichan)) / fDeltaTime - fBCM_Offset[bcm_ind]) /
fBCM_Gain[bcm_ind];
dvarsFirst[ivar] =
dvarsFirst[ivar] +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(dvars[ivar] - fBCM_SatOffset[bcm_ind], 0.0), 2.);
}
if (bcm_ind == fbcm_Current_Threshold_Index)
scal_current = dvarsFirst[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
Double_t cur_temp =
((scalers[idx]->GetData(ichan)) / fDeltaTime - fBCM_Offset[bcm_ind]) /
fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.);
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
dvarsFirst[ivar] += fBCM_delta_charge[bcm_ind];
}
}
}
if (fDebugFile)
*fDebugFile << " dvarsFirst " << scalerloc[ivar]->ikind << " " << dvarsFirst[ivar]
<< endl;
}
} else {
cout << "THcScalerEvtHandler:: ERROR:: incorrect index " << ivar << " " << idx << " "
<< ichan << endl;
}
}else{ // evcount != 0
if (fDebugFile) *fDebugFile << "Debug dvars "<<i<<" "<<ivar<<" "<<idx<<" "<<ichan<<endl;
if ((ivar < scalerloc.size()) &&
(idx < scalers.size()) &&
(ichan < MAXCHAN)) {
if (scalerloc[ivar]->ikind == ICOUNT) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
if(scaldata < scal_prev_read[nscal]) {
scal_overflows[nscal]++;
}
dvars[ivar] = scaldata + (1+((Double_t)kMaxUInt))*scal_overflows[nscal]
-dvarsFirst[ivar];
scal_present_read[nscal]=scaldata;
nscal++;
}
if (scalerloc[ivar]->ikind == ITIME) {
dvars[ivar] = fTotalTime;
}
if (scalerloc[ivar]->ikind == IRATE) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if(scaldata < scal_prev_read[nscal-1]) {
diff = (kMaxUInt-(scal_prev_read[nscal-1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal-1];
}
dvars[ivar] = diff/fDeltaTime;
// printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan));//checks
}
if(scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE)
{
Int_t bcm_ind=-1;
for(Int_t itemp =0; itemp<fNumBCMs;itemp++)
{
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match!=string::npos)
{
bcm_ind=itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvars[ivar]=0;
if (bcm_ind != -1) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if(scaldata < scal_prev_read[nscal-1]) {
diff = (kMaxUInt-(scal_prev_read[nscal-1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal-1];
}
dvars[ivar]=0.;
if (fDeltaTime>0) dvars[ivar]=(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
}
if (bcm_ind == fbcm_Current_Threshold_Index) scal_current= dvars[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if(scaldata < scal_prev_read[nscal-1]) {
diff = (kMaxUInt-(scal_prev_read[nscal-1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal-1];
}
fBCM_delta_charge[bcm_ind]=0;
if (fDeltaTime>0) fBCM_delta_charge[bcm_ind]=fDeltaTime*(diff/fDeltaTime-fBCM_Offset[bcm_ind])/fBCM_Gain[bcm_ind];
dvars[ivar]+=fBCM_delta_charge[bcm_ind];
}
}
// printf("event %i index %i fBCMname %s scalerloc %s offset %f gain %f computed %f\n",evcount, bcm_ind, fBCM_Name[bcm_ind],scalerloc[ivar]->name.Data(),fBCM_Offset[bcm_ind],fBCM_Gain[bcm_ind],dvars[ivar]);
}
if (fDebugFile) *fDebugFile << " dvars "<<scalerloc[ivar]->ikind<<" "<<dvars[ivar]<<endl;
} else { // evcount != 0
if (fDebugFile)
*fDebugFile << "Debug dvars " << i << " " << ivar << " " << idx << " " << ichan << endl;
if ((ivar < scalerloc.size()) && (idx < scalers.size()) && (ichan < MAXCHAN)) {
if (scalerloc[ivar]->ikind == ICOUNT) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
if (scaldata < scal_prev_read[nscal]) {
scal_overflows[nscal]++;
}
dvars[ivar] =
scaldata + (1 + ((Double_t)kMaxUInt)) * scal_overflows[nscal] - dvarsFirst[ivar];
scal_present_read[nscal] = scaldata;
nscal++;
}
if (scalerloc[ivar]->ikind == ITIME) {
dvars[ivar] = fTotalTime;
}
if (scalerloc[ivar]->ikind == IRATE) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if (scaldata < scal_prev_read[nscal - 1]) {
diff = (kMaxUInt - (scal_prev_read[nscal - 1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal - 1];
}
dvars[ivar] = diff / fDeltaTime;
// printf("%s %f\n",scalerloc[ivar]->name.Data(),scalers[idx]->GetRate(ichan));//checks
}
if (scalerloc[ivar]->ikind == ICURRENT || scalerloc[ivar]->ikind == ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scalerloc[ivar]->ikind == ICURRENT) {
dvars[ivar] = 0;
if (bcm_ind != -1) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if (scaldata < scal_prev_read[nscal - 1]) {
diff = (kMaxUInt - (scal_prev_read[nscal - 1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal - 1];
}
dvars[ivar] = 0.;
if (fDeltaTime > 0) {
Double_t cur_temp = (diff / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.);
dvars[ivar] = cur_temp;
}
}
if (bcm_ind == fbcm_Current_Threshold_Index)
scal_current = dvars[ivar];
}
if (scalerloc[ivar]->ikind == ICHARGE) {
if (bcm_ind != -1) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
UInt_t diff;
if (scaldata < scal_prev_read[nscal - 1]) {
diff = (kMaxUInt - (scal_prev_read[nscal - 1] - 1)) + scaldata;
} else {
diff = scaldata - scal_prev_read[nscal - 1];
}
fBCM_delta_charge[bcm_ind] = 0;
if (fDeltaTime > 0) {
Double_t cur_temp = (diff / fDeltaTime - fBCM_Offset[bcm_ind]) / fBCM_Gain[bcm_ind];
cur_temp =
cur_temp +
fBCM_SatQuadratic[bcm_ind] *
TMath::Power(TMath::Max(cur_temp - fBCM_SatOffset[bcm_ind], 0.0), 2.);
fBCM_delta_charge[bcm_ind] = fDeltaTime * cur_temp;
}
dvars[ivar] += fBCM_delta_charge[bcm_ind];
}
}
// printf("event %i index %i fBCMname %s scalerloc %s offset %f gain %f computed
//%f\n",evcount, bcm_ind,
// fBCM_Name[bcm_ind],scalerloc[ivar]->name.Data(),fBCM_Offset[bcm_ind],fBCM_Gain[bcm_ind],dvars[ivar]);
}
if (fDebugFile)
*fDebugFile << " dvars " << scalerloc[ivar]->ikind << " " << dvars[ivar] << endl;
} else {
cout << "THcScalerEvtHandler:: ERROR:: incorrect index "<<ivar<<" "<<idx<<" "<<ichan<<endl;
cout << "THcScalerEvtHandler:: ERROR:: incorrect index " << ivar << " " << idx << " "
<< ichan << endl;
}
}
}
//
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t idx = scalerloc[i]->index;
for (size_t i = 0; i < scalerloc.size(); i++) {
size_t ivar = scalerloc[i]->ivar;
size_t idx = scalerloc[i]->index;
size_t ichan = scalerloc[i]->ichan;
if (scalerloc[ivar]->ikind == ICUT+ICOUNT){
if (scalerloc[ivar]->ikind == ICUT + ICOUNT) {
UInt_t scaldata = scalers[idx]->GetData(ichan);
if ( scal_current > fbcm_Current_Threshold) {
UInt_t diff;
if(scaldata < dvars_prev_read[ivar]) {
diff = (kMaxUInt-(dvars_prev_read[ivar] - 1)) + scaldata;
} else {
diff = scaldata - dvars_prev_read[ivar];
}
dvars[ivar] += diff;
}
if (scal_current > fbcm_Current_Threshold) {
UInt_t diff;
if (scaldata < dvars_prev_read[ivar]) {
diff = (kMaxUInt - (dvars_prev_read[ivar] - 1)) + scaldata;
} else {
diff = scaldata - dvars_prev_read[ivar];
}
dvars[ivar] += diff;
}
dvars_prev_read[ivar] = scaldata;
}
if (scalerloc[ivar]->ikind == ICUT+ICHARGE){
Int_t bcm_ind=-1;
for(Int_t itemp =0; itemp<fNumBCMs;itemp++)
{
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match!=string::npos)
{
bcm_ind=itemp;
}
}
if ( scal_current > fbcm_Current_Threshold && bcm_ind != -1) {
dvars[ivar] += fBCM_delta_charge[bcm_ind];
}
if (scalerloc[ivar]->ikind == ICUT + ICHARGE) {
Int_t bcm_ind = -1;
for (Int_t itemp = 0; itemp < fNumBCMs; itemp++) {
size_t match = string(scalerloc[ivar]->name.Data()).find(string(fBCM_Name[itemp]));
if (match != string::npos) {
bcm_ind = itemp;
}
}
if (scal_current > fbcm_Current_Threshold && bcm_ind != -1) {
dvars[ivar] += fBCM_delta_charge[bcm_ind];
}
}
if (scalerloc[ivar]->ikind == ICUT+ITIME){
if ( scal_current > fbcm_Current_Threshold) {
dvars[ivar] += fDeltaTime;
}
if (scalerloc[ivar]->ikind == ICUT + ITIME) {
if (scal_current > fbcm_Current_Threshold) {
dvars[ivar] += fDeltaTime;
}
}
}
//
evcount = evcount + 1;
evcount = evcount + 1;
evcountR = evcount;
//
for (size_t j=0; j<scal_prev_read.size(); j++) scal_prev_read[j]=scal_present_read[j];
//
for (size_t j=0; j<scalers.size(); j++) scalers[j]->Clear("");
for (size_t j = 0; j < scal_prev_read.size(); j++)
scal_prev_read[j] = scal_present_read[j];
//
for (size_t j = 0; j < scalers.size(); j++)
scalers[j]->Clear("");
return 1;
}
THaAnalysisObject::EStatus THcScalerEvtHandler::Init(const TDatime& date)
{
THaAnalysisObject::EStatus THcScalerEvtHandler::Init(const TDatime& date) {
//
ReadDatabase(date);
const int LEN = 200;
char cbuf[LEN];
char cbuf[LEN];
fStatus = kOK;
fStatus = kOK;
fNormIdx = -1;
for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
it != fDelayedEvents.end(); ++it )
delete [] *it;
for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
delete[] * it;
fDelayedEvents.clear();
//cout << "Howdy ! We are initializing THcScalerEvtHandler !! name = "
// cout << "Howdy ! We are initializing THcScalerEvtHandler !! name = "
// << fName << endl;
if(eventtypes.size()==0) {
eventtypes.push_back(0); // Default Event Type
if (eventtypes.size() == 0) {
eventtypes.push_back(0); // Default Event Type
}
TString dfile;
dfile = fName + "scaler.txt";
// Parse the map file which defines what scalers exist and the global variables.
// Parse the map file which defines what scalers exist and the global variables.
TString sname0 = "Scalevt";
TString sname;
sname = fName+sname0;
sname = fName + sname0;
FILE *fi = OpenFile(sname.Data(), date);
if ( !fi ) {
cout << "Cannot find db file for "<<fName<<" scaler event handler"<<endl;
FILE* fi = OpenFile(sname.Data(), date);
if (!fi) {
cout << "Cannot find db file for " << fName << " scaler event handler" << endl;
return kFileError;
}
size_t minus1 = -1;
size_t pos1;
string scomment = "#";
string svariable = "variable";
string smap = "map";
size_t minus1 = -1;
size_t pos1;
string scomment = "#";
string svariable = "variable";
string smap = "map";
vector<string> dbline;
while( fgets(cbuf, LEN, fi) != NULL) {
while (fgets(cbuf, LEN, fi) != NULL) {
std::string sin(cbuf);
std::string sinput(sin.substr(0,sin.find_first_of("#")));
if (fDebugFile) *fDebugFile << "string input "<<sinput<<endl;
std::string sinput(sin.substr(0, sin.find_first_of("#")));
if (fDebugFile)
*fDebugFile << "string input " << sinput << endl;
dbline = vsplit(sinput);
if (dbline.size() > 0) {
pos1 = FindNoCase(dbline[0],scomment);
if (pos1 != minus1) continue;
pos1 = FindNoCase(dbline[0],svariable);
if (pos1 != minus1 && dbline.size()>4) {
string sdesc = "";
for (size_t j=5; j<dbline.size(); j++) sdesc = sdesc+" "+dbline[j];
UInt_t islot = atoi(dbline[1].c_str());
UInt_t ichan = atoi(dbline[2].c_str());
UInt_t ikind = atoi(dbline[3].c_str());
if (fDebugFile)
*fDebugFile << "add var "<<dbline[1]<<" desc = "<<sdesc<<" islot= "<<islot<<" "<<ichan<<" "<<ikind<<endl;
TString tsname(dbline[4].c_str());
TString tsdesc(sdesc.c_str());
AddVars(tsname,tsdesc,islot,ichan,ikind);
// add extra scaler which is cut on the current
if (ikind == ICOUNT ||ikind == ITIME ||ikind == ICHARGE ) {
tsname=tsname+"Cut";
AddVars(tsname,tsdesc,islot,ichan,ICUT+ikind);
}
pos1 = FindNoCase(dbline[0], scomment);
if (pos1 != minus1)
continue;
pos1 = FindNoCase(dbline[0], svariable);
if (pos1 != minus1 && dbline.size() > 4) {
string sdesc = "";
for (size_t j = 5; j < dbline.size(); j++)
sdesc = sdesc + " " + dbline[j];
UInt_t islot = atoi(dbline[1].c_str());
UInt_t ichan = atoi(dbline[2].c_str());
UInt_t ikind = atoi(dbline[3].c_str());
if (fDebugFile)
*fDebugFile << "add var " << dbline[1] << " desc = " << sdesc << " islot= " << islot
<< " " << ichan << " " << ikind << endl;
TString tsname(dbline[4].c_str());
TString tsdesc(sdesc.c_str());
AddVars(tsname, tsdesc, islot, ichan, ikind);
// add extra scaler which is cut on the current
if (ikind == ICOUNT || ikind == ITIME || ikind == ICHARGE) {
tsname = tsname + "Cut";
AddVars(tsname, tsdesc, islot, ichan, ICUT + ikind);
}
}
pos1 = FindNoCase(dbline[0],smap);
if (fDebugFile) *fDebugFile << "map ? "<<dbline[0]<<" "<<smap<<" "<<pos1<<" "<<dbline.size()<<endl;
if (pos1 != minus1 && dbline.size()>6) {
Int_t imodel, icrate, islot, inorm;
UInt_t header, mask;
char cdum[20];
sscanf(sinput.c_str(),"%s %d %d %d %x %x %d \n",cdum,&imodel,&icrate,&islot, &header, &mask, &inorm);
if ((fNormSlot >= 0) && (fNormSlot != inorm)) cout << "THcScalerEvtHandler::WARN: contradictory norm slot "<<fNormSlot<<" "<<inorm<<endl;
fNormSlot = inorm; // slot number used for normalization. This variable is not used but is checked.
Int_t clkchan = -1;
Double_t clkfreq = 1;
if (dbline.size()>8) {
clkchan = atoi(dbline[7].c_str());
clkfreq = 1.0*atoi(dbline[8].c_str());
fClockChan=clkchan;
fClockFreq=clkfreq;
}
if (fDebugFile) {
*fDebugFile << "map line "<<dec<<imodel<<" "<<icrate<<" "<<islot<<endl;
*fDebugFile <<" header 0x"<<hex<<header<<" 0x"<<mask<<dec<<" "<<inorm<<" "<<clkchan<<" "<<clkfreq<<endl;
}
switch (imodel) {
case 560:
scalers.push_back(new Scaler560(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 1151:
scalers.push_back(new Scaler1151(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 3800:
scalers.push_back(new Scaler3800(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 3801:
scalers.push_back(new Scaler3801(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 9001: // TI Scalers
scalers.push_back(new Scaler9001(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 9250: // FADC250 Scalers
scalers.push_back(new Scaler9250(icrate, islot));
if(!fOnlyBanks) fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
}
if (scalers.size() > 0) {
UInt_t idx = scalers.size()-1;
// Headers must be unique over whole event, not
// just within a ROC
scalers[idx]->SetHeader(header, mask);
// The normalization slot has the clock in it, so we automatically recognize it.
// fNormIdx is the index in scaler[] and
// fNormSlot is the slot#, checked for consistency
if (clkchan >= 0) {
scalers[idx]->SetClock(defaultDT, clkchan, clkfreq);
_param_logger->info("Setting scaler clock ... channel = {} ... freq = {} ",
clkchan, clkfreq);
if (fDebugFile) *fDebugFile <<"Setting scaler clock ... channel = "<<clkchan<<" ... freq = "<<clkfreq<<endl;
fNormIdx = idx;
if (islot != fNormSlot) cout << "THcScalerEvtHandler:: WARN: contradictory norm slot ! "<<islot<<endl;
}
}
pos1 = FindNoCase(dbline[0], smap);
if (fDebugFile)
*fDebugFile << "map ? " << dbline[0] << " " << smap << " " << pos1 << " "
<< dbline.size() << endl;
if (pos1 != minus1 && dbline.size() > 6) {
Int_t imodel, icrate, islot, inorm;
UInt_t header, mask;
char cdum[20];
sscanf(sinput.c_str(), "%s %d %d %d %x %x %d \n", cdum, &imodel, &icrate, &islot, &header,
&mask, &inorm);
if ((fNormSlot >= 0) && (fNormSlot != inorm))
cout << "THcScalerEvtHandler::WARN: contradictory norm slot " << fNormSlot << " "
<< inorm << endl;
fNormSlot =
inorm; // slot number used for normalization. This variable is not used but is checked.
Int_t clkchan = -1;
Double_t clkfreq = 1;
if (dbline.size() > 8) {
clkchan = atoi(dbline[7].c_str());
clkfreq = 1.0 * atoi(dbline[8].c_str());
fClockChan = clkchan;
fClockFreq = clkfreq;
}
if (fDebugFile) {
*fDebugFile << "map line " << dec << imodel << " " << icrate << " " << islot << endl;
*fDebugFile << " header 0x" << hex << header << " 0x" << mask << dec << " " << inorm
<< " " << clkchan << " " << clkfreq << endl;
}
switch (imodel) {
case 560:
scalers.push_back(new Scaler560(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 1151:
scalers.push_back(new Scaler1151(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 3800:
scalers.push_back(new Scaler3800(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 3801:
scalers.push_back(new Scaler3801(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 9001: // TI Scalers
scalers.push_back(new Scaler9001(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
case 9250: // FADC250 Scalers
scalers.push_back(new Scaler9250(icrate, islot));
if (!fOnlyBanks)
fRocSet.insert(icrate);
fModuleSet.insert(imodel);
break;
}
if (scalers.size() > 0) {
UInt_t idx = scalers.size() - 1;
// Headers must be unique over whole event, not
// just within a ROC
scalers[idx]->SetHeader(header, mask);
// The normalization slot has the clock in it, so we automatically recognize it.
// fNormIdx is the index in scaler[] and
// fNormSlot is the slot#, checked for consistency
if (clkchan >= 0) {
scalers[idx]->SetClock(defaultDT, clkchan, clkfreq);
_param_logger->info("Setting scaler clock ... channel = {} ... freq = {} ", clkchan,
clkfreq);
if (fDebugFile)
*fDebugFile << "Setting scaler clock ... channel = " << clkchan
<< " ... freq = " << clkfreq << endl;
fNormIdx = idx;
if (islot != fNormSlot)
cout << "THcScalerEvtHandler:: WARN: contradictory norm slot ! " << islot << endl;
}
}
}
}
}
// can't compare UInt_t to Int_t (compiler warning), so do this
nscalers=0;
for (size_t i=0; i<scalers.size(); i++) nscalers++;
nscalers = 0;
for (size_t i = 0; i < scalers.size(); i++)
nscalers++;
// need to do LoadNormScaler after scalers created and if fNormIdx found
if (fDebugFile) *fDebugFile <<"fNormIdx = "<<fNormIdx<<endl;
if (fDebugFile)
*fDebugFile << "fNormIdx = " << fNormIdx << endl;
if ((fNormIdx >= 0) && fNormIdx < nscalers) {
for (Int_t i = 0; i < nscalers; i++) {
if (i==fNormIdx) continue;
if (i == fNormIdx)
continue;
scalers[i]->LoadNormScaler(scalers[fNormIdx]);
}
}
......@@ -757,27 +846,26 @@ THaAnalysisObject::EStatus THcScalerEvtHandler::Init(const TDatime& date)
// This code is superseded by the parsing of a map file above. It's another way ...
if (fName == "Left") {
AddVars("TSbcmu1", "BCM x1 counts", 1, 4, ICOUNT);
AddVars("TSbcmu1r","BCM x1 rate", 1, 4, IRATE);
AddVars("TSbcmu1r", "BCM x1 rate", 1, 4, IRATE);
AddVars("TSbcmu3", "BCM u3 counts", 1, 5, ICOUNT);
AddVars("TSbcmu3r", "BCM u3 rate", 1, 5, IRATE);
AddVars("TSbcmu3r", "BCM u3 rate", 1, 5, IRATE);
} else {
AddVars("TSbcmu1", "BCM x1 counts", 0, 4, ICOUNT);
AddVars("TSbcmu1r","BCM x1 rate", 0, 4, IRATE);
AddVars("TSbcmu1r", "BCM x1 rate", 0, 4, IRATE);
AddVars("TSbcmu3", "BCM u3 counts", 0, 5, ICOUNT);
AddVars("TSbcmu3r", "BCM u3 rate", 0, 5, IRATE);
AddVars("TSbcmu3r", "BCM u3 rate", 0, 5, IRATE);
}
#endif
DefVars();
#ifdef HARDCODED
// This code is superseded by the parsing of a map file above. It's another way ...
if (fName == "Left") {
scalers.push_back(new Scaler1151(1,0));
scalers.push_back(new Scaler3800(1,1));
scalers.push_back(new Scaler3800(1,2));
scalers.push_back(new Scaler3800(1,3));
scalers.push_back(new Scaler1151(1, 0));
scalers.push_back(new Scaler3800(1, 1));
scalers.push_back(new Scaler3800(1, 2));
scalers.push_back(new Scaler3800(1, 3));
scalers[0]->SetHeader(0xabc00000, 0xffff0000);
scalers[1]->SetHeader(0xabc10000, 0xffff0000);
scalers[2]->SetHeader(0xabc20000, 0xffff0000);
......@@ -787,10 +875,10 @@ THaAnalysisObject::EStatus THcScalerEvtHandler::Init(const TDatime& date)
scalers[2]->LoadNormScaler(scalers[1]);
scalers[3]->LoadNormScaler(scalers[1]);
} else {
scalers.push_back(new Scaler3800(2,0));
scalers.push_back(new Scaler3800(2,0));
scalers.push_back(new Scaler1151(2,1));
scalers.push_back(new Scaler1151(2,2));
scalers.push_back(new Scaler3800(2, 0));
scalers.push_back(new Scaler3800(2, 0));
scalers.push_back(new Scaler1151(2, 1));
scalers.push_back(new Scaler1151(2, 2));
scalers[0]->SetHeader(0xceb00000, 0xffff0000);
scalers[1]->SetHeader(0xceb10000, 0xffff0000);
scalers[2]->SetHeader(0xceb20000, 0xffff0000);
......@@ -803,88 +891,87 @@ THaAnalysisObject::EStatus THcScalerEvtHandler::Init(const TDatime& date)
#endif
// Verify that the slots are not defined twice
for (UInt_t i1=0; i1 < scalers.size()-1; i1++) {
for (UInt_t i2=i1+1; i2 < scalers.size(); i2++) {
if (scalers[i1]->GetSlot()==scalers[i2]->GetSlot())
cout << "THcScalerEvtHandler:: WARN: same slot defined twice"<<endl;
for (UInt_t i1 = 0; i1 < scalers.size() - 1; i1++) {
for (UInt_t i2 = i1 + 1; i2 < scalers.size(); i2++) {
if (scalers[i1]->GetSlot() == scalers[i2]->GetSlot())
cout << "THcScalerEvtHandler:: WARN: same slot defined twice" << endl;
}
}
// Identify indices of scalers[] vector to variables.
for (UInt_t i=0; i < scalers.size(); i++) {
for (UInt_t i = 0; i < scalers.size(); i++) {
for (UInt_t j = 0; j < scalerloc.size(); j++) {
if (scalerloc[j]->islot==static_cast<UInt_t>(scalers[i]->GetSlot()))
scalerloc[j]->index = i;
if (scalerloc[j]->islot == static_cast<UInt_t>(scalers[i]->GetSlot()))
scalerloc[j]->index = i;
}
}
if(fDebugFile) *fDebugFile << "THcScalerEvtHandler:: Name of scaler bank "<<fName<<endl;
for (size_t i=0; i<scalers.size(); i++) {
if(fDebugFile) {
*fDebugFile << "Scaler # "<<i<<endl;
if (fDebugFile)
*fDebugFile << "THcScalerEvtHandler:: Name of scaler bank " << fName << endl;
for (size_t i = 0; i < scalers.size(); i++) {
if (fDebugFile) {
*fDebugFile << "Scaler # " << i << endl;
scalers[i]->SetDebugFile(fDebugFile);
scalers[i]->DebugPrint(fDebugFile);
}
}
//
return kOK;
}
void THcScalerEvtHandler::AddVars(TString name, TString desc, UInt_t islot,
UInt_t ichan, UInt_t ikind)
{
void THcScalerEvtHandler::AddVars(TString name, TString desc, UInt_t islot, UInt_t ichan,
UInt_t ikind) {
// need to add fName here to make it a unique variable. (Left vs Right HRS, for example)
TString name1 = fName + name;
TString desc1 = fName + desc;
// We don't yet know the correspondence between index of scalers[] and slots.
// Will put that in later.
scalerloc.push_back( new HCScalerLoc(name1, desc1, 0, islot, ichan, ikind,
scalerloc.size()) );
scalerloc.push_back(new HCScalerLoc(name1, desc1, 0, islot, ichan, ikind, scalerloc.size()));
}
void THcScalerEvtHandler::DefVars()
{
void THcScalerEvtHandler::DefVars() {
// called after AddVars has finished being called.
Nvars = scalerloc.size();
if (Nvars == 0) return;
dvars_prev_read = new UInt_t[Nvars]; // dvars_prev_read is a member of this class
dvars = new Double_t[Nvars]; // dvars is a member of this class
dvarsFirst = new Double_t[Nvars]; // dvarsFirst is a member of this class
memset(dvars, 0, Nvars*sizeof(Double_t));
memset(dvars_prev_read, 0, Nvars*sizeof(UInt_t));
memset(dvarsFirst, 0, Nvars*sizeof(Double_t));
if (Nvars == 0)
return;
dvars_prev_read = new UInt_t[Nvars]; // dvars_prev_read is a member of this class
dvars = new Double_t[Nvars]; // dvars is a member of this class
dvarsFirst = new Double_t[Nvars]; // dvarsFirst is a member of this class
memset(dvars, 0, Nvars * sizeof(Double_t));
memset(dvars_prev_read, 0, Nvars * sizeof(UInt_t));
memset(dvarsFirst, 0, Nvars * sizeof(Double_t));
if (gHaVars) {
if(fDebugFile) *fDebugFile << "THcScalerEVtHandler:: Have gHaVars "<<gHaVars<<endl;
if (fDebugFile)
*fDebugFile << "THcScalerEVtHandler:: Have gHaVars " << gHaVars << endl;
} else {
cout << "No gHaVars ?! Well, that's a problem !!"<<endl;
cout << "No gHaVars ?! Well, that's a problem !!" << endl;
return;
}
if(fDebugFile) *fDebugFile << "THcScalerEvtHandler:: scalerloc size "<<scalerloc.size()<<endl;
if (fDebugFile)
*fDebugFile << "THcScalerEvtHandler:: scalerloc size " << scalerloc.size() << endl;
const Int_t* count = 0;
for (size_t i = 0; i < scalerloc.size(); i++) {
gHaVars->DefineByType(scalerloc[i]->name.Data(), scalerloc[i]->description.Data(),
&dvars[i], kDouble, count);
//gHaVars->DefineByType(scalerloc[i]->name.Data(), scalerloc[i]->description.Data(),
gHaVars->DefineByType(scalerloc[i]->name.Data(), scalerloc[i]->description.Data(), &dvars[i],
kDouble, count);
// gHaVars->DefineByType(scalerloc[i]->name.Data(), scalerloc[i]->description.Data(),
// &dvarsFirst[i], kDouble, count);
}
}
size_t THcScalerEvtHandler::FindNoCase(const string& sdata, const string& skey)
{
size_t THcScalerEvtHandler::FindNoCase(const string& sdata, const string& skey) {
// Find iterator of word "sdata" where "skey" starts. Case insensitive.
string sdatalc, skeylc;
sdatalc = ""; skeylc = "";
for (string::const_iterator p =
sdata.begin(); p != sdata.end(); ++p) {
sdatalc = "";
skeylc = "";
for (string::const_iterator p = sdata.begin(); p != sdata.end(); ++p) {
sdatalc += tolower(*p);
}
for (string::const_iterator p =
skey.begin(); p != skey.end(); ++p) {
for (string::const_iterator p = skey.begin(); p != skey.end(); ++p) {
skeylc += tolower(*p);
}
if (sdatalc.find(skeylc,0) == string::npos) return -1;
return sdatalc.find(skeylc,0);
if (sdatalc.find(skeylc, 0) == string::npos)
return -1;
return sdatalc.find(skeylc, 0);
};
ClassImp(THcScalerEvtHandler)
src/THcScalerEvtHandler.h
View file @ e974befb
......@@ -59,6 +59,8 @@ public:
Int_t fNumBCMs;
Double_t *fBCM_Gain;
Double_t *fBCM_Offset;
Double_t *fBCM_SatOffset;
Double_t *fBCM_SatQuadratic;
Double_t *fBCM_delta_charge;
Double_t fTotalTime;
Double_t fDeltaTime;
......
src/THcScintillatorPlane.cxx
View file @ e974befb
......@@ -534,11 +534,21 @@ Int_t THcScintillatorPlane::DefineVariables( EMode mode )
{"DiffDisTrackCorr", "TW Corrected Dist Difference between track and scintillator position (cm)", "fGoodDiffDistTrack"},
{"TrackXPos", "Track X position at plane (cm)", "fTrackXPosition"},
{"TrackYPos", "Track Y position at plane (cm)", "fTrackYPosition"},
{"ScinXPos", "Scint Average Y position at plane (cm)", "fScinXPos"},
{"ScinYPos", "Scint Average Xposition at plane (cm)", "fScinYPos"},
{"NumClus", "Number of clusters", "fNumberClusters"},
{"Clus.Pos", "Position of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterPosition()"},
{"Clus.Size", "Size of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterSize()"},
{"Clus.Flag", "Flag of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterFlag()"},
{"Clus.UsedFlag", "USed Flag of each paddle clusters", "fCluster.THcScintPlaneCluster.GetClusterUsedFlag()"},
{"PosTdcRefTime", "Reference time of Pos TDC", "fPosTdcRefTime"},
{"NegTdcRefTime", "Reference time of Neg TDC", "fNegTdcRefTime"},
{"PosAdcRefTime", "Reference time of Pos ADC", "fPosAdcRefTime"},
{"NegAdcRefTime", "Reference time of Neg aDC", "fNegAdcRefTime"},
{"PosTdcRefDiffTime", "Reference Diff time of Pos TDC", "fPosTdcRefDiffTime"},
{"NegTdcRefDiffTime", "Reference Diff time of Neg TDC", "fNegTdcRefDiffTime"},
{"PosAdcRefDiffTime", "Reference Diff time of Pos ADC", "fPosAdcRefDiffTime"},
{"NegAdcRefDiffTime", "Reference Diff time of Neg aDC", "fNegAdcRefDiffTime"},
//{"ngoodhits", "Number of paddle hits (passed tof tolerance and used to determine the focal plane time )", "GetNGoodHits() "},
{ 0 }
};
......@@ -643,9 +653,19 @@ void THcScintillatorPlane::Clear( Option_t* )
fpTime = -1.e4;
fHitDistance = kBig;
fScinYPos = kBig;
fScinXPos = kBig;
fTrackXPosition = kBig;
fTrackYPosition = kBig;
fNumberClusters=0;
fPosTdcRefTime = kBig;
fPosAdcRefTime = kBig;
fNegTdcRefTime = kBig;
fNegAdcRefTime = kBig;
fPosTdcRefDiffTime = kBig;
fPosAdcRefDiffTime = kBig;
fNegTdcRefDiffTime = kBig;
fNegAdcRefDiffTime = kBig;
}
//_____________________________________________________________________________
......@@ -700,6 +720,14 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
*
*/
//raw
fPosTdcRefTime = kBig;
fPosAdcRefTime = kBig;
fNegTdcRefTime = kBig;
fNegAdcRefTime = kBig;
fPosTdcRefDiffTime = kBig;
fPosAdcRefDiffTime = kBig;
fNegTdcRefDiffTime = kBig;
fNegAdcRefDiffTime = kBig;
Int_t nrPosTDCHits=0;
Int_t nrNegTDCHits=0;
Int_t nrPosADCHits=0;
......@@ -774,6 +802,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
// might include multiple hits if it uses a multihit tdc.
// Use "ihit" as the index over THcRawHodoHit objects. Use
// "thit" to index over multiple tdc hits within an "ihit".
Bool_t problem_flag=kFALSE;
while(ihit < nrawhits) {
THcRawHodoHit* hit = (THcRawHodoHit *) rawhits->At(ihit);
if(hit->fPlane > fPlaneNum) {
......@@ -784,8 +813,17 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
Int_t index=padnum-1;
THcRawTdcHit& rawPosTdcHit = hit->GetRawTdcHitPos();
if (rawPosTdcHit.GetNHits() >0 && rawPosTdcHit.HasRefTime()) {
if (fPosTdcRefTime == kBig) {
fPosTdcRefTime=rawPosTdcHit.GetRefTime() ;
fPosTdcRefDiffTime=rawPosTdcHit.GetRefDiffTime() ;
}
if (fPosTdcRefTime != rawPosTdcHit.GetRefTime()) {
cout << "THcScintillatorPlane: " << GetName() << " reftime problem at paddle num = " << padnum << " TDC pos hits = " << rawPosTdcHit.GetNHits() << endl;
problem_flag=kTRUE;
}
}
for (UInt_t thit=0; thit<rawPosTdcHit.GetNHits(); ++thit) {
((THcSignalHit*) frPosTdcTimeRaw->ConstructedAt(nrPosTdcHits))->Set(padnum, rawPosTdcHit.GetTimeRaw(thit));
((THcSignalHit*) frPosTdcTime->ConstructedAt(nrPosTdcHits))->Set(padnum, rawPosTdcHit.GetTime(thit));
......@@ -794,6 +832,17 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fTotNumPosTdcHits++;
}
THcRawTdcHit& rawNegTdcHit = hit->GetRawTdcHitNeg();
if (rawNegTdcHit.GetNHits() >0 && rawNegTdcHit.HasRefTime()) {
if (fNegTdcRefTime == kBig) {
fNegTdcRefTime=rawNegTdcHit.GetRefTime() ;
fNegTdcRefDiffTime=rawNegTdcHit.GetRefDiffTime() ;
}
if (fNegTdcRefTime != rawNegTdcHit.GetRefTime()) {
cout << "THcScintillatorPlane: " << GetName()<< " Neg TDC reftime problem at paddle num = " << padnum << " TDC neg hits = " << rawNegTdcHit.GetNHits() << endl;
problem_flag=kTRUE;
}
}
// cout << " paddle num = " << padnum << " TDC Neg hits = " << rawNegTdcHit.GetNHits() << endl;
for (UInt_t thit=0; thit<rawNegTdcHit.GetNHits(); ++thit) {
((THcSignalHit*) frNegTdcTimeRaw->ConstructedAt(nrNegTdcHits))->Set(padnum, rawNegTdcHit.GetTimeRaw(thit));
((THcSignalHit*) frNegTdcTime->ConstructedAt(nrNegTdcHits))->Set(padnum, rawNegTdcHit.GetTime(thit));
......@@ -802,6 +851,17 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fTotNumNegTdcHits++;
}
THcRawAdcHit& rawPosAdcHit = hit->GetRawAdcHitPos();
if (rawPosAdcHit.GetNPulses() >0 && rawPosAdcHit.HasRefTime()) {
if (fPosAdcRefTime == kBig ) {
fPosAdcRefTime=rawPosAdcHit.GetRefTime() ;
fPosAdcRefDiffTime=rawPosAdcHit.GetRefDiffTime() ;
}
if (fPosAdcRefTime != rawPosAdcHit.GetRefTime()) {
cout << "THcScintillatorPlane: " << GetName()<< " Pos ADC reftime problem at paddle num = " << padnum << " ADC pos hits = " << rawPosAdcHit.GetNPulses() << endl;
problem_flag=kTRUE;
}
}
// cout << " paddle num = " << padnum << " ADC Pos hits = " << rawPosAdcHit.GetNPulses() << endl;
for (UInt_t thit=0; thit<rawPosAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPedRaw());
((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPed());
......@@ -824,6 +884,17 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fTotNumPosAdcHits++;
}
THcRawAdcHit& rawNegAdcHit = hit->GetRawAdcHitNeg();
if (rawNegAdcHit.GetNPulses()>0 && rawNegAdcHit.HasRefTime()) {
if (fNegAdcRefTime == kBig) {
fNegAdcRefTime=rawNegAdcHit.GetRefTime() ;
fNegAdcRefDiffTime=rawNegAdcHit.GetRefDiffTime() ;
}
if (fNegAdcRefTime != rawNegAdcHit.GetRefTime()) {
cout << "THcScintillatorPlane: " << GetName()<< " Neg ADC reftime problem at paddle num = " << padnum << " TDC pos hits = " << rawNegAdcHit.GetNPulses() << endl;
problem_flag=kTRUE;
}
}
// cout << " paddle num = " << padnum << " ADC Neg hits = " << rawNegAdcHit.GetNPulses() << endl;
for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPedRaw());
((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPed());
......@@ -882,7 +953,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(tdc_pos >= fScinTdcMin && tdc_pos <= fScinTdcMax) {
btdcraw_pos = kTRUE;
good_ielem_postdc = thit;
break;
break;
}
}
for(UInt_t thit=0; thit<hit->GetRawTdcHitNeg().GetNHits(); thit++) {
......@@ -890,58 +961,86 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if(tdc_neg >= fScinTdcMin && tdc_neg <= fScinTdcMax) {
btdcraw_neg = kTRUE;
good_ielem_negtdc = thit;
break;
break;
}
}
//
if(fADCMode == kADCDynamicPedestal) {
Int_t good_ielem_negadc_test2=-1;
Int_t good_ielem_posadc_test2=-1;
//Loop Here over all hits per event for neg side of plane
if (good_ielem_negtdc != -1) {
Double_t max_adcamp_test=-1000.;
Double_t max_adctdcdiff_test=1000.;
for (UInt_t ielem=0;ielem<rawNegAdcHit.GetNPulses();ielem++) {
Double_t pulsePed = rawNegAdcHit.GetPed();
Double_t pulseInt = rawNegAdcHit.GetPulseInt(ielem);
Double_t pulseAmp = rawNegAdcHit.GetPulseAmp(ielem);
Double_t pulseTime = rawNegAdcHit.GetPulseTime(ielem)+fAdcTdcOffset;
Bool_t errorflag = 0 ;
Double_t TdcAdcTimeDiff = tdc_neg*fScinTdcToTime-pulseTime;
if (rawNegAdcHit.GetPulseAmpRaw(ielem) <= 0) errorflag=1;
if (rawNegAdcHit.GetPulseAmpRaw(ielem) <= 0)pulseAmp= 200.;
Bool_t pulseTimeCut =( TdcAdcTimeDiff > fHodoNegAdcTimeWindowMin[index]) && (TdcAdcTimeDiff < fHodoNegAdcTimeWindowMax[index]);
if (!errorflag && pulseTimeCut && adcint_neg == -999) {
adcped_neg = pulsePed;
adcmult_neg = rawNegAdcHit.GetNPulses();
adchitused_neg = ielem+1;
adcint_neg = pulseInt;
adcamp_neg = pulseAmp;
adctime_neg = pulseTime;
badcraw_neg = kTRUE;
if (pulseTimeCut && pulseAmp>max_adcamp_test) {
good_ielem_negadc = ielem;
adctdcdifftime_neg=TdcAdcTimeDiff;
max_adcamp_test=pulseAmp;
}
if (abs(TdcAdcTimeDiff) < max_adctdcdiff_test) {
good_ielem_negadc_test2 = ielem;
max_adctdcdiff_test=abs(TdcAdcTimeDiff);
}
}
}
//
if ( good_ielem_negadc == -1 && good_ielem_negadc_test2 != -1) good_ielem_negadc=good_ielem_negadc_test2;
if ( good_ielem_negadc == -1 && good_ielem_negadc_test2 == -1 && rawNegAdcHit.GetNPulses()>0) {
good_ielem_negadc=0;
}
//
if (good_ielem_negadc != -1 && good_ielem_negadc<rawNegAdcHit.GetNPulses()) {
adcped_neg = rawNegAdcHit.GetPed();
adcmult_neg = rawNegAdcHit.GetNPulses();
adchitused_neg = good_ielem_negadc+1;
adcint_neg = rawNegAdcHit.GetPulseInt(good_ielem_negadc);
adcamp_neg = rawNegAdcHit.GetPulseAmp(good_ielem_negadc);
if (rawNegAdcHit.GetPulseAmpRaw(good_ielem_negadc) <= 0) adcamp_neg= 200.;
adctime_neg = rawNegAdcHit.GetPulseTime(good_ielem_negadc)+fAdcTdcOffset;
badcraw_neg = kTRUE;
adctdcdifftime_neg=tdc_neg*fScinTdcToTime-adctime_neg;
}
//
//Loop Here over all hits per event for pos side of plane
if (good_ielem_postdc != -1) {
Double_t max_adcamp_test=-1000.;
Double_t max_adctdcdiff_test=1000.;
//
for (UInt_t ielem=0;ielem<rawPosAdcHit.GetNPulses();ielem++) {
Double_t pulsePed = rawPosAdcHit.GetPed();
Double_t pulseInt = rawPosAdcHit.GetPulseInt(ielem);
Double_t pulseAmp = rawPosAdcHit.GetPulseAmp(ielem);
Double_t pulseTime = rawPosAdcHit.GetPulseTime(ielem)+fAdcTdcOffset;
Bool_t errorflag = 0 ;
Double_t TdcAdcTimeDiff = tdc_pos*fScinTdcToTime-pulseTime;
if (rawPosAdcHit.GetPulseAmpRaw(ielem) <= 0) errorflag=1;
Bool_t pulseTimeCut =( TdcAdcTimeDiff > fHodoPosAdcTimeWindowMin[index]) && (TdcAdcTimeDiff < fHodoPosAdcTimeWindowMax[index]);
if (!errorflag && pulseTimeCut && adcint_pos == -999) {
adcped_pos = pulsePed;
adcmult_pos = rawPosAdcHit.GetNPulses();
adchitused_pos = ielem+1;
adcint_pos = pulseInt;
adcamp_pos = pulseAmp;
adctime_pos = pulseTime;
badcraw_pos = kTRUE;
if (rawPosAdcHit.GetPulseAmpRaw(ielem) <= 0)pulseAmp= 200.;
if (pulseTimeCut && pulseAmp>max_adcamp_test) {
good_ielem_posadc = ielem;
adctdcdifftime_pos=TdcAdcTimeDiff;
max_adcamp_test=pulseAmp;
}
if (abs(TdcAdcTimeDiff) < max_adctdcdiff_test) {
good_ielem_posadc_test2 = ielem;
max_adctdcdiff_test=abs(TdcAdcTimeDiff);
}
}
} //
if ( good_ielem_posadc == -1 && good_ielem_posadc_test2 != -1) good_ielem_posadc=good_ielem_posadc_test2;
if ( good_ielem_posadc == -1 && good_ielem_posadc_test2 == -1 && rawPosAdcHit.GetNPulses()>0) good_ielem_posadc=0;
if (good_ielem_posadc != -1 && good_ielem_posadc<rawPosAdcHit.GetNPulses()) {
adcped_pos = rawPosAdcHit.GetPed();
adcmult_pos = rawPosAdcHit.GetNPulses();
adchitused_pos = good_ielem_posadc+1;
adcint_pos = rawPosAdcHit.GetPulseInt(good_ielem_posadc);
adcamp_pos = rawPosAdcHit.GetPulseAmp(good_ielem_posadc);
if (rawPosAdcHit.GetPulseAmpRaw(good_ielem_posadc) <= 0) adcamp_pos= 200.;
adctime_pos = rawPosAdcHit.GetPulseTime(good_ielem_posadc)+fAdcTdcOffset;
badcraw_pos = kTRUE;
adctdcdifftime_pos=tdc_pos*fScinTdcToTime-adctime_pos;
//
}
} else if (fADCMode == kADCSampleIntegral) {
adcint_pos = hit->GetRawAdcHitPos().GetSampleIntRaw() - fPosPed[index];
......@@ -972,7 +1071,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
//
if((btdcraw_pos && badcraw_pos) || (btdcraw_neg && badcraw_neg )) {
if (good_ielem_posadc != -1) {
if (good_ielem_posadc != -1) {
//good adc multiplicity
fTotNumGoodPosAdcHits++;
fTotNumGoodAdcHits++;
......@@ -1045,9 +1144,12 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
}
// Do corrections if valid TDC on both ends of bar
if( (btdcraw_pos && btdcraw_neg) && (badcraw_pos && badcraw_neg) ) {
// Do the pulse height correction to the time. (Position dependent corrections later)
Double_t adc_timec_pos= adctime_pos;
Double_t adc_timec_neg= adctime_neg;
Double_t timec_pos, timec_neg;
if(fTofUsingInvAdc) {
timec_pos = tdc_pos*fScinTdcToTime
......@@ -1059,6 +1161,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else { // FADC style
timec_pos = tdc_pos*fScinTdcToTime -tw_corr_pos + fHodo_LCoeff[index];
timec_neg = tdc_neg*fScinTdcToTime -tw_corr_neg- 2*fHodoCableFit[index] + fHodo_LCoeff[index];
adc_timec_pos = adc_timec_pos -tw_corr_pos + fHodo_LCoeff[index];
adc_timec_neg = adc_timec_neg -tw_corr_neg- 2*fHodoCableFit[index] + fHodo_LCoeff[index];
}
Double_t TWCorrDiff = fGoodNegTdcTimeWalkCorr.at(padnum-1) - 2*fHodoCableFit[index] - fGoodPosTdcTimeWalkCorr.at(padnum-1);
......@@ -1077,6 +1181,8 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
hit_position=TMath::Min(hit_position,fPosLeft);
hit_position=TMath::Max(hit_position,fPosRight);
Double_t scin_corrected_time, postime, negtime;
Double_t adc_postime=adc_timec_pos;
Double_t adc_negtime=adc_timec_neg;
if(fTofUsingInvAdc) {
timec_pos -= (fPosLeft-hit_position)/
fHodoPosInvAdcLinear[index];
......@@ -1091,23 +1197,32 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
negtime = timec_neg - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
}
} else {
scin_corrected_time = 0.5*(timec_neg+timec_pos); // add constants for each paddle, 25ns, 25 + zpos, . . . //remove propagation time
scin_corrected_time = 0.5*(timec_neg+timec_pos);
timec_pos= scin_corrected_time;
timec_neg= scin_corrected_time;
timec_neg= scin_corrected_time;
Double_t adc_time_corrected = 0.5*(adc_timec_pos+adc_timec_neg);
if (fCosmicFlag) {
postime = timec_pos + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
negtime = timec_neg + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
adc_postime = adc_time_corrected + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
adc_negtime = adc_time_corrected + (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
} else {
postime = timec_pos - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
negtime = timec_neg - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
adc_postime = adc_time_corrected - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
adc_negtime = adc_time_corrected - (fZpos+(index%2)*fDzpos)/(29.979*fBetaNominal);
}
}
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPaddleCenter(fPosCenter[index]);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCorrectedTimes(timec_pos,timec_neg,
postime, negtime,
scin_corrected_time);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adcamp_pos); // need for new TWCOrr
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adcamp_neg); // need for new TWCOrr
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adcamp_pos);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adcamp_neg);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCCorrtime(adc_postime);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCCorrtime(adc_negtime);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCalcPosition(fHitDistCorr); //
fGoodPosTdcTimeCorr.at(padnum-1) = timec_pos;
fGoodNegTdcTimeCorr.at(padnum-1) = timec_neg;
fGoodPosTdcTimeTOFCorr.at(padnum-1) = postime;
......@@ -1138,14 +1253,24 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
if (badcraw_neg) adc_neg=adcamp_neg;
if (badcraw_pos) adc_pos=adcamp_pos;
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPaddleCenter(fPosCenter[index]);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCorrectedTimes(timec_pos,timec_neg,
timec_pos,timec_neg,0.0);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCorrectedTimes(timec_pos,timec_neg);
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCpeak(adc_neg); // needed for new TWCOrr
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCpeak(adc_pos); // needed for new TWCOrr
if (badcraw_neg) {
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCtime(adctime_neg);
} else {
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetNegADCtime(-999.);
}
if (badcraw_pos) {
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCtime(adctime_pos);
} else {
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetPosADCtime(-999.);
}
((THcHodoHit*) fHodoHits->At(fNScinHits))->SetCalcPosition(kBig); //
fGoodPosTdcTimeCorr.at(padnum-1) = timec_pos;
fGoodNegTdcTimeCorr.at(padnum-1) = timec_neg;
fGoodPosTdcTimeTOFCorr.at(padnum-1) = timec_pos;
fGoodNegTdcTimeTOFCorr.at(padnum-1) = timec_neg;
fGoodPosTdcTimeTOFCorr.at(padnum-1) = kBig;
fGoodNegTdcTimeTOFCorr.at(padnum-1) = kBig;
}
// if ( ((THcHodoHit*) fHodoHits->At(fNScinHits))->GetPosTOFCorrectedTime() != ((THcHodoHit*) fHodoHits->At(fNScinHits))->GetPosTOFCorrectedTime()) cout << " ihit = " << ihit << " scinhit = " << fNScinHits << " plane = " << fPlaneNum << " padnum = " << padnum << " " << tdc_pos<< " "<< tdc_neg<< " " << ((THcHodoHit*) fHodoHits->At(fNScinHits))->GetPosTOFCorrectedTime() << endl;
fNScinHits++; // One or more good time counter
......@@ -1154,7 +1279,10 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
ihit++; // Raw hit counter
}
// cout << "THcScintillatorPlane: ihit = " << ihit << endl;
if (problem_flag) {
cout << "THcScintillatorPlane::ProcessHits " << fPlaneNum << " " << nexthit << "/" << nrawhits << endl;
cout << " Ref problem end *******" << endl;
}
return(ihit);
}
......
src/THcScintillatorPlane.h
View file @ e974befb
......@@ -63,6 +63,8 @@ class THcScintillatorPlane : public THaSubDetector {
void SetFpTime(Double_t f) {fFptime=f;};
void SetNGoodHits(Int_t ng) {fNGoodHits=ng;};
void SetHitDistance(Double_t f) {fHitDistance=f;}; // Distance between track and hit paddle
void SetScinYPos(Double_t f) {fScinYPos=f;}; // Scint Average Y position at plane (cm)
void SetScinXPos(Double_t f) {fScinXPos=f;}; // Scint Average X position at plane (cm)
void SetTrackXPosition(Double_t f) {fTrackXPosition=f;}; // Distance track X position at plane
void SetTrackYPosition(Double_t f) {fTrackYPosition=f;}; // Distance track Y position at plane
void SetNumberClusters(Int_t nclus) {fNumberClusters=nclus;}; // number of paddle
......@@ -181,7 +183,17 @@ class THcScintillatorPlane : public THaSubDetector {
vector<Double_t> fGoodDiffDistTrack;
Int_t fDebugAdc;
Double_t fPosTdcRefTime;
Double_t fPosAdcRefTime;
Double_t fNegTdcRefTime;
Double_t fNegAdcRefTime;
Double_t fPosTdcRefDiffTime;
Double_t fPosAdcRefDiffTime;
Double_t fNegTdcRefDiffTime;
Double_t fNegAdcRefDiffTime;
Double_t fHitDistance;
Double_t fScinXPos;
Double_t fScinYPos;
Double_t fTrackXPosition;
Double_t fTrackYPosition;
Int_t fCosmicFlag; //
......@@ -272,7 +284,10 @@ class THcScintillatorPlane : public THaSubDetector {
Double_t *fNegPed;
Double_t *fNegSig;
Double_t *fNegThresh;
//
Int_t fEventType;
Int_t fEventNum;
//
Int_t fNScinGoodHits; // number of hits for which both ends of the paddle fired in time!
Double_t fpTime; // the original code only has one fpTime per plane!
......