diff --git a/CALIBRATION/shms_dc_calib/scripts/DC_calib.C b/CALIBRATION/shms_dc_calib/scripts/DC_calib.C
index aa0f9a00f17df2a6b2a131da5932e53576eb5fe7..33de7e7743b777421272dca6503951992ac83146 100644
--- a/CALIBRATION/shms_dc_calib/scripts/DC_calib.C
+++ b/CALIBRATION/shms_dc_calib/scripts/DC_calib.C
@@ -29,6 +29,7 @@ DC_calib::DC_calib(TString a, TString b, const int c, Long64_t d, TString e)
entries = NULL;
t_zero = NULL;
t_zero_err = NULL;
+ t_zero_final = NULL;
bin_max = NULL;
bin_maxContent = NULL;
time_max = NULL;
@@ -67,6 +68,7 @@ DC_calib::~DC_calib()
delete [] entries[ip];
delete [] t_zero[ip];
delete [] t_zero_err[ip];
+ delete [] t_zero_final[ip];
delete [] cell_dt[ip];
delete [] cell_dt_corr[ip];
delete [] fitted_cell_dt[ip];
@@ -80,6 +82,7 @@ DC_calib::~DC_calib()
delete [] entries; entries = NULL;
delete [] t_zero; t_zero = NULL;
delete [] t_zero_err; t_zero_err = NULL;
+ delete [] t_zero_final; t_zero_final = NULL;
delete [] cell_dt; cell_dt = NULL;
delete [] cell_dt_corr; cell_dt_corr = NULL;
delete [] fitted_cell_dt; fitted_cell_dt = NULL;
@@ -175,10 +178,10 @@ void DC_calib::GetDCLeafs()
if (num_evts > nentries)
{
- cout << "Number of entries exceeds: " << nentries << endl;
- cout << "Please input a value that is <= " << nentries << " entries" << endl;
- cout << "Exiting NOW! " << endl;
- exit(EXIT_SUCCESS);
+ cout << "Number of entries entered: " << num_evts << " exeeds MAX number of entries: " << nentries << endl;
+ cout << "Setting the number of entries to: " << nentries << endl;
+
+ num_evts = nentries;
}
@@ -200,25 +203,90 @@ void DC_calib::GetDCLeafs()
}
+
if (spec=="SHMS")
{
- cer_npe_name = "P.ngcer.npeSum";
- EL_CLEAN_name = "T.shms.pEL_CLEAN_tdcTime";
- //EL_CLEAN_name = "T.coin.pEL_CLEAN_ROC2_tdcTime";
+ cal_etotnorm_leaf = "P.cal.etotnorm";
+ cer_npe_leaf = "P.ngcer.npeSum";
+ EL_CLEAN_leaf = "T.shms.pEL_CLEAN_tdcTime";
+ //EL_CLEAN_leaf = "T.coin.pEL_CLEAN_ROC2_tdcTime";
+
+ //Check Branch Status
+ status_cal = tree->GetBranchStatus(cal_etotnorm_leaf);
+ status_cer = tree->GetBranchStatus(cer_npe_leaf);
+ status_EL_clean = tree->GetBranchStatus(EL_CLEAN_leaf);
- tree->SetBranchAddress(cer_npe_name, &cer_npe);
- tree->SetBranchAddress(EL_CLEAN_name, &EL_CLEAN);
+
+ if ((!status_cal || !status_cer || !status_EL_clean) && (pid=="pid_elec" || pid=="pid_prot" || pid=="bkg_cut"))
+ {
+ cout << "*************ATTENTION!**************" << endl;
+ cout << "" << endl;
+ cout << " One or more of the following leafs " << endl;
+ cout << " is *NOT* present in current ROOTfile. " << endl;
+ cout << "1) " << cal_etotnorm_leaf<< endl;
+ cout << "2) " << cer_npe_leaf << endl;
+ cout << "3) " << EL_CLEAN_leaf << endl;
+ cout << "" << endl;
+ cout << "Please add them if you want to make " << endl;
+ cout << "any cuts during calibration." << endl;
+ cout << "" << endl;
+ cout << "OR, set the pid flag in main_calib.C " << endl;
+ cout << "to pid_kFALSE" << endl;
+ cout << " Exiting NOW! " << endl;
+ cout << "*************************************" << endl;
+
+ exit (EXIT_SUCCESS);
+ }
+
+ else
+ {
+ tree->SetBranchAddress(cal_etotnorm_leaf, &cal_etot_norm);
+ tree->SetBranchAddress(cer_npe_leaf, &cer_npe);
+ tree->SetBranchAddress(EL_CLEAN_leaf, &EL_CLEAN);
+ }
+
}
else if (spec=="HMS")
{
- cer_npe_name = "H.cer.npeSum";
- EL_CLEAN_name = "T.hms.hEL_CLEAN_tdcTime";
- //EL_CLEAN_name = "T.coin.hEL_CLEAN_ROC2_tdcTime";
+ cal_etotnorm_leaf = "H.cal.etotnorm";
+ cer_npe_leaf = "H.cer.npeSum";
+ EL_CLEAN_leaf = "T.hms.hEL_CLEAN_tdcTime";
+ //EL_CLEAN_leaf = "T.coin.hEL_CLEAN_ROC2_tdcTime";
+
+ //Check Branch Status with Boolean
+ status_cal = tree->GetBranchStatus(cal_etotnorm_leaf);
+ status_cer = tree->GetBranchStatus(cer_npe_leaf);
+ status_EL_clean = tree->GetBranchStatus(EL_CLEAN_leaf);
+
+ if ((!status_cal || !status_cer || !status_EL_clean) && (pid=="pid_elec" || pid=="pid_prot" || pid=="bkg_cut"))
+ {
+ cout << "*************ATTENTION!**************" << endl;
+ cout << "" << endl;
+ cout << " One or more of the following leafs " << endl;
+ cout << " is *NOT* present in current ROOTfile. " << endl;
+ cout << "1) " << cal_etotnorm_leaf<< endl;
+ cout << "2) " << cer_npe_leaf << endl;
+ cout << "3) " << EL_CLEAN_leaf << endl;
+ cout << "" << endl;
+ cout << "Please add them if you want to make " << endl;
+ cout << "any cuts during calibration." << endl;
+ cout << "" << endl;
+ cout << "OR, set the pid flag in main_calib.C "<< endl;
+ cout << "to: pid_kFALSE" << endl;
+ cout << " Exiting NOW! " << endl;
+ cout << "*************************************" << endl;
+ exit (EXIT_SUCCESS);
+ }
- tree->SetBranchAddress(cer_npe_name, &cer_npe);
- tree->SetBranchAddress(EL_CLEAN_name, &EL_CLEAN);
-
+ else
+ {
+ tree->SetBranchAddress(cal_etotnorm_leaf, &cal_etot_norm);
+ tree->SetBranchAddress(cer_npe_leaf, &cer_npe);
+ tree->SetBranchAddress(EL_CLEAN_leaf, &EL_CLEAN);
+ }
+
+
}
@@ -238,6 +306,7 @@ void DC_calib::AllocateDynamicArrays()
entries = new Int_t*[NPLANES];
t_zero = new Double_t*[NPLANES];
t_zero_err = new Double_t*[NPLANES];
+ t_zero_final = new Double_t*[NPLANES];
cell_dt = new TH1F*[NPLANES]; /*create array to store cell drift times*/
cell_dt_corr = new TH1F*[NPLANES];
fitted_cell_dt = new TH1F*[NPLANES]; /*create array to store cell drift times*/
@@ -253,6 +322,7 @@ void DC_calib::AllocateDynamicArrays()
entries[ip] = new Int_t[nwires[ip]];
t_zero[ip] = new Double_t[nwires[ip]];
t_zero_err[ip] = new Double_t[nwires[ip]];
+ t_zero_final[ip] = new Double_t[nwires[ip]];
cell_dt[ip] = new TH1F[nwires[ip]];
cell_dt_corr[ip] = new TH1F[nwires[ip]];
fitted_cell_dt[ip] = new TH1F[nwires[ip]];
@@ -360,32 +430,46 @@ void DC_calib::EventLoop()
//------READ USER 'pid' input to determine particle type to calibrate----------
- //NO PID Cut, Set Bool_t to always kTRUE
+ //NO PID Cut,
if(pid=="pid_kFALSE")
{
+ //cal_elec = 1;
cer_elec = 1;
elec_clean = 1;
}
- //PID Cut, Set Bool_t to actual value, and see if it passes cut
+ //PID Cut, Set Bool_t to actual leaf value, and see if it passes cut
else if (pid=="pid_elec")
{
- cer_elec = cer_npe>1.0;
- elec_clean = EL_CLEAN>0; //tdcTime>0
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (bkg cleanup)
+ cer_elec = cer_npe>1.0; //number of photoelec. > 1 (electrons)
+ elec_clean = EL_CLEAN>0; //tdcTime>0 (reduce bkg events)
}
//PID Cut, hadron, Set Bool_t to actual value, and see if it passes cut
else if (pid=="pid_hadron")
{
- cer_elec = cer_npe<1.0;
- elec_clean = EL_CLEAN==0; //tdcTime==0
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (bkg cleanup)
+ cer_elec = cer_npe<1.0; //number of photoelec. < 1 (hadrons)
+ elec_clean = EL_CLEAN>0; //tdcTime>0 (reduce bkg events)
}
+ //PID Cut, BACKGROUND, Set Bool_t to actual value, and see if it passes cut
+ else if (pid=="pid_bkg")
+ {
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (bkg cleanup)
+ cer_elec = 1; //Set to always true
+ elec_clean = EL_CLEAN>0; //tdcTime>0 (reduce bkg events)
+ }
+
+
+
else
{
cout << "Enter which particle to calibrate: " << endl;
cout << "For electrons: 'pid_elec' " << endl;
cout << "For hadrons: 'pid_hadron' " << endl;
+ cout << "For background cut (Only EL-CLEAN trigger): 'pid_bkg' " << endl;
cout << "NO PID Cuts: 'pid_KFALSE' " << endl;
}
@@ -399,95 +483,96 @@ void DC_calib::EventLoop()
// cout << "PLANE: " << ip << endl;
//-----------------------------------------------------------------------------------------
-
- //Loop over number of hits for each trigger in each DC plane
- for(Int_t j = 0, k = 0; j < ndata_time[ip], k < ndata_wirenum[ip]; j++, k++)
- {
-
-
- //get wire hit for ith event in 'ip' plane
- wire = int(wire_num[ip][k]);
- //cout << "event: " << i << " ::pEL-Clean: " << pEL_CLEAN << endl;
-
-
-
- //Fill uncorrected plane drift times (from: get_pdc_time_histo.C )
- if (ip== 0 && wire > 80)
- {
- plane_dt[ip].Fill(drift_time[ip][j] - offset);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- }
-
- else if (ip== 6 && wire > 80)
- {
- plane_dt[ip].Fill(drift_time[ip][j] - offset);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- }
-
- else if (ip == 0 && wire <= 80)
- {
- plane_dt[ip].Fill(drift_time[ip][j]);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- }
- else if (ip == 6 && wire <= 80)
- {
- plane_dt[ip].Fill(drift_time[ip][j]);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- }
-
-
-
- else if ((ip == 4 || ip == 10) && wire > 48 && wire <65)
- {
- plane_dt[ip].Fill(drift_time[ip][j] - offset);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
- }
- else if (ip == 4 && (wire <=48 || wire >= 65))
- {
- plane_dt[ip].Fill(drift_time[ip][j]);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- }
- else if (ip == 10 && (wire <= 48 || wire >= 65))
- {
- plane_dt[ip].Fill(drift_time[ip][j]);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- }
-
- else if (ip!=4 || ip!=0 || ip!=6 || ip!=10)
- {
- plane_dt[ip].Fill(drift_time[ip][j]);
- dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
- cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
- }
-
- } // end TRG Cut
+
+ //Require number of chamber hits per event per plane to be UNITY.
+ if (ndata_time[ip]==1 && ndata_wirenum[ip]==1)
+ {
+ //Loop over number of hits for each trigger in each DC plane
+ for(Int_t j = 0, k = 0; j < ndata_time[ip], k < ndata_wirenum[ip]; j++, k++)
+ {
+
+
+ //get wire hit for ith event in 'ip' plane
+ wire = int(wire_num[ip][k]);
+ //cout << "event: " << i << " ::pEL-Clean: " << pEL_CLEAN << endl;
+
+
+
+ //Fill uncorrected plane drift times (from: get_pdc_time_histo.C )
+ if (ip== 0 && wire > 80)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j] - offset);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ }
+
+ else if (ip== 6 && wire > 80)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j] - offset);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ }
+
+ else if (ip == 0 && wire <= 80)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j]);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ }
+ else if (ip == 6 && wire <= 80)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j]);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ }
+
+
+
+ else if ((ip == 4 || ip == 10) && wire > 48 && wire <65)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j] - offset);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j] - offset);
+ }
+ else if (ip == 4 && (wire <=48 || wire >= 65))
+ {
+ plane_dt[ip].Fill(drift_time[ip][j]);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ }
+ else if (ip == 10 && (wire <= 48 || wire >= 65))
+ {
+ plane_dt[ip].Fill(drift_time[ip][j]);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ }
+
+ else if (ip!=4 || ip!=0 || ip!=6 || ip!=10)
+ {
+ plane_dt[ip].Fill(drift_time[ip][j]);
+ dt_vs_wire[ip].Fill(wire_num[ip][k], drift_time[ip][j]);
+ cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ fitted_cell_dt[ip][wire-1].Fill(drift_time[ip][j]);
+ }
+ } //end loop over hits
+
+ } // end if statement (require ONLY 1 hits per event)
//-----------------------------------------------------------------------------------------
- } //END loop over hits
-
-
-
+ } //END loop over planes
- } //end loop over planes
+ } //end loop over pid cuts
} //end loop over events
-} // end event loop
+} // end event loop method
//_________________________________________________________________________
@@ -903,18 +988,18 @@ void DC_calib::WriteTZeroParam()
//cout << "Plane: " << ip << endl;
for (wire=0; wire<nwires[ip]; wire++)
{
- // cout << "wire: " << wire+1 << " :: " << "tzero: "<< t_zero[ip][wire] << endl;
+ //cout << "Plane: " << ip << " :: wire: " << wire+1 << " :: " << "tzero_final: "<< t_zero_final[ip][wire] << endl;
if (wire <= 10)
{
- out_txtFILE << setprecision(6) << t_zero[ip][wire] << fixed << ",";
+ out_txtFILE << setprecision(6) << t_zero_final[ip][wire] << fixed << ",";
}
else if (wire>10 && wire <(nwires[ip]-1))
{
- out_txtFILE << setprecision(6) << t_zero[ip][wire] << ((wire+1) % 16 ? ", " : "\n") << fixed;
+ out_txtFILE << setprecision(6) << t_zero_final[ip][wire] << ((wire+1) % 16 ? ", " : "\n") << fixed;
}
else if (wire==nwires[ip]-1)
{
- out_txtFILE << setprecision(6) << t_zero[ip][wire] << fixed << endl;
+ out_txtFILE << setprecision(6) << t_zero_final[ip][wire] << fixed << endl;
}
} //END LOOP OVER WIRES
@@ -945,7 +1030,7 @@ void DC_calib::ApplyTZeroCorrection()
//PID Cut, Set Bool_t to always kTRUE
if(pid=="pid_kFALSE")
{
-
+ //cal_elec = 1;
cer_elec = 1;
elec_clean = 1;
@@ -954,26 +1039,37 @@ void DC_calib::ApplyTZeroCorrection()
//PID Cut, Set Bool_t to actual value, and see if it passes cut
else if (pid=="pid_elec")
{
-
- cer_elec = cer_npe>1.0;
- elec_clean = EL_CLEAN>0; //tdcTime>0
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (reduce bkg events)
+ cer_elec = cer_npe>1.0; //number of photoelec. > 1 (electrons)
+ elec_clean = EL_CLEAN>0.; //tdcTime>0 (reduce bkg events)
}
//PID Cut, hadron, Set Bool_t to actual value, and see if it passes cut
else if (pid=="pid_hadron")
{
-
- cer_elec = cer_npe<1.0;
- elec_clean = EL_CLEAN==0; //tdcTime==0
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (reduce bkg events)
+ cer_elec = cer_npe<1.0; //number of photoelec. < 1 (hadrons)
+ elec_clean = EL_CLEAN>0.; //tdcTime>0 (reduce bkg events)
}
+ //PID Cut, hadron, Set Bool_t to actual value, and see if it passes cut
+ else if (pid=="pid_bkg")
+ {
+ //cal_elec = cal_etot_norm>0.1; //normalize energy > 0.1 (reduce bkg events)
+ cer_elec = 1; //set to always true
+ elec_clean = EL_CLEAN>0.; //tdcTime>0 (reduce bkg events)
+
+ }
+
+
else
{
cout << "Enter which particle to calibrate: " << endl;
cout << "For electrons: 'pid_elec' " << endl;
cout << "For hadrons: 'pid_hadron' " << endl;
+ cout << "For background cut (Only EL-CLEAN trigger): 'pid_bkg' " << endl;
cout << "NO PID Cuts: 'pid_KFALSE' " << endl;
}
@@ -987,195 +1083,227 @@ void DC_calib::ApplyTZeroCorrection()
{
//cout << "ApplyTZeroCorr: " << weighted_avg[ip] << endl;
//Loop over number of hits for each trigger in each DC plane
-
- for(Int_t j = 0, k = 0; j < ndata_time[ip], k < ndata_wirenum[ip]; j++, k++)
+
+ //Require number of chamber hits per event per plane to be UNITY.
+ if (ndata_time[ip]==1 && ndata_wirenum[ip]==1)
{
- //get wire hit for ith event in 'ip' plane
- wire = int(wire_num[ip][j]);
-
- //Apply general offsets to account for TDC shift
- if (ip== 0 && wire > 80)
+ for(Int_t j = 0, k = 0; j < ndata_time[ip], k < ndata_wirenum[ip]; j++, k++)
{
- //apply weighted average corr to low stats wires
- if ( t_zero[ip][wire-1] == weighted_avg[ip])
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
- }
+ //get wire hit for ith event in 'ip' plane
+ wire = int(wire_num[ip][j]);
- else
+ //Apply general offsets to account for TDC shift
+ if (ip== 0 && wire > 80)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ //apply weighted average corr to low stats wires
+ if ( t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = offset + weighted_avg[ip];
+ //cout << "offset: " << offset << endl;
+ //cout << "weighted avg: " << weighted_avg[ip] << endl;
+ //cout << "t_zero_final: " << t_zero_final[ip][wire-1] << endl;
+ }
+
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = offset + t_zero[ip][wire-1];
+
+ }
}
- }
-
- else if ( ip == 0 && wire <=80)
- {
- if (t_zero[ip][wire-1] == weighted_avg[ip])
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
- }
- else
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
- }
-
- }
-
- //apply tdc offset to plane 6
- else if (ip== 6 && wire > 80)
- {
- if ( t_zero[ip][wire-1] == weighted_avg[ip])
+ else if ( ip == 0 && wire <=80)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
+ if (t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = weighted_avg[ip];
+
+ }
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = t_zero[ip][wire-1];
+ }
+
}
- else
+ //apply tdc offset to plane 6
+ else if (ip== 6 && wire > 80)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ if ( t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = offset + weighted_avg[ip];
+
+ }
+
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = offset + t_zero[ip][wire-1];
+ }
}
- }
-
- else if (ip == 6 && wire <= 80)
- {
- if ( t_zero[ip][wire-1] == weighted_avg[ip])
+ else if (ip == 6 && wire <= 80)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+
+ if ( t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = weighted_avg[ip];
+
+ }
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = t_zero[ip][wire-1];
+ }
+
}
- else
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
- }
- }
-
- else if ((ip == 4 || ip == 10) && wire > 48 && wire <65)
- {
- if (t_zero[ip][wire-1] == weighted_avg[ip])
+ else if ((ip == 4 || ip == 10) && wire > 48 && wire <65)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
- }
- else
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] -offset - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ if (t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = offset + weighted_avg[ip];
+
+ }
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] -offset - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - offset - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = offset + t_zero[ip][wire-1];
+ }
+
}
- }
-
- else if (ip == 4 && (wire <=48 || wire >= 65 ))
- {
- if (t_zero[ip][wire-1] == weighted_avg[ip])
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
- }
- else
+ else if (ip == 4 && (wire <=48 || wire >= 65 ))
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ if (t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = weighted_avg[ip];
+
+ }
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = t_zero[ip][wire-1];
+ }
+
}
- }
-
-
- else if (ip == 10 && (wire <= 48 || wire >= 65))
- {
- if (t_zero[ip][wire-1] == weighted_avg[ip])
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
- }
- else
+ else if (ip == 10 && (wire <= 48 || wire >= 65))
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ if (t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = weighted_avg[ip];
+
+ }
+
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = t_zero[ip][wire-1];
+ }
+
}
- }
-
- else if (ip!=4 || ip!=0 || ip!=6 || ip!=10)
- {
-
- if (t_zero[ip][wire-1] == weighted_avg[ip])
+ else if (ip!=4 || ip!=0 || ip!=6 || ip!=10)
{
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
- }
- else
- {
- //Fill corrected plane drift times
- plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
- dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+
+ if (t_zero[ip][wire-1] == weighted_avg[ip])
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - weighted_avg[ip]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - weighted_avg[ip]);
+ t_zero_final[ip][wire-1] = weighted_avg[ip];
+
+ }
+ else
+ {
+ //Fill corrected plane drift times
+ plane_dt_corr[ip].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ cell_dt_corr[ip][wire-1].Fill(drift_time[ip][j] - t_zero[ip][wire-1]);
+ dt_vs_wire_corr[ip].Fill(wire_num[ip][k], drift_time[ip][j] - t_zero[ip][wire-1]);
+ t_zero_final[ip][wire-1] = t_zero[ip][wire-1];
+ }
+
}
- }
+ } // end LOOP over number of hits
- } // end EL-CLEAN TRG Cut
-
- } //end loop over hits
+ } //end if statement(require #of hits to be UNITY)
- } //end loop over planes
-
+ } //end loop over planes
+
+ } //end loop over PID cuts
+
} //end loop over events
-}
+
+} // end loop over method
//_________________________________________________________________________________
void DC_calib::WriteLookUpTable()
{
otxtfile_name = "./"+spectre+"dc_calib_"+std::to_string(run_NUM)+".param";
out_txtFILE.open(otxtfile_name);
- Double_t t_offset_firstbin = -19.;
+ Double_t t_offset_firstbin = 0.0;
//Set headers for subsequent columns of data
out_txtFILE << Form("; Lookup Table: RUN %d", run_NUM) << "\n";
out_txtFILE << "; number of bins in time to distance lookup table" << "\n";
out_txtFILE << Form(spectre+"driftbins = %d", TOTAL_BINS+1) << "\n";
out_txtFILE << "; number of 1st bin in table in ns" << "\n";
- out_txtFILE << spectre+Form("drift1stbin=%f", t_offset_firstbin) << "\n";
+ out_txtFILE << spectre+Form("drift1stbin=%f", t_offset_firstbin + 1.) << "\n";
out_txtFILE << "; bin size in ns" << "\n";
out_txtFILE << spectre+"driftbinsz=1" << "\n";
@@ -1183,9 +1311,9 @@ void DC_calib::WriteLookUpTable()
for (int ip=0; ip<NPLANES; ip++){
-
+
//Get bin corresponding to t0 = 0 ns
- bin_t0[ip] = plane_dt_corr[ip].GetXaxis()->FindBin(t_offset_firstbin -1.);
+ bin_t0[ip] = plane_dt_corr[ip].GetXaxis()->FindBin(t_offset_firstbin);
//Get final bin
bin_final[ip] = bin_t0[ip] + TOTAL_BINS;
diff --git a/CALIBRATION/shms_dc_calib/scripts/DC_calib.h b/CALIBRATION/shms_dc_calib/scripts/DC_calib.h
index 153d9b022bb881aa96f9b2e8c3485f6c11d0db6c..2ed5e01a211c812896213ef8238901255de30e9d 100644
--- a/CALIBRATION/shms_dc_calib/scripts/DC_calib.h
+++ b/CALIBRATION/shms_dc_calib/scripts/DC_calib.h
@@ -4,9 +4,9 @@
#define NPLANES 12
#define NBINS 400
-#define MINBIN -50.5
-#define MAXBIN 349.5
-#define TOTAL_BINS 274
+#define MINBIN -50.0
+#define MAXBIN 350.0
+#define TOTAL_BINS 189
class DC_calib
{
public:
@@ -56,19 +56,28 @@ class DC_calib
TString drifttime;
TString wirenum;
- TString cer_npe_name;
- TString EL_CLEAN_name;
+ TString cal_etotnorm_leaf;
+ TString cer_npe_leaf;
+ TString EL_CLEAN_leaf;
- Double_t cer_npe;
- Double_t EL_CLEAN;
+ Double_t cal_etot_norm; //calorimeter normalized energy
+ Double_t cer_npe; //cerenkon photoelectron Sum
+ Double_t EL_CLEAN; //electron clean trigger
Double_t hcer_npe;
+ //Boolean for checking if TBranch exists
+ Bool_t status_cal;
+ Bool_t status_cer;
+ Bool_t status_EL_clean;
//Boolean for PID cuts
- Bool_t cer_elec; //hms cerenkov cut
+ Bool_t cal_elec; //calorimeter normalized energy cut
+ Bool_t cer_elec; //cerenkov cut
Bool_t elec_clean; //e- clean trigger tdctime cut
+
+
Int_t wire;
@@ -155,7 +164,11 @@ class DC_calib
Double_t std_dev;
Double_t **t_zero;
Double_t **t_zero_err;
-
+
+ //tzero with tdc offsets taken into account,
+ //to be written into tzero param file
+ Double_t **t_zero_final;
+
//declare variables to make plot of tzero v. wire number
Double_t weighted_avg[NPLANES];
diff --git a/CALIBRATION/shms_dc_calib/scripts/main_calib.C b/CALIBRATION/shms_dc_calib/scripts/main_calib.C
index baa58aa954a19df573d689c95daf88474386487a..47d3ced161696575292c63a5ce41201ba93fc43b 100644
--- a/CALIBRATION/shms_dc_calib/scripts/main_calib.C
+++ b/CALIBRATION/shms_dc_calib/scripts/main_calib.C
@@ -16,11 +16,12 @@ int main_calib()
clock_t cl;
cl = clock();
- //pid_elec, pid_hadron, pid_kFALSE (no PID cuts)
+ //pid_elec, pid_hadron, pid_bkg pid_kFALSE (no PID cuts)
// |
// v
- DC_calib obj("HMS", "../../../ROOTfiles/coin_replay_production_1866_-1_dcuncalib.root", 1866,1400000, "pid_hadron");
-
+ DC_calib obj("HMS", "../../../ROOTfiles/hms_replay_production_all_1856_dcuncal.root", 1856,2000000, "pid_kFALSE");
+// DC_calib obj("SHMS", "../../../ROOTfiles/shms_replay_production_all_2071_-1_dcuncalib.root", 2071, 3000000, "pid_bkg");
+
obj.printInitVar();
obj.SetPlaneNames();
@@ -29,12 +30,11 @@ int main_calib()
obj.CreateHistoNames();
obj.EventLoop();
obj.Calculate_tZero();
- obj.WriteTZeroParam();
obj.ApplyTZeroCorrection();
+ obj.WriteTZeroParam();
obj.WriteLookUpTable();
obj.WriteToFile(1); //set argument to (1) for debugging
-
-
+
//stop clock
cl = clock() - cl;