diff --git a/cfg/moduleCharacterization.cfg b/cfg/moduleCharacterization.cfg
index 6fd6dacb..5197c27f 100644
--- a/cfg/moduleCharacterization.cfg
+++ b/cfg/moduleCharacterization.cfg
@@ -1,19 +1,42 @@
inputDir /data/TOFPET2/reco/
fileBaseName run
-runs 2436
+runs 2672
maxEntries -1
-step1FileName /home/cmsdaq/martina/Lab5015Analysis/plots/moduleCharacterization_step1_run2436.root
+step1FileName /data/ALIO/plots/moduleCharacterization_step1_run2672.root
+#step1FileName /data/ALIO/plots/moduleCharacterizationCoinNa22_step1_run2634_prova.root
+#step1FileName /home/cmsdaq/guglielmi/Lab5015Analysis/plots/analyzeTOFPET2_wirelessBarInArray_HDR2_UVlaser_tune93_step1_run2652.root
usePedestals 0
-sourceName Na22
+sourceName Co60
+#sourceName Co60SumPeak
+#sourceName Na22
+#sourceName Na22SingleBar
+#sourceName Laser
+SiPM HDR2
+#SiPM FBK
@@ -23,11 +46,22 @@ tResMax 250
tResMode 2
Vov 3 5 7
-energyBins 800 600 400
+energyBinsNa22 800 600 400
+energyBinsCo60 1000 750 500
+energyBinsLaser 800 600 400
energyMin 0
-energyMax 40
-
+energyMaxNa22 40
+energyMaxCo60 50
+energyMaxLaser 35
refBar 7
refVov 7
refTh 20
+
+
+status no
+chL 226
+chR 244
+peak511eBin 5
+
+
diff --git a/cfg/moduleCharacterization_step3.cfg b/cfg/moduleCharacterization_step3.cfg
index 37a678b1..eddf79c2 100644
--- a/cfg/moduleCharacterization_step3.cfg
+++ b/cfg/moduleCharacterization_step3.cfg
@@ -1,19 +1,55 @@
-inputDir /home/cmsdaq/martina/Lab5015Analysis/plots/
+inputDir1 /data/ALIO/plots/
+inputDir /data/ALIO/plots/
fileBaseName1 moduleCharacterization_step
+#fileBaseName1 moduleCharacterizationCoinNa22_step
fileBaseName2 run
-#runs 2428-2453
+runs 2428-2453
#runs 2548-2574
-runs 2375-2401
+#runs 2375-2401
+#runs 2608-2634
+#runs 2670-2672
+#runs 2674-2678
+#runs 2313-2315
+#runs 2316-2318
+#runs 2304-2312
+#runs 2326-2334
+#sourceName Co60
+#sourceName Co60SumPeak
+sourceName Na22
+#sourceName SingleBarNa22
+#sourceName SingleBarNa22_coinc
+
+
@@ -21,6 +57,8 @@ tResMin 0
tResMax 250
tResMode 2
+
+
Vov 3 5 7
energyBins 800 600 400
energyMin 0
@@ -29,6 +67,18 @@ energyMax 40
refBar 7
refVov 7
refTh 20
+#refTh_HDR2
+#refTh3 15
+#refTh5 15
+#refTh7 20
+#refTh_SingleHDR2
+#refTh3 20
+#refTh5 20
+#refTh7 20
+#refTh_FBK
+refTh3 10
+refTh5 15
+refTh7 20
diff --git a/interface/Co60SpectrumAnalyzer_2Peaks.h b/interface/Co60SpectrumAnalyzer_2Peaks.h
new file mode 100644
index 00000000..5935ed23
--- /dev/null
+++ b/interface/Co60SpectrumAnalyzer_2Peaks.h
@@ -0,0 +1,18 @@
+#ifndef CO60SPECTRUMANALYZER_2PEAKS_H
+#define CO60SPECTRUMANALYZER_2PEAKS_H
+
+#include
+#include
+
+#include "TH1F.h"
+#include "TF1.h"
+#include "TLine.h"
+#include "TSpectrum.h"
+#include "TStyle.h"
+
+
+
+std::map > Co60SpectrumAnalyzer_2Peaks(TH1F* histo,
+ std::vector* ranges = NULL);
+
+#endif
diff --git a/interface/Na22SpectrumAnalyzerSingleBar.h b/interface/Na22SpectrumAnalyzerSingleBar.h
new file mode 100644
index 00000000..bdd7a010
--- /dev/null
+++ b/interface/Na22SpectrumAnalyzerSingleBar.h
@@ -0,0 +1,18 @@
+#ifndef NA22SPECTRUMANALYZERSINGLEBAR_H
+#define NA22SPECTRUMANALYZERSINGLEBAR_H
+
+#include
+#include
+
+#include "TH1F.h"
+#include "TF1.h"
+#include "TLine.h"
+#include "TSpectrum.h"
+#include "TStyle.h"
+
+
+
+std::map > Na22SpectrumAnalyzerSingleBar(TH1F* histo,
+ std::vector* ranges = NULL);
+
+#endif
diff --git a/interface/SiPM_HDR2.h b/interface/SiPM_HDR2.h
index 928a3856..9884d115 100644
--- a/interface/SiPM_HDR2.h
+++ b/interface/SiPM_HDR2.h
@@ -1,30 +1,36 @@
-// HDR2 SiPM parameters
-
-/*
-float PDE_vs_OV(const float& ov)
+float PDE_vs_OV(const float& ov, const std::string& type = "HDR2")
{
- return 39.4321 * ( 1. - exp(-1.*0.738063*ov) );
+ if( type == "HDR2" )
+ return 0.394321 * ( 1. - exp(-1.*0.738063*ov) );
+
+ if( type == "FBK_W7S" )
+ return 0.391748 * ( 1. - exp(-1.*0.288262*ov) );
+
+ return 0.;
}
-float Gain_vs_OV(const float& ov)
+float Gain_vs_OV(const float& ov, const std::string& type = "HDR2")
{
- return 36890.225 + 97602.904*ov;
+ if( type == "HDR2" )
+ return 36890.225 + 97602.904*ov;
+
+ if( type == "FBK_W7S" )
+ //return 29131.588 + 104567.57 * ov;
+ //return 54174.9 + 95169.2 * ov;
+ return 50738.5 + 95149 * ov;
+
+
+ return 0.;
}
-*/
-
-
-// FBK SiPM parameters
-float PDE_vs_OV(const float& ov)
+float ECF_vs_OV(const float& ov, const std::string& type = "HDR2")
{
-
- //return ( -1.5389053 + 22.812361*ov - 10.091946*ov*ov + 3.6867230 *ov*ov*ov -0.73282017*ov*ov*ov*ov + 0.056861423*ov*ov*ov*ov*ov ) ;
-
- return 45.5275 * ( 1. - exp(-1.*0.36229*ov) );
-}
-
-float Gain_vs_OV(const float& ov)
-{
- return 29131.588 + 104567.57 * ov;
-
+ if( type == "HDR2" )
+ return 1. + 0.0030812903 * ov + 0.0015623938 * ov * ov;
+
+ if( type == "FBK_W7S" )
+ //return 1. - 0.006356446 * ov + 0.0057041163 * ov * ov;
+ return 1.02366 - 0.00168754 * ov + 0.00284079 * ov * ov;
+
+ return 0.;
}
diff --git a/macros/plotLaser_guglielmi.py b/macros/plotLaser_guglielmi.py
new file mode 100644
index 00000000..f026b26c
--- /dev/null
+++ b/macros/plotLaser_guglielmi.py
@@ -0,0 +1,1091 @@
+#! /usr/bin/env python
+import os
+import shutil
+import glob
+import math
+import array
+import sys
+import time
+import numpy as np
+
+
+import ROOT
+from ROOT import *
+import CMS_lumi, tdrstyle
+
+from ROOT import TLatex
+
+#set the tdr style
+tdrstyle.setTDRStyle()
+ROOT.gStyle.SetOptStat(0)
+ROOT.gStyle.SetOptFit(1)
+ROOT.gStyle.SetOptTitle(0)
+ROOT.gStyle.SetLabelSize(0.04,'X')
+ROOT.gStyle.SetLabelSize(0.04,'Y')
+ROOT.gStyle.SetTitleSize(0.04,'X')
+ROOT.gStyle.SetTitleSize(0.04,'Y')
+ROOT.gStyle.SetTitleOffset(1.1,'X')
+ROOT.gStyle.SetTitleOffset(1.2,'Y')
+ROOT.gStyle.SetImageScaling(2.5);
+
+
+
+def GetLinearizedEnergy(x, par):
+ xx = x[0]
+ if( xx < par[0] ): return ( par[1]*(math.exp(par[2]*xx)-1) + par[3] );
+ else: return ( par[1]*(math.exp(par[2]*par[0])-1)+par[3])/(math.exp(par[4]*par[0]))*math.exp(par[4]*xx );
+
+
+def GetSaturationCorrection(Ncells, Edep, LY, PDE, LCE):
+ Npe = Edep * LY * PDE * LCE
+ Nfired = Ncells * (1 - math.exp(-Npe/Ncells))
+ k = Nfired/Npe
+ #print Edep, LY, PDE, LCE, ECF
+ #print 'Ncells = %d, Npe = %f, Nfired = %f, k = %f'%(Ncells, Npe, Nfired, k)
+ return k
+
+
+
+
+
+#to run python plotLaser_guglielmi.py /home/cmsdaq/guglielmi/Lab5015Analysis/ /data/ALIO/plots/ /home/cmsdaq/guglielmi/Lab5015Analysis/plots/
+
+###
+outdir = '/var/www/html/ModuleCharacterization/guglielmi/Laser/new_funcTot/'
+path = os.path.join(outdir)
+#os.mkdir(path)
+plotdir = sys.argv[1] + 'plots/'
+step1dir = sys.argv[2]
+otherdir = sys.argv[3]
+
+Vov = 3.0
+LY = 40000. # gamma/MeV
+LCE = 0.15
+Ncells = 40000
+PDE = 0.351
+
+
+
+#laserTunes = [10,20,30,40,50,60,70,80,83,85,86,87,88,89,90]
+laserTunes = [20,30,40,50,60,70,75,80,85,87,90,91,92,93,94]
+run =[2665, 2664, 2663, 2662, 2661, 2660, 2659, 2658, 2655, 2657, 2654, 2653, 2656, 2652, 2651]
+
+f_3_20_Na22 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2611.root'%step1dir)
+f_5_20_Na22 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2620.root'%step1dir)
+f_7_20_Na22 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2629.root'%step1dir)
+
+h_en_Vov3_th20_Na22 = f_3_20_Na22.Get('h1_energy_bar00L-R_Vov3.0_th20')
+h_en_Vov5_th20_Na22 = f_5_20_Na22.Get('h1_energy_bar00L-R_Vov5.0_th20')
+h_en_Vov7_th20_Na22 = f_7_20_Na22.Get('h1_energy_bar00L-R_Vov7.0_th20')
+
+f_3_20_Co60 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2670.root'%step1dir)
+f_5_20_Co60 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2671.root'%step1dir)
+f_7_20_Co60 = ROOT.TFile.Open('%smoduleCharacterization_step1_run2672.root'%step1dir)
+
+h_en_Vov3_th20_Co60 = f_3_20_Co60.Get('h1_energy_bar00L-R_Vov3.0_th20')
+h_en_Vov5_th20_Co60 = f_5_20_Co60.Get('h1_energy_bar00L-R_Vov5.0_th20')
+h_en_Vov7_th20_Co60 = f_7_20_Co60.Get('h1_energy_bar00L-R_Vov7.0_th20')
+
+
+
+energy = {}
+linearizedEnergy = {}
+f = {}
+fstep2 = {}
+h1_energy = {}
+g_tRes_vs_threshold = {}
+g_enPeak_vs_bar_Vov3_th20 = {}
+
+g_tRes_vs_energy = ROOT.TGraphErrors()
+g_tRes_vs_linearizedEnergy = ROOT.TGraphErrors()
+g_tRes_vs_energy_bestTh = ROOT.TGraphErrors()
+g_tRes_vs_linearizedEnergy_bestTh = ROOT.TGraphErrors()
+g_energy_vs_laserTune = ROOT.TGraphErrors()
+
+tRes_err_sys = 2.
+energy_err_sys = 0.015 # 1.5%
+
+
+# outfile root
+outFileName = ('%splotLaser.root' %plotdir)
+outFile = ROOT.TFile.Open(str(outFileName),"RECREATE");
+
+
+
+
+# from Andrea
+funcLinearizedEnergy = ROOT.TF1('funcLinearizedEnergy',GetLinearizedEnergy, 0, 100,5)
+funcLinearizedEnergy.SetParameter(0,1.79617e+01)
+funcLinearizedEnergy.SetParameter(1,7.07013e+00)
+funcLinearizedEnergy.SetParameter(2,1.87292e-02)
+funcLinearizedEnergy.SetParameter(3,-1.70555e-02)
+funcLinearizedEnergy.SetParameter(4,6.59613e-02)
+
+funcLinearizedEnergyCo60 = ROOT.TF1('funcLinearizedEnergyCo60',GetLinearizedEnergy, 0, 100,5)
+funcLinearizedEnergyCo60.SetParameter(0,1.79617e+01)
+funcLinearizedEnergyCo60.SetParameter(1,7.07013e+00)
+funcLinearizedEnergyCo60.SetParameter(2,1.87292e-02)
+funcLinearizedEnergyCo60.SetParameter(3,-1.70555e-02)
+funcLinearizedEnergyCo60.SetParameter(4,6.59613e-02)
+
+
+
+
+
+
+#Linearization
+f1 = ROOT.TFile.Open('%smoduleCharacterization_single_HPK_HDR2_step3.root'%plotdir)
+f2 = ROOT.TFile.Open('%smoduleCharacterization_single_HPK_HDR2_Co60_step3.root'%plotdir)
+
+'''funcLinearizedEnergy = TF1('f_linearity_Vov%.01f' %Vov, '[0]*(exp([1]*x)-1) + [2]' ,0.,200.)
+funcLinearizedEnergy.SetParameters(6.92,0.0158,-0.0639)
+funcLinearizedEnergyCo60 = TF1('f_linearity_Vov%.01f_Co60' %Vov, '[0]*(exp([1]*x)-1) + [2]' ,0.,200.)
+funcLinearizedEnergyCo60.SetParameters(6.92,0.0158,-0.0639)'''
+
+
+
+g_Linearization_Na22 = ROOT.TGraphErrors()
+g_Linearization_Na22 = f1.Get('g_linearity_Vov%.01f' %Vov)
+g_Linearization_Co60 = ROOT.TGraphErrors()
+g_Linearization_Co60 = f2.Get('g_linearity_Vov%.01f' %Vov)
+g_Vov3_Na22 = ROOT.TGraphErrors()
+g_Vov3_Na22 = f1.Get('g_ov3_Vov%.01f' %Vov)
+g_Vov5_Na22 = ROOT.TGraphErrors()
+g_Vov5_Na22 = f1.Get('g_ov5_Vov%.01f' %Vov)
+g_Vov7_Na22 = ROOT.TGraphErrors()
+g_Vov7_Na22 = f1.Get('g_ov7_Vov%.01f' %Vov)
+g_Vov3_Co60 = ROOT.TGraphErrors()
+g_Vov3_Co60 = f2.Get('g_ov3_Vov%.01f' %Vov)
+g_Vov5_Co60 = ROOT.TGraphErrors()
+g_Vov5_Co60 = f2.Get('g_ov5_Vov%.01f' %Vov)
+g_Vov7_Co60 = ROOT.TGraphErrors()
+g_Vov7_Co60 = f2.Get('g_ov7_Vov%.01f' %Vov)
+
+
+
+g_Linearization_Co60_rescaled = ROOT.TGraphErrors()
+g_Vov3_Co60_rescaled = ROOT.TGraphErrors()
+g_Vov5_Co60_rescaled = ROOT.TGraphErrors()
+g_Vov7_Co60_rescaled = ROOT.TGraphErrors()
+for i in range(g_Linearization_Co60.GetN()):
+ g_Linearization_Co60_rescaled.SetPoint(i,g_Linearization_Co60.GetPointX(i)/0.93, g_Linearization_Co60.GetPointY(i))
+ #print(g_Linearization_Co60.GetPointX(i)/0.93)
+for i in range(g_Vov3_Co60.GetN()):
+ g_Vov3_Co60_rescaled.SetPoint(i,g_Vov3_Co60.GetPointX(i)/0.93, g_Vov3_Co60.GetPointY(i))
+ g_Vov5_Co60_rescaled.SetPoint(i,g_Vov5_Co60.GetPointX(i)/0.93, g_Vov5_Co60.GetPointY(i))
+ g_Vov7_Co60_rescaled.SetPoint(i,g_Vov7_Co60.GetPointX(i)/0.93, g_Vov7_Co60.GetPointY(i))
+
+'''g_Vov3_Co60_rescaled.SetPoint(0,7.57231769017414, g_Vov3_Co60.GetPointY(0))
+g_Vov3_Co60_rescaled.SetPoint(1,8.59874438475017, g_Vov3_Co60.GetPointY(1))
+g_Vov3_Co60_rescaled.SetPoint(2,16.1710620749243, g_Vov3_Co60.GetPointY(2))
+#g_Vov3_Co60_rescaled.SetPoint(2,g_Vov3_Co60.GetPointX(2), g_Vov3_Co60.GetPointY(2))
+g_Linearization_Co60_rescaled.SetPoint(0,7.57231769017414, g_Vov3_Co60.GetPointY(0))
+g_Linearization_Co60_rescaled.SetPoint(1,8.59874438475017, g_Vov3_Co60.GetPointY(1))
+g_Linearization_Co60_rescaled.SetPoint(2,16.1710620749243, g_Vov3_Co60.GetPointY(2))
+#g_Linearization_Co60_rescaled.SetPoint(2,g_Vov3_Co60.GetPointX(2), g_Vov3_Co60.GetPointY(2))
+
+g_Vov5_Co60_rescaled.SetPoint(0,13.2141438415333, g_Vov5_Co60.GetPointY(0))
+g_Vov5_Co60_rescaled.SetPoint(1,15.0053193494649, g_Vov5_Co60.GetPointY(1))
+g_Vov5_Co60_rescaled.SetPoint(2,28.2194631909982, g_Vov5_Co60.GetPointY(2))
+#g_Vov5_Co60_rescaled.SetPoint(2,g_Vov5_Co60.GetPointX(2), g_Vov5_Co60.GetPointY(2))
+g_Linearization_Co60_rescaled.SetPoint(3,13.2141438415333, g_Vov5_Co60.GetPointY(0))
+g_Linearization_Co60_rescaled.SetPoint(4,15.0053193494649, g_Vov5_Co60.GetPointY(1))
+g_Linearization_Co60_rescaled.SetPoint(5,28.2194631909982, g_Vov5_Co60.GetPointY(2))
+g_Linearization_Co60_rescaled.SetPoint(5,g_Vov5_Co60.GetPointX(2), g_Vov5_Co60.GetPointY(2))
+
+g_Vov7_Co60_rescaled.SetPoint(0,18.033908680328, g_Vov7_Co60.GetPointY(0))
+g_Vov7_Co60_rescaled.SetPoint(1,20.4784026958201, g_Vov7_Co60.GetPointY(1))
+g_Vov7_Co60_rescaled.SetPoint(2,38.5123113761481, g_Vov7_Co60.GetPointY(2))
+#g_Vov7_Co60_rescaled.SetPoint(2,g_Vov7_Co60.GetPointX(2), g_Vov7_Co60.GetPointY(2))
+g_Linearization_Co60_rescaled.SetPoint(6,18.033908680328, g_Vov7_Co60.GetPointY(0))
+g_Linearization_Co60_rescaled.SetPoint(7,20.4784026958201, g_Vov7_Co60.GetPointY(1))
+g_Linearization_Co60_rescaled.SetPoint(8,38.5123113761481, g_Vov7_Co60.GetPointY(2))
+#g_Linearization_Co60_rescaled.SetPoint(8,g_Vov7_Co60.GetPointX(2), g_Vov7_Co60.GetPointY(2))'''
+
+
+c_colour_rescaled = ROOT.TCanvas('Lin coloured scaled','Lin coloured scaled',700,500)
+hPad_lin10 = ROOT.gPad.DrawFrame(0.,0.,30.,10.)
+hPad_lin10.GetYaxis().SetTitleSize(0.04)
+hPad_lin10.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hPad_lin10.GetYaxis().SetTitle('E_{MeV} #times G #times PDE #times ECF #times k_{sat}/[G #times PDE #times ECF]_{%.1f V}' %Vov)
+hPad_lin10.Draw()
+gPad.SetGridx()
+gPad.SetGridy()
+g_Vov3_Na22.SetMarkerStyle(20)
+g_Vov3_Na22.SetMarkerSize(1.)
+g_Vov3_Na22.SetMarkerColor(kGreen)
+g_Vov3_Na22.Draw("P")
+g_Vov5_Na22.SetMarkerStyle(20)
+g_Vov5_Na22.SetMarkerSize(1.)
+g_Vov5_Na22.SetMarkerColor(kRed)
+g_Vov5_Na22.Draw("Psame")
+g_Vov7_Na22.SetMarkerStyle(20)
+g_Vov7_Na22.SetMarkerSize(1.)
+g_Vov7_Na22.SetMarkerColor(kBlue)
+g_Vov7_Na22.Draw("Psame")
+g_Vov3_Co60_rescaled.SetMarkerStyle(20)
+g_Vov3_Co60_rescaled.SetMarkerSize(1.)
+g_Vov3_Co60_rescaled.SetMarkerColor(kYellow)
+g_Vov3_Co60_rescaled.Draw("Psame")
+g_Vov5_Co60_rescaled.SetMarkerStyle(20)
+g_Vov5_Co60_rescaled.SetMarkerSize(1.)
+g_Vov5_Co60_rescaled.SetMarkerColor(kViolet)
+g_Vov5_Co60_rescaled.Draw("Psame")
+g_Vov7_Co60_rescaled.SetMarkerStyle(20)
+g_Vov7_Co60_rescaled.SetMarkerSize(1.)
+g_Vov7_Co60_rescaled.SetMarkerColor(kBlack)
+g_Vov7_Co60_rescaled.Draw("Psame")
+'''g_Linearization_Co60_rescaled.SetMarkerStyle(20)
+g_Linearization_Co60_rescaled.SetMarkerSize(1.)
+g_Linearization_Co60_rescaled.SetMarkerColor(kYellow)
+g_Linearization_Co60_rescaled.Draw("Psame")'''
+leg_rescaled = ROOT.TLegend(0.33,0.73,0.53,0.85)
+leg_rescaled.SetFillStyle(0)
+leg_rescaled.SetBorderSize(0)
+leg_rescaled.AddEntry(g_Vov3_Na22,'Na22 Vov3','PL')
+leg_rescaled.AddEntry(g_Vov5_Na22,'Na22 Vov5','PL')
+leg_rescaled.AddEntry(g_Vov7_Na22,'Na22 Vov7','PL')
+leg_rescaled.AddEntry(g_Vov3_Co60_rescaled,'Co60 Vov3','PL')
+leg_rescaled.AddEntry(g_Vov5_Co60_rescaled,'Co60 Vov5','PL')
+leg_rescaled.AddEntry(g_Vov7_Co60_rescaled,'Co60 Vov7','PL')
+leg_rescaled.Draw('same')
+
+
+
+
+x = []
+y = []
+n = g_Linearization_Na22.GetN() + g_Linearization_Co60.GetN()
+for i in range( n ):
+ if i %d'%(g_tRes_vs_energy_Na22.GetPointX(eBin-1), laserTune[eBin])
+
+#find laser tune that gives energy closest to sumPeak (2.5MeV) Co60
+for eBin in refPeaks_Co60:
+#for eBin in range(1, 13):
+ d = 999999
+ for i,tune in enumerate(laserTunes):
+ if (abs(g_tRes_vs_energy_Co60.GetPointX(eBin-1)-energy[tune]) < d):
+ laserTune[eBin]= tune
+ d = abs(g_tRes_vs_energy_Co60.GetPointX(eBin-1)-energy[tune])
+ print 'laser tune for %d a.u. Co60 --> %d'%(g_tRes_vs_energy_Co60.GetPointX(eBin-1), laserTune[eBin])
+
+# check diff in quadrature of Na22 - laser vs threshold
+#g_diff_vs_threshold = {}
+#for eBin in refPeaks:
+# g_diff_vs_threshold[eBin] = ROOT.TGraphErrors()
+# for i in range(0, g_tRes_vs_threshold[eBin].GetN()):
+# resNa22 = ROOT.Double(0)
+# th = ROOT.Double(0)
+# g_tRes_vs_threshold[eBin].GetPoint(i, th, resNa22)
+# resLaser = g_tRes_vs_threshold[laserTune[eBin]].Eval(th)
+# diff = -1
+# if (resNa22 > resLaser):
+# diff = math.sqrt(resNa22*resNa22-resLaser*resLaser)
+# g_diff_vs_threshold[eBin].SetPoint(g_diff_vs_threshold[eBin].GetN(), th, diff)
+
+### PLOTS
+
+tl = ROOT.TLatex( 0.16, 0.93, 'LYSO:Ce 3x3x57 mm^{2} - HDR2 - TOFPET2' )
+tl.SetNDC()
+tl.SetTextSize(0.030)
+
+
+# plot laser energy
+canvas1 = ROOT.TCanvas('c_energy_laser','c_energy_laser')
+canvas1.SetGridy()
+canvas1.SetGridx()
+hdummy1 = ROOT.TH2F('hdummy1','', 100, 0, 30, 100, 0.1, 4000)
+hdummy1.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy1.GetYaxis().SetTitle('counts')
+hdummy1.GetYaxis().SetTitleOffset(1.5)
+hdummy1.Draw()
+for i,tune in enumerate(laserTunes):
+ h1_energy[tune].SetLineWidth(2)
+ h1_energy[tune].SetLineColor(51 + i*3)
+ h1_energy[tune].Draw('histo same')
+tl.Draw()
+
+
+
+# plot laser tune vs energy at th20
+canvas4 = ROOT.TCanvas('c_energy_vs_laserTune','c_energy_vs_laserTune')
+canvas4.SetGridy()
+canvas4.SetGridx()
+hdummy4 = ROOT.TH2F('hdummy4','', 100, 0, 100, 100, 0, 30)
+hdummy4.GetYaxis().SetTitle('energy')
+hdummy4.GetXaxis().SetTitle('laser tune')
+hdummy4.Draw()
+g_energy_vs_laserTune.Draw('plsame')
+tl.Draw()
+
+# plot time resolution vs threshold for laser and Na22 and Co60
+canvas3 = ROOT.TCanvas('c_tRes_vs_threshold','c_tRes_vs_threshold')
+canvas3.Divide(2)
+#canvas3.SetGridy()
+#canvas3.SetGridx()
+hdummy3 = ROOT.TH2F('hdummy3','', 100, 0, 120, 250, 0, 700)
+hdummy3.GetXaxis().SetTitle('threshold (DAC)')
+hdummy3.GetYaxis().SetTitle('#sigma_{bar} (ps)')
+
+legend1 = ROOT.TLegend(0.75,0.2,0.95,0.9)
+legend1.SetBorderSize(0)
+legend2 = ROOT.TLegend(0.75,0.2,0.95,0.9)
+legend2.SetBorderSize(0)
+canvas3.cd(1)
+hPad1 = ROOT.gPad.DrawFrame(0.,0.,30.,7.);
+hPad1.Draw();
+gPad.SetGridx();
+gPad.SetGridy();
+hdummy3.Draw()
+for i,tune in enumerate(laserTunes):
+ g_tRes_vs_threshold[tune].SetMarkerColor(i+1)
+ g_tRes_vs_threshold[tune].SetLineColor(i+1)
+ if i == 9:
+ g_tRes_vs_threshold[tune].SetMarkerColor(len(laserTunes)+1)
+ g_tRes_vs_threshold[tune].SetLineColor(len(laserTunes)+1)
+ g_tRes_vs_threshold[tune].SetMarkerStyle(26)
+ g_tRes_vs_threshold[tune].Draw('plsame')
+ legend1.AddEntry(g_tRes_vs_threshold[tune], 'laser tune %d'%tune, 'pl')
+legend1.Draw()
+tl.Draw()
+canvas3.cd(2)
+hPad2 = ROOT.gPad.DrawFrame(0.,0.,2.,500.);
+hPad2.Draw();
+gPad.SetGridx();
+gPad.SetGridy();
+hdummy3.Draw()
+#g_tRes_vs_threshold[5].SetMarkerColor(eBin)
+#g_tRes_vs_threshold[5].SetLineColor(eBin)
+#g_tRes_vs_threshold[5].SetMarkerStyle(24)
+#g_tRes_vs_threshold[5].Draw('plsame')
+for eBin in range(2,13):
+ g_tRes_vs_threshold[eBin].SetMarkerColor(eBin)
+ g_tRes_vs_threshold[eBin].SetLineColor(eBin)
+ g_tRes_vs_threshold[eBin].SetMarkerStyle(24)
+ g_tRes_vs_threshold[eBin].Draw('plsame')
+ legend2.AddEntry(g_tRes_vs_threshold[eBin], 'Na22 enBin %d'%eBin, 'pl')
+g_tRes_vs_threshold_Co60.SetMarkerColor(kOrange)
+g_tRes_vs_threshold_Co60.SetLineColor(kOrange)
+g_tRes_vs_threshold_Co60.SetMarkerStyle(24)
+g_tRes_vs_threshold_Co60.Draw('plsame')
+legend2.AddEntry(g_tRes_vs_threshold_Co60, 'Co60 enBin SumPeak', 'pl')
+
+'''g_tRes_vs_threshold_Na22_coinc.SetMarkerColor(kViolet)
+g_tRes_vs_threshold_Na22_coinc.SetLineColor(kViolet)
+g_tRes_vs_threshold_Na22_coinc.SetMarkerStyle(24)
+g_tRes_vs_threshold_Na22_coinc.Draw('plsame')
+legend2.AddEntry(g_tRes_vs_threshold_Na22_coinc, 'Na22 enBin 5 Coinc', 'pl')'''
+
+legend2.Draw()
+tl.Draw()
+
+# plot time resolution vs threshold for laser and Na22
+#canvas7 = ROOT.TCanvas('c_diff_vs_threshold','c_diff_vs_threshold')
+#canvas7.SetGridy()
+#canvas7.SetGridx()
+#hdummy7 = ROOT.TH2F('hdummy7','', 100, 0, 64, 250, 0, 250)
+#hdummy7.GetXaxis().SetTitle('threshold (DAC)')
+#hdummy7.GetYaxis().SetTitle('diff (ps)')
+#hdummy7.Draw()
+#for i,peak in enumerate(refPeaks):
+# g_diff_vs_threshold[eBin].SetMarkerColor(i+1)
+# g_diff_vs_threshold[eBin].SetLineColor(i+1)
+# g_diff_vs_threshold[eBin].SetMarkerStyle(21)
+# g_diff_vs_threshold[eBin].Draw('plsame')
+#tl.Draw()
+
+
+# plot time resolution vs energy th20
+canvas2 = ROOT.TCanvas('c_tRes_vs_energy','c_tRes_vs_energy')
+canvas2.SetGridy()
+canvas2.SetGridx()
+hdummy2 = ROOT.TH2F('hdummy2','', 100, 0, 30, 250, 0, 500)
+hdummy2.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy2.GetYaxis().SetTitle('#sigma_{bar} (ps)')
+hdummy2.Draw()
+g_tRes_vs_energy.SetLineWidth(2)
+g_tRes_vs_energy.Draw('plsame')
+g_tRes_vs_energy_Na22.SetMarkerStyle(21)
+g_tRes_vs_energy_Na22.SetMarkerSize(0.7)
+g_tRes_vs_energy_Na22.SetMarkerColor(ROOT.kRed)
+g_tRes_vs_energy_Na22.SetLineColor(ROOT.kRed)
+g_tRes_vs_energy_Na22.Draw('psame')
+g_tRes_vs_energy_Na22_coinc.SetMarkerStyle(21)
+g_tRes_vs_energy_Na22_coinc.SetMarkerSize(0.7)
+g_tRes_vs_energy_Na22_coinc.SetMarkerColor(ROOT.kBlue)
+g_tRes_vs_energy_Na22_coinc.SetLineColor(ROOT.kBlue)
+g_tRes_vs_energy_Na22_coinc.Draw('psame')
+g_tRes_vs_energy_Co60.SetMarkerStyle(21)
+g_tRes_vs_energy_Co60.SetMarkerSize(0.7)
+g_tRes_vs_energy_Co60.SetMarkerColor(ROOT.kGreen)
+g_tRes_vs_energy_Co60.SetLineColor(ROOT.kGreen)
+g_tRes_vs_energy_Co60.Draw('psame')
+
+#leg = ROOT.TLegend(0.15,0.13,0.33,0.28)
+leg = ROOT.TLegend(0.33,0.73,0.53,0.85)
+leg.SetFillStyle(0)
+leg.SetBorderSize(0)
+leg.AddEntry(g_tRes_vs_energy,'UV laser','PL')
+leg.AddEntry(g_tRes_vs_energy_Na22,'^{22}Na','PL')
+leg.AddEntry(g_tRes_vs_energy_Na22_coinc,'^{22}Na coincidence','PL')
+leg.AddEntry(g_tRes_vs_energy_Co60,'^{60}Co','PL')
+leg.Draw('same')
+tl.Draw()
+
+#plot function for linearized energy
+canvas6 = ROOT.TCanvas('c_function_linearity','c_function_linearity')
+canvas6.SetGridy()
+canvas6.SetGridx()
+hdummy6 = ROOT.TH2F('hdummy6','', 30, 0, 30, 6, 0, 6)
+hdummy6.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy6.GetYaxis().SetTitle('energy (MeV)')
+hdummy6.Draw()
+funcLinearizedEnergy.SetLineWidth(2)
+funcLinearizedEnergy.Draw('same')
+
+# plot time resolution vs linearized energy at th20
+canvas5 = ROOT.TCanvas('c_tRes_vs_linearizedEnergy','c_tRes_vs_linearizedEnergy')
+canvas5.SetGridy()
+canvas5.SetGridx()
+hdummy5 = ROOT.TH2F('hdummy5','', 100, 0, 7.0, 250, 0, 500)
+hdummy5.GetXaxis().SetTitle('energy (MeV)')
+hdummy5.GetYaxis().SetTitle('#sigma_{bar} (ps)')
+hdummy5.Draw()
+g_tRes_vs_linearizedEnergy.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy.Draw('psame')
+'''fitFunc = ROOT.TF1('fitFunc','sqrt([0]*[0] + [1]*[1]/pow(x,[2])/pow(x,[2]))', 0.0,10.0)
+fitFunc.SetRange(0.2,6.7)
+fitFunc.SetLineWidth(2)
+fitFunc.SetParName(0,'c')
+fitFunc.SetParName(1,'s')
+fitFunc.SetParName(2,'#alpha')
+fitFunc.SetLineColor(1)
+fitFunc.SetParameter(0, 15)
+fitFunc.SetParameter(1, 50)
+fitFunc.SetParameter(2, 0.5)
+g_tRes_vs_linearizedEnergy.Fit('fitFunc','QRS')'''
+fitFunc = ROOT.TF1('fitFunc','sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))',0.0,10.0)
+fitFunc.SetRange(0.1, 7.0)
+fitFunc.SetParName(0,'n')
+fitFunc.SetParName(1,'s')
+fitFunc.SetParName(2,'c')
+fitFunc.SetLineColor(kBlack)
+fitFunc.SetParameter(0, 50)
+fitFunc.SetParameter(1, 100)
+fitFunc.FixParameter(2, 15)
+fitFunc.SetParLimits(0, 0.00001, 1000)
+fitFunc.SetParLimits(1, 0.00001, 1000)
+
+fitFunc2 = ROOT.TF1('fitFunc2','sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))',0.0,10.0)
+fitFunc2.SetRange(0.1, 7.0)
+fitFunc2.SetParName(0,'n')
+fitFunc2.SetParName(1,'s')
+fitFunc2.SetParName(2,'c')
+fitFunc2.SetLineColor(kBlack)
+fitFunc2.SetParameter(0, 50)
+fitFunc2.SetParameter(1, 100)
+fitFunc2.FixParameter(2, 15)
+fitFunc2.SetParLimits(0, 0.00001, 1000)
+fitFunc2.SetParLimits(1, 0.00001, 1000)
+#fitFunc.SetParLimits(2, 0.00001, 1000)
+g_tRes_vs_linearizedEnergy.Fit('fitFunc','QRS')
+fitFunc2.SetLineColor(kRed)
+g_tRes_vs_linearizedEnergy_Na22.Fit('fitFunc2','QRS+')
+canvas5.Modified()
+canvas5.Update()
+gPad.Update()
+st = g_tRes_vs_linearizedEnergy.GetListOfFunctions().FindObject('stats')
+st.SetX1NDC(0.62)
+st.SetX2NDC(0.98)
+st.SetY1NDC(0.59)
+st.SetY2NDC(0.79)
+g_tRes_vs_linearizedEnergy_Na22.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Na22.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Na22.SetMarkerColor(ROOT.kRed)
+g_tRes_vs_linearizedEnergy_Na22.SetLineColor(ROOT.kRed)
+g_tRes_vs_linearizedEnergy_Na22.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Na22.Draw('psame')
+g_tRes_vs_linearizedEnergy_Na22_coinc.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Na22_coinc.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Na22_coinc.SetMarkerColor(ROOT.kBlue)
+g_tRes_vs_linearizedEnergy_Na22_coinc.SetLineColor(ROOT.kBlue)
+g_tRes_vs_linearizedEnergy_Na22_coinc.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Na22_coinc.Draw('psame')
+g_tRes_vs_linearizedEnergy_Co60.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerColor(kGreen)
+g_tRes_vs_linearizedEnergy_Co60.Draw('psame')
+canvas5.Modified()
+canvas5.Update()
+leg.Draw('same')
+tl.Draw()
+
+print 'Time resolution at 4.2 MeV = %.01f ps'%fitFunc.Eval(4.2)
+vline = ROOT.TLine(4.2,0,4.2,fitFunc.Eval(4.2))
+vline.SetLineStyle(2)
+vline.Draw('same')
+
+hline = ROOT.TLine(0,fitFunc.Eval(4.2),4.2,fitFunc.Eval(4.2))
+hline.SetLineStyle(2)
+hline.Draw('same')
+
+
+
+
+
+# print difference in quadrature laser vs Na22
+g_diff_vs_linearizedEnergy = ROOT.TGraphErrors()
+
+for ipeak,eBin in enumerate(refPeaks):
+ diff = -1
+ diff = math.sqrt(fabs(g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1))*g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)) - g_tRes_vs_energy.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1))* g_tRes_vs_energy.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)) ))
+ print '%d a.u. --> laser = %.02f ps Na22 = %.02f ps diff = %.02f ps '%(g_tRes_vs_energy_Na22.GetPointX(eBin-1), g_tRes_vs_energy.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)), g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)), diff)
+ g_diff_vs_linearizedEnergy.SetPoint(g_diff_vs_linearizedEnergy.GetN(), g_tRes_vs_linearizedEnergy_Na22.GetPointX(eBin-1), diff)
+ diff_err = g_tRes_vs_energy_Na22.GetEY()[ipeak]
+ g_diff_vs_linearizedEnergy.SetPointError(g_diff_vs_linearizedEnergy.GetN()-1, 0, diff_err)
+
+# print difference in quadrature Na22 vs Na22coincidence
+g_diff_vs_linearizedEnergy_Na22Coinc = ROOT.TGraphErrors()
+
+refPeaksCoin = [5]
+
+for ipeak,eBin in enumerate(refPeaksCoin):
+ diff = -1
+ diff = math.sqrt(fabs(g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1))*g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)) - g_tRes_vs_energy_Na22_coinc.Eval(g_tRes_vs_energy_Na22_coinc.GetPointX(eBin-1))* g_tRes_vs_energy_Na22_coinc.Eval(g_tRes_vs_energy_Na22_coinc.GetPointX(eBin-1)) ))
+ print '%d a.u. --> Na22Coinc = %.02f ps Na22 = %.02f ps diff = %.02f ps '%(g_tRes_vs_energy_Na22.GetPointX(eBin-1), g_tRes_vs_energy_Na22_coinc.Eval(g_tRes_vs_energy_Na22_coinc.GetPointX(eBin-1)), g_tRes_vs_energy_Na22.Eval(g_tRes_vs_energy_Na22.GetPointX(eBin-1)), diff)
+
+
+
+
+# print difference in quadrature laser vs Co60
+g_diff_vs_linearizedEnergy_Co60 = ROOT.TGraphErrors()
+for ipeak,eBin in enumerate(refPeaks_Co60):
+#for eBin in refPeaks:
+ diff = -1
+ #if ( g_tRes_vs_energy_Co60.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1)) > g_tRes_vs_energy.Eval(energy[laserTune[eBin]])):
+ diff = math.sqrt(fabs(g_tRes_vs_energy_Co60.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1))*g_tRes_vs_energy_Co60.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1)) - g_tRes_vs_energy.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1))* g_tRes_vs_energy.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1)) ))
+ print '%d a.u. --> laser = %.02f ps Co60 = %.02f ps diff = %.02f ps '%(g_tRes_vs_energy_Co60.GetPointX(eBin-1), g_tRes_vs_energy.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1)), g_tRes_vs_energy_Co60.Eval(g_tRes_vs_energy_Co60.GetPointX(eBin-1)), diff)
+ g_diff_vs_linearizedEnergy_Co60.SetPoint(g_diff_vs_linearizedEnergy_Co60.GetN(), g_tRes_vs_linearizedEnergy_Co60.GetPointX(eBin-1), diff)
+ diff_err = g_tRes_vs_energy_Co60.GetEY()[ipeak]
+ g_diff_vs_linearizedEnergy_Co60.SetPointError(g_diff_vs_linearizedEnergy_Co60.GetN()-1, 0, diff_err)
+
+
+# plot time resolution diff vs linearized energy
+canvas8 = ROOT.TCanvas('c_diff_vs_linearizedEnergy','c_diff_vs_linearizedEnergy')
+canvas8.SetGridy()
+canvas8.SetGridx()
+hdummy8 = ROOT.TH2F('hdummy8','', 100, 0, 6.5, 250, 0, 150)
+hdummy8.GetXaxis().SetTitle('energy (MeV)')
+hdummy8.GetYaxis().SetTitle('#sigma_{bar}^{source} (-) #sigma_{bar}^{laser} (ps)')
+hdummy8.Draw()
+g_diff_vs_linearizedEnergy.SetLineWidth(2)
+g_diff_vs_linearizedEnergy.SetMarkerSize(0.8)
+g_diff_vs_linearizedEnergy.SetMarkerStyle(20)
+g_diff_vs_linearizedEnergy.SetLineColor(kRed)
+g_diff_vs_linearizedEnergy.Draw('psame')
+g_diff_vs_linearizedEnergy_Co60.SetLineWidth(2)
+g_diff_vs_linearizedEnergy_Co60.SetMarkerSize(0.8)
+g_diff_vs_linearizedEnergy_Co60.SetMarkerStyle(20)
+g_diff_vs_linearizedEnergy_Co60.SetLineColor(kGreen)
+g_diff_vs_linearizedEnergy_Co60.Draw('psame')
+
+
+# plot time resolution vs energy at bestTh
+canvas9 = ROOT.TCanvas('c_tRes_vs_energy_bestTh','c_tRes_vs_energy_bestTh')
+canvas9.SetGridy()
+canvas9.SetGridx()
+hdummy9 = ROOT.TH2F('hdummy2','', 100, 0, 30, 250, 0, 500)
+hdummy9.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy9.GetYaxis().SetTitle('#sigma_{bar} (ps)')
+hdummy9.Draw()
+g_tRes_vs_energy_bestTh.SetLineWidth(2)
+g_tRes_vs_energy_bestTh.Draw('plsame')
+g_tRes_vs_energy_Na22_bestTh.SetMarkerStyle(21)
+g_tRes_vs_energy_Na22_bestTh.SetMarkerSize(0.7)
+g_tRes_vs_energy_Na22_bestTh.SetMarkerColor(ROOT.kRed)
+g_tRes_vs_energy_Na22_bestTh.SetLineColor(ROOT.kRed)
+g_tRes_vs_energy_Na22_bestTh.Draw('psame')
+g_tRes_vs_energy_Na22_coinc_bestTh.SetMarkerStyle(21)
+g_tRes_vs_energy_Na22_coinc_bestTh.SetMarkerSize(0.7)
+g_tRes_vs_energy_Na22_coinc_bestTh.SetMarkerColor(ROOT.kBlue)
+g_tRes_vs_energy_Na22_coinc_bestTh.SetLineColor(ROOT.kBlue)
+g_tRes_vs_energy_Na22_coinc_bestTh.Draw('psame')
+g_tRes_vs_energy_Co60.SetMarkerStyle(21)
+g_tRes_vs_energy_Co60.SetMarkerSize(0.7)
+g_tRes_vs_energy_Co60.SetMarkerColor(ROOT.kGreen)
+g_tRes_vs_energy_Co60.SetLineColor(ROOT.kGreen)
+g_tRes_vs_energy_Co60.Draw('psame')
+
+#leg = ROOT.TLegend(0.15,0.13,0.33,0.28)
+leg9 = ROOT.TLegend(0.33,0.73,0.53,0.85)
+leg9.SetFillStyle(0)
+leg9.SetBorderSize(0)
+leg9.AddEntry(g_tRes_vs_energy_bestTh,'UV laser','PL')
+leg9.AddEntry(g_tRes_vs_energy_Na22_bestTh,'^{22}Na','PL')
+leg9.AddEntry(g_tRes_vs_energy_Na22_coinc_bestTh,'^{22}Na coincidence','PL')
+leg9.AddEntry(g_tRes_vs_energy_Co60,'^{60}Co','PL')
+leg9.Draw('same')
+tl.Draw()
+
+# plot time resolution vs linearized energy at bestTh
+canvas10 = ROOT.TCanvas('c_tRes_vs_linearizedEnergy_bestTh','c_tRes_vs_linearizedEnergy_bestTh')
+canvas10.SetGridy()
+canvas10.SetGridx()
+hdummy10 = ROOT.TH2F('hdummy10','', 100, 0, 7.0, 250, 0, 500)
+hdummy10.GetXaxis().SetTitle('energy (MeV)')
+hdummy10.GetYaxis().SetTitle('#sigma_{bar} (ps)')
+hdummy10.Draw()
+g_tRes_vs_linearizedEnergy_bestTh.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_bestTh.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_bestTh.Draw('psame')
+fitFunc = ROOT.TF1('fitFunc','sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))',0.0,10.0)
+fitFunc.SetRange(0.25,6.7)
+fitFunc.SetParName(0,'n')
+fitFunc.SetParName(1,'s')
+fitFunc.SetParName(2,'c')
+fitFunc.SetLineColor(kBlack)
+fitFunc.SetParameter(0, 50)
+fitFunc.SetParameter(1, 100)
+fitFunc.SetParameter(2, 50)
+g_tRes_vs_linearizedEnergy_bestTh.Fit('fitFunc','QRS')
+'''fitFunc1 = ROOT.TF1('fitFunc1','sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))',0.0,10.0)
+fitFunc1.SetRange(0.25,6.7)
+fitFunc1.SetParName(0,'n')
+fitFunc1.SetParName(1,'s')
+fitFunc1.SetParName(2,'c')
+fitFunc1.SetLineColor(kBlack)
+fitFunc1.SetParameter(0, 50)
+fitFunc1.SetParameter(1, 100)
+fitFunc1.SetParameter(2, 50)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.Fit('fitFunc1','QRS')'''
+canvas10.Modified()
+canvas10.Update()
+st1 = g_tRes_vs_linearizedEnergy_bestTh.GetListOfFunctions().FindObject('stats')
+st1.SetX1NDC(0.6)
+st1.SetX2NDC(0.90)
+st1.SetY1NDC(0.7)
+st1.SetY2NDC(0.89)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.SetMarkerColor(ROOT.kRed)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.SetLineColor(ROOT.kRed)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Na22_bestTh.Draw('psame')
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.SetMarkerColor(ROOT.kBlue)
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.SetLineColor(ROOT.kBlue)
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Na22_coinc_bestTh.Draw('psame')
+g_tRes_vs_linearizedEnergy_Co60.SetLineWidth(2)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerSize(0.7)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerStyle(21)
+g_tRes_vs_linearizedEnergy_Co60.SetMarkerColor(kGreen)
+g_tRes_vs_linearizedEnergy_Co60.Draw('psame')
+canvas10.Modified()
+canvas10.Update()
+leg9.Draw('same')
+tl.Draw()
+
+#energy Na22, Co60 at th20
+c_en_Na22 = ROOT.TCanvas('c_en_Na22','c_en_Na22')
+c_en_Na22.SetGridy()
+c_en_Na22.SetGridx()
+hdummy_en_Na22 = ROOT.TH2F('hdummy_en_Na22','', 100, 0, 30, 100, 1, 140000)
+hdummy_en_Na22.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy_en_Na22.GetYaxis().SetTitle('entries')
+hdummy_en_Na22.Draw()
+h_en_Vov3_th20_Na22.SetLineColor(kRed)
+h_en_Vov5_th20_Na22.SetLineColor(kBlue)
+h_en_Vov7_th20_Na22.SetLineColor(kGreen)
+h_en_Vov3_th20_Na22.Draw('same')
+h_en_Vov5_th20_Na22.Draw('same')
+h_en_Vov7_th20_Na22.Draw('same')
+c_en_Na22.Modified()
+c_en_Na22.Update()
+tl.Draw()
+
+c_en_Co60 = ROOT.TCanvas('c_en_Co60','c_en_Co60')
+c_en_Co60.SetGridy()
+c_en_Co60.SetGridx()
+c_en_Co60.cd()
+hdummy_en_Co60 = ROOT.TH2F('hdummy_en_Co60','', 100, 0, 30, 100, 1, 200000)
+hdummy_en_Co60.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy_en_Co60.GetYaxis().SetTitle('entries')
+hdummy_en_Co60.Draw()
+h_en_Vov3_th20_Co60.SetLineColor(kRed)
+h_en_Vov5_th20_Co60.SetLineColor(kBlue)
+h_en_Vov7_th20_Co60.SetLineColor(kGreen)
+h_en_Vov3_th20_Co60.Draw('same')
+h_en_Vov5_th20_Co60.Draw('same')
+h_en_Vov7_th20_Co60.Draw('same')
+c_en_Co60.Modified()
+c_en_Co60.Update()
+tl.Draw()
+
+# plot source energy
+canvas11 = ROOT.TCanvas('c_energy_sources','c_energy_sources')
+canvas11.SetGridy()
+canvas11.SetGridx()
+canvas11.SetLogy()
+hdummy11 = ROOT.TH2F('hdummy1','', 100, 0, 30, 100, 1, 500000)
+hdummy11.GetXaxis().SetTitle('TOFPET2 amp (a.u.)')
+hdummy11.GetYaxis().SetTitle('counts')
+hdummy11.Draw()
+h_en_Vov3_th20_Na22.SetLineColor(kRed)
+h_en_Vov3_th20_Na22.SetMarkerColor(kRed)
+h_en_Vov3_th20_Na22.SetLineWidth(2)
+h_en_Vov3_th20_Na22.Draw('histo same')
+#h1_energy[1173].Scale(300./600.)
+h_en_Vov3_th20_Co60.SetLineColor(kBlack)
+h_en_Vov3_th20_Co60.SetLineWidth(2)
+h_en_Vov3_th20_Co60.Draw('histo same')
+leg_source = ROOT.TLegend(0.73,0.73,0.85,0.85)
+leg_source.SetFillStyle(0)
+leg_source.SetBorderSize(0)
+leg_source.AddEntry(h_en_Vov3_th20_Na22,'Na22','L')
+leg_source.AddEntry(h_en_Vov3_th20_Co60,'Co60','L')
+leg_source.Draw('same')
+tl.Draw()
+
+
+
+
+# save plots
+print 'Saving plots in %s'%outdir
+os.system('mkdir %s'%outdir)
+shutil.copy('/var/www/html/index.php', '%s'%outdir)
+
+
+for c in [canvas1, canvas2, canvas3, canvas4, canvas5, canvas6, canvas8, canvas9, canvas10, canvas11, c_lin, c_colour, c_colour_rescaled, c_en_Na22, c_en_Co60]:
+ c.SaveAs(outdir+'/'+c.GetName()+'.png')
+ c.SaveAs(outdir+'/'+c.GetName()+'.pdf')
+ c.SaveAs(outdir+'/'+c.GetName()+'.C')
+
+outFile.Close()
diff --git a/main/moduleCharacterization_step1.cpp b/main/moduleCharacterization_step1.cpp
index 293fb479..719b57bc 100644
--- a/main/moduleCharacterization_step1.cpp
+++ b/main/moduleCharacterization_step1.cpp
@@ -1,6 +1,8 @@
#include "interface/AnalysisUtils.h"
#include "interface/FitUtils.h"
#include "interface/SetTDRStyle.h"
+#include "interface/Na22SpectrumAnalyzer.h"
+#include "interface/Na22SpectrumAnalyzerSingleBar.h"
#include "CfgManager/interface/CfgManager.h"
#include "CfgManager/interface/CfgManagerT.h"
@@ -45,7 +47,6 @@ int main(int argc, char** argv)
//--- parse the config file
CfgManager opts;
opts.ParseConfigFile(argv[1]);
-
int debugMode = 0;
if( argc > 2 ) debugMode = atoi(argv[2]);
@@ -58,6 +59,7 @@ int main(int argc, char** argv)
int maxEntries = opts.GetOpt("Input.maxEntries");
int usePedestals = opts.GetOpt("Input.usePedestals");
TChain* tree = new TChain("data","data");
+
std::stringstream ss(runs);
std::string token;
@@ -89,13 +91,63 @@ int main(int argc, char** argv)
}
}
-
+
+
+
+std::string source = opts.GetOpt("Input.sourceName");
+std::string Na22 = "Na22";
+std::string Na22SingleBar = "Na22SingleBar";
+std::string Co60 = "Co60";
+std::string Co60SumPeak = "Co60SumPeak";
+std::string Laser = "Laser";
+
//--- define channels
- int chsL[16];
- int chsR[16];
+
+int chsL[16];
+int chsR[16];
+
+if(!source.compare(Na22SingleBar) || !source.compare(Co60) || !source.compare(Co60SumPeak) || !source.compare(Laser)){
+ chsL[0] = 249;
+ //chsR[0] = 4;
+ chsR[0] = 49;
+ chsL[1] = 242;
+ chsR[1] = 6;
+ chsL[2] = 251;
+ chsR[2] = 15;
+ chsL[3] = 253;
+ chsR[3] = 8;
+ chsL[4] = 250;
+ chsR[4] = 13;
+ chsL[5] = 254;
+ chsR[5] = 2;
+ chsL[6] = 201;
+ chsR[6] = 11;
+ chsL[7] = 230;
+ chsR[7] = 27;
+ chsL[8] = 226;
+ chsR[8] = 32;
+ chsL[9] = 225;
+ chsR[9] = 31;
+ chsL[10] = 228;
+ chsR[10] = 30;
+ chsL[11] = 227;
+ chsR[11] = 29;
+ chsL[12] = 229;
+ chsR[12] = 28;
+ chsL[13] = 231;
+ chsR[13] = 26;
+ chsL[14] = 233;
+ chsR[14] = 24;
+ chsL[15] = 232;
+ chsR[15] = 25;
+}
+
+
+if(!source.compare(Na22)){
chsL[0] = 249;
chsR[0] = 4;
+ //chsR[0] = 49;
chsL[1] = 242;
chsR[1] = 6;
chsL[2] = 251;
@@ -125,8 +177,9 @@ int main(int argc, char** argv)
chsL[14] = 233;
chsR[14] = 24;
chsL[15] = 232;
- chsR[15] = 25;
-
+ chsR[15] = 25;
+}
+
//--- define branches
@@ -147,13 +200,37 @@ int main(int argc, char** argv)
//--- get plot settings
std::vector Vov = opts.GetOpt >("Plots.Vov");
- std::vector energyBins = opts.GetOpt >("Plots.energyBins");
+ std::vector energyBins;
+ if(!opts.GetOpt("Input.sourceName").compare("Na22") || !opts.GetOpt("Input.sourceName").compare("Na22SingleBar")){
+ energyBins = opts.GetOpt >("Plots.energyBinsNa22");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Co60")){
+ energyBins = opts.GetOpt >("Plots.energyBinsCo60");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Co60SumPeak")){
+ energyBins = opts.GetOpt >("Plots.energyBinsCo60");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Laser")){
+ energyBins = opts.GetOpt >("Plots.energyBinsLaser");
+ }
std::map map_energyBins;
for(unsigned int ii = 0; ii < Vov.size(); ++ii)
map_energyBins[Vov[ii]] = energyBins[ii];
float energyMin = opts.GetOpt("Plots.energyMin");
- float energyMax = opts.GetOpt("Plots.energyMax");
+ float energyMax;
+ if(!opts.GetOpt("Input.sourceName").compare("Na22") || !opts.GetOpt("Input.sourceName").compare("Na22SingleBar")){
+ energyMax = opts.GetOpt("Plots.energyMaxNa22");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Co60")){
+ energyMax = opts.GetOpt("Plots.energyMaxCo60");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Co60SumPeak")){
+ energyMax = opts.GetOpt("Plots.energyMaxCo60");
+ }
+ if(!opts.GetOpt("Input.sourceName").compare("Laser")){
+ energyMax = opts.GetOpt("Plots.energyMaxLaser");
+ }
@@ -171,11 +248,222 @@ int main(int argc, char** argv)
std::map h1_energyL;
std::map h1_energyR;
std::map h1_energyLR;
+ std::map h1_energyLR_ext;
+ std::map c;
+ std::map*> rangesLR;
+ std::map > > peaksLR;
+ std::map acceptEvent;
+ //Coincidence pre loop
+ if(!opts.GetOpt("Coincidence.status").compare("yes") && (!opts.GetOpt("Input.sourceName").compare("Na22SingleBar") || !opts.GetOpt("Input.sourceName").compare("Na22"))){
+
+ float qfineL_ext;
+ float qfineR_ext;
+ float energyL_ext;
+ float energyR_ext;
+ int chL_ext = opts.GetOpt("Coincidence.chL");
+ int chR_ext = opts.GetOpt("Coincidence.chR");
+
+ int nEntries = tree->GetEntries();
+ if( maxEntries > 0 ) nEntries = maxEntries;
+ for(int entry = 0; entry < nEntries; ++entry){
+ tree -> GetEntry(entry);
+ acceptEvent[entry] = false;
+ if( entry%200000 == 0 ){
+ std::cout << "\n>>> external bar loop: reading entry " << entry << " / " << nEntries << " (" << 100.*entry/nEntries << "%)" << std::endl;
+ TrackProcess(cpu, mem, vsz, rss);
+ }
+
+ float Vov = step1;
+ float vth1 = float(int(step2/10000)-1);;
+
+ qfineL_ext = qfine[chL_ext];
+ qfineR_ext = qfine[chR_ext];
+
+ energyL_ext = energy[chL_ext];
+ energyR_ext = energy[chR_ext];
+
+ if (qfineL_ext < 13) continue;
+ if (qfineR_ext < 13) continue;
+
+ int index( (10000*int(Vov*100.)) + (100*vth1) + 99 );
+
+ //--- create histograms, if needed
+ if( h1_energyLR_ext[index] == NULL ){
+ c[index] = new TCanvas(Form("c1_Vov%.1f_th%02.0f",Vov,vth1), Form("c1_Vov%.1f_th%02.0f",Vov,vth1));
+ c[index] -> cd();
+ h1_energyLR_ext[index] = new TH1F(Form("h1_energy_external_barL-R_Vov%.1f_th%02.0f",Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax);
+ }
+ acceptEvent[entry] = true;
+ h1_energyLR_ext[index] -> Fill(0.5*(energyL_ext + energyR_ext));
+ }
+ std::cout << std::endl;
+
+
+ for( auto index : h1_energyLR_ext){
+ rangesLR[index.first] = new std::vector;
+ if(!opts.GetOpt("Input.sourceName").compare("Na22SingleBar")) peaksLR[index.first] = Na22SpectrumAnalyzerSingleBar(index.second,rangesLR[index.first]);
+ if(!opts.GetOpt("Input.sourceName").compare("Na22")) peaksLR[index.first] = Na22SpectrumAnalyzer(index.second,rangesLR[index.first]);
+ }
+ }
+
+
+ //------------------------
+ //--- 1st loop over events
+
+if (!opts.GetOpt("Input.sourceName").compare("Na22") || !opts.GetOpt("Input.sourceName").compare("Na22SingleBar") || !opts.GetOpt("Input.sourceName").compare("Laser")){
+
+
+ ModuleEventClass anEvent;
+
+ float qfineL[16];
+ float qfineR[16];
+ float totL[16];
+ float totR[16];
+ long long int timeL[16];
+ long long int timeR[16];
+ float energyL[16];
+ float energyR[16];
+ int chL[16];
+ int chR[16];
+
+
+ /*float qfineL[1];
+ float qfineR[1];
+ float totL[1];
+ float totR[1];
+ long long int timeL[1];
+ long long int timeR[1];
+ float energyL[1];
+ float energyR[1];
+ int chL[1];
+ int chR[1];*/
+
+ int nEntries = tree->GetEntries();
+ if( maxEntries > 0 ) nEntries = maxEntries;
+ for(int entry = 0; entry < nEntries; ++entry)
+ {
+ tree -> GetEntry(entry);
+ if( entry%200000 == 0 )
+ {
+ std::cout << "\n>>> 1st loop: reading entry " << entry << " / " << nEntries << " (" << 100.*entry/nEntries << "%)" << std::endl;
+ TrackProcess(cpu, mem, vsz, rss);
+ }
+
+ float Vov = step1;
+ float vth1 = float(int(step2/10000)-1);;
+ // float vth2 = float(int((step2-10000*vth1)/100)-1);
+ // float vthe = float(int((step2-10000*vth1-step2-100*vth2)/1)-1);
+
+ if(!opts.GetOpt("Coincidence.status").compare("yes")){
+ if( !acceptEvent[entry] ) continue;
+ }
+ for(int iBar = 0; iBar < 16; ++iBar)
+ {
+
+ chL[iBar]=chsL[iBar];
+ chR[iBar]=chsR[iBar];
+
+ qfineL[iBar]=qfine[chL[iBar]];
+ qfineR[iBar]=qfine[chR[iBar]];
+ totL[iBar]=0.001*tot[chL[iBar]];
+ totR[iBar]=0.001*tot[chR[iBar]];
+
+ energyL[iBar]=energy[chL[iBar]];
+ energyR[iBar]=energy[chR[iBar]];
+ timeL[iBar]=time[chL[iBar]];
+ timeR[iBar]=time[chR[iBar]];
+
+ }
+
+ int maxEn=0;
+ int maxBar=0;
+
+ if(!opts.GetOpt("Coincidence.status").compare("yes")){
+ int chL_ext = opts.GetOpt("Coincidence.chL");
+ int chR_ext = opts.GetOpt("Coincidence.chR");
+ float energyL_ext = energy[chL_ext];
+ float energyR_ext = energy[chR_ext];
+
+ int label = (10000*int(Vov*100.)) + (100*vth1) + 99;
+ int eBin = opts.GetOpt("Coincidence.peak511eBin");
+ float avEn = 0.5 * ( energyL_ext + energyR_ext);
+ //std::cout<<"rangesLR[label]-> at(eBin-1)"< at(eBin-1)< at(eBin)"< at(eBin)< at(eBin-1)) continue;
+ if ( avEn > rangesLR[label]-> at(eBin)) continue;
+ }
+
+
+ for(int iBar = 0; iBar < 16; ++iBar)
+ {
+
+ if (qfineL[iBar]>13 && qfineR[iBar]>13){
+ float energyMean=(energyL[iBar]+energyR[iBar])/2;
+ if(energyMean>maxEn){
+ maxEn = energyMean;
+ maxBar = iBar;
+ }
+ }
+
+ int index( (10000*int(Vov*100.)) + (100*vth1) + iBar );
+
+ //--- create histograms, if needed
+ if( h1_totL[index] == NULL )
+ {
+ h1_qfineL[index] = new TH1F(Form("h1_qfine_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5);
+ h1_qfineR[index] = new TH1F(Form("h1_qfine_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",512,-0.5,511.5);
+
+ h1_totL[index] = new TH1F(Form("h1_tot_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.);
+ h1_totR[index] = new TH1F(Form("h1_tot_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",2000,0.,1000.);
+
+ h1_energyL[index] = new TH1F(Form("h1_energy_bar%02dL_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax);
+ h1_energyR[index] = new TH1F(Form("h1_energy_bar%02dR_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax);
+
+ outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1));
+ outTrees[index] -> Branch("event",&anEvent);
+
+ h1_energyLR[index] = new TH1F(Form("h1_energy_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),"",map_energyBins[Vov],energyMin,energyMax);
+ }
+
+ }
+ int index( (10000*int(Vov*100.)) + (100*vth1) + maxBar );
+
+ //--- fill histograms
+ h1_qfineL[index] -> Fill( qfineL[maxBar] );
+ h1_totL[index] -> Fill( totL[maxBar] );
+ h1_energyL[index] -> Fill( energyL[maxBar] );
+
+ h1_qfineR[index] -> Fill( qfineR[maxBar] );
+ h1_totR[index] -> Fill( totR[maxBar] );
+ h1_energyR[index] -> Fill( energyR[maxBar] );
+
+ h1_energyLR[index] -> Fill(0.5*(energyL[maxBar]+energyR[maxBar]));
+
+ anEvent.barID = maxBar;
+ anEvent.Vov = Vov;
+ anEvent.vth1 = vth1;
+ anEvent.energyL = energyL[maxBar];
+ anEvent.energyR = energyR[maxBar];
+ anEvent.timeL = timeL[maxBar];
+ anEvent.timeR = timeR[maxBar];
+ outTrees[index] -> Fill();
+
+ }
+ std::cout << std::endl;
- //------------------------
- //--- 1st loop over events
+ int bytes = outFile -> Write();
+ std::cout << "============================================" << std::endl;
+ std::cout << "nr of B written: " << int(bytes) << std::endl;
+ std::cout << "nr of KB written: " << int(bytes/1024.) << std::endl;
+ std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl;
+ std::cout << "============================================" << std::endl;
+
+}
+
+if (!opts.GetOpt("Input.sourceName").compare("Co60") || !opts.GetOpt("Input.sourceName").compare("Co60SumPeak")){
ModuleEventClass anEvent;
int nEntries = tree->GetEntries();
@@ -268,3 +556,8 @@ int main(int argc, char** argv)
std::cout << "nr of MB written: " << int(bytes/1024./1024.) << std::endl;
std::cout << "============================================" << std::endl;
}
+
+
+
+
+}
diff --git a/main/moduleCharacterization_step2.cpp b/main/moduleCharacterization_step2.cpp
index 33c12e47..22aee81c 100644
--- a/main/moduleCharacterization_step2.cpp
+++ b/main/moduleCharacterization_step2.cpp
@@ -1,6 +1,7 @@
#include "interface/AnalysisUtils.h"
#include "interface/Na22SpectrumAnalyzer.h"
-#include "interface/Co60SpectrumAnalyzer.h"
+#include "interface/Na22SpectrumAnalyzerSingleBar.h"
+#include "interface/Co60SpectrumAnalyzer_2Peaks.h"
#include "interface/FitUtils.h"
#include "interface/SetTDRStyle.h"
#include "CfgManager/interface/CfgManager.h"
@@ -63,6 +64,7 @@ int main(int argc, char** argv)
system(Form("mkdir -p %s/CTR/",plotDir.c_str()));
system(Form("mkdir -p %s/CTR_energyRatioCorr/",plotDir.c_str()));
+
//copia il file .php che si trova in una cartella fuori plotDir in plotDir definita nel config file
system(Form("cp %s/../index.php %s",plotDir.c_str(),plotDir.c_str()));
system(Form("cp %s/../index.php %s/qfine/",plotDir.c_str(),plotDir.c_str()));
@@ -73,6 +75,7 @@ int main(int argc, char** argv)
system(Form("cp %s/../index.php %s/CTR_energyRatioCorr/",plotDir.c_str(),plotDir.c_str()));
+
std::vector LRLabels;
LRLabels.push_back("L");
LRLabels.push_back("R");
@@ -167,10 +170,17 @@ int main(int argc, char** argv)
std::map outTrees2;
float energy511LR;
float energy1275LR;
+ float energy1786LR;
float energy511R;
float energy1275R;
+ float energy1786R;
float energy511L;
float energy1275L;
+ float energy1786L;
+ float entries511L;
+ float entries511R;
+ float entries1275L;
+ float entries1275R;
float theIndex;
float timeRes;
@@ -184,6 +194,7 @@ int main(int argc, char** argv)
std::map h1_deltaT_energyRatioCorr;
std::map*> rangesLR;
+ std::map*> rangesLR_new;
std::map > > peaksLR;
std::map > energyBinLR;
@@ -202,11 +213,13 @@ int main(int argc, char** argv)
float* vals = new float[6];
TLatex* latex;
TH1F* histo;
+ TH1F* histoL;
TProfile* prof;
TF1* func;
-
+ int x = 0;
+ int z = 0;
@@ -214,7 +227,11 @@ int main(int argc, char** argv)
//--- draw 1st plots
std::string source = opts.GetOpt("Input.sourceName");
std::string Na22 = "Na22";
+ std::string Na22SingleBar = "Na22SingleBar";
std::string Co60 = "Co60";
+ std::string Co60SumPeak = "Co60SumPeak";
+ std::string Laser = "Laser";
+
for(auto stepLabel : stepLabels)
{
float Vov = map_Vovs[stepLabel];
@@ -233,11 +250,19 @@ int main(int argc, char** argv)
outTrees[index] = new TTree(Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1),Form("data_bar%02dL-R_Vov%.1f_th%02.0f",iBar,Vov,vth1));
outTrees[index] -> Branch("energyPeak511LR",&energy511LR);
outTrees[index] -> Branch("energyPeak1275LR",&energy1275LR);
+ outTrees[index] -> Branch("energyPeak1786LR",&energy1786LR);
outTrees[index] -> Branch("energyPeak511L",&energy511L);
outTrees[index] -> Branch("energyPeak1275L",&energy1275L);
+ outTrees[index] -> Branch("energyPeak1786L",&energy1786L);
outTrees[index] -> Branch("energyPeak511R",&energy511R);
outTrees[index] -> Branch("energyPeak1275R",&energy1275R);
+ outTrees[index] -> Branch("energyPeak1786R",&energy1786R);
+ outTrees[index] -> Branch("entriesPeak511R",&entries511R);
+ outTrees[index] -> Branch("entriesPeak511L",&entries511L);
+ outTrees[index] -> Branch("entriesPeak511R",&entries1275R);
+ outTrees[index] -> Branch("entriesPeak511L",&entries1275L);
outTrees[index] -> Branch("indexID",&theIndex);
+
for(auto LRLabel : LRLabels )
@@ -300,201 +325,453 @@ int main(int argc, char** argv)
histo -> SetLineColor(kRed);
histo -> SetLineWidth(2);
histo -> Draw();
+
+ if( source.compare(Na22) && source.compare(Na22SingleBar) && source.compare(Co60) && source.compare(Co60SumPeak) && source.compare(Laser)){
+ std::cout << " Missspelled radioactive source " << std::endl;
+ return(0);
+ }
+
+ if( LRLabel.compare("L")){
+
+ if(!source.compare(Na22)){
+ rangesL[index] = new std::vector;
+ peaksL[index] = Na22SpectrumAnalyzer(histo,rangesL[index]);
+
+ if (peaksL[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
+
+ if (peaksL[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40);
+ break;
+ }
+ peaksL.erase(index);
+ rangesL.erase(index);
+ peaksL[index]["0.511 MeV"].first = -10;
+ peaksL[index]["1.275 MeV"].first = -10;
+
+ }
+
+
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksL[index]["0.511 MeV"].first != -10){
+ for(int i = 1; i < rangesL[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinL[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
+ }
- if( LRLabel.compare("L")){
- if(!source.compare(Na22)){
- rangesL[index] = new std::vector;
- peaksL[index] = Na22SpectrumAnalyzer(histo,rangesL[index]);
- if (peaksL[index]["0.511 MeV"].first== -9999){
- peaksL.erase(index);
- rangesL.erase(index);
- peaksL[index]["0.511 MeV"].first = -10;
- peaksL[index]["1.275 MeV"].first = -10;
+ if(!source.compare(Na22SingleBar)){
+ rangesL[index] = new std::vector;
+ peaksL[index] = Na22SpectrumAnalyzerSingleBar(histo,rangesL[index]);
- }
+ if (peaksL[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
+ if (peaksL[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40);
+ break;
+ }
+ peaksL.erase(index);
+ rangesL.erase(index);
+ peaksL[index]["0.511 MeV"].first = -10;
+ peaksL[index]["1.275 MeV"].first = -10;
+ peaksL[index]["1.786 MeV"].first = -10;
+
+ }
+
+
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksL[index]["0.511 MeV"].first != -10){
+ for(int i = 1; i < rangesL[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinL[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
+ }
- energy511L = peaksL[index]["0.511 MeV"].first;
- energy1275L = peaksL[index]["1.275 MeV"].first;
- theIndex = index;
- outTrees[index] -> Fill();
-
-
-
- for(auto peak : peaksL[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
- }
-
-
- histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
- if(peaksL[index]["0.511 MeV"].first != -10){
- for(int i = 1; i < rangesL[index]->size(); i++){
- TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
- int j = 1;
- for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
- binHisto -> SetBinContent( j, histo->GetBinContent(bin));
- j++;
+
+ if(!source.compare(Co60)){
+ rangesL[index] = new std::vector;
+ peaksL[index] = Co60SpectrumAnalyzer_2Peaks(histo,rangesL[index]);
+ if (peaksL[index]["1.173 MeV"].first > 0){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
}
- energyBinL[index][i] = binHisto -> GetMean();
- binHisto -> Delete();
- }
- }
-
- }
- if(!source.compare(Co60)){
- rangesL[index] = new std::vector;
- peaksL[index] = Co60SpectrumAnalyzer(histo,rangesL[index]);
-
+ if (peaksL[index]["1.173 MeV"].first== -9999){
+ for(auto peak : peaksL[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ break;
+ }
+ peaksL.erase(index);
+ rangesL.erase(index);
+ peaksL[index]["1.173 MeV"].first = -10;
+ peaksL[index]["1.332 MeV"].first = -10;
+
+ }
- if (peaksL[index]["1.173 MeV"].first== -9999){
- peaksL.erase(index);
- rangesL.erase(index);
- peaksL[index]["1.173 MeV"].first = -10;
- peaksL[index]["1.332 MeV"].first = -10;
-
- }
+
- energy511L = peaksL[index]["1.173 MeV"].first;
- energy1275L = peaksL[index]["1.332 MeV"].first;
- theIndex = index;
- outTrees[index] -> Fill();
-
-
-
- for(auto peak : peaksL[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
- }
-
-
- histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
- if(peaksL[index]["1.173 MeV"].first != -10){
- for(int i = 1; i < rangesL[index]->size(); i++){
- TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
- int j = 1;
- for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
- binHisto -> SetBinContent( j, histo->GetBinContent(bin));
- j++;
+
+
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksL[index]["1.173 MeV"].first != -10){
+ for(int i = 1; i < rangesL[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinL[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
}
- energyBinL[index][i] = binHisto -> GetMean();
- binHisto -> Delete();
+
}
- }
-
- }
+
+
+
+
+ if(!source.compare(Co60SumPeak) ){
+
+ rangesL[index] = new std::vector;
+ TF1* fitFunc = new TF1 ( Form("funcL_%d",index), "gaus(0)", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+
+ //fitFunc -> SetParameter (1, histo -> GetBinCenter(histo->GetMaximumBin()));
+ //fitFunc -> SetParameter (1, histo -> GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+
+ rangesL[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5);
+ rangesL[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ energy511L = fitFunc->GetParameter(1);
+
+ for(auto range: (*rangesL[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesL[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinL[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
+ if( !source.compare(Laser)){
+
+ rangesL[index] = new std::vector;
+ TF1* fitFunc = new TF1 ( Form("funcL_%d",index), "gaus", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
- if( source.compare(Na22) && source.compare(Co60)){
- std::cout << " Missspelled radioactive source " << std::endl;
- return(0);
- }
- }
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *10, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *10);
+
+ if ((fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5)>0) rangesL[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5);
+ else rangesL[index] -> push_back(0);
+ rangesL[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ energy511L = fitFunc->GetParameter(1);
+
+ for(auto range: (*rangesL[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesL[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesL[index]->at(i)) - histo->FindBin(rangesL[index]-> at(i-1)), rangesL[index]-> at(i-1), rangesL[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesL[index]->at(i-1)) ; bin < histo -> FindBin(rangesL[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinL[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
- if( LRLabel.compare("R")){
- if(!source.compare(Na22)){
- rangesR[index] = new std::vector;
- peaksR[index] = Na22SpectrumAnalyzer(histo,rangesR[index]);
- if (peaksR[index]["0.511 MeV"].first== -9999){
- peaksR.erase(index);
- rangesR.erase(index);
- peaksR[index]["0.511 MeV"].first = -10;
- peaksR[index]["1.275 MeV"].first = -10;
-
- }
- energy511R = peaksR[index]["0.511 MeV"].first;
- energy1275R = peaksR[index]["1.275 MeV"].first;
- theIndex = index;
- outTrees[index] -> Fill();
-
-
-
- for(auto peak : peaksR[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
}
+ if( LRLabel.compare("R")){
+ if(!source.compare(Na22)){
+ rangesR[index] = new std::vector;
+ peaksR[index] = Na22SpectrumAnalyzer(histo,rangesR[index]);
- histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
- if(peaksR[index]["0.511 MeV"].first != -10){
- for(int i = 1; i < rangesR[index]->size(); i++){
- TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
- int j = 1;
- for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
- binHisto -> SetBinContent( j, histo->GetBinContent(bin));
- j++;
+ if (peaksR[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
+
+ if (peaksR[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ break;
+ }
+ peaksR.erase(index);
+ rangesR.erase(index);
+ peaksR[index]["0.511 MeV"].first = -10;
+ peaksR[index]["1.275 MeV"].first = -10;
+
+ }
+
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksR[index]["0.511 MeV"].first != -10){
+ for(int i = 1; i < rangesR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinR[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
}
- energyBinR[index][i] = binHisto -> GetMean();
- binHisto -> Delete();
- }
- }
-
- }
- if(!source.compare(Co60)){
- rangesR[index] = new std::vector;
- peaksR[index] = Co60SpectrumAnalyzer(histo,rangesR[index]);
+ if(!source.compare(Na22SingleBar)){
+ rangesR[index] = new std::vector;
+ peaksR[index] = Na22SpectrumAnalyzerSingleBar(histo,rangesR[index]);
+ if (peaksR[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
- if (peaksR[index]["1.173 MeV"].first== -9999){
- peaksR.erase(index);
- rangesR.erase(index);
- peaksR[index]["1.173 MeV"].first = -10;
- peaksR[index]["1.332 MeV"].first = -10;
+ if (peaksR[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ break;
+ }
+ peaksR.erase(index);
+ rangesR.erase(index);
+ peaksR[index]["0.511 MeV"].first = -10;
+ peaksR[index]["1.275 MeV"].first = -10;
+ peaksR[index]["1.786 MeV"].first = -10;
+
+ }
+
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksR[index]["0.511 MeV"].first != -10){
+ for(int i = 1; i < rangesR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinR[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
- }
+ }
- energy511R = peaksR[index]["1.173 MeV"].first;
- energy1275R = peaksR[index]["1.332 MeV"].first;
- theIndex = index;
- outTrees[index] -> Fill();
-
-
-
- for(auto peak : peaksR[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
- }
-
-
- histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
- if(peaksR[index]["1.173 MeV"].first != -10){
- for(int i = 1; i < rangesR[index]->size(); i++){
- TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
- int j = 1;
- for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
- binHisto -> SetBinContent( j, histo->GetBinContent(bin));
- j++;
+ if(!source.compare(Co60)){
+ rangesR[index] = new std::vector;
+ peaksR[index] = Co60SpectrumAnalyzer_2Peaks(histo,rangesR[index]);
+
+ if (peaksR[index]["1.173 MeV"].first > 0){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
+
+ if (peaksR[index]["1.173 MeV"].first== -9999){
+ for(auto peak : peaksR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ break;
+ }
+ peaksR.erase(index);
+ rangesR.erase(index);
+ peaksR[index]["1.173 MeV"].first = -10;
+ peaksR[index]["1.332 MeV"].first = -10;
+
+ }
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksR[index]["1.173 MeV"].first != -10){
+ for(int i = 1; i < rangesR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinR[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
}
- energyBinR[index][i] = binHisto -> GetMean();
- binHisto -> Delete();
+
}
- }
-
- }
+
+ if(!source.compare(Co60SumPeak)){
+ rangesR[index] = new std::vector;
-
- if( source.compare(Na22) && source.compare(Co60)){
- std::cout << " Missspelled radioactive source " << std::endl;
- return(0);
- }
- }
+ TF1* fitFunc = new TF1 ( Form("funcR_%d",index), "gaus(0)", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ rangesR[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(0) *5);
+ rangesR[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(0) *5);
+ energy511R = fitFunc->GetParameter(1);
+
+ for(auto range: (*rangesR[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinR[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
+
+
+ if(!source.compare(Laser)){
+ rangesR[index] = new std::vector;
+
+ TF1* fitFunc = new TF1 ( Form("funcR_%d",index), "gaus(0)", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *10, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *10);
+ if ((fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5)>0) rangesR[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5);
+ else rangesR[index] -> push_back(0);
+ rangesR[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ energy511R = fitFunc->GetParameter(1);
+
+ for(auto range: (*rangesR[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesR[index]->at(i)) - histo->FindBin(rangesR[index]-> at(i-1)), rangesR[index]-> at(i-1), rangesR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesR[index]->at(i-1)) ; bin < histo -> FindBin(rangesR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+ }
+ energyBinR[index][i] = binHisto -> GetMean();
+ binHisto -> Delete();
+ }
+ }
+
+
+
+
+
+ }
+
+
+
+
- latex -> Draw("same");
- histo -> Write();
+ latex -> Draw("same");
+ outFile -> cd();
+ histo->Write();
c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str()));
c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str()));
delete c;
@@ -527,9 +804,24 @@ int main(int argc, char** argv)
if(!source.compare(Na22)){
rangesLR[index] = new std::vector;
+ rangesLR_new[index] = new std::vector;
peaksLR[index] = Na22SpectrumAnalyzer(histo,rangesLR[index]);
+
+ if (peaksLR[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
if (peaksLR[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ std::cout<<"bar"<Integral(histo->FindBin(energy511LR-0.1*energy511LR),histo->FindBin(energy511LR));
+ entries511R = histo->Integral(histo->FindBin(energy511LR)+1, histo->FindBin(energy511LR+0.1*energy511LR));
+
+
+
+
+
+
outTrees[index] -> Fill();
-
- for(auto peak : peaksLR[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
+
+ histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
+ if(peaksLR[index]["0.511 MeV"].first != -10){
+ for(int i = 1; i < rangesLR[index]->size(); i++){
+ TH1F *binHisto = new TH1F ( "binHisto", "binHisto", histo -> FindBin(rangesLR[index]->at(i)) - histo->FindBin(rangesLR[index]-> at(i-1)), rangesLR[index]-> at(i-1), rangesLR[index]->at(i));
+ int j = 1;
+ for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){
+ binHisto -> SetBinContent( j, histo->GetBinContent(bin));
+ j++;
+
+ }
+ energyBinLR[index][i] = binHisto -> GetMean();
+
+ binHisto -> Delete();
+ }
}
+
+ }
+
+
+ gStyle->SetOptFit(1111);
+
+ if(!source.compare(Na22SingleBar)){
+ rangesLR[index] = new std::vector;
+ rangesLR_new[index] = new std::vector;
+ peaksLR[index] = Na22SpectrumAnalyzerSingleBar(histo,rangesLR[index]);
+
+ if (peaksLR[index]["0.511 MeV"].first > 0){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
+
+ if (peaksLR[index]["0.511 MeV"].first== -9999){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40.);
+ std::cout<<"bar"<Integral(histo->FindBin(energy511LR-0.1*energy511LR),histo->FindBin(energy511LR));
+ entries511R = histo->Integral(histo->FindBin(energy511LR)+1, histo->FindBin(energy511LR+0.1*energy511LR));
+
+ outTrees[index] -> Fill();
+
histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
@@ -560,42 +938,71 @@ int main(int argc, char** argv)
for (int bin = histo->FindBin(rangesLR[index]->at(i-1)) ; bin < histo -> FindBin(rangesLR[index]->at(i)) +1 ; bin++){
binHisto -> SetBinContent( j, histo->GetBinContent(bin));
j++;
+
}
energyBinLR[index][i] = binHisto -> GetMean();
+
binHisto -> Delete();
}
}
}
- if(!source.compare(Co60)){
+
+
+
+
+
+
+
+
+
+
+
+ if(!source.compare(Co60)){
+
rangesLR[index] = new std::vector;
- peaksLR[index] = Co60SpectrumAnalyzer(histo,rangesLR[index]);
+ peaksLR[index] = Co60SpectrumAnalyzer_2Peaks(histo,rangesLR[index]);
+ if (peaksLR[index]["1.173 MeV"].first > 0){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
+ break;
+ }
+ }
if (peaksLR[index]["1.173 MeV"].first== -9999){
+ for(auto peak : peaksLR[index] )
+ {
+ histo -> GetXaxis() -> SetRangeUser(0.,40);
+ break;
+ }
peaksLR.erase(index);
rangesLR.erase(index);
peaksLR[index]["1.173 MeV"].first = -10;
peaksLR[index]["1.332 MeV"].first = -10;
+ peaksLR[index]["2.505 MeV"].first = -10;
}
+
+
+
+
+ energy511R = peaksR[index]["1.173 MeV"].first;
+ energy1275R = peaksR[index]["1.332 MeV"].first;
+ energy511L = peaksL[index]["1.173 MeV"].first;
+ energy1275L = peaksL[index]["1.332 MeV"].first;
energy511LR = peaksLR[index]["1.173 MeV"].first;
energy1275LR = peaksLR[index]["1.332 MeV"].first;
+ energy1786LR = peaksLR[index]["2.505 MeV"].first;
theIndex = index;
outTrees[index] -> Fill();
-
-
- for(auto peak : peaksLR[index] )
- {
- histo -> GetXaxis() -> SetRangeUser(0.,5.*peak.second.first);
- break;
- }
-
+
histo -> GetYaxis() -> SetRangeUser(3.,5.*histo->GetMaximum());
if(peaksLR[index]["1.173 MeV"].first != -10){
for(int i = 1; i < rangesLR[index]->size(); i++){
@@ -610,21 +1017,88 @@ int main(int argc, char** argv)
}
}
- }
-
-
- if( source.compare(Na22) && source.compare(Co60)){
- std::cout << " Missspelled radioactive source " << std::endl;
- return(0);
}
-
-
+
+ if(!source.compare(Co60SumPeak)){
+ rangesLR[index] = new std::vector;
+
+ TF1* fitFunc = new TF1 ( Form("funcLR_%d",index), "gaus(0)", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ rangesLR[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5);
+ rangesLR[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ energy511LR = fitFunc->GetParameter(1);
+
+
+ energy1275LR = -10;
+ energy1275L = -10;
+ energy1275R = -10;
+ theIndex = index;
+ outTrees[index] -> Fill();
+
+ for(auto range: (*rangesLR[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesLR[index]->size(); i++){
+ energyBinLR[index][i] = energy511LR;
+ }
+ }
+
+
+
+
+
+ if(!source.compare(Laser)){
+ rangesLR[index] = new std::vector;
+
+ TF1* fitFunc = new TF1 ( Form("funcLR_%d",index), "gaus", histo -> GetBinCenter(histo->GetMaximumBin()) - 2*histo->GetRMS(), histo -> GetBinCenter(histo->GetMaximumBin()) + 2*histo->GetRMS());
+ histo -> Fit( fitFunc, "NQR");
+ fitFunc -> SetParameter(0, fitFunc->GetParameter(0));
+ fitFunc -> SetParameter(1, fitFunc->GetParameter(1));
+ fitFunc -> SetParameter(2, fitFunc->GetParameter(2));
+ fitFunc -> SetLineColor(kBlack);
+ fitFunc -> SetLineWidth(2);
+ histo -> Fit ( fitFunc, "QRS+", "", fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *10, fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *10);
+ if ((fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5)>0) rangesLR[index] -> push_back( fitFunc->GetParameter(1) - fitFunc->GetParameter(2) *5);
+ else rangesLR[index] -> push_back(0);
+ rangesLR[index] -> push_back( fitFunc->GetParameter(1) + fitFunc->GetParameter(2) *5);
+ energy511LR = fitFunc->GetParameter(1);
+
+
+ energy1275LR = -10;
+ energy1275L = -10;
+ energy1275R = -10;
+ theIndex = index;
+ outTrees[index] -> Fill();
+
+ for(auto range: (*rangesLR[index])){
+ float yval = std::max(10., histo->GetBinContent(histo->FindBin(range)));
+ TLine* line = new TLine(range,3.,range, yval);
+ line -> SetLineWidth(1);
+ line -> SetLineStyle(7);
+ line -> Draw("same");
+ }
+
+ for(int i = 1; i < rangesLR[index]->size(); i++){
+ energyBinLR[index][i] = energy511LR;
+ }
+ }
-
latex -> Draw("same");
- histo -> Write();
+ outFile -> cd();
+ histo->Write();
c -> Print(Form("%s/energy/c_energy__%s.png",plotDir.c_str(),label.c_str()));
c -> Print(Form("%s/energy/c_energy__%s.pdf",plotDir.c_str(),label.c_str()));
delete c;
@@ -639,9 +1113,7 @@ int main(int argc, char** argv)
// end 1st plots
-
-
-
+
//------------------------
@@ -663,8 +1135,8 @@ int main(int argc, char** argv)
accept[index1][entry] = false;
-
- if( !rangesLR[index1] ) continue;
+ if(!rangesLR[index1] ) continue;
+
int energyBinAverage = FindBin(0.5*(anEvent->energyL+anEvent->energyR),rangesLR[index1])+1;
if( energyBinAverage < 1 ) continue;
@@ -673,21 +1145,27 @@ int main(int argc, char** argv)
accept[index1][entry] = true;
-
+
if( h1_energyRatio[index2] == NULL )
{
- std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
-
- h1_energyRatio[index2] = new TH1F(Form("h1_energyRatio_%s",labelLR_energyBin.c_str()),"",1000,0.,5.);
+ std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%02d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
+
+
+
+ h1_energyRatio[index2] = new TH1F(Form("h1_energyRatio_%s",labelLR_energyBin.c_str()),"",1000,0.,5.);
+
+
h1_deltaT_raw[index2] = new TH1F(Form("h1_deltaT_raw_%s",labelLR_energyBin.c_str()),"",1250,-5000.,5000.);
}
-
- h1_energyRatio[index2] -> Fill( anEvent->energyR / anEvent->energyL );
- h1_deltaT_raw[index2] -> Fill( anEvent->timeR-anEvent->timeL );
+
+ if ((anEvent->energyR / anEvent->energyL >0) && (anEvent->energyR / anEvent->energyL <5)){
+ h1_energyRatio[index2] -> Fill( anEvent->energyR / anEvent->energyL );
+ h1_deltaT_raw[index2] -> Fill( anEvent->timeR-anEvent->timeL );
+ }
}
std::cout << std::endl;
}
-
+
//------------------
//--- draw 2nd plots
@@ -695,6 +1173,7 @@ int main(int argc, char** argv)
std::map CTRSigmas;
std::map fitFunc_energyRatio;
+ std::map fitFunc_energyRatio2;
for(auto mapIt : h1_deltaT_raw)
{
@@ -721,23 +1200,26 @@ int main(int argc, char** argv)
int index1( (10000*int(Vov*100.)) + (100*vth1) + iBar );
std::string labelLR(Form("bar%02dL-R_%s",iBar,stepLabel.c_str()));
-
if( !rangesLR[index1] ) continue;
int nEnergyBins = rangesLR[index1]->size()-1;
for(int iEnergyBin = 1; iEnergyBin <= nEnergyBins; ++iEnergyBin)
{
+ //if (rangesLR[index1]->at(iEnergyBin)<0) continue;
double index2( 10000000*iEnergyBin+index1 );
+
+
if (!h1_energyRatio[index2]) continue;
+
+ std::string labelLR_energyBin(Form("%s_energyBin%02d",labelLR.c_str(),iEnergyBin));
+
- std::string labelLR_energyBin(Form("%s_energyBin%d",labelLR.c_str(),iEnergyBin));
-
- c = new TCanvas(Form("c_energyRatio_%s",labelLR_energyBin.c_str()),Form("c_energyRatio_%s",labelLR_energyBin.c_str()));
+ c = new TCanvas(Form("c_energyRatio_%s",labelLR_energyBin.c_str()),Form("c_energyRatio_%s",labelLR_energyBin.c_str()));
histo = h1_energyRatio[index2];
histo -> GetXaxis() -> SetRangeUser(histo->GetMean()-5.*histo->GetRMS(),histo->GetMean()+5.*histo->GetRMS());
- histo -> SetMaximum(1.25*histo->GetBinContent(histo->FindBin(FindXMaximum(histo,histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()))));
+ histo -> SetMaximum(1.25*histo->GetBinContent(histo->FindBin(FindXMaximum(histo,histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS()))));
histo -> SetTitle(Form(";energy_{right} / energy_{left};entries"));
histo -> SetLineColor(kRed);
histo -> SetLineWidth(2);
@@ -745,12 +1227,32 @@ int main(int argc, char** argv)
histo -> Write();
- fitFunc_energyRatio[index2] = new TF1(Form("fitFunc_energyRatio_%s",labelLR_energyBin.c_str()),"gaus",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS());
+ /*fitFunc_energyRatio[index2] = new TF1(Form("fitFunc_energyRatio_%s",labelLR_energyBin.c_str()),"gaus",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS());
histo -> Fit(fitFunc_energyRatio[index2],"QNRS");
histo -> Fit(fitFunc_energyRatio[index2],"QS+","",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS());
fitFunc_energyRatio[index2] -> SetLineColor(kBlack);
fitFunc_energyRatio[index2] -> SetLineWidth(2);
- fitFunc_energyRatio[index2] -> Draw("same");
+ fitFunc_energyRatio[index2] -> Draw("same");*/
+
+
+
+ fitFunc_energyRatio[index2] = new TF1(Form("fitFunc_energyRatio_%s",labelLR_energyBin.c_str()),"gaus",histo->GetMean()-2.*histo->GetRMS(),histo->GetMean()+2.*histo->GetRMS());
+ histo -> Fit(fitFunc_energyRatio[index2],"QNRS");
+ histo -> Fit(fitFunc_energyRatio[index2],"QS+","",fitFunc_energyRatio[index2]->GetParameter(1)-2.*fitFunc_energyRatio[index2]->GetParameter(2),fitFunc_energyRatio[index2]->GetParameter(1)+2.*fitFunc_energyRatio[index2]->GetParameter(2));
+
+
+ double min =fitFunc_energyRatio[index2]->GetParameter(1)-2.*fitFunc_energyRatio[index2]->GetParameter(2);
+ double max =fitFunc_energyRatio[index2]->GetParameter(1)+2.*fitFunc_energyRatio[index2]->GetParameter(2);
+ fitFunc_energyRatio2[index2] = new TF1(Form("fitFunc_energyRatio2_%s",labelLR_energyBin.c_str()),"gaus",min,max);
+ fitFunc_energyRatio2[index2]->SetParameter(0, fitFunc_energyRatio[index2]->GetParameter(0));
+ fitFunc_energyRatio2[index2]->SetParameter(1, fitFunc_energyRatio[index2]->GetParameter(1));
+ fitFunc_energyRatio2[index2]->SetParameter(2, fitFunc_energyRatio[index2]->GetParameter(2));
+ histo -> Fit(fitFunc_energyRatio2[index2],"NQRS");
+
+
+ fitFunc_energyRatio2[index2] -> SetLineColor(kBlack);
+ fitFunc_energyRatio2[index2] -> SetLineWidth(2);
+ fitFunc_energyRatio2[index2] -> Draw("same");
@@ -771,7 +1273,7 @@ int main(int argc, char** argv)
}
-
+
//------------------------
@@ -797,21 +1299,31 @@ int main(int argc, char** argv)
- float energyRatioMean = fitFunc_energyRatio[index2]->GetParameter(1);
- float energyRatioSigma = fitFunc_energyRatio[index2]->GetParameter(2);
+ float energyRatioMean = fitFunc_energyRatio2[index2]->GetParameter(1);
+ float energyRatioSigma = fitFunc_energyRatio2[index2]->GetParameter(2);
+
+
if( fabs(anEvent->energyR/anEvent->energyL-energyRatioMean) > 2.*energyRatioSigma )
{
accept[index1][entry] = false;
continue;
}
+
+
+ if( (anEvent->energyR/anEvent->energyL)>5 || (anEvent->energyR/anEvent->energyL)<0)
+ {
+ accept[index1][entry] = false;
+ continue;
+ }
+
long long deltaT = anEvent->timeR - anEvent->timeL;
if( h1_deltaT[index2] == NULL )
{
- std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
+ std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%02d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
h1_deltaT[index2] = new TH1F(Form("h1_deltaT_%s",labelLR_energyBin.c_str()),"",1000,-5000,5000.);
@@ -846,15 +1358,16 @@ int main(int argc, char** argv)
std::string labelLR(Form("bar%02dL-R_%s",iBar,stepLabel.c_str()));
if( !rangesLR[index1] ) continue;
+
int nEnergyBins = rangesLR[index1]->size()-1;
for(int iEnergyBin = 1; iEnergyBin <= nEnergyBins; ++iEnergyBin)
{
+ //if (rangesLR[index1]->at(iEnergyBin)<0) continue;
double index2( 10000000*iEnergyBin+index1 );
if(!p1_deltaT_vs_energyRatio[index2]) continue;
-
- std::string labelLR_energyBin(Form("%s_energyBin%d",labelLR.c_str(),iEnergyBin));
+ std::string labelLR_energyBin(Form("%s_energyBin%02d",labelLR.c_str(),iEnergyBin));
c = new TCanvas(Form("c_deltaT_vs_energyRatio_%s",labelLR_energyBin.c_str()),Form("c_deltaT_vs_energyRatio_%s",labelLR_energyBin.c_str()));
@@ -871,7 +1384,7 @@ int main(int argc, char** argv)
latex -> SetTextColor(kRed);
latex -> Draw("same");
- func = fitFunc_energyRatio[index2];
+ func = fitFunc_energyRatio2[index2];
float fitXMin = func->GetParameter(1) - 2.*func->GetParameter(2);
float fitXMax = func->GetParameter(1) + 2.*func->GetParameter(2);
@@ -892,6 +1405,9 @@ int main(int argc, char** argv)
//------------------------
//--- 4th loop over events
+
+ gStyle->SetOptFit(1111);
+
for(auto mapIt : trees)
{
ModuleEventClass* anEvent = new ModuleEventClass();
@@ -905,10 +1421,10 @@ int main(int argc, char** argv)
mapIt.second -> GetEntry(entry);
int index1( (10000*int(anEvent->Vov*100.)) + (100*anEvent->vth1) + anEvent->barID );
+
if( !accept[index1][entry] ) continue;
-
int energyBinAverage = FindBin(0.5*(anEvent->energyL+anEvent->energyR),rangesLR[index1])+1;
double index2( 10000000*energyBinAverage+index1 );
@@ -920,15 +1436,17 @@ int main(int argc, char** argv)
energyRatioCorr = fitFunc_energyRatioCorr[index2]->Eval(anEvent->energyR/anEvent->energyL) -
- fitFunc_energyRatioCorr[index2]->Eval(fitFunc_energyRatio[index2]->GetParameter(1));
+ fitFunc_energyRatioCorr[index2]->Eval(fitFunc_energyRatio2[index2]->GetParameter(1));
if( h1_deltaT_energyRatioCorr[index2] == NULL )
{
- std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
+ std::string labelLR_energyBin(Form("bar%02dL-R_Vov%.1f_th%02d_energyBin%02d",anEvent->barID,anEvent->Vov,anEvent->vth1,energyBinAverage));
+
h1_deltaT_energyRatioCorr[index2] = new TH1F(Form("h1_deltaT_energyRatioCorr_%s",labelLR_energyBin.c_str()),"",1000,-5000.,5000.);
+
}
h1_deltaT_energyRatioCorr[index2] -> Fill( deltaT - energyRatioCorr );
@@ -940,6 +1458,11 @@ int main(int argc, char** argv)
double theIndex2;
float enBin;
+ float theDeltaTMean;
+ float theEntriesCTR;
+ float errTimeRes;
+ float timeResEffSigma;
+
//------------------
//--- draw 4th plots
@@ -961,18 +1484,28 @@ int main(int argc, char** argv)
for(int iEnergyBin = 1; iEnergyBin <= nEnergyBins; ++iEnergyBin)
{
- double index2( 10000000*iEnergyBin + index1 );
+
+ double index2( 10000000*iEnergyBin + index1 );
+ double Vov_iBar_enBin_ID =10000000*iEnergyBin + 10000*int(Vov*100.) + iBar;
outTrees2[index2] = new TTree(Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin),Form("dataRes_bar%02dL-R_Vov%.1f_th%02.0f_enBin%02d",iBar,Vov,vth1,iEnergyBin));
- outTrees2[index2] -> Branch("energyBin", &enBin);
+ outTrees2[index2] -> Branch("energyBin", &enBin);
outTrees2[index2] -> Branch("timeResolution",&timeRes);
+ outTrees2[index2] -> Branch("effSigma",&timeResEffSigma);
+ outTrees2[index2] -> Branch("errTimeResolution",&errTimeRes);
outTrees2[index2] -> Branch("indexID2",&theIndex2);
+ outTrees2[index2] -> Branch("DeltaTMean",&theDeltaTMean);
+ outTrees2[index2] -> Branch("EntriesCTR",&theEntriesCTR);
+
+
+
+
if(!h1_deltaT_energyRatioCorr[index2]) continue;
- std::string labelLR_energyBin(Form("%s_energyBin%d",labelLR.c_str(),iEnergyBin));
+ std::string labelLR_energyBin(Form("%s_energyBin%02d",labelLR.c_str(),iEnergyBin));
c = new TCanvas(Form("c_deltaT_energyRatioCorr_%s",labelLR_energyBin.c_str()),Form("c_deltaT_energyRatioCorr_%s",labelLR_energyBin.c_str()));
@@ -984,48 +1517,127 @@ int main(int argc, char** argv)
histo -> SetLineWidth(2);
histo -> SetLineColor(kBlue);
histo -> SetMarkerColor(kBlue);
- histo -> RebinX(2);
+ histo -> RebinX(2);
+
+
histo -> Draw("");
+
+
+ FindSmallestInterval(vals,histo,0.68);
+ float min = vals[4];
+ float max = vals[5];
+ float delta = max-min;
+ float sigma = 0.5*delta;
+ float effSigma = sigma;
+
float fitXMin = histo->GetBinCenter(histo->GetMaximumBin()) - 200.;
float fitXMax = histo->GetBinCenter(histo->GetMaximumBin()) + 200.;
+
+
+
TF1* fitFunc = new TF1(Form("fitFunc_energyCorr_%s",labelLR_energyBin.c_str()),"gaus",fitXMin,fitXMax);
fitFunc -> SetParameters(1,histo->GetMean(),histo->GetRMS());
histo -> Fit(fitFunc,"QNRSL+","");
histo -> Fit(fitFunc,"QNRSL+","",fitFunc->GetParameter(1)-fitFunc->GetParameter(2),fitFunc->GetParameter(1)+fitFunc->GetParameter(2));
-
- float fitXMin2 = fitFunc->GetParameter(1)-1.5*fitFunc->GetParameter(2);
- float fitXMax2 = fitFunc->GetParameter(1)+1.5*fitFunc->GetParameter(2);
+
+ float fitXMin2 = fitFunc->GetParameter(1)-1.5*fitFunc->GetParameter(2);
+ float fitXMax2 = fitFunc->GetParameter(1)+1.5*fitFunc->GetParameter(2);
+ //if(!SiPM_type.compare(HDR2) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3) && !source.compare(Na22)) fitXMin2 = fitFunc->GetParameter(1)-1.1*fitFunc->GetParameter(2);
+ //if(!source.compare(Na22SingleBar) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3)) fitXMin2 = fitFunc->GetParameter(1)-1.*fitFunc->GetParameter(2);
+ //if(!source.compare(Na22SingleBar) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3)) fitXMax2 = fitFunc->GetParameter(1)+1.1*fitFunc->GetParameter(2);
TF1* fitFunc2 = new TF1(Form("fitFunc2_energyCorr_%s",labelLR_energyBin.c_str()),"gaus(0)",fitXMin2,fitXMax2);
fitFunc2 -> SetParameter(0,fitFunc->GetParameter(0));
fitFunc2 -> SetParameter(1,fitFunc->GetParameter(1));
fitFunc2 -> SetParameter(2,fitFunc->GetParameter(2));
histo -> Fit(fitFunc2,"QRSL+");
+ histo -> Fit(fitFunc2,"QRSL+");
- fitFunc2 -> SetLineColor(kBlue+1);
- fitFunc2 -> SetLineWidth(3);
- fitFunc2 -> Draw("same");
-
- FindSmallestInterval(vals,histo,0.68);
- float min = vals[4];
- float max = vals[5];
- float delta = max-min;
- float sigma = 0.5*delta;
- float effSigma = sigma;
+ /*if (effSigma>1.1*fabs(fitFunc2->GetParameter(2)) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3) && !source.compare(Na22)){
+ FindSmallestInterval(vals,histo,0.68);
+ float min = vals[4];
+ float max = vals[5];
+ float delta = max-min;
+ float sigma = 0.5*delta;
+ float effSigma = sigma;
+
+ float fitXMin = histo->GetBinCenter(histo->GetMaximumBin()) - 2000.;
+ float fitXMax = histo->GetBinCenter(histo->GetMaximumBin()) + 2000.;
+
+ if (!SiPM_type.compare(HDR2) && histo->GetBinCenter(histo->GetMaximumBin())<-3000){
+ fitXMin = 0. -2000.;
+ fitXMax = 0. +2000.;
+ }
+
+ fitFunc = new TF1(Form("fitFunc_energyCorr_%s",labelLR_energyBin.c_str()),"gaus",fitXMin,fitXMax);
+ fitFunc -> SetParameters(1,histo->GetMean(),histo->GetRMS());
+ histo -> Fit(fitFunc,"QNRSL+","");
+ histo -> Fit(fitFunc,"QNRSL+","",fitFunc->GetParameter(1)-fitFunc->GetParameter(2),fitFunc->GetParameter(1)+fitFunc->GetParameter(2));
+
+ float fitXMin2 = fitFunc->GetParameter(1)-1.5*fitFunc->GetParameter(2);
+ if(!SiPM_type.compare(HDR2)) fitXMin2 = fitFunc->GetParameter(1)-1.1*fitFunc->GetParameter(2);
+ //if(th==40 || th==50 || th==60) fitXMin2 = fitFunc->GetParameter(1)-1.*fitFunc->GetParameter(2);
+ float fitXMax2 = fitFunc->GetParameter(1)+1.5*fitFunc->GetParameter(2);
+ fitFunc2 = new TF1(Form("fitFunc2_energyCorr_%s",labelLR_energyBin.c_str()),"gaus(0)",fitXMin2,fitXMax2);
+ fitFunc2 -> SetParameter(0,fitFunc->GetParameter(0));
+ fitFunc2 -> SetParameter(1,fitFunc->GetParameter(1));
+ fitFunc2 -> SetParameter(2,fitFunc->GetParameter(2));
+ histo -> Fit(fitFunc2,"QRSL+");
+ histo -> Fit(fitFunc2,"QRSL+");
+ x++;
+
+ }*/
+
+
+ fitFunc2 -> SetLineColor(kGreen+1);
+ fitFunc2 -> SetLineWidth(3);
+ fitFunc2 -> Draw("same");
+
+
+ outFile -> cd();
+ histo -> Write();
+
+
+
+
enBin = energyBinLR[index1][iEnergyBin];
theIndex2 = index2;
timeRes = fabs(fitFunc2->GetParameter(2));
+ timeResEffSigma = effSigma;
+
+ theDeltaTMean = fitFunc2->GetParameter(1);
+ theEntriesCTR = fitFunc2->GetParameter(0);
+ errTimeRes = fabs(fitFunc2->GetParError(2));
+
+ int contr = 0;
+ if(!source.compare("Na22") && fitFunc2->GetParameter(0)<20){
+ timeRes = 0;
+ effSigma = 0;
+ contr = 1;
+ }
+
+ /*if(!source.compare("Na22SingleBar") && iEnergyBin==1 && fitFunc2->GetParameter(0)<100){
+ timeRes = 0;
+ effSigma = 0;
+ contr = 1;
+ }*/
+
+
outTrees2[index2]->Fill();
- histo -> GetXaxis() -> SetRangeUser(fitFunc2->GetParameter(1)-10.*fitFunc2->GetParameter(2),
+ if (!source.compare("Laser")) histo -> GetXaxis() -> SetRangeUser(fitFunc2->GetParameter(1)-10.*fitFunc2->GetParameter(2),
fitFunc2->GetParameter(1)+10.*fitFunc2->GetParameter(2));
+ histo -> SetMaximum(histo->GetMaximum()+0.1*histo->GetMaximum());
+
+
latex = new TLatex(0.55,0.85,Form("#splitline{#sigma_{corr.}^{eff} = %.0f ps}{#sigma_{corr.}^{gaus} = %.0f ps}",effSigma,fabs(fitFunc2->GetParameter(2))));
+ if(contr==1) latex = new TLatex(0.55,0.85,Form("#splitline{#sigma_{corr.}^{eff} = %.0f ps}{#sigma_{corr.}^{gaus} = %.0f ps}",effSigma,timeRes));
latex -> SetNDC();
latex -> SetTextFont(42);
latex -> SetTextSize(0.04);
@@ -1033,40 +1645,90 @@ int main(int argc, char** argv)
latex -> Draw("same");
// -- raw delta T
- c->cd();
+ c->cd();
histo = h1_deltaT[index2];
histo -> SetLineWidth(2);
histo -> SetLineColor(kRed);
histo -> SetMarkerColor(kRed);
histo -> RebinX(2);
+
histo -> Draw("same");
+
+
+ FindSmallestInterval(vals,histo,0.68);
+ min = vals[4];
+ max = vals[5];
+ delta = max-min;
+ sigma = 0.5*delta;
+ effSigma = sigma;
fitXMin = histo->GetBinCenter(histo->GetMaximumBin()) - 200.;
fitXMax = histo->GetBinCenter(histo->GetMaximumBin()) + 200.;
+
+ /*if (!SiPM_type.compare(HDR2) && histo->GetBinCenter(histo->GetMaximumBin())<-3000){
+ fitXMin = 0. -200.;
+ fitXMax = 0. +200.;
+ }*/
+
fitFunc = new TF1(Form("fitFunc_%s",labelLR_energyBin.c_str()),"gaus",fitXMin,fitXMax);
fitFunc -> SetParameters(1,histo->GetMean(),histo->GetRMS());
histo -> Fit(fitFunc,"QNRSL+","");
histo -> Fit(fitFunc,"QNRSL+","",fitFunc->GetParameter(1)-fitFunc->GetParameter(2),fitFunc->GetParameter(1)+fitFunc->GetParameter(2));
fitXMin2 = fitFunc->GetParameter(1)-fabs(1.5*fitFunc->GetParameter(2));
- fitXMax2 = fitFunc->GetParameter(1)+fabs(1.5*fitFunc->GetParameter(2));
+ fitXMax2 = fitFunc->GetParameter(1)+fabs(1.5*fitFunc->GetParameter(2));
+ //if(!SiPM_type.compare(HDR2) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3) && !source.compare(Na22)) fitXMin2 = fitFunc->GetParameter(1)-fabs(1.3*fitFunc->GetParameter(2));
+ //if((!source.compare(Na22SingleBar)) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3)) fitXMin2 = fitFunc->GetParameter(1)-1.*fitFunc->GetParameter(2);
+ //if(!source.compare(Na22SingleBar) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3)) fitXMax2 = fitFunc->GetParameter(1)+1.1*fitFunc->GetParameter(2);
fitFunc2 = new TF1(Form("fitFunc2_energyCorr_%s",labelLR_energyBin.c_str()),"gaus(0)",fitXMin2,fitXMax2);
fitFunc2 -> SetParameter(0,fitFunc->GetParameter(0));
fitFunc2 -> SetParameter(1,fitFunc->GetParameter(1));
fitFunc2 -> SetParameter(2,fitFunc->GetParameter(2));
histo -> Fit(fitFunc2,"QRSL+");
histo -> Fit(fitFunc2,"QRSL+");
-
- fitFunc2 -> SetLineColor(kRed+1);
- fitFunc2 -> SetLineWidth(3);
- fitFunc2 -> Draw("same");
-
- FindSmallestInterval(vals,histo,0.68);
+
+
+ /*if (effSigma>1.1*fabs(fitFunc2->GetParameter(2)) && (iEnergyBin==1 || iEnergyBin==2 || iEnergyBin==3) && !source.compare(Na22)){
+
+
+ FindSmallestInterval(vals,histo,0.68);
min = vals[4];
max = vals[5];
delta = max-min;
sigma = 0.5*delta;
effSigma = sigma;
+
+
+ fitXMin = histo->GetBinCenter(histo->GetMaximumBin()) - 2000.;
+ fitXMax = histo->GetBinCenter(histo->GetMaximumBin()) + 2000.;
+
+ if (!SiPM_type.compare(HDR2) && histo->GetBinCenter(histo->GetMaximumBin())<-3000){
+ fitXMin = 0. -2000.;
+ fitXMax = 0. +2000.;
+ }
+
+ fitFunc = new TF1(Form("fitFunc_%s",labelLR_energyBin.c_str()),"gaus",fitXMin,fitXMax);
+ fitFunc -> SetParameters(1,histo->GetMean(),histo->GetRMS());
+ histo -> Fit(fitFunc,"QNRSL+","");
+ histo -> Fit(fitFunc,"QNRSL+","",fitFunc->GetParameter(1)-fitFunc->GetParameter(2),fitFunc->GetParameter(1)+fitFunc->GetParameter(2));
+
+ fitXMin2 = fitFunc->GetParameter(1)-1.5*fitFunc->GetParameter(2);
+ if(!SiPM_type.compare(HDR2)) fitXMin2 = fitFunc->GetParameter(1)-1.1*fitFunc->GetParameter(2);
+ fitXMax2 = fitFunc->GetParameter(1)+fabs(1.5*fitFunc->GetParameter(2));
+ fitFunc2 = new TF1(Form("fitFunc2_energyCorr_%s",labelLR_energyBin.c_str()),"gaus(0)",fitXMin2,fitXMax2);
+ fitFunc2 -> SetParameter(0,fitFunc->GetParameter(0));
+ fitFunc2 -> SetParameter(1,fitFunc->GetParameter(1));
+ fitFunc2 -> SetParameter(2,fitFunc->GetParameter(2));
+ histo -> Fit(fitFunc2,"QRSL+");
+ histo -> Fit(fitFunc2,"QRSL+");
+ z++;
+ }*/
+
+ fitFunc2 -> SetLineColor(kRed+1);
+ fitFunc2 -> SetLineWidth(3);
+ fitFunc2 -> Draw("same");
+
+
latex = new TLatex(0.20,0.85,Form("#splitline{#sigma_{raw}^{eff} = %.0f ps}{#sigma_{raw}^{gaus} = %.0f ps}",effSigma,fabs(fitFunc2->GetParameter(2))));
latex -> SetNDC();
@@ -1074,6 +1736,8 @@ int main(int argc, char** argv)
latex -> SetTextSize(0.04);
latex -> SetTextColor(kRed);
latex -> Draw("same");
+ outFile -> cd();
+ histo -> Write();
c -> Print(Form("%s/CTR_energyRatioCorr/c_deltaT_energyRatioCorr__%s.pdf",plotDir.c_str(),labelLR_energyBin.c_str()));
c -> Print(Form("%s/CTR_energyRatioCorr/c_deltaT_energyRatioCorr__%s.png",plotDir.c_str(),labelLR_energyBin.c_str()));
delete c;
@@ -1081,6 +1745,7 @@ int main(int argc, char** argv)
}
}
+
int bytes = outFile -> Write();
diff --git a/main/moduleCharacterization_step3.cpp b/main/moduleCharacterization_step3.cpp
index e4063cce..d7ad841d 100644
--- a/main/moduleCharacterization_step3.cpp
+++ b/main/moduleCharacterization_step3.cpp
@@ -27,7 +27,7 @@
#include "TLegend.h"
#include "TSpectrum.h"
-
+#include "interface/SiPM_HDR2.h"
int main(int argc, char** argv)
{
@@ -51,8 +51,8 @@ int main(int argc, char** argv)
typedef std::numeric_limits dbl;
std::cout.precision(dbl::max_digits10);
-
-//--- get parameters
+
+ //--- get parameters
std::string plotDir = opts.GetOpt("Output.plotDirStep3");
//system(Form("rm -r %s", plotDir.c_str()));
system(Form("mkdir -p %s",plotDir.c_str()));
@@ -60,22 +60,30 @@ int main(int argc, char** argv)
system(Form("mkdir -p %s/summaryPlots/tot/",plotDir.c_str()));
system(Form("mkdir -p %s/summaryPlots/energy/",plotDir.c_str()));
system(Form("mkdir -p %s/summaryPlots/timeResolution/",plotDir.c_str()));
-
- //copia il file .php che si trova in una cartella fuori plotDir in plotDir definita nel config file
- system(Form("cp %s/../index.php %s",plotDir.c_str(),plotDir.c_str()));
- system(Form("cp %s/../index.php %s/summaryPlots/",plotDir.c_str(),plotDir.c_str()));
- system(Form("cp %s/../index.php %s/summaryPlots/tot/",plotDir.c_str(),plotDir.c_str()));
- system(Form("cp %s/../index.php %s/summaryPlots/energy/",plotDir.c_str(),plotDir.c_str()));
- system(Form("cp %s/../index.php %s/summaryPlots/timeResolution/",plotDir.c_str(),plotDir.c_str()));
+ system(Form("mkdir -p %s/summaryPlots/DeltaTMean/",plotDir.c_str()));
+ system(Form("mkdir -p %s/summaryPlots/Entries_CTR/",plotDir.c_str()));
+ system(Form("mkdir -p %s/summaryPlots/Energy_linearization/",plotDir.c_str()));
//--- open files and make the tree chain
//--- Prima di eseguire cambia il numero del run nel config file
+ std::string inputDir1 = opts.GetOpt("Input.inputDir1");
std::string inputDir = opts.GetOpt("Input.inputDir");
std::string fileBaseName1 = opts.GetOpt("Input.fileBaseName1");
std::string fileBaseName2 = opts.GetOpt("Input.fileBaseName2");
std::string runs = opts.GetOpt("Input.runs");
+ std::string source = opts.GetOpt("Input.sourceName");
+ std::string Na22 = "Na22";
+ std::string SingleBarNa22 = "SingleBarNa22";
+ std::string SingleBarNa22_coinc = "SingleBarNa22_coinc";
+ std::string Co60 = "Co60";
+ std::string Co60SumPeak = "Co60SumPeak";
+ std::string Laser = "Laser";
+
+
+
+
//Define variables
std::vector listStep1;
@@ -100,7 +108,7 @@ int main(int argc, char** argv)
{
std::string fileNameStep1;
std::string fileNameStep2;
- fileNameStep1 = Form("%s%s1_%s%04d.root",inputDir.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run);
+ fileNameStep1 = Form("%s%s1_%s%04d.root",inputDir1.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run);
fileNameStep2 = Form("%s%s2_%s%04d.root",inputDir.c_str(),fileBaseName1.c_str(),fileBaseName2.c_str(),run);
std::cout << ">>> Adding file " << fileNameStep1 << std::endl;
std::cout << ">>> Adding file " << fileNameStep2 << std::endl;
@@ -113,6 +121,7 @@ int main(int argc, char** argv)
//--- define output root file
std::string outFileName = opts.GetOpt("Output.outFileNameStep3");
TFile* outFile = TFile::Open(outFileName.c_str(),"RECREATE");
+
//--- get plot settings
TCanvas* c;
@@ -123,7 +132,14 @@ int main(int argc, char** argv)
std::vector LRLabels;
LRLabels.push_back("L");
LRLabels.push_back("R");
+
+ float refTh = opts.GetOpt("Plots.refTh");
+ float refTh3 = opts.GetOpt("Plots.refTh3");
+ float refTh5 = opts.GetOpt("Plots.refTh5");
+ float refTh7 = opts.GetOpt("Plots.refTh7");
+ float refVov = opts.GetOpt("Plots.refVov");
+
std::map trees;
std::map trees2;
std::map VovLabels;
@@ -133,11 +149,14 @@ int main(int argc, char** argv)
std::map map_Vovs;
std::map map_ths;
std::map map_EnBin;
+
+ std::vector vec_th1;
+ std::vector vec_Vov;
//define TGraph
std::map g_tot_vs_th;
std::map g_tot_vs_Vov;
- std::map g_tot_vs_iBar;
+ std::map g_tot_vs_bar;
for(auto file: listStep1){
@@ -169,11 +188,19 @@ int main(int argc, char** argv)
std::string string_th = tokens[4];
string_th.erase(0,2);
map_ths[stepLabel] = atof(string_th.c_str());
+ float vth1 = map_ths[stepLabel];
+ vec_th1.push_back(vth1);
+ float Vov = map_Vovs[stepLabel];
+ vec_Vov.push_back(Vov);
}
}
std::sort(stepLabels.begin(),stepLabels.end());
stepLabels.erase(std::unique(stepLabels.begin(),stepLabels.end()),stepLabels.end());
-
+ std::sort(vec_th1.begin(),vec_th1.end());
+ vec_th1.erase(std::unique(vec_th1.begin(),vec_th1.end()),vec_th1.end());
+
+
+
for(auto stepLabel : stepLabels)
@@ -193,6 +220,7 @@ int main(int argc, char** argv)
for(auto LRLabel : LRLabels )
{
+
//label histo
std::string label(Form("bar%02d%s_%s",iBar,LRLabel.c_str(),stepLabel.c_str()));
@@ -222,11 +250,13 @@ int main(int argc, char** argv)
histo -> Fit(fitFunc1,"QNS+", " ", fitFunc1->GetParameter(1)-fitFunc1->GetParameter(2), fitFunc1->GetParameter(1)+fitFunc1->GetParameter(2));
fitFunc1 -> SetLineColor(kBlack);
fitFunc1 -> SetLineWidth(3);
- fitFunc1 -> Draw("same");
+ fitFunc1 -> Draw("same");
+ //outFile -> cd();
//histo -> Write();
- c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.png",plotDir.c_str(),label.c_str()));
+
+ /*c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.png",plotDir.c_str(),label.c_str()));
c -> Print(Form("%s/summaryPlots/tot/c_tot__%s.pdf",plotDir.c_str(),label.c_str()));
- delete c;
+ delete c;*/
if(LRLabel=="L"){
@@ -247,13 +277,13 @@ int main(int argc, char** argv)
g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())] -> SetPoint(g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())]->GetN(),Vov,fitFunc1->GetMaximumX());
g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())] -> SetPointError(g_tot_vs_Vov[Form("bar%02dL_%s",iBar,thLabel.c_str())]->GetN()-1,0,0);
- //g_tot_vs_iBar
+ //g_tot_vs_bar
- if( g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){
- g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors();
+ if( g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){
+ g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors();
}
- g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX());
- g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_iBar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0);
+ g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX());
+ g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_bar[Form("L_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0);
}
if(LRLabel=="R"){
@@ -274,19 +304,20 @@ int main(int argc, char** argv)
g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())] -> SetPoint(g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())]->GetN(),Vov,fitFunc1->GetMaximumX());
g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())] -> SetPointError(g_tot_vs_Vov[Form("bar%02dR_%s",iBar,thLabel.c_str())]->GetN()-1,0,0);
- //g_tot_vs_iBar
+ //g_tot_vs_bar
- if( g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){
- g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors();
+ if( g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] == NULL ){
+ g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] = new TGraphErrors();
}
- g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX());
- g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_iBar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0);
+ g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPoint(g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN(),iBar,fitFunc1->GetMaximumX());
+ g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())] -> SetPointError(g_tot_vs_bar[Form("R_%s_%s",VovLabel.c_str(),thLabel.c_str())]->GetN()-1,0,0);
}
}
}
}
- }
-
+ }
+
+
//------------------
//--- draw summary plots
@@ -306,13 +337,13 @@ int main(int argc, char** argv)
std::string label2(Form("bar%02dR_%s",iBar,mapIt.first.c_str()));
if(! g_tot_vs_th[Form("bar%02dL_%s",iBar,mapIt.first.c_str())]){
- std::cout<<"NON ESISTE barL "< SetMarkerColor(1+iter);
g_totL -> SetMarkerStyle(20);
g_totL -> Draw("PL,same");
+ outFile -> cd();
+ g_totL -> Write(Form("g_tot_vs_th_bar%02dL_%s",iBar,mapIt.first.c_str()));
g_totR -> SetLineColor(1+iter);
g_totR -> SetMarkerColor(1+iter);
g_totR -> SetMarkerStyle(25);
g_totR -> Draw("PL,same");
+ outFile -> cd();
+ g_totR -> Write(Form("g_tot_vs_th_bar%02dR_%s",iBar,mapIt.first.c_str()));
latex = new TLatex(0.55,0.85-0.04*iter,Form("%s",mapIt.first.c_str()));
latex -> SetNDC();
@@ -336,7 +371,6 @@ int main(int argc, char** argv)
latex -> Draw("same");
++iter;
- //std::cout<<"iter "< Print(Form("%s/summaryPlots/tot/c_tot_vs_th_bar%02d.png",plotDir.c_str(),iBar));
@@ -358,13 +392,13 @@ int main(int argc, char** argv)
std::string label2(Form("bar%02dR_%s",iBar,mapIt.first.c_str()));
if(! g_tot_vs_Vov[Form("bar%02dL_%s",iBar,mapIt.first.c_str())]){
- std::cout<<"NON ESISTE barL "< SetMarkerColor(1+iter);
g_totL -> SetMarkerStyle(20);
g_totL -> Draw("PL,same");
+ outFile -> cd();
+ g_totL -> Write(Form("g_tot_vs_Vov_bar%02dL_%s",iBar,mapIt.first.c_str()));
g_totR -> SetLineColor(1+iter);
g_totR -> SetMarkerColor(1+iter);
g_totR -> SetMarkerStyle(25);
g_totR -> Draw("PL,same");
+ outFile -> cd();
+ g_totR -> Write(Form("g_tot_vs_Vov_bar%02dR_%s",iBar,mapIt.first.c_str()));
latex = new TLatex(0.55,0.85-0.04*iter,Form("%s",mapIt.first.c_str()));
latex -> SetNDC();
@@ -388,7 +426,6 @@ int main(int argc, char** argv)
latex -> Draw("same");
++iter;
- //std::cout<<"iter "< Print(Form("%s/summaryPlots/tot/c_tot_vs_Vov_bar%02d.png",plotDir.c_str(),iBar));
@@ -398,6 +435,7 @@ int main(int argc, char** argv)
for(auto mapIt1 : thLabels){
+ if( mapIt1.first != Form("th%d",int(refTh))) continue;
for(auto mapIt2 : VovLabels){
c1 = new TCanvas(Form("c_tot_vs_bar_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()),Form("c_tot_vs_bar_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()));
@@ -410,56 +448,56 @@ int main(int argc, char** argv)
std::string label1(Form("L_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()));
std::string label2(Form("R_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()));
- if( g_tot_vs_iBar[Form("L_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){
- TGraphErrors* g_totL = g_tot_vs_iBar[label1];
- if(g_totL->GetN()>0){
- g_totL -> SetLineColor(kBlue);
- g_totL -> SetMarkerColor(kBlue);
- g_totL -> SetMarkerStyle(20);
- g_totL -> Draw("P");
-
- std::string label_prova = "Blue: L Red: R";
- latex = new TLatex(0.55,0.85-0.04,label_prova.c_str());
- latex -> SetNDC();
- latex -> SetTextFont(42);
- latex -> SetTextSize(0.04);
- latex -> SetTextColor(kBlack);
- latex -> Draw("same");
-
- c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
- c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.pdf",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
- }
+ if( g_tot_vs_bar[Form("L_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){
+ TGraphErrors* g_totL = g_tot_vs_bar[label1];
+ if(g_totL->GetN()>0){
+ g_totL -> SetLineColor(kBlue);
+ g_totL -> SetMarkerColor(kBlue);
+ g_totL -> SetMarkerStyle(20);
+ g_totL -> Draw("P");
+ outFile -> cd();
+ g_totL -> Write(Form("g_tot_vs_barL_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()));
+
+ TLegend* legenda = new TLegend(0.55,0.80,0.70,0.85);
+ legenda->AddEntry(g_totL, "Left", "P");
+ legenda -> SetBorderSize(0);
+ legenda->Draw();
+
+
+ c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
+ c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.pdf",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
+ }
}
- if(g_tot_vs_iBar[Form("R_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){
- TGraphErrors* g_totR = g_tot_vs_iBar[label2];
- if(g_totR->GetN()>0){
- g_totR -> SetLineColor(kRed);
- g_totR -> SetMarkerColor(kRed);
- g_totR -> SetMarkerStyle(20);
- g_totR -> Draw("P,same");
+ if(g_tot_vs_bar[Form("R_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str())]){
+ TGraphErrors* g_totR = g_tot_vs_bar[label2];
+ if(g_totR->GetN()>0){
+ g_totR -> SetLineColor(kRed);
+ g_totR -> SetMarkerColor(kRed);
+ g_totR -> SetMarkerStyle(20);
+ g_totR -> Draw("P,same");
+ outFile -> cd();
+ g_totR -> Write(Form("g_tot_vs_barR_%s_%s",mapIt2.first.c_str(),mapIt1.first.c_str()));
-
-
-
- c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
- c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.pdf",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
- }
- }
-
-
-
- }
- }
+ TLegend* legenda = new TLegend(0.55,0.85,0.70,0.89);
+ legenda->AddEntry(g_totR, "Right", "P");
+ legenda -> SetBorderSize(0);
+ legenda->Draw();
+
+ c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.png",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
+ c1 -> Print(Form("%s/summaryPlots/tot/c_tot_vs_bar_%s_%s.pdf",plotDir.c_str(),mapIt2.first.c_str(),mapIt1.first.c_str()));
+ }
+ }
+ }
+ }
//int d=0;
- std::vector vec_Vov;
+
std::vector vec_th;
for(auto file: listStep2){
- //int i=0;
TFile* inFile = TFile::Open(file.c_str(),"READ");
TList* list = inFile -> GetListOfKeys();
@@ -470,9 +508,9 @@ int main(int argc, char** argv)
std::string name(object->GetName());
std::vector tokens = GetTokens(name,'_');
std::size_t found;
- std::size_t found1;
+ std::size_t found1;
found = name.find("data_");
- found1 = name.find("dataRes_");
+ found1 = name.find("dataRes_");
//tree
if( found!=std::string::npos )
{
@@ -488,7 +526,6 @@ int main(int argc, char** argv)
std::string string_th = tokens[3];
string_th.erase(0,2);
map_ths[stepLabel] = atof(string_th.c_str());
- //i++;
}
std::sort(stepLabels.begin(),stepLabels.end());
@@ -508,21 +545,17 @@ int main(int argc, char** argv)
std::string string_th = tokens[3];
string_th.erase(0,2);
map_ths[stepLabel] = atof(string_th.c_str());
- float Vov = map_Vovs[stepLabel];
+ float Vov = map_Vovs[stepLabel];
float vth1 = map_ths[stepLabel];
- vec_Vov.push_back(Vov);
+ vec_Vov.push_back(Vov);
vec_th.push_back(vth1);
EnBinLabels[tokens[4]] += 1;
std::string string_EnBin = tokens[4];
string_EnBin.erase(0,9);
map_EnBin[stepLabel] = atof(string_EnBin.c_str());
- //d++;
}
}
- //std::cout<<"i="< map_en511;
std::map map_en1275;
+ std::map map_en1786;
std::map map_enRatio;
std::map map_en511R;
std::map map_en1275R;
+ std::map map_en1786R;
std::map map_enRatioR;
std::map map_en511L;
std::map map_en1275L;
+ std::map map_en1786L;
std::map map_enRatioL;
+
+
std::map g_en511_vs_th;
std::map g_en511_vs_Vov;
- std::map g_en511_vs_iBar;
+ std::map g_en511_vs_bar;
std::map g_en1275_vs_th;
std::map g_en1275_vs_Vov;
- std::map g_en1275_vs_iBar;
+ std::map g_en1275_vs_bar;
+ std::map g_en1786_vs_th;
+ std::map g_en1786_vs_Vov;
+ std::map g_en1786_vs_bar;
std::map g_enRatio_vs_th;
std::map g_enRatio_vs_Vov;
- std::map g_enRatio_vs_iBar;
+ std::map g_enRatio_vs_bar;
std::map g_en511_vs_th_L;
std::map g_en511_vs_Vov_L;
- std::map g_en511_vs_iBar_L;
+ std::map g_en511_vs_bar_L;
std::map g_en1275_vs_th_L;
std::map g_en1275_vs_Vov_L;
- std::map g_en1275_vs_iBar_L;
+ std::map g_en1275_vs_bar_L;
std::map g_enRatio_vs_th_L;
std::map g_enRatio_vs_Vov_L;
- std::map g_enRatio_vs_iBar_L;
+ std::map g_enRatio_vs_bar_L;
std::map g_en511_vs_th_R;
std::map g_en511_vs_Vov_R;
- std::map g_en511_vs_iBar_R;
+ std::map g_en511_vs_bar_R;
std::map g_en1275_vs_th_R;
std::map g_en1275_vs_Vov_R;
- std::map g_en1275_vs_iBar_R;
+ std::map g_en1275_vs_bar_R;
+ std::map g_en1786_vs_bar_R;
+ std::map g_en1786_vs_bar_L;
std::map g_enRatio_vs_th_R;
std::map g_enRatio_vs_Vov_R;
- std::map g_enRatio_vs_iBar_R;
+ std::map g_enRatio_vs_bar_R;
- //int index( (10000*int(Vov*100.)) + (100*vth1) + iBar );
-
for (auto mapIt : trees){
float energy511 = -1;
float energy1275 = -1;
+ float energy1786 = -1;
float energy511L = -1;
float energy1275L = -1;
+ float energy1786L = -1;
float energy511R = -1;
float energy1275R = -1;
+ float energy1786R = -1;
float theIndex = -1;
mapIt.second -> SetBranchStatus("*",0);
mapIt.second -> SetBranchStatus("energyPeak511LR", 1); mapIt.second -> SetBranchAddress("energyPeak511LR", &energy511);
@@ -590,6 +634,17 @@ int main(int argc, char** argv)
mapIt.second -> SetBranchStatus("energyPeak511R", 1); mapIt.second -> SetBranchAddress("energyPeak511R", &energy511R);
mapIt.second -> SetBranchStatus("energyPeak1275R", 1); mapIt.second -> SetBranchAddress("energyPeak1275R", &energy1275R);
mapIt.second -> SetBranchStatus("indexID", 1); mapIt.second -> SetBranchAddress("indexID", &theIndex);
+ if(!source.compare(SingleBarNa22)){
+ mapIt.second -> SetBranchStatus("energyPeak1786LR", 1); mapIt.second -> SetBranchAddress("energyPeak1786LR", &energy1786);
+ mapIt.second -> SetBranchStatus("energyPeak1786L", 1); mapIt.second -> SetBranchAddress("energyPeak1786L", &energy1786L);
+ mapIt.second -> SetBranchStatus("energyPeak1786R", 1); mapIt.second -> SetBranchAddress("energyPeak1786R", &energy1786R);
+ }
+ if(!source.compare(Co60)){
+ mapIt.second -> SetBranchStatus("energyPeak1786LR", 1); mapIt.second -> SetBranchAddress("energyPeak1786LR", &energy1786);
+ //mapIt.second -> SetBranchStatus("energyPeak1786L", 1); mapIt.second -> SetBranchAddress("energyPeak1786L", &energy1786L);
+ //mapIt.second -> SetBranchStatus("energyPeak1786R", 1); mapIt.second -> SetBranchAddress("energyPeak1786R", &energy1786R);
+ }
+
int nEntries = mapIt.second->GetEntries();
for(int entry = 0; entry < nEntries; ++entry)
{
@@ -601,6 +656,7 @@ int main(int argc, char** argv)
if(energy1275 > 0.1){
map_en1275[theIndex] = energy1275;
}
+
if(energy511 > 0.1 && energy1275 > 0.1){
map_enRatio[theIndex] = energy1275/energy511;
}
@@ -625,16 +681,41 @@ int main(int argc, char** argv)
map_enRatioR[theIndex] = energy1275R/energy511R;
}
+ if(!source.compare(SingleBarNa22)){
+ if(energy1786 > 0.1){
+ map_en1786[theIndex] = energy1786;
+ }
+ if(energy1786L > 0.1){
+ map_en1786L[theIndex] = energy1786L;
+ }
+ if(energy1786R > 0.1){
+ map_en1786R[theIndex] = energy1786R;
+ }
+ }
+
+ if(!source.compare(Co60)){
+ if(energy1786 > 0.1){
+ map_en1786[theIndex] = energy1786;
+ }
+ }
+
+
+
}
+
+
+
//energy Peak 511 vs th and vs Vov and vs iBar (LR, L, R)
for(std::map::iterator index = map_en511.begin(); index != map_en511.end(); index++){
//std::cout<<"index"<first<first/10000))/100);
+
+
+ Vov = float(int(index->first/10000))/100;
th = int((index->first - Vov*1000000)/100);
iBar = int(index->first - Vov*1000000 - th*100);
int Vov_iBar_ID;
@@ -644,34 +725,24 @@ int main(int argc, char** argv)
th_iBar_ID = 100*th + iBar;
Vov_th_ID =(10000*int(Vov*100.)) + (100*th);
+
+
if( g_en511_vs_th[Vov_iBar_ID] == NULL ){
g_en511_vs_th[Vov_iBar_ID] = new TGraphErrors();
}
if( g_en511_vs_Vov[th_iBar_ID] == NULL ){
g_en511_vs_Vov[th_iBar_ID] = new TGraphErrors();
}
- if( g_en511_vs_iBar[Vov_th_ID] == NULL ){
- g_en511_vs_iBar[Vov_th_ID] = new TGraphErrors();
+ if( g_en511_vs_bar[Vov_th_ID] == NULL ){
+ g_en511_vs_bar[Vov_th_ID] = new TGraphErrors();
}
- /*if( g_en511_vs_th_L[Vov_iBar_ID] == NULL ){
- g_en511_vs_th_L[Vov_iBar_ID] = new TGraphErrors();
+ if( g_en511_vs_bar_L[Vov_th_ID] == NULL ){
+ g_en511_vs_bar_L[Vov_th_ID] = new TGraphErrors();
}
- if( g_en511_vs_Vov_L[th_iBar_ID] == NULL ){
- g_en511_vs_Vov_L[th_iBar_ID] = new TGraphErrors();
- }*/
- if( g_en511_vs_iBar_L[Vov_th_ID] == NULL ){
- g_en511_vs_iBar_L[Vov_th_ID] = new TGraphErrors();
- }
-
- /*if( g_en511_vs_th_R[Vov_iBar_ID] == NULL ){
- g_en511_vs_th_R[Vov_iBar_ID] = new TGraphErrors();
- }
- if( g_en511_vs_Vov_R[th_iBar_ID] == NULL ){
- g_en511_vs_Vov_R[th_iBar_ID] = new TGraphErrors();
- }*/
- if( g_en511_vs_iBar_R[Vov_th_ID] == NULL ){
- g_en511_vs_iBar_R[Vov_th_ID] = new TGraphErrors();
+
+ if( g_en511_vs_bar_R[Vov_th_ID] == NULL ){
+ g_en511_vs_bar_R[Vov_th_ID] = new TGraphErrors();
}
if (index->second>0.1){
@@ -684,9 +755,9 @@ int main(int argc, char** argv)
//std::cout<<"N"<GetN()<<"Vov"<second<<"thIbarID"<SetPointError(g_en511_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.);
- g_en511_vs_iBar[Vov_th_ID]->SetPoint(g_en511_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second);
- //std::cout<<"N"<GetN()<<"Vov"<second<<"thVovID"<SetPointError(g_en511_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en511_vs_bar[Vov_th_ID]->SetPoint(g_en511_vs_bar[Vov_th_ID]->GetN(), iBar, index->second);
+ //std::cout<<"N"<GetN()<<"Vov"<second<<"thVovID"<SetPointError(g_en511_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.);
}
@@ -698,8 +769,8 @@ int main(int argc, char** argv)
g_en511_vs_Vov_L[th_iBar_ID]->SetPoint(g_en511_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_en511L[index->first]);
g_en511_vs_Vov_L[th_iBar_ID]->SetPointError(g_en511_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_en511_vs_iBar_L[Vov_th_ID]->SetPoint(g_en511_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_en511L[index->first]);
- g_en511_vs_iBar_L[Vov_th_ID]->SetPointError(g_en511_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en511_vs_bar_L[Vov_th_ID]->SetPoint(g_en511_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_en511L[index->first]);
+ g_en511_vs_bar_L[Vov_th_ID]->SetPointError(g_en511_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.);
}
@@ -711,21 +782,24 @@ int main(int argc, char** argv)
g_en511_vs_Vov_R[th_iBar_ID]->SetPoint(g_en511_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_en511R[index->first]);
g_en511_vs_Vov_R[th_iBar_ID]->SetPointError(g_en511_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_en511_vs_iBar_R[Vov_th_ID]->SetPoint(g_en511_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_en511R[index->first]);
- g_en511_vs_iBar_R[Vov_th_ID]->SetPointError(g_en511_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en511_vs_bar_R[Vov_th_ID]->SetPoint(g_en511_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_en511R[index->first]);
+ g_en511_vs_bar_R[Vov_th_ID]->SetPointError(g_en511_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.);
}
}
+
//energy Peak 1275 vs th and vs Vov and vs iBar (LR, L, R)
for(std::map::iterator index = map_en1275.begin(); index != map_en1275.end(); index++){
+ if(!source.compare(SingleBarNa22_coinc)) continue;
//std::cout<<"index"<first<first/10000))/100);
+
+ Vov = float(int(index->first/10000))/100;
th = int((index->first - Vov*1000000)/100);
iBar = int(index->first - Vov*1000000 - th*100);
int Vov_iBar_ID;
@@ -733,7 +807,9 @@ int main(int argc, char** argv)
int Vov_th_ID;
Vov_iBar_ID = 10000*int(Vov*100.) + iBar;
th_iBar_ID = 100*th + iBar;
- Vov_th_ID = 10000*int(Vov*100.) + 100*th;
+ Vov_th_ID =(10000*int(Vov*100.)) + (100*th);
+
+
if( g_en1275_vs_th[Vov_iBar_ID] == NULL ){
g_en1275_vs_th[Vov_iBar_ID] = new TGraphErrors();
@@ -741,26 +817,16 @@ int main(int argc, char** argv)
if( g_en1275_vs_Vov[th_iBar_ID] == NULL ){
g_en1275_vs_Vov[th_iBar_ID] = new TGraphErrors();
}
- if( g_en1275_vs_iBar[Vov_th_ID] == NULL ){
- g_en1275_vs_iBar[Vov_th_ID] = new TGraphErrors();
- }
- /*if( g_en1275_vs_th_L[Vov_iBar_ID] == NULL ){
- g_en1275_vs_th_L[Vov_iBar_ID] = new TGraphErrors();
- }
- if( g_en1275_vs_Vov_L[th_iBar_ID] == NULL ){
- g_en1275_vs_Vov_L[th_iBar_ID] = new TGraphErrors();
- }*/
- if( g_en1275_vs_iBar_L[Vov_th_ID] == NULL ){
- g_en1275_vs_iBar_L[Vov_th_ID] = new TGraphErrors();
+ if( g_en1275_vs_bar[Vov_th_ID] == NULL ){
+ g_en1275_vs_bar[Vov_th_ID] = new TGraphErrors();
}
- /*if( g_en1275_vs_th_R[Vov_iBar_ID] == NULL ){
- g_en1275_vs_th_R[Vov_iBar_ID] = new TGraphErrors();
+
+ if( g_en1275_vs_bar_L[Vov_th_ID] == NULL ){
+ g_en1275_vs_bar_L[Vov_th_ID] = new TGraphErrors();
}
- if( g_en1275_vs_Vov_R[th_iBar_ID] == NULL ){
- g_en1275_vs_Vov_R[th_iBar_ID] = new TGraphErrors();
- }*/
- if( g_en1275_vs_iBar_R[Vov_th_ID] == NULL ){
- g_en1275_vs_iBar_R[Vov_th_ID] = new TGraphErrors();
+
+ if( g_en1275_vs_bar_R[Vov_th_ID] == NULL ){
+ g_en1275_vs_bar_R[Vov_th_ID] = new TGraphErrors();
}
if (index->second>0.1){
@@ -771,8 +837,8 @@ int main(int argc, char** argv)
g_en1275_vs_Vov[th_iBar_ID]->SetPoint(g_en1275_vs_Vov[th_iBar_ID]->GetN(), Vov, index->second);
g_en1275_vs_Vov[th_iBar_ID]->SetPointError(g_en1275_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.);
- g_en1275_vs_iBar[Vov_th_ID]->SetPoint(g_en1275_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second);
- g_en1275_vs_iBar[Vov_th_ID]->SetPointError(g_en1275_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en1275_vs_bar[Vov_th_ID]->SetPoint(g_en1275_vs_bar[Vov_th_ID]->GetN(), iBar, index->second);
+ g_en1275_vs_bar[Vov_th_ID]->SetPointError(g_en1275_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.);
}
if (map_en1275R[index->first]>0.1){
@@ -783,8 +849,8 @@ int main(int argc, char** argv)
g_en1275_vs_Vov_R[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_en1275R[index->first]);
g_en1275_vs_Vov_R[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_en1275_vs_iBar_R[Vov_th_ID]->SetPoint(g_en1275_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_en1275R[index->first]);
- g_en1275_vs_iBar_R[Vov_th_ID]->SetPointError(g_en1275_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en1275_vs_bar_R[Vov_th_ID]->SetPoint(g_en1275_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_en1275R[index->first]);
+ g_en1275_vs_bar_R[Vov_th_ID]->SetPointError(g_en1275_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.);
}
if (map_en1275L[index->first]>0.1){
@@ -795,19 +861,22 @@ int main(int argc, char** argv)
g_en1275_vs_Vov_L[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_en1275L[index->first]);
g_en1275_vs_Vov_L[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_en1275_vs_iBar_L[Vov_th_ID]->SetPoint(g_en1275_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_en1275L[index->first]);
- g_en1275_vs_iBar_L[Vov_th_ID]->SetPointError(g_en1275_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_en1275_vs_bar_L[Vov_th_ID]->SetPoint(g_en1275_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_en1275L[index->first]);
+ g_en1275_vs_bar_L[Vov_th_ID]->SetPointError(g_en1275_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.);
}
}
//energy Peak Ratio vs th and vs Vov vs iBar (LR, L, R)
for(std::map::iterator index = map_enRatio.begin(); index != map_enRatio.end(); index++){
+ if(!source.compare(SingleBarNa22_coinc)) continue;
//std::cout<<"index"<first<first/10000))/100);
+
+ Vov = float(int(index->first/10000))/100;
+
th = int((index->first - Vov*1000000)/100);
iBar = int(index->first - Vov*1000000 - th*100);
int Vov_iBar_ID;
@@ -817,6 +886,8 @@ int main(int argc, char** argv)
th_iBar_ID = 100*th + iBar;
Vov_th_ID = 10000*int(Vov*100.) + 100*th;
+
+
if( g_enRatio_vs_th[Vov_iBar_ID] == NULL ){
g_enRatio_vs_th[Vov_iBar_ID] = new TGraphErrors();
}
@@ -824,30 +895,16 @@ int main(int argc, char** argv)
g_enRatio_vs_Vov[th_iBar_ID] = new TGraphErrors();
}
- if( g_enRatio_vs_iBar[Vov_th_ID] == NULL ){
- g_enRatio_vs_iBar[Vov_th_ID] = new TGraphErrors();
+ if( g_enRatio_vs_bar[Vov_th_ID] == NULL ){
+ g_enRatio_vs_bar[Vov_th_ID] = new TGraphErrors();
}
- /*if( g_enRatio_vs_th_L[Vov_iBar_ID] == NULL ){
- g_enRatio_vs_th_L[Vov_iBar_ID] = new TGraphErrors();
- }
- if( g_enRatio_vs_Vov_L[th_iBar_ID] == NULL ){
- g_enRatio_vs_Vov_L[th_iBar_ID] = new TGraphErrors();
- }*/
-
- if( g_enRatio_vs_iBar_L[Vov_th_ID] == NULL ){
- g_enRatio_vs_iBar_L[Vov_th_ID] = new TGraphErrors();
+ if( g_enRatio_vs_bar_L[Vov_th_ID] == NULL ){
+ g_enRatio_vs_bar_L[Vov_th_ID] = new TGraphErrors();
}
-
- /* if( g_enRatio_vs_th_R[Vov_iBar_ID] == NULL ){
- g_enRatio_vs_th_R[Vov_iBar_ID] = new TGraphErrors();
- }
- if( g_enRatio_vs_Vov_R[th_iBar_ID] == NULL ){
- g_enRatio_vs_Vov_R[th_iBar_ID] = new TGraphErrors();
- }*/
-
- if( g_enRatio_vs_iBar_R[Vov_th_ID] == NULL ){
- g_enRatio_vs_iBar_R[Vov_th_ID] = new TGraphErrors();
+
+ if( g_enRatio_vs_bar_R[Vov_th_ID] == NULL ){
+ g_enRatio_vs_bar_R[Vov_th_ID] = new TGraphErrors();
}
if (index->second>0.1){
@@ -858,8 +915,8 @@ int main(int argc, char** argv)
g_enRatio_vs_Vov[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov[th_iBar_ID]->GetN(), Vov, index->second);
g_enRatio_vs_Vov[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.);
- g_enRatio_vs_iBar[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar[Vov_th_ID]->GetN(), iBar, index->second);
- g_enRatio_vs_iBar[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_enRatio_vs_bar[Vov_th_ID]->SetPoint(g_enRatio_vs_bar[Vov_th_ID]->GetN(), iBar, index->second);
+ g_enRatio_vs_bar[Vov_th_ID]->SetPointError(g_enRatio_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.);
}
if (map_enRatioR[index->first]>0.1){
@@ -870,8 +927,8 @@ int main(int argc, char** argv)
g_enRatio_vs_Vov_R[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_enRatioR[index->first]);
g_enRatio_vs_Vov_R[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_enRatio_vs_iBar_R[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar_R[Vov_th_ID]->GetN(), iBar, map_enRatioR[index->first]);
- g_enRatio_vs_iBar_R[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar_R[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_enRatio_vs_bar_R[Vov_th_ID]->SetPoint(g_enRatio_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_enRatioR[index->first]);
+ g_enRatio_vs_bar_R[Vov_th_ID]->SetPointError(g_enRatio_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.);
}
if (map_enRatioL[index->first]>0.1){
@@ -882,46 +939,137 @@ int main(int argc, char** argv)
g_enRatio_vs_Vov_L[th_iBar_ID]->SetPoint(g_enRatio_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_enRatioL[index->first]);
g_enRatio_vs_Vov_L[th_iBar_ID]->SetPointError(g_enRatio_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/
- g_enRatio_vs_iBar_L[Vov_th_ID]->SetPoint(g_enRatio_vs_iBar_L[Vov_th_ID]->GetN(), iBar, map_enRatioL[index->first]);
- g_enRatio_vs_iBar_L[Vov_th_ID]->SetPointError(g_enRatio_vs_iBar_L[Vov_th_ID]->GetN()-1, 0., 0.);
+ g_enRatio_vs_bar_L[Vov_th_ID]->SetPoint(g_enRatio_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_enRatioL[index->first]);
+ g_enRatio_vs_bar_L[Vov_th_ID]->SetPointError(g_enRatio_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.);
}
}
+
+
+//energy Peak 1786 vs th and vs Vov and vs iBar (LR, L, R)
+ if(!source.compare(SingleBarNa22) || !source.compare(Co60)){
+ for(std::map::iterator index = map_en1786.begin(); index != map_en1786.end(); index++){
+ if(!source.compare(SingleBarNa22_coinc)) continue;
+ //std::cout<<"index"<first<first/10000))/100;
+
+ th = int((index->first - Vov*1000000)/100);
+ iBar = int(index->first - Vov*1000000 - th*100);
+ int Vov_iBar_ID;
+ int th_iBar_ID;
+ int Vov_th_ID;
+ Vov_iBar_ID = 10000*int(Vov*100.) + iBar;
+ th_iBar_ID = 100*th + iBar;
+ Vov_th_ID = 10000*int(Vov*100.) + 100*th;
+
+ if( g_en1786_vs_th[Vov_iBar_ID] == NULL ){
+ g_en1786_vs_th[Vov_iBar_ID] = new TGraphErrors();
+ }
+ if( g_en1786_vs_Vov[th_iBar_ID] == NULL ){
+ g_en1786_vs_Vov[th_iBar_ID] = new TGraphErrors();
+ }
+ if( g_en1786_vs_bar[Vov_th_ID] == NULL ){
+ g_en1786_vs_bar[Vov_th_ID] = new TGraphErrors();
+ }
+
+ if( g_en1786_vs_bar_L[Vov_th_ID] == NULL ){
+ g_en1786_vs_bar_L[Vov_th_ID] = new TGraphErrors();
+ }
+
+ if( g_en1786_vs_bar_R[Vov_th_ID] == NULL ){
+ g_en1786_vs_bar_R[Vov_th_ID] = new TGraphErrors();
+ }
+
+ if (index->second>0.1){
+ //1275 Peak LR
+ g_en1786_vs_th[Vov_iBar_ID]->SetPoint(g_en1786_vs_th[Vov_iBar_ID]->GetN(), th, index->second);
+ g_en1786_vs_th[Vov_iBar_ID]->SetPointError(g_en1786_vs_th[Vov_iBar_ID]->GetN()-1, 0., 0.);
+
+ g_en1786_vs_Vov[th_iBar_ID]->SetPoint(g_en1786_vs_Vov[th_iBar_ID]->GetN(), Vov, index->second);
+ g_en1786_vs_Vov[th_iBar_ID]->SetPointError(g_en1786_vs_Vov[th_iBar_ID]->GetN()-1, 0., 0.);
+
+ g_en1786_vs_bar[Vov_th_ID]->SetPoint(g_en1786_vs_bar[Vov_th_ID]->GetN(), iBar, index->second);
+ g_en1786_vs_bar[Vov_th_ID]->SetPointError(g_en1786_vs_bar[Vov_th_ID]->GetN()-1, 0., 0.);
+ }
+
+ if(!source.compare(SingleBarNa22)){
+
+ if (map_en1786R[index->first]>0.1){
+ //1275 Peak Right
+ /* g_en1275_vs_th_R[Vov_iBar_ID]->SetPoint(g_en1275_vs_th_R[Vov_iBar_ID]->GetN(), th, map_en1275R[index->first]);
+ g_en1275_vs_th_R[Vov_iBar_ID]->SetPointError(g_en1275_vs_th_R[Vov_iBar_ID]->GetN()-1, 0., 0.);
+
+ g_en1275_vs_Vov_R[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_R[th_iBar_ID]->GetN(), Vov, map_en1275R[index->first]);
+ g_en1275_vs_Vov_R[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_R[th_iBar_ID]->GetN()-1, 0., 0.);*/
+ g_en1786_vs_bar_R[Vov_th_ID]->SetPoint(g_en1786_vs_bar_R[Vov_th_ID]->GetN(), iBar, map_en1786R[index->first]);
+ g_en1786_vs_bar_R[Vov_th_ID]->SetPointError(g_en1786_vs_bar_R[Vov_th_ID]->GetN()-1, 0., 0.);
+ }
+
+ if (map_en1786L[index->first]>0.1){
+ //1275 Peak Left
+ /* g_en1275_vs_th_L[Vov_iBar_ID]->SetPoint(g_en1275_vs_th_L[Vov_iBar_ID]->GetN(), th, map_en1275L[index->first]);
+ g_en1275_vs_th_L[Vov_iBar_ID]->SetPointError(g_en1275_vs_th_L[Vov_iBar_ID]->GetN()-1, 0., 0.);
+
+ g_en1275_vs_Vov_L[th_iBar_ID]->SetPoint(g_en1275_vs_Vov_L[th_iBar_ID]->GetN(), Vov, map_en1275L[index->first]);
+ g_en1275_vs_Vov_L[th_iBar_ID]->SetPointError(g_en1275_vs_Vov_L[th_iBar_ID]->GetN()-1, 0., 0.);*/
+
+ g_en1786_vs_bar_L[Vov_th_ID]->SetPoint(g_en1786_vs_bar_L[Vov_th_ID]->GetN(), iBar, map_en1786L[index->first]);
+ g_en1786_vs_bar_L[Vov_th_ID]->SetPointError(g_en1786_vs_bar_L[Vov_th_ID]->GetN()-1, 0., 0.);
+ }
+ }
+ }
+ }
+
+
std::map c_en511_vs_th;
std::map c_en511_vs_Vov;
- std::map c_en511_vs_iBar;
+ std::map c_en511_vs_bar;
std::map c_en1275_vs_th;
std::map c_en1275_vs_Vov;
- std::map c_en1275_vs_iBar;
+ std::map c_en1275_vs_bar;
std::map c_enRatio_vs_th;
std::map c_enRatio_vs_Vov;
- std::map c_enRatio_vs_iBar;
+ std::map c_enRatio_vs_bar;
+ std::map c_en1786_vs_th;
+ std::map c_en1786_vs_Vov;
+ std::map c_en1786_vs_bar;
+
std::map c_en511_vs_th_L;
std::map c_en511_vs_Vov_L;
- std::map c_en511_vs_iBar_L;
+ std::map c_en511_vs_bar_L;
std::map c_en1275_vs_th_L;
std::map c_en1275_vs_Vov_L;
- std::map c_en1275_vs_iBar_L;
+ std::map c_en1275_vs_bar_L;
+ std::map c_en1786_vs_th_L;
+ std::map c_en1786_vs_Vov_L;
+ std::map c_en1786_vs_bar_L;
std::map c_enRatio_vs_th_L;
std::map c_enRatio_vs_Vov_L;
- std::map c_enRatio_vs_iBar_L;
+ std::map c_enRatio_vs_bar_L;
std::map c_en511_vs_th_R;
std::map c_en511_vs_Vov_R;
- std::map c_en511_vs_iBar_R;
+ std::map c_en511_vs_bar_R;
std::map c_en1275_vs_th_R;
std::map c_en1275_vs_Vov_R;
- std::map c_en1275_vs_iBar_R;
+ std::map c_en1275_vs_bar_R;
+ std::map c_en1786_vs_th_R;
+ std::map c_en1786_vs_Vov_R;
+ std::map c_en1786_vs_bar_R;
std::map c_enRatio_vs_th_R;
std::map c_enRatio_vs_Vov_R;
- std::map