Ignore:
Timestamp:
Nov 5, 2010, 4:08:39 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/advanced/microbeam/plot.C

    r1230 r1342  
    11// -------------------------------------------------------------------
    2 // $Id: plot.C,v 1.5 2009/04/30 10:23:57 sincerti Exp $
     2// $Id: plot.C,v 1.6 2010/10/07 14:03:11 sincerti Exp $
    33// -------------------------------------------------------------------
    44//
    … …  
    2424gROOT->SetStyle("Plain");
    2525Double_t scale;
    26        
     26
     27
    2728c1 = new TCanvas ("c1","",20,20,1200,900);
    2829c1.Divide(4,3);
    29 
    30 FILE * fp = fopen("dose.txt","r");
    31 Float_t nD,cD;
    32 Int_t ncols=0;
    33 Int_t nlines = 0;
    34 
    35 TH1F *h1  = new TH1F("Absorbed dose distribution in Nucleus","Dose distribution in Nucleus",100,0.001,0.5);
    36 TH1F *h10 = new TH1F("Absorbed dose distribution in Cytoplasm","Dose distribution in Cytoplasm",100,0.001,0.2);
    37 
    38 while (1)
    39    {
    40       ncols = fscanf(fp,"%f %f",&nD,&cD);
    41       if (ncols < 0) break;
    42       h1->Fill(nD);
    43       h10->Fill(cD);
    44       nlines++;
    45    }
    46 fclose(fp);
    47 
    48 c1.cd(2);
    49         scale = 1/h1->Integral();
    50         h1->Scale(scale);
    51         h1->Draw();
    52         h1->GetXaxis()->SetLabelSize(0.025);
    53         h1->GetYaxis()->SetLabelSize(0.025);
    54         h1->GetXaxis()->SetTitleSize(0.035);
    55         h1->GetYaxis()->SetTitleSize(0.035);
    56         h1->GetXaxis()->SetTitleOffset(1.4);
    57         h1->GetYaxis()->SetTitleOffset(1.4);
    58         h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
    59         h1->GetYaxis()->SetTitle("Fraction of events");
    60         h1->SetLineColor(3);
    61         h1->SetFillColor(3);
    62 
    63 //*****************
    64 // DOSE IN CYTOPLASM
    65 //*****************
    66 
    67 c1.cd(6);
    68         scale = 1/h10->Integral();
    69         h10->Scale(scale);
    70         h10->Draw();
    71         h10->GetXaxis()->SetLabelSize(0.025);
    72         h10->GetYaxis()->SetLabelSize(0.025);
    73         h10->GetXaxis()->SetTitleSize(0.035);
    74         h10->GetYaxis()->SetTitleSize(0.035);
    75         h10->GetXaxis()->SetTitleOffset(1.4);
    76         h10->GetYaxis()->SetTitleOffset(1.4);
    77         h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
    78         h10->GetYaxis()->SetTitle("Fraction of events");
    79         h10->SetLineColor(2);
    80         h10->SetFillColor(2);
    81 
    82 //********************************
    83 // STOPPING POWER AT CELL ENTRANCE
    84 //********************************
    85  
    86 gStyle->SetOptStat(0000);
    87 gStyle->SetOptFit();
    88 gStyle->SetPalette(1);
    89 gROOT->SetStyle("Plain");
    90 
    91 Float_t d;
    92 FILE * fp = fopen("stoppingPower.txt","r");
    93 
    94 TH1F *h2 = new TH1F("Beam stopping Power at cell entrance","h1",200,0,300);
    95 while (1)
    96    {
    97       ncols = fscanf(fp,"%f",&d);
    98       if (ncols < 0) break;
    99       h2->Fill(d);
    100       nlines++;
    101    }
    102 fclose(fp);
    103    
    104 c1.cd(9);
    105         scale = 1/h2->Integral();
    106         h2->Scale(scale);
    107         h2->Draw();
    108         h2->GetXaxis()->SetLabelSize(0.025);
    109         h2->GetYaxis()->SetLabelSize(0.025);
    110         h2->GetXaxis()->SetTitleSize(0.035);
    111         h2->GetYaxis()->SetTitleSize(0.035);
    112         h2->GetXaxis()->SetTitleOffset(1.4);
    113         h2->GetYaxis()->SetTitleOffset(1.4);
    114         h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
    115         h2->GetYaxis()->SetTitle("Fraction of events");
    116         h2->SetTitle("dE/dx at cell entrance");
    117         h2->SetFillColor(4);
    118         h2->SetLineColor(4);
    119         h2->Fit("gaus");
    120         gaus->SetLineColor(6);
    121         h2->Fit("gaus");
    122 
    123 //**************
    124 // RANGE IN CELL
    125 //**************
    126 
    127 Float_t Xc,Zc,X1,Y1,Z1,x,z,X2,Y2,Z2;
    128 Float_t d;
    129 
    130 // X position of target in World
    131 Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
    132 
    133 // Z position of target in World
    134 Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
    135 
    136 // Line alignment (cf MicrobeamEMField.cc)
    137 Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
    138 Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
    139 
    140 FILE * fp = fopen("range.txt","r");
    141 
    142 TNtuple *ntuple = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
    143 
    144 while (1)
    145    {
    146       ncols = fscanf(fp,"%f %f %f",&X1,&Y1,&Z1);
    147       if (ncols < 0) break;
    148       x = X1-Xc;
    149       z = Z1-Zc;
    150       Z2 = z*cos(10*TMath::Pi()/180)-x*sin(10*TMath::Pi()/180);
    151       X2 = z*sin(10*TMath::Pi()/180)+x*cos(10*TMath::Pi()/180);
    152       Y2 = Y1;
    153      
    154       ntuple->Fill(Z2,Y2,X2);
    155       nlines++;
    156    }
    157 fclose(fp);
    158      
    159 c1.cd(10);
    160   ntuple->Draw("X2:Z2","abs(X2)<50","surf3");
    161   gPad->SetLogz();
    162   htemp->GetXaxis()->SetLabelSize(0.025);
    163   htemp->GetYaxis()->SetLabelSize(0.025);
    164   htemp->GetZaxis()->SetLabelSize(0.025);
    165   htemp->GetXaxis()->SetTitleSize(0.035);
    166   htemp->GetYaxis()->SetTitleSize(0.035);
    167   htemp->GetXaxis()->SetTitleOffset(1.4);
    168   htemp->GetYaxis()->SetTitleOffset(1.4);
    169   htemp->GetXaxis()->SetTitle("Z (µm)");
    170   htemp->GetYaxis()->SetTitle("X (µm)");
    171   htemp->SetTitle("Range in cell");
    17230
    17331//*********************
    17432// INTENSITY HISTOGRAMS
    17533//*********************
     34jump:
    17635
    17736FILE * fp = fopen("phantom.dat","r");
    … …  
    18948
    19049TNtuple *ntupleYXN = new TNtuple("NUCLEUS","ntuple","Y:X:vox");
    191 TNtuple *ntupleZX = new TNtuple("CYTOPASM","ntuple","Z:X:vox");
    192 TNtuple *ntupleYX = new TNtuple("CYTOPASM","ntuple","Y:X:vox");
     50TNtuple *ntupleZX = new TNtuple("CYTOPLASM","ntuple","Z:X:vox");
     51TNtuple *ntupleYX = new TNtuple("CYTOPLASM","ntuple","Y:X:vox");
    19352
    19453nlines=0;
    … …  
    308167  hist->SetTitle("Nucleus intensity on transverse section");
    309168
    310 //****************
    311 // ENERGY DEPOSITS
    312 //****************
    313 
     169//
     170
     171TFile f("microbeam.root");
     172
     173TNtuple* ntuple0;
     174TNtuple* ntuple1;
     175TNtuple* ntuple2;
     176TNtuple* ntuple3;
     177TNtuple* ntuple4;
     178
     179ntuple0 = (TNtuple*)f->Get("ntuple0");
     180ntuple1 = (TNtuple*)f->Get("ntuple1");
     181ntuple2 = (TNtuple*)f->Get("ntuple2");
     182ntuple3 = (TNtuple*)f->Get("ntuple3");
     183ntuple4 = (TNtuple*)f->Get("ntuple4");
     184
     185TH1F *h1  = new TH1F("h1","Dose distribution in Nucleus",100,0.001,1.);
     186TH1F *h10 = new TH1F("h10","Dose distribution in Cytoplasm",100,0.001,.2);
     187
     188c1.cd(2);
     189
     190        ntuple3->Project("h1","doseN");
     191        scale = 1/h1->Integral();
     192        h1->Scale(scale);
     193        h1->Draw();
     194        h1->GetXaxis()->SetLabelSize(0.025);
     195        h1->GetYaxis()->SetLabelSize(0.025);
     196        h1->GetXaxis()->SetTitleSize(0.035);
     197        h1->GetYaxis()->SetTitleSize(0.035);
     198        h1->GetXaxis()->SetTitleOffset(1.4);
     199        h1->GetYaxis()->SetTitleOffset(1.4);
     200        h1->GetXaxis()->SetTitle("Absorbed dose (Gy)");
     201        h1->GetYaxis()->SetTitle("Fraction of events");
     202        h1->SetLineColor(3);
     203        h1->SetFillColor(3);
     204
     205//*****************
     206// DOSE IN CYTOPLASM
     207//*****************
     208
     209c1.cd(6);
     210        ntuple3->Project("h10","doseC");
     211        scale = 1/h10->Integral();
     212        h10->Scale(scale);
     213        h10->Draw();
     214        h10->GetXaxis()->SetLabelSize(0.025);
     215        h10->GetYaxis()->SetLabelSize(0.025);
     216        h10->GetXaxis()->SetTitleSize(0.035);
     217        h10->GetYaxis()->SetTitleSize(0.035);
     218        h10->GetXaxis()->SetTitleOffset(1.4);
     219        h10->GetYaxis()->SetTitleOffset(1.4);
     220        h10->GetXaxis()->SetTitle("Absorbed dose (Gy)");
     221        h10->GetYaxis()->SetTitle("Fraction of events");
     222        h10->SetLineColor(2);
     223        h10->SetFillColor(2);
     224
     225//********************************
     226// STOPPING POWER AT CELL ENTRANCE
     227//********************************
     228 
    314229gStyle->SetOptStat(0000);
    315230gStyle->SetOptFit();
    … …  
    317232gROOT->SetStyle("Plain");
    318233
    319 FILE * fp = fopen("3DDose.txt","r");
    320 
    321 TNtuple *ntuple2 = new TNtuple("CELL","ntuple","yVox:xVox:dose");
    322 TNtuple *ntuple3 = new TNtuple("CELL","ntuple","xVox:zVox:dose");
    323 TNtuple *ntuplezyx = new TNtuple("DOSE","ntuple","zVox:yVox:xVox:dose");
    324 
    325 while (1)
    326    {
    327       ncols = fscanf(fp,"%f %f %f %f",&xVox, &yVox, &zVox, &dose);
    328       if (ncols < 0) break;
    329      
    330       xVox= xVox*voxelSizeX;
    331       yVox= yVox*voxelSizeY;
    332       zVox= zVox*voxelSizeZ;
    333      
    334       ntuple2->Fill(yVox,xVox,dose);
    335       ntuple3->Fill(xVox,zVox,dose);
    336       ntuplezyx->Fill(zVox,yVox,xVox,dose);
    337    }
    338 fclose(fp);
     234Float_t d;
     235
     236TH1F *h2 = new TH1F("h2","Beam stopping power at cell entrance",200,0,300);
     237   
     238c1.cd(9);
     239        ntuple0->Project("h2","sp");
     240        scale = 1/h2->Integral();
     241        h2->Scale(scale);
     242        h2->Draw();
     243        h2->GetXaxis()->SetLabelSize(0.025);
     244        h2->GetYaxis()->SetLabelSize(0.025);
     245        h2->GetXaxis()->SetTitleSize(0.035);
     246        h2->GetYaxis()->SetTitleSize(0.035);
     247        h2->GetXaxis()->SetTitleOffset(1.4);
     248        h2->GetYaxis()->SetTitleOffset(1.4);
     249        h2->GetXaxis()->SetTitle("dE/dx (keV/µm)");
     250        h2->GetYaxis()->SetTitle("Fraction of events");
     251        h2->SetTitle("dE/dx at cell entrance");
     252        h2->SetFillColor(4);
     253        h2->SetLineColor(4);
     254        h2->Fit("gaus");
     255        gaus->SetLineColor(6);
     256        h2->Fit("gaus");
     257
     258
     259//**************
     260// RANGE IN CELL
     261//**************
     262
     263Double_t Xc,Zc,X1,Y1,Z1,X2,Y2,Z2;
     264
     265// X position of target in World
     266Xc = -1295.59e3 - 955e3*sin(10*TMath::Pi()/180);
     267
     268// Z position of target in World
     269Zc = -1327e3 + 955e3*cos(10*TMath::Pi()/180);
     270
     271// Line alignment (cf MicrobeamEMField.cc)
     272Xc = Xc + 5.24*cos(10*TMath::Pi()/180);
     273Zc = Zc + 5.24*sin(10*TMath::Pi()/180);
     274
     275TNtuple *ntupleR = new TNtuple("Rmax","ntuple","Z2:Y2:X2");
     276Double_t x,y,z,xx,zz;
     277ntuple2->SetBranchAddress("x",&x);
     278ntuple2->SetBranchAddress("y",&y);
     279ntuple2->SetBranchAddress("z",&z);
     280Int_t nentries = (Int_t)ntuple2->GetEntries();
     281for (Int_t i=0;i<nentries;i++)
     282{
     283      ntuple2->GetEntry(i);
     284      X1=x;
     285      Y1=y;
     286      Z1=z;
     287      xx = X1-Xc;
     288      zz = Z1-Zc;
     289      Z2 = zz*cos(10*TMath::Pi()/180)-xx*sin(10*TMath::Pi()/180);
     290      X2 = zz*sin(10*TMath::Pi()/180)+xx*cos(10*TMath::Pi()/180);
     291      Y2 = Y1;   
     292      ntupleR->Fill(Z2,Y2,X2);
     293}
     294     
     295c1.cd(10);
     296  ntupleR->Draw("X2:Z2","abs(X2)<50","surf3");
     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  */
     310
     311//****************
     312// ENERGY DEPOSITS
     313//****************
     314
     315
     316gStyle->SetOptStat(0000);
     317gStyle->SetOptFit();
     318gStyle->SetPalette(1);
     319gROOT->SetStyle("Plain");
    339320
    340321c1.cd(11);
    341322  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
    342   ntuple2->Draw("yVox:xVox>>hist","dose","contz");
     323  ntuple4->Draw("y*0.359060:x*0.359060>>hist","doseV","contz");
    343324  gPad->SetLogz();
    344325  hist->Draw("contz");
    … …  
    353334  hist->GetYaxis()->SetTitle("X (µm)");
    354335  hist->SetTitle("Mean energy deposit -transverse- (z axis in eV)");
     336
    355337 
    356338c1.cd(12);
    357339  TH2F *hist = new TH2F("hist","hist",50,-20,20,50,-20,20);
    358   ntuple3->Draw("xVox:zVox>>hist","dose","contz");
     340  ntuple4->Draw("x*0.359060:z*0.162810>>hist","doseV","contz");
    359341  gPad->SetLogz();
    360342  hist->Draw("contz");
    … …  
    379361gROOT->SetStyle("Plain");
    380362
    381 Float_t bx, by;
    382 FILE * fp = fopen("beamPosition.txt","r");
    383 
    384 TH1F *h77 = new TH1F("Beam X transverse position at cell entrance","h1",200,-10,10);
    385 TH1F *h88 = new TH1F("Beam Y transverse position at cell entrance","h1",200,-10,10);
    386 while (1)
    387    {
    388       ncols = fscanf(fp,"%f %f",&bx, &by);
    389       if (ncols < 0) break;
    390       h77->Fill(bx);
    391       h88->Fill(by);
    392       nlines++;
    393     }
    394 fclose(fp);
    395    
     363TH1F *h77 = new TH1F("hx","h1",200,-10,10);
     364TH1F *h88 = new TH1F("hy","h1",200,-10,10);
     365   
    396366c1.cd(4);
    397         scale = 1/h77->Integral();
     367        ntuple1->Project("hx","x");
     368        scale = 1/h77->Integral();
    398369        h77->Scale(scale);
    399370        h77->Draw();
    … …  
    413384
    414385c1.cd(8);
    415         scale = 1/h88->Integral();
     386        ntuple1->Project("hy","y");
     387        scale = 1/h88->Integral();
    416388        h88->Scale(scale);
    417389        h88->Draw();
Note: See TracChangeset for help on using the changeset viewer.