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"
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#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),300000,100000,400000);
#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);
// --- 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;
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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
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
TFile * fsave = new TFile(Form("out_K11_alpha240Pu_window%ius_%i.root",AlphaDecay_TimeWindow,fileno),"recreate");
h1_Qmax->Write();
fsave->Close();
}