// Edited by Dhruval Shah (Summer 2021) // // function prototypes void drawDash (int counter); // draws Dashes (-) void convertAxisAndDrawX (TH1I* hist_in); // convert X-axis from channel number to Energy // main function void Si22DecayProton5 () { // global constants const Int_t peaksFound = 32; // total number of peaks const Int_t peaksFound_g = 32; // finds directory of file TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName()); dir.ReplaceAll("Si22DecayProton5.C",""); dir.ReplaceAll("/./","/"); ifstream inData; TString inputFileName = Form("%sRawData/Si22/corr.dmax.2!decay.dat", dir.Data()); inData.open(inputFileName); if(!(inData.is_open())) { cout << "!!!Error!!! File Opening Failed" << endl; } TCanvas* c1 = new TCanvas("c1", "22Si Protons", 1300, 1000); c1->ToggleEventStatus(); //c1->ToggleEditor(); //c1->Divide(1,2); /* Load Histogram Data */ //c1->cd(1); TH1I* h1 = new TH1I("h1", "22Si Proton Spectrum; Channel Number; Number of Counts", 256, 0.5, 256.5); int channelNumber; int loopCounter = 1000; double countNumber; // inputs data to histogram while(inData) { inData >> channelNumber >> countNumber; // corrects for last channel duplication if(channelNumber == loopCounter) { break; } loopCounter = channelNumber; h1->SetBinContent(channelNumber, countNumber); } inData.close(); h1->SetAxisRange(0.5, 256.5); h1->Draw(); /* Find Peaks and Fit */ //c1->cd(2); TH1I* h2 = (TH1I*)h1->Clone("h2"); // fit results variables double constants[peaksFound], error_constants[peaksFound]; double means[peaksFound], error_means[peaksFound]; double gaussianArea[peaksFound]; TF1* skew[peaksFound]; double peakPosition[peaksFound] = {17,24, 30,35,30,30, 45,50, 60,60,60, 80,80,80, 80,80,80, 90,90,90, 100,100, 105, 110, 115, 120, 120, 130, 130, 140, 154, 154}; double leftRange[peaksFound] = {10,19.5, 10,20,20,20, 35,40, 45,50,50, 50,50,50, 50,50,50, 75,75,75, 80,80, 90, 100, 100, 100, 110, 100, 120, 130, 135, 140}; double rightRange[peaksFound] = {30,33.5, 50,50,50,50, 55,60, 70,70,70, 90,90,90, 90,90,90, 105,105,105, 110,110, 120, 120, 130, 140, 140, 140, 140, 155, 165, 165}; for(int i = 0; i < peaksFound; i++) { // skewed gaussians and default parameters // we are guessing the function here to fit the data skew[i] = new TF1("skew", "[0] * exp(-0.5*((x-[1])/([2]+(x<[1])*[3]*(x-[1])))^2)"); // (funcName, funcFormula) skew[i]->SetParNames("Height", "Centroid", "Sigma", "Skewness"); // 4 free parametersf skew[i]->SetParLimits(0, 1, 7500); skew[i]->SetParameter(0, 3500); skew[i]->SetParLimits(1, 10, 160); skew[i]->SetParameter(1, 100); skew[i]->SetParLimits(2, 1, 3); skew[i]->FixParameter(2, 2.03945); skew[i]->SetParLimits(3, 0.1, 0.3); skew[i]->FixParameter(3, 0.199132); // Fit Each peak appropriately skew[i]->SetParameter(1, peakPosition[i]); // if(i == 0) // Actual Fix // { // skew[i]->FixParameter(1, 17.6875); // skew[i]->FixParameter(0, 436.632); // } if(i == 0) // Actual Fix // 419.5 { skew[i]->FixParameter(1, 17.6875); skew[i]->FixParameter(0, 405); //skew[i]->FixParameter(2, 3); skew[i]->FixParameter(3, 0.45); } if(i == 1) { skew[i]->FixParameter(1, 24.603); // Actual Fix skew[i]->FixParameter(0, 1072.16); } if(i == 2) { skew[i]->FixParameter(1, 30.0); skew[i]->FixParameter(0, 130); } if(i == 3) { skew[i]->FixParameter(1, 34.9); skew[i]->FixParameter(0, 120); } if(i == 4) { skew[i]->FixParameter(1, 38.8); // 36 skew[i]->FixParameter(0, 180); //140 } if(i == 5) { skew[i]->FixParameter(1, 41.8); skew[i]->FixParameter(0, 520); //480 } if(i == 6) // Actual Fix { skew[i]->FixParameter(1, 45.2752); skew[i]->FixParameter(0, 7050); } if(i == 7) { skew[i]->FixParameter(1, 49.5); skew[i]->FixParameter(0, 2050); } if(i == 8) { skew[i]->FixParameter(1, 55); skew[i]->FixParameter(0, 260); } if(i == 9) { skew[i]->FixParameter(1, 59.2); skew[i]->FixParameter(0, 135); } if(i == 10) { skew[i]->FixParameter(1, 63); skew[i]->FixParameter(0, 120); } if(i == 11) { skew[i]->FixParameter(1, 66.7); skew[i]->FixParameter(0, 240); } if(i == 12) { skew[i]->FixParameter(1, 69.5); skew[i]->FixParameter(0, 220); } if(i == 13) { skew[i]->FixParameter(1, 73.7); skew[i]->FixParameter(0, 300); } if(i == 14) { skew[i]->FixParameter(1, 77.5); skew[i]->FixParameter(0, 330); } if(i == 15) { skew[i]->FixParameter(1, 80.5); skew[i]->FixParameter(0, 270); } if(i == 16) { skew[i]->FixParameter(1, 83.7); skew[i]->FixParameter(0, 215); } if(i == 17) { skew[i]->FixParameter(1, 87.4); skew[i]->FixParameter(0, 250); } if(i == 18) { skew[i]->FixParameter(1, 91); skew[i]->FixParameter(0, 130); } if(i == 19) { skew[i]->FixParameter(1, 93); skew[i]->FixParameter(0, 140); } if(i == 20) { skew[i]->FixParameter(1, 96.2); skew[i]->FixParameter(0, 130); } if(i == 21) { skew[i]->FixParameter(1, 99.4); skew[i]->FixParameter(0, 90); } if(i == 22) { skew[i]->FixParameter(1, 103.993); skew[i]->FixParameter(0, 320.802); } if(i == 23) { skew[i]->FixParameter(1, 109); skew[i]->FixParameter(0, 80); } if(i == 24) { skew[i]->FixParameter(1, 113.7); skew[i]->FixParameter(0, 43); } if(i == 25) { skew[i]->FixParameter(1, 118.61); skew[i]->FixParameter(0, 64); } if(i == 26) { skew[i]->FixParameter(1, 123.5); skew[i]->FixParameter(0, 20); } if(i == 27) { skew[i]->FixParameter(1, 128.221); skew[i]->FixParameter(0, 30); } if(i == 28) { skew[i]->FixParameter(1, 132); skew[i]->FixParameter(0, 15); } if(i == 29) { skew[i]->FixParameter(1, 138); skew[i]->FixParameter(0, 5); } if(i == 30) { skew[i]->FixParameter(1, 147); skew[i]->FixParameter(0, 4); } if(i == 31) { skew[i]->FixParameter(1, 152.5); skew[i]->FixParameter(0, 7.5); } h2->Fit("skew", "0", "", leftRange[i], rightRange[i]); // "0" Do not plot the result of the fit // Fit results TF1* fit = h2->GetFunction("skew"); constants[i] = fit->GetParameter(0); error_constants[i] = fit->GetParError(0); means[i] = fit->GetParameter(1); error_means[i] = fit->GetParError(1); gaussianArea[i] = fit->Integral(means[i] - leftRange[i], means[i] + rightRange[i]); } ///////////////////////////////////////////////////////////////////////////// for(int i = 0; i < peaksFound; i++) { // skewed gaussians and default parameters skew[i] = new TF1("skew", "[0] * exp(-0.5*((x-[1])/([2]+(x<[1])*[3]*(x-[1])))^2)"); skew[i]->SetParNames("Height", "Centroid", "Sigma", "Skewness"); skew[i]->SetParLimits(0, 1, 7500); skew[i]->SetParameter(0, 3500); skew[i]->SetParLimits(1, 10, 160); skew[i]->SetParameter(1, 100); skew[i]->SetParLimits(2, 1, 3); skew[i]->FixParameter(2, 2.03945); skew[i]->SetParLimits(3, 0.1, 0.3); skew[i]->FixParameter(3, 0.199132); // Fit Each peak appropratelty skew[i]->SetParameter(1, peakPosition[i]); // if(i == 0) // Actual Fix // { // skew[i]->FixParameter(1, 17.6875); // skew[i]->FixParameter(0, 436.632); // } if(i == 0) // Actual Fix 419.5 { skew[i]->FixParameter(1, 17.6875); skew[i]->FixParameter(0, 405); //skew[i]->FixParameter(2, 3); skew[i]->FixParameter(3, 0.45); } if(i == 1) { skew[i]->FixParameter(1, 24.603); // Actual Fix skew[i]->FixParameter(0, 1072.16); } if(i == 2) { skew[i]->FixParameter(1, 30.0); skew[i]->FixParameter(0, 130); } if(i == 3) { skew[i]->FixParameter(1, 34.9); skew[i]->FixParameter(0, 120); } if(i == 4) { skew[i]->FixParameter(1, 38.8); // 36 skew[i]->FixParameter(0, 180); //140 } if(i == 5) { skew[i]->FixParameter(1, 41.8); skew[i]->FixParameter(0, 520); //480 } if(i == 6) // Actual Fix { skew[i]->FixParameter(1, 45.2752); skew[i]->FixParameter(0, 7050); } if(i == 7) { skew[i]->FixParameter(1, 49.5); skew[i]->FixParameter(0, 2050); } if(i == 8) { skew[i]->FixParameter(1, 55); skew[i]->FixParameter(0, 260); } if(i == 9) { skew[i]->FixParameter(1, 59.2); skew[i]->FixParameter(0, 135); } if(i == 10) { skew[i]->FixParameter(1, 63); skew[i]->FixParameter(0, 120); } if(i == 11) { skew[i]->FixParameter(1, 66.7); skew[i]->FixParameter(0, 240); } if(i == 12) { skew[i]->FixParameter(1, 69.5); skew[i]->FixParameter(0, 220); } if(i == 13) { skew[i]->FixParameter(1, 73.7); skew[i]->FixParameter(0, 300); } if(i == 14) { skew[i]->FixParameter(1, 77.5); skew[i]->FixParameter(0, 330); } if(i == 15) { skew[i]->FixParameter(1, 80.5); skew[i]->FixParameter(0, 270); } if(i == 16) { skew[i]->FixParameter(1, 83.7); skew[i]->FixParameter(0, 215); } if(i == 17) { skew[i]->FixParameter(1, 87.4); skew[i]->FixParameter(0, 250); } if(i == 18) { skew[i]->FixParameter(1, 91); skew[i]->FixParameter(0, 130); } if(i == 19) { skew[i]->FixParameter(1, 93); skew[i]->FixParameter(0, 140); } if(i == 20) { skew[i]->FixParameter(1, 96.2); skew[i]->FixParameter(0, 130); } if(i == 21) { skew[i]->FixParameter(1, 99.4); skew[i]->FixParameter(0, 90); } if(i == 22) { skew[i]->FixParameter(1, 103.993); skew[i]->FixParameter(0, 320.802); } if(i == 23) { skew[i]->FixParameter(1, 109); skew[i]->FixParameter(0, 80); } if(i == 24) { skew[i]->FixParameter(1, 113.7); skew[i]->FixParameter(0, 43); } if(i == 25) { skew[i]->FixParameter(1, 118.61); skew[i]->FixParameter(0, 64); } if(i == 26) { skew[i]->FixParameter(1, 123.5); skew[i]->FixParameter(0, 20); } if(i == 27) { skew[i]->FixParameter(1, 128.221); skew[i]->FixParameter(0, 30); } if(i == 28) { skew[i]->FixParameter(1, 132); skew[i]->FixParameter(0, 15); } if(i == 29) { skew[i]->FixParameter(1, 138); skew[i]->FixParameter(0, 5); } if(i == 30) { skew[i]->FixParameter(1, 147); skew[i]->FixParameter(0, 4); } if(i == 31) { skew[i]->FixParameter(1, 152.5); skew[i]->FixParameter(0, 7.5); } skew[i]->SetLineColor(kGreen+1); h2->Fit("skew", "Q+", "", leftRange[i], rightRange[i]); } // gui settings h2->SetTitle("22Si Proton Spectrum"); h2->SetLineColor(kBlue); h2->SetAxisRange(10.5, 160.5); // whole spectrum //h2->SetAxisRange(10.5, 35.5); // low E //h2->SetAxisRange(34.5, 58.5); // doublet //h2->SetAxisRange(60.5, 110.5); // mid E //h2->SetAxisRange(107.5, 160.5); // high E gPad->SetLogy(0); h2->Draw(); // Draw global fit // A single fit for the whole plot TGraph* g1 = new TGraph(); for(double i = 0; i < 256; i+= 0.25) // increase i by 0.25 { double sum = 0.0; for(int j = 0; j < peaksFound; j++) { sum += skew[j]->Eval(i); } g1->SetPoint(i*4, i, sum); } g1->SetLineColor(kRed); g1->SetLineWidth(2); g1->Draw("same"); cout << endl; cout << "********************************************************************" << endl; cout << endl; // display results for quick reading of the fit parameters double energies[peaksFound]; cout << "Peak# Height Centroid Area Energy Peak?" << endl; //drawDash(50); for(int i = 0; i < peaksFound; i ++) { energies[i] = means[i] * 61.2144 - 792.425; cout << i << '\t' << constants[i] << '\t' << means[i] << '\t' << gaussianArea[i] << '\t' << energies[i]; if( i == 2 || i == 6 || i == 7 || i == 11 || i == 12 || i == 22 || i == 25 || i == 27 || i == 31 ) { cout << '\t' << "<---- ***"; } cout << endl; //drawDash(50); } cout << endl; cout << "********************************************************************" << endl; cout << endl; /* Fit Statistics */ bool isSignificant[peaksFound]; double ggchi2 = 0; // for Red global fit for(int i = 0; i < 256; i++) { if(g1->Eval(i) < 1) { ggchi2 += pow(h2->GetBinContent(i) - g1->Eval(i), 2); } else { ggchi2 += pow(h2->GetBinContent(i) - g1->Eval(i), 2) / g1->Eval(i); // this is the chi2 formula } } //drawDash(40); cout << "| Global Histogram Area: " << h2->Integral() << " |" << endl; cout << "| Global Fit Function Area: " << g1->Integral() << " |" << endl; cout << "| Global Fit Chi2: " << ggchi2 << " |" << endl; //drawDash(40); cout << endl; // this is an extra endl on top of the default one in drawDash cout << "Peak# Global Local~ Local Global Local Peak (gr * lr / area > 0.6)" << endl; cout << " Chi2 Chi2 Chi2 Ratio Ratio Cutoff" << endl; //drawDash(70); for(int i = 0; i < peaksFound; i++) { double gchi2 = 0; double llchi2 = 0; double lchi2 = 0; double peakTest; // subtract current peak for(int j = leftRange[i]; j < rightRange[i] + 1; j++) { g1->SetPoint(j, j, g1->Eval(j) - skew[i]->Eval(j)); } // calculate local chi2, without current peak for(int k = means[i] - 10; k < means[i] + 10; k++) { if(g1->Eval(k) < 1) { lchi2 += pow(h2->GetBinContent(k) - g1->Eval(k), 2); } else { lchi2 += pow(h2->GetBinContent(k) - g1->Eval(k), 2) / g1->Eval(k); } } // calculate global chi2, without current peak for(int l = 0; l < 256; l++) { if(g1->Eval(l) < 1) { gchi2 += pow(h2->GetBinContent(l) - g1->Eval(l), 2); } else { gchi2 += pow(h2->GetBinContent(l) - g1->Eval(l), 2) / g1->Eval(l); } } // re-add current peak for(int m = leftRange[i]; m < rightRange[i] + 1; m++) { g1->SetPoint(m, m, g1->Eval(m) + skew[i]->Eval(m)); } // calculate local chi2, with current peak for(int n = means[i] - 10; n < means[i] + 10; n++) { if(g1->Eval(n) < 1) { llchi2 += pow(h2->GetBinContent(n) - g1->Eval(n), 2); } else { llchi2 += pow(h2->GetBinContent(n) - g1->Eval(n), 2) / g1->Eval(n); } } // output peakTest = (gchi2 / ggchi2) * pow(lchi2 / llchi2, 1) / gaussianArea[i]; cout << i << '\t' << gchi2 << '\t' << llchi2 << '\t' << lchi2 << '\t' << gchi2 / ggchi2 << '\t' << lchi2 / llchi2 << '\t' << peakTest; if(peakTest > 0.6 && i != 0) { isSignificant[i] = true; cout << " <----"; } else { isSignificant[i] = false; } cout << endl; //drawDash(70); } // Draw global fit for(double i = 0; i < 256; i+= 0.25) { double sum = 0.0; for(int j = 0; j < peaksFound; j++) { sum += skew[j]->Eval(i); } g1->SetPoint(i*4, i, sum); } g1->Draw("same"); // creates energy labels for each peak for(int i = 0; i < peaksFound; i++) { string number1 = std::to_string(i); const char* number2 = number1.c_str(); if (i < 22 && i > 4) { TText* t1 = new TText(means[i] - 0.5, constants[i] + 150, number2); t1->SetTextSize(0.015); if(isSignificant[i]) { t1->SetTextColor(6); } t1->Draw("+"); } else if (i > 22) { TText* t3 = new TText(means[i] - 0.5, constants[i] + 30, number2); t3->SetTextSize(0.015); if(isSignificant[i]) { t3->SetTextColor(6); } t3->Draw("+"); } else { TText* t2 = new TText(means[i] - 0.5, constants[i] + 60, number2); t2->SetTextSize(0.015); if(isSignificant[i]) { t2->SetTextColor(6); } t2->Draw("+"); } } // Convert X-axis to Energy convertAxisAndDrawX (h2); // make legend for the different colored lines TLegend* legend = new TLegend(0.7, 0.7, 0.89, 0.89, NULL, "brNDC"); TLegendEntry* entry2 = legend->AddEntry("h1", "Raw Data"); entry2->SetMarkerStyle(0); entry2->SetLineColor(4); entry2->SetMarkerColor(4); TLegendEntry* entry = legend->AddEntry("g1","Global Fitting Function"); entry->SetMarkerStyle(0); entry->SetLineColor(2); entry->SetMarkerColor(2); TLegendEntry* entry1 = legend->AddEntry("skew","Skewed Gaussians"); entry1->SetMarkerStyle(0); entry1->SetLineColor(kGreen+1); entry1->SetMarkerColor(kGreen+1); legend->Draw(); cout << endl; cout << "********************************************************************" << endl; cout << endl; cout << "The Proton Peaks of Si22 Decay:" << endl << endl; double total_BR = 0; double resolved = 0; for(int i = 0; i < peaksFound; i++) { if(isSignificant[i]) { total_BR += gaussianArea[i] / (g1->Integral() - gaussianArea[0]) * 100; cout << energies[i] << " KeV " << "at " << gaussianArea[i] / (g1->Integral() - gaussianArea[0]) * 100 << " % BR" << endl; cout << "---------------------------------" << endl; resolved += gaussianArea[i] / (g1->Integral() - gaussianArea[0]) * 100; } } cout << "Unresolved peaks at " << 100 - resolved << " % BR" << endl; cout << "---------------------------------" << endl; //cout << "Total of all BR: " << total_BR + (100 - resolved) << endl; //cout << "---------------------------------" << endl; for(int i = 0; i < peaksFound; i++) { cout << (gaussianArea[i] / (h2->Integral() - gaussianArea[0]) * 100) - (gaussianArea[i] / (g1->Integral() - gaussianArea[0]) * 100) << endl; } h2->FitPanel(); cout << endl; cout << "********************************************************************" << endl; cout << endl; cout << "Coincidence Check" << endl << "-------------------" << endl; cout << "206KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 146 && (energies[j] - energies[i]) < 266) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "332KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 272 && (energies[j]) - (energies[i]) < 392) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "596KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 536 && (energies[j] - energies[i]) < 656) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "717KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 657 && (energies[j] - energies[i]) < 777) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "880KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 820 && (energies[j] - energies[i]) < 940) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "1633KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 1573 && (energies[j] - energies[i]) < 1693) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "1088KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 1028 && (energies[j] - energies[i]) < 1148) { cout << "(" << i << "," << j << ")"; } } } cout << endl; //drawDash(150); cout << "3243KeV Difference: "; for(int i = 1; i < peaksFound; i++) { for(int j = i; j < peaksFound; j++) { if ((energies[j] - energies[i]) > 3183 && (energies[j] - energies[i]) < 3303) { cout << "(" << i << "," << j << ")"; } } } cout << endl; cout << endl; cout << "********************************************************************" << endl; cout << endl; /* Gated Spectrum */ ifstream inData_g; TString inputFileName_g = Form("%sRawData/Si22/corr.dmax.2!decay_g206.dat", dir.Data()); inData_g.open(inputFileName_g); if(!(inData_g.is_open())) cout << "File Opening Failed" << endl; TCanvas* c1_g = new TCanvas("c1_g", "22Si Protons", 1300, 1000); c1_g->ToggleEventStatus(); //c1->Divide(1,2); /* Load Histogram Data */ //c1->cd(1); TH1I* h1_g = new TH1I("h1_g", "22Si - Gate on gamma; Channel Number; Number of Counts", 256, 0.5, 256.5); int channelNumber_g, nextChannel_g = 100000; double countNumber_g; // inputs data to histogram while(inData_g) { inData_g >> channelNumber_g >> countNumber_g; // corrects for last channel duplication if(channelNumber_g == nextChannel_g) break; nextChannel_g = channelNumber_g; h1_g->SetBinContent(channelNumber_g, countNumber_g); } inData_g.close(); h1_g->SetAxisRange(0.5, 256.5); h1_g->Draw(); // proper scale to match ungated spectrum int min_g = h1_g->GetMinimum(); int max_g = h1_g->GetMaximum(); double scale = h1->GetMaximum() / (max_g - min_g); min_g = 0; double scale2 = h1->GetMaximum() / max_g; // draw on ungated spectrum TH1I* h3_g = (TH1I*)h1_g->Clone("h3_g"); for(int i = 0; i < 256; i++) { h3_g->SetBinContent(i, (h1_g->GetBinContent(i) - min_g) * scale2); //h3_g->SetBinContent(i, h1_g->GetBinContent(i) * 5.10989011); } h3_g->SetLineColor(kGray+2); c1->cd(); //h3_g->Draw("SAME"); c1_g->cd(); // creates energy labels for each peak for(int i = 0; i < peaksFound_g; i++) { string number1_g = std::to_string(i); const char* number2_g = number1_g.c_str(); TText* t1_g = new TText(means[i] - 0.5, h1_g->GetBinContent(means[i]) + 2, number2_g); t1_g->SetTextSize(0.015); if(isSignificant[i]) t1_g->SetTextColor(6); t1_g->Draw("+"); } /* Decay Scheme Graph */ TCanvas* c2 = new TCanvas("c2", "22Si Decay Graph", 1300, 1000); double branchingRatio[9] = {7.38, 62.66 + 7.38, 3.17 + 62.66 + 7.38, 4.34 + 3.17 + 62.66 + 7.38, 1.72 + 4.34 + 3.17 + 62.66 + 7.38, 0.44 + 1.72 + 4.34 + 3.17 + 62.66 + 7.38, 0.61 + 0.44 + 1.72 + 4.34 + 3.17 + 62.66 + 7.38, 0.05 + 0.61 + 0.44 + 1.72 + 4.34 + 3.17 + 62.66 + 7.38, 2.51 + 0.05 + 0.61 + 0.44 + 1.72 + 4.34 + 3.17 + 62.66 + 7.38}; double energyLevels[9] = {740, 2250, 3480, 3970, 4580, 6490, 7300, 8560, 9410}; double x[100000], y[100000]; for(int i = 0; i < 100000; i++) { double j = i; x[i] = j / 1000; //cout << x[i] << endl; if(x[i] < branchingRatio[0]) y[i] = energyLevels[0]; else if(x[i] < branchingRatio[1]) y[i] = energyLevels[1]; else if(x[i] < branchingRatio[2]) y[i] = energyLevels[2]; else if(x[i] < branchingRatio[3]) y[i] = energyLevels[3]; else if(x[i] < branchingRatio[4]) y[i] = energyLevels[4]; else if(x[i] < branchingRatio[5]) y[i] = energyLevels[5]; else if(x[i] < branchingRatio[6]) y[i] = energyLevels[6]; else if(x[i] < branchingRatio[7]) y[i] = energyLevels[7]; else if(x[i] < branchingRatio[8]) y[i] = energyLevels[8]; else if(x[i] < 100) y[i] = 0; } TGraph* g2 = new TGraph(100000, x, y); g2->SetTitle("Decay Scheme Graph; Branching Ratio (%); Energy Level (KeV)"); //g2->SetMarkerStyle(20); g2->SetMarkerColor(kBlue); g2->SetLineColor(kBlue); g2->Draw(); /* Ground Transitions */ if(h1_g->GetMaximum() > 1000) { TCanvas* c3 = new TCanvas("c3", "22Si Ground Transitions", 1300, 1000); TH1I* h4_g = (TH1I*)h1_g->Clone("h4_g"); TH1I* h3 = (TH1I*)h1->Clone("h3"); h3->Add(h4_g, -1/(0.140331695575063)); h3->Draw(); for(int i = 0; i < peaksFound_g; i++) { string number1_g = std::to_string(i); const char* number2_g = number1_g.c_str(); TText* t1_g = new TText(means[i] - 0.5, h3->GetBinContent(means[i]) + 50, number2_g); t1_g->SetTextSize(0.015); if(isSignificant[i]) t1_g->SetTextColor(6); t1_g->Draw("+"); } } } // end of option 5 // *********************************************************************************** // function Implementations // *********************************************************************************** void drawDash (int counter) { for (int a = 0; a < counter; a++) { cout << "-"; } cout << endl; } void convertAxisAndDrawX (TH1I* hist_in) { // Convert X-axis to Energy hist_in->GetXaxis()->SetBinLabel((0 + 792.425) / 61.2144, "0"); hist_in->GetXaxis()->SetBinLabel((1000 + 792.425) / 61.2144, "1000"); hist_in->GetXaxis()->SetBinLabel((2000 + 792.425) / 61.2144, "2000"); hist_in->GetXaxis()->SetBinLabel((3000 + 792.425) / 61.2144, "3000"); hist_in->GetXaxis()->SetBinLabel((4000 + 792.425) / 61.2144, "4000"); hist_in->GetXaxis()->SetBinLabel((5000 + 792.425) / 61.2144, "5000"); hist_in->GetXaxis()->SetBinLabel((6000 + 792.425) / 61.2144, "6000"); hist_in->GetXaxis()->SetBinLabel((7000 + 792.425) / 61.2144, "7000"); hist_in->GetXaxis()->SetBinLabel((8000 + 792.425) / 61.2144, "8000"); hist_in->GetXaxis()->SetBinLabel((9000 + 792.425) / 61.2144, "9000"); hist_in->GetXaxis()->SetBinLabel((10000 + 792.425) / 61.2144, "10000"); hist_in->GetXaxis()->SetBinLabel((500 + 792.425) / 61.2144, "500"); hist_in->GetXaxis()->SetBinLabel((1500 + 792.425) / 61.2144, "1500"); hist_in->GetXaxis()->SetBinLabel((2500 + 792.425) / 61.2144, "2500"); hist_in->GetXaxis()->SetBinLabel((3500 + 792.425) / 61.2144, "3500"); hist_in->GetXaxis()->SetBinLabel((4500 + 792.425) / 61.2144, "4500"); hist_in->GetXaxis()->SetBinLabel((5500 + 792.425) / 61.2144, "5500"); hist_in->GetXaxis()->SetBinLabel((6500 + 792.425) / 61.2144, "6500"); hist_in->GetXaxis()->SetBinLabel((7500 + 792.425) / 61.2144, "7500"); hist_in->GetXaxis()->SetBinLabel((8500 + 792.425) / 61.2144, "8500"); hist_in->GetXaxis()->SetBinLabel((9500 + 792.425) / 61.2144, "9500"); hist_in->GetXaxis()->LabelsOption("h"); // draw labels horizontal hist_in->GetXaxis()->SetLabelSize(0.03); hist_in->GetXaxis()->SetTickLength(0.00); hist_in->GetXaxis()->SetTitle("Energy (KeV)"); }