Geant4 Cross Reference |
1 //******************************************** 1 2 // Spectrum_proton.C 3 // Root command file 4 // Type: root Spectrum_proton.C 5 // 6 // It visualizes the spectrum of protons and p 7 // simulation result ProtonAtExit.dat 8 // 9 // More information is available in UserGuide 10 // Created by Z.LI LP2i Bordeaux 2022 11 //******************************************** 12 13 #include <math.h> 14 #include <stdint.h> 15 #include <stdio.h> 16 #include <string.h> 17 18 #include <vector> 19 // using namespace std; 20 21 // Define a structure to read and write each e 22 struct RunInfo 23 { 24 // uint_16t 25 uint8_t projectionIndex; // 1 byte 26 uint16_t sliceIndex; // 27 uint16_t pixelIndex; 28 uint32_t nbParticle; // 4 bytes int 29 }; 30 31 struct ParticleInfo 32 { 33 float energy_keV; 34 float mx; 35 float my; 36 float mz; 37 }; 38 39 void Plot(vector<double>& energies, int bin, d 40 { 41 auto mycanvas = new TCanvas("canvas", "canva 42 gPad->SetLeftMargin(0.15); 43 44 // unit is in keV 45 auto hist = new TH1D("hist (keV)", "Spectrum 46 47 for (int i = 0; i < energies.size(); ++i) { 48 hist->Fill(energies[i]); 49 } 50 51 hist->Draw(); 52 hist->GetXaxis()->SetTitle("Energy (keV)"); 53 hist->GetYaxis()->SetTitle("Counts"); 54 hist->GetXaxis()->CenterTitle(); 55 hist->GetYaxis()->CenterTitle(); 56 57 mycanvas->Print("spectrum_proton.png"); 58 } 59 60 void Spectrum_proton() 61 { 62 FILE* input = fopen("../build/ProtonAtExit.d 63 if (input == NULL) { 64 printf("error for opening the input file\n 65 return; 66 } 67 68 //****************************************** 69 //**************************Selection parame 70 //****************************************** 71 72 const int nbProjection = 10; 73 const int nbSlice = 128; 74 const int nbPixel = 20; 75 76 int projection_index_begin = 0; // starter 77 int projection_index_end = 0; // end of the 78 79 int slice_index_begin = 64; // starter of t 80 int slice_index_end = 64; // end of the sli 81 82 //********************Parameters for spectru 83 int bin = 100; 84 double eMin = 0; // keV 85 double eMax = 0; // keV 86 87 //****************************************** 88 //**************************Selection parame 89 //****************************************** 90 91 RunInfo runInfo; 92 vector<double> energies; 93 int runID = -1; // index of simulations, na 94 // while(!feof(input)) //if not the end, rea 95 while (fread(&runInfo, sizeof(RunInfo), 1, i 96 runID++; 97 int nbParticle = runInfo.nbParticle; 98 99 // ***********the following codes are used 100 // if************************************* 101 // ***********the index of projection, sli 102 // configured in the simulation 103 runInfo.projectionIndex = runID / (nbSlice 104 int remain = runID % (nbSlice * nbPixel); 105 runInfo.sliceIndex = remain / nbPixel; 106 runInfo.pixelIndex = remain % nbPixel; 107 //**************************************** 108 109 if (!nbParticle) continue; 110 std::vector<ParticleInfo> proton(nbParticl 111 fread(&proton[0], sizeof(ParticleInfo), nb 112 113 if (runInfo.projectionIndex >= projection_ 114 && runInfo.projectionIndex <= projecti 115 { 116 if (runInfo.sliceIndex >= slice_index_be 117 for (int i = 0; i < nbParticle; ++i) { 118 energies.push_back(proton[i].energy_ 119 if (proton[i].energy_keV > eMax) eMa 120 } 121 } 122 } 123 else 124 break; 125 } 126 127 fclose(input); 128 Plot(energies, bin, eMin, eMax + 10); 129 } 130