00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "rootAngleHist.h"
00010 #include <math.h>
00011 #include <iostream>
00012 #include <iomanip>
00013 #include <fstream>
00014
00015
00016 rootAngleHist::rootAngleHist(int bins) : num_bins(bins)
00017 {
00018 theta_hist = new rootHist(bins);
00019 theta_sigma_hist = new rootHist(bins);
00020 theta_raw_hist = new rootHist(bins);
00021 phi_hist = new rootHist(bins);
00022 phi_sigma_hist = new rootHist(bins);
00023 phi_raw_hist = new rootHist(bins);
00024 currentType = linear;
00025 graphTitle = "No Title";
00026 PhiXLabel = "";
00027 PhiYLabel = "";
00028 ThetaXLabel = "";
00029 ThetaYLabel = "";
00030 }
00031
00032 rootAngleHist::~rootAngleHist(void)
00033 {
00034 delete theta_hist;
00035 delete theta_sigma_hist;
00036 delete theta_raw_hist;
00037 delete phi_hist;
00038 delete phi_sigma_hist;
00039 delete phi_raw_hist;
00040 }
00041
00042 rootAngleHist::rootAngleHist(const rootAngleHist& oldHist) : num_bins(oldHist.num_bins)
00043 {
00044 theta_hist = new rootHist(num_bins);
00045 theta_sigma_hist = new rootHist(num_bins);
00046 theta_raw_hist = new rootHist(num_bins);
00047 phi_hist = new rootHist(num_bins);
00048 phi_sigma_hist = new rootHist(num_bins);
00049 phi_raw_hist = new rootHist(num_bins);
00050
00051 for(int i = 0; i < num_bins; i++)
00052 {
00053 double oldCount = (oldHist.theta_hist)->retrieveBin(i);
00054 theta_hist->updateBin(i, oldCount);
00055
00056 oldCount = (oldHist.theta_sigma_hist)->retrieveBin(i);
00057 theta_sigma_hist->updateBin(i, oldCount);
00058
00059 oldCount = (oldHist.theta_raw_hist)->retrieveBin(i);
00060 theta_raw_hist->updateBin(i, oldCount);
00061
00062 oldCount = (oldHist.phi_hist)->retrieveBin(i);
00063 phi_hist->updateBin(i, oldCount);
00064
00065 oldCount = (oldHist.phi_sigma_hist)->retrieveBin(i);
00066 phi_sigma_hist->updateBin(i, oldCount);
00067
00068 oldCount = (oldHist.phi_raw_hist)->retrieveBin(i);
00069 phi_raw_hist->updateBin(i, oldCount);
00070 }
00071 currentType = oldHist.currentType;
00072 graphTitle = oldHist.graphTitle;
00073 ThetaXLabel = oldHist.ThetaXLabel;
00074 ThetaYLabel = oldHist.ThetaYLabel;
00075 PhiXLabel = oldHist.PhiXLabel;
00076 PhiYLabel = oldHist.PhiYLabel;
00077 }
00078
00079 void rootAngleHist::setTitle(std::string title)
00080 {
00081 graphTitle = title;
00082 }
00083
00084 void rootAngleHist::setPhiXLabel(std::string label)
00085 {
00086 PhiXLabel = label;
00087 }
00088
00089 void rootAngleHist::setPhiYLabel(std::string label)
00090 {
00091 PhiYLabel = label;
00092 }
00093
00094 void rootAngleHist::setThetaXLabel(std::string label)
00095 {
00096 ThetaXLabel = label;
00097 }
00098
00099 void rootAngleHist::setThetaYLabel(std::string label)
00100 {
00101 ThetaYLabel = label;
00102 }
00103
00104 void rootAngleHist::setGraphType(const char *graph_type)
00105 {
00106 if(0 == strcmp(graph_type,"linear") || 0 == strcmp(graph_type,"semilogx"))
00107 currentType = linear;
00108 else if(0 == strcmp(graph_type,"log") || 0 == strcmp(graph_type,"semilogy"))
00109 currentType = log;
00110 else
00111 std::cerr << "ERROR: Invalid Graph Type" << std::endl;
00112 }
00113
00114 void rootAngleHist::storeTheta(double cos_theta_value)
00115 {
00116 cos_theta_value += 1;
00117
00118
00119 int binNumber = floor((cos_theta_value)/2 * num_bins);
00120 double count = theta_raw_hist->retrieveBin(binNumber);
00121 theta_raw_hist->updateBin(binNumber,++count);
00122 }
00123
00124 void rootAngleHist::storePhi(double phi_value)
00125 {
00126 int binNumber = floor(phi_value / (2*M_PI) * num_bins);
00127 double count = phi_raw_hist->retrieveBin(binNumber);
00128 phi_raw_hist->updateBin(binNumber,++count);
00129 }
00130
00131 double rootAngleHist::retrievePhiCount(int binNumber)
00132 {
00133 return phi_raw_hist->retrieveBin(binNumber);
00134 }
00135
00136 double rootAngleHist::retrievePhiFlux(int binNumber)
00137 {
00138 return phi_hist->retrieveBin(binNumber);
00139 }
00140
00141 double rootAngleHist::retrievePhiFluxUncertainty(int binNumber)
00142 {
00143 return phi_sigma_hist->retrieveBin(binNumber);
00144 }
00145
00146 double rootAngleHist::retrieveThetaCount(int binNumber)
00147 {
00148 return theta_raw_hist->retrieveBin(binNumber);
00149 }
00150
00151 double rootAngleHist::retrieveThetaFlux(int binNumber)
00152 {
00153 return theta_hist->retrieveBin(binNumber);
00154 }
00155
00156 double rootAngleHist::retrieveThetaFluxUncertainty(int binNumber)
00157 {
00158 return theta_sigma_hist->retrieveBin(binNumber);
00159 }
00160
00161 double rootAngleHist::retrieveThetaAngle(int binNumber)
00162 {
00163 return (binNumber+0.5) / num_bins * 2 - 1;
00164 }
00165
00166 double rootAngleHist::retrievePhiAngle(int binNumber)
00167 {
00168 return ((binNumber+0.5) / num_bins * 2 * M_PI);
00169 }
00170
00171
00172 void rootAngleHist::apply(double scale_factor)
00173 {
00174 for(int i = 0; i < num_bins; ++i)
00175 {
00176 double current_flux = theta_hist->retrieveBin(i);
00177 double current_count = theta_raw_hist->retrieveBin(i);
00178 current_flux += (current_count * scale_factor);
00179 theta_hist->updateBin(i,current_flux);
00180
00181
00182 current_flux = sqrt(pow(theta_sigma_hist->retrieveBin(i),2)+current_count) * scale_factor;
00183 theta_sigma_hist->updateBin(i,current_flux);
00184
00185
00186 theta_raw_hist->updateBin(i,0);
00187
00188 current_flux = phi_hist->retrieveBin(i);
00189 current_count = phi_raw_hist->retrieveBin(i);
00190 current_flux += (current_count * scale_factor);
00191 phi_hist->updateBin(i,current_flux);
00192
00193
00194 current_flux = sqrt(pow(phi_sigma_hist->retrieveBin(i),2)+current_count) * scale_factor;
00195 phi_sigma_hist->updateBin(i,current_flux);
00196
00197
00198 phi_raw_hist->updateBin(i,0);
00199 }
00200 }
00201
00202 void rootAngleHist::reset(void)
00203 {
00204 for(int i = 0; i < num_bins; ++i)
00205 {
00206 phi_hist->updateBin(i,0);
00207 phi_sigma_hist->updateBin(i,0);
00208 phi_raw_hist->updateBin(i,0);
00209 theta_hist->updateBin(i,0);
00210 theta_sigma_hist->updateBin(i,0);
00211 theta_raw_hist->updateBin(i,0);
00212 }
00213
00214 graphTitle = "No Title";
00215 PhiXLabel = "";
00216 PhiYLabel = "";
00217 ThetaXLabel = "";
00218 ThetaYLabel = "";
00219 }
00220
00221 void rootAngleHist::draw(double scale_factor, std::string mode, int current_plot, int total_plots)
00222 {
00223 char *window_title = "Graph Window";
00224
00225 if(current_plot >= total_plots)
00226 {
00227 std::cerr << "Error: Invalid plot number" << std::endl;
00228 exit(0);
00229 }
00230
00231 if(current_plot == 0)
00232 {
00233 std::ofstream out_file;
00234
00235 if(mode != "end")
00236 {
00237 out_file.open("graph.cxx");
00238
00239 if(false == out_file.is_open())
00240 {
00241 std::cerr << "Unable to open temporary file for writing." << std::endl;
00242 exit(0);
00243 }
00244
00245 out_file <<
00246 "{\n"
00247 " gROOT->Reset();\n";
00248 }
00249 else
00250 {
00251 out_file.open("graph.cxx", std::ios::app);
00252 }
00253
00254 out_file <<
00255 " char *angle_window_title = \"" << window_title << "\";\n"
00256 " char *theta_graph_title = \"" << "Particle Flux vs. Zenith Angle" << "\";\n"
00257 " char *phi_graph_title = \"" << "Particle Flux vs. Azimuth Angle" << "\";\n"
00258 " char *theta_y_label = \"" << ThetaYLabel << "\";\n"
00259 " char *theta_x_label = \"" << ThetaXLabel << "\";\n"
00260 " char *phi_y_label = \"" << PhiYLabel << "\";\n"
00261 " char *phi_x_label = \"" << PhiXLabel << "\";\n"
00262 " int num_bins = "<< num_bins << ";\n"
00263 " double theta_angle[num_bins];\n"
00264 " double e_theta_angle[num_bins];\n"
00265 " double phi_angle[num_bins];\n"
00266 " double e_phi_angle[num_bins];\n"
00267
00268 " for(int i = 0; i < num_bins; i++) {\n"
00269 " theta_angle[i] = (i+0.5) / num_bins * 2 - 1;\n"
00270 " e_theta_angle[i] = 0;\n"
00271 " phi_angle[i] = (double(i)+0.5)/double(num_bins)*360;\n"
00272 " e_phi_angle[i] = 0;\n"
00273 " }\n"
00274
00275 " // name, title, wtopx, wtopy, ww, wh\n"
00276 " c2 = new TCanvas(\"c2\",angle_window_title, 200, 200, 700, 500);\n"
00277 " c2->SetGrid();\n"
00278 " c2->GetFrame()->SetFillColor(21);\n"
00279 " c2->GetFrame()->SetBorderSize(12);\n"
00280 " c3 = new TCanvas(\"c3\",angle_window_title, 200, 200, 700, 500);\n"
00281 " c3->SetGrid();\n"
00282 " c3->GetFrame()->SetFillColor(21);\n"
00283 " c3->GetFrame()->SetBorderSize(12);\n";
00284
00285 out_file <<
00286 " theta_leg = new TLegend(0.73,0.83,0.99,0.99);\n"
00287 " phi_leg = new TLegend(0.73,0.83,0.99,0.99);\n"
00288 " double scale_factor" << current_plot << " = " << scale_factor << ";\n"
00289 " double theta_count" << current_plot << "[] = {\n";
00290 {for(int i = 0; i < num_bins; i++)
00291 out_file << std::setw(5) << retrieveThetaFlux(i) << (i%5==4? ",\n" : ",");
00292 }
00293
00294 out_file <<
00295 " };\n"
00296 " double phi_count" << current_plot << "[] = {\n";
00297 {for(int i = 0; i < num_bins; i++)
00298 out_file << std::setw(5) << retrievePhiFlux(i) << (i%5==4? ",\n" : ",");
00299 }
00300
00301 out_file <<
00302 " };\n"
00303 " double upper_angle;"
00304 " double lower_angle;"
00305 " double e_theta_count" << current_plot << "[] = {\n";
00306 {for(int i = 0; i < num_bins; i++)
00307 out_file << std::setw(5) << retrieveThetaFluxUncertainty(i) << (i%5==4? ",\n" : ",");
00308 }
00309 out_file <<
00310 " };\n"
00311 " double e_phi_count" << current_plot << "[] = {\n";
00312 {for(int i = 0; i < num_bins; i++)
00313 out_file << std::setw(5) << retrievePhiFluxUncertainty(i) << (i%5==4? ",\n" : ",");
00314 }
00315 out_file <<
00316 " };\n"
00317 " for(int i = 0; i < num_bins; i++) {\n"
00318 " e_theta_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00319 " e_phi_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00320 " theta_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00321 " phi_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00322 " }\n"
00323 " c2->cd();"
00324 " theta_gr" << current_plot << " = new TGraphErrors(num_bins,theta_angle,theta_count" << current_plot
00325 << ",e_theta_angle,e_theta_count" << current_plot << ");\n"
00326 " theta_gr" << current_plot << "->SetTitle(theta_graph_title);\n"
00327 " theta_gr" << current_plot << "->SetMarkerColor(2);\n"
00328 " theta_gr" << current_plot << "->SetMarkerStyle(21);\n"
00329 " theta_leg->AddEntry(theta_gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n"
00330 " c3->cd();"
00331 " phi_gr" << current_plot << " = new TGraphErrors(num_bins,phi_angle,phi_count" << current_plot
00332 << ",e_phi_angle,e_phi_count" << current_plot << ");\n"
00333 " phi_gr" << current_plot << "->SetTitle(phi_graph_title);\n"
00334 " phi_gr" << current_plot << "->SetMarkerColor(2);\n"
00335 " phi_gr" << current_plot << "->SetMarkerStyle(21);\n"
00336 " phi_leg->AddEntry(phi_gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n";
00337
00338 out_file.close();
00339 }
00340 else if(current_plot <= total_plots - 1)
00341 {
00342 std::ofstream out_file("graph.cxx", std::ios::app);
00343 out_file <<
00344 " double scale_factor" << current_plot << " = " << scale_factor << ";\n"
00345 " double theta_count" << current_plot << "[] = {\n";
00346 {for(int i = 0; i < num_bins; i++)
00347 out_file << std::setw(5) << retrieveThetaFlux(i) << (i%5==4? ",\n" : ",");
00348 }
00349
00350 out_file <<
00351 " };\n"
00352 " double phi_count" << current_plot << "[] = {\n";
00353 {for(int i = 0; i < num_bins; i++)
00354 out_file << std::setw(5) << retrievePhiFlux(i) << (i%5==4? ",\n" : ",");
00355 }
00356
00357 out_file <<
00358 " };\n"
00359 " double e_theta_count" << current_plot << "[] = {\n";
00360 {for(int i = 0; i < num_bins; i++)
00361 out_file << std::setw(5) << retrieveThetaFluxUncertainty(i) << (i%5==4? ",\n" : ",");
00362 }
00363 out_file <<
00364 " };\n"
00365 " double e_phi_count" << current_plot << "[] = {\n";
00366 {for(int i = 0; i < num_bins; i++)
00367 out_file << std::setw(5) << retrievePhiFluxUncertainty(i) << (i%5==4? ",\n" : ",");
00368 }
00369 out_file <<
00370 " };\n"
00371 " for(i = 0; i < num_bins; i++) {\n"
00372 " e_theta_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00373 " e_phi_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00374 " theta_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00375 " phi_count" << current_plot << "[i] *= scale_factor" << current_plot << ";\n"
00376 " }\n"
00377 " c2->cd();"
00378 " theta_gr" << current_plot << " = new TGraphErrors(num_bins,theta_angle,theta_count" << current_plot
00379 << ",e_theta_angle,e_theta_count" << current_plot << ");\n"
00380 " theta_gr" << current_plot << "->SetMarkerColor(" << 2+current_plot << ");\n"
00381 " theta_gr" << current_plot << "->SetMarkerStyle(" << 21 << ");\n"
00382 " theta_leg->AddEntry(theta_gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n"
00383 " c3->cd();"
00384 " phi_gr" << current_plot << " = new TGraphErrors(num_bins,phi_angle,phi_count" << current_plot
00385 << ",e_phi_angle,e_phi_count" << current_plot << ");\n"
00386 " phi_gr" << current_plot << "->SetMarkerColor(" << 2+current_plot << ");\n"
00387 " phi_gr" << current_plot << "->SetMarkerStyle(" << 21 << ");\n"
00388 " phi_leg->AddEntry(phi_gr" << current_plot << ",\"" << graphTitle << "\",\"P\");\n";
00389 out_file.close();
00390 }
00391
00392 if(current_plot >= total_plots - 1)
00393 {
00394 std::ofstream out_file("graph.cxx", std::ios::app);
00395
00396 out_file <<
00397 " c2->cd();"
00398 " double theta_max_count = 0;\n"
00399 " double theta_min_count = 1e12;\n"
00400 " for(int i = 0; i < num_bins; i++)\n"
00401 " {\n";
00402 {for(int plot = 0; plot < total_plots; plot++) {
00403 out_file <<
00404 " if(theta_count" << plot << "[i] > theta_max_count)\n"
00405 " theta_max_count = theta_count" << plot << "[i];\n"
00406 " if(theta_count" << plot << "[i] < theta_min_count && theta_count" << plot << "[i] > 0)\n"
00407 " theta_min_count = theta_count" << plot << "[i];\n";
00408 }}
00409 out_file <<
00410 " }\n"
00411 " double theta_angle_limits[] = {-1,1};\n"
00412 " double theta_count_limits[] = {0,theta_max_count};\n"
00413 " int theta_bin_limits = 2;\n"
00414 " theta_graph0 = new TGraph(theta_bin_limits,theta_angle_limits,theta_count_limits);\n"
00415 " theta_graph0->SetTitle(theta_graph_title);\n"
00416 " theta_graph0->Draw(\"AP\");\n"
00417 " TAxis *theta_ax = theta_graph0->GetXaxis();\n"
00418 " TAxis *theta_ay = theta_graph0->GetYaxis();\n"
00419 " theta_ay->SetTitle(theta_y_label); theta_ay->CenterTitle(1);\n"
00420 " theta_ax->SetLimits(-1,1);\n"
00421 " theta_ax->SetTitle(theta_x_label); theta_ax->CenterTitle(1); \n"
00422 " theta_ax->SetTitleOffset(1.2);\n";
00423
00424 {for(int plot = 0; plot < total_plots; plot++) {
00425 out_file <<
00426 " theta_gr" << plot << "->Draw(\"P\");\n";
00427 }}
00428
00429 out_file <<
00430 " theta_leg->Draw();\n"
00431 " c2->Modified();\n"
00432 " c2->Update();\n";
00433
00434 out_file <<
00435 " c3->cd();"
00436 " double phi_max_count = 0;\n"
00437 " double phi_min_count = 1e12;\n"
00438 " for(int i = 0; i < num_bins; i++)\n"
00439 " {\n";
00440 {for(int plot = 0; plot < total_plots; plot++) {
00441 out_file <<
00442 " if(phi_count" << plot << "[i] > phi_max_count)\n"
00443 " phi_max_count = phi_count" << plot << "[i];\n"
00444 " if(phi_count" << plot << "[i] < phi_min_count && phi_count" << plot << "[i] > 0)\n"
00445 " phi_min_count = phi_count" << plot << "[i];\n";
00446 }}
00447 out_file <<
00448 " }\n"
00449 " double phi_angle_limits[] = {0,360};\n"
00450 " double phi_count_limits[] = {0,phi_max_count};\n"
00451 " int phi_bin_limits = 2;\n"
00452 " phi_graph0 = new TGraph(phi_bin_limits,phi_angle_limits,phi_count_limits);\n"
00453 " phi_graph0->SetTitle(phi_graph_title);\n"
00454 " phi_graph0->Draw(\"AP\");\n"
00455 " TAxis *phi_ax = phi_graph0->GetXaxis();\n"
00456 " TAxis *phi_ay = phi_graph0->GetYaxis();\n"
00457 " phi_ay->SetTitle(phi_y_label); phi_ay->CenterTitle(1);\n"
00458 " phi_ax->SetLimits(0,360);\n"
00459 " phi_ax->SetTitle(phi_x_label); phi_ax->CenterTitle(1); \n"
00460 " phi_ax->SetTitleOffset(1.2);\n";
00461
00462 {for(int plot = 0; plot < total_plots; plot++) {
00463 out_file <<
00464 " phi_gr" << plot << "->Draw(\"P\");\n";
00465 }}
00466
00467 out_file <<
00468 " phi_leg->Draw();\n"
00469 " c3->Modified();\n"
00470 " c3->Update();\n";
00471
00472 if(mode != "begin")
00473 out_file << "\n}\n";
00474
00475 out_file.close();
00476
00477 if(mode != "begin")
00478 system("root -l graph.cxx");
00479 }
00480 }
00481
00482 void rootAngleHist::writeFile(double scale_factor, std::ostream& out_file)
00483 {
00484
00485 }
