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 ;
int nbins_QvsT = 114;
int nbins_QvsT_external = 114;
int nbins_convol_external;
// --- 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);
TH1D * h1_QmaxExternal = new TH1D("QmaxExternalAnode_per_GEFevent","QmaxExternalAnode_per_GEFevent",100000,0,100000);
TCanvas * can1D_QvsT = new TCanvas("1D_QvsT","One Event",0,0,2500,2500);
can1D_QvsT->Divide(1,2);
#endif
// === ============================================================================
// === Loop over the entries
unsigned int nentries = t->GetEntries();
std::cout << " nentries = " << nentries << std::endl;
for (unsigned int Entry = 201; Entry < 1000; Entry++) {
for (unsigned int Entry = 0; Entry < nentries; 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) {
137
138
139
140
141
142
143
144
145
146
147
148
149
150
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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
// nbins_QvsT = 0;
// nbins_QvsT_external = 0;
// for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
// h1_AnodeNumber->Fill(m_Data->GetAnodeNbr(pmult));
// int nsteps = m_Data->GetAnodeNbr(pmult);
// if(m_Data->GetAnodeNbr(pmult)<100){
// mult_internal++;
// if(nsteps > 114) cout << "Attention !!! nsteps_internal = " << nsteps << " larger than maximum 114 steps " << endl;
// for(int smult=0; smult<nsteps; smult++){
// int nbins_current = ceil(m_Data->GetTimeCreationElectronsPerStep(pmult).at(smult) / time_bin_width_ns) + 114;
// if(nbins_QvsT < nbins_current) nbins_QvsT = nbins_current;
// }
// }
// else{
// mult_external++;
// if(nsteps > 114) cout << "Attention !!! nsteps_external = " << nsteps << " larger than maximum 114 steps " << endl;
// for(int smult=0; smult<nsteps; smult++){
// int nbins_current = ceil(m_Data->GetTimeCreationElectronsPerStep(pmult).at(smult) / time_bin_width_ns) + 114;
// if(nbins_QvsT_external < nbins_current) nbins_QvsT_external = nbins_current;
// }
// }
// }
// nbins_convol = nbins_QvsT + 2*nbins_gauss;
// vector<double> influence(nbins_convol,0);
// nbins_convol_external = nbins_QvsT_external + 2*nbins_gauss;
// vector<double> influence_external(nbins_convol_external,0);
//
//#if DEBUG
// TH1D * h1int_QvsT = new TH1D("InternalAnode_influence_vs_time_1D","InternalAnode_influence_vs_time_1D",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
// h1int_QvsT->GetXaxis()->SetTitle("Absolute time [ps] 200ps/bin");
// h1int_QvsT->SetLineColor(kBlue);
// h1int_QvsT->SetDirectory(0);
// TH1D * h1int_convolution_QvsT = new TH1D("InternalAnode_QvsT_convol","InternalAnode_QvsT_convol",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
// h1int_convolution_QvsT->SetLineColor(kRed);
//
// TH1D * h1ext_QvsT = new TH1D("ExternalAnode_influence_vs_time_1D","ExternalAnode_influence_vs_time_1D",nbins_convol_external,-2*range_gauss,nbins_QvsT_external*1000*time_bin_width_ns+2*range_gauss);
// h1ext_QvsT->GetXaxis()->SetTitle("Absolute time [ps] 200ps/bin");
// h1ext_QvsT->SetLineColor(kBlue);
// h1ext_QvsT->SetDirectory(0);
// TH1D * h1ext_convolution_QvsT = new TH1D("ExternalAnode_QvsT_convol","ExternalAnode_QvsT_convol",nbins_convol_external,-2*range_gauss,nbins_QvsT_external*1000*time_bin_width_ns+2*range_gauss);
// h1ext_convolution_QvsT->SetLineColor(kRed);
//#endif
//
// for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
// int nsteps = m_Data->GetNumberOfSteps(pmult);
// for(int smult=0; smult<nsteps; smult++){
// double t_Step = m_Data->GetTimeCreationElectronsPerStep(pmult).at(smult) ; // [ns NB. Vff >> Vdrift]
// int n_Step = m_Data->GetNumElectronsPerStep(pmult).at(smult);
// int b_Step = t_Step/time_bin_width_ns;
// for(int bin=smult; bin<114; bin++){
// if(m_Data->GetAnodeNbr(pmult)<100){
// if(nbins_QvsT < b_Step+bin-smult) {
// cout << "=== nbins_QvsT = " << nbins_QvsT << " but need to put data at index " << b_Step+bin-smult << endl;
// }
// influence[b_Step+bin-smult+nbins_gauss] += n_Step * 22. / 2500. ;
// }
// else{
// if(nbins_QvsT_external < b_Step+bin-smult) {
// cout << "=== nbins_QvsT_external = " << nbins_QvsT_external << " but need to put data at index " << b_Step+bin-smult << endl;
// }
// influence_external[b_Step+bin-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 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
// h1int_QvsT->SetBinContent(i+1,influence[i]);
// h1int_QvsT->SetBinError(i+1,0);
// h1int_convolution_QvsT->SetBinContent(i+1+0.5*nbins_gauss,convolution);
//#endif
// if (convolution>Qmax) Qmax=convolution;
// }
// h1_QmaxInternal->Fill(Qmax);
// influence.clear();
//
// // if data @ external proceed
// if(nbins_QvsT_external>0){
// Qmax = 0;
// for(int i=0; i<nbins_QvsT_external+nbins_gauss; i++){
// double convolution = 0;
// for (int j = 0; j < nbins_gauss; j++) convolution += influence_external[i+j] * gaussian_distribution[j];
//#if DEBUG
// h1ext_QvsT->SetBinContent(i+1,influence_external[i]);
// h1ext_QvsT->SetBinError(i+1,0);
// h1ext_convolution_QvsT->SetBinContent(i+1+0.5*nbins_gauss,convolution);
//#endif
// if (convolution>Qmax) Qmax=convolution;
// }
// h1_QmaxExternal->Fill(Qmax);
// }
// influence_external.clear();
//#if DEBUG
// can1D_QvsT->cd(1); h1int_QvsT->Draw(); h1int_convolution_QvsT->Draw("same");
// can1D_QvsT->cd(2); h1ext_QvsT->Draw(); h1ext_convolution_QvsT->Draw("same");
//
// can1D_QvsT->Update();
// can1D_QvsT->WaitPrimitive();
//
// h1int_QvsT->Reset(); h1int_convolution_QvsT->Reset(); delete h1int_convolution_QvsT;
// h1ext_QvsT->Reset(); h1ext_convolution_QvsT->Reset(); delete h1ext_convolution_QvsT;
//#endif
}// end of if else mult>1
// if mult==1, only internal anode
int nsteps = m_Data->GetNumberOfSteps(0);
if(nsteps > 114) cout << "Attention !!! nsteps = " << nsteps << " larger than maximum 114 steps " << endl;
nbins_QvsT = 114;
for(int smult=0; smult<nsteps; smult++){
int nbins_current = ceil(m_Data->GetTimeCreationElectronsPerStep(0).at(smult) / time_bin_width_ns) + 114;
if(nbins_QvsT < nbins_current) nbins_QvsT = nbins_current;
}
nbins_convol = nbins_QvsT + 2*nbins_gauss;
vector<double> influence(nbins_convol,0);
TH1D * h1_QvsT = new TH1D("Mult1_InternalAnode_influence_vs_time_1D","Mult1_InternalAnode_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("Mult1_InternalAnode_QvsT_convol","Mult1_InternalAnode_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] NB. Vff >> Vdrift
int n_Step = m_Data->GetNumElectronsPerStep(0).at(smult);
int b_Step = ceil(t_Step/time_bin_width_ns);
for(int bin=smult; bin<114; bin++){
if(nbins_QvsT < b_Step+bin-smult) {
cout << "--- nbins_QvsT = " << nbins_QvsT << " but need to put data at index " << b_Step+bin-smult << endl;
}
else
influence[b_Step+bin-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];
h1_QvsT->SetBinContent(i+1,influence[i]);
h1_QvsT->SetBinError(i+1,0);
h1_convolution_QvsT->SetBinContent(i+1+0.5*nbins_gauss,convolution);
if (convolution>Qmax) Qmax=convolution;
}
h1_AnodeNumber->Fill(m_Data->GetAnodeNbr(0));
if(m_Data->GetAnodeNbr(0)<100){
mult_internal++;
h1_QmaxInternal->Fill(Qmax);
}
else{
cout << "****** ******* ******** WARNING: MULT_TOT = 1 (" << m_Data->GetMultiplicity() << ") : INTERNAL ANODE ONLY BUT ANODE #" << m_Data->GetAnodeNbr(0) << endl;
mult_external++;
}
influence.clear();
can1D_QvsT->cd(1); h1_QvsT->Draw(); h1_convolution_QvsT->Draw("same");
can1D_QvsT->cd(2);
can1D_QvsT->Update();
can1D_QvsT->WaitPrimitive();
h1_QvsT->Reset(); delete h1_QvsT;
h1_convolution_QvsT->Reset(); delete h1_convolution_QvsT;
}// 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());
cout << endl;
// === 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_Qmax = new TCanvas("Qmax","Qmax",0,0,2000,1500);
can_Qmax->Divide(2,1);
can_Qmax->cd(1); h1_QmaxInternal->Draw();
can_Qmax->cd(1); h1_QmaxExternal->Draw();
// TFile * fsave = new TFile(Form("/media/audrey/DATA1/EPIC/simulation/out_K11_GEF240Pu.root",AlphaDecay_TimeWindow,fileno),"recreate");
// h1_Qmax->Write();
// fsave->Close();