exoticphysics/saxs/scattAnalysis.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 scattAnalysis()
0029 {
0030     gROOT->Reset();
0031
0032     //----------------------------------input--------------------------------

0033     //open a simulation result file

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

0037     Double_t thetalimit = 0.5;      //deg

0038     Double_t Nbin = 120.;
0039     Double_t thetaMin = 0.;         //deg

0040     Double_t thetaMax = 30.;        //deg

0041
0042     //options

0043     Bool_t IwantTotalScatterig = false;
0044     Bool_t IWantPlotComptonScattering = false;
0045     Bool_t IWantPlotPlotComparison = false;
0046     Bool_t IWantPlotProcessDistrib = false;
0047
0048     //export

0049     Bool_t IwantExportScattDistrib = false;
0050     Char_t ExportFileName[128];
0051     sprintf(ExportFileName,"SimResults.txt");
0052     //-----------------------------------------------------------------------

0053
0054     //definitions

0055     Char_t cutproc[256];
0056     Double_t val, angle;
0057
0058     //cut definition

0059     if (IwantTotalScatterig) {
0060         sprintf(cutproc,"(processID == 1 || processID == 2)");
0061     } else {
0062         sprintf(cutproc,"processID == 1");
0063     }
0064
0065     //define a tree from the ntuple contained in the input file

0066     TTree* scatt = (TTree*)file->Get("scatt");
0067     Int_t N = scatt->GetEntries();
0068
0069     //Rayleigh/Total scattering

0070     TH1D* thetaDistr = new TH1D("thetaDistr", "", Nbin, thetaMin, thetaMax);
0071     scatt->Project("thetaDistr", "theta", cutproc);
0072
0073     Double_t Norig = thetaDistr->Integral();
0074
0075     for (Int_t i=1; i<=Nbin; i++) {
0076         angle = thetaDistr->GetBinCenter(i);
0077         if (angle >= thetalimit) {
0078             val = thetaDistr->GetBinContent(i);
0079             thetaDistr->SetBinContent(i,val/TMath::Sin(angle*pi/180));
0080         }
0081     }
0082
0083     Double_t Ndiv = thetaDistr->Integral();
0084     thetaDistr->Scale(Norig/Ndiv);
0085
0086     cout << endl << "Scattering events: " << thetaDistr->Integral() << endl << endl;
0087
0088     thetaDistr->SetLineColor(kRed);
0089     thetaDistr->SetLineWidth(2);
0090     thetaDistr->SetXTitle("#theta (deg)");
0091     thetaDistr->GetXaxis()->CenterTitle();
0092     thetaDistr->GetXaxis()->SetTitleOffset(1.1);
0093     thetaDistr->GetXaxis()->SetRangeUser(0., thetaMax);
0094     thetaDistr->SetYTitle("Counts (a.u.)");
0095     thetaDistr->GetYaxis()->CenterTitle();
0096     thetaDistr->GetYaxis()->SetTitleOffset(1.5);
0097     thetaDistr->SetTitle("");
0098     TCanvas* c1 = new TCanvas("c1","",1200,800);
0099     c1->cd();
0100     gStyle->SetOptStat(kFALSE);
0101     thetaDistr->Draw("HIST");
0102
0103
0104     //Compton scattering

0105     if (IWantPlotComptonScattering) {
0106         TH1D* thetaDistrCompt = new TH1D("thetaDistrCompt", "", Nbin, thetaMin, thetaMax);
0107         scatt->Project("thetaDistrCompt", "theta", "processID == 2");
0108
0109         Double_t NorigC = thetaDistrCompt->Integral();
0110
0111         for (Int_t i=1; i<=Nbin; i++) {
0112             angle = thetaDistrCompt->GetBinCenter(i);
0113             if (angle >= thetalimit) {
0114                 val = thetaDistrCompt->GetBinContent(i);
0115                 thetaDistrCompt->SetBinContent(i,val/TMath::Sin(angle*pi/180));
0116             }
0117         }
0118
0119         Double_t NdivC = thetaDistrCompt->Integral();
0120         thetaDistrCompt->Scale(NorigC/NdivC);
0121
0122         thetaDistrCompt->SetLineColor(kBlue);
0123         thetaDistrCompt->SetLineWidth(2);
0124
0125         if (IWantPlotPlotComparison) {
0126             thetaDistrCompt->Draw("HIST SAME");
0127             TLegend* leg = new TLegend(0.7,0.75,0.85,0.85);
0128             leg->SetBorderSize(0);
0129             leg->AddEntry(thetaDistr,"Rayleigh","lp");
0130             leg->AddEntry(thetaDistrCompt,"Compton","lp");
0131             leg->Draw();
0132         } else {
0133             thetaDistrCompt->SetXTitle("#theta (deg)");
0134             thetaDistrCompt->GetXaxis()->CenterTitle();
0135             thetaDistrCompt->GetXaxis()->SetTitleOffset(1.1);
0136             thetaDistrCompt->GetXaxis()->SetRangeUser(0., thetaMax);
0137             thetaDistr->SetYTitle("Counts (a.u.)");
0138             thetaDistrCompt->GetYaxis()->CenterTitle();
0139             thetaDistrCompt->GetYaxis()->SetTitleOffset(1.5);
0140             thetaDistrCompt->SetTitle("");
0141             TCanvas* c2 = new TCanvas("c2","",1200,800);
0142             c2->cd();
0143             thetaDistrCompt->Draw("HIST SAME");
0144         }
0145
0146     }
0147
0148
0149     //process distribution (0->transport, 1->Rayleigh, 2->Compton, 3->Photoel,

0150     //                      4->Pair prod, 5->Nuclear, 6->Diffraction)

0151     if (IWantPlotProcessDistrib) {
0152         TH1D* procDistr = new TH1D("procDistr", "process distribution", 70, 0, 7);
0153         scatt->Project("procDistr", "processID");
0154         procDistr->SetLineColor(kBlack);
0155         procDistr->SetLineWidth(2);
0156         procDistr->SetXTitle("process");
0157         procDistr->GetXaxis()->CenterTitle();
0158         procDistr->GetXaxis()->SetTitleOffset(1.1);
0159         procDistr->SetYTitle("Counts");
0160         procDistr->GetYaxis()->CenterTitle();
0161         procDistr->GetYaxis()->SetTitleOffset(1.2);
0162         TCanvas* c3 = new TCanvas("c3","",1200,800);
0163         c3->cd();
0164         procDistr->Draw();
0165     }
0166
0167
0168     //export scattering data

0169     if (IwantExportScattDistrib) {
0170         ofstream f(ExportFileName);
0171         if (!f) {
0172             cout << "Error opening the file!";
0173             return;
0174         }
0175         for (Int_t i=1; i<=Nbin; i++) {
0176             angle = thetaDistr->GetBinCenter(i);
0177             if (angle >= thetalimit) {
0178                 val = thetaDistr->GetBinContent(i);
0179                 f << angle << " " << val << endl;
0180             }
0181         }
0182         f.close();
0183         cout << "writing successful!" << endl << endl;
0184     }
0185
0186 }
0187