Geant4 Cross Reference |
1 // ------------------------------------------- 1 // ------------------------------------------------------------------- >> 2 // $Id$ 2 // ------------------------------------------- 3 // ------------------------------------------------------------------- 3 // 4 // 4 // ******************************************* 5 // ********************************************************************* 5 // To execute this macro under ROOT, 6 // To execute this macro under ROOT, 6 // 1 - launch ROOT (usually type 'root' at y 7 // 1 - launch ROOT (usually type 'root' at your machine's prompt) 7 // 2 - type '.X plot.C' at the ROOT session 8 // 2 - type '.X plot.C' at the ROOT session prompt 8 // This macro needs five files : dose.txt, sto 9 // This macro needs five files : dose.txt, stoppingPower.txt, range.txt, 9 // 3DDose.txt and beamPosition.txt 10 // 3DDose.txt and beamPosition.txt 10 // written by S. Incerti and O. Boissonnade, 1 11 // written by S. Incerti and O. Boissonnade, 10/04/2006 11 // ******************************************* 12 // ********************************************************************* 12 { 13 { 13 14 14 gROOT->Reset(); 15 gROOT->Reset(); 15 16 16 //**************** 17 //**************** 17 // DOSE IN NUCLEUS 18 // DOSE IN NUCLEUS 18 //**************** 19 //**************** 19 20 20 gStyle->SetOptStat(0000); 21 gStyle->SetOptStat(0000); 21 gStyle->SetOptFit(); 22 gStyle->SetOptFit(); 22 gStyle->SetPalette(1); 23 gStyle->SetPalette(1); 23 gROOT->SetStyle("Plain"); 24 gROOT->SetStyle("Plain"); 24 Double_t scale; 25 Double_t scale; 25 26 26 27 27 c1 = new TCanvas ("c1","",20,20,1200,900); 28 c1 = new TCanvas ("c1","",20,20,1200,900); 28 c1->Divide(4,3); << 29 c1.Divide(4,3); 29 30 30 //********************* 31 //********************* 31 // INTENSITY HISTOGRAMS 32 // INTENSITY HISTOGRAMS 32 //********************* 33 //********************* >> 34 jump: 33 35 34 FILE * fp = fopen("phantom.dat","r"); 36 FILE * fp = fopen("phantom.dat","r"); 35 Float_t xVox, yVox, zVox, tmp, den, dose; 37 Float_t xVox, yVox, zVox, tmp, den, dose; 36 38 37 Float_t X,Y,Z; 39 Float_t X,Y,Z; 38 Float_t vox = 0, mat = 0; 40 Float_t vox = 0, mat = 0; 39 Float_t voxelSizeX, voxelSizeY, voxelSizeZ; 41 Float_t voxelSizeX, voxelSizeY, voxelSizeZ; 40 42 41 TH1F *h1 = new TH1F("h1","Nucleus marker inte 43 TH1F *h1 = new TH1F("h1","Nucleus marker intensity",100,1,300); 42 TH1F *h11 = new TH1F("h11 ","",100,1,300); 44 TH1F *h11 = new TH1F("h11 ","",100,1,300); 43 45 44 TH1F *h2 = new TH1F("h2","Cytoplasm marker in 46 TH1F *h2 = new TH1F("h2","Cytoplasm marker intensity",100,1,300); 45 TH1F *h20 = new TH1F("h20 ","",100,1,300); 47 TH1F *h20 = new TH1F("h20 ","",100,1,300); 46 48 47 TNtuple *ntupleYXN = new TNtuple("NUCLEUS","nt 49 TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox"); 48 TNtuple *ntupleZX = new TNtuple("CYTOPLASM","n 50 TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox"); 49 TNtuple *ntupleYX = new TNtuple("CYTOPLASM","n 51 TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox"); 50 52 51 Int_t nlines=0; << 53 nlines=0; 52 Int_t ncols=0; << 54 ncols=0; 53 55 54 while (1) 56 while (1) 55 { 57 { 56 if ( nlines == 0 ) ncols = fscanf(fp,"%f 58 if ( nlines == 0 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp); 57 if ( nlines == 1 ) ncols = fscanf(fp,"%f 59 if ( nlines == 1 ) ncols = fscanf(fp,"%f %f %f",&voxelSizeX,&voxelSizeY,&voxelSizeZ); 58 if ( nlines == 2 ) ncols = fscanf(fp,"%f 60 if ( nlines == 2 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp); 59 if ( nlines >= 3 ) ncols = fscanf(fp,"%f 61 if ( nlines >= 3 ) ncols = fscanf(fp,"%f %f %f %f %f %f", &X, &Y, &Z, &mat, &den, &vox); 60 if (ncols < 0) break; 62 if (ncols < 0) break; 61 63 62 X= X*voxelSizeX; 64 X= X*voxelSizeX; 63 Y= Y*voxelSizeY; 65 Y= Y*voxelSizeY; 64 Z= Z*voxelSizeZ; 66 Z= Z*voxelSizeZ; 65 67 66 if ( mat == 2 ) // noyau 68 if ( mat == 2 ) // noyau 67 { 69 { 68 if (den==1) h1->Fill( vox ); 70 if (den==1) h1->Fill( vox ); 69 if (den==2) h11->Fill( vox ); 71 if (den==2) h11->Fill( vox ); 70 ntupleYXN->Fill(Y,X,vox); 72 ntupleYXN->Fill(Y,X,vox); 71 } 73 } 72 if ( mat == 1 ) // cytoplasm 74 if ( mat == 1 ) // cytoplasm 73 { 75 { 74 if (den==1) h2->Fill( vox ); 76 if (den==1) h2->Fill( vox ); 75 if (den==2) h20->Fill( vox ); 77 if (den==2) h20->Fill( vox ); 76 ntupleZX->Fill(Z,X,vox); 78 ntupleZX->Fill(Z,X,vox); 77 ntupleYX->Fill(Y,X,vox); 79 ntupleYX->Fill(Y,X,vox); 78 } 80 } 79 nlines++; 81 nlines++; 80 82 81 } 83 } 82 fclose(fp); 84 fclose(fp); 83 85 84 // HISTO NUCLEUS 86 // HISTO NUCLEUS 85 87 86 c1->cd(1); << 88 c1.cd(1); 87 h1->Draw(); 89 h1->Draw(); 88 h1->GetXaxis()->SetLabelSize(0.025); 90 h1->GetXaxis()->SetLabelSize(0.025); 89 h1->GetYaxis()->SetLabelSize(0.025); 91 h1->GetYaxis()->SetLabelSize(0.025); 90 h1->GetXaxis()->SetTitleSize(0.035); 92 h1->GetXaxis()->SetTitleSize(0.035); 91 h1->GetYaxis()->SetTitleSize(0.035); 93 h1->GetYaxis()->SetTitleSize(0.035); 92 h1->GetXaxis()->SetTitleOffset(1.4); 94 h1->GetXaxis()->SetTitleOffset(1.4); 93 h1->GetYaxis()->SetTitleOffset(1.4); 95 h1->GetYaxis()->SetTitleOffset(1.4); 94 h1->GetXaxis()->SetTitle("Voxel intensity (0 96 h1->GetXaxis()->SetTitle("Voxel intensity (0-255)"); 95 h1->GetYaxis()->SetTitle("Number of events") 97 h1->GetYaxis()->SetTitle("Number of events"); 96 h1->SetLineColor(3); 98 h1->SetLineColor(3); 97 h1->SetFillColor(3); // green 99 h1->SetFillColor(3); // green 98 100 99 h11->SetLineColor(8); 101 h11->SetLineColor(8); 100 h11->SetFillColor(8); // dark green 102 h11->SetFillColor(8); // dark green 101 h11->Draw("same"); 103 h11->Draw("same"); 102 104 103 // HISTO CYTOPLASM 105 // HISTO CYTOPLASM 104 106 105 c1->cd(5); << 107 c1.cd(5); 106 h2->Draw(); 108 h2->Draw(); 107 h2->GetXaxis()->SetLabelSize(0.025); 109 h2->GetXaxis()->SetLabelSize(0.025); 108 h2->GetYaxis()->SetLabelSize(0.025); 110 h2->GetYaxis()->SetLabelSize(0.025); 109 h2->GetXaxis()->SetTitleSize(0.035); 111 h2->GetXaxis()->SetTitleSize(0.035); 110 h2->GetYaxis()->SetTitleSize(0.035); 112 h2->GetYaxis()->SetTitleSize(0.035); 111 h2->GetXaxis()->SetTitleOffset(1.4); 113 h2->GetXaxis()->SetTitleOffset(1.4); 112 h2->GetYaxis()->SetTitleOffset(1.4); 114 h2->GetYaxis()->SetTitleOffset(1.4); 113 h2->GetXaxis()->SetTitle("Voxel intensity (0 115 h2->GetXaxis()->SetTitle("Voxel intensity (0-255)"); 114 h2->GetYaxis()->SetTitle("Number of events") 116 h2->GetYaxis()->SetTitle("Number of events"); 115 h2->SetLineColor(2); 117 h2->SetLineColor(2); 116 h2->SetFillColor(2); // red 118 h2->SetFillColor(2); // red 117 119 118 h20->SetLineColor(5); 120 h20->SetLineColor(5); 119 h20->SetFillColor(5); // yellow (nucleoli) 121 h20->SetFillColor(5); // yellow (nucleoli) 120 h20->Draw("same"); 122 h20->Draw("same"); 121 123 122 //************************* 124 //************************* 123 // CUMULATED CELL INTENSITY 125 // CUMULATED CELL INTENSITY 124 //************************* 126 //************************* 125 127 126 gStyle->SetOptStat(0000); 128 gStyle->SetOptStat(0000); 127 gStyle->SetOptFit(); 129 gStyle->SetOptFit(); 128 gStyle->SetPalette(1); 130 gStyle->SetPalette(1); 129 gROOT->SetStyle("Plain"); 131 gROOT->SetStyle("Plain"); 130 132 131 //CYTOPLASM 133 //CYTOPLASM 132 134 133 c1->cd(7); // axe YX << 135 c1.cd(7); // axe YX 134 TH2F *hist = new TH2F("hist","hist",50,-20,2 136 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20); 135 ntupleYX->Draw("Y:X>>hist","vox","colz"); 137 ntupleYX->Draw("Y:X>>hist","vox","colz"); 136 gPad->SetLogz(); 138 gPad->SetLogz(); 137 hist->Draw("colz"); 139 hist->Draw("colz"); 138 hist->GetXaxis()->SetLabelSize(0.025); 140 hist->GetXaxis()->SetLabelSize(0.025); 139 hist->GetYaxis()->SetLabelSize(0.025); 141 hist->GetYaxis()->SetLabelSize(0.025); 140 hist->GetZaxis()->SetLabelSize(0.025); 142 hist->GetZaxis()->SetLabelSize(0.025); 141 hist->GetXaxis()->SetTitleSize(0.035); 143 hist->GetXaxis()->SetTitleSize(0.035); 142 hist->GetYaxis()->SetTitleSize(0.035); 144 hist->GetYaxis()->SetTitleSize(0.035); 143 hist->GetXaxis()->SetTitleOffset(1.4); 145 hist->GetXaxis()->SetTitleOffset(1.4); 144 hist->GetYaxis()->SetTitleOffset(1.4); 146 hist->GetYaxis()->SetTitleOffset(1.4); 145 hist->GetXaxis()->SetTitle("Y (um)"); << 147 hist->GetXaxis()->SetTitle("Y (µm)"); 146 hist->GetYaxis()->SetTitle("X (um)"); << 148 hist->GetYaxis()->SetTitle("X (µm)"); 147 hist->SetTitle("Cytoplasm intensity on trans 149 hist->SetTitle("Cytoplasm intensity on transverse section"); 148 150 149 //NUCLEUS 151 //NUCLEUS 150 152 151 c1->cd(3); // axe YX << 153 c1.cd(3); // axe YX 152 TH2F *hist2 = new TH2F("hist2","hist2",50,-2 << 154 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20); 153 ntupleYXN->Draw("Y:X>>hist2","vox","colz"); << 155 ntupleYXN->Draw("Y:X>>hist","vox","colz"); 154 gPad->SetLogz(); 156 gPad->SetLogz(); 155 hist2->Draw("colz"); << 157 hist->Draw("colz"); 156 hist2->GetXaxis()->SetLabelSize(0.025); << 158 hist->GetXaxis()->SetLabelSize(0.025); 157 hist2->GetYaxis()->SetLabelSize(0.025); << 159 hist->GetYaxis()->SetLabelSize(0.025); 158 hist2->GetZaxis()->SetLabelSize(0.025); << 160 hist->GetZaxis()->SetLabelSize(0.025); 159 hist2->GetXaxis()->SetTitleSize(0.035); << 161 hist->GetXaxis()->SetTitleSize(0.035); 160 hist2->GetYaxis()->SetTitleSize(0.035); << 162 hist->GetYaxis()->SetTitleSize(0.035); 161 hist2->GetXaxis()->SetTitleOffset(1.4); << 163 hist->GetXaxis()->SetTitleOffset(1.4); 162 hist2->GetYaxis()->SetTitleOffset(1.4); << 164 hist->GetYaxis()->SetTitleOffset(1.4); 163 hist2->GetXaxis()->SetTitle("Y (um)"); << 165 hist->GetXaxis()->SetTitle("Y (µm)"); 164 hist2->GetYaxis()->SetTitle("X (um)"); << 166 hist->GetYaxis()->SetTitle("X (µm)"); 165 hist2->SetTitle("Nucleus intensity on transv << 167 hist->SetTitle("Nucleus intensity on transverse section"); 166 168 167 // 169 // 168 170 169 system ("rm -rf microbeam.root"); << 170 system ("hadd -O microbeam.root microbeam_*.ro << 171 << 172 TFile f("microbeam.root"); 171 TFile f("microbeam.root"); 173 172 174 TNtuple* ntuple0; 173 TNtuple* ntuple0; 175 TNtuple* ntuple1; 174 TNtuple* ntuple1; 176 TNtuple* ntuple2; 175 TNtuple* ntuple2; 177 TNtuple* ntuple3; 176 TNtuple* ntuple3; 178 TNtuple* ntuple4; 177 TNtuple* ntuple4; 179 178 180 ntuple0 = (TNtuple*)f.Get("ntuple0"); << 179 ntuple0 = (TNtuple*)f->Get("ntuple0"); 181 ntuple1 = (TNtuple*)f.Get("ntuple1"); << 180 ntuple1 = (TNtuple*)f->Get("ntuple1"); 182 ntuple2 = (TNtuple*)f.Get("ntuple2"); << 181 ntuple2 = (TNtuple*)f->Get("ntuple2"); 183 ntuple3 = (TNtuple*)f.Get("ntuple3"); << 182 ntuple3 = (TNtuple*)f->Get("ntuple3"); 184 ntuple4 = (TNtuple*)f.Get("ntuple4"); << 183 ntuple4 = (TNtuple*)f->Get("ntuple4"); 185 << 184 186 TH1F *h1bis = new TH1F("h1bis","Dose distribu << 185 TH1F *h1 = new TH1F("h1","Dose distribution in Nucleus",100,0.001,1.); 187 TH1F *h10 = new TH1F("h10bis","Dose distributi << 186 TH1F *h10 = new TH1F("h10","Dose distribution in Cytoplasm",100,0.001,.2); 188 << 187 189 c1->cd(2); << 188 c1.cd(2); 190 << 189 191 ntuple3->Project("h1bis","doseN"); << 190 ntuple3->Project("h1","doseN"); 192 scale = 1/h1bis->Integral(); << 191 scale = 1/h1->Integral(); 193 h1bis->Scale(scale); << 192 h1->Scale(scale); 194 h1bis->Draw(); << 193 h1->Draw(); 195 h1bis->GetXaxis()->SetLabelSize(0.025); << 194 h1->GetXaxis()->SetLabelSize(0.025); 196 h1bis->GetYaxis()->SetLabelSize(0.025); << 195 h1->GetYaxis()->SetLabelSize(0.025); 197 h1bis->GetXaxis()->SetTitleSize(0.035); << 196 h1->GetXaxis()->SetTitleSize(0.035); 198 h1bis->GetYaxis()->SetTitleSize(0.035); << 197 h1->GetYaxis()->SetTitleSize(0.035); 199 h1bis->GetXaxis()->SetTitleOffset(1.4); << 198 h1->GetXaxis()->SetTitleOffset(1.4); 200 h1bis->GetYaxis()->SetTitleOffset(1.4); << 199 h1->GetYaxis()->SetTitleOffset(1.4); 201 h1bis->GetXaxis()->SetTitle("Absorbed dose ( << 200 h1->GetXaxis()->SetTitle("Absorbed dose (Gy)"); 202 h1bis->GetYaxis()->SetTitle("Fraction of eve << 201 h1->GetYaxis()->SetTitle("Fraction of events"); 203 h1bis->SetLineColor(3); << 202 h1->SetLineColor(3); 204 h1bis->SetFillColor(3); << 203 h1->SetFillColor(3); 205 204 206 //***************** 205 //***************** 207 // DOSE IN CYTOPLASM 206 // DOSE IN CYTOPLASM 208 //***************** 207 //***************** 209 208 210 c1->cd(6); << 209 c1.cd(6); 211 ntuple3->Project("h10bis","doseC"); << 210 ntuple3->Project("h10","doseC"); 212 scale = 1/h10->Integral(); 211 scale = 1/h10->Integral(); 213 h10->Scale(scale); 212 h10->Scale(scale); 214 h10->Draw(); 213 h10->Draw(); 215 h10->GetXaxis()->SetLabelSize(0.025); 214 h10->GetXaxis()->SetLabelSize(0.025); 216 h10->GetYaxis()->SetLabelSize(0.025); 215 h10->GetYaxis()->SetLabelSize(0.025); 217 h10->GetXaxis()->SetTitleSize(0.035); 216 h10->GetXaxis()->SetTitleSize(0.035); 218 h10->GetYaxis()->SetTitleSize(0.035); 217 h10->GetYaxis()->SetTitleSize(0.035); 219 h10->GetXaxis()->SetTitleOffset(1.4); 218 h10->GetXaxis()->SetTitleOffset(1.4); 220 h10->GetYaxis()->SetTitleOffset(1.4); 219 h10->GetYaxis()->SetTitleOffset(1.4); 221 h10->GetXaxis()->SetTitle("Absorbed dose (Gy 220 h10->GetXaxis()->SetTitle("Absorbed dose (Gy)"); 222 h10->GetYaxis()->SetTitle("Fraction of event 221 h10->GetYaxis()->SetTitle("Fraction of events"); 223 h10->SetLineColor(2); 222 h10->SetLineColor(2); 224 h10->SetFillColor(2); 223 h10->SetFillColor(2); 225 224 226 //******************************** 225 //******************************** 227 // STOPPING POWER AT CELL ENTRANCE 226 // STOPPING POWER AT CELL ENTRANCE 228 //******************************** 227 //******************************** 229 228 230 gStyle->SetOptStat(0000); 229 gStyle->SetOptStat(0000); 231 gStyle->SetOptFit(); 230 gStyle->SetOptFit(); 232 gStyle->SetPalette(1); 231 gStyle->SetPalette(1); 233 gROOT->SetStyle("Plain"); 232 gROOT->SetStyle("Plain"); 234 233 235 Float_t d; 234 Float_t d; 236 235 237 TH1F *h2bis = new TH1F("h2bis","Beam stopping << 236 TH1F *h2 = new TH1F("h2","Beam stopping power at cell entrance",200,0,300); 238 237 239 c1->cd(9); << 238 c1.cd(9); 240 ntuple0->Project("h2bis","sp"); << 239 ntuple0->Project("h2","sp"); 241 scale = 1/h2bis->Integral(); << 240 scale = 1/h2->Integral(); 242 h2bis->Scale(scale); << 241 h2->Scale(scale); 243 h2bis->Draw(); << 242 h2->Draw(); 244 h2bis->GetXaxis()->SetLabelSize(0.025); << 243 h2->GetXaxis()->SetLabelSize(0.025); 245 h2bis->GetYaxis()->SetLabelSize(0.025); << 244 h2->GetYaxis()->SetLabelSize(0.025); 246 h2bis->GetXaxis()->SetTitleSize(0.035); << 245 h2->GetXaxis()->SetTitleSize(0.035); 247 h2bis->GetYaxis()->SetTitleSize(0.035); << 246 h2->GetYaxis()->SetTitleSize(0.035); 248 h2bis->GetXaxis()->SetTitleOffset(1.4); << 247 h2->GetXaxis()->SetTitleOffset(1.4); 249 h2bis->GetYaxis()->SetTitleOffset(1.4); << 248 h2->GetYaxis()->SetTitleOffset(1.4); 250 h2bis->GetXaxis()->SetTitle("dE/dx (keV/um << 249 h2->GetXaxis()->SetTitle("dE/dx (keV/µm)"); 251 h2bis->GetYaxis()->SetTitle("Fraction of eve << 250 h2->GetYaxis()->SetTitle("Fraction of events"); 252 h2bis->SetTitle("dE/dx at cell entrance"); << 251 h2->SetTitle("dE/dx at cell entrance"); 253 h2bis->SetFillColor(4); << 252 h2->SetFillColor(4); 254 h2bis->SetLineColor(4); << 253 h2->SetLineColor(4); 255 h2bis->Fit("gaus"); << 254 h2->Fit("gaus"); 256 gaus->SetLineColor(6); 255 gaus->SetLineColor(6); 257 h2bis->Fit("gaus"); << 256 h2->Fit("gaus"); >> 257 258 258 259 //************** 259 //************** 260 // RANGE IN CELL 260 // RANGE IN CELL 261 //************** 261 //************** 262 262 263 Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2; 263 Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2; 264 264 265 // X position of target in World 265 // X position of target in World 266 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180 266 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180); 267 267 268 // Z position of target in World 268 // Z position of target in World 269 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180); 269 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180); 270 270 271 // Line alignment (cf MicrobeamEMField.cc) 271 // Line alignment (cf MicrobeamEMField.cc) 272 Xc = Xc + 5.24*cos(10*TMath::Pi()/180); 272 Xc = Xc + 5.24*cos(10*TMath::Pi()/180); 273 Zc = Zc + 5.24*sin(10*TMath::Pi()/180); 273 Zc = Zc + 5.24*sin(10*TMath::Pi()/180); 274 274 275 TNtuple *ntupleR = new TNtuple("Rmax","ntuple" 275 TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2"); 276 Double_t x,y,z,xx,zz; 276 Double_t x,y,z,xx,zz; 277 ntuple2->SetBranchAddress("x",&x); 277 ntuple2->SetBranchAddress("x",&x); 278 ntuple2->SetBranchAddress("y",&y); 278 ntuple2->SetBranchAddress("y",&y); 279 ntuple2->SetBranchAddress("z",&z); 279 ntuple2->SetBranchAddress("z",&z); 280 Int_t nentries = (Int_t)ntuple2->GetEntries(); 280 Int_t nentries = (Int_t)ntuple2->GetEntries(); 281 for (Int_t i=0;i<nentries;i++) 281 for (Int_t i=0;i<nentries;i++) 282 { 282 { 283 ntuple2->GetEntry(i); 283 ntuple2->GetEntry(i); 284 X1=x; 284 X1=x; 285 Y1=y; 285 Y1=y; 286 Z1=z; 286 Z1=z; 287 xx = X1-Xc; 287 xx = X1-Xc; 288 zz = Z1-Zc; 288 zz = Z1-Zc; 289 Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(1 289 Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(10*TMath::Pi()/180); 290 X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(1 290 X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(10*TMath::Pi()/180); 291 Y2 = Y1; 291 Y2 = Y1; 292 ntupleR->Fill(Z2,Y2,X2); 292 ntupleR->Fill(Z2,Y2,X2); 293 } 293 } 294 294 295 c1->cd(10); << 295 c1.cd(10); 296 ntupleR->Draw("X2:Z2","abs(X2)<50","surf3"); 296 ntupleR->Draw("X2:Z2","abs(X2)<50","surf3"); 297 gPad->SetLogz(); 297 gPad->SetLogz(); >> 298 /* >> 299 htemp->GetXaxis()->SetLabelSize(0.025); >> 300 htemp->GetYaxis()->SetLabelSize(0.025); >> 301 htemp->GetZaxis()->SetLabelSize(0.025); >> 302 htemp->GetXaxis()->SetTitleSize(0.035); >> 303 htemp->GetYaxis()->SetTitleSize(0.035); >> 304 htemp->GetXaxis()->SetTitleOffset(1.4); >> 305 htemp->GetYaxis()->SetTitleOffset(1.4); >> 306 htemp->GetXaxis()->SetTitle("Z (µm)"); >> 307 htemp->GetYaxis()->SetTitle("X (µm)"); >> 308 htemp->SetTitle("Range in cell"); >> 309 */ 298 310 299 //**************** 311 //**************** 300 // ENERGY DEPOSITS 312 // ENERGY DEPOSITS 301 //**************** 313 //**************** 302 314 >> 315 303 gStyle->SetOptStat(0000); 316 gStyle->SetOptStat(0000); 304 gStyle->SetOptFit(); 317 gStyle->SetOptFit(); 305 gStyle->SetPalette(1); 318 gStyle->SetPalette(1); 306 gROOT->SetStyle("Plain"); 319 gROOT->SetStyle("Plain"); 307 320 308 c1->cd(11); << 321 c1.cd(11); 309 TH2F *histbis = new TH2F("histbis","histbis" << 322 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20); 310 ntuple4->Draw("y*0.359060:x*0.359060>>histbi << 323 ntuple4->Draw("y*0.359060:x*0.359060>>hist","doseV","contz"); 311 gPad->SetLogz(); 324 gPad->SetLogz(); 312 histbis->Draw("contz"); << 325 hist->Draw("contz"); 313 histbis->GetXaxis()->SetLabelSize(0.025); << 326 hist->GetXaxis()->SetLabelSize(0.025); 314 histbis->GetYaxis()->SetLabelSize(0.025); << 327 hist->GetYaxis()->SetLabelSize(0.025); 315 histbis->GetZaxis()->SetLabelSize(0.025); << 328 hist->GetZaxis()->SetLabelSize(0.025); 316 histbis->GetXaxis()->SetTitleSize(0.035); << 329 hist->GetXaxis()->SetTitleSize(0.035); 317 histbis->GetYaxis()->SetTitleSize(0.035); << 330 hist->GetYaxis()->SetTitleSize(0.035); 318 histbis->GetXaxis()->SetTitleOffset(1.4); << 331 hist->GetXaxis()->SetTitleOffset(1.4); 319 histbis->GetYaxis()->SetTitleOffset(1.4); << 332 hist->GetYaxis()->SetTitleOffset(1.4); 320 histbis->GetXaxis()->SetTitle("Y (um)"); << 333 hist->GetXaxis()->SetTitle("Y (µm)"); 321 histbis->GetYaxis()->SetTitle("X (um)"); << 334 hist->GetYaxis()->SetTitle("X (µm)"); 322 histbis->SetTitle("Energy deposit -transvers << 335 hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)"); 323 << 336 324 c1->cd(12); << 337 325 TH2F *histter = new TH2F("histter","histter" << 338 c1.cd(12); 326 ntuple4->Draw("x*0.359060:(z+1500/0.162810+2 << 339 TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20); >> 340 ntuple4->Draw("x*0.359060:z*0.162810>>hist","doseV","contz"); 327 gPad->SetLogz(); 341 gPad->SetLogz(); 328 histter->Draw("contz"); << 342 hist->Draw("contz"); 329 histter->GetXaxis()->SetLabelSize(0.025); << 343 hist->GetXaxis()->SetLabelSize(0.025); 330 histter->GetYaxis()->SetLabelSize(0.025); << 344 hist->GetYaxis()->SetLabelSize(0.025); 331 histter->GetZaxis()->SetLabelSize(0.025); << 345 hist->GetZaxis()->SetLabelSize(0.025); 332 histter->GetXaxis()->SetTitleSize(0.035); << 346 hist->GetXaxis()->SetTitleSize(0.035); 333 histter->GetYaxis()->SetTitleSize(0.035); << 347 hist->GetYaxis()->SetTitleSize(0.035); 334 histter->GetXaxis()->SetTitleOffset(1.4); << 348 hist->GetXaxis()->SetTitleOffset(1.4); 335 histter->GetYaxis()->SetTitleOffset(1.4); << 349 hist->GetYaxis()->SetTitleOffset(1.4); 336 histter->GetXaxis()->SetTitle("Z (um)"); << 350 hist->GetXaxis()->SetTitle("Z (µm)"); 337 histter->GetYaxis()->SetTitle("X (um)"); << 351 hist->GetYaxis()->SetTitle("X (µm)"); 338 histter->SetTitle("Energy deposit -longitudi << 352 hist->SetTitle("Mean energy deposit -longitudinal- (z axis in eV)"); 339 353 340 //******************************* 354 //******************************* 341 // BEAM POSITION AT CELL ENTRANCE 355 // BEAM POSITION AT CELL ENTRANCE 342 //******************************* 356 //******************************* 343 357 344 gStyle->SetOptStat(0000); 358 gStyle->SetOptStat(0000); 345 gStyle->SetOptFit(); 359 gStyle->SetOptFit(); 346 gStyle->SetPalette(1); 360 gStyle->SetPalette(1); 347 gROOT->SetStyle("Plain"); 361 gROOT->SetStyle("Plain"); 348 362 349 TH1F *h77 = new TH1F("hx","h1",200,-10,10); 363 TH1F *h77 = new TH1F("hx","h1",200,-10,10); 350 TH1F *h88 = new TH1F("hy","h1",200,-10,10); 364 TH1F *h88 = new TH1F("hy","h1",200,-10,10); 351 365 352 c1->cd(4); << 366 c1.cd(4); 353 ntuple1->Project("hx","x"); 367 ntuple1->Project("hx","x"); 354 scale = 1/h77->Integral(); 368 scale = 1/h77->Integral(); 355 h77->Scale(scale); 369 h77->Scale(scale); 356 h77->Draw(); 370 h77->Draw(); 357 h77->GetXaxis()->SetLabelSize(0.025); 371 h77->GetXaxis()->SetLabelSize(0.025); 358 h77->GetYaxis()->SetLabelSize(0.025); 372 h77->GetYaxis()->SetLabelSize(0.025); 359 h77->GetXaxis()->SetTitleSize(0.035); 373 h77->GetXaxis()->SetTitleSize(0.035); 360 h77->GetYaxis()->SetTitleSize(0.035); 374 h77->GetYaxis()->SetTitleSize(0.035); 361 h77->GetXaxis()->SetTitleOffset(1.4); 375 h77->GetXaxis()->SetTitleOffset(1.4); 362 h77->GetYaxis()->SetTitleOffset(1.4); 376 h77->GetYaxis()->SetTitleOffset(1.4); 363 h77->GetXaxis()->SetTitle("Position (um)") << 377 h77->GetXaxis()->SetTitle("Position (µm)"); 364 h77->GetYaxis()->SetTitle("Fraction of event 378 h77->GetYaxis()->SetTitle("Fraction of events"); 365 h77->SetTitle("Beam X position on cell"); 379 h77->SetTitle("Beam X position on cell"); 366 h77->SetFillColor(4); 380 h77->SetFillColor(4); 367 h77->SetLineColor(4); 381 h77->SetLineColor(4); >> 382 gaus->SetLineColor(6); 368 h77->Fit("gaus"); 383 h77->Fit("gaus"); 369 384 370 c1->cd(8); << 385 c1.cd(8); 371 ntuple1->Project("hy","y"); 386 ntuple1->Project("hy","y"); 372 scale = 1/h88->Integral(); 387 scale = 1/h88->Integral(); 373 h88->Scale(scale); 388 h88->Scale(scale); 374 h88->Draw(); 389 h88->Draw(); 375 h88->GetXaxis()->SetLabelSize(0.025); 390 h88->GetXaxis()->SetLabelSize(0.025); 376 h88->GetYaxis()->SetLabelSize(0.025); 391 h88->GetYaxis()->SetLabelSize(0.025); 377 h88->GetXaxis()->SetTitleSize(0.035); 392 h88->GetXaxis()->SetTitleSize(0.035); 378 h88->GetYaxis()->SetTitleSize(0.035); 393 h88->GetYaxis()->SetTitleSize(0.035); 379 h88->GetXaxis()->SetTitleOffset(1.4); 394 h88->GetXaxis()->SetTitleOffset(1.4); 380 h88->GetYaxis()->SetTitleOffset(1.4); 395 h88->GetYaxis()->SetTitleOffset(1.4); 381 h88->GetXaxis()->SetTitle("Position (um)") << 396 h88->GetXaxis()->SetTitle("Position (µm)"); 382 h88->GetYaxis()->SetTitle("Fraction of event 397 h88->GetYaxis()->SetTitle("Fraction of events"); 383 h88->SetTitle("Beam Y position on cell"); 398 h88->SetTitle("Beam Y position on cell"); 384 h88->SetFillColor(4); 399 h88->SetFillColor(4); 385 h88->SetLineColor(4); 400 h88->SetLineColor(4); >> 401 gaus->SetLineColor(6); 386 h88->Fit("gaus"); 402 h88->Fit("gaus"); >> 403 387 } 404 } 388 405