advanced/exp_microdosimetry/ProcessMicro.C

0001 {
0002 ////////////////////////////////////////////
0003     gROOT->Reset();
0004
0005     bool saveImage = true;
0006     string imageFileName = "microdosimetricSpectra.png";
0007
0008     bool saveLinearEdep = true;
0009     bool saveLogEdep = true;
0010     bool saveydySpec = true;
0011
0012 ////////////////////////////////////////////
0013 //Open File where data has been stored
0014 ////////////////////////////////////////////
0015     TFile* f = new TFile("radioprotection_t0.root");
0016
0017     TDirectory* dir = (TDirectory*)f->Get("radioprotection_ntuple");
0018     TTree * nEdep = (TTree*)dir->Get("102");
0019     Double_t edep;
0020     nEdep->SetBranchAddress("edep", &edep);
0021
0022     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)
0023 ////////////////////////////////////////////
0024 //Put hits in microdosimeter into linear binned space histogram (typical of experimental setup)
0025 ////////////////////////////////////////////
0026     int numberLinearBins = 4096;
0027     double minLinealEnergy = 0;
0028     double maxLinealEnergy = 1000;
0029     TH1D*  histoLin = new TH1D("h1", "f(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
0030
0031     TH1D*  histoYFY = new TH1D("h2", "yf(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
0032     TH1D*  histoDY = new TH1D("h3", "d(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
0033
0034 ////////////////////////////////////////////
0035 //Create log bin for microdosimetric plotting
0036 ////////////////////////////////////////////
0037     //Log bin stuff
0038     int B = 20; //the number of increments in each decade
0039     //eg if B = 10: 1 2 3 ... 9 10 20 30 ... 90 100 200 ...
0040     double yMin = 0.01; // y = lineal energy
0041     double yMax = 1000;
0042
0043     double range = yMax / yMin;
0044     double fooRange = range;
0045     double fooMag = 10; //any old silly no > 1
0046     int order = 0;
0047     while (fooMag > 1)
0048     {
0049         fooMag = fooRange / 10 ;
0050         fooRange = fooRange / 10 ;
0051         order = order + 1 ;
0052     }
0053
0054     int binsPerDecade = B;
0055     const int numberOfLogBins = binsPerDecade * order +1; //get the total number of bins
0056
0057     vector<double> xbins;
0058     for (int b = 0; b < numberOfLogBins; b++)
0059     {
0060         xbins.push_back((yMin)*pow(10, ((b) / double(B))));
0061         //  cout << b << " " << xbins[b] << endl;
0062     }
0063
0064     TH1D*  histoLog = new TH1D("l1", "edepLin", (numberOfLogBins-1), &xbins[0]);
0065
0066
0067 ////////////////////////////////////////////
0068 //loop through the hits in the detector stored in the ntuple
0069 ////////////////////////////////////////////
0070     int numberHits = nEdep->GetEntries();
0071     cout << "Number of hits in detector = " << numberHits << endl;
0072
0073     for (int i = 0; i < numberHits; i++)
0074     {
0075         nEdep->GetEntry(i);
0076         if (edep <= 0)
0077         {
0078             cout << "Edep = 0" << endl;
0079         }
0080
0081         histoLin->Fill(edep/meanChordLength);
0082         histoLog->Fill(edep/meanChordLength);
0083     }
0084
0085 ////////////////////////////////////////////
0086 //Normalise the log histogram based on the bin width
0087 ////////////////////////////////////////////
0088     for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0089     {
0090         histoLog->SetBinContent(bb, (histoLog->GetBinContent(bb) / histoLog->GetBinWidth(bb)));
0091     }
0092
0093
0094 ////////////////////////////////////////////
0095 //---------Save values to file--------------
0096 ////////////////////////////////////////////
0097     ofstream outFile;
0098     if (saveLinearEdep)
0099     {
0100         outFile.open("linEdepSpec.txt");
0101         for (int i = 1 ; i <= numberLinearBins; i++)
0102         {
0103             outFile << histoLin->GetBinCenter(i) * meanChordLength << "\t" << histoLin->GetBinContent(i) << endl;
0104         }
0105         outFile.close();
0106     }
0107     if (saveLogEdep)
0108     {
0109         outFile.open("logEdepSpec.txt");
0110         for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0111         {
0112             outFile << histoLog->GetBinCenter(bb) * meanChordLength << "\t" << histoLog->GetBinContent(bb) << endl;
0113         }
0114         outFile.close();
0115     }
0116
0117 ////////////////////////////////////////////
0118 //----------Normalise Histo to get f(y)----
0119 ////////////////////////////////////////////
0120     double intOfHist = histoLin->Integral("width") ; //gets counts in each histogram and multiplies by bin width
0121     histoLin->Scale(1./intOfHist);
0122     //histoLin->Sumw2();
0123
0124     double intOfLogHist = histoLog->Integral("width") ;
0125     histoLog->Scale(1./intOfLogHist);
0126     //histoLog->Sumw2();
0127
0128
0129 //////////////////////////////////////////////////////////////
0130 //-------Calculate RBE using the modified MK model-----------
0131 /////////////////////////////////////////////////////////////
0132 //Technically the use of the modified MKM is applied in
0133 //therapeutic ion beams and not used for radioprotection applications
0134 //though it is included here for convience.
0135 //NOTE: the form used here is for "heavy ions" and not applicable
0136 //to the more dose depedent protons
0137
0138     //---------------------------------
0139     //Constants
0140     //---------------------------------
0141     const double satPar = 150.; //saturation paramater obv. in keV/um
0142     double satParSqr = satPar * satPar;
0143     const double pi = 3.14159265359;
0144     const double SiToMuscleCorrection = 0.57;
0145
0146     //Cell reltated numbers
0147     double alpha0 = 0.13;//0.164;//0.13; //Gy^-1 //0.13 for carbon
0148     double beta = 0.05; //Gy^-2
0149     double rho = 1.; //g/cm^3
0150     double rSV = 0.42; //um, the radius of the HSG cell
0151     //-----------------------------------------------------------
0152     //In order for consistent units, which in their raw form are:
0153     //(1/(Gy^2))(keV/g)(cm^3/um^3)
0154     // (1/Gy^2)(keV/um)(cm^3/g)(1/um^2)
0155     //------------------------------------------------------------
0156     //So convert (keV/g) into Gy it is scaled by (1.602E-16/0.001)(J/kg)
0157     //And to get cm^3 into um^3 just scale by 10E12
0158     //------------------------------------------------------------
0159     double scaleGray = 1.602E-16/0.001;
0160     double scaleLength = 1.E12;
0161
0162 //-------Calculate y*--------
0163     double topIntegral = 0.;
0164     double botIntegral = 0.;
0165
0166     //top integral
0167     for (Int_t i = 1; i <= numberLinearBins; i++)
0168     {
0169         topIntegral += (1 - exp(-(histoLin->GetBinCenter(i) * histoLin->GetBinCenter(i)) / satParSqr))*histoLin->GetBinContent(i) * histoLin->GetBinWidth(i);
0170         //              (1 - exp(-y^2/y_{0}^2)) * f(y) * dy
0171     }
0172     //bottom integral
0173     for (Int_t i = 1; i <= numberLinearBins; i++)
0174     {
0175         botIntegral += (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i));
0176         //              y (size of each being time bin number) * f(y) * dy (bin size)
0177     }
0178
0179     //cout << "Top: " << topIntegral << endl;
0180     //cout << "Bot: " << botIntegral << endl;
0181
0182     //saturation corrected dose averaged linel energy value (y*)
0183     double satCor = satParSqr * topIntegral / botIntegral;
0184     //(y_{0})^2*top/bottom
0185     //cout << "--------------------------------" << endl;
0186     //cout << "y* = " << satCor << endl;
0187     //cout << "--------------------------------" << endl;
0188
0189
0190 //-----------alpha--------
0191     double alpha = alpha0 + (beta / (rho*pi*rSV*rSV))*satCor*scaleGray*scaleLength;
0192     //cout << "alpha = " << alpha << endl;
0193
0194 //-----------finally RBE itself--------
0195     double ln01 = -2.302585093;
0196     double RBE10;//= (2*beta*5.)((sqrt(alpha*alpha - 4*beta*ln01) - alpha));
0197
0198     double topRBE = (2 * beta*5.);
0199     double botRBE = (sqrt(alpha*alpha - 4 * beta*ln01) - alpha);
0200     RBE10 = topRBE / botRBE;
0201
0202     //cout << "--------------------------------" << endl;
0203     cout << "RBE10 = " << RBE10 << endl;
0204     //cout << "--------------------------------" << endl;
0205
0206 ////////////////////////////////////////////
0207 //-----------Calculate yF-------------------
0208 ////////////////////////////////////////////
0209 //yF = int(y * f(y) dy)
0210     double linyF = 0.;
0211     for (int i = 1 ; i <= numberLinearBins; i++)
0212     {
0213         linyF += histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i) ;
0214     }
0215
0216     cout << "yF = " << linyF << endl;
0217
0218     double logyF = 0;
0219     for (int i = 1 ; i <= numberOfLogBins; i++)
0220     {
0221         logyF += histoLog->GetBinCenter(i) * histoLog->GetBinContent(i) * histoLog->GetBinWidth(i) ;
0222     }
0223     cout << "Log yF = " << logyF << endl;
0224
0225
0226 ////////////////////////////////////////////
0227 //--------------Calculate d(y)--------------
0228 ////////////////////////////////////////////
0229 //d(y) = y*f(y) / yF
0230     for (int i = 1 ; i <= numberLinearBins; i++)
0231     {
0232         histoLin->SetBinContent(i, (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) / linyF));
0233     }
0234
0235     //check that it's normalised to 1
0236     double intOfHistDy = histoLin->Integral("width") ;
0237     cout << "Int of d(y) = " << intOfHistDy << endl;
0238
0239     for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0240     {
0241         histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb) / logyF));
0242     }
0243
0244     double intOfHistLogDy = histoLog->Integral("width") ;
0245     cout << "Int of log d(y) = " << intOfHistLogDy << endl;
0246
0247 ////////////////////////////////////////////
0248 //-------Calculate Average Q(y)-----------
0249 ////////////////////////////////////////////
0250 //Q(y) = (a1/y)[1-exp(-(a2)y*y - (a3)y*y*y)]
0251
0252     //Quality factor values from ICRP 60 report
0253     double qualityCoeff1 = 5510.;
0254     double qualityCoeff2 = 5.E-5;
0255     double qualityCoeff3 = 2.E-7;
0256
0257     double aveQuality = 0.;
0258     //AveQuality = int(Q(y)d(y)dy)
0259
0260     for (int i = 1 ; i <= numberLinearBins; i++)
0261     {
0262         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) ;
0263     }
0264     cout << "<Q> = " << aveQuality << endl;
0265
0266
0267 ////////////////////////////////////////////
0268 //---------Calculate yD--------------------
0269 ////////////////////////////////////////////
0270 //yD = int (y * d(y) dy)
0271     double linyD = 0.;
0272
0273     for (int i = 1 ; i <= numberLinearBins; i++)
0274     {
0275         linyD +=  histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin -> GetBinWidth(i) ;
0276     }
0277
0278     cout << "yD = " << linyD << endl;
0279
0280     double logyD = 0.;
0281
0282     for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0283     {
0284         logyD += histoLog->GetBinCenter(bb) * histoLog->GetBinContent(bb) * histoLog -> GetBinWidth(bb) ;
0285     }
0286     cout << "log yD = " << logyD << endl;
0287
0288 ////////////////////////////////////////////
0289 //---------Create yd(y)-------------------
0290 ////////////////////////////////////////////
0291     for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0292     {
0293         histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb)) );
0294     }
0295     //Scale bins by log(10)/B
0296     //See Appendix B of the ICRU 36 report for more details
0297     histoLog->Scale((log(10) / B));
0298
0299     if (saveydySpec)
0300     {
0301         outFile.open("logYDYspec.txt");
0302         for (int bb = 1; bb <= numberOfLogBins-1; bb++)
0303         {
0304             outFile << histoLog->GetBinCenter(bb) << "\t" << histoLog->GetBinContent(bb) << endl;
0305         }
0306         outFile.close();
0307     }
0308
0309     if (saveImage)
0310     {
0311
0312         TCanvas *c1 = new TCanvas("c1", "", 1000, 800);
0313         //gStyle->SetOptStat(0); //un-comment to remove stat box on top right
0314
0315         histoLog->SetTitle("");
0316         histoLog->GetXaxis()->SetTitle("y (keV/#mum)");
0317         histoLog->GetYaxis()->SetTitle("yd(y)");
0318
0319         histoLog->GetYaxis()->SetLabelFont(42);
0320         histoLog->GetXaxis()->SetLabelFont(42);
0321         histoLog->GetYaxis()->SetTitleFont(42);
0322         histoLog->GetXaxis()->SetTitleFont(42);
0323         histoLog->GetYaxis()->CenterTitle(1);
0324         histoLog->GetXaxis()->CenterTitle(1);
0325         histoLog->GetXaxis()->SetTitleSize(0.045);
0326         histoLog->GetYaxis()->SetTitleSize(0.045);
0327         //histoLog->GetYaxis()->SetRangeUser(0., 1.8);
0328         //histoLog -> GetXaxis()->SetRangeUser(0.01., 1000.);
0329         histoLog->GetYaxis()->SetLabelSize(0.04);
0330         histoLog->GetXaxis()->SetLabelSize(0.04);
0331         histoLog->SetLineColor(1);
0332         histoLog->SetLineWidth(4);
0333         c1->SetLogx();
0334         //c1->SetLogy();
0335         //Plotting misbehaves in root6.XX
0336         histoLog->Draw("");
0337
0338         c1->SaveAs(imageFileName.c_str());
0339     }
0340 }