Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/microbeam/plot.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/microbeam/plot.C (Version 11.3.0) and /examples/advanced/microbeam/plot.C (Version 11.0)


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