Gabriele Cosmo
341 lines
11 KiB
C++
341 lines
11 KiB
C++
{
|
|
////////////////////////////////////////////
|
|
gROOT->Reset();
|
|
|
|
bool saveImage = true;
|
|
string imageFileName = "microdosimetricSpectra.png";
|
|
|
|
bool saveLinearEdep = true;
|
|
bool saveLogEdep = true;
|
|
bool saveydySpec = true;
|
|
|
|
////////////////////////////////////////////
|
|
//Open File where data has been stored
|
|
////////////////////////////////////////////
|
|
TFile* f = new TFile("radioprotection_t0.root");
|
|
|
|
TDirectory* dir = (TDirectory*)f->Get("radioprotection_ntuple");
|
|
TTree * nEdep = (TTree*)dir->Get("102");
|
|
Double_t edep;
|
|
nEdep->SetBranchAddress("edep", &edep);
|
|
|
|
double meanChordLength = 1; //change based on the thickness of the sensitive volume (make sure to apply a conversion factor if converting to a tissue equivalent material)
|
|
////////////////////////////////////////////
|
|
//Put hits in microdosimeter into linear binned space histogram (typical of experimental setup)
|
|
////////////////////////////////////////////
|
|
int numberLinearBins = 4096;
|
|
double minLinealEnergy = 0;
|
|
double maxLinealEnergy = 1000;
|
|
TH1D* histoLin = new TH1D("h1", "f(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
|
|
|
|
TH1D* histoYFY = new TH1D("h2", "yf(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
|
|
TH1D* histoDY = new TH1D("h3", "d(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
|
|
|
|
////////////////////////////////////////////
|
|
//Create log bin for microdosimetric plotting
|
|
////////////////////////////////////////////
|
|
//Log bin stuff
|
|
int B = 20; //the number of increments in each decade
|
|
//eg if B = 10: 1 2 3 ... 9 10 20 30 ... 90 100 200 ...
|
|
double yMin = 0.01; // y = lineal energy
|
|
double yMax = 1000;
|
|
|
|
double range = yMax / yMin;
|
|
double fooRange = range;
|
|
double fooMag = 10; //any old silly no > 1
|
|
int order = 0;
|
|
while (fooMag > 1)
|
|
{
|
|
fooMag = fooRange / 10 ;
|
|
fooRange = fooRange / 10 ;
|
|
order = order + 1 ;
|
|
}
|
|
|
|
int binsPerDecade = B;
|
|
const int numberOfLogBins = binsPerDecade * order +1; //get the total number of bins
|
|
|
|
vector<double> xbins;
|
|
for (int b = 0; b < numberOfLogBins; b++)
|
|
{
|
|
xbins.push_back((yMin)*pow(10, ((b) / double(B))));
|
|
// cout << b << " " << xbins[b] << endl;
|
|
}
|
|
|
|
TH1D* histoLog = new TH1D("l1", "edepLin", (numberOfLogBins-1), &xbins[0]);
|
|
|
|
|
|
////////////////////////////////////////////
|
|
//loop through the hits in the detector stored in the ntuple
|
|
////////////////////////////////////////////
|
|
int numberHits = nEdep->GetEntries();
|
|
cout << "Number of hits in detector = " << numberHits << endl;
|
|
|
|
for (int i = 0; i < numberHits; i++)
|
|
{
|
|
nEdep->GetEntry(i);
|
|
if (edep <= 0)
|
|
{
|
|
cout << "Edep = 0" << endl;
|
|
}
|
|
|
|
histoLin->Fill(edep/meanChordLength);
|
|
histoLog->Fill(edep/meanChordLength);
|
|
}
|
|
|
|
////////////////////////////////////////////
|
|
//Normalise the log histogram based on the bin width
|
|
////////////////////////////////////////////
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
histoLog->SetBinContent(bb, (histoLog->GetBinContent(bb) / histoLog->GetBinWidth(bb)));
|
|
}
|
|
|
|
|
|
////////////////////////////////////////////
|
|
//---------Save values to file--------------
|
|
////////////////////////////////////////////
|
|
ofstream outFile;
|
|
if (saveLinearEdep)
|
|
{
|
|
outFile.open("linEdepSpec.txt");
|
|
for (int i = 1 ; i <= numberLinearBins; i++)
|
|
{
|
|
outFile << histoLin->GetBinCenter(i) * meanChordLength << "\t" << histoLin->GetBinContent(i) << endl;
|
|
}
|
|
outFile.close();
|
|
}
|
|
if (saveLogEdep)
|
|
{
|
|
outFile.open("logEdepSpec.txt");
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
outFile << histoLog->GetBinCenter(bb) * meanChordLength << "\t" << histoLog->GetBinContent(bb) << endl;
|
|
}
|
|
outFile.close();
|
|
}
|
|
|
|
////////////////////////////////////////////
|
|
//----------Normalise Histo to get f(y)----
|
|
////////////////////////////////////////////
|
|
double intOfHist = histoLin->Integral("width") ; //gets counts in each histogram and multiplies by bin width
|
|
histoLin->Scale(1./intOfHist);
|
|
//histoLin->Sumw2();
|
|
|
|
double intOfLogHist = histoLog->Integral("width") ;
|
|
histoLog->Scale(1./intOfLogHist);
|
|
//histoLog->Sumw2();
|
|
|
|
|
|
//////////////////////////////////////////////////////////////
|
|
//-------Calculate RBE using the modified MK model-----------
|
|
/////////////////////////////////////////////////////////////
|
|
//Technically the use of the modified MKM is applied in
|
|
//therapeutic ion beams and not used for radioprotection applications
|
|
//though it is included here for convience.
|
|
//NOTE: the form used here is for "heavy ions" and not applicable
|
|
//to the more dose depedent protons
|
|
|
|
//---------------------------------
|
|
//Constants
|
|
//---------------------------------
|
|
const double satPar = 150.; //saturation paramater obv. in keV/um
|
|
double satParSqr = satPar * satPar;
|
|
const double pi = 3.14159265359;
|
|
const double SiToMuscleCorrection = 0.57;
|
|
|
|
//Cell reltated numbers
|
|
double alpha0 = 0.13;//0.164;//0.13; //Gy^-1 //0.13 for carbon
|
|
double beta = 0.05; //Gy^-2
|
|
double rho = 1.; //g/cm^3
|
|
double rSV = 0.42; //um, the radius of the HSG cell
|
|
//-----------------------------------------------------------
|
|
//In order for consistent units, which in their raw form are:
|
|
//(1/(Gy^2))(keV/g)(cm^3/um^3)
|
|
// (1/Gy^2)(keV/um)(cm^3/g)(1/um^2)
|
|
//------------------------------------------------------------
|
|
//So convert (keV/g) into Gy it is scaled by (1.602E-16/0.001)(J/kg)
|
|
//And to get cm^3 into um^3 just scale by 10E12
|
|
//------------------------------------------------------------
|
|
double scaleGray = 1.602E-16/0.001;
|
|
double scaleLength = 1.E12;
|
|
|
|
//-------Calculate y*--------
|
|
double topIntegral = 0.;
|
|
double botIntegral = 0.;
|
|
|
|
//top integral
|
|
for (Int_t i = 1; i <= numberLinearBins; i++)
|
|
{
|
|
topIntegral += (1 - exp(-(histoLin->GetBinCenter(i) * histoLin->GetBinCenter(i)) / satParSqr))*histoLin->GetBinContent(i) * histoLin->GetBinWidth(i);
|
|
// (1 - exp(-y^2/y_{0}^2)) * f(y) * dy
|
|
}
|
|
//bottom integral
|
|
for (Int_t i = 1; i <= numberLinearBins; i++)
|
|
{
|
|
botIntegral += (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i));
|
|
// y (size of each being time bin number) * f(y) * dy (bin size)
|
|
}
|
|
|
|
//cout << "Top: " << topIntegral << endl;
|
|
//cout << "Bot: " << botIntegral << endl;
|
|
|
|
//saturation corrected dose averaged linel energy value (y*)
|
|
double satCor = satParSqr * topIntegral / botIntegral;
|
|
//(y_{0})^2*top/bottom
|
|
//cout << "--------------------------------" << endl;
|
|
//cout << "y* = " << satCor << endl;
|
|
//cout << "--------------------------------" << endl;
|
|
|
|
|
|
//-----------alpha--------
|
|
double alpha = alpha0 + (beta / (rho*pi*rSV*rSV))*satCor*scaleGray*scaleLength;
|
|
//cout << "alpha = " << alpha << endl;
|
|
|
|
//-----------finally RBE itself--------
|
|
double ln01 = -2.302585093;
|
|
double RBE10;//= (2*beta*5.)((sqrt(alpha*alpha - 4*beta*ln01) - alpha));
|
|
|
|
double topRBE = (2 * beta*5.);
|
|
double botRBE = (sqrt(alpha*alpha - 4 * beta*ln01) - alpha);
|
|
RBE10 = topRBE / botRBE;
|
|
|
|
//cout << "--------------------------------" << endl;
|
|
cout << "RBE10 = " << RBE10 << endl;
|
|
//cout << "--------------------------------" << endl;
|
|
|
|
////////////////////////////////////////////
|
|
//-----------Calculate yF-------------------
|
|
////////////////////////////////////////////
|
|
//yF = int(y * f(y) dy)
|
|
double linyF = 0.;
|
|
for (int i = 1 ; i <= numberLinearBins; i++)
|
|
{
|
|
linyF += histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i) ;
|
|
}
|
|
|
|
cout << "yF = " << linyF << endl;
|
|
|
|
double logyF = 0;
|
|
for (int i = 1 ; i <= numberOfLogBins; i++)
|
|
{
|
|
logyF += histoLog->GetBinCenter(i) * histoLog->GetBinContent(i) * histoLog->GetBinWidth(i) ;
|
|
}
|
|
cout << "Log yF = " << logyF << endl;
|
|
|
|
|
|
////////////////////////////////////////////
|
|
//--------------Calculate d(y)--------------
|
|
////////////////////////////////////////////
|
|
//d(y) = y*f(y) / yF
|
|
for (int i = 1 ; i <= numberLinearBins; i++)
|
|
{
|
|
histoLin->SetBinContent(i, (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) / linyF));
|
|
}
|
|
|
|
//check that it's normalised to 1
|
|
double intOfHistDy = histoLin->Integral("width") ;
|
|
cout << "Int of d(y) = " << intOfHistDy << endl;
|
|
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb) / logyF));
|
|
}
|
|
|
|
double intOfHistLogDy = histoLog->Integral("width") ;
|
|
cout << "Int of log d(y) = " << intOfHistLogDy << endl;
|
|
|
|
////////////////////////////////////////////
|
|
//-------Calculate Average Q(y)-----------
|
|
////////////////////////////////////////////
|
|
//Q(y) = (a1/y)[1-exp(-(a2)y*y - (a3)y*y*y)]
|
|
|
|
//Quality factor values from ICRP 60 report
|
|
double qualityCoeff1 = 5510.;
|
|
double qualityCoeff2 = 5.E-5;
|
|
double qualityCoeff3 = 2.E-7;
|
|
|
|
double aveQuality = 0.;
|
|
//AveQuality = int(Q(y)d(y)dy)
|
|
|
|
for (int i = 1 ; i <= numberLinearBins; i++)
|
|
{
|
|
aveQuality += (qualityCoeff1/histoLin->GetBinCenter(i))* (1 - exp(-qualityCoeff2*pow(histoLin->GetBinCenter(i), 2.)-qualityCoeff3*pow(histoLin->GetBinCenter(i), 3.)))* histoLin->GetBinContent(i)* histoLin->GetBinWidth(i) ;
|
|
}
|
|
cout << "<Q> = " << aveQuality << endl;
|
|
|
|
|
|
////////////////////////////////////////////
|
|
//---------Calculate yD--------------------
|
|
////////////////////////////////////////////
|
|
//yD = int (y * d(y) dy)
|
|
double linyD = 0.;
|
|
|
|
for (int i = 1 ; i <= numberLinearBins; i++)
|
|
{
|
|
linyD += histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin -> GetBinWidth(i) ;
|
|
}
|
|
|
|
cout << "yD = " << linyD << endl;
|
|
|
|
double logyD = 0.;
|
|
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
logyD += histoLog->GetBinCenter(bb) * histoLog->GetBinContent(bb) * histoLog -> GetBinWidth(bb) ;
|
|
}
|
|
cout << "log yD = " << logyD << endl;
|
|
|
|
////////////////////////////////////////////
|
|
//---------Create yd(y)-------------------
|
|
////////////////////////////////////////////
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb)) );
|
|
}
|
|
//Scale bins by log(10)/B
|
|
//See Appendix B of the ICRU 36 report for more details
|
|
histoLog->Scale((log(10) / B));
|
|
|
|
if (saveydySpec)
|
|
{
|
|
outFile.open("logYDYspec.txt");
|
|
for (int bb = 1; bb <= numberOfLogBins-1; bb++)
|
|
{
|
|
outFile << histoLog->GetBinCenter(bb) << "\t" << histoLog->GetBinContent(bb) << endl;
|
|
}
|
|
outFile.close();
|
|
}
|
|
|
|
if (saveImage)
|
|
{
|
|
|
|
TCanvas *c1 = new TCanvas("c1", "", 1000, 800);
|
|
//gStyle->SetOptStat(0); //un-comment to remove stat box on top right
|
|
|
|
histoLog->SetTitle("");
|
|
histoLog->GetXaxis()->SetTitle("y (keV/#mum)");
|
|
histoLog->GetYaxis()->SetTitle("yd(y)");
|
|
|
|
histoLog->GetYaxis()->SetLabelFont(42);
|
|
histoLog->GetXaxis()->SetLabelFont(42);
|
|
histoLog->GetYaxis()->SetTitleFont(42);
|
|
histoLog->GetXaxis()->SetTitleFont(42);
|
|
histoLog->GetYaxis()->CenterTitle(1);
|
|
histoLog->GetXaxis()->CenterTitle(1);
|
|
histoLog->GetXaxis()->SetTitleSize(0.045);
|
|
histoLog->GetYaxis()->SetTitleSize(0.045);
|
|
//histoLog->GetYaxis()->SetRangeUser(0., 1.8);
|
|
//histoLog -> GetXaxis()->SetRangeUser(0.01., 1000.);
|
|
histoLog->GetYaxis()->SetLabelSize(0.04);
|
|
histoLog->GetXaxis()->SetLabelSize(0.04);
|
|
histoLog->SetLineColor(1);
|
|
histoLog->SetLineWidth(4);
|
|
c1->SetLogx();
|
|
//c1->SetLogy();
|
|
//Plotting misbehaves in root6.XX
|
|
histoLog->Draw("");
|
|
|
|
c1->SaveAs(imageFileName.c_str());
|
|
}
|
|
}
|