|
| 1 | +#include <utility> |
| 2 | +#include <vector> |
| 3 | +#include <TString.h> |
| 4 | +#include "TCut.h" |
| 5 | +#include <list> |
| 6 | +#include <dirent.h> |
| 7 | +#include <unistd.h> |
| 8 | +#include <fstream> |
| 9 | +#include <iostream> |
| 10 | +#include <sstream> |
| 11 | +#include <string> |
| 12 | +#include <algorithm> |
| 13 | +#include <TSystem.h> |
| 14 | +#include <TChain.h> |
| 15 | +#include <TFile.h> |
| 16 | +#include <TCanvas.h> |
| 17 | +#include <TROOT.h> |
| 18 | +#include <TH1.h> |
| 19 | +#include <TH2.h> |
| 20 | +void plotDetPlaneCuts(Double_t r_outer_cut=0, Int_t n_steps = 15, Double_t r_step = 10.0, Double_t z_cut = 27000.0, Int_t detector = 28, TString type = "moller", Int_t radialCut = 1, TString fname = "NULL"){ |
| 21 | + /* I want to take a given ROOT output file and plot a specific detector's |
| 22 | + * hit.x/y/r/theta, etc. distributions (weighted by rate) |
| 23 | + * and using cuts that are set at the z_cut position, cutting out all part.tjr |
| 24 | + * (calculated radii) tracks that are > r_outer_cut |
| 25 | + * |
| 26 | + * I would also like to have this loop over r_outer_cut from some initial |
| 27 | + * to final in fine steps and make a plot of the number of cut out/in epelastic/moller |
| 28 | + * particles (weighted by rate) vs. that cut radius |
| 29 | + * |
| 30 | + * I would also like to do a similar scan in z maybe... or at least compare at multiple |
| 31 | + * z cut places */ |
| 32 | + |
| 33 | + if (r_outer_cut==0) { |
| 34 | + std::cout<<"Usage: Load Data Root File to gDirectory, root> .x plotDetPlaneCuts(r_outer_cut=0, n_steps = 15, r_step = 10.0, z_cut = 27000.0, detector = 28, type = \"moller\")\n\n"<<std::endl; |
| 35 | + return; |
| 36 | + } |
| 37 | + |
| 38 | + TChain * T = new TChain("T"); |
| 39 | + if (fname == "NULL"){ |
| 40 | + T = (TChain*)gDirectory->Get("T"); |
| 41 | + } |
| 42 | + else { |
| 43 | + T->Add(fname); |
| 44 | + } |
| 45 | + FILE *fp; |
| 46 | + fp = fopen(Form("Output_cuts.txt"),"a"); |
| 47 | + |
| 48 | + Int_t n_events = 100000; |
| 49 | + Int_t start_event_n = 1; |
| 50 | + Double_t z_cut_space = 1.0; // mm |
| 51 | + |
| 52 | + const Int_t n_canvases = 100; |
| 53 | + const Int_t n_plots = 1; |
| 54 | + Double_t norm_rate_weight_all[n_canvases*n_plots]; |
| 55 | + Double_t norm_rate_weight_in[n_canvases*n_plots]; |
| 56 | + Double_t norm_rate_weight_out[n_canvases*n_plots]; |
| 57 | + TCanvas *canvas[n_canvases]; |
| 58 | + TH1F *totalRate[n_canvases*n_plots]; |
| 59 | + TH1F *r_inside[n_canvases*n_plots]; |
| 60 | + TH1F *r_outside[n_canvases*n_plots]; |
| 61 | + TH2F *xy_inside[n_canvases*n_plots]; |
| 62 | + TH2F *xy_outside[n_canvases*n_plots]; |
| 63 | + |
| 64 | + for (Int_t n_canvas = 0; n_canvas < n_steps; n_canvas++) { |
| 65 | + r_outer_cut = r_outer_cut - r_step; |
| 66 | + canvas[n_canvas] = new TCanvas(Form("Plot of %s in detector %d hits for z plane %f, r > %f",type.Data(),detector,z_cut,r_outer_cut),Form("Plot of %s in detector %d hits for z plane %f, r > %f",type.Data(),detector,z_cut,r_outer_cut)); |
| 67 | + canvas[n_canvas]->Divide(n_plots*2,2); |
| 68 | + for (Int_t n_plot = 0; n_plot<n_plots; n_plot++){ |
| 69 | + std::cout<<"Plot number "<<n_canvas*n_plots+n_plot<<", n_canvas "<<n_canvas<<", n_plot "<<n_plot<<std::endl; |
| 70 | + |
| 71 | + canvas[n_canvas]->cd(n_plot+1); |
| 72 | + T->Draw(Form("hit.r[0]>>h_totalRate_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f)",detector,z_cut,z_cut_space,z_cut,z_cut_space),"",n_events,start_event_n); |
| 73 | + totalRate[n_canvas*n_plots + n_plot] = (TH1F*)gROOT->FindObject(Form("h_totalRate_%d",n_canvas*n_plots + n_plot)); |
| 74 | + norm_rate_weight_all[n_canvas*n_plots + n_plot] = totalRate[n_canvas*n_plots + n_plot]->Integral(); |
| 75 | + Printf("Normalization total rate weight = %f",norm_rate_weight_all[n_canvas*n_plots + n_plot]); |
| 76 | + |
| 77 | + // Radial |
| 78 | + if (radialCut == 1){ |
| 79 | + T->Draw(Form("hit.r[0]>>h_r_inside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && sqrt((part.tjx*part.tjx)+(part.tjy*part.tjy)) < %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 80 | + } |
| 81 | + else { |
| 82 | + // X cut |
| 83 | + T->Draw(Form("hit.r[0]>>h_r_inside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && part.tjx < %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 84 | + } |
| 85 | + r_inside[n_canvas*n_plots + n_plot] = (TH1F*)gROOT->FindObject(Form("h_r_inside_%d",n_canvas*n_plots + n_plot)); |
| 86 | + norm_rate_weight_in[n_canvas*n_plots + n_plot] = r_inside[n_canvas*n_plots + n_plot]->Integral(); |
| 87 | + |
| 88 | + canvas[n_canvas]->cd(n_plot+2); |
| 89 | + // Radial |
| 90 | + if (radialCut == 1){ |
| 91 | + T->Draw(Form("hit.r[0]>>h_r_outside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && sqrt((part.tjx*part.tjx)+(part.tjy*part.tjy)) > %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 92 | + // X cut |
| 93 | + } |
| 94 | + else { |
| 95 | + T->Draw(Form("hit.r[0]>>h_r_outside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && part.tjx > %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 96 | + } |
| 97 | + r_outside[n_canvas*n_plots + n_plot] = (TH1F*)gROOT->FindObject(Form("h_r_outside_%d",n_canvas*n_plots + n_plot)); |
| 98 | + norm_rate_weight_out[n_canvas*n_plots + n_plot] = r_outside[n_canvas*n_plots + n_plot]->Integral(); |
| 99 | + |
| 100 | + canvas[n_canvas]->cd(n_plot+3); |
| 101 | + // Radial |
| 102 | + if (radialCut == 1){ |
| 103 | + T->Draw(Form("hit.y[0]:hit.x[0]>>h_xy_inside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && sqrt((part.tjx*part.tjx)+(part.tjy*part.tjy)) < %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 104 | + } |
| 105 | + else { |
| 106 | + // X cut |
| 107 | + T->Draw(Form("hit.y[0]:hit.x[0]>>h_xy_inside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && part.tjx < %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 108 | + } |
| 109 | + xy_inside[n_canvas*n_plots + n_plot] = (TH2F*)gROOT->FindObject(Form("h_xy_inside_%d",n_canvas*n_plots + n_plot)); |
| 110 | + |
| 111 | + canvas[n_canvas]->cd(n_plot+4); |
| 112 | + // Radial |
| 113 | + if (radialCut == 1){ |
| 114 | + T->Draw(Form("hit.y[0]:hit.x[0]>>h_xy_outside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && sqrt((part.tjx*part.tjx)+(part.tjy*part.tjy)) > %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 115 | + } |
| 116 | + // X cut |
| 117 | + else { |
| 118 | + T->Draw(Form("hit.y[0]:hit.x[0]>>h_xy_outside_%d",n_canvas*n_plots + n_plot),Form("rate*(hit.det[0] == %d && part.tjz > %f-%f && part.tjz < %f+%f && part.tjx > %f)",detector,z_cut,z_cut_space,z_cut,z_cut_space,r_outer_cut),"",n_events,start_event_n); |
| 119 | + } |
| 120 | + xy_outside[n_canvas*n_plots + n_plot] = (TH2F*)gROOT->FindObject(Form("h_xy_outside_%d",n_canvas*n_plots + n_plot)); |
| 121 | + |
| 122 | + // CAMGUIN feature: writeFile_h takes the name of the variable, the value to set it to, the "run number" which is a unique identifier, the split number (which is another unique identifier), and the number of runs (another one) |
| 123 | + //writeFile_h("r_outer_cut", r_outer_cut, r_outer_cut, z_cut, detector*1.0); |
| 124 | + //writeFile_h("z_cut", z_cut, r_outer_cut, z_cut, detector*1.0); |
| 125 | + //writeFile_h("detector", detector*1.0, r_outer_cut, z_cut, detector*1.0); |
| 126 | + //writeFile_h(Form("%s_norm_rate_weight_all",type.Data()), norm_rate_weight_all[n_canvas*n_plots + n_plot], r_outer_cut, z_cut, detector*1.0); |
| 127 | + //writeFile_h(Form("%s_norm_rate_weight_in",type.Data()), norm_rate_weight_in[n_canvas*n_plots + n_plot], r_outer_cut, z_cut, detector*1.0); |
| 128 | + //writeFile_h(Form("%s_norm_rate_weight_out",type.Data()), norm_rate_weight_out[n_canvas*n_plots + n_plot], r_outer_cut, z_cut, detector*1.0); |
| 129 | + |
| 130 | + if(fp!=NULL){ |
| 131 | + fprintf(fp,"r_outer_cut = %f, z_cut = %f, detector = %d, %s_norm_rate_weight_all = %f, %s_norm_rate_weight_in = %f, %s_norm_rate_weight_out = %f\n",r_outer_cut,z_cut,detector,type.Data(),norm_rate_weight_all[n_canvas*n_plots + n_plot],type.Data(),norm_rate_weight_in[n_canvas*n_plots + n_plot],type.Data(),norm_rate_weight_out[n_canvas*n_plots + n_plot]); |
| 132 | + } |
| 133 | + } |
| 134 | + if (n_canvas == 0 && n_steps != 1){ |
| 135 | + canvas[n_canvas]->SaveAs(Form("%s_norm_rate_weight.pdf(",type.Data())); |
| 136 | + } |
| 137 | + else if (n_canvas < n_steps || n_steps==1) { |
| 138 | + canvas[n_canvas]->SaveAs(Form("%s_norm_rate_weight.pdf",type.Data())); |
| 139 | + } |
| 140 | + else { |
| 141 | + canvas[n_canvas]->SaveAs(Form("%s_norm_rate_weight.pdf)",type.Data())); |
| 142 | + } |
| 143 | + } |
| 144 | + |
| 145 | + fclose(fp); |
| 146 | +} |
0 commit comments