Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 09:20:24 

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 // This macro needs five files : dose.txt, stoppingPower.txt, range.txt,
0009 // 3DDose.txt and beamPosition.txt
0010 // written by S. Incerti and O. Boissonnade, 10/04/2006
0011 // *********************************************************************
0012 {
0013 
0014 gROOT->Reset();
0015 
0016 //****************
0017 // DOSE IN NUCLEUS
0018 //****************
0019 
0020 gStyle->SetOptStat(0000);
0021 gStyle->SetOptFit();
0022 gStyle->SetPalette(1);
0023 gROOT->SetStyle("Plain");
0024 Double_t scale;
0025 
0026 
0027 c1 = new TCanvas ("c1","",20,20,1200,900);
0028 c1->Divide(4,3);
0029 
0030 //*********************
0031 // INTENSITY HISTOGRAMS 
0032 //*********************
0033 
0034 FILE * fp = fopen("phantom.dat","r");
0035 Float_t xVox, yVox, zVox, tmp, den, dose;
0036 
0037 Float_t X,Y,Z;
0038 Float_t vox = 0,  mat = 0;
0039 Float_t voxelSizeX, voxelSizeY, voxelSizeZ;
0040 
0041 TH1F *h1  = new TH1F("h1","Nucleus marker intensity",100,1,300);
0042 TH1F *h11 = new TH1F("h11 ","",100,1,300);
0043 
0044 TH1F *h2  = new TH1F("h2","Cytoplasm marker intensity",100,1,300);
0045 TH1F *h20 = new TH1F("h20 ","",100,1,300);
0046 
0047 TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
0048 TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox");
0049 TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox");
0050 
0051 Int_t nlines=0;
0052 Int_t ncols=0;
0053 
0054 while (1) 
0055    {
0056       if ( nlines == 0 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
0057       if ( nlines == 1 ) ncols = fscanf(fp,"%f %f %f",&voxelSizeX,&voxelSizeY,&voxelSizeZ);
0058       if ( nlines == 2 ) ncols = fscanf(fp,"%f %f %f",&tmp,&tmp,&tmp);
0059       if ( nlines >= 3 ) ncols = fscanf(fp,"%f %f %f %f %f %f", &X, &Y, &Z, &mat, &den, &vox);
0060       if (ncols < 0) break;
0061 
0062       X= X*voxelSizeX;
0063       Y= Y*voxelSizeY;
0064       Z= Z*voxelSizeZ;
0065   
0066       if ( mat == 2 )  // noyau
0067          {
0068       if (den==1) h1->Fill( vox );
0069       if (den==2) h11->Fill( vox );
0070       ntupleYXN->Fill(Y,X,vox);
0071      }
0072       if ( mat == 1 ) // cytoplasm
0073          {
0074       if (den==1) h2->Fill( vox );
0075       if (den==2) h20->Fill( vox );
0076       ntupleZX->Fill(Z,X,vox);
0077       ntupleYX->Fill(Y,X,vox);
0078      }
0079       nlines++;
0080       
0081    }
0082 fclose(fp);
0083 
0084 // HISTO NUCLEUS
0085 
0086 c1->cd(1);
0087     h1->Draw();
0088     h1->GetXaxis()->SetLabelSize(0.025);
0089     h1->GetYaxis()->SetLabelSize(0.025);
0090     h1->GetXaxis()->SetTitleSize(0.035);
0091     h1->GetYaxis()->SetTitleSize(0.035);
0092     h1->GetXaxis()->SetTitleOffset(1.4);
0093     h1->GetYaxis()->SetTitleOffset(1.4);
0094     h1->GetXaxis()->SetTitle("Voxel intensity (0-255)");
0095     h1->GetYaxis()->SetTitle("Number of events");  
0096     h1->SetLineColor(3);
0097     h1->SetFillColor(3);   // green
0098 
0099     h11->SetLineColor(8);
0100     h11->SetFillColor(8);  // dark green
0101     h11->Draw("same");
0102 
0103 // HISTO CYTOPLASM
0104 
0105 c1->cd(5);
0106     h2->Draw();
0107     h2->GetXaxis()->SetLabelSize(0.025);
0108     h2->GetYaxis()->SetLabelSize(0.025);
0109     h2->GetXaxis()->SetTitleSize(0.035);
0110     h2->GetYaxis()->SetTitleSize(0.035);
0111     h2->GetXaxis()->SetTitleOffset(1.4);
0112     h2->GetYaxis()->SetTitleOffset(1.4);
0113     h2->GetXaxis()->SetTitle("Voxel intensity (0-255)");
0114     h2->GetYaxis()->SetTitle("Number of events");  
0115     h2->SetLineColor(2);
0116     h2->SetFillColor(2);   // red
0117     
0118     h20->SetLineColor(5);
0119     h20->SetFillColor(5);  // yellow (nucleoli)
0120     h20->Draw("same");
0121 
0122 //*************************
0123 // CUMULATED CELL INTENSITY 
0124 //*************************
0125 
0126 gStyle->SetOptStat(0000);
0127 gStyle->SetOptFit();
0128 gStyle->SetPalette(1);
0129 gROOT->SetStyle("Plain");
0130 
0131 //CYTOPLASM
0132 
0133 c1->cd(7);  // axe YX
0134   TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
0135   ntupleYX->Draw("Y:X>>hist","vox","colz");
0136   gPad->SetLogz();
0137   hist->Draw("colz");
0138   hist->GetXaxis()->SetLabelSize(0.025);
0139   hist->GetYaxis()->SetLabelSize(0.025);
0140   hist->GetZaxis()->SetLabelSize(0.025);
0141   hist->GetXaxis()->SetTitleSize(0.035);
0142   hist->GetYaxis()->SetTitleSize(0.035);
0143   hist->GetXaxis()->SetTitleOffset(1.4);
0144   hist->GetYaxis()->SetTitleOffset(1.4);
0145   hist->GetXaxis()->SetTitle("Y (um)");
0146   hist->GetYaxis()->SetTitle("X (um)");
0147   hist->SetTitle("Cytoplasm intensity on transverse section");
0148 
0149 //NUCLEUS
0150 
0151 c1->cd(3);  // axe YX
0152   TH2F *hist2 = new TH2F("hist2","hist2",50,-20,20,50,-20,20);
0153   ntupleYXN->Draw("Y:X>>hist2","vox","colz");
0154   gPad->SetLogz();
0155   hist2->Draw("colz");
0156   hist2->GetXaxis()->SetLabelSize(0.025);
0157   hist2->GetYaxis()->SetLabelSize(0.025);
0158   hist2->GetZaxis()->SetLabelSize(0.025);
0159   hist2->GetXaxis()->SetTitleSize(0.035);
0160   hist2->GetYaxis()->SetTitleSize(0.035);
0161   hist2->GetXaxis()->SetTitleOffset(1.4);
0162   hist2->GetYaxis()->SetTitleOffset(1.4);
0163   hist2->GetXaxis()->SetTitle("Y (um)");
0164   hist2->GetYaxis()->SetTitle("X (um)");
0165   hist2->SetTitle("Nucleus intensity on transverse section");
0166 
0167 //
0168 
0169 system ("rm -rf microbeam.root");
0170 system ("hadd -O microbeam.root microbeam_*.root");
0171 
0172 TFile f("microbeam.root"); 
0173 
0174 TNtuple* ntuple0;
0175 TNtuple* ntuple1;
0176 TNtuple* ntuple2;
0177 TNtuple* ntuple3;
0178 TNtuple* ntuple4;
0179 
0180 ntuple0 = (TNtuple*)f.Get("ntuple0"); 
0181 ntuple1 = (TNtuple*)f.Get("ntuple1"); 
0182 ntuple2 = (TNtuple*)f.Get("ntuple2"); 
0183 ntuple3 = (TNtuple*)f.Get("ntuple3"); 
0184 ntuple4 = (TNtuple*)f.Get("ntuple4"); 
0185 
0186 TH1F *h1bis  = new TH1F("h1bis","Dose distribution in Nucleus",100,0.001,1.);
0187 TH1F *h10 = new TH1F("h10bis","Dose distribution in Cytoplasm",100,0.001,.2);
0188 
0189 c1->cd(2);
0190 
0191         ntuple3->Project("h1bis","doseN");
0192     scale = 1/h1bis->Integral();
0193     h1bis->Scale(scale);
0194     h1bis->Draw();
0195     h1bis->GetXaxis()->SetLabelSize(0.025);
0196     h1bis->GetYaxis()->SetLabelSize(0.025);
0197     h1bis->GetXaxis()->SetTitleSize(0.035);
0198     h1bis->GetYaxis()->SetTitleSize(0.035);
0199     h1bis->GetXaxis()->SetTitleOffset(1.4);
0200     h1bis->GetYaxis()->SetTitleOffset(1.4);
0201     h1bis->GetXaxis()->SetTitle("Absorbed dose (Gy)");
0202     h1bis->GetYaxis()->SetTitle("Fraction of events");
0203     h1bis->SetLineColor(3);
0204     h1bis->SetFillColor(3);
0205 
0206 //*****************
0207 // DOSE IN CYTOPLASM
0208 //*****************
0209 
0210 c1->cd(6);
0211         ntuple3->Project("h10bis","doseC");
0212         scale = 1/h10->Integral();
0213     h10->Scale(scale);
0214     h10->Draw();
0215     h10->GetXaxis()->SetLabelSize(0.025);
0216     h10->GetYaxis()->SetLabelSize(0.025);
0217     h10->GetXaxis()->SetTitleSize(0.035);
0218     h10->GetYaxis()->SetTitleSize(0.035);
0219     h10->GetXaxis()->SetTitleOffset(1.4);
0220     h10->GetYaxis()->SetTitleOffset(1.4);
0221     h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
0222     h10->GetYaxis()->SetTitle("Fraction of events");
0223     h10->SetLineColor(2);
0224     h10->SetFillColor(2);
0225 
0226 //********************************
0227 // STOPPING POWER AT CELL ENTRANCE
0228 //********************************
0229  
0230 gStyle->SetOptStat(0000);
0231 gStyle->SetOptFit();
0232 gStyle->SetPalette(1);
0233 gROOT->SetStyle("Plain");
0234 
0235 Float_t d;
0236 
0237 TH1F *h2bis = new TH1F("h2bis","Beam stopping power at cell entrance",200,0,300); 
0238     
0239 c1->cd(9);
0240         ntuple0->Project("h2bis","sp");
0241     scale = 1/h2bis->Integral();
0242     h2bis->Scale(scale);
0243     h2bis->Draw();
0244     h2bis->GetXaxis()->SetLabelSize(0.025);
0245     h2bis->GetYaxis()->SetLabelSize(0.025);
0246     h2bis->GetXaxis()->SetTitleSize(0.035);
0247     h2bis->GetYaxis()->SetTitleSize(0.035);
0248     h2bis->GetXaxis()->SetTitleOffset(1.4);
0249     h2bis->GetYaxis()->SetTitleOffset(1.4);
0250     h2bis->GetXaxis()->SetTitle("dE/dx (keV/um)");
0251     h2bis->GetYaxis()->SetTitle("Fraction of events");
0252     h2bis->SetTitle("dE/dx at cell entrance");
0253     h2bis->SetFillColor(4);
0254     h2bis->SetLineColor(4);
0255     h2bis->Fit("gaus");
0256     gaus->SetLineColor(6);
0257     h2bis->Fit("gaus");
0258 
0259 //**************
0260 // RANGE IN CELL
0261 //**************
0262 
0263 Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2;
0264 
0265 // X position of target in World
0266 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180); 
0267 
0268 // Z position of target in World
0269 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180); 
0270 
0271 // Line alignment (cf MicrobeamEMField.cc)
0272 Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
0273 Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
0274 
0275 TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
0276 Double_t x,y,z,xx,zz;
0277 ntuple2->SetBranchAddress("x",&x);
0278 ntuple2->SetBranchAddress("y",&y);
0279 ntuple2->SetBranchAddress("z",&z);
0280 Int_t nentries = (Int_t)ntuple2->GetEntries();
0281 for (Int_t i=0;i<nentries;i++) 
0282 {
0283       ntuple2->GetEntry(i);
0284       X1=x;
0285       Y1=y;
0286       Z1=z;
0287       xx = X1-Xc;
0288       zz = Z1-Zc;
0289       Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(10*TMath::Pi()/180);
0290       X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(10*TMath::Pi()/180);
0291       Y2 = Y1;    
0292       ntupleR->Fill(Z2,Y2,X2);
0293 }
0294      
0295 c1->cd(10);
0296   ntupleR->Draw("X2:Z2","abs(X2)<50","surf3");
0297   gPad->SetLogz();
0298 
0299 //****************
0300 // ENERGY DEPOSITS 
0301 //****************
0302 
0303 gStyle->SetOptStat(0000);
0304 gStyle->SetOptFit();
0305 gStyle->SetPalette(1);
0306 gROOT->SetStyle("Plain");
0307 
0308 c1->cd(11);
0309   TH2F *histbis = new TH2F("histbis","histbis",50,-20,20,50,-20,20);
0310   ntuple4->Draw("y*0.359060:x*0.359060>>histbis","doseV","contz");
0311   gPad->SetLogz();
0312   histbis->Draw("contz");
0313   histbis->GetXaxis()->SetLabelSize(0.025);
0314   histbis->GetYaxis()->SetLabelSize(0.025);
0315   histbis->GetZaxis()->SetLabelSize(0.025);
0316   histbis->GetXaxis()->SetTitleSize(0.035);
0317   histbis->GetYaxis()->SetTitleSize(0.035);
0318   histbis->GetXaxis()->SetTitleOffset(1.4);
0319   histbis->GetYaxis()->SetTitleOffset(1.4);
0320   histbis->GetXaxis()->SetTitle("Y (um)");
0321   histbis->GetYaxis()->SetTitle("X (um)");
0322   histbis->SetTitle("Energy deposit -transverse- (z axis in eV)");
0323 
0324 c1->cd(12);
0325   TH2F *histter = new TH2F("histter","histter",50,-20,20,50,-20,20);
0326   ntuple4->Draw("x*0.359060:(z+1500/0.162810+21)*0.162810>>histter","doseV","contz");
0327   gPad->SetLogz();
0328   histter->Draw("contz");
0329   histter->GetXaxis()->SetLabelSize(0.025);
0330   histter->GetYaxis()->SetLabelSize(0.025);
0331   histter->GetZaxis()->SetLabelSize(0.025);
0332   histter->GetXaxis()->SetTitleSize(0.035);
0333   histter->GetYaxis()->SetTitleSize(0.035);
0334   histter->GetXaxis()->SetTitleOffset(1.4);
0335   histter->GetYaxis()->SetTitleOffset(1.4);
0336   histter->GetXaxis()->SetTitle("Z (um)");
0337   histter->GetYaxis()->SetTitle("X (um)");
0338   histter->SetTitle("Energy deposit -longitudinal- (z axis in eV)");
0339 
0340 //*******************************
0341 // BEAM POSITION AT CELL ENTRANCE
0342 //*******************************
0343  
0344 gStyle->SetOptStat(0000);
0345 gStyle->SetOptFit();
0346 gStyle->SetPalette(1);
0347 gROOT->SetStyle("Plain");
0348 
0349 TH1F *h77 = new TH1F("hx","h1",200,-10,10); 
0350 TH1F *h88 = new TH1F("hy","h1",200,-10,10); 
0351    
0352 c1->cd(4);
0353         ntuple1->Project("hx","x");
0354     scale = 1/h77->Integral();
0355     h77->Scale(scale);
0356     h77->Draw();
0357     h77->GetXaxis()->SetLabelSize(0.025);
0358     h77->GetYaxis()->SetLabelSize(0.025);
0359     h77->GetXaxis()->SetTitleSize(0.035);
0360     h77->GetYaxis()->SetTitleSize(0.035);
0361     h77->GetXaxis()->SetTitleOffset(1.4);
0362     h77->GetYaxis()->SetTitleOffset(1.4);
0363     h77->GetXaxis()->SetTitle("Position (um)");
0364     h77->GetYaxis()->SetTitle("Fraction of events");
0365     h77->SetTitle("Beam X position on cell");
0366     h77->SetFillColor(4);
0367     h77->SetLineColor(4);
0368     h77->Fit("gaus");
0369 
0370 c1->cd(8);
0371         ntuple1->Project("hy","y");
0372     scale = 1/h88->Integral();
0373     h88->Scale(scale);
0374     h88->Draw();
0375     h88->GetXaxis()->SetLabelSize(0.025);
0376     h88->GetYaxis()->SetLabelSize(0.025);
0377     h88->GetXaxis()->SetTitleSize(0.035);
0378     h88->GetYaxis()->SetTitleSize(0.035);
0379     h88->GetXaxis()->SetTitleOffset(1.4);
0380     h88->GetYaxis()->SetTitleOffset(1.4);
0381     h88->GetXaxis()->SetTitle("Position (um)");
0382     h88->GetYaxis()->SetTitle("Fraction of events");
0383     h88->SetTitle("Beam Y position on cell");
0384     h88->SetFillColor(4);
0385     h88->SetLineColor(4);
0386     h88->Fit("gaus");
0387 }
[ 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!