Newer
Older
#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 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
{
// === ==========
// === input data
// === ==========
TChain * t = new TChain("SimulatedTree");
t->Add("/media/audrey/DATA1/EPIC/simulation/K11_GEF_240Pfns.root");
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
// === ==========
TH1D * h1_Qmax = new TH1D("Qmax_per_GEFevent","Qmax_per_GEFevent",100000,0,100000);
//// === DRAW
TCanvas * can1D_QvsT = new TCanvas("1D_QvsT","One Event",0,0,2500,2500);
// === ============================================================================
// === 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++) {
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
cout << endl ;
cout << "***********************************************" << endl;
cout << "@ Entry = " << Entry << ", mult = " << m_Data->GetMultiplicity() << endl;
#endif
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
if(m_Data->GetMultiplicity()==0) continue;
else if(m_Data->GetMultiplicity()!=1) {
cout << "Particle Multiplicity = " << m_Data->GetMultiplicity() << endl;
time_window_ns = 0;
// FIXME : LOOP
// loop over the particle multiplicity
for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
}
}// end of if else mult>1
else{
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;
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 steps up to anodes
}// end of loop over the steps
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 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>Qmax) Qmax=convolution;
h1_Qmax->Fill(Qmax);
influence.clear();
can1D_QvsT->cd();
h1_QvsT->Draw();
h1_convolution_QvsT->Draw("same");
can1D_QvsT->Update();
can1D_QvsT->WaitPrimitive();
h1_QvsT->Reset();
h1_convolution_QvsT->Reset();
}// end of else mult=1
// //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_Qmax = new TCanvas("Qmax","Qmax",0,0,2000,1500);
can_Qmax->cd();
h1_Qmax->Draw();
// TFile * fsave = new TFile(Form("/media/audrey/DATA1/EPIC/simulation/out_K11_alpha240Pu_window%ius_%i.root",AlphaDecay_TimeWindow,fileno),"recreate");
// h1_Qmax->Write();
// fsave->Close();