xray_SiliconPoreOptics/analysis/analysis.C

0001 #include "TROOT.h"
0002 #include "TTree.h"
0003 #include "TCanvas.h"
0004 #include "TBranch.h"
0005 #include <TFile.h>
0006 #include <TH1D.h>
0007 #include <iostream>
0008 #include <fstream>
0009 #include <string>
0010 #include <sstream>
0011 #include <vector>
0012 #include <stdio.h>
0013 #include <stdlib.h>
0014 #include <iomanip>
0015 #include <math.h>
0016 #include "TLegend.h"
0017 #include "TApplication.h"
0018 #include "TString.h"
0019 #include <algorithm>
0020 #include "TPie.h"
0021
0022 using namespace std;
0023 TApplication myapp("app",NULL,NULL);
0024
0025 const double pi = 3.14159;
0026 double to_deg(double angle)
0027 {
0028   return angle * 180 / pi;
0029 }
0030
0031 //////////////////////////////////////////////////////
0032
0033 int main(int argc, char *argv[]){
0034   // input = ASCII data filename
0035   string ifilename = "";
0036
0037   if(argc == 1)
0038   {
0039     cout << "No input file selected. Please use: ./read_tree_spectrum [input_file.root]" << endl;
0040     return 0;
0041   }
0042   else
0043   {
0044     ifilename = argv[1];
0045   }
0046   TString ifilenameforroot = ifilename;
0047
0048   // create the file, the Tree and a few branches
0049   TFile *f = new TFile(ifilenameforroot, "READ");
0050   TTree *mytree = (TTree*)f->Get("XraySPO");  // name of the ntuple in the rootfile
0051
0052   // create variables to store column values
0053   int eventID = 0;
0054   int trackID = 0;
0055   double x = 0;
0056   double y = 0;
0057   double z = 0;
0058   double theta = 0;
0059   double phi = 0;
0060   int parentID = 0;
0061   char vol_name[500];
0062   char proc_name[500];
0063   int num_reflections = 0;
0064
0065   mytree->SetBranchAddress("eventID", &eventID);
0066   mytree->SetBranchAddress("vol_name", &vol_name);
0067   mytree->SetBranchAddress("trackID", &trackID);
0068   mytree->SetBranchAddress("x", &x);
0069   mytree->SetBranchAddress("y", &y);
0070   mytree->SetBranchAddress("z", &z);
0071   mytree->SetBranchAddress("theta", &theta);
0072   mytree->SetBranchAddress("phi", &phi);
0073   mytree->SetBranchAddress("proc_name", &proc_name);
0074   mytree->SetBranchAddress("parentID", &parentID);
0075   mytree->SetBranchAddress("num_reflections", &num_reflections);
0076
0077   int bin_theta = 50;
0078   int bin_phi = 100;
0079
0080   int min_val_theta = 0;
0081   int max_val_theta = 5;
0082   int min_val_phi = -100;
0083   int max_val_phi = 100;
0084
0085   double bin_width_theta = (max_val_theta-min_val_theta)/(double)bin_theta;
0086   double bin_width_phi = (max_val_phi-min_val_phi)/(double)bin_phi;
0087
0088   TH1D *angle_theta = new TH1D ("theta", "Efficiency - theta [deg]", bin_theta, min_val_theta, max_val_theta);
0089   TH1D *angle_phi = new TH1D ("phi", "Efficiency - phi [deg]", bin_phi, min_val_phi, max_val_phi);
0090
0091   TH1I *reflections = new TH1I ("reflections", "Number of reflections", 8, -1.5, 6.5);
0092
0093   // Read all entries and fill the histograms
0094   Long64_t nentries = mytree->GetEntries();
0095   cout << "Reading " << nentries << " from the TTree." << endl;
0096
0097   int trackID_old = 0;
0098
0099   bool passed_entrance = false;
0100   bool passed_exit = false;
0101   int eventID_old = 0;
0102
0103   int check = -999;
0104
0105   int count_enter = 0;
0106   int count_exit = 0;
0107
0108   double theta_exit = 0 ;
0109   double phi_exit = 0;
0110
0111   // std::string new_particle_event = "";
0112   for (Long64_t i=0; i<nentries; i++)
0113   {
0114     mytree->GetEntry(i);
0115     if (eventID != eventID_old || eventID == 0)
0116     {
0117       passed_entrance = false;
0118       passed_exit = false;
0119
0120       if ((string)vol_name == "pDummyEntrance")
0121       {
0122         passed_entrance = true;
0123         passed_exit = false;
0124         count_enter += 1;
0125       }
0126     }
0127     else
0128     {
0129       if ((string)vol_name == "pDummyExit" && passed_entrance == true)
0130       {
0131         // Plot the angle at the
0132         theta_exit = theta;
0133         phi_exit = phi;
0134         passed_exit = true;
0135       }
0136
0137       if (passed_exit)
0138       {
0139         if ((string)vol_name == "pDummySphere")
0140         {
0141           passed_entrance = false;
0142           passed_exit = false;
0143
0144           if (num_reflections < 1) num_reflections = -1;
0145           reflections->Fill(num_reflections);
0146
0147           // Check if the eventID has not already been counted before
0148           if (eventID == check) cout << "Evento uguale! " << endl;
0149
0150           theta_exit = pi-theta_exit;
0151           if (phi_exit <= 0) phi_exit = -pi-phi_exit;
0152           else phi_exit = pi-phi_exit;
0153
0154           if (theta_exit >= min_val_theta && theta_exit <= max_val_theta)
0155           {
0156             if (phi_exit >= min_val_phi && phi_exit <= max_val_phi)
0157             {
0158               count_exit+=1;
0159               angle_theta->Fill(to_deg(theta_exit));
0160               angle_phi->Fill(to_deg(phi_exit));
0161             }
0162           }
0163           check = eventID;
0164         }
0165       }
0166     }
0167     eventID_old = eventID;
0168   }
0169
0170   // pie chart
0171   Float_t val0 = reflections->GetBinContent(1);  // number of entries with 0 reflections
0172   Float_t val2 = reflections->GetBinContent(3);  // number of entries with 1 reflection
0173   Float_t val3 = reflections->GetBinContent(4);
0174   Float_t val4 = reflections->GetBinContent(5);
0175   Float_t val5 = reflections->GetBinContent(6) + reflections->GetBinContent(7); // value for 4+ reflections
0176   cout << "Values are: " << val0 << " " << val2 << " " << val3 << " " << val4 << " " << val5 << " " << endl;
0177
0178
0179   Float_t vals[] = {val0,val2,val3,val4,val5};
0180   Int_t colors[] = {0,1,2,3,4};
0181   Int_t nvals = sizeof(vals)/sizeof(vals[0]);
0182
0183   TCanvas *cpie = new TCanvas("cpie", "cpie", 0, 0, 900, 900);
0184   TPie *pie1 = new TPie("pie1", "100 keV, 1deg., SS model",nvals,vals,colors);
0185   pie1->SetLabelsOffset(.01);
0186   pie1->SetRadius(.2);
0187
0188   pie1->SetEntryLabel(0, "0");
0189   pie1->SetEntryLabel(1, "1");
0190   pie1->SetEntryLabel(2, "2");
0191   pie1->SetEntryLabel(3, "3");
0192   pie1->SetEntryLabel(4, ">3");
0193   pie1->SetLabelFormat("%txt (%perc)");
0194   pie1->Draw();
0195
0196   // Normalization
0197   cout << "Count particle: <entered: " << count_enter << "> <exited: " << count_exit << ">" << endl;
0198   angle_theta->Scale(1.0/((double)count_enter));
0199   angle_phi->Scale(1.0/((double)count_enter));
0200   angle_theta->Scale(1.0/(bin_width_theta));
0201   angle_phi->Scale(1.0/(bin_width_phi));
0202
0203   TCanvas *c1 = new TCanvas("c1", "c1", 0, 50, 900, 900);
0204   int c1_id = c1->GetCanvasID();
0205   cout << "The first canvas ID is: " << c1_id << endl;
0206   c1->SetLogy();
0207   angle_theta->SetFillColor(kBlack);
0208   angle_theta->SetMarkerColor(kBlack);
0209   angle_theta->Draw("");
0210   angle_theta->GetXaxis()->SetTitle("Theta [deg]");
0211   angle_theta->GetYaxis()->SetTitle("Normalized efficiency [deg^-1]");
0212   angle_theta->GetXaxis()->SetTitleSize(0.03);
0213   angle_theta->GetYaxis()->SetTitleSize(0.03);
0214   cout << "integral of angle_theta from histo:" << angle_theta->Integral() << " and from fraction: " << count_exit/count_enter << endl;
0215
0216   TCanvas *c2 = new TCanvas("c2", "c2", 0, 100, 900, 900);
0217   int c2_id = c2->GetCanvasID();
0218   cout << "The second canvas ID is: " << c2_id << endl;
0219   c2->SetLogy();
0220   angle_phi->SetFillColor(kBlack);
0221   angle_phi->SetMarkerColor(kBlack);
0222   angle_phi->Draw("");
0223   angle_phi->GetXaxis()->SetTitle("Phi [deg]");
0224   angle_phi->GetYaxis()->SetTitle("Normalized efficiency [deg^-1]");
0225   angle_phi->GetXaxis()->SetTitleSize(0.03);
0226   angle_phi->GetYaxis()->SetTitleSize(0.03);
0227   cout << "integral of angle_phi is:" << angle_theta->Integral() << " and from fraction: " << count_exit/count_enter << endl;
0228
0229   myapp.Run(true);
0230
0231   f->Close();
0232   return 0;
0233 }