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");
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);
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++) {
//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++){
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)
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
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);
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
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];
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();
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
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();