exoticphysics/saxs/ADXRD.C

0001
0002 #include "TFile.h"
0003 #include "TTree.h"
0004 #include "TH1D.h"
0005 #include "TH2D.h"
0006 #include "TH1I.h"
0007 #include "TF1.h"
0008 #include "TProfile2D.h"
0009 #include "TGraph.h"
0010 #include "TGraph2D.h"
0011 #include "TGraphErrors.h"
0012 #include "TGaxis.h"
0013 #include "TAxis.h"
0014 #include "TMath.h"
0015 #include "TRandom3.h"
0016 #include "TCanvas.h"
0017 #include "TLegend.h"
0018 #include <iostream>
0019 #include <fstream>
0020 #include <vector>
0021 #include "string.h"
0022 #include <sstream>
0023
0024 #define pi 3.1415926535
0025
0026 using namespace std;
0027
0028 void ADXRD()
0029 {
0030     gROOT->Reset();
0031
0032     //----------------------------------input--------------------------------

0033     //open a simulation result file

0034     TFile* file = TFile::Open("output.root");
0035
0036     //histogram parameters

0037     Double_t binXsize = 0.5;        //mm

0038     Double_t binXStart = -200.;     //mm

0039     Double_t binXEnd = -binXStart;
0040     Int_t nbinsX = binXEnd*2./binXsize;
0041
0042     Double_t binYsize = 0.5;        //mm

0043     Double_t binYStart = -200.;     //mm

0044     Double_t binYEnd = binXEnd;
0045     Int_t nbinsY = binYEnd*2./binYsize;
0046
0047     Double_t Nbin = 92;
0048     Double_t thetaMin = 1;          //deg

0049     Double_t thetaMax = 27;         //deg

0050
0051     //cuts

0052     Double_t Emin = 0.;             //keV

0053     Double_t Emax = 140.;           //keV

0054     Double_t angCut = 0.;           //deg

0055     Bool_t IWantOnlyScatt = true;
0056     Bool_t IWantOnlyRayleighScatt = false;
0057     Bool_t IWantOnlyComptonScatt = false;
0058     Bool_t IWantOnlyDiffraction = false;
0059
0060     //plot options

0061     Bool_t IWantThetaDistribution = true;
0062     Bool_t IWantEdistrib = true;
0063
0064     Bool_t IWantPlot2D = true;
0065     Bool_t IWantProfiles = false;
0066
0067     Bool_t IWantBoxAnalysis = false;
0068     Bool_t IWantPlotEtheta = true;
0069
0070     Int_t nbp = 2;
0071     gStyle->SetOptStat(kFALSE);
0072     //gStyle->SetPalette(52);   //52->gray, 53->hot

0073
0074     //Energy Distribution Analysis

0075     Double_t NbinE = 140;
0076     Double_t EMin = 0;              //keV

0077     Double_t EMax = 140;            //keV

0078     Double_t Angle = 2.58;          //deg

0079     Double_t DeltaAngle = 2.;       //deg

0080     Bool_t ApplyThetaCut = false;
0081
0082     //export

0083     Bool_t IwantExportScattering = false;
0084     Char_t ExportScattFileName[128];
0085     sprintf(ExportScattFileName,"scatt.txt");
0086
0087     Bool_t IwantExportVariables = false;
0088     Char_t ExportVarFileName[128];
0089     sprintf(ExportVarFileName,"var.txt");
0090
0091     Bool_t IwantExportImage = false;
0092     Char_t ExportImageFileName[128];
0093     sprintf(ExportImageFileName,"image.txt");
0094     //-----------------------------------------------------------------------

0095
0096     //tree definition

0097     Double_t x,y,energy;
0098     Int_t kind,ID,nri,nci,ndi;
0099     TTree* t1 = (TTree*)file->Get("part");
0100     t1->SetBranchAddress("e",&energy);
0101     t1->SetBranchAddress("posx",&x);
0102     t1->SetBranchAddress("posy",&y);
0103     t1->SetBranchAddress("type",&kind);
0104     t1->SetBranchAddress("trackID",&ID);
0105     t1->SetBranchAddress("NRi",&nri);
0106     t1->SetBranchAddress("NCi",&nci);
0107     t1->SetBranchAddress("NDi",&ndi);
0108     Int_t N = t1->GetEntries();
0109     const Int_t Nval = N;
0110
0111     //scatt cut

0112     Char_t cutScatt[128];
0113     if (IWantOnlyScatt) {
0114         if (IWantOnlyRayleighScatt) {
0115             sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi>=1 & NCi==0 & NDi==0",Emin,Emax);
0116         } else if (IWantOnlyComptonScatt) {
0117             sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi>=1 & NDi==0",Emin,Emax);
0118         } else if (IWantOnlyDiffraction) {
0119             sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & NRi==0 & NCi==0 & NDi>1",Emin,Emax);
0120         } else {
0121             sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f & (NRi+NCi+NDi)>=1",Emin,Emax);
0122         }
0123     } else {
0124         sprintf(cutScatt,"type==0 & trackID==1 & e>=%f & e<=%f",Emin,Emax);
0125     }
0126
0127
0128     //spatial distribution

0129     TH2D* SpatialDistribution = new TH2D("SpatialDistribution", "Spatial Distribution", nbinsX, binXStart, binXEnd, nbinsY, binYStart, binYEnd);
0130     t1->Project("SpatialDistribution","posy:posx",cutScatt);
0131     if (IWantPlot2D) {
0132         SpatialDistribution->SetXTitle("X (mm)");
0133         SpatialDistribution->GetXaxis()->CenterTitle();
0134         SpatialDistribution->GetXaxis()->SetTitleOffset(1.2);
0135         SpatialDistribution->SetYTitle("Y (mm)");
0136         SpatialDistribution->GetYaxis()->CenterTitle();
0137         SpatialDistribution->GetYaxis()->SetTitleOffset(1.3);
0138         TCanvas* c1 = new TCanvas("c1","",1000,1100);
0139         c1->cd();
0140         SpatialDistribution->SetContour(50);
0141         SpatialDistribution->Draw("colz");
0142     }
0143
0144     //profiles (projections transformed in profiles by myself)

0145     if (IWantProfiles) {
0146         TCanvas* c2 = new TCanvas("c2","",1200,900);
0147         c2->Divide(2,1);
0148
0149         Int_t bcx = Int_t(nbinsX/2.);
0150         Int_t bcy = Int_t(nbinsY/2.);
0151         Double_t sf = 1./(2.*nbp+1.);
0152
0153         TH1D* profileX = SpatialDistribution->ProjectionX("px", bcx-nbp, bcx+nbp);
0154         profileX->Scale(sf);
0155         profileX->SetTitle("X profile");
0156         c2->cd(1);
0157         //profileX->Fit("gaus");

0158         profileX->Draw();
0159
0160         TH1D* profileY = SpatialDistribution->ProjectionY("py", bcy-nbp, bcy+nbp);
0161         profileY->Scale(sf);
0162         profileY->SetTitle("Y profile");
0163         c2->cd(2);
0164         //profileY->Fit("gaus");

0165         profileY->Draw();
0166
0167         cout << endl;
0168         cout << "profileX FWHM: " << profileX->GetRMS()*2.35 << " mm" << endl;
0169         cout << "profileY FWHM: " << profileY->GetRMS()*2.35 << " mm" << endl;
0170         cout << endl;
0171     }
0172
0173
0174     //theta distribution

0175     TH1D* thetaDistr = new TH1D("thetaDistr", "", Nbin, thetaMin, thetaMax);
0176     t1->Project("thetaDistr","acos(momz)*180/acos(-1)",cutScatt);
0177
0178     Double_t Norig = thetaDistr->Integral();
0179
0180     Double_t val, angle;
0181     for (Int_t i=1; i<=Nbin; i++) {
0182         angle = thetaDistr->GetBinCenter(i)*acos(-1)/180;
0183         val = thetaDistr->GetBinContent(i);
0184         thetaDistr->SetBinContent(i,val/TMath::Sin(angle));
0185     }
0186
0187     Double_t Ndiv = thetaDistr->Integral();
0188     thetaDistr->Scale(Norig/Ndiv);
0189
0190     cout << endl << "total counts of thetaDistr: " << thetaDistr->Integral() << endl;
0191
0192     if (IWantThetaDistribution) {
0193         thetaDistr->SetLineColor(kBlack);
0194         thetaDistr->SetLineWidth(2);
0195         thetaDistr->SetXTitle("#theta (degree)");
0196         thetaDistr->GetXaxis()->CenterTitle();
0197         thetaDistr->GetXaxis()->SetTitleOffset(1.1);
0198         //thetaDistr->GetXaxis()->SetRangeUser(2., 12.);

0199         thetaDistr->SetYTitle("Counts");
0200         thetaDistr->GetYaxis()->CenterTitle();
0201         thetaDistr->GetYaxis()->SetTitleOffset(1.2);
0202         TCanvas* c3 = new TCanvas("c3","",1100,1100);
0203         c3->cd();
0204         thetaDistr->Draw("HIST");
0205     }
0206
0207
0208     //Energy Distribution Analysis

0209     if (IWantEdistrib) {
0210         //Theta cut

0211         Char_t cutTheta[256];
0212
0213         if (ApplyThetaCut) {
0214             sprintf(cutTheta,"type==0 & TMath::Abs(acos(momz)-%f) <= %f",Angle*0.017453293,DeltaAngle*0.017453293);
0215         } else {
0216             sprintf(cutTheta,"");
0217         }
0218
0219         //Energy distribution of scattered photons

0220         TH1D* edistrib = new TH1D("edistrib", "", NbinE, EMin, EMax);
0221         t1->Project("edistrib", "e", cutTheta);
0222         edistrib->SetLineColor(kRed);
0223         edistrib->SetLineWidth(2);
0224         edistrib->SetXTitle("E (keV)");
0225         edistrib->GetXaxis()->CenterTitle();
0226         edistrib->GetXaxis()->SetTitleOffset(1.1);
0227         edistrib->SetYTitle("Counts (a. u.)");
0228         edistrib->GetYaxis()->CenterTitle();
0229         edistrib->GetYaxis()->SetTitleOffset(1.3);
0230         edistrib->SetTitle("Energy spectrum of the particles impinging on the detector");
0231         TCanvas* c4 = new TCanvas("c4","",1200,800);
0232         c4->cd();
0233         gStyle->SetOptStat(kFALSE);
0234         edistrib->Draw();
0235
0236         Int_t Ntot = edistrib->Integral();
0237         cout << "total counts of edistrib: " << Ntot << endl;
0238     }
0239
0240
0241     //Box Score Analysis

0242     if (IWantBoxAnalysis) {
0243         Int_t xBins = 40;
0244         Int_t yBins = 40;
0245         Int_t zBins = 20;
0246
0247         Double_t xlow = -25;            //mm

0248         Double_t xup = 25;
0249         Double_t ylow = -25;
0250         Double_t yup = 25;
0251         Double_t zlow = -10;
0252         Double_t zup = 10;
0253
0254         Double_t rngadd = 1.;           //Box plot range margin (mm)

0255
0256         TH3D* XYZ = new TH3D("XYZ","Hot spots", xBins, xlow-rngadd, xup+rngadd, zBins, zlow-rngadd, zup+rngadd, yBins, ylow-rngadd, yup+rngadd);
0257         t1->Project("XYZ", "posy:posz:posx",cutScatt);
0258         XYZ->SetFillColor(2);
0259         XYZ->SetXTitle("x (mm)");
0260         XYZ->GetXaxis()->CenterTitle();
0261         XYZ->GetXaxis()->SetTitleOffset(1.8);
0262         XYZ->SetYTitle("z (mm)");
0263         XYZ->GetYaxis()->CenterTitle();
0264         XYZ->GetYaxis()->SetTitleOffset(2.4);
0265         XYZ->SetZTitle("y (mm)");
0266         XYZ->GetZaxis()->CenterTitle();
0267         XYZ->GetZaxis()->SetTitleOffset(1.8);
0268         gStyle->SetCanvasPreferGL(kTRUE);
0269         TCanvas* c5 = new TCanvas("c5","Box Score Analysis",1200,800);
0270         c5->cd();
0271         XYZ->Draw("glbox");
0272     }
0273
0274
0275     //energy-angle correlation of impinging photons

0276     if (IWantPlotEtheta) {
0277         TH2D* Etheta = new TH2D("Etheta", "", NbinE, EMin, EMax, Nbin, thetaMin, thetaMax);
0278         t1->Project("Etheta","acos(momz)*180/acos(-1):e");
0279         Etheta->SetXTitle("E (keV)");
0280         Etheta->GetXaxis()->CenterTitle();
0281         Etheta->GetXaxis()->SetTitleOffset(1.2);
0282         Etheta->SetYTitle("#theta (degree)");
0283         Etheta->GetYaxis()->CenterTitle();
0284         Etheta->GetYaxis()->SetTitleOffset(1.3);
0285         TCanvas* c6 = new TCanvas("c6","",1000,1100);
0286         c6->cd();
0287         Etheta->SetContour(50);
0288         Etheta->Draw("colz");
0289     }
0290
0291
0292     //export scattering data

0293     if (IwantExportScattering) {
0294         //open a txt file

0295         ofstream f(ExportScattFileName);
0296         if(!f) {
0297             cout << "Error opening the file!";
0298             return;
0299         }
0300         //variables extraction from tree

0301         for (Int_t i=1; i<=Nbin; i++) {
0302             Double_t ang = thetaDistr->GetBinCenter(i);
0303             Double_t scatt = thetaDistr->GetBinContent(i);
0304             if (ang >= angCut) {
0305                 f << ang << " " << scatt << endl;
0306             }
0307         }
0308         //close the txt file

0309         f.close();
0310         cout << "writing the file with the scattering data successful!" << endl << endl;
0311     }
0312
0313
0314     //export (x,y) variables

0315     if (IwantExportVariables) {
0316         //open a txt file

0317         ofstream f(ExportVarFileName);
0318         if(!f) {
0319             cout << "Error opening the file!";
0320             return;
0321         }
0322         for (Int_t i=0; i<Nval; i++) {
0323             t1->GetEntry(i);
0324             if (IWantOnlyScatt) {
0325                 if (kind==0 && ID==1 && energy>=Emin && energy<=Emax && nri+nci+ndi>=1) {
0326                     f << x << " " << y << endl;
0327                 }
0328             } else {
0329                 f << x << " " << y << endl;
0330             }
0331         }
0332         //close the txt file

0333         f.close();
0334         cout << "writing the file with the (x,y) variables successful!" << endl << endl;
0335     }
0336
0337
0338     //export Image

0339     if (IwantExportImage) {
0340         //open a txt file

0341         ofstream fout(ExportImageFileName);
0342         if (!fout) {
0343             cout << "Error opening the file!";
0344             return;
0345         }
0346         //variables extraction from histogram

0347         Double_t counts = 0.;
0348         for (Int_t i=1; i<=nbinsY; i++) {
0349             for (Int_t j=1; j<=nbinsX; j++) {
0350                 counts = SpatialDistribution->GetBinContent(j,i);
0351                 fout << counts << " ";
0352             }
0353             fout << endl;
0354         }
0355         //close the txt file

0356         fout.close();
0357         cout << "writing text image file successful!" << endl << endl;
0358     }
0359
0360 }
0361