xray_TESdetector/analysis/analysis.C

0001 #include <TROOT.h>
0002 #include <TTree.h>
0003 #include <TCanvas.h>
0004 #include <TApplication.h>
0005 #include <TString.h>
0006 #include <TBranch.h>
0007 #include <TStyle.h>
0008 #include <TPie.h>
0009 #include <TFile.h>
0010 #include <TLegend.h>
0011 #include <TH1D.h>
0012 #include <TH2D.h>
0013 #include <iostream>
0014 #include <fstream>
0015 #include <string>
0016 #include <sstream>
0017 #include <vector>
0018 #include <stdio.h>
0019 #include <stdlib.h>
0020 #include <iomanip>
0021 #include <algorithm>
0022
0023 using namespace std;
0024 TApplication myapp("app",NULL,NULL);
0025
0026 //////////////////////////////////////////////////////
0027
0028 int main(int argc, char *argv[]){
0029   string ifilename = "";
0030   int no_events = 10000000; // number of primaries
0031
0032   if(argc == 1)
0033   {
0034       cout << "No input file selected. Please use: ./read_tree_spectrum [input_file.root]" << endl;
0035       return 0;
0036   }
0037   else
0038   {
0039       ifilename = argv[1];
0040       if (argc == 3) no_events = atoi(argv[2]);
0041   }
0042
0043   TString ifilenameforroot = ifilename;
0044   TFile *f = new TFile(ifilenameforroot, "READ");
0045
0046   int no_pixel = 317;
0047
0048   // Histogram definition
0049   //  Detector deposits
0050   TH1D *pixel_dep = new TH1D ("pixel", "Avg. energy deposit per pixel [MeV]", no_pixel, 0, no_pixel);  //(nbinx, xdown, xup, nbiny, ylow, yup) for drawing the hexagon
0051   TH1D *part_per_pixel = new TH1D ("Part per pixel", "No. particles per pixel", no_pixel, 0, no_pixel);
0052
0053   // Spectra for energies (initial, incident and deposited)
0054   int bin_number = 500;
0055   int max_val = 10000;
0056   int min_val = 0;
0057   double bin_width = (max_val-min_val)/(double)bin_number;
0058   TH1D *initial_energies = new TH1D ("energies4", "GCR Protons impacting spectra on the detector", bin_number, min_val, max_val);
0059   TH1D *incident_energies = new TH1D ("energies2", " ", bin_number, min_val, max_val);
0060   TH1D *dep_energies = new TH1D ("energies3", " ", bin_number, min_val, max_val);
0061
0062   // Spectra for energies (initial, incident and deposited) - for primaries (protons)
0063   TH1D *initial_energies_P = new TH1D ("energies4_P", " ", bin_number, min_val, max_val);
0064   TH1D *incident_energies_P = new TH1D ("energies2_P", " ", bin_number, min_val, max_val);
0065   TH1D *dep_energies_P = new TH1D ("energies3_P", "Primary GCR Protons impacting spectra on the detector", bin_number, min_val, max_val);
0066
0067   // Incident particles
0068   TH1D *inc_particles = new TH1D ("energies6", "GCR Protons particle fluxes on the detector composition", 5, 0, 5);
0069
0070   // Variables to store the tuples' values
0071   int eventID = 0;
0072   char vol_name[500];
0073   int trackID = 0;
0074   double x = 0;
0075   double y = 0;
0076   double z = 0;
0077   double theta = 0;
0078   double phi = 0;
0079   int parentID = 0;
0080   int pixel_number = 0;
0081   double step_energy_dep = 0;
0082   int step_number = 0;
0083   double init_kinetic_energy = 0;
0084   double kinetic_energy = 0;
0085   char particle_name[500];
0086   char pre_step_name[500];
0087   char post_step_name[500];
0088
0089   // Define tuple elements
0090   TTree *mytree = (TTree*)f->Get("TES_Tuple");
0091   Long64_t nentries = mytree->GetEntries();
0092   mytree->SetBranchAddress("eventID", &eventID);
0093   mytree->SetBranchAddress("vol_name", &vol_name);
0094   mytree->SetBranchAddress("trackID", &trackID);
0095   mytree->SetBranchAddress("x", &x);
0096   mytree->SetBranchAddress("y", &y);
0097   mytree->SetBranchAddress("z", &z);
0098   mytree->SetBranchAddress("theta", &theta);
0099   mytree->SetBranchAddress("phi", &phi);
0100   mytree->SetBranchAddress("parentID", &parentID);
0101   mytree->SetBranchAddress("pixel_number", &pixel_number);
0102   mytree->SetBranchAddress("step_energy_dep", &step_energy_dep);
0103   mytree->SetBranchAddress("step_number", &step_number);
0104   mytree->SetBranchAddress("init_kinetic_energy", &init_kinetic_energy);
0105   mytree->SetBranchAddress("kinetic_energy", &kinetic_energy);
0106   mytree->SetBranchAddress("particle_name", &particle_name);
0107   mytree->SetBranchAddress("pre_step_name", &pre_step_name);
0108   mytree->SetBranchAddress("post_step_name", &post_step_name);
0109
0110   cout << "Reading " << nentries << " from the TTree." << endl;
0111   string vol_name_old = "";
0112
0113   int eventID_old = 0;
0114   int trackID_old = 0;
0115   double total_step_dep = 0;
0116   double total_step_dep_P = 0;
0117   double kine = 0;        // kinetic energy
0118   double kine_P = 0;      // kinetic energy for protons
0119   double in_kine = 0;     // initial kinetic energy
0120   double in_kine_P = 0;   // initial kinetic energy for protons
0121   bool entered = false;
0122   bool entered_P = false;
0123   char arr[50];
0124   string str = "";
0125   int pixel_number_old = 0;
0126
0127   for (Long64_t i=0; i<nentries; i++)
0128   {
0129     mytree->GetEntry(i);
0130
0131     if (entered || entered_P)
0132     {
0133       if (eventID != eventID_old || (eventID == eventID_old && trackID != trackID_old))
0134       {
0135         if (entered && total_step_dep > 0)
0136         {
0137           initial_energies->Fill(in_kine);
0138           incident_energies->Fill(kine);
0139           dep_energies->Fill(total_step_dep);
0140           inc_particles->Fill(arr, 1);
0141           total_step_dep = 0;
0142         }
0143         if (entered_P && total_step_dep_P > 0)
0144         {
0145           initial_energies_P->Fill(in_kine_P);
0146           incident_energies_P->Fill(kine_P);
0147           dep_energies_P->Fill(total_step_dep_P);
0148           total_step_dep_P = 0;
0149         }
0150         entered = false;
0151         entered_P = false;
0152       }
0153     }
0154     str = (string)particle_name;
0155
0156     if((string)vol_name == "Bipxl")
0157     {
0158       if (vol_name_old != "Bipxl")
0159       {
0160         if (!entered)
0161         {
0162           entered = true;
0163           kine = kinetic_energy;
0164           in_kine = init_kinetic_energy;
0165
0166           if (str != "proton" && str != "gamma" && str != "e+" && str != "e-")
0167           {
0168             str = "Other particles";
0169           }
0170           strcpy(arr, str.c_str());
0171
0172           // Same but onyl for primaries
0173           if (!entered_P)
0174           {
0175             if (str == "proton" && parentID == 0)
0176             {
0177               entered_P = true;
0178               kine_P = kinetic_energy;
0179               in_kine_P = init_kinetic_energy;
0180             }
0181           }
0182         }
0183       }
0184
0185       // Sum the total energy deposit from all particles
0186       total_step_dep += step_energy_dep;
0187
0188       // Sum the total energy deposit from the primaries
0189       if (str == "proton" && parentID == 0) total_step_dep_P += step_energy_dep;
0190       // hist_det_count->Fill(x, y, step_energy_dep);
0191
0192       if (pixel_number != pixel_number_old)
0193       {
0194         part_per_pixel->Fill(-pixel_number);
0195       }
0196
0197       pixel_dep->Fill(-pixel_number, step_energy_dep);
0198     }
0199     eventID_old = eventID;
0200     trackID_old = trackID;
0201     vol_name_old = (string)vol_name;
0202     pixel_number_old = pixel_number;
0203   }
0204
0205   cout << "Starts plotting..." << endl;
0206
0207   // Draw histograms
0208   TCanvas *c1 = new TCanvas ("c1", "Energy dep", 0, 0, 1000, 900);
0209   double radius = 27.0;
0210   double pi = 3.1415;
0211   double T = no_events/(0.407*4*pi*pi*radius*radius);     //equivalent time
0212   double S = 2.1; //cm2
0213   dep_energies->Scale(1.0/(T*S*bin_width));               //scale based on the equivalent time and detector surface
0214   initial_energies->Scale(1.0/(T*S*bin_width));
0215   incident_energies->Scale(1.0/(T*S*bin_width));
0216   gStyle->SetOptStat(0);
0217   c1->SetLogx();
0218   c1->SetLogy();
0219   dep_energies->SetLineColor(kRed);
0220   initial_energies->Draw("hist");
0221   initial_energies->SetLineColor(kBlack);
0222   incident_energies->Draw("histsame");
0223   dep_energies->Draw("histsame");
0224   incident_energies->SetLineColor(kBlue);
0225   initial_energies->GetYaxis()->SetRangeUser(1.0e-6, 1.0);
0226   initial_energies->GetXaxis()->SetTitle("Energy [MeV]");
0227   initial_energies->GetYaxis()->SetTitle("Counts/cm2/s/MeV");
0228   auto legend1 = new TLegend(0.75,0.8,0.9,0.9);
0229   legend1->AddEntry(dep_energies,"Edep");
0230   legend1->AddEntry(incident_energies, "Einc");
0231   legend1->AddEntry(initial_energies, "Ei");
0232   legend1->Draw("same");
0233
0234   TCanvas *c2 = new TCanvas ("c2", "Proton energy dep", 50, 0, 1000, 900);
0235   dep_energies_P->Scale(1.0/(T*S*bin_width));
0236   initial_energies_P->Scale(1.0/(T*S*bin_width));
0237   incident_energies_P->Scale(1.0/(T*S*bin_width));
0238   gStyle->SetOptStat(0);
0239   c2->SetLogx();
0240   c2->SetLogy();
0241   dep_energies_P->Draw("hist");
0242   dep_energies_P->SetLineColor(kRed);
0243   initial_energies_P->Draw("histsame");
0244   initial_energies_P->SetLineColor(kBlack);
0245   incident_energies_P->Draw("histsame");
0246   incident_energies_P->SetLineColor(kBlue);
0247   dep_energies_P->GetYaxis()->SetRangeUser(1.0e-6, 1.0);
0248   dep_energies_P->GetXaxis()->SetTitle("Energy [MeV]");
0249   dep_energies_P->GetYaxis()->SetTitle("Counts/cm2/s/MeV");
0250   auto legend2 = new TLegend(0.75,0.8,0.9,0.9);
0251   legend2->AddEntry(dep_energies,"Edep");
0252   legend2->AddEntry(incident_energies, "Einc");
0253   legend2->AddEntry(initial_energies, "Ei");
0254   legend2->Draw("same");
0255
0256   TCanvas *cpie = new TCanvas("Particless distribution", "Particles distribution", 100, 0, 1000, 1000);
0257   string label1 = inc_particles->GetXaxis()->GetBinLabel(1);
0258   string label2 = inc_particles->GetXaxis()->GetBinLabel(2);
0259   string label3 = inc_particles->GetXaxis()->GetBinLabel(3);
0260   string label4 = inc_particles->GetXaxis()->GetBinLabel(4);
0261   string label5 = inc_particles->GetXaxis()->GetBinLabel(5);
0262   Float_t val1 = inc_particles->GetBinContent(1);
0263   Float_t val2 = inc_particles->GetBinContent(2);
0264   Float_t val3 = inc_particles->GetBinContent(3);
0265   Float_t val4 = inc_particles->GetBinContent(4);
0266   Float_t val5 = inc_particles->GetBinContent(5);
0267   Float_t vals[] = {val1,val2,val3,val4,val5};
0268   Int_t colors[] = {1,2,3,4,5};
0269   Int_t nvals = sizeof(vals)/sizeof(vals[0]);
0270   TPie *pie1 = new TPie("pie1", "Particles distribution",nvals,vals,colors);
0271   pie1->SetLabelsOffset(.01);
0272   pie1->SetRadius(.2);
0273   pie1->SetEntryLabel(0, label1.c_str());
0274   pie1->SetEntryLabel(1, label2.c_str());
0275   pie1->SetEntryLabel(2, label3.c_str());
0276   pie1->SetEntryLabel(3, label4.c_str());
0277   pie1->SetEntryLabel(4, label5.c_str());
0278   pie1->SetLabelFormat("%txt (%perc)");
0279   pie1->Draw();
0280
0281   TCanvas *c4_num = new TCanvas ("c4_num", "Particle count per pixel", 150, 0, 1000, 900);
0282   part_per_pixel->GetXaxis()->SetTitle("Pixel number");
0283   part_per_pixel->GetYaxis()->SetTitle("Counts");
0284   part_per_pixel->SetFillColor(kBlack);
0285   part_per_pixel->SetMarkerColor(kBlack);
0286   part_per_pixel->Draw();
0287
0288   // Draw the histogram of the single pixel deposition
0289   TCanvas *c5 = new TCanvas ("c5", "Avg. energy deposit per pixel", 250, 0, 1000, 900);
0290   pixel_dep->Divide(part_per_pixel);
0291   pixel_dep->GetXaxis()->SetTitle("Pixel number");
0292   pixel_dep->GetYaxis()->SetTitle("Avg. energy deposit [MeV]");
0293   pixel_dep->SetFillColor(kBlack);
0294   pixel_dep->SetMarkerColor(kBlack);
0295   pixel_dep->Draw("");
0296
0297   myapp.Run(true);
0298   f->Close();
0299   return 0;
0300 }