#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 Convolution_Gaussian_Sigma_ns 10
#define Convolution_Gaussian_Range_nSigma 3
void run()
{
// === ==========
// === input data
// === ==========
TChain * t = new TChain("SimulatedTree");
for (int file=1; file<=10; file++){
if (file != 3)
t->Add(Form("/media/audrey/DATA1/EPIC/simulation/K11_240Pu_nf_GEF_part%i.root",file));
}
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;
int nbins_QvsT;
int nbins_convol;
// --- Gaussian Distribution for convolution
double range_gauss = Convolution_Gaussian_Sigma_ns * Convolution_Gaussian_Range_nSigma * 1000; // n sigma in ps
int nbins_gauss = (2*range_gauss) / (1000*time_bin_width_ns);
cout << "nbins_gauss = " << nbins_gauss << " to cover +/- " << Convolution_Gaussian_Range_nSigma << " sigma" << endl;
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;
// === ==========
// === Histograms
// === ==========
TH1I * h1_AnodeNumber = new TH1I("AnodeNumber","AnodeNumber",120,-1.5, 119.5);
TH1I * h1_MultInternal = new TH1I("MultInternal","MultInternal",6,-1.5,4.5);
TH1I * h1_MultExternal = new TH1I("MultExternal","MultExternal",6,-1.5,4.5);
TH1I * h1_MultTot = new TH1I("MultTot","MultTot",7,-1.5,5.5);
TH2I * h2_MultTot_vs_MultInternal = new TH2I("MultTotal_vs_MultInternal","MultTotal_vs_MultInternal",6,-1.5,4.5,7,-1.5,5.5);
TH2I * h2_MultTot_vs_MultExternal = new TH2I("MultTotal_vs_MultExternal","MultTotal_vs_MultExternal",6,-1.5,4.5,7,-1.5,5.5);
TH1D * h1_QmaxInternal = new TH1D("QmaxInternalAnode_per_GEFevent","QmaxInternalAnode_per_GEFevent",100000,0,100000);
//// === DRAW
#if DEBUG
TCanvas * can1D_QvsT = new TCanvas("1D_QvsT","One Event",0,0,2500,2500);
#endif
// === ============================================================================
// === Loop over the entries
unsigned int nentries = t->GetEntries();
std::cout << "nentries = " << nentries << std::endl;
#if DEBUG
for (unsigned int Entry = 0; Entry < 100; Entry++) {
#else
//for (unsigned int Entry = 0; Entry < nentries; Entry++) {
for (unsigned int Entry = 0; Entry < 1000000; Entry++) {
if((Entry%100000)==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
cout << endl ;
cout << "***********************************************" << endl;
cout << "@ Entry = " << Entry << ", mult = " << m_Data->GetMultiplicity() << endl;
#endif
int mult_internal = 0;
int mult_external = 0;
h1_MultTot->Fill(m_Data->GetMultiplicity());
if(m_Data->GetMultiplicity()==0) {
h1_MultInternal->Fill(0);
h1_MultExternal->Fill(0);
}
else if(m_Data->GetMultiplicity()!=1) {
#if DEBUG
cout << "CASE OF 2FF : m_Data->GetMultiplicity() = " << m_Data->GetMultiplicity() << endl;
#endif
double time_window_ns_max = 0;
for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
h1_AnodeNumber->Fill(m_Data->GetAnodeNbr(pmult));
if(m_Data->GetAnodeNbr(pmult)<100) mult_internal++;
else mult_external++;
// FIXME: consider the external anode in an independant variable
if(m_Data->GetAnodeNbr(pmult)>=100) continue;
int nsteps = m_Data->GetNumberOfSteps(pmult);
double tfirst = m_Data->GetTimeCreationElectronsPerStep(pmult).at(0);
double tlast = m_Data->GetTimeCreationElectronsPerStep(pmult).at(nsteps-1);
time_window_ns = max(tfirst+time_KtoA_ns,tlast+(114-nsteps)*time_bin_width_ns);
if(time_window_ns>time_window_ns_max) time_window_ns_max = time_window_ns;
}
nbins_QvsT = ceil( time_window_ns_max / time_bin_width_ns );
nbins_convol = nbins_QvsT + 2*nbins_gauss;
vector<double> influence(nbins_convol,0);
#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);
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
for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
// FIXME: consider the external anode in an independant variable
if(m_Data->GetAnodeNbr(pmult)>=100) continue;
int nsteps = m_Data->GetNumberOfSteps(pmult);
for(int smult=0; smult<nsteps; 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;
// FIXME
//for(int step=smult; step<114; step++){
for(int step=b_Step; step<114; step++){
if(nbins_QvsT < b_Step+step-smult){
cout << endl;
cout << "=== size of the influence vector: " << nbins_QvsT
<< " but index to fill vector is " << b_Step+step-smult << endl;
cout << " ANODE " << m_Data->GetAnodeNbr(pmult) << ", mult_tot = " << m_Data->GetMultiplicity()
<< ", smult / nsteps : " << smult << " / " << nsteps << " : b_Step = " << b_Step << endl;
}
else {
influence[b_Step+step-smult+nbins_gauss] += n_Step * 22. / 2500. ;
}
}// end of loop over the steps up to anodes
}// end of loop over the steps
}// end of loop over the particle multiplicity (if 2 FF)
double QmaxInternal = 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 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 (convolution>QmaxInternal) QmaxInternal=convolution;
}
h1_QmaxInternal->Fill(QmaxInternal);
influence.clear();
#if DEBUG
can1D_QvsT->cd();
h1_QvsT->Draw();
h1_convolution_QvsT->Draw("same");
can1D_QvsT->Update();
can1D_QvsT->WaitPrimitive();
h1_QvsT->Reset();
h1_convolution_QvsT->Reset();
#endif
}// end of if else mult>1
else{
h1_AnodeNumber->Fill(m_Data->GetAnodeNbr(0));
if(m_Data->GetAnodeNbr(0)<100) mult_internal++;
else mult_external++;
// FIXME: consider the external anode in an independant variable
if(m_Data->GetAnodeNbr(0)<100) {
int nsteps = m_Data->GetNumberOfSteps(0);
double tfirst = m_Data->GetTimeCreationElectronsPerStep(0).at(0);
double tlast = m_Data->GetTimeCreationElectronsPerStep(0).at(nsteps-1);
time_window_ns = max(tfirst+time_KtoA_ns,tlast+(114-nsteps)*time_bin_width_ns);
nbins_QvsT = ceil( time_window_ns / time_bin_width_ns );
nbins_convol = nbins_QvsT + 2*nbins_gauss;
vector<double> influence(nbins_convol,0);
#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);
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
for(int smult=0; smult<nsteps; smult++){
double t_Step = m_Data->GetTimeCreationElectronsPerStep(0).at(smult) ; // [ns]
int n_Step = m_Data->GetNumElectronsPerStep(0).at(smult);
int b_Step = t_Step/time_bin_width_ns;
// FIXME
//for(int step=smult; step<114; step++){
for(int step=b_Step; step<114; step++){
if(nbins_QvsT < b_Step+step-smult){
cout << endl;
cout << "--- size of the vector: " << nbins_QvsT
<< " but index to fill vector is " << b_Step+step-smult << endl;
cout << " ANODE " << m_Data->GetAnodeNbr(0) << ", mult_tot = " << m_Data->GetMultiplicity()
<< ", smult / nsteps : " << smult << " / " << nsteps << " : b_Step = " << b_Step << endl;
}
else {
influence[b_Step+step-smult+nbins_gauss] += n_Step * 22. / 2500. ;
}
}// end of loop over the steps up to anodes
}// end of loop over the steps
double QmaxInternal = 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 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 (convolution>QmaxInternal) QmaxInternal=convolution;
}
h1_QmaxInternal->Fill(QmaxInternal);
influence.clear();
#if DEBUG
can1D_QvsT->cd();
h1_QvsT->Draw();
h1_convolution_QvsT->Draw("same");
can1D_QvsT->Update();
can1D_QvsT->WaitPrimitive();
h1_QvsT->Reset();
h1_convolution_QvsT->Reset();
#endif
}
}// end of else mult=1
h1_MultInternal->Fill(mult_internal);
h1_MultExternal->Fill(mult_external);
h2_MultTot_vs_MultInternal->Fill(mult_internal,m_Data->GetMultiplicity());
h2_MultTot_vs_MultExternal->Fill(mult_external,m_Data->GetMultiplicity());
// //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;
// // }
// //}
//
}
// === End of loop over the entries
// === ============================================================================
TCanvas * can_mult = new TCanvas("mult","mult",0,0,2000,2000);
can_mult->Divide(3,2);
can_mult->cd(1); gPad->SetLogy(); h1_AnodeNumber->Draw();
can_mult->cd(2); gPad->SetLogy(); h1_MultInternal->Draw();
can_mult->cd(3); gPad->SetLogy(); h1_MultExternal->Draw();
can_mult->cd(4); gPad->SetLogy(); h1_MultTot->Draw();
can_mult->cd(5); h2_MultTot_vs_MultInternal->Draw("COL");
can_mult->cd(6); h2_MultTot_vs_MultExternal->Draw("COL");
TCanvas * can_QmaxInternal = new TCanvas("QmaxInternal","QmaxInternal",0,0,2000,1500);
can_QmaxInternal->cd();
h1_QmaxInternal->Draw();
// TFile * fsave = new TFile(Form("/media/audrey/DATA1/EPIC/simulation/out_K11_GEF240Pu.root",AlphaDecay_TimeWindow,fileno),"recreate");
// h1_Qmax->Write();
// fsave->Close();
}