ampcomparison.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <vector>
4 #include <string>
5 #include <TCanvas.h>
6 #include <TH2F.h>
7 #include <TMath.h>
8 #include <TFile.h>
9 #include <TGraph.h>
10 #include <TROOT.h>
11 #include <TStyle.h>
12 #include <TColor.h>
13 #include <TTimeStamp.h>
14 #define AUGER_RADIO_ENABLED
15 #include "RecEventFile.h"
16 #include "RecEvent.h"
17 #include "RdEvent.h"
18 #include "DetectorGeometry.h"
19 #include "SDEvent.h"
20 #include "DetectorGeometry.h"
21 #include "SdRecStation.h"
22 #include "ShowerRRecDataQuantities.h"
23 #include "StationRRecDataQuantities.h"
24 #include "PolarHist.h"
25 #include "TLegend.h"
26 #include "TPaletteAxis.h"
27 
28 //Added
29 #include <TVector3.h>
30 #include <TProfile.h>
31 #include <TLatex.h>
32 #include <TPaveStats.h>
33 
34 #include <RdTrace.h>
35 
36 using namespace std;
37 
38 class AntennaData {
39  public:
40  int runnumber;
41  int eventid;
42  int stationid;
43  double peakamplitude;
44  double theta;
45  double energy;
46  bool contained;
47  double axis_distance;
48 
49  bool operator==(const AntennaData &other);
50  bool operator!=(const AntennaData &other);
51 };
52 
53 
54 bool AntennaData::operator==(const AntennaData &other){
55  return (this->runnumber == other.runnumber) && (this->eventid == other.eventid) && (this->stationid == other.stationid);
56 }
57 
58 bool AntennaData::operator!=(const AntennaData &other){
59  return !(*this == other);
60 }
61 
62 vector<AntennaData> readAmplitudeFile(string filename){
63  int run, event, station;
64  double amp, theta, energy, dist;
65  bool contained;
66 
67  vector<AntennaData> entries;
68  ifstream fin (filename.c_str());
69  while (fin >> run >> event >> station >> amp >> theta >> energy >> contained >> dist){
70  AntennaData a;
71  a.runnumber = run;
72  a.eventid = event;
73  a.stationid = station;
74  a.peakamplitude = amp;
75  a.theta = theta;
76  a.energy = energy;
77  a.contained = contained;
78  a.axis_distance = dist;
79 
80  entries.push_back(a);
81  }
82  cout << "read " << entries.size() << " entries from " << filename << endl;
83  return entries;
84 }
85 
86 int main() {
87  TH2F* Amplitude_comparison_log = new TH2F("Amplitude_comparison_log", "Amplitude_comparison_log",
88  100, 0, 5, 100, 0, 5);
89  TProfile* Amplitude_comparison_profile_log = new TProfile("Amplitude_comparison_profile_log",
90  "Amplitude_comparison_profile_log",
91  25, 0, 5, 0, 5, "S");
92 
93  TH1F* Hist_Deviations = new TH1F("Hist_Amp_Deviations", "Hist_Amp_Deviations",
94  40, -2., 2.);
95  TH1F* Hist_Deviations_Log = new TH1F("Hist_Amp_Deviations_Log", "Hist_Amp_Deviations_Log",
96  40, -1., 1.);
97  TH1F* Hist_ClassicalDeviations = new TH1F("Hist_Amp_ClassicalDeviations", "Hist_Amp_ClassicalDeviations",
98  40, -1., 3.);
99 
100  TH1F* Hist_DiffDiag=new TH1F("Hist_DiffDiag", "Hist_DiffDiag", 100, -1.5, 1.5);
101 
102  vector<AntennaData> sim = readAmplitudeFile("Amplitude_sim.txt");
103  vector<AntennaData> data = readAmplitudeFile("Amplitude_data.txt");
104 
105  for(int simid=0; simid<sim.size(); simid++) {
106  for(int dataid=0; dataid<data.size(); dataid++) {
107  if(sim[simid] == data[dataid]){
108  //cout << sim[simid].runnumber << "." << sim[simid].eventid << endl;
109  double simamp = sim[simid].peakamplitude;
110  double dataamp = data[dataid].peakamplitude;
111 
112  Amplitude_comparison_log->Fill(log10(dataamp),log10(simamp));
113  Amplitude_comparison_profile_log->Fill(log10(dataamp),log10(simamp));
114  Hist_Deviations->Fill(2.*(simamp-dataamp)/(simamp+dataamp));
115  Hist_ClassicalDeviations->Fill((simamp-dataamp)/dataamp);
116  Hist_Deviations_Log->Fill(log10(simamp/dataamp));
117 
118  double distanceToDiagonal = (simamp-dataamp)/sqrt(2);
119  double mean = (simamp+dataamp)/2;
120  Hist_DiffDiag->Fill(distanceToDiagonal/mean);
121  }
122  }
123  }
124 
125  //Drawing:
126  gStyle->SetFitFormat("7.2g");
127  gStyle->SetStatFormat("7.2g");
128 
129  TLatex Tl;
130  float y, x1, x2;
131  y = 0.90;
132  x1 = 0.07;
133  x2 = x1+0.2;
134  cout<<Tl.DrawLatex(x2, y, "#mu")<<endl;
135 
136  TLatex* text = new TLatex();
137  text->DrawLatex(0.45,0.05," #mu (GeV) ");
138 
139  TF1 *f1 = new TF1("f1","x",0,10000);
140  f1->SetLineColor(kGray+2);
141  f1->SetLineStyle(2);
142 
143  TCanvas *c2 = new TCanvas("c2", "", 1450, 1200);
144  c2->SetLeftMargin(0.18);
145  c2->SetRightMargin(0.18);
146  c2->SetTopMargin(0.06);
147  c2->SetBottomMargin(0.14);
148  c2->SetFixedAspectRatio();
149 
150 // Amplitude_comparison_log->SetTitle("Amplitude Comparison: #theta: 62#circ - 80#circ, E: 10^{17.4} eV - 10^{19.5} eV");
151  Amplitude_comparison_log->SetTitle("");
152 // Amplitude_comparison_log->SetTitleSize(0.04);
153  Amplitude_comparison_log->GetXaxis()->SetTitle("Data: log_{10}(Amplitude / (#muV/m))");
154  Amplitude_comparison_log->GetXaxis()->SetLabelFont(42);
155  Amplitude_comparison_log->GetXaxis()->SetLabelOffset(0.015);
156  Amplitude_comparison_log->GetXaxis()->SetLabelSize(0.05);
157  Amplitude_comparison_log->GetXaxis()->SetTitleSize(0.055);
158  Amplitude_comparison_log->GetXaxis()->SetTitleOffset(1.15);
159  Amplitude_comparison_log->GetXaxis()->SetTitleFont(42);
160  Amplitude_comparison_log->GetXaxis()->SetRangeUser(1.5,3.9);
161  Amplitude_comparison_log->GetXaxis()->SetNdivisions(505);
162  Amplitude_comparison_log->GetYaxis()->SetTitle("Simulation: log_{10}(Amplitude / (#muV/m))");
163  Amplitude_comparison_log->GetYaxis()->SetLabelFont(42);
164  Amplitude_comparison_log->GetYaxis()->SetLabelOffset(0.01);
165  Amplitude_comparison_log->GetYaxis()->SetLabelSize(0.05);
166  Amplitude_comparison_log->GetYaxis()->SetTitleSize(0.055);
167  Amplitude_comparison_log->GetYaxis()->SetTitleOffset(1.25);
168  Amplitude_comparison_log->GetYaxis()->SetTitleFont(42);
169  Amplitude_comparison_log->GetYaxis()->SetRangeUser(1.5,3.9);
170  Amplitude_comparison_log->GetYaxis()->SetNdivisions(505);
171  Amplitude_comparison_log->GetZaxis()->SetTitle("Entries");
172  Amplitude_comparison_log->GetZaxis()->SetTitleFont(42);
173  Amplitude_comparison_log->GetZaxis()->SetTitleSize(0.055);
174  Amplitude_comparison_log->GetZaxis()->SetLabelFont(42);
175  Amplitude_comparison_log->GetZaxis()->SetLabelOffset(0.01);
176  Amplitude_comparison_log->GetZaxis()->SetLabelSize(0.05);
177  Amplitude_comparison_log->SetMarkerSize(0.3);
178  Amplitude_comparison_log->SetMarkerStyle(20);
179  Amplitude_comparison_log->Draw("COLZ");
180  Amplitude_comparison_profile_log->SetLineColor(kRed);
181  Amplitude_comparison_profile_log->SetLineWidth(4.0);
182  Amplitude_comparison_log->SetStats(kFALSE);
183 
184 // Amplitude_comparison_profile_log->Draw("P, SAME");
185 
186  f1->SetLineWidth(3.);
187  f1->Draw("SAME");
188 
189  TF1 *f2 = new TF1("f2","x+0.08",0,10000);
190  f2->SetLineColor(kGray+2);
191  f2->SetLineWidth(3.);
192  f2->SetLineStyle(10);
193 
194  TF1 *f3 = new TF1("f3","x-0.1",0,10000);
195  f3->SetLineColor(kGray+2);
196  f3->SetLineStyle(10);
197  f3->SetLineWidth(3.);
198 
199  c2->SaveAs("Amplitude_comparison_log_proton.pdf");
200 
201  TFile *f = new TFile("ampdeviations.root","RECREATE");
202  Hist_Deviations->Write();
203  f->Close();
204 
205  TCanvas Histo("Histo","Histo",1800,1000);
206  Histo.SetLeftMargin(0.12);
207  Histo.SetRightMargin(0.04);
208  Histo.SetTopMargin(0.06);
209  Histo.SetBottomMargin(0.2);
210 
211  Hist_Deviations_Log->GetXaxis()->SetTitle("log10(Amplitude_sim/Amplitude_data)");
212  Hist_Deviations_Log->GetYaxis()->SetTitle("Number of entries");
213  //Hist_Energy_sd->SetMarkerStyle(20);
214  //Hist_Energy_sd->SetMarkerSize(0.4);
215  //Hist_Energy_sd->SetMarkerColor(kBlack);
216  //Hist_Energy_sd->GetXaxis()->SetTitleOffset(1.2);
217  //Hist_Energy_sd->GetYaxis()->SetTitleOffset(1.2);
218  Hist_Deviations_Log->SetTitle("");
219  Hist_Deviations_Log->SetFillColor(4);
220  Hist_Deviations_Log->Draw("HIST");
221  Histo.SaveAs("Amp_Deviations_Log.pdf");
222 
223  Hist_Deviations->GetXaxis()->SetTitle("2#upoint(A#lower[0.5]{#scale[0.8]{sim}} - A#lower[0.5]{#scale[0.8]{data}})/(A#lower[0.5]{#scale[0.8]{sim}} + A#lower[0.5]{#scale[0.8]{data}})");
224  Hist_Deviations->SetStats(kFALSE);
225  Hist_Deviations->GetXaxis()->SetLabelFont(42);
226  Hist_Deviations->GetXaxis()->SetLabelOffset(0.01);
227  Hist_Deviations->GetXaxis()->SetLabelSize(0.075);
228  Hist_Deviations->GetXaxis()->SetTitleSize(0.08);
229  Hist_Deviations->GetXaxis()->SetTitleOffset(1.15);
230  Hist_Deviations->GetXaxis()->SetTitleFont(42);
231  Hist_Deviations->GetYaxis()->SetTitle("Entries");
232  Hist_Deviations->GetYaxis()->SetLabelFont(42);
233  Hist_Deviations->GetYaxis()->SetLabelOffset(0.01);
234  Hist_Deviations->GetYaxis()->SetLabelSize(0.075);
235  Hist_Deviations->GetYaxis()->SetTitleSize(0.08);
236  Hist_Deviations->GetYaxis()->SetTitleOffset(0.8);
237  Hist_Deviations->GetYaxis()->SetTitleFont(42);
238  Hist_Deviations->SetTitle("");
239  Hist_Deviations->SetFillColor(4);
240  Hist_Deviations->Draw("HIST");
241  Histo.SaveAs("Amp_Deviations.pdf");
242 
243  cout << "number of entries in antenna comparison: "
244  << Hist_Deviations->GetEntries() << endl;
245 
246  Hist_ClassicalDeviations->GetXaxis()->SetTitle("(A_sim-A_data)/A_data");
247  Hist_ClassicalDeviations->GetYaxis()->SetTitle("Number of entries");
248  Hist_ClassicalDeviations->SetTitle("");
249  Hist_ClassicalDeviations->SetFillColor(4);
250  Hist_ClassicalDeviations->Draw("HIST");
251  Histo.SaveAs("Amp_ClassicalDeviations.pdf");
252 
253  Hist_DiffDiag->GetXaxis()->SetTitle("Distance to diagonal");
254  Hist_DiffDiag->GetYaxis()->SetTitle("Number of entries");
255  Hist_DiffDiag->SetTitle("");
256  Hist_DiffDiag->SetFillColor(4);
257  Hist_DiffDiag->Draw("HIST");
258  Histo.SaveAs("Amp_DiffDiag.pdf");
259 }
AntennaData::peakamplitude
double peakamplitude
Definition: ampcomparison.cc:43
main
int main(int argc, char *argv[])
Definition: DBSync.cc:58
utl::operator==
bool operator==(const TimeStamp &ts, const TimeRange &tr)
Definition: TimeRange.h:97
run
fwk::VModule::ResultFlag(fwk::VModule::* run)(evt::Event &)
Definition: fwkPython.cc:59
readAmplitudeFile
vector< AntennaData > readAmplitudeFile(string filename)
Definition: ampcomparison.cc:62
data
uint16_t * data
Definition: dump1090.h:228
AntennaData::operator!=
bool operator!=(const AntennaData &other)
Definition: ampcomparison.cc:58
AntennaData::axis_distance
double axis_distance
Definition: ampcomparison.cc:47
filename
char * filename
Definition: dump1090.h:266
AntennaData::operator==
bool operator==(const AntennaData &other)
Definition: ampcomparison.cc:54
utl::operator!=
bool operator!=(const TimeStamp &ts, const TimeRange &tr)
Definition: TimeRange.h:100
, generated on Tue Sep 26 2023.