Gabriele Cosmo
362 lines
11 KiB
C++
362 lines
11 KiB
C++
|
|
#include "TFile.h"
|
|
#include "TTree.h"
|
|
#include "TH1D.h"
|
|
#include "TH2D.h"
|
|
#include "TH1I.h"
|
|
#include "TF1.h"
|
|
#include "TProfile2D.h"
|
|
#include "TGraph.h"
|
|
#include "TGraph2D.h"
|
|
#include "TGraphErrors.h"
|
|
#include "TGaxis.h"
|
|
#include "TAxis.h"
|
|
#include "TMath.h"
|
|
#include "TRandom3.h"
|
|
#include "TCanvas.h"
|
|
#include "TLegend.h"
|
|
#include <iostream>
|
|
#include <fstream>
|
|
#include <vector>
|
|
#include "string.h"
|
|
#include <sstream>
|
|
|
|
#define pi 3.1415926535
|
|
|
|
using namespace std;
|
|
|
|
void ADXRD()
|
|
{
|
|
gROOT->Reset();
|
|
|
|
//----------------------------------input--------------------------------
|
|
//open a simulation result file
|
|
TFile* file = TFile::Open("output.root");
|
|
|
|
//histogram parameters
|
|
Double_t binXsize = 0.5; //mm
|
|
Double_t binXStart = -200.; //mm
|
|
Double_t binXEnd = -binXStart;
|
|
Int_t nbinsX = binXEnd*2./binXsize;
|
|
|
|
Double_t binYsize = 0.5; //mm
|
|
Double_t binYStart = -200.; //mm
|
|
Double_t binYEnd = binXEnd;
|
|
Int_t nbinsY = binYEnd*2./binYsize;
|
|
|
|
Double_t Nbin = 92;
|
|
Double_t thetaMin = 1; //deg
|
|
Double_t thetaMax = 27; //deg
|
|
|
|
//cuts
|
|
Double_t Emin = 0.; //keV
|
|
Double_t Emax = 140.; //keV
|
|
Double_t angCut = 0.; //deg
|
|
Bool_t IWantOnlyScatt = true;
|
|
Bool_t IWantOnlyRayleighScatt = false;
|
|
Bool_t IWantOnlyComptonScatt = false;
|
|
Bool_t IWantOnlyDiffraction = false;
|
|
|
|
//plot options
|
|
Bool_t IWantThetaDistribution = true;
|
|
Bool_t IWantEdistrib = true;
|
|
|
|
Bool_t IWantPlot2D = true;
|
|
Bool_t IWantProfiles = false;
|
|
|
|
Bool_t IWantBoxAnalysis = false;
|
|
Bool_t IWantPlotEtheta = true;
|
|
|
|
Int_t nbp = 2;
|
|
gStyle->SetOptStat(kFALSE);
|
|
//gStyle->SetPalette(52); //52->gray, 53->hot
|
|
|
|
//Energy Distribution Analysis
|
|
Double_t NbinE = 140;
|
|
Double_t EMin = 0; //keV
|
|
Double_t EMax = 140; //keV
|
|
Double_t Angle = 2.58; //deg
|
|
Double_t DeltaAngle = 2.; //deg
|
|
Bool_t ApplyThetaCut = false;
|
|
|
|
//export
|
|
Bool_t IwantExportScattering = false;
|
|
Char_t ExportScattFileName[128];
|
|
sprintf(ExportScattFileName,"scatt.txt");
|
|
|
|
Bool_t IwantExportVariables = false;
|
|
Char_t ExportVarFileName[128];
|
|
sprintf(ExportVarFileName,"var.txt");
|
|
|
|
Bool_t IwantExportImage = false;
|
|
Char_t ExportImageFileName[128];
|
|
sprintf(ExportImageFileName,"image.txt");
|
|
//-----------------------------------------------------------------------
|
|
|
|
//tree definition
|
|
Double_t x,y,energy;
|
|
Int_t kind,ID,nri,nci,ndi;
|
|
TTree* t1 = (TTree*)file->Get("part");
|
|
t1->SetBranchAddress("e",&energy);
|
|
t1->SetBranchAddress("posx",&x);
|
|
t1->SetBranchAddress("posy",&y);
|
|
t1->SetBranchAddress("type",&kind);
|
|
t1->SetBranchAddress("trackID",&ID);
|
|
t1->SetBranchAddress("NRi",&nri);
|
|
t1->SetBranchAddress("NCi",&nci);
|
|
t1->SetBranchAddress("NDi",&ndi);
|
|
Int_t N = t1->GetEntries();
|
|
const Int_t Nval = N;
|
|
|
|
//scatt cut
|
|
Char_t cutScatt[128];
|
|
if (IWantOnlyScatt) {
|
|
if (IWantOnlyRayleighScatt) {
|
|
sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi>=1 & NCi==0 & NDi==0",Emin,Emax);
|
|
} else if (IWantOnlyComptonScatt) {
|
|
sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi>=1 & NDi==0",Emin,Emax);
|
|
} else if (IWantOnlyDiffraction) {
|
|
sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi==0 & NDi>1",Emin,Emax);
|
|
} else {
|
|
sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & (NRi+NCi+NDi)>=1",Emin,Emax);
|
|
}
|
|
} else {
|
|
sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f",Emin,Emax);
|
|
}
|
|
|
|
|
|
//spatial distribution
|
|
TH2D* SpatialDistribution = new TH2D("SpatialDistribution", "Spatial Distribution", nbinsX, binXStart, binXEnd, nbinsY, binYStart, binYEnd);
|
|
t1->Project("SpatialDistribution","posy:posx",cutScatt);
|
|
if (IWantPlot2D) {
|
|
SpatialDistribution->SetXTitle("X (mm)");
|
|
SpatialDistribution->GetXaxis()->CenterTitle();
|
|
SpatialDistribution->GetXaxis()->SetTitleOffset(1.2);
|
|
SpatialDistribution->SetYTitle("Y (mm)");
|
|
SpatialDistribution->GetYaxis()->CenterTitle();
|
|
SpatialDistribution->GetYaxis()->SetTitleOffset(1.3);
|
|
TCanvas* c1 = new TCanvas("c1","",1000,1100);
|
|
c1->cd();
|
|
SpatialDistribution->SetContour(50);
|
|
SpatialDistribution->Draw("colz");
|
|
}
|
|
|
|
//profiles (projections transformed in profiles by myself)
|
|
if (IWantProfiles) {
|
|
TCanvas* c2 = new TCanvas("c2","",1200,900);
|
|
c2->Divide(2,1);
|
|
|
|
Int_t bcx = Int_t(nbinsX/2.);
|
|
Int_t bcy = Int_t(nbinsY/2.);
|
|
Double_t sf = 1./(2.*nbp+1.);
|
|
|
|
TH1D* profileX = SpatialDistribution->ProjectionX("px", bcx-nbp, bcx+nbp);
|
|
profileX->Scale(sf);
|
|
profileX->SetTitle("X profile");
|
|
c2->cd(1);
|
|
//profileX->Fit("gaus");
|
|
profileX->Draw();
|
|
|
|
TH1D* profileY = SpatialDistribution->ProjectionY("py", bcy-nbp, bcy+nbp);
|
|
profileY->Scale(sf);
|
|
profileY->SetTitle("Y profile");
|
|
c2->cd(2);
|
|
//profileY->Fit("gaus");
|
|
profileY->Draw();
|
|
|
|
cout << endl;
|
|
cout << "profileX FWHM: " << profileX->GetRMS()*2.35 << " mm" << endl;
|
|
cout << "profileY FWHM: " << profileY->GetRMS()*2.35 << " mm" << endl;
|
|
cout << endl;
|
|
}
|
|
|
|
|
|
//theta distribution
|
|
TH1D* thetaDistr = new TH1D("thetaDistr", "", Nbin, thetaMin, thetaMax);
|
|
t1->Project("thetaDistr","acos(momz)*180/acos(-1)",cutScatt);
|
|
|
|
Double_t Norig = thetaDistr->Integral();
|
|
|
|
Double_t val, angle;
|
|
for (Int_t i=1; i<=Nbin; i++) {
|
|
angle = thetaDistr->GetBinCenter(i)*acos(-1)/180;
|
|
val = thetaDistr->GetBinContent(i);
|
|
thetaDistr->SetBinContent(i,val/TMath::Sin(angle));
|
|
}
|
|
|
|
Double_t Ndiv = thetaDistr->Integral();
|
|
thetaDistr->Scale(Norig/Ndiv);
|
|
|
|
cout << endl << "total counts of thetaDistr: " << thetaDistr->Integral() << endl;
|
|
|
|
if (IWantThetaDistribution) {
|
|
thetaDistr->SetLineColor(kBlack);
|
|
thetaDistr->SetLineWidth(2);
|
|
thetaDistr->SetXTitle("#theta (degree)");
|
|
thetaDistr->GetXaxis()->CenterTitle();
|
|
thetaDistr->GetXaxis()->SetTitleOffset(1.1);
|
|
//thetaDistr->GetXaxis()->SetRangeUser(2., 12.);
|
|
thetaDistr->SetYTitle("Counts");
|
|
thetaDistr->GetYaxis()->CenterTitle();
|
|
thetaDistr->GetYaxis()->SetTitleOffset(1.2);
|
|
TCanvas* c3 = new TCanvas("c3","",1100,1100);
|
|
c3->cd();
|
|
thetaDistr->Draw("HIST");
|
|
}
|
|
|
|
|
|
//Energy Distribution Analysis
|
|
if (IWantEdistrib) {
|
|
//Theta cut
|
|
Char_t cutTheta[256];
|
|
|
|
if (ApplyThetaCut) {
|
|
sprintf(cutTheta,"type==0 & TMath::Abs(acos(momz)-%f) <= %f",Angle*0.017453293,DeltaAngle*0.017453293);
|
|
} else {
|
|
sprintf(cutTheta,"");
|
|
}
|
|
|
|
//Energy distribution of scattered photons
|
|
TH1D* edistrib = new TH1D("edistrib", "", NbinE, EMin, EMax);
|
|
t1->Project("edistrib", "e", cutTheta);
|
|
edistrib->SetLineColor(kRed);
|
|
edistrib->SetLineWidth(2);
|
|
edistrib->SetXTitle("E (keV)");
|
|
edistrib->GetXaxis()->CenterTitle();
|
|
edistrib->GetXaxis()->SetTitleOffset(1.1);
|
|
edistrib->SetYTitle("Counts (a. u.)");
|
|
edistrib->GetYaxis()->CenterTitle();
|
|
edistrib->GetYaxis()->SetTitleOffset(1.3);
|
|
edistrib->SetTitle("Energy spectrum of the particles impinging on the detector");
|
|
TCanvas* c4 = new TCanvas("c4","",1200,800);
|
|
c4->cd();
|
|
gStyle->SetOptStat(kFALSE);
|
|
edistrib->Draw();
|
|
|
|
Int_t Ntot = edistrib->Integral();
|
|
cout << "total counts of edistrib: " << Ntot << endl;
|
|
}
|
|
|
|
|
|
//Box Score Analysis
|
|
if (IWantBoxAnalysis) {
|
|
Int_t xBins = 40;
|
|
Int_t yBins = 40;
|
|
Int_t zBins = 20;
|
|
|
|
Double_t xlow = -25; //mm
|
|
Double_t xup = 25;
|
|
Double_t ylow = -25;
|
|
Double_t yup = 25;
|
|
Double_t zlow = -10;
|
|
Double_t zup = 10;
|
|
|
|
Double_t rngadd = 1.; //Box plot range margin (mm)
|
|
|
|
TH3D* XYZ = new TH3D("XYZ","Hot spots", xBins, xlow-rngadd, xup+rngadd, zBins, zlow-rngadd, zup+rngadd, yBins, ylow-rngadd, yup+rngadd);
|
|
t1->Project("XYZ", "posy:posz:posx",cutScatt);
|
|
XYZ->SetFillColor(2);
|
|
XYZ->SetXTitle("x (mm)");
|
|
XYZ->GetXaxis()->CenterTitle();
|
|
XYZ->GetXaxis()->SetTitleOffset(1.8);
|
|
XYZ->SetYTitle("z (mm)");
|
|
XYZ->GetYaxis()->CenterTitle();
|
|
XYZ->GetYaxis()->SetTitleOffset(2.4);
|
|
XYZ->SetZTitle("y (mm)");
|
|
XYZ->GetZaxis()->CenterTitle();
|
|
XYZ->GetZaxis()->SetTitleOffset(1.8);
|
|
gStyle->SetCanvasPreferGL(kTRUE);
|
|
TCanvas* c5 = new TCanvas("c5","Box Score Analysis",1200,800);
|
|
c5->cd();
|
|
XYZ->Draw("glbox");
|
|
}
|
|
|
|
|
|
//energy-angle correlation of impinging photons
|
|
if (IWantPlotEtheta) {
|
|
TH2D* Etheta = new TH2D("Etheta", "", NbinE, EMin, EMax, Nbin, thetaMin, thetaMax);
|
|
t1->Project("Etheta","acos(momz)*180/acos(-1):e");
|
|
Etheta->SetXTitle("E (keV)");
|
|
Etheta->GetXaxis()->CenterTitle();
|
|
Etheta->GetXaxis()->SetTitleOffset(1.2);
|
|
Etheta->SetYTitle("#theta (degree)");
|
|
Etheta->GetYaxis()->CenterTitle();
|
|
Etheta->GetYaxis()->SetTitleOffset(1.3);
|
|
TCanvas* c6 = new TCanvas("c6","",1000,1100);
|
|
c6->cd();
|
|
Etheta->SetContour(50);
|
|
Etheta->Draw("colz");
|
|
}
|
|
|
|
|
|
//export scattering data
|
|
if (IwantExportScattering) {
|
|
//open a txt file
|
|
ofstream f(ExportScattFileName);
|
|
if(!f) {
|
|
cout << "Error opening the file!";
|
|
return;
|
|
}
|
|
//variables extraction from tree
|
|
for (Int_t i=1; i<=Nbin; i++) {
|
|
Double_t ang = thetaDistr->GetBinCenter(i);
|
|
Double_t scatt = thetaDistr->GetBinContent(i);
|
|
if (ang >= angCut) {
|
|
f << ang << " " << scatt << endl;
|
|
}
|
|
}
|
|
//close the txt file
|
|
f.close();
|
|
cout << "writing the file with the scattering data successful!" << endl << endl;
|
|
}
|
|
|
|
|
|
//export (x,y) variables
|
|
if (IwantExportVariables) {
|
|
//open a txt file
|
|
ofstream f(ExportVarFileName);
|
|
if(!f) {
|
|
cout << "Error opening the file!";
|
|
return;
|
|
}
|
|
for (Int_t i=0; i<Nval; i++) {
|
|
t1->GetEntry(i);
|
|
if (IWantOnlyScatt) {
|
|
if (kind==0 && ID==1 && energy>=Emin && energy<=Emax && nri+nci+ndi>=1) {
|
|
f << x << " " << y << endl;
|
|
}
|
|
} else {
|
|
f << x << " " << y << endl;
|
|
}
|
|
}
|
|
//close the txt file
|
|
f.close();
|
|
cout << "writing the file with the (x,y) variables successful!" << endl << endl;
|
|
}
|
|
|
|
|
|
//export Image
|
|
if (IwantExportImage) {
|
|
//open a txt file
|
|
ofstream fout(ExportImageFileName);
|
|
if (!fout) {
|
|
cout << "Error opening the file!";
|
|
return;
|
|
}
|
|
//variables extraction from histogram
|
|
Double_t counts = 0.;
|
|
for (Int_t i=1; i<=nbinsY; i++) {
|
|
for (Int_t j=1; j<=nbinsX; j++) {
|
|
counts = SpatialDistribution->GetBinContent(j,i);
|
|
fout << counts << " ";
|
|
}
|
|
fout << endl;
|
|
}
|
|
//close the txt file
|
|
fout.close();
|
|
cout << "writing text image file successful!" << endl << endl;
|
|
}
|
|
|
|
}
|
|
|