Skip to content
Snippets Groups Projects
Unverified Commit e131a4d9 authored by Eric Pooser's avatar Eric Pooser Committed by GitHub
Browse files

Hodo time walk calib (#462) 

* Initial commit of time-walk calibration script

* Updates for safe keeping

* Monior tweaks

* Minor edits

* Add time walk function

* Minor tweaks to histo macro, first hack at fit macro

* Add functions for readability

* Added more functionality

* Still hacking out min/max algo

* More tweaks

* Modify cuts for more recent (12/5/17) cosmic running

* Add time walk scripts for HMS

* Update fit ranges and histogram ranges
parent d15e91b9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CALIBRATION/hms_hodo_calib/timeWalkCalib.C 0 → 100644
+ 296
− 0
View file @ e131a4d9
// Macro to perform time-walk fits and extract the calibration parameters
// Author: Eric Pooser, pooser@jlab.org
// Declare ROOT files
TFile *histoFile, *outFile;
// Declare constants
static const UInt_t nPlanes = 4;
static const UInt_t nSides = 2;
static const UInt_t nBarsMax = 21;
static const UInt_t nTwFitPars = 3;
static const Double_t tdcThresh = 120.0; // 30 mV in units of FADC channels
static const Double_t twFitRangeLow = 10.0;
static const Double_t twFitRangeHigh = 600.0;
static const Double_t c0twParInit = -15.0;
static const Double_t c1twParInit = 5.0;
static const Double_t c2twParInit = 0.25;
static const Double_t fontSize = 0.05;
static const Double_t yTitleOffset = 0.75;
static const Double_t markerSize = 2.0;
static const Double_t minScale = 0.75;
static const Double_t maxScale = 0.75;
static const UInt_t lineWidth = 4;
static const UInt_t lineStyle = 7;
static const UInt_t nbars[nPlanes] = {16, 10, 16, 10};
static const TString planeNames[nPlanes] = {"1x", "1y", "2x", "2y"};
static const TString sideNames[nSides] = {"pos", "neg"};
static const TString twFitParNames[nTwFitPars] = {"c_{0}", "c_{1}", "c_{2}"};
// Declare directories
TDirectory *dataDir, *planeDir[nPlanes], *sideDir[nPlanes][nSides];
TDirectory *twDir[nPlanes][nSides];
// Declare fits
TF1 *twFit[nPlanes][nSides][nBarsMax], *avgParFit[nPlanes][nSides][nTwFitPars];
// Declare arrays
Double_t paddleIndex[nPlanes][nSides][nBarsMax];
Double_t twFitPar[nPlanes][nSides][nTwFitPars][nBarsMax], twFitParErr[nPlanes][nSides][nTwFitPars][nBarsMax];
Double_t avgPar[nPlanes][nSides][nTwFitPars];
Double_t minPar[nPlanes][nSides][nTwFitPars], maxPar[nPlanes][nSides][nTwFitPars];
// Declare canvases
TCanvas *twFitCan[nPlanes][nSides], *twFitParCan[nTwFitPars];
// Declare histos
// 2D histos
TH2F *h2_adcTdcTimeDiffWalk[nPlanes][nSides][nBarsMax];
// Declare graphs
TGraph *twFitParGraph[nPlanes][nSides][nTwFitPars];
// Declare multi-graphs
TMultiGraph *twFitParMultiGraph[nPlanes][nTwFitPars];
// Declare legends
TLegend *twFitParLegend[nPlanes][nTwFitPars], *avgParLegend[nPlanes][nTwFitPars];
// Declare text objects
TPaveText *twFitParText[nPlanes][nSides][nBarsMax];
// Functions to extract calibration parameters
//=:=:=:=:=:=:
//=: Level 1
//=:=:=:=:=:=:
// Add color to fit lines
void addColorToFitLine(UInt_t style, UInt_t width, UInt_t color, TF1 *fit1D) {
fit1D->SetLineStyle(lineStyle);
fit1D->SetLineWidth(lineWidth);
fit1D->SetLineColor(color);
return;
} // addColorToFitLine()
// Add color to fit graphs
void addColorToGraph(UInt_t style, UInt_t size, UInt_t color, TGraph *graph) {
graph->SetMarkerStyle(style);
graph->SetMarkerSize(size);
graph->SetMarkerColor(color);
return;
} // addColorToFitGraph()
// Add centered titles to multigraphs
void modifyMultiGraph(TMultiGraph *mg, TString xtitle, TString ytitle) {
mg->GetYaxis()->SetTitle(ytitle);
mg->GetYaxis()->CenterTitle();
mg->GetXaxis()->SetTitle(xtitle);
mg->GetXaxis()->CenterTitle();
return;
} // modifyMultiGraph()
// Create canvases
TCanvas *makeCan(UInt_t numColumns, UInt_t numRows, UInt_t winWidth, UInt_t winHeight, TCanvas *can, TString name, TString title) {
can = new TCanvas(name, title, winWidth, winHeight);
can->Divide(numColumns, numRows);
return can;
} // makeCan()
// Define the time-walk fit function
Double_t twFitFunc(Double_t *a, Double_t *c) {
Double_t twFitVal = c[0] + c[1]/(TMath::Power((a[0]/tdcThresh), c[2]));
return twFitVal;
} // twFitFunc()
// Locate min or max value from input array
Double_t calcMinOrMax(Double_t *array, UInt_t iplane, TString minOrmax) {
auto result = minmax_element(array, array+nbars[iplane]);
if (minOrmax == "min") return *result.first;
else if (minOrmax == "max") return *result.second;
else return 0.0;
} // calcMinOrMax()
//=:=:=:=:=:=:
//=: Level 2
//=:=:=:=:=:=:
// Perform the timw-walk fits
void doTwFits(UInt_t iplane, UInt_t iside, UInt_t ipaddle) {
// Draw fits on canvas
twFitCan[iplane][iside]->cd(ipaddle+1);
gPad->SetLogz();
h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->Draw("COLZ");
gPad->Modified(); gPad->Update();
// Declare and initialize the fits
twFit[iplane][iside][ipaddle] = new TF1("twFit", twFitFunc, twFitRangeLow, twFitRangeHigh, nTwFitPars);
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFit[iplane][iside][ipaddle]->SetParName(ipar, twFitParNames[ipar]);
twFit[iplane][iside][ipaddle]->SetParameters(c0twParInit, c1twParInit, c2twParInit);
// Perform the fits and scream if it failed
Int_t twFitStatus = h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->Fit("twFit", "REQ");
if (twFitStatus != 0)
cout << "ERROR: Time Walk Fit Failed!!! " << "Status = " << twFitStatus << " For Plane: " << planeNames[iplane] << " Side: " << sideNames[iside] << " Paddle: " << ipaddle+1 << endl;
// Create text box to display fir parameters
twFitParText[iplane][iside][ipaddle] = new TPaveText(0.4, 0.6, 0.895, 0.895, "NBNDC");
twFitParText[iplane][iside][ipaddle]->AddText(Form("Entries = %.0f", h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->GetEntries()));
// Obtain the fit parameters and associated errors
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
twFitPar[iplane][iside][ipar][ipaddle] = twFit[iplane][iside][ipaddle]->GetParameter(ipar);
twFitParErr[iplane][iside][ipar][ipaddle] = twFit[iplane][iside][ipaddle]->GetParError(ipar);
twFitParText[iplane][iside][ipaddle]->AddText(Form(twFitParNames[ipar]+" = %.2f #pm %.2f", twFitPar[iplane][iside][ipar][ipaddle], twFitParErr[iplane][iside][ipar][ipaddle]));
} // Parameter loop
// Draw the fit parameter text
twFitParText[iplane][iside][ipaddle]->SetFillColor(kWhite);
twFitParText[iplane][iside][ipaddle]->SetTextAlign(12);
twFitParText[iplane][iside][ipaddle]->SetFillStyle(3003);
twFitParText[iplane][iside][ipaddle]->Draw();
gPad->Modified(); gPad->Update();
return;
} // doTwFits()
// Calculate the averege of the time-walk fit parameters
void calcParAvg(UInt_t iplane, UInt_t iside) {
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
// Calculate the weighted average while ignoring fit errors provided by Minuit
avgParFit[iplane][iside][ipar] = new TF1("avgParFit", "pol0", 1, nbars[iplane]);
avgParFit[iplane][iside][ipar]->SetParName(0, "#bar{"+twFitParNames[ipar]+"}");
// Add color to fit lines
if (iside == 0) addColorToFitLine(lineStyle, lineWidth, kRed, avgParFit[iplane][iside][ipar]);
if (iside == 1) addColorToFitLine(lineStyle, lineWidth, kBlue, avgParFit[iplane][iside][ipar]);
// Initialize the parameters
if (ipar == 0) avgParFit[iplane][iside][ipar]->SetParameter(0, c0twParInit);
if (ipar == 1) avgParFit[iplane][iside][ipar]->SetParameter(0, c1twParInit);
if (ipar == 2) avgParFit[iplane][iside][ipar]->SetParameter(0, c2twParInit);
// Perform the least squares fit
Int_t avgParFitStatus = twFitParGraph[iplane][iside][ipar]->Fit("avgParFit", "REQ");
if (avgParFitStatus != 0)
cout << "ERROR: Parameter Fit Failed!!! " << "Status = " << avgParFitStatus << " For Plane: " << planeNames[iplane] << " Side: " << sideNames[iside] << endl;
// Store the fit parameters
avgPar[iplane][iside][ipar] = avgParFit[iplane][iside][ipar]->GetParameter(0);
// Add graphs to multi graph
twFitParMultiGraph[iplane][ipar]->Add(twFitParGraph[iplane][iside][ipar]);
// Calculate the min and max value of each parameter
minPar[iplane][iside][ipar] = calcMinOrMax(twFitPar[iplane][iside][ipar], iplane, "min");
maxPar[iplane][iside][ipar] = calcMinOrMax(twFitPar[iplane][iside][ipar], iplane, "max");
} // Parameter loop
return;
} // calcParAvg()
// Draw the time-walk fit parameter multi-graph
void drawParams(UInt_t iplane) {
// Populate the parameter canvases
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
twFitParCan[ipar]->cd(iplane+1);
twFitParMultiGraph[iplane][ipar]->Draw("AP");
gPad->Modified(); gPad->Update();
// Add aesthetics
modifyMultiGraph(twFitParMultiGraph[iplane][ipar], "Paddle Number", "Fit Parameter "+twFitParNames[ipar]);
twFitParMultiGraph[iplane][ipar]->GetXaxis()->SetNdivisions(21);
twFitParMultiGraph[iplane][ipar]->Draw("AP");
// Draw legends
twFitParLegend[iplane][ipar] = new TLegend(0.1, 0.7, 0.3, 0.9);
avgParLegend[iplane][ipar] = new TLegend(0.6, 0.7, 0.9, 0.9);
for(UInt_t iside = 0; iside < nSides; iside++) {
if (iside == 0) twFitParLegend[iplane][ipar]->AddEntry(twFitParGraph[iplane][iside][ipar], "Positive Side", "p");
if (iside == 1) twFitParLegend[iplane][ipar]->AddEntry(twFitParGraph[iplane][iside][ipar], "Negative Side", "p");
avgParLegend[iplane][ipar]->AddEntry(avgParFit[iplane][iside][ipar], "#bar{"+twFitParNames[ipar]+"}"+Form(" = %.2f", avgPar[iplane][iside][ipar]),"l");
// Define the time-walk fit parameter multigraph y-axis range
if (iside == 0) {
// Set y-axis range min
if (minPar[iplane][iside][ipar] < 0.0 || minPar[iplane][iside+1][ipar] < 0.0) {
if (minPar[iplane][iside][ipar] < minPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside][ipar] + minScale*minPar[iplane][iside][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside+1][ipar] + minScale*minPar[iplane][iside+1][ipar]);
}
if (minPar[iplane][iside][ipar] > 0.0 || minPar[iplane][iside+1][ipar] > 0.0) {
if (minPar[iplane][iside][ipar] < minPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside][ipar] - minScale*minPar[iplane][iside][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside+1][ipar] - minScale*minPar[iplane][iside+1][ipar]);
}
// Set y-axis range max
if (maxPar[iplane][iside][ipar] < 0.0 || maxPar[iplane][iside+1][ipar] < 0.0) {
if (maxPar[iplane][iside][ipar] < maxPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside+1][ipar] - maxScale*maxPar[iplane][iside+1][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside][ipar] - maxScale*maxPar[iplane][iside][ipar]);
}
if (maxPar[iplane][iside][ipar] > 0.0 || maxPar[iplane][iside+1][ipar] > 0.0) {
if (maxPar[iplane][iside][ipar] < maxPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside+1][ipar] + maxScale*maxPar[iplane][iside+1][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside][ipar] + maxScale*maxPar[iplane][iside][ipar]);
}
} // Side index = 0 condition
} // Side loop
twFitParLegend[iplane][ipar]->Draw();
avgParLegend[iplane][ipar]->Draw();
} // Parameter loop
return;
} // drawParams
//=:=:=:=:=
//=: Main
//=:=:=:=:=
void timeWalkCalib() {
// ROOT settings
gStyle->SetTitleFontSize(fontSize);
gStyle->SetLabelSize(fontSize, "XY");
gStyle->SetTitleSize(fontSize, "XY");
gStyle->SetTitleOffset(yTitleOffset, "Y");
gStyle->SetStatFormat(".2f");
gStyle->SetOptFit(0);
gStyle->SetOptStat(0);
// Read the ROOT file containing the time-walk histos
histoFile = new TFile("timeWalkHistos.root", "READ");
// Obtain the top level directory
dataDir = dynamic_cast <TDirectory*> (histoFile->FindObjectAny("hodoUncalib"));
// Create the parameter canvases
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFitParCan[ipar] = makeCan(2, 2, 1600, 800, twFitParCan[ipar], twFitParNames[ipar]+"FitParCan", "Parameter "+twFitParNames[ipar]+" Canvas");
// Loop over the planes
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
//for(UInt_t iplane = 0; iplane < 1; iplane++) {
// Obtain the plane directory
planeDir[iplane] = dynamic_cast <TDirectory*> (dataDir->FindObjectAny(planeNames[iplane]));
// Create multigraphs
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFitParMultiGraph[iplane][ipar] = new TMultiGraph(planeNames[iplane]+"_"+twFitParNames[ipar]+"_Multigraph", "Plane "+planeNames[iplane]+" Parameter "+twFitParNames[ipar]);
// Loop over the sides
for(UInt_t iside = 0; iside < nSides; iside++) {
// Obtain the side and time walk directories
sideDir[iplane][iside] = dynamic_cast <TDirectory*> (planeDir[iplane]->FindObjectAny(sideNames[iside]));
twDir[iplane][iside] = dynamic_cast <TDirectory*> (sideDir[iplane][iside]->FindObjectAny("adcTdcTimeDiffWalk"));
// Create the time-walk histo and fit canvases
if (planeNames[iplane] != "1y" || planeNames[iplane] != "2y") twFitCan[iplane][iside] = makeCan(4, 4, 1600, 800, twFitCan[iplane][iside], planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan", planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan");
if (planeNames[iplane] == "1y" || planeNames[iplane] == "2y") twFitCan[iplane][iside] = makeCan(4, 3, 1600, 800, twFitCan[iplane][iside], planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan", planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan");
// Loop over the paddles
for(UInt_t ipaddle = 0; ipaddle < nbars[iplane]; ipaddle++) {
// Populate the paddle index arrays
paddleIndex[iplane][iside][ipaddle] = Double_t (ipaddle + 1);
// Obtain the time-walk histos
h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle] = dynamic_cast <TH2F*> (twDir[iplane][iside]->FindObjectAny(Form("h2_adcTdcTimeDiffWalk_paddle_%d", ipaddle+1)));
// Perform the time-walk fits while ignoring the S2Y plane (for now)
// if (planeNames[iplane] == "2y") continue;
doTwFits(iplane, iside, ipaddle);
} // Paddle loop
// Produce the time-walk fit parameter graphs
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
// Populate graphs and multi-graphs
twFitParGraph[iplane][iside][ipar] = new TGraph(nbars[iplane], paddleIndex[iplane][iside], twFitPar[iplane][iside][ipar]);
if (iside == 0) addColorToGraph(22, markerSize, kRed, twFitParGraph[iplane][iside][ipar]);
if (iside == 1) addColorToGraph(23, markerSize, kBlue, twFitParGraph[iplane][iside][ipar]);
// // Add graphs to multi graph
// twFitParMultiGraph[iplane][ipar]->Add(twFitParGraph[iplane][iside][ipar]);
} // Parameter loop
// Calculate the average of the time-walk fit parameters
calcParAvg(iplane, iside);
} // Side loop
// Draw the time-walk parameter graphs
drawParams(iplane);
} // Plane loop
return;
} // timeWalkFits()
CALIBRATION/hms_hodo_calib/timeWalkHistos.C 0 → 100644
+ 427
− 0
View file @ e131a4d9
// Macro to generate time-walk histograms for the hodoscopes
// Author: Eric Pooser, pooser@jlab.org
#include <time.h>
// Declare replay data file and output file
TFile *replayFile, *outFile;
// Declare data tree
TTree *rawDataTree;
// Declare constants
static const UInt_t nPlanes = 4;
static const UInt_t nSides = 2;
static const UInt_t nAdcSignals = 3;
static const UInt_t nTdcSignals = 1;
static const UInt_t nSignals = nAdcSignals + nTdcSignals;
static const UInt_t maxTdcHits = 128;
static const UInt_t maxAdcHits = 4;
static const UInt_t nMaxBars = 16;
static const TString spec = "H";
static const TString detec = "hod";
static const TString trig = "T";
static const TString hms = "hms";
static const UInt_t nbars[nPlanes] = {16, 10, 16, 10};
static const TString planeNames[nPlanes] = {"1x", "1y", "2x", "2y"};
static const TString sideNames[nSides] = {"pos", "neg"};
static const TString signalNames[nSignals] = {"Adc", "Tdc"};
static const TString adcData[nAdcSignals] = {"ErrorFlag", "PulseTimeRaw", "PulseAmp"};
static const TString tdcData[nTdcSignals] = {"TimeRaw"};
static const Double_t tdcChanToTime = 0.100; // Units of ns
static const Double_t adcChanToTime = 0.0625; // Units of ns
static const Double_t adcDynamicRange = 1000.0; // Units of mV
static const Double_t nAdcChan = 4096.0; // Units of ADC channels
static const Double_t adcChanTomV = adcDynamicRange/nAdcChan; // Units of mV/ADC Chan
static const Double_t hodoPulseAmpCutLow = 10.0; // Units of mV
static const Double_t hodoPulseAmpCutHigh = 1000.0; // Units of mV
static const Double_t refAdcPulseAmpCutLow = 50.0; // Units of mV
static const Double_t refAdcPulseAmpCutHigh = 60.0; // Units of mV
static const Double_t refAdcPulseTimeCutLow = 210.0; // Units of ns
static const Double_t refAdcPulseTimeCutHigh = 225.0; // Units of ns
static const Double_t adcTdcTimeDiffCutLow = -100.0; // Units of ns
static const Double_t adcTdcTimeDiffCutHigh = 100.0; // Units of ns
static const Double_t calEtotCutVal = 0.100; // Units of GeV
static const Double_t cerNpeSumCutVal = 1.5; // Units of NPE
// static const Double_t adcTdcTimeDiffCutLow = -6000.0; // Units of ns
// static const Double_t adcTdcTimeDiffCutHigh = 1000.0; // Units of ns
static const Double_t fontSize = 0.04;
// Declare variables to obtain from data tree
// Number of ADC & TDC hits
Int_t adcHits[nPlanes][nSides][nSignals], tdcHits[nPlanes][nSides][nSignals];
// Paddle number which fired
Double_t adcPaddle[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t tdcPaddle[nPlanes][nSides][nSignals][nMaxBars*maxTdcHits];
// FADC data
Double_t hodoAdcPulseTimeRaw[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t hodoAdcPulseAmp[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t hodoAdcErrorFlag[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
// TDC Data
Double_t hodoTdcTimeRaw[nPlanes][nSides][nSignals][nMaxBars*maxTdcHits];
// Trigger apparatus data
Double_t refAdcPulseTimeRaw, refAdcPulseAmp, refAdcMultiplicity;
Double_t refT1TdcTimeRaw, refT2TdcTimeRaw;
// Declare variables to steer data to arrays
TString *adcBaseName, *adcNdataName, *adcPaddleName;
TString *adcErrorFlagName, *adcPulseTimeRawName, *adcPulseAmpName;
TString *tdcBaseName, *tdcNdataName, *tdcPaddleName;
TString *tdcTimeRawName;
// Declare directories for output ROOT file
TDirectory *trigRawDir, *hodoRawDir, *hodoUncalibDir;
TDirectory *planeRawDir[nPlanes], *sideRawDir[nPlanes][nSides], *planeUncalibDir[nPlanes], *sideUncalibDir[nPlanes][nSides];
TDirectory *adcTimeWalkDir[nPlanes][nSides], *tdcTimeWalkDir[nPlanes][nSides], *adcTdcTimeDiffWalkDir[nPlanes][nSides];
// Declare histos
// 1D histos
TH1F *h1_refAdcPulseTimeRaw, *h1_refAdcPulseAmp, *h1_refAdcMultiplicity;
TH1F *h1_refT1TdcTimeRaw, *h1_refT2TdcTimeRaw;
// 2D histos
TH2F *h2_adcErrorFlag[nPlanes][nSides];
TH2F *h2_adcPulseTimeRaw[nPlanes][nSides];
TH2F *h2_adcPulseAmp[nPlanes][nSides];
TH2F *h2_tdcTimeRaw[nPlanes][nSides];
TH2F *h2_adcPulseTime[nPlanes][nSides];
TH2F *h2_tdcTime[nPlanes][nSides];
TH2F *h2_adcTdcTimeDiff[nPlanes][nSides];
TH2F *h2_adcPulseAmpCuts[nPlanes][nSides];
TH2F *h2_adcTimeWalk[nPlanes][nSides][nMaxBars];
TH2F *h2_tdcTimeWalk[nPlanes][nSides][nMaxBars];
TH2F *h2_adcTdcTimeDiffWalk[nPlanes][nSides][nMaxBars];
// Declare event data objects
UInt_t nentries, ievent, adcPaddleNum, tdcPaddleNum;
Int_t numAdcHits, numTdcHits;
Double_t adcErrorFlag, adcPulseTimeRaw, adcPulseTime, adcPulseAmp;
Double_t tdcTimeRaw, tdcTime, adcTdcTimeDiff;
Double_t calEtot, cerNpeSum;
Bool_t adcRefMultiplicityCut, adcRefPulseAmpCut, adcRefPulseTimeCut;
Bool_t edtmCut, adcErrorFlagCut, adcAndTdcHitCut;
Bool_t adcPulseAmpCut, adcTdcTimeDiffCut;
Bool_t calEtotCut, cerNpeSumCut;
void generatePlots(UInt_t iplane, UInt_t iside, UInt_t ipaddle) {
// Create trigger aparatus directory
trigRawDir = dynamic_cast <TDirectory*> (outFile->Get("trigAppRaw"));
if (!trigRawDir) {trigRawDir = outFile->mkdir("trigAppRaw"); trigRawDir->cd();}
else outFile->cd("trigAppRaw");
// FADC reference
if (!h1_refAdcPulseTimeRaw) h1_refAdcPulseTimeRaw = new TH1F("h1_refAdcPulseTimeRaw", "ROC1 Raw FADC Reference Pulse Time; Raw FADC Pulse Time (ns); Number of Entries / 100 ps", 4000, 0, 400);
if (!h1_refAdcPulseAmp) h1_refAdcPulseAmp = new TH1F("h1_refAdcPulseAmp", "ROC1 FADC Reference Pulse Amplitude; FADC Pulse Amplitude (mV); Number of Entries / 1 mV", 1000, 0, 1000);
if (!h1_refAdcMultiplicity) h1_refAdcMultiplicity = new TH1F("h1_refAdcMultiplicity", "ROC1 FADC Reference Multiplicity; Raw FADC Multiplicity; Number of Entries", 5, -0.5, 4.5);
// TDC reference
if (!h1_refT1TdcTimeRaw) h1_refT1TdcTimeRaw = new TH1F("h1_refT1TdcTimeRaw", "ROC1 T1 Raw TDC Reference TDC Time (Slot 20, Channel 15); Raw TDC Time (ns); Number of Entries / 100 ps", 6000, 0, 600);
if (!h1_refT2TdcTimeRaw) h1_refT2TdcTimeRaw = new TH1F("h1_refT2TdcTimeRaw", "ROC1 T2 Raw TDC Reference TDC Time (Slot 19, Channel 31); Raw TDC Time (ns); Number of Entries / 100 ps", 6000, 0, 600);
// Create plane directory for raw hodoscope quantities
hodoRawDir = dynamic_cast <TDirectory*> (outFile->Get("hodoRaw"));
if (!hodoRawDir) {hodoRawDir = outFile->mkdir("hodoRaw"); hodoRawDir->cd();}
else outFile->cd("hodoRaw");
planeRawDir[iplane] = dynamic_cast <TDirectory*> (hodoRawDir->Get(planeNames[iplane]));
if (!planeRawDir[iplane]) {planeRawDir[iplane] = hodoRawDir->mkdir(planeNames[iplane]); planeRawDir[iplane]->cd();}
else hodoRawDir->cd(planeNames[iplane]);
// Create side directory
sideRawDir[iplane][iside] = dynamic_cast <TDirectory*> (planeRawDir[iplane]->Get(sideNames[iside]));
if (!sideRawDir[iplane][iside]) {sideRawDir[iplane][iside] = planeRawDir[iplane]->mkdir(sideNames[iside]); sideRawDir[iplane][iside]->cd();}
else outFile->cd("hodoRaw/"+planeNames[iplane]+"/"+sideNames[iside]);
// Book histos
if (!h2_adcErrorFlag[iplane][iside]) h2_adcErrorFlag[iplane][iside] = new TH2F("h2_adcErrorFlag", "FADC Error Flag Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Error Flag", nbars[iplane], 0.5, nbars[iplane]+0.5, 2, -0.5, 1.5);
if (!h2_adcPulseTimeRaw[iplane][iside]) h2_adcPulseTimeRaw[iplane][iside] = new TH2F("h2_adcPulseTimeRaw", "Raw FADC Pulse Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; Raw FADC Pulse Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 4000, 0, 400);
if (!h2_adcPulseAmp[iplane][iside]) h2_adcPulseAmp[iplane][iside] = new TH2F("h2_adcPulseAmp", "FADC Pulse Amplitude Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Amplitude (mV) / 1 mV", nbars[iplane], 0.5, nbars[iplane]+0.5, 1000, 0, 1000);
if (!h2_tdcTimeRaw[iplane][iside]) h2_tdcTimeRaw[iplane][iside] = new TH2F("h2_tdcTimeRaw", "Raw TDC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; Raw TDC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 4000, 0, 400);
// Create plane directory for uncalibrated hodoscope quantities
hodoUncalibDir = dynamic_cast <TDirectory*> (outFile->Get("hodoUncalib"));
if (!hodoUncalibDir) {hodoUncalibDir = outFile->mkdir("hodoUncalib"); hodoUncalibDir->cd();}
else outFile->cd("hodoUncalib");
planeUncalibDir[iplane] = dynamic_cast <TDirectory*> (hodoUncalibDir->Get(planeNames[iplane]));
if (!planeUncalibDir[iplane]) {planeUncalibDir[iplane] = hodoUncalibDir->mkdir(planeNames[iplane]); planeUncalibDir[iplane]->cd();}
else hodoUncalibDir->cd(planeNames[iplane]);
// Create side directory
sideUncalibDir[iplane][iside] = dynamic_cast <TDirectory*> (planeUncalibDir[iplane]->Get(sideNames[iside]));
if (!sideUncalibDir[iplane][iside]) {sideUncalibDir[iplane][iside] = planeUncalibDir[iplane]->mkdir(sideNames[iside]); sideUncalibDir[iplane][iside]->cd();}
else outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]);
// Book histos
if (!h2_adcPulseTime[iplane][iside]) h2_adcPulseTime[iplane][iside] = new TH2F("h2_adcPulseTime", "FADC Pulse Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 8000, -400, 400);
if (!h2_tdcTime[iplane][iside]) h2_tdcTime[iplane][iside] = new TH2F("h2_tdcTime", "TDC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; TDC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 8000, -400, 400);
if (!h2_adcTdcTimeDiff[iplane][iside]) h2_adcTdcTimeDiff[iplane][iside] = new TH2F("h2_adcTdcTimeDiff", "TDC-ADC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; TDC - ADC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 16000, -800, 800);
if (!h2_adcPulseAmpCuts[iplane][iside]) h2_adcPulseAmpCuts[iplane][iside] = new TH2F("h2_adcPulseAmpCuts", "FADC Pulse Amplitude Cuts Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Amplitude (mV) / 1 mV", nbars[iplane], 0.5, nbars[iplane]+0.5, 1000, 0, 1000);
// Create ADC time-walk directory
adcTimeWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("adcTimeWalk"));
if (!adcTimeWalkDir[iplane][iside]) {adcTimeWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("adcTimeWalk"); adcTimeWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/adcTimeWalk"));
// Book histos
if (!h2_adcTimeWalk[iplane][iside][ipaddle]) h2_adcTimeWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_adcTimeWalk_paddle_%d", ipaddle+1), "ADC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; ADC Time (ns) / 100 ps", 1000, 0, 1000, 4000, -200, 200);
// Create TDC time-walk directory
tdcTimeWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("tdcTimeWalk"));
if (!tdcTimeWalkDir[iplane][iside]) {tdcTimeWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("tdcTimeWalk"); tdcTimeWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/tdcTimeWalk"));
// Book histos
if (!h2_tdcTimeWalk[iplane][iside][ipaddle]) h2_tdcTimeWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_tdcTimeWalk_paddle_%d", ipaddle+1), "TDC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; TDC Time (ns) / 100 ps", 1000, 0, 1000, 4000, -300, 100);
// Create TDC-ADC time-walk directory
adcTdcTimeDiffWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("adcTdcTimeDiffWalk"));
if (!adcTdcTimeDiffWalkDir[iplane][iside]) {adcTdcTimeDiffWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("adcTdcTimeDiffWalk"); adcTdcTimeDiffWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/adcTdcTimeDiffWalk"));
// Book histos
if (!h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]) h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_adcTdcTimeDiffWalk_paddle_%d", ipaddle+1), "TDC-ADC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; TDC-ADC Time (ns) / 100 ps", 1000, 0, 1000, 150, -20, -5);
} // generatePlots()
void timeWalkHistos(Int_t runNum) {
// Global ROOT settings
gStyle->SetOptFit();
gStyle->SetStatFormat(".2f");
gStyle->SetTitleFontSize(fontSize);
gStyle->SetLabelSize(fontSize, "XY");
gStyle->SetTitleSize(fontSize, "XY");
// Initialize the analysis clock
clock_t t;
t = clock();
// Obtain the replay data file and create new output ROOT file
// replayFile = new TFile("ROOTfiles/hms_replay_production_all_1267_-1.root", "READ");
// replayFile = new TFile("ROOTfiles/hms_replay_production_all_1268_-1.root", "READ");
replayFile = new TFile("ROOTfiles/hms_replay_production_all_1267_1268_-1.root", "READ");
// replayFile = new TFile(Form("ROOTfiles/hms_replay_production_all_%d_-1.root", runNum), "READ");
// replayFile = new TFile(Form("ROOTfiles/hms_coin_replay_production_%d_-1.root", runNum), "READ");
outFile = new TFile("timeWalkHistos_temp.root", "RECREATE");
// Obtain the tree
rawDataTree = dynamic_cast <TTree*> (replayFile->Get("T"));
// Acquire the trigger apparatus data
rawDataTree->SetBranchAddress("T.hms.hFADC_TREF_ROC1_adcPulseTimeRaw", &refAdcPulseTimeRaw);
rawDataTree->SetBranchAddress("T.hms.hFADC_TREF_ROC1_adcPulseAmp", &refAdcPulseAmp);
rawDataTree->SetBranchAddress("T.hms.hFADC_TREF_ROC1_adcMultiplicity", &refAdcMultiplicity);
rawDataTree->SetBranchAddress("T.hms.hT1_tdcTimeRaw", &refT1TdcTimeRaw);
rawDataTree->SetBranchAddress("T.hms.hT2_tdcTimeRaw", &refT2TdcTimeRaw);
rawDataTree->SetBranchAddress("H.cal.etot", &calEtot);
rawDataTree->SetBranchAddress("H.cer.npeSum", &cerNpeSum);
// Loop over the planes, sides, signals, leafs, and fill data arrays
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
for(UInt_t iside = 0; iside < nSides; iside++) {
// Generate directory structure and histograms
for(UInt_t ipaddle = 0; ipaddle < nbars[iplane]; ipaddle++)
generatePlots(iplane, iside, ipaddle);
for(UInt_t isignal = 0; isignal < nSignals; isignal++) {
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Define the leaf strings of interest
adcBaseName = new TString(spec+"."+detec+"."+planeNames[iplane]+"."+
sideNames[iside]+signalNames[isignal]);
adcNdataName = new TString("Ndata."+*adcBaseName+"Counter");
adcPaddleName = new TString(*adcBaseName+"Counter");
// Fill the leafs
rawDataTree->SetBranchAddress(*adcNdataName, &adcHits[iplane][iside][isignal]);
rawDataTree->SetBranchAddress(*adcPaddleName, &adcPaddle[iplane][iside][isignal]);
// Loop over the adc data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
if (adcData[iadcdata] == "ErrorFlag") {
adcErrorFlagName = new TString(*adcBaseName+"ErrorFlag");
rawDataTree->SetBranchAddress(*adcErrorFlagName, &hodoAdcErrorFlag[iplane][iside][isignal]);
} // Error flag leaf
if (adcData[iadcdata] == "PulseTimeRaw") {
adcPulseTimeRawName = new TString(*adcBaseName+"PulseTimeRaw");
rawDataTree->SetBranchAddress(*adcPulseTimeRawName, &hodoAdcPulseTimeRaw[iplane][iside][isignal]);
} // Raw pulse time leaf
if (adcData[iadcdata] == "PulseAmp") {
adcPulseAmpName = new TString(*adcBaseName+"PulseAmp");
rawDataTree->SetBranchAddress(*adcPulseAmpName, &hodoAdcPulseAmp[iplane][iside][isignal]);
} // Pedestal subtracted pulse amplitude leaf
} // ADC signal data loop
} // ADC signal
// Acquire the hodoscope TDC data objects
if(signalNames[isignal] == "Tdc") {
// Define the leaf strings of interest
tdcBaseName = new TString(spec+"."+detec+"."+planeNames[iplane]+"."+
sideNames[iside]+signalNames[isignal]);
tdcNdataName = new TString("Ndata."+*tdcBaseName+"Counter");
tdcPaddleName = new TString(*tdcBaseName+"Counter");
// Fill the leafs
rawDataTree->SetBranchAddress(*tdcNdataName, &tdcHits[iplane][iside][isignal]);
rawDataTree->SetBranchAddress(*tdcPaddleName, &tdcPaddle[iplane][iside][isignal]);
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
if (tdcData[itdcdata] == "TimeRaw") {
tdcTimeRawName = new TString(*tdcBaseName+"TimeRaw");
rawDataTree->SetBranchAddress(*tdcTimeRawName, &hodoTdcTimeRaw[iplane][iside][isignal]);
} // Raw TDC time leaf
} // TDC signal data loop
} // TDC signal
} // Signal loop
} // Side loop
} // Plane loop
// Loop over the events and fill histograms
nentries = rawDataTree->GetEntries();
// nentries = 100000;
cout << "\n******************************************" << endl;
cout << nentries << " Events Will Be Processed" << endl;
cout << "******************************************\n" << endl;
for (ievent = 0; ievent < nentries; ievent++) {
// Acquire the event from the data tree
rawDataTree->GetEntry(ievent);
// Fiducial PID cuts
calEtotCut = (calEtot < calEtotCutVal);
cerNpeSumCut = (cerNpeSum < cerNpeSumCutVal);
if (calEtotCut || cerNpeSumCut) continue;
// Fill trigger apparatus histos
h1_refAdcPulseTimeRaw->Fill(refAdcPulseTimeRaw*adcChanToTime);
h1_refAdcPulseAmp->Fill(refAdcPulseAmp);
h1_refAdcMultiplicity->Fill(refAdcMultiplicity);
h1_refT1TdcTimeRaw->Fill(refT1TdcTimeRaw*tdcChanToTime);
h1_refT2TdcTimeRaw->Fill(refT2TdcTimeRaw*tdcChanToTime);
// Loop over the planes, sides, signals, data arrays, and fill hodoscope histograms
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
for(UInt_t iside = 0; iside < nSides; iside++) {
for(UInt_t isignal = 0; isignal < nSignals; isignal++) {
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Loop over the ADC data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
numAdcHits = adcHits[iplane][iside][isignal];
for (Int_t iadchit = 0; iadchit < numAdcHits; iadchit++) {
// Obtain variables
adcPaddleNum = adcPaddle[iplane][iside][isignal][iadchit];
adcErrorFlag = hodoAdcErrorFlag[iplane][iside][isignal][iadchit];
adcPulseTimeRaw = hodoAdcPulseTimeRaw[iplane][iside][isignal][iadchit]*adcChanToTime;
adcPulseTime = adcPulseTimeRaw - refAdcPulseTimeRaw*adcChanToTime;
adcPulseAmp = hodoAdcPulseAmp[iplane][iside][isignal][iadchit];
// Fill histos
if (adcData[iadcdata] == "ErrorFlag")
h2_adcErrorFlag[iplane][iside]->Fill(adcPaddleNum, adcErrorFlag);
if (adcData[iadcdata] == "PulseTimeRaw") {
h2_adcPulseTimeRaw[iplane][iside]->Fill(adcPaddleNum, adcPulseTimeRaw);
h2_adcPulseTime[iplane][iside]->Fill(adcPaddleNum, adcPulseTime);
}
if (adcData[iadcdata] == "PulseAmp")
h2_adcPulseAmp[iplane][iside]->Fill(adcPaddleNum, adcPulseAmp);
} // ADC hit loop
} // ADC signal data loop
} // ADC signal
// Acquire the hodoscope TDC data objects
if(signalNames[isignal] == "Tdc") {
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
numTdcHits = tdcHits[iplane][iside][isignal];
// Define cuts
// edtmCut = (numTdcHits == nbars[0] || numTdcHits == nbars[2] || numTdcHits == nbars[3]);
// Implement cuts
// if (edtmCut) continue;
for (Int_t itdchit = 0; itdchit < numTdcHits; itdchit++) {
// Obtain variables
tdcPaddleNum = tdcPaddle[iplane][iside][isignal][itdchit];
tdcTimeRaw = hodoTdcTimeRaw[iplane][iside][isignal][itdchit]*tdcChanToTime;
tdcTime = tdcTimeRaw - refT2TdcTimeRaw*tdcChanToTime;
if (tdcData[itdcdata] == "TimeRaw") {
h2_tdcTimeRaw[iplane][iside]->Fill(tdcPaddleNum, tdcTimeRaw);
h2_tdcTime[iplane][iside]->Fill(tdcPaddleNum, tdcTime);
} // RDC raw time signal
} // TDC hit loop
} // TDC signal data loop
} // TDC signal
// Define cuts
adcRefMultiplicityCut = (refAdcMultiplicity != 1.0);
adcRefPulseAmpCut = (refAdcPulseAmp < refAdcPulseAmpCutLow || refAdcPulseAmp > refAdcPulseAmpCutHigh);
adcRefPulseTimeCut = (refAdcPulseTimeRaw*adcChanToTime < refAdcPulseTimeCutLow || refAdcPulseTimeRaw*adcChanToTime > refAdcPulseTimeCutHigh);
// Implement cuts
// if (adcRefMultiplicityCut || adcRefPulseAmpCut || adcRefPulseTimeCut) continue;
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Loop over the signals again
for(UInt_t jsignal = 0; jsignal < nSignals; jsignal++) {
// Acquire othe hodoscope TDC data objects
if(signalNames[jsignal] == "Tdc") {
// Loop over the ADC data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
if (adcData[iadcdata] != "PulseTimeRaw") continue;
// Acquire the number of ADC hits
numAdcHits = adcHits[iplane][iside][isignal];
for (Int_t iadchit = 0; iadchit < numAdcHits; iadchit++) {
// Obtain variables
adcPaddleNum = UInt_t (adcPaddle[iplane][iside][isignal][iadchit]);
adcErrorFlag = hodoAdcErrorFlag[iplane][iside][isignal][iadchit];
adcPulseTimeRaw = hodoAdcPulseTimeRaw[iplane][iside][isignal][iadchit]*adcChanToTime;
adcPulseTime = adcPulseTimeRaw - refAdcPulseTimeRaw*adcChanToTime;
adcPulseAmp = hodoAdcPulseAmp[iplane][iside][isignal][iadchit];
// Define cuts
adcErrorFlagCut = (adcErrorFlag != 0);
adcPulseAmpCut = (adcPulseAmp < hodoPulseAmpCutLow || adcPulseAmp > hodoPulseAmpCutHigh);
// Implement cuts
if (adcErrorFlagCut || adcPulseAmpCut) continue;
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
if (tdcData[itdcdata] != "TimeRaw") continue;
// Acquire the number of TDC hits
numTdcHits = tdcHits[iplane][iside][jsignal];
// Define cuts
//edtmCut = (numTdcHits == nbars[0] || numTdcHits == nbars[2] || numTdcHits == nbars[3]);
// Implement cuts
//if (edtmCut) continue;
for (Int_t itdchit = 0; itdchit < numTdcHits; itdchit++) {
// Obtain variables
tdcPaddleNum = UInt_t (tdcPaddle[iplane][iside][jsignal][itdchit]);
tdcTimeRaw = hodoTdcTimeRaw[iplane][iside][jsignal][itdchit]*tdcChanToTime;
tdcTime = tdcTimeRaw - refT2TdcTimeRaw*tdcChanToTime;
adcTdcTimeDiff = tdcTime - adcPulseTime;
// Define cuts
adcAndTdcHitCut = (adcPaddleNum != tdcPaddleNum);
adcTdcTimeDiffCut = (adcTdcTimeDiff < adcTdcTimeDiffCutLow || adcTdcTimeDiff > adcTdcTimeDiffCutHigh);
// Implement cuts
if (adcAndTdcHitCut || adcTdcTimeDiffCut) continue;
h2_adcPulseAmpCuts[iplane][iside]->Fill(tdcPaddleNum, adcPulseAmp);
h2_adcTdcTimeDiff[iplane][iside]->Fill(tdcPaddleNum, adcTdcTimeDiff);
h2_adcTimeWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, adcPulseTime);
h2_tdcTimeWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, tdcTime);
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, adcTdcTimeDiff);
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->GetXaxis()->CenterTitle();
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->GetYaxis()->CenterTitle();
} // TDC hit loop
} // TDC signal loop
} // ADC hit loop
} // ADC signal loop
} // TDC signal
} // Nested signal loop
} // ADC signal
} // Signal loop
} // Side loop
} // Plane loop
if (ievent % 100000 == 0 && ievent != 0)
cout << ievent << " Events Have Been Processed..." << endl;
} // rawDataTree event loop
cout << "\n***************************************" << endl;
cout << ievent << " Events Were Processed" << endl;
cout << "***************************************\n" << endl;
// Calculate the analysis rate
t = clock() - t;
printf ("The Analysis Took %.1f seconds \n", ((float) t) / CLOCKS_PER_SEC);
printf ("The Analysis Event Rate Was %.3f kHz \n", (ievent + 1) / (((float) t) / CLOCKS_PER_SEC*1000.));
outFile->Write();
//outFile->Close();
return 0;
} // time_walk_calib()
CALIBRATION/shms_hodo_calib/timeWalkCalib.C 0 → 100644
+ 296
− 0
View file @ e131a4d9
// Macro to perform time-walk fits and extract the calibration parameters
// Author: Eric Pooser, pooser@jlab.org
// Declare ROOT files
TFile *histoFile, *outFile;
// Declare constants
static const UInt_t nPlanes = 4;
static const UInt_t nSides = 2;
static const UInt_t nBarsMax = 21;
static const UInt_t nTwFitPars = 3;
static const Double_t tdcThresh = 120.0; // 30 mV in units of FADC channels
static const Double_t twFitRangeLow = 20.0;
static const Double_t twFitRangeHigh = 400.0;
static const Double_t c0twParInit = -35.0;
static const Double_t c1twParInit = 5.0;
static const Double_t c2twParInit = 0.25;
static const Double_t fontSize = 0.05;
static const Double_t yTitleOffset = 0.75;
static const Double_t markerSize = 2.0;
static const Double_t minScale = 0.75;
static const Double_t maxScale = 0.75;
static const UInt_t lineWidth = 4;
static const UInt_t lineStyle = 7;
static const UInt_t nbars[nPlanes] = {13, 13, 14, 21};
static const TString planeNames[nPlanes] = {"1x", "1y", "2x", "2y"};
static const TString sideNames[nSides] = {"pos", "neg"};
static const TString twFitParNames[nTwFitPars] = {"c_{0}", "c_{1}", "c_{2}"};
// Declare directories
TDirectory *dataDir, *planeDir[nPlanes], *sideDir[nPlanes][nSides];
TDirectory *twDir[nPlanes][nSides];
// Declare fits
TF1 *twFit[nPlanes][nSides][nBarsMax], *avgParFit[nPlanes][nSides][nTwFitPars];
// Declare arrays
Double_t paddleIndex[nPlanes][nSides][nBarsMax];
Double_t twFitPar[nPlanes][nSides][nTwFitPars][nBarsMax], twFitParErr[nPlanes][nSides][nTwFitPars][nBarsMax];
Double_t avgPar[nPlanes][nSides][nTwFitPars];
Double_t minPar[nPlanes][nSides][nTwFitPars], maxPar[nPlanes][nSides][nTwFitPars];
// Declare canvases
TCanvas *twFitCan[nPlanes][nSides], *twFitParCan[nTwFitPars];
// Declare histos
// 2D histos
TH2F *h2_adcTdcTimeDiffWalk[nPlanes][nSides][nBarsMax];
// Declare graphs
TGraph *twFitParGraph[nPlanes][nSides][nTwFitPars];
// Declare multi-graphs
TMultiGraph *twFitParMultiGraph[nPlanes][nTwFitPars];
// Declare legends
TLegend *twFitParLegend[nPlanes][nTwFitPars], *avgParLegend[nPlanes][nTwFitPars];
// Declare text objects
TPaveText *twFitParText[nPlanes][nSides][nBarsMax];
// Functions to extract calibration parameters
//=:=:=:=:=:=:
//=: Level 1
//=:=:=:=:=:=:
// Add color to fit lines
void addColorToFitLine(UInt_t style, UInt_t width, UInt_t color, TF1 *fit1D) {
fit1D->SetLineStyle(lineStyle);
fit1D->SetLineWidth(lineWidth);
fit1D->SetLineColor(color);
return;
} // addColorToFitLine()
// Add color to fit graphs
void addColorToGraph(UInt_t style, UInt_t size, UInt_t color, TGraph *graph) {
graph->SetMarkerStyle(style);
graph->SetMarkerSize(size);
graph->SetMarkerColor(color);
return;
} // addColorToFitGraph()
// Add centered titles to multigraphs
void modifyMultiGraph(TMultiGraph *mg, TString xtitle, TString ytitle) {
mg->GetYaxis()->SetTitle(ytitle);
mg->GetYaxis()->CenterTitle();
mg->GetXaxis()->SetTitle(xtitle);
mg->GetXaxis()->CenterTitle();
return;
} // modifyMultiGraph()
// Create canvases
TCanvas *makeCan(UInt_t numColumns, UInt_t numRows, UInt_t winWidth, UInt_t winHeight, TCanvas *can, TString name, TString title) {
can = new TCanvas(name, title, winWidth, winHeight);
can->Divide(numColumns, numRows);
return can;
} // makeCan()
// Define the time-walk fit function
Double_t twFitFunc(Double_t *a, Double_t *c) {
Double_t twFitVal = c[0] + c[1]/(TMath::Power((a[0]/tdcThresh), c[2]));
return twFitVal;
} // twFitFunc()
// Locate min or max value from input array
Double_t calcMinOrMax(Double_t *array, UInt_t iplane, TString minOrmax) {
auto result = minmax_element(array, array+nbars[iplane]);
if (minOrmax == "min") return *result.first;
else if (minOrmax == "max") return *result.second;
else return 0.0;
} // calcMinOrMax()
//=:=:=:=:=:=:
//=: Level 2
//=:=:=:=:=:=:
// Perform the timw-walk fits
void doTwFits(UInt_t iplane, UInt_t iside, UInt_t ipaddle) {
// Draw fits on canvas
twFitCan[iplane][iside]->cd(ipaddle+1);
gPad->SetLogz();
h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->Draw("COLZ");
gPad->Modified(); gPad->Update();
// Declare and initialize the fits
twFit[iplane][iside][ipaddle] = new TF1("twFit", twFitFunc, twFitRangeLow, twFitRangeHigh, nTwFitPars);
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFit[iplane][iside][ipaddle]->SetParName(ipar, twFitParNames[ipar]);
twFit[iplane][iside][ipaddle]->SetParameters(c0twParInit, c1twParInit, c2twParInit);
// Perform the fits and scream if it failed
Int_t twFitStatus = h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->Fit("twFit", "REQ");
if (twFitStatus != 0)
cout << "ERROR: Time Walk Fit Failed!!! " << "Status = " << twFitStatus << " For Plane: " << planeNames[iplane] << " Side: " << sideNames[iside] << " Paddle: " << ipaddle+1 << endl;
// Create text box to display fir parameters
twFitParText[iplane][iside][ipaddle] = new TPaveText(0.4, 0.6, 0.895, 0.895, "NBNDC");
twFitParText[iplane][iside][ipaddle]->AddText(Form("Entries = %.0f", h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]->GetEntries()));
// Obtain the fit parameters and associated errors
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
twFitPar[iplane][iside][ipar][ipaddle] = twFit[iplane][iside][ipaddle]->GetParameter(ipar);
twFitParErr[iplane][iside][ipar][ipaddle] = twFit[iplane][iside][ipaddle]->GetParError(ipar);
twFitParText[iplane][iside][ipaddle]->AddText(Form(twFitParNames[ipar]+" = %.2f #pm %.2f", twFitPar[iplane][iside][ipar][ipaddle], twFitParErr[iplane][iside][ipar][ipaddle]));
} // Parameter loop
// Draw the fit parameter text
twFitParText[iplane][iside][ipaddle]->SetFillColor(kWhite);
twFitParText[iplane][iside][ipaddle]->SetTextAlign(12);
twFitParText[iplane][iside][ipaddle]->SetFillStyle(3003);
twFitParText[iplane][iside][ipaddle]->Draw();
gPad->Modified(); gPad->Update();
return;
} // doTwFits()
// Calculate the averege of the time-walk fit parameters
void calcParAvg(UInt_t iplane, UInt_t iside) {
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
// Calculate the weighted average while ignoring fit errors provided by Minuit
avgParFit[iplane][iside][ipar] = new TF1("avgParFit", "pol0", 1, nbars[iplane]);
avgParFit[iplane][iside][ipar]->SetParName(0, "#bar{"+twFitParNames[ipar]+"}");
// Add color to fit lines
if (iside == 0) addColorToFitLine(lineStyle, lineWidth, kRed, avgParFit[iplane][iside][ipar]);
if (iside == 1) addColorToFitLine(lineStyle, lineWidth, kBlue, avgParFit[iplane][iside][ipar]);
// Initialize the parameters
if (ipar == 0) avgParFit[iplane][iside][ipar]->SetParameter(0, c0twParInit);
if (ipar == 1) avgParFit[iplane][iside][ipar]->SetParameter(0, c1twParInit);
if (ipar == 2) avgParFit[iplane][iside][ipar]->SetParameter(0, c2twParInit);
// Perform the least squares fit
Int_t avgParFitStatus = twFitParGraph[iplane][iside][ipar]->Fit("avgParFit", "REQ");
if (avgParFitStatus != 0)
cout << "ERROR: Parameter Fit Failed!!! " << "Status = " << avgParFitStatus << " For Plane: " << planeNames[iplane] << " Side: " << sideNames[iside] << endl;
// Store the fit parameters
avgPar[iplane][iside][ipar] = avgParFit[iplane][iside][ipar]->GetParameter(0);
// Add graphs to multi graph
twFitParMultiGraph[iplane][ipar]->Add(twFitParGraph[iplane][iside][ipar]);
// Calculate the min and max value of each parameter
minPar[iplane][iside][ipar] = calcMinOrMax(twFitPar[iplane][iside][ipar], iplane, "min");
maxPar[iplane][iside][ipar] = calcMinOrMax(twFitPar[iplane][iside][ipar], iplane, "max");
} // Parameter loop
return;
} // calcParAvg()
// Draw the time-walk fit parameter multi-graph
void drawParams(UInt_t iplane) {
// Populate the parameter canvases
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
twFitParCan[ipar]->cd(iplane+1);
twFitParMultiGraph[iplane][ipar]->Draw("AP");
gPad->Modified(); gPad->Update();
// Add aesthetics
modifyMultiGraph(twFitParMultiGraph[iplane][ipar], "Paddle Number", "Fit Parameter "+twFitParNames[ipar]);
twFitParMultiGraph[iplane][ipar]->GetXaxis()->SetNdivisions(21);
twFitParMultiGraph[iplane][ipar]->Draw("AP");
// Draw legends
twFitParLegend[iplane][ipar] = new TLegend(0.1, 0.7, 0.3, 0.9);
avgParLegend[iplane][ipar] = new TLegend(0.6, 0.7, 0.9, 0.9);
for(UInt_t iside = 0; iside < nSides; iside++) {
if (iside == 0) twFitParLegend[iplane][ipar]->AddEntry(twFitParGraph[iplane][iside][ipar], "Positive Side", "p");
if (iside == 1) twFitParLegend[iplane][ipar]->AddEntry(twFitParGraph[iplane][iside][ipar], "Negative Side", "p");
avgParLegend[iplane][ipar]->AddEntry(avgParFit[iplane][iside][ipar], "#bar{"+twFitParNames[ipar]+"}"+Form(" = %.2f", avgPar[iplane][iside][ipar]),"l");
// Define the time-walk fit parameter multigraph y-axis range
if (iside == 0) {
// Set y-axis range min
if (minPar[iplane][iside][ipar] < 0.0 || minPar[iplane][iside+1][ipar] < 0.0) {
if (minPar[iplane][iside][ipar] < minPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside][ipar] + minScale*minPar[iplane][iside][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside+1][ipar] + minScale*minPar[iplane][iside+1][ipar]);
}
if (minPar[iplane][iside][ipar] > 0.0 || minPar[iplane][iside+1][ipar] > 0.0) {
if (minPar[iplane][iside][ipar] < minPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside][ipar] - minScale*minPar[iplane][iside][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMinimum(minPar[iplane][iside+1][ipar] - minScale*minPar[iplane][iside+1][ipar]);
}
// Set y-axis range max
if (maxPar[iplane][iside][ipar] < 0.0 || maxPar[iplane][iside+1][ipar] < 0.0) {
if (maxPar[iplane][iside][ipar] < maxPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside+1][ipar] - maxScale*maxPar[iplane][iside+1][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside][ipar] - maxScale*maxPar[iplane][iside][ipar]);
}
if (maxPar[iplane][iside][ipar] > 0.0 || maxPar[iplane][iside+1][ipar] > 0.0) {
if (maxPar[iplane][iside][ipar] < maxPar[iplane][iside+1][ipar])
twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside+1][ipar] + maxScale*maxPar[iplane][iside+1][ipar]);
else twFitParMultiGraph[iplane][ipar]->SetMaximum(maxPar[iplane][iside][ipar] + maxScale*maxPar[iplane][iside][ipar]);
}
} // Side index = 0 condition
} // Side loop
twFitParLegend[iplane][ipar]->Draw();
avgParLegend[iplane][ipar]->Draw();
} // Parameter loop
return;
} // drawParams
//=:=:=:=:=
//=: Main
//=:=:=:=:=
void timeWalkCalib() {
// ROOT settings
gStyle->SetTitleFontSize(fontSize);
gStyle->SetLabelSize(fontSize, "XY");
gStyle->SetTitleSize(fontSize, "XY");
gStyle->SetTitleOffset(yTitleOffset, "Y");
gStyle->SetStatFormat(".2f");
gStyle->SetOptFit(0);
gStyle->SetOptStat(0);
// Read the ROOT file containing the time-walk histos
histoFile = new TFile("timeWalkHistos.root", "READ");
// Obtain the top level directory
dataDir = dynamic_cast <TDirectory*> (histoFile->FindObjectAny("hodoUncalib"));
// Create the parameter canvases
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFitParCan[ipar] = makeCan(2, 2, 1600, 800, twFitParCan[ipar], twFitParNames[ipar]+"FitParCan", "Parameter "+twFitParNames[ipar]+" Canvas");
// Loop over the planes
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
//for(UInt_t iplane = 0; iplane < 1; iplane++) {
// Obtain the plane directory
planeDir[iplane] = dynamic_cast <TDirectory*> (dataDir->FindObjectAny(planeNames[iplane]));
// Create multigraphs
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++)
twFitParMultiGraph[iplane][ipar] = new TMultiGraph(planeNames[iplane]+"_"+twFitParNames[ipar]+"_Multigraph", "Plane "+planeNames[iplane]+" Parameter "+twFitParNames[ipar]);
// Loop over the sides
for(UInt_t iside = 0; iside < nSides; iside++) {
// Obtain the side and time walk directories
sideDir[iplane][iside] = dynamic_cast <TDirectory*> (planeDir[iplane]->FindObjectAny(sideNames[iside]));
twDir[iplane][iside] = dynamic_cast <TDirectory*> (sideDir[iplane][iside]->FindObjectAny("adcTdcTimeDiffWalk"));
// Create the time-walk histo and fit canvases
if (planeNames[iplane] != "2y") twFitCan[iplane][iside] = makeCan(5, 3, 1600, 800, twFitCan[iplane][iside], planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan", planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan");
if (planeNames[iplane] == "2y") twFitCan[iplane][iside] = makeCan(6, 4, 1600, 800, twFitCan[iplane][iside], planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan", planeNames[iplane]+"_"+sideNames[iside]+"_twFitCan");
// Loop over the paddles
for(UInt_t ipaddle = 0; ipaddle < nbars[iplane]; ipaddle++) {
// Populate the paddle index arrays
paddleIndex[iplane][iside][ipaddle] = Double_t (ipaddle + 1);
// Obtain the time-walk histos
h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle] = dynamic_cast <TH2F*> (twDir[iplane][iside]->FindObjectAny(Form("h2_adcTdcTimeDiffWalk_paddle_%d", ipaddle+1)));
// Perform the time-walk fits while ignoring the S2Y plane (for now)
// if (planeNames[iplane] == "2y") continue;
doTwFits(iplane, iside, ipaddle);
} // Paddle loop
// Produce the time-walk fit parameter graphs
for (UInt_t ipar = 0; ipar < nTwFitPars; ipar++) {
// Populate graphs and multi-graphs
twFitParGraph[iplane][iside][ipar] = new TGraph(nbars[iplane], paddleIndex[iplane][iside], twFitPar[iplane][iside][ipar]);
if (iside == 0) addColorToGraph(22, markerSize, kRed, twFitParGraph[iplane][iside][ipar]);
if (iside == 1) addColorToGraph(23, markerSize, kBlue, twFitParGraph[iplane][iside][ipar]);
// // Add graphs to multi graph
// twFitParMultiGraph[iplane][ipar]->Add(twFitParGraph[iplane][iside][ipar]);
} // Parameter loop
// Calculate the average of the time-walk fit parameters
calcParAvg(iplane, iside);
} // Side loop
// Draw the time-walk parameter graphs
drawParams(iplane);
} // Plane loop
return;
} // timeWalkFits()
CALIBRATION/shms_hodo_calib/timeWalkHistos.C 0 → 100644
+ 429
− 0
View file @ e131a4d9
// Macro to generate time-walk histograms for the hodoscopes
// Author: Eric Pooser, pooser@jlab.org
#include <time.h>
// Declare replay data file and output file
TFile *replayFile, *outFile;
// Declare data tree
TTree *rawDataTree;
// Declare constants
static const UInt_t nPlanes = 4;
static const UInt_t nSides = 2;
static const UInt_t nAdcSignals = 3;
static const UInt_t nTdcSignals = 1;
static const UInt_t nSignals = nAdcSignals + nTdcSignals;
static const UInt_t maxTdcHits = 128;
static const UInt_t maxAdcHits = 4;
static const UInt_t nMaxBars = 21;
static const TString spec = "P";
static const TString detec = "hod";
static const TString trig = "T";
static const TString shms = "shms";
static const UInt_t nbars[nPlanes] = {13, 13, 14, 21};
static const TString planeNames[nPlanes] = {"1x", "1y", "2x", "2y"};
static const TString sideNames[nSides] = {"pos", "neg"};
static const TString signalNames[nSignals] = {"Adc", "Tdc"};
static const TString adcData[nAdcSignals] = {"ErrorFlag", "PulseTimeRaw", "PulseAmp"};
static const TString tdcData[nTdcSignals] = {"TimeRaw"};
static const Double_t tdcChanToTime = 0.100; // Units of ns
static const Double_t adcChanToTime = 0.0625; // Units of ns
static const Double_t adcDynamicRange = 1000.0; // Units of mV
static const Double_t nAdcChan = 4096.0; // Units of ADC channels
static const Double_t adcChanTomV = adcDynamicRange/nAdcChan; // Units of mV/ADC Chan
static const Double_t hodoPulseAmpCutLow = 10.0; // Units of mV
static const Double_t hodoPulseAmpCutHigh = 1000.0; // Units of mV
static const Double_t refAdcPulseAmpCutLow = 50.0; // Units of mV
static const Double_t refAdcPulseAmpCutHigh = 60.0; // Units of mV
static const Double_t refAdcPulseTimeCutLow = 210.0; // Units of ns
static const Double_t refAdcPulseTimeCutHigh = 225.0; // Units of ns
static const Double_t adcTdcTimeDiffCutLow = -100.0; // Units of ns
static const Double_t adcTdcTimeDiffCutHigh = 100.0; // Units of ns
static const Double_t calEtotCutVal = 0.100; // Units of GeV
static const Double_t cerNpeSumCutVal = 1.5; // Units of NPE
// static const Double_t adcTdcTimeDiffCutLow = -6000.0; // Units of ns
// static const Double_t adcTdcTimeDiffCutHigh = 1000.0; // Units of ns
static const Double_t fontSize = 0.04;
// Declare variables to obtain from data tree
// Number of ADC & TDC hits
Int_t adcHits[nPlanes][nSides][nSignals], tdcHits[nPlanes][nSides][nSignals];
// Paddle number which fired
Double_t adcPaddle[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t tdcPaddle[nPlanes][nSides][nSignals][nMaxBars*maxTdcHits];
// FADC data
Double_t hodoAdcPulseTimeRaw[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t hodoAdcPulseAmp[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
Double_t hodoAdcErrorFlag[nPlanes][nSides][nSignals][nMaxBars*maxAdcHits];
// TDC Data
Double_t hodoTdcTimeRaw[nPlanes][nSides][nSignals][nMaxBars*maxTdcHits];
// Trigger apparatus data
Double_t refAdcPulseTimeRaw, refAdcPulseAmp, refAdcMultiplicity;
Double_t refT1TdcTimeRaw, refT2TdcTimeRaw, refT3TdcTimeRaw;
// Declare variables to steer data to arrays
TString *adcBaseName, *adcNdataName, *adcPaddleName;
TString *adcErrorFlagName, *adcPulseTimeRawName, *adcPulseAmpName;
TString *tdcBaseName, *tdcNdataName, *tdcPaddleName;
TString *tdcTimeRawName;
// Declare directories for output ROOT file
TDirectory *trigRawDir, *hodoRawDir, *hodoUncalibDir;
TDirectory *planeRawDir[nPlanes], *sideRawDir[nPlanes][nSides], *planeUncalibDir[nPlanes], *sideUncalibDir[nPlanes][nSides];
TDirectory *adcTimeWalkDir[nPlanes][nSides], *tdcTimeWalkDir[nPlanes][nSides], *adcTdcTimeDiffWalkDir[nPlanes][nSides];
// Declare histos
// 1D histos
TH1F *h1_refAdcPulseTimeRaw, *h1_refAdcPulseAmp, *h1_refAdcMultiplicity;
TH1F *h1_refT1TdcTimeRaw, *h1_refT2TdcTimeRaw, *h1_refT3TdcTimeRaw;
// 2D histos
TH2F *h2_adcErrorFlag[nPlanes][nSides];
TH2F *h2_adcPulseTimeRaw[nPlanes][nSides];
TH2F *h2_adcPulseAmp[nPlanes][nSides];
TH2F *h2_tdcTimeRaw[nPlanes][nSides];
TH2F *h2_adcPulseTime[nPlanes][nSides];
TH2F *h2_tdcTime[nPlanes][nSides];
TH2F *h2_adcTdcTimeDiff[nPlanes][nSides];
TH2F *h2_adcPulseAmpCuts[nPlanes][nSides];
TH2F *h2_adcTimeWalk[nPlanes][nSides][nMaxBars];
TH2F *h2_tdcTimeWalk[nPlanes][nSides][nMaxBars];
TH2F *h2_adcTdcTimeDiffWalk[nPlanes][nSides][nMaxBars];
// Declare event data objects
UInt_t nentries, ievent, adcPaddleNum, tdcPaddleNum;
Int_t numAdcHits, numTdcHits;
Double_t adcErrorFlag, adcPulseTimeRaw, adcPulseTime, adcPulseAmp;
Double_t tdcTimeRaw, tdcTime, adcTdcTimeDiff;
Double_t calEtot, cerNpeSum;
Bool_t adcRefMultiplicityCut, adcRefPulseAmpCut, adcRefPulseTimeCut;
Bool_t edtmCut, adcErrorFlagCut, adcAndTdcHitCut;
Bool_t adcPulseAmpCut, adcTdcTimeDiffCut;
Bool_t calEtotCut, cerNpeSumCut;
void generatePlots(UInt_t iplane, UInt_t iside, UInt_t ipaddle) {
// Create trigger aparatus directory
trigRawDir = dynamic_cast <TDirectory*> (outFile->Get("trigAppRaw"));
if (!trigRawDir) {trigRawDir = outFile->mkdir("trigAppRaw"); trigRawDir->cd();}
else outFile->cd("trigAppRaw");
// FADC reference
if (!h1_refAdcPulseTimeRaw) h1_refAdcPulseTimeRaw = new TH1F("h1_refAdcPulseTimeRaw", "ROC2 Raw FADC Reference Pulse Time; Raw FADC Pulse Time (ns); Number of Entries / 100 ps", 4000, 0, 400);
if (!h1_refAdcPulseAmp) h1_refAdcPulseAmp = new TH1F("h1_refAdcPulseAmp", "ROC2 FADC Reference Pulse Amplitude; FADC Pulse Amplitude (mV); Number of Entries / 1 mV", 1000, 0, 1000);
if (!h1_refAdcMultiplicity) h1_refAdcMultiplicity = new TH1F("h1_refAdcMultiplicity", "ROC2 FADC Reference Multiplicity; Raw FADC Multiplicity; Number of Entries", 5, -0.5, 4.5);
// TDC reference
if (!h1_refT1TdcTimeRaw) h1_refT1TdcTimeRaw = new TH1F("h1_refT1TdcTimeRaw", "ROC2 T1 Raw TDC Reference TDC Time (Slot 20, Channel 15); Raw TDC Time (ns); Number of Entries / 100 ps", 6000, 0, 600);
if (!h1_refT2TdcTimeRaw) h1_refT2TdcTimeRaw = new TH1F("h1_refT2TdcTimeRaw", "ROC2 T2 Raw TDC Reference TDC Time (Slot 19, Channel 31); Raw TDC Time (ns); Number of Entries / 100 ps", 6000, 0, 600);
if (!h1_refT3TdcTimeRaw) h1_refT3TdcTimeRaw = new TH1F("h1_refT3TdcTimeRaw", "ROC2 T3 Raw TDC Reference TDC Time (Slot 19, Channel 38); Raw TDC Time (ns); Number of Entries / 100 ps", 6000, 0, 600);
// Create plane directory for raw hodoscope quantities
hodoRawDir = dynamic_cast <TDirectory*> (outFile->Get("hodoRaw"));
if (!hodoRawDir) {hodoRawDir = outFile->mkdir("hodoRaw"); hodoRawDir->cd();}
else outFile->cd("hodoRaw");
planeRawDir[iplane] = dynamic_cast <TDirectory*> (hodoRawDir->Get(planeNames[iplane]));
if (!planeRawDir[iplane]) {planeRawDir[iplane] = hodoRawDir->mkdir(planeNames[iplane]); planeRawDir[iplane]->cd();}
else hodoRawDir->cd(planeNames[iplane]);
// Create side directory
sideRawDir[iplane][iside] = dynamic_cast <TDirectory*> (planeRawDir[iplane]->Get(sideNames[iside]));
if (!sideRawDir[iplane][iside]) {sideRawDir[iplane][iside] = planeRawDir[iplane]->mkdir(sideNames[iside]); sideRawDir[iplane][iside]->cd();}
else outFile->cd("hodoRaw/"+planeNames[iplane]+"/"+sideNames[iside]);
// Book histos
if (!h2_adcErrorFlag[iplane][iside]) h2_adcErrorFlag[iplane][iside] = new TH2F("h2_adcErrorFlag", "FADC Error Flag Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Error Flag", nbars[iplane], 0.5, nbars[iplane]+0.5, 2, -0.5, 1.5);
if (!h2_adcPulseTimeRaw[iplane][iside]) h2_adcPulseTimeRaw[iplane][iside] = new TH2F("h2_adcPulseTimeRaw", "Raw FADC Pulse Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; Raw FADC Pulse Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 4000, 0, 400);
if (!h2_adcPulseAmp[iplane][iside]) h2_adcPulseAmp[iplane][iside] = new TH2F("h2_adcPulseAmp", "FADC Pulse Amplitude Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Amplitude (mV) / 1 mV", nbars[iplane], 0.5, nbars[iplane]+0.5, 1000, 0, 1000);
if (!h2_tdcTimeRaw[iplane][iside]) h2_tdcTimeRaw[iplane][iside] = new TH2F("h2_tdcTimeRaw", "Raw TDC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; Raw TDC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 4000, 0, 400);
// Create plane directory for uncalibrated hodoscope quantities
hodoUncalibDir = dynamic_cast <TDirectory*> (outFile->Get("hodoUncalib"));
if (!hodoUncalibDir) {hodoUncalibDir = outFile->mkdir("hodoUncalib"); hodoUncalibDir->cd();}
else outFile->cd("hodoUncalib");
planeUncalibDir[iplane] = dynamic_cast <TDirectory*> (hodoUncalibDir->Get(planeNames[iplane]));
if (!planeUncalibDir[iplane]) {planeUncalibDir[iplane] = hodoUncalibDir->mkdir(planeNames[iplane]); planeUncalibDir[iplane]->cd();}
else hodoUncalibDir->cd(planeNames[iplane]);
// Create side directory
sideUncalibDir[iplane][iside] = dynamic_cast <TDirectory*> (planeUncalibDir[iplane]->Get(sideNames[iside]));
if (!sideUncalibDir[iplane][iside]) {sideUncalibDir[iplane][iside] = planeUncalibDir[iplane]->mkdir(sideNames[iside]); sideUncalibDir[iplane][iside]->cd();}
else outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]);
// Book histos
if (!h2_adcPulseTime[iplane][iside]) h2_adcPulseTime[iplane][iside] = new TH2F("h2_adcPulseTime", "FADC Pulse Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 8000, -400, 400);
if (!h2_tdcTime[iplane][iside]) h2_tdcTime[iplane][iside] = new TH2F("h2_tdcTime", "TDC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; TDC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 8000, -400, 400);
if (!h2_adcTdcTimeDiff[iplane][iside]) h2_adcTdcTimeDiff[iplane][iside] = new TH2F("h2_adcTdcTimeDiff", "TDC-ADC Time Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; TDC - ADC Time (ns) / 100 ps", nbars[iplane], 0.5, nbars[iplane]+0.5, 16000, -800, 800);
if (!h2_adcPulseAmpCuts[iplane][iside]) h2_adcPulseAmpCuts[iplane][iside] = new TH2F("h2_adcPulseAmpCuts", "FADC Pulse Amplitude Cuts Plane "+planeNames[iplane]+" Side "+sideNames[iside]+"; PMT Number; FADC Pulse Amplitude (mV) / 1 mV", nbars[iplane], 0.5, nbars[iplane]+0.5, 1000, 0, 1000);
// Create ADC time-walk directory
adcTimeWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("adcTimeWalk"));
if (!adcTimeWalkDir[iplane][iside]) {adcTimeWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("adcTimeWalk"); adcTimeWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/adcTimeWalk"));
// Book histos
if (!h2_adcTimeWalk[iplane][iside][ipaddle]) h2_adcTimeWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_adcTimeWalk_paddle_%d", ipaddle+1), "ADC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; ADC Time (ns) / 100 ps", 1000, 0, 1000, 4000, -200, 200);
// Create TDC time-walk directory
tdcTimeWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("tdcTimeWalk"));
if (!tdcTimeWalkDir[iplane][iside]) {tdcTimeWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("tdcTimeWalk"); tdcTimeWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/tdcTimeWalk"));
// Book histos
if (!h2_tdcTimeWalk[iplane][iside][ipaddle]) h2_tdcTimeWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_tdcTimeWalk_paddle_%d", ipaddle+1), "TDC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; TDC Time (ns) / 100 ps", 1000, 0, 1000, 4000, -300, 100);
// Create TDC-ADC time-walk directory
adcTdcTimeDiffWalkDir[iplane][iside] = dynamic_cast <TDirectory*> (sideUncalibDir[iplane][iside]->Get("adcTdcTimeDiffWalk"));
if (!adcTdcTimeDiffWalkDir[iplane][iside]) {adcTdcTimeDiffWalkDir[iplane][iside] = sideUncalibDir[iplane][iside]->mkdir("adcTdcTimeDiffWalk"); adcTdcTimeDiffWalkDir[iplane][iside]->cd();}
else (outFile->cd("hodoUncalib/"+planeNames[iplane]+"/"+sideNames[iside]+"/adcTdcTimeDiffWalk"));
// Book histos
if (!h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle]) h2_adcTdcTimeDiffWalk[iplane][iside][ipaddle] = new TH2F(Form("h2_adcTdcTimeDiffWalk_paddle_%d", ipaddle+1), "TDC-ADC Time vs. Pulse Amp Plane "+planeNames[iplane]+" Side "+sideNames[iside]+Form(" Paddle %d", ipaddle+1)+"; Pulse Amplitude (mV) / 1 mV; TDC-ADC Time (ns) / 100 ps", 1000, 0, 1000, 500, -20, 30);
} // generatePlots()
void timeWalkHistos(Int_t runNum) {
// Global ROOT settings
gStyle->SetOptFit();
gStyle->SetStatFormat(".2f");
gStyle->SetTitleFontSize(fontSize);
gStyle->SetLabelSize(fontSize, "XY");
gStyle->SetTitleSize(fontSize, "XY");
// Initialize the analysis clock
clock_t t;
t = clock();
// Obtain the replay data file and create new output ROOT file
replayFile = new TFile(Form("ROOTfiles/shms_replay_production_all_%d_-1.root", runNum), "READ");
// replayFile = new TFile(Form("ROOTfiles/shms_coin_replay_production_%d_-1.root", runNum), "READ");
outFile = new TFile("timeWalkHistos_temp.root", "RECREATE");
// Obtain the tree
rawDataTree = dynamic_cast <TTree*> (replayFile->Get("T"));
// Acquire the trigger apparatus data
rawDataTree->SetBranchAddress("T.shms.pFADC_TREF_ROC2_adcPulseTimeRaw", &refAdcPulseTimeRaw);
rawDataTree->SetBranchAddress("T.shms.pFADC_TREF_ROC2_adcPulseAmp", &refAdcPulseAmp);
rawDataTree->SetBranchAddress("T.shms.pFADC_TREF_ROC2_adcMultiplicity", &refAdcMultiplicity);
rawDataTree->SetBranchAddress("T.shms.pT1_tdcTimeRaw", &refT1TdcTimeRaw);
rawDataTree->SetBranchAddress("T.shms.pT2_tdcTimeRaw", &refT2TdcTimeRaw);
rawDataTree->SetBranchAddress("T.shms.pT3_tdcTimeRaw", &refT3TdcTimeRaw);
rawDataTree->SetBranchAddress("P.cal.etot", &calEtot);
rawDataTree->SetBranchAddress("P.ngcer.npeSum", &cerNpeSum);
// Loop over the planes, sides, signals, leafs, and fill data arrays
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
for(UInt_t iside = 0; iside < nSides; iside++) {
// Generate directory structure and histograms
for(UInt_t ipaddle = 0; ipaddle < nbars[iplane]; ipaddle++)
generatePlots(iplane, iside, ipaddle);
for(UInt_t isignal = 0; isignal < nSignals; isignal++) {
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Define the leaf strings of interest
adcBaseName = new TString(spec+"."+detec+"."+planeNames[iplane]+"."+
sideNames[iside]+signalNames[isignal]);
adcNdataName = new TString("Ndata."+*adcBaseName+"Counter");
adcPaddleName = new TString(*adcBaseName+"Counter");
// Fill the leafs
rawDataTree->SetBranchAddress(*adcNdataName, &adcHits[iplane][iside][isignal]);
rawDataTree->SetBranchAddress(*adcPaddleName, &adcPaddle[iplane][iside][isignal]);
// Loop over the adc data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
if (adcData[iadcdata] == "ErrorFlag") {
adcErrorFlagName = new TString(*adcBaseName+"ErrorFlag");
rawDataTree->SetBranchAddress(*adcErrorFlagName, &hodoAdcErrorFlag[iplane][iside][isignal]);
} // Error flag leaf
if (adcData[iadcdata] == "PulseTimeRaw") {
adcPulseTimeRawName = new TString(*adcBaseName+"PulseTimeRaw");
rawDataTree->SetBranchAddress(*adcPulseTimeRawName, &hodoAdcPulseTimeRaw[iplane][iside][isignal]);
} // Raw pulse time leaf
if (adcData[iadcdata] == "PulseAmp") {
adcPulseAmpName = new TString(*adcBaseName+"PulseAmp");
rawDataTree->SetBranchAddress(*adcPulseAmpName, &hodoAdcPulseAmp[iplane][iside][isignal]);
} // Pedestal subtracted pulse amplitude leaf
} // ADC signal data loop
} // ADC signal
// Acquire the hodoscope TDC data objects
if(signalNames[isignal] == "Tdc") {
// Define the leaf strings of interest
tdcBaseName = new TString(spec+"."+detec+"."+planeNames[iplane]+"."+
sideNames[iside]+signalNames[isignal]);
tdcNdataName = new TString("Ndata."+*tdcBaseName+"Counter");
tdcPaddleName = new TString(*tdcBaseName+"Counter");
// Fill the leafs
rawDataTree->SetBranchAddress(*tdcNdataName, &tdcHits[iplane][iside][isignal]);
rawDataTree->SetBranchAddress(*tdcPaddleName, &tdcPaddle[iplane][iside][isignal]);
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
if (tdcData[itdcdata] == "TimeRaw") {
tdcTimeRawName = new TString(*tdcBaseName+"TimeRaw");
rawDataTree->SetBranchAddress(*tdcTimeRawName, &hodoTdcTimeRaw[iplane][iside][isignal]);
} // Raw TDC time leaf
} // TDC signal data loop
} // TDC signal
} // Signal loop
} // Side loop
} // Plane loop
// Loop over the events and fill histograms
// nentries = rawDataTree->GetEntries();
nentries = 100000;
cout << "\n******************************************" << endl;
cout << nentries << " Events Will Be Processed" << endl;
cout << "******************************************\n" << endl;
for (ievent = 0; ievent < nentries; ievent++) {
// Acquire the event from the data tree
rawDataTree->GetEntry(ievent);
// Fiducial PID cuts
calEtotCut = (calEtot < calEtotCutVal);
cerNpeSumCut = (cerNpeSum < cerNpeSumCutVal);
if (calEtotCut || cerNpeSumCut) continue;
// Fill trigger apparatus histos
h1_refAdcPulseTimeRaw->Fill(refAdcPulseTimeRaw*adcChanToTime);
h1_refAdcPulseAmp->Fill(refAdcPulseAmp);
h1_refAdcMultiplicity->Fill(refAdcMultiplicity);
h1_refT1TdcTimeRaw->Fill(refT1TdcTimeRaw*tdcChanToTime);
h1_refT2TdcTimeRaw->Fill(refT2TdcTimeRaw*tdcChanToTime);
h1_refT3TdcTimeRaw->Fill(refT3TdcTimeRaw*tdcChanToTime);
// Loop over the planes, sides, signals, data arrays, and fill hodoscope histograms
for(UInt_t iplane = 0; iplane < nPlanes; iplane++) {
for(UInt_t iside = 0; iside < nSides; iside++) {
for(UInt_t isignal = 0; isignal < nSignals; isignal++) {
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Loop over the ADC data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
numAdcHits = adcHits[iplane][iside][isignal];
for (Int_t iadchit = 0; iadchit < numAdcHits; iadchit++) {
// Obtain variables
adcPaddleNum = adcPaddle[iplane][iside][isignal][iadchit];
adcErrorFlag = hodoAdcErrorFlag[iplane][iside][isignal][iadchit];
adcPulseTimeRaw = hodoAdcPulseTimeRaw[iplane][iside][isignal][iadchit]*adcChanToTime;
adcPulseTime = adcPulseTimeRaw - refAdcPulseTimeRaw*adcChanToTime;
adcPulseAmp = hodoAdcPulseAmp[iplane][iside][isignal][iadchit];
// Fill histos
if (adcData[iadcdata] == "ErrorFlag")
h2_adcErrorFlag[iplane][iside]->Fill(adcPaddleNum, adcErrorFlag);
if (adcData[iadcdata] == "PulseTimeRaw") {
h2_adcPulseTimeRaw[iplane][iside]->Fill(adcPaddleNum, adcPulseTimeRaw);
h2_adcPulseTime[iplane][iside]->Fill(adcPaddleNum, adcPulseTime);
}
if (adcData[iadcdata] == "PulseAmp")
h2_adcPulseAmp[iplane][iside]->Fill(adcPaddleNum, adcPulseAmp);
} // ADC hit loop
} // ADC signal data loop
} // ADC signal
// Acquire the hodoscope TDC data objects
if(signalNames[isignal] == "Tdc") {
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
numTdcHits = tdcHits[iplane][iside][isignal];
// Define cuts
//edtmCut = (numTdcHits == nbars[0] || numTdcHits == nbars[2] || numTdcHits == nbars[3]);
// Implement cuts
//if (edtmCut) continue;
for (Int_t itdchit = 0; itdchit < numTdcHits; itdchit++) {
// Obtain variables
tdcPaddleNum = tdcPaddle[iplane][iside][isignal][itdchit];
tdcTimeRaw = hodoTdcTimeRaw[iplane][iside][isignal][itdchit]*tdcChanToTime;
if (iplane == 3 && iside == 0) tdcTime = tdcTimeRaw - refT2TdcTimeRaw*tdcChanToTime;
else tdcTime = tdcTimeRaw - refT1TdcTimeRaw*tdcChanToTime;
if (tdcData[itdcdata] == "TimeRaw") {
h2_tdcTimeRaw[iplane][iside]->Fill(tdcPaddleNum, tdcTimeRaw);
h2_tdcTime[iplane][iside]->Fill(tdcPaddleNum, tdcTime);
} // RDC raw time signal
} // TDC hit loop
} // TDC signal data loop
} // TDC signal
// Define cuts
adcRefMultiplicityCut = (refAdcMultiplicity != 1.0);
adcRefPulseAmpCut = (refAdcPulseAmp < refAdcPulseAmpCutLow || refAdcPulseAmp > refAdcPulseAmpCutHigh);
adcRefPulseTimeCut = (refAdcPulseTimeRaw*adcChanToTime < refAdcPulseTimeCutLow || refAdcPulseTimeRaw*adcChanToTime > refAdcPulseTimeCutHigh);
// Implement cuts
//if (adcRefMultiplicityCut || adcRefPulseAmpCut || adcRefPulseTimeCut) continue;
// Acquire the hodoscope ADC data objects
if(signalNames[isignal] == "Adc") {
// Loop over the signals again
for(UInt_t jsignal = 0; jsignal < nSignals; jsignal++) {
// Acquire othe hodoscope TDC data objects
if(signalNames[jsignal] == "Tdc") {
// Loop over the ADC data objects
for(UInt_t iadcdata = 0; iadcdata < nAdcSignals; iadcdata++) {
if (adcData[iadcdata] != "PulseTimeRaw") continue;
// Acquire the number of ADC hits
numAdcHits = adcHits[iplane][iside][isignal];
for (Int_t iadchit = 0; iadchit < numAdcHits; iadchit++) {
// Obtain variables
adcPaddleNum = UInt_t (adcPaddle[iplane][iside][isignal][iadchit]);
adcErrorFlag = hodoAdcErrorFlag[iplane][iside][isignal][iadchit];
adcPulseTimeRaw = hodoAdcPulseTimeRaw[iplane][iside][isignal][iadchit]*adcChanToTime;
adcPulseTime = adcPulseTimeRaw - refAdcPulseTimeRaw*adcChanToTime;
adcPulseAmp = hodoAdcPulseAmp[iplane][iside][isignal][iadchit];
// Define cuts
adcErrorFlagCut = (adcErrorFlag != 0);
adcPulseAmpCut = (adcPulseAmp < hodoPulseAmpCutLow || adcPulseAmp > hodoPulseAmpCutHigh);
// Implement cuts
if (adcErrorFlagCut || adcPulseAmpCut) continue;
// Loop over the TDC data objects
for(UInt_t itdcdata = 0; itdcdata < nTdcSignals; itdcdata++) {
if (tdcData[itdcdata] != "TimeRaw") continue;
// Acquire the number of TDC hits
numTdcHits = tdcHits[iplane][iside][jsignal];
// Define cuts
//edtmCut = (numTdcHits == nbars[0] || numTdcHits == nbars[2] || numTdcHits == nbars[3]);
// Implement cuts
//if (edtmCut) continue;
for (Int_t itdchit = 0; itdchit < numTdcHits; itdchit++) {
// Obtain variables
tdcPaddleNum = UInt_t (tdcPaddle[iplane][iside][jsignal][itdchit]);
tdcTimeRaw = hodoTdcTimeRaw[iplane][iside][jsignal][itdchit]*tdcChanToTime;
if (iplane == 3 && iside == 0) tdcTime = tdcTimeRaw - refT2TdcTimeRaw*tdcChanToTime;
else tdcTime = tdcTimeRaw - refT1TdcTimeRaw*tdcChanToTime;
adcTdcTimeDiff = tdcTime - adcPulseTime;
// Define cuts
adcAndTdcHitCut = (adcPaddleNum != tdcPaddleNum);
adcTdcTimeDiffCut = (adcTdcTimeDiff < adcTdcTimeDiffCutLow || adcTdcTimeDiff > adcTdcTimeDiffCutHigh);
// Implement cuts
if (adcAndTdcHitCut || adcTdcTimeDiffCut) continue;
h2_adcPulseAmpCuts[iplane][iside]->Fill(tdcPaddleNum, adcPulseAmp);
h2_adcTdcTimeDiff[iplane][iside]->Fill(tdcPaddleNum, adcTdcTimeDiff);
h2_adcTimeWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, adcPulseTime);
h2_tdcTimeWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, tdcTime);
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->Fill(adcPulseAmp, adcTdcTimeDiff);
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->GetXaxis()->CenterTitle();
h2_adcTdcTimeDiffWalk[iplane][iside][tdcPaddleNum-1]->GetYaxis()->CenterTitle();
} // TDC hit loop
} // TDC signal loop
} // ADC hit loop
} // ADC signal loop
} // TDC signal
} // Nested signal loop
} // ADC signal
} // Signal loop
} // Side loop
} // Plane loop
if (ievent % 100000 == 0 && ievent != 0)
cout << ievent << " Events Have Been Processed..." << endl;
} // rawDataTree event loop
cout << "\n***************************************" << endl;
cout << ievent << " Events Were Processed" << endl;
cout << "***************************************\n" << endl;
// Calculate the analysis rate
t = clock() - t;
printf ("The Analysis Took %.1f seconds \n", ((float) t) / CLOCKS_PER_SEC);
printf ("The Analysis Event Rate Was %.3f kHz \n", (ievent + 1) / (((float) t) / CLOCKS_PER_SEC*1000.));
outFile->Write();
//outFile->Close();
return 0;
} // time_walk_calib()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment