1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
| #include <iomanip>
#include <time.h>
#include <signal.h>
#include <math.h>
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THStack.h>
#include <TLegend.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TArrayD.h>
#include <TMatrixD.h>
#include <TVectorD.h>
#include <TBranch.h>
#include <TColor.h>
#include <TRandom3.h>
#include <TFile.h>
#include <TApplication.h>
#include <TProfile2D.h>
#include <TRandom3.h>
using namespace std;
int main(int argc, char **argv)
{
//#########################################################
// root handling
// ./main -b : batch mode (no display of the figures)
//#########################################################
TApplication theApp("App", &argc, argv);
// Factor of statistics reduction (to speed up the program execution)
Double_t StatFactor = 1;
//root colz style (smooth 2D palette)
const Int_t NRGBs = 5;
const Int_t NCont = 255;
Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
gStyle->SetNumberContours(NCont);
//General settings
static const Int_t radius = 100; //phantom (cylinder) radius in mm
Double_t I0 = 500.; //photons number per bin for each projection
Double_t pixelSize = 1.; //in mm
static const Int_t Nproj = 180; //projection number
Double_t theta[Nproj] = {0.};
Double_t AngularStep = 1.; // Angular step bewteen the radiography
for (int div=0; div<Nproj; div++) theta[div] = AngularStep*div*TMath::Pi()/180.; //RunID = projection angle
// Int_t iter = 3; //iteration number
static const Int_t idiam = 2*radius;
// Histograms
// TH2F *hEnergyZ = new TH2F("hEnergyZ","hEnergyZ",200,-100,100,100,0,100);
TH2F *hYZ = new TH2F("hYZ","hYZ",200,-100,100,21,-10,10);
//CT Tables x=angle, y=position on projection axis
Double_t TabCT[Nproj][idiam] = {{0.,0.}};
//Get tree results (Gate output)
TFile *f = new TFile("RootFile/PhaseSpace.root","read");
TTree *tree = (TTree *) f->Get("PhaseSpace");
//declare tree variables
Float_t Y, Z, Ekine; // position in mm and Ekine in eV
Int_t RunID;
//Get tree variables
tree->SetBranchAddress("Z",&Z);
tree->SetBranchAddress("Y",&Y);
tree->SetBranchAddress("Ekine",&Ekine);
tree->SetBranchAddress("RunID", &RunID );
//loop over tree entries
Int_t nentries = (Int_t)tree->GetEntries();
printf("Number of entries = %d\n",nentries);
cout << "Number of read entries: " << nentries/StatFactor << endl;
for (Int_t n=0; n < nentries/StatFactor; n++)
{
//Get tree entry n
tree->GetEntry(n);
// Energy E(keV) and position (mm) on the detector
// hEnergyZ->Fill(Z,Ekine*1000.,1.);
hYZ->Fill(Y,Z,1.);
}
gStyle->SetOptStat(0000);
TCanvas *c0 = new TCanvas("hYZ","hYZ",500,500);
hYZ->Draw("colz");
hYZ->GetXaxis()->SetTitle("Position Y (mm)");
hYZ->GetYaxis()->SetTitle("Position Z (mm)");
c0->SetLogz();
// TCanvas *c1 = new TCanvas("hEnergyZ","hEnergyZ",500,500);
// hEnergyZ->Draw("colz");
// hEnergyZ->GetXaxis()->SetTitle("Position Z (mm)");
// hEnergyZ->GetYaxis()->SetTitle("Energie (keV)");
// c1->SetLogz();
// Save the figures
c0-> SaveAs("Figures/hZY.pdf");
// c1-> SaveAs("Figures/hEnergyZ.pdf");
// Display of the figures if "interactive" mode (not batch)
if(!gROOT->IsBatch()) theApp.Run();
return 0;
} |
Partager