#include
#include
#include
#include
//#include "/usr/include/root/TH2.h"
#include "TH2.h"
#include
#include "TCanvas.h"
#include "TFile.h"
#include
#include "TLegend.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TLine.h"
#include "TLatex.h"
#include "TStyle.h"
//#define draw_strips main
//this requires ROOT, I generally just run root draw_strips.cc+, but there are probably better options
//This prints histograms of the covariance map based on number of hits, the expected covariance map based on diagonal values, and a map of how many stdevs observed is away from expected
//This also prints 1d histogram distributions of said number of stdevs between expected and observed
void draw_strips()
//int main(void)
{
// typical making a canvas and setting the color scheme, palette can be switched. A picture of what this palette makes is attached in the email
TCanvas *c = new TCanvas("pad1","pad1",150,10,990,660);;
gStyle->SetPalette(55);
//FIXME below is the file name of the csv made in covariance.cc, make sure it matches up. I just replace and recompile because I can't get inputs with argv to work with ROOT. You're welcome to try though
std::ifstream in("/home/rutca/c++/test/covariance_dac11_new_big.csv");
// the following are 4 2d histograms I make for each file, 240 by 240
TH2F * output = new TH2F("covariance","covariance",240, 0.5,240.5, 240, 0.5, 240.5);
// output is the data collected
TH2F * expected = new TH2F("expected","expected",240, 0.5,240.5, 240, 0.5, 240.5);
// expected is the extrapolation from the diagonals
// e.g the number of hits at 1,2 should be the number at 1,1 * 2,2 divided by (total events collected * 4). The 4 is because diagonals are double counted
TH2F * difference = new TH2F("difference","difference",240, 0.5,240.5, 240, 0.5, 240.5);
// difference is the number of standard deviations between expected and experimental (poisson distribution so stdev is sqrt of output)
// Be careful with high dac values, since some positions end with very low chances of having a hit, so if they get 1, it ends up being hundreds of deviations out from expected
TH2F * strange = new TH2F("outliers","outliers",240, 0.5,240.5, 240, 0.5, 240.5);
// this is a map of points with high difference values, so this would in theory map positions with high or low covariance
gStyle->SetPalette(55); //? not sure why this is here, could probably be removed
// the following are 1d
TH1F * ssa1 = new TH1F("ssa1 distribution","ssa1 distribution",100,-10,10);
// This is a histogram of the difference (the 2d hist) for the first ssa. This should look like a gaussian if things go right
TH1F * ssa2 = new TH1F("ssa2 distribution","ssa2 distribution",100,-10,10);
// Same as above but for ssa2
TH1F * overlap = new TH1F("overlap distribution","overlap distribution",100,-10,10);
// same as above but for the overlap region
TH1F * overall = new TH1F("overall distribution","overall distribution",100,-10,10);
// same as above but for the entire 240 x 240 difference histogram
// getting rid of stat boxes
output->SetBit(TH1::kNoStats);
ssa1->SetBit(TH1::kNoStats);
ssa2->SetBit(TH1::kNoStats);
overlap->SetBit(TH1::kNoStats);
overall->SetBit(TH1::kNoStats);
// expected->SetBit(TH1::kNoStats);
difference->SetBit(TH2::kNoStats);
//initializing to read the input csv
std::string line;
int comma;
std::string temp;
int hits;
int x = 0;
int y = 0;
std::vector diagonals;
// below fills in output and a vector of diagonals for use in expected
for(int j = 1; j <= 240; j++)
{
getline(in, line);
for(int i = 1; i <=240; i++)
{
comma = line.find(",");
temp = line.substr(0,comma);
// std::cout << temp << ",";
hits = std::stoi(temp);
if(i == j)
{
diagonals.push_back(hits);
}
if(hits!=0)
{
output->Fill(j,i,hits);
}
line.erase(0,comma+1);
}
// std::cout << std::endl;
}
for(int i = 0; i <240; i++) // debug purposes
{
// std::cout <GetBinContent(120,i) << ",";
}
std::cout << std::endl;
// below fills in expected, difference and outlier
for(int i = 1; i <=240; i ++)
{
for(int j = 1; j<=240;j++)
{
expected->Fill(i,j,diagonals[i-1]*diagonals[j-1]/20000000); // 20,000,000 since 5,000,000 events, and each diagonal is counted twice
if(abs(i-j) > 1)
{
float diff = output->GetBinContent(i,j)-expected->GetBinContent(i,j);
float dev = sqrt(output->GetBinContent(i,j));
difference->Fill(i,j,diff/dev);
if(i == 40 && j == 45) // testing purposes, see if things are being calculated properly or why a position has weird behavior
{
std::cout << "diff:" << diff << std::endl;
std::cout << "dev:" << dev << std::endl;
std::cout << "expected:" << expected->GetBinContent(i,j)<< std::endl;
std::cout << "observed:" << output->GetBinContent(i,j) << std::endl;
std::cout << "xvalue:" << diagonals[i-1]< 10) // 10 stdev away, outside the range of the 1d distribution histograms
{
strange->Fill(i,j);
}
overall->Fill(difference->GetBinContent(i,j)); // filling overall
// std::cout << difference->GetBinContent(i,j) << ",";
}
}
// std::cout << std::endl;
}
// below fills the 1d distribution histograms
for(int i = 1; i <=120; i++)
{
for(int j = 1; j <=120; j++)
{
ssa1->Fill(difference->GetBinContent(i,j));
overlap->Fill(difference->GetBinContent(i,j+120));
ssa2->Fill(difference->GetBinContent(i+120,j+120));
}
}
// drawing all the relevant pieces
difference->Draw("COLZ");
TCanvas *c2 = new TCanvas("pad2","pad2",150,10,990,660);
ssa1->Draw("HIST");
TCanvas *c3 = new TCanvas("pad3","pad3",150,10,990,660);
ssa2->Draw("HIST");
TCanvas *c4 = new TCanvas("pad4","pad4",150,10,990,660);
overlap->Draw("HIST");
TCanvas *c5 = new TCanvas("pad5","pad5",150,10,990,660);
overall->Draw("HIST");
TCanvas *c6 = new TCanvas("pad6","pad6",150,10,990,660);
strange->Draw("COLZ");
}