Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-13 08:24:57 

0001 // -------------------------------------------------------------------
0002 // -------------------------------------------------------------------
0003 //
0004 // *********************************************************************
0005 // To execute this macro under ROOT,
0006 //   1 - launch ROOT (usually type 'root' at your machine's prompt)
0007 //   2 - type '.X plot.C' at the ROOT session prompt
0008 // Written by S. Incerti, 25/01/2025
0009 // *********************************************************************
0010 {
0011 gROOT->Reset();
0012 gROOT->SetStyle("Plain");
0013 gStyle->SetOptStat(0000);
0014 gStyle->SetPalette(1);
0015 
0016 auto c1 = new TCanvas ("c1","",20,20,1200,900);
0017 c1->Divide(4,3);
0018 
0019 //------------------------------
0020 // Original phantom file view
0021 //------------------------------
0022 
0023 FILE * fp = fopen("phantoms/phantom.dat","r");
0024 
0025 Double_t X, Y, Z, mat, tmp;
0026 char unit[100];
0027 Double_t voxelSizeX, voxelSizeY, voxelSizeZ;
0028 Long_t numberVoxTot, numberVoxRed, numberVoxGreen, numberVoxBlue;
0029 
0030 TNtuple *ntuplePhantom = new TNtuple("PHANTOM","ntuple","X:Y:Z:mat");
0031 
0032 Long_t nlines=0;
0033 Long_t ncols=0;
0034 
0035 while (1)
0036    {
0037       if ( nlines == 0 ) ncols = fscanf(fp,"%ld %ld %ld %ld",&numberVoxTot,&numberVoxRed,&numberVoxGreen,&numberVoxBlue);
0038       if ( nlines == 1 ) ncols = fscanf(fp,"%lf %lf %lf %s",&tmp,&tmp,&tmp,unit);
0039       if ( nlines == 2 ) ncols = fscanf(fp,"%lf %lf %lf %s",&voxelSizeX,&voxelSizeY,&voxelSizeZ, unit);
0040       if ( nlines >= 3 ) ncols = fscanf(fp,"%lf %lf %lf %lf", &X, &Y, &Z, &mat);
0041       //cout << X << " " << Y << " " << Z << " " << mat << endl;
0042       if (ncols < 0) break;
0043       ntuplePhantom->Fill(X,Y,Z,mat);
0044       nlines++;
0045    }
0046 fclose(fp);
0047 
0048 c1->cd(1);
0049 
0050 ntuplePhantom->SetMarkerColor(1);
0051 ntuplePhantom->Draw("Y:X");
0052 // RED
0053 ntuplePhantom->SetMarkerColor(2);
0054 ntuplePhantom->Draw("Y:X","mat==1","same");
0055 // GREEN
0056 ntuplePhantom->SetMarkerColor(3);
0057 ntuplePhantom->Draw("Y:X","mat==2","same");
0058 // BLUE
0059 ntuplePhantom->SetMarkerColor(4);
0060 ntuplePhantom->Draw("Y:X","mat==3","same");
0061 //
0062 TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
0063 htemp->GetXaxis()->SetTitle("X (microns)");
0064 htemp->GetYaxis()->SetTitle("Y (mirons)");
0065 htemp->GetXaxis()->SetLabelSize(0.025);
0066 htemp->GetYaxis()->SetLabelSize(0.025);
0067 htemp->GetXaxis()->SetTitleSize(0.035);
0068 htemp->GetYaxis()->SetTitleSize(0.035);
0069 htemp->GetXaxis()->SetTitleOffset(1.4);
0070 htemp->GetYaxis()->SetTitleOffset(1.4);
0071 htemp->SetTitle("RGB phantom YX view");
0072 
0073 c1->cd(5);
0074 
0075 ntuplePhantom->SetMarkerColor(1);
0076 ntuplePhantom->Draw("Y:Z");
0077 // RED
0078 ntuplePhantom->SetMarkerColor(2);
0079 ntuplePhantom->Draw("Y:Z","mat==1","same");
0080 // GREEN
0081 ntuplePhantom->SetMarkerColor(3);
0082 ntuplePhantom->Draw("Y:Z","mat==2","same");
0083 // BLUE
0084 ntuplePhantom->SetMarkerColor(4);
0085 ntuplePhantom->Draw("Y:Z","mat==3","same");
0086 //
0087 TH2F *htempBis = (TH2F*)gPad->GetPrimitive("htemp");
0088 htempBis->GetXaxis()->SetTitle("Z (microns)");
0089 htempBis->GetYaxis()->SetTitle("Y (mirons)");
0090 htempBis->GetXaxis()->SetLabelSize(0.025);
0091 htempBis->GetYaxis()->SetLabelSize(0.025);
0092 htempBis->GetXaxis()->SetTitleSize(0.035);
0093 htempBis->GetYaxis()->SetTitleSize(0.035);
0094 htempBis->GetXaxis()->SetTitleOffset(1.4);
0095 htempBis->GetYaxis()->SetTitleOffset(1.4);
0096 htempBis->SetTitle("RGB phantom YZ view");
0097 
0098 c1->cd(9);
0099 
0100 ntuplePhantom->SetMarkerColor(1);
0101 ntuplePhantom->Draw("X:Z");
0102 // RED
0103 ntuplePhantom->SetMarkerColor(2);
0104 ntuplePhantom->Draw("X:Z","mat==1","same");
0105 // GREEN
0106 ntuplePhantom->SetMarkerColor(3);
0107 ntuplePhantom->Draw("X:Z","mat==2","same");
0108 // BLUE
0109 ntuplePhantom->SetMarkerColor(4);
0110 ntuplePhantom->Draw("X:Z","mat==3","same");
0111 //
0112 TH2F *htempTer = (TH2F*)gPad->GetPrimitive("htemp");
0113 htempTer->GetXaxis()->SetTitle("Z (microns)");
0114 htempTer->GetYaxis()->SetTitle("X (mirons)");
0115 htempTer->GetXaxis()->SetLabelSize(0.025);
0116 htempTer->GetYaxis()->SetLabelSize(0.025);
0117 htempTer->GetXaxis()->SetTitleSize(0.035);
0118 htempTer->GetYaxis()->SetTitleSize(0.035);
0119 htempTer->GetXaxis()->SetTitleOffset(1.4);
0120 htempTer->GetYaxis()->SetTitleOffset(1.4);
0121 htempTer->SetTitle("RGB phantom XZ view");
0122 
0123 //------------------
0124 // Read ROOT file
0125 //------------------
0126 
0127 TFile *f = new TFile ("phantom.root");
0128 
0129 TNtuple* ntuple1;
0130 TNtuple* ntuple2;
0131 TNtuple* ntuple3;
0132 
0133 ntuple1 = (TNtuple*)f->Get("ntuple1");
0134 ntuple2 = (TNtuple*)f->Get("ntuple2");
0135 ntuple3 = (TNtuple*)f->Get("ntuple3");
0136 
0137 Double_t x, y, z, energy, dose;
0138 Int_t voxelID;
0139 
0140 ntuple1->SetBranchAddress("x",&x);
0141 ntuple1->SetBranchAddress("y",&y);
0142 ntuple1->SetBranchAddress("z",&z);
0143 ntuple1->SetBranchAddress("energy",&energy);
0144 ntuple1->SetBranchAddress("dose",&dose);
0145 ntuple1->SetBranchAddress("voxelID",&voxelID);
0146 
0147 ntuple2->SetBranchAddress("x",&x);
0148 ntuple2->SetBranchAddress("y",&y);
0149 ntuple2->SetBranchAddress("z",&z);
0150 ntuple2->SetBranchAddress("energy",&energy);
0151 ntuple2->SetBranchAddress("dose",&dose);
0152 ntuple2->SetBranchAddress("voxelID",&voxelID);
0153 
0154 ntuple3->SetBranchAddress("x",&x);
0155 ntuple3->SetBranchAddress("y",&y);
0156 ntuple3->SetBranchAddress("z",&z);
0157 ntuple3->SetBranchAddress("energy",&energy);
0158 ntuple3->SetBranchAddress("dose",&dose);
0159 ntuple3->SetBranchAddress("voxelID",&voxelID);
0160 
0161 //---------------------------------
0162 // Absorbed energy distributions
0163 //---------------------------------
0164 
0165 c1->cd(2);
0166 gPad->SetLogy();
0167 ntuple1->Draw("energy","energy>0");
0168 TH1F *htemp2 = (TH1F*)gPad->GetPrimitive("htemp");
0169 htemp2->GetXaxis()->SetTitle("Energy (keV)");
0170 htemp2->GetXaxis()->SetLabelSize(0.025);
0171 htemp2->GetXaxis()->SetTitleSize(0.035);
0172 htemp2->GetXaxis()->SetTitleOffset(1.4);
0173 htemp2->SetTitle("RED voxel energy");
0174 htemp2->SetFillStyle(1001);
0175 htemp2->SetFillColor(2);
0176 
0177 c1->cd(6);
0178 gPad->SetLogy();
0179 ntuple2->Draw("energy","energy>0");
0180 TH1F *htemp3 = (TH1F*)gPad->GetPrimitive("htemp");
0181 htemp3->GetXaxis()->SetTitle("Energy (keV)");
0182 htemp3->GetXaxis()->SetLabelSize(0.025);
0183 htemp3->GetXaxis()->SetTitleSize(0.035);
0184 htemp3->GetXaxis()->SetTitleOffset(1.4);
0185 htemp3->SetTitle("GREEN voxel energy");
0186 htemp3->SetFillStyle(1001);
0187 htemp3->SetFillColor(3);
0188 
0189 c1->cd(10);
0190 gPad->SetLogy();
0191 ntuple3->Draw("energy","energy>0");
0192 TH1F *htemp4 = (TH1F*)gPad->GetPrimitive("htemp");
0193 htemp4->GetXaxis()->SetTitle("Energy (keV)");
0194 htemp4->GetXaxis()->SetLabelSize(0.025);
0195 htemp4->GetXaxis()->SetTitleSize(0.035);
0196 htemp4->GetXaxis()->SetTitleOffset(1.4);
0197 htemp4->SetTitle("BLUE voxel energy");
0198 htemp4->SetFillStyle(1001);
0199 htemp4->SetFillColor(4);
0200 
0201 //------------------------------
0202 // Map of energy distribution
0203 //------------------------------
0204 
0205 c1->cd(3);
0206 TH2F *histNrjRed = new TH2F("histNrjRed","histNrjRed",100,0,800,100,0,800);
0207 ntuple1->Draw("y:x>>histNrjRed","energy","contz");
0208 gPad->SetLogz();
0209 histNrjRed->Draw("contz");
0210 histNrjRed->GetXaxis()->SetTitle("X (microns)");
0211 histNrjRed->GetYaxis()->SetTitle("Y (mirons)");
0212 histNrjRed->GetZaxis()->SetTitle("Energy (keV)");
0213 histNrjRed->GetXaxis()->SetLabelSize(0.025);
0214 histNrjRed->GetYaxis()->SetLabelSize(0.025);
0215 histNrjRed->GetZaxis()->SetLabelSize(0.025);
0216 histNrjRed->GetXaxis()->SetTitleSize(0.035);
0217 histNrjRed->GetYaxis()->SetTitleSize(0.035);
0218 histNrjRed->GetZaxis()->SetTitleSize(0.035);
0219 histNrjRed->GetXaxis()->SetTitleOffset(1.4);
0220 histNrjRed->GetYaxis()->SetTitleOffset(1.4);
0221 histNrjRed->GetZaxis()->SetTitleOffset(.6);
0222 histNrjRed->SetTitle("Energy map for RED voxels");
0223 
0224 c1->cd(7);
0225 TH2F *histNrjGreen = new TH2F("histNrjGreen","histNrjGreen",100,0,800,100,0,800);
0226 ntuple2->Draw("y:x>>histNrjGreen","energy","contz");
0227 gPad->SetLogz();
0228 histNrjGreen->Draw("contz");
0229 histNrjGreen->GetXaxis()->SetTitle("X (microns)");
0230 histNrjGreen->GetYaxis()->SetTitle("Y (mirons)");
0231 histNrjGreen->GetZaxis()->SetTitle("Energy (keV)");
0232 histNrjGreen->GetXaxis()->SetLabelSize(0.025);
0233 histNrjGreen->GetYaxis()->SetLabelSize(0.025);
0234 histNrjGreen->GetZaxis()->SetLabelSize(0.025);
0235 histNrjGreen->GetXaxis()->SetTitleSize(0.035);
0236 histNrjGreen->GetYaxis()->SetTitleSize(0.035);
0237 histNrjGreen->GetZaxis()->SetTitleSize(0.035);
0238 histNrjGreen->GetXaxis()->SetTitleOffset(1.4);
0239 histNrjGreen->GetYaxis()->SetTitleOffset(1.4);
0240 histNrjGreen->GetZaxis()->SetTitleOffset(.6);
0241 histNrjGreen->SetTitle("Energy map for GREEN voxels");
0242 
0243 c1->cd(11);
0244 TH2F *histNrjBlue = new TH2F("histNrjBlue","histNrjBlue",100,0,800,100,0,800);
0245 ntuple3->Draw("y:x>>histNrjBlue","energy","contz");
0246 gPad->SetLogz();
0247 histNrjBlue->Draw("contz");
0248 histNrjBlue->GetXaxis()->SetTitle("X (microns)");
0249 histNrjBlue->GetYaxis()->SetTitle("Y (mirons)");
0250 histNrjBlue->GetZaxis()->SetTitle("Energy (keV)");
0251 histNrjBlue->GetXaxis()->SetLabelSize(0.025);
0252 histNrjBlue->GetYaxis()->SetLabelSize(0.025);
0253 histNrjBlue->GetZaxis()->SetLabelSize(0.025);
0254 histNrjBlue->GetXaxis()->SetTitleSize(0.035);
0255 histNrjBlue->GetYaxis()->SetTitleSize(0.035);
0256 histNrjBlue->GetZaxis()->SetTitleSize(0.035);
0257 histNrjBlue->GetXaxis()->SetTitleOffset(1.4);
0258 histNrjBlue->GetYaxis()->SetTitleOffset(1.4);
0259 histNrjBlue->GetZaxis()->SetTitleOffset(.6);
0260 histNrjBlue->SetTitle("Energy map for BLUE voxels");
0261 
0262 //----------------------------
0263 // Map of dose distribution
0264 //----------------------------
0265 
0266 c1->cd(4);
0267 TH2F *histDoseRed = new TH2F("histDoseRed","histDoseRed",100,0,800,100,0,800);
0268 // WARNING : dose scaling to mGy
0269 ntuple1->Draw("y:x>>histDoseRed","dose/1000","contz");
0270 //gPad->SetLogz();
0271 histDoseRed->Draw("contz");
0272 histDoseRed->GetXaxis()->SetTitle("X (microns)");
0273 histDoseRed->GetYaxis()->SetTitle("Y (mirons)");
0274 histDoseRed->GetZaxis()->SetTitle("Dose (mGy)");
0275 histDoseRed->GetXaxis()->SetLabelSize(0.025);
0276 histDoseRed->GetYaxis()->SetLabelSize(0.025);
0277 histDoseRed->GetZaxis()->SetLabelSize(0.025);
0278 histDoseRed->GetXaxis()->SetTitleSize(0.035);
0279 histDoseRed->GetYaxis()->SetTitleSize(0.035);
0280 histDoseRed->GetZaxis()->SetTitleSize(0.035);
0281 histDoseRed->GetXaxis()->SetTitleOffset(1.4);
0282 histDoseRed->GetYaxis()->SetTitleOffset(1.4);
0283 histDoseRed->GetZaxis()->SetTitleOffset(.6);
0284 histDoseRed->SetTitle("Dose map for RED voxels");
0285 
0286 c1->cd(8);
0287 TH2F *histDoseGreen = new TH2F("histDoseGreen","histDoseGreen",100,0,800,100,0,800);
0288 // WARNING : dose scaling to mGy
0289 ntuple2->Draw("y:x>>histDoseGreen","dose/1000","contz");
0290 //gPad->SetLogz();
0291 histDoseGreen->Draw("contz");
0292 histDoseGreen->GetXaxis()->SetTitle("X (microns)");
0293 histDoseGreen->GetYaxis()->SetTitle("Y (mirons)");
0294 histDoseGreen->GetZaxis()->SetTitle("Dose (mGy)");
0295 histDoseGreen->GetXaxis()->SetLabelSize(0.025);
0296 histDoseGreen->GetYaxis()->SetLabelSize(0.025);
0297 histDoseGreen->GetZaxis()->SetLabelSize(0.025);
0298 histDoseGreen->GetXaxis()->SetTitleSize(0.035);
0299 histDoseGreen->GetYaxis()->SetTitleSize(0.035);
0300 histDoseGreen->GetZaxis()->SetTitleSize(0.035);
0301 histDoseGreen->GetXaxis()->SetTitleOffset(1.4);
0302 histDoseGreen->GetYaxis()->SetTitleOffset(1.4);
0303 histDoseGreen->GetZaxis()->SetTitleOffset(.6);
0304 histDoseGreen->SetTitle("Dose map for GREEN voxels");
0305 
0306 c1->cd(12);
0307 TH2F *histDoseBlue = new TH2F("histDoseBlue","histDoseBlue",100,0,800,100,0,800);
0308 // WARNING : dose scaling to mGy
0309 ntuple3->Draw("y:x>>histDoseBlue","dose/1000","contz");
0310 //gPad->SetLogz();
0311 histDoseBlue->Draw("contz");
0312 histDoseBlue->GetXaxis()->SetTitle("X (microns)");
0313 histDoseBlue->GetYaxis()->SetTitle("Y (mirons)");
0314 histDoseBlue->GetZaxis()->SetTitle("Dose (mGy)");
0315 histDoseBlue->GetXaxis()->SetLabelSize(0.025);
0316 histDoseBlue->GetYaxis()->SetLabelSize(0.025);
0317 histDoseBlue->GetZaxis()->SetLabelSize(0.025);
0318 histDoseBlue->GetXaxis()->SetTitleSize(0.035);
0319 histDoseBlue->GetYaxis()->SetTitleSize(0.035);
0320 histDoseBlue->GetZaxis()->SetTitleSize(0.035);
0321 histDoseBlue->GetXaxis()->SetTitleOffset(1.4);
0322 histDoseBlue->GetYaxis()->SetTitleOffset(1.4);
0323 histDoseBlue->GetZaxis()->SetTitleOffset(.6);
0324 histDoseBlue->SetTitle("Dose map for BLUE voxels");
0325 
0326 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!