#include <iostream>
#include <Riostream.h>
#include <TROOT.h>
#include <TObject.h>
#include <TApplication.h>
#include <TBranch.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TClass.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TTree.h>
#include <TVirtualFFT.h>
#include "TEpicData.h"
#define DEBUG 0
#define z_bin_length 22.0 // [um]
#define Electron_Vdrift 110.0 // [um/ns]
#define AlphaDecay_TimeWindow 100 // us
#define Convolution_Gaussian_Sigma_ns 10
#define Convolution_Gaussian_Range_nSigma 3
void run(int fileno=1)
{
// === ==========
// === input data
// === ==========
TChain * t = new TChain("SimulatedTree");
t->Add(Form("/media/audrey/DATA1/EPIC/simulation/K11_alpha240Pu_window%ius_%i.root",AlphaDecay_TimeWindow,fileno));
t->ls();
// === ===================
// === Initialize the data
// === ===================
TClass *epicClass = TClass::GetClass("TEpicData");
if (!epicClass) {
std::cerr << "ERROR: TEpicData Class is not loaded" << std::endl;
return;
}
TEpicData *m_Data = nullptr;
t->SetBranchAddress("Epic", &m_Data);
// === =========
// === Variables
// === =========
double time_bin_width_ns = z_bin_length / Electron_Vdrift ;
double time_KtoA_ns = 2500 / Electron_Vdrift ;
double time_window_ns = AlphaDecay_TimeWindow * 1000 ;
double range_gauss = Convolution_Gaussian_Sigma_ns * Convolution_Gaussian_Range_nSigma * 1000; // n sigma in ps
int nbins_QvsT = ceil( (time_window_ns+time_KtoA_ns) / time_bin_width_ns ) ;
int nbins_gauss = (2*range_gauss) / (1000*time_bin_width_ns);
int nbins_convol = nbins_QvsT + 2*nbins_gauss;
cout << "nbins_gauss = " << nbins_gauss << " to cover +/- " << Convolution_Gaussian_Range_nSigma << " sigma" << endl;
// === ==========
// === Histograms
// === ==========
TH1D * h1_Qmax = new TH1D(Form("Qmax_per_AlphaDecayWindowOf%ius",AlphaDecay_TimeWindow),
Form("Qmax_per_AlphaDecayWindowOf%ius",AlphaDecay_TimeWindow),100000,0,100000);
#if DEBUG
TH1D * h1_QvsT = new TH1D("influence_vs_time_1D","influence_vs_time_1D",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
h1_QvsT->GetXaxis()->SetTitle("Absolute time [ps] 200ps/bin");
h1_QvsT->SetLineColor(kBlue);
h1_QvsT->SetDirectory(0);
#endif
// --- Gaussian Distribution for convolution
vector<double> gaussian_distribution(nbins_gauss,0.);
double normalization_factor = 0;
for (int i = 0; i < nbins_gauss; i++) {
gaussian_distribution[i] = TMath::Gaus(-range_gauss+i*time_bin_width_ns*1000, 0, Convolution_Gaussian_Sigma_ns*1000, kTRUE);
normalization_factor += gaussian_distribution[i];
}
for(int i=0; i<nbins_gauss; i++) gaussian_distribution[i] /= normalization_factor;
#if DEBUG
TH1D * h1_convolution_QvsT = new TH1D("QvsT_convol","QvsT_convol",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
h1_convolution_QvsT->SetLineColor(kRed);
#endif
// === ============================================================================
// === Loop over the entries
unsigned int nentries = t->GetEntries();
std::cout << "nentries = " << nentries << std::endl;
#if DEBUG
for (unsigned int Entry = 0; Entry < 1; Entry++) {
#else
for (unsigned int Entry = 0; Entry < nentries; Entry++) {
if((Entry%250)==0) cout << "\r === Entry = " << Entry << " ===" << flush;
#endif
int bytesRead = t->GetEntry(Entry);
if (bytesRead <= 0) {
std::cerr << "Error : GetEntry failed at Entry #" << Entry << std::endl;
continue;
}
if (!m_Data) {
std::cerr << "Error: m_Data is NULL at Entry #" << Entry << std::endl;
continue;
}
#if DEBUG
if(m_Data->GetMultiplicity()==0) continue;
cout << endl ;
cout << "***********************************************" << endl;
cout << "@ Entry = " << Entry << ", mult = " << m_Data->GetMultiplicity() << endl;
#endif
vector<double> influence(nbins_convol,0);
// loop over the particle multiplicity
for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
for(int smult=0; smult<m_Data->GetNumberOfSteps(pmult); smult++){
double t_Step = m_Data->GetTimeCreationElectronsPerStep(pmult).at(smult) ; // [ns]
int n_Step = m_Data->GetNumElectronsPerStep(pmult).at(smult);
int b_Step = t_Step/time_bin_width_ns;
for(int step=smult; step<114; step++){
if(nbins_QvsT < b_Step+step-smult){
cout << "size of the vector: " << nbins_QvsT
<< " but index to fill vector is " << b_Step+step-smult << endl;
}
else {
influence[b_Step+step-smult+nbins_gauss] += n_Step * 22. / 2500. ;
}
}
}
}// end of loop over the EventMult one mult per track
double Qmax = 0;
for(int i=0; i<nbins_QvsT+nbins_gauss; i++){
double convolution = 0;
for (int j = 0; j < nbins_gauss; j++) convolution += influence[i+j] * gaussian_distribution[j];
if(convolution>Qmax) Qmax=convolution;
#if DEBUG
h1_QvsT->SetBinContent(i+1,influence[i]);
h1_QvsT->SetBinError(i+1,0);
h1_convolution_QvsT->SetBinContent(i+1+0.5*nbins_gauss,convolution);
#endif
}
if(Qmax>h1_Qmax->GetBinCenter(h1_Qmax->GetNbinsX())) cout << "Qmax out of range" << endl;
h1_Qmax->Fill(Qmax);
//bool kProcess = kFALSE;
//for(int bin=1; bin<=h1_QvsT->GetNbinsX(); bin++){
// if( kProcess==kFALSE && h1_QvsT->GetBinContent(bin)>0 ){
// kProcess = kTRUE ;
// }
// if(kProcess==kTRUE){
// // loop over the 50 ns
// // one bin is 200 ps
// // 50 ns is 250 bins
// double Q1_50ns = 0;
// for(int i=bin; i<=bin+250;i++) Q1_50ns += h1_QvsT->GetBinContent(i);
// h1_Q1_50ns->Fill(Q1_50ns);
// // reinit
// bin += 249;
// Q1_50ns = 0 ;
// kProcess = kFALSE;
// }
//}
influence.clear();
}
// === End of loop over the entries
// === ============================================================================
//// === DRAW
#if DEBUG
TCanvas * can1D_QvsT = new TCanvas("1D_QvsT_100us_TimeWindow","One Event QvsT_100us_TimeWindow",0,0,2500,2500);
h1_QvsT->Draw();
h1_convolution_QvsT->Draw("same");
#endif
TCanvas * can_Qmax = new TCanvas("Qmax","Qmax",0,0,2000,1500);
can_Qmax->cd();
h1_Qmax->Draw();
TFile * fsave = new TFile(Form("out_K11_alpha240Pu_window%ius_%i.root",AlphaDecay_TimeWindow,fileno),"recreate");
h1_Qmax->Write();
fsave->Close();
}