// Vijay Kumar, Univerity of Regina - 24/07/20
// vijay36361@gmail.com
#ifndef calibration_h
#define calibration_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>
#include <TH1.h>
#include <TH2.h>
#include <TNtuple.h>
const Int_t fhgc_pmts = 4;
const Double_t fhgc_zpos = 156.27;
class calibration : public TSelector {
public :
TTreeReader fReader;
TTree *fChain = 0; //!pointer to the analyzed TTree or TChain
Bool_t fFullRead;
Bool_t fFullShow;
Bool_t fTrack;
Bool_t fCut;
Bool_t fPions;
Bool_t frun;
// Declaration of histograms
TH1F **fPulseInt;
TH1F **fPulseInt_poiss;
TH1F ***fPulseInt_quad;
TH1F *fBeta_Cut;
TH1F *fBeta_Full;
TH1F *fTiming_Full;
TH1F *fTim1;
TH1F *fTim2;
TH1F *fTim3;
TH1F *fTim4;
TH1F *fTim1_full;
TH1F *fTim2_full;
TH1F *fTim3_full;
TH1F *fTim4_full;
// Declaration of histograms used in fitting/analysis
TH1F *scaled_clone;
TH1F *fscaled[4];
TH1F *fscaled_nobackground[4];
TH1F *fscaled_mk2[4];
TH1F *fscaled_mk2_nobackground[4];
TH1F *fscaled_temp[4];
TH1F *fscaled_combined[4];
TH1F *fscaled_total;
TH1F *fscaled_temp_mk2[4];
TH1F *fscaled_combined_mk2[4];
TH1F *fscaled_total_mk2;
// Declaration of canvases used in fitting/analysis
TCanvas *t;
TCanvas *quad_cuts_ipmt;
TCanvas *quad_cuts[4];
TCanvas *low_stats_ipmt;
TCanvas *background_ipmt;
TCanvas *Full_zoom_fit_ipmt;
TCanvas *final_spectra_ipmt;
TCanvas *background_mk2_ipmt;
TCanvas *final_spectra_mk2_ipmt;
TCanvas *final_spectra_combined;
TCanvas *final_spectra_combined_mk2;
TCanvas *scaled_poisson;
TCanvas *scaled_total;
// Declaration of preprocessing quantities
Double_t timing_mean[4];
Double_t timing_std[4];
Double_t x;
// Readers to access the data
// These leaves MUST all be present in your replay file for this scrip to run!
TTreeReaderValue<Int_t> Ndata_P_tr_beta = {fReader, "Ndata.P.tr.beta"};
TTreeReaderArray<Double_t> P_tr_beta = {fReader, "P.tr.beta"};
TTreeReaderArray<Double_t> P_hgcer_goodAdcTdcDiffTime = {fReader, "P.hgcer.goodAdcTdcDiffTime"};
TTreeReaderArray<Double_t> P_hgcer_goodAdcPulseInt = {fReader, "P.hgcer.goodAdcPulseInt"};
TTreeReaderArray<Double_t> P_hgcer_goodAdcPulseAmp = {fReader, "P.hgcer.goodAdcPulseAmp"};
TTreeReaderArray<Double_t> P_hgcer_numTracksFired = {fReader, "P.hgcer.numTracksFired"};
TTreeReaderValue<Double_t> P_cal_fly_earray = {fReader, "P.cal.fly.earray"};
TTreeReaderValue<Double_t> P_cal_pr_eplane = {fReader, "P.cal.pr.eplane"};
TTreeReaderValue<Double_t> P_cal_etotnorm = {fReader, "P.cal.etotnorm"};
TTreeReaderValue<Double_t> P_gtr_p = {fReader, "P.gtr.p"};
TTreeReaderArray<Double_t> P_gtr_dp = {fReader, "P.gtr.dp"};
TTreeReaderArray<Double_t> P_tr_x = {fReader, "P.tr.x"};
TTreeReaderArray<Double_t> P_tr_ph = {fReader, "P.tr.ph"};
TTreeReaderArray<Double_t> P_tr_y = {fReader, "P.tr.y"};
TTreeReaderArray<Double_t> P_tr_th = {fReader, "P.tr.th"};
TTreeReaderArray<Double_t> P_hgcer_xAtCer = {fReader, "P.hgcer.xAtCer"};
TTreeReaderArray<Double_t> P_hgcer_yAtCer = {fReader, "P.hgcer.yAtCer"};
calibration(TTree * /*tree*/ =0) : fChain(0) {fPulseInt = 0, fPulseInt_poiss = 0, fPulseInt_quad = 0, fBeta_Cut = 0, fBeta_Full = 0, fTiming_Full = 0,fTim1 =0, fTim1_full = 0,fTim2 =0, fTim2_full = 0, fTim3 = 0, fTim3_full = 0, fTim4 = 0, fTim4_full = 0, fFullRead = kFALSE, fFullShow = kFALSE, fTrack = kFALSE, fCut = kFALSE, fPions = kFALSE;}
virtual ~calibration() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
virtual void Init(TTree *tree);
virtual Bool_t Notify();
virtual Bool_t Process(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
virtual void SetOption(const char *option) { fOption = option; }
virtual void SetObject(TObject *obj) { fObject = obj; }
virtual void SetInputList(TList *input) { fInput = input; }
virtual TList *GetOutputList() const { return fOutput; }
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(calibration,0);
};
#endif
#ifdef calibration_cxx
void calibration::Init(TTree *tree)
{
fReader.SetTree(tree);
}
Bool_t calibration::Notify()
{
return kTRUE;
}
//Poisson distribution is used to remove background from larger values of NPE
Double_t poisson(Double_t *x, Double_t *par)
{
Double_t PoissFit1 = (par[1]*pow(par[0],x[0])*exp(-par[0]))/(tgamma(x[0]+1));
return PoissFit1;
}
//Gaussian distribution is used to find the mean of the SPE and determine spacing between subsequent NPE
Double_t gauss(Double_t *x, Double_t *par)
{
Double_t GaussFit1 = par[0]*exp((-0.5)*(pow((x[0] - par[1]),2)/pow((par[2]),2)));
Double_t GaussFit2 = par[3]*exp((-0.5)*(pow((x[0] - par[4]),2)/pow((par[5]),2)));
return GaussFit1 + GaussFit2;
}
// Function used for quality contron of the calibration
Double_t fun_4gauss_2poisson(Double_t *x, Double_t *par)
{
Double_t GaussFit1 = par[0]*exp((-0.5)*(pow((x[0] - par[1]),2)/pow((par[2]),2)));
Double_t GaussFit2 = par[3]*exp((-0.5)*(pow((x[0] - par[4]),2)/pow((par[5]),2)));
Double_t GaussFit3 = par[6]*exp((-0.5)*(pow((x[0] - par[7]),2)/pow((par[8]),2)));
Double_t GaussFit4 = par[9]*exp((-0.5)*(pow((x[0] - par[10]),2)/pow((par[11]),2)));
Double_t PoissFit1 = (par[13]*pow(par[12],x[0])*exp(-par[12]))/(tgamma(x[0]+1));
Double_t PoissFit2 = (par[15]*pow(par[14],x[0])*exp(-par[14]))/(tgamma(x[0]+1));
return GaussFit1 + GaussFit2 + GaussFit3 + GaussFit4 + PoissFit1 + PoissFit2;
}
//A simple linear equation is used to determine how linear the means of the NPE are
Double_t linear(Double_t *x, Double_t *par)
{
Double_t LinFit1 = par[0]*x[0] + par[1];
return LinFit1;
}
#endif // #ifdef calibration_cxx