Here is the whole code:
#include <cmath>
#include <vector>
#include <fstream>
{
// allpackets-SPBEUSO-ACQUISITION-20160930-050017-001.001--bg.root
// allpackets-SPBEUSO-ACQUISITION-20160930-051250-001.001--10Hz_50mus_1500mV.root
//---- define structures -----------//
struct Squares
{
double Radius[100][100]; // matrix to represent the distance from the center of the square to the center of each 1cm*1cm squares
double Incident_angle[100][100]; //matrix of incident angles of light to each 1cm by 1cm squares
double Distance[100][100]; // matrix of the dstances from LED to each 1cm by 1cm squares
double Incident_light[100][100]; // matrix to store the amount of incident light to each 1cm by 1cm square
double Reflected_light[100][100]; // matrix to store the reflected light from each 1cm by 1 cm square
};
/////////////////////////////////////////////////////////////////////////////////////
struct EUSO
{
Int_t count_entry[128] ={0};
vector<Int_t> ENTRY;
};
//----- open the file and read the tree -----------//
TFile f1("allpackets-SPBEUSO-ACQUISITION-20160930-051250-001.001--10Hz_50mus_1500mV.root");
TTree* tevent; f1.GetObject("tevent",tevent);
UChar_t photon_count_utah[1][1][48][48];
tevent->SetBranchAddress("photon_count_data",photon_count_utah);
//-----------define some stuffs ---------//
TCanvas* c2 = new TCanvas("c1","",700,700);
//c2->Divide(1,4);
c2->SetGrid(1,1);
gStyle->SetPalette(1);
gStyle->SetOptStat(0);
// ---- read the first branch ----//
// ---- read the first branch ----//
std:ifstream input;
input.open("relative.txt");
if(!input) {cout<<"wrong"<<endl;}
const int N=48;
const int Entries = 25600;
EUSO back_entry[200];
EUSO LED_entry[200];
EUSO en_back[200];
EUSO LED[200];
double Average_Count[Entries];
double Entry_x[Entries];
double Background[N][N];
double Back_sum[N][N]={0};
double Av_count[N][N] ={0};
double LED_BACK[N][N] = {0};
double Calib[N][N]={0};
double Re_LED[N][N];
double count = 0;
int count1 = 0;
int s=0;
//-----------Define 2D histograms-------//
TH2F* Back_ground = new TH2F("Average Background count", "",48,0,48,48,0,48);
Back_ground->SetTitle("Average Background count; Pixel x; Pixel y");
TH2F* LED_average = new TH2F("Average LED count", "",48,0,48,48,0,48);
LED_average->SetTitle("Average LED count; Pixel x; Pixel y");
LED_average->GetXaxis()->SetNdivisions(-706); LED_average->GetYaxis()->SetNdivisions(-706);
TH2F* LED_Back = new TH2F("Average LED-Background count", "",48,0,48,48,0,48);
LED_Back->SetTitle("Background Substracted Signal on (Utah); Pixel x; Pixel y");
LED_Back->GetXaxis()->SetNdivisions(-706); LED_Back->GetYaxis()->SetNdivisions(-706);
TH2F* Illumination = new TH2F("Illumination of the PDM", "",48,0,48,48,0,48);
Illumination->SetTitle("Illumination of the PDM; Pixel x; Pixel y");
Illumination->GetXaxis()->SetNdivisions(-706); Illumination->GetYaxis()->SetNdivisions(-706);
TH2F* CALIB = new TH2F("Illumination divided on the count", "",48,0,48,48,0,48);
CALIB->SetTitle("Illumination Divided by the average count; Pixel x; Pixel y");
CALIB->GetXaxis()->SetNdivisions(-706); CALIB->GetYaxis()->SetNdivisions(-706);
TH2F* re_led = new TH2F("Illumination divided on the count", "",48,0,48,48,0,48);
re_led->SetTitle("Relative LED Signal; Pixel x; Pixel y");
re_led->GetXaxis()->SetNdivisions(-706); re_led->GetYaxis()->SetNdivisions(-706);
TH1F* Poisson = new TH1F("# Count","",48,0,48);
Poisson->SetTitle("Background Average Number of Count; Bin number; number of count");
Double_t av_avB=0;
Double_t av_avL=0;
Double_t cnt_b=0;
Double_t cnt_L=0;
Double_t sigma_B =0;
Double_t sigma_L=0;
for(int r =0; r<200;r++)
{
for( s=(r*128);s<(128*(r+1));++s)
{
double sum = 0;
tevent->GetEntry(s);
for(int i =0; i<N;++i)
{
for(int j =0;j<N;++j)
{
Background[i][j] = static_cast<double>(photon_count_utah[0][0][i][j]);
sum+=Background[i][j];
}
}
Average_Count[s]=sum/(N*N); //cout<<s<<"\t"<<Average_Count[s]<<endl;
if(Average_Count[s] <2.9) {sigma_B +=(Average_Count[s]-2.75176)*(Average_Count[s]-2.75176)/21875; cnt_b++;}
if(Average_Count[s] >3) {sigma_L +=(Average_Count[s]-3.13761)*(Average_Count[s]-3.13761)/3613; cnt_L++;}
Entry_x[s] = s+1;
if(Average_Count[s] <2.8438342 && Average_Count[s] > 2.6596858 ) // it was zero
{
back_entry[r].count_entry[(s-(r*128))] = s;
}
else if(Average_Count[s] >3.036382 && Average_Count[s] < 3.388398)
{
LED_entry[r].count_entry[(s-(r*128))] = s;
}
}
}
cout<<av_avB/cnt_b<<"\t"<<cnt_b<<"\t"<<TMath::Sqrt(sigma_B)<<endl;
cout<<av_avL/cnt_L<<"\t"<<cnt_L<<"\t"<<TMath::Sqrt(sigma_L)<<endl;
//TGraph* PDM_count = new TGraph(Entries,Entry_x,Average_Count);
//PDM_count->SetTitle("PDM Mean Count Light Curve (Utah);GTU[2.5#mus] ; Measured Signal [PDM counts/pixel/GTU]");
//PDM_count->Draw();
//cout<<count<<endl;
/////////////////////////////////////////////
for(int i = 0;i<200;++i)
{
for(int j = 0;j<128;++j)
{
if(back_entry[i].count_entry[j] != 0)
(en_back[i].ENTRY).push_back(back_entry[i].count_entry[j]);
if(LED_entry[i].count_entry[j] !=0)
(LED[i].ENTRY).push_back(LED_entry[i].count_entry[j]);
}
}
//-----average background pixels -------//
for(int i = 0; i <200;++i)
{
Int_t z = round((en_back[i].ENTRY).size()*0.1);
for(int j=z;j<(((en_back[i].ENTRY).size())-z);++j)
{
tevent->GetEntry(en_back[i].ENTRY[j]);
for(int l =0; l<N;++l)
{
for(int k =0;k<N;++k)
{
Back_sum[l][k]+= static_cast<double>(photon_count_utah[0][0][l][k]);
}
}
count++;
}
}
//------------------Fill The Background--------------------------------//
for(int i =0; i<N;++i)
{
for(int j =0;j<N;++j)
{
Back_ground->SetBinContent(i+1,j+1,Back_sum[i][j]/count);
//Poisson->Fill(Back_sum[i][j]/count);
}
}
//- Fill the LED On -----//
for(int i = 0; i <200;++i)
{
int z = round((LED[i].ENTRY).size()*0.1);
for(int j=z;j<(((LED[i].ENTRY).size())-z);++j)
{
tevent->GetEntry(LED[i].ENTRY[j]);
for(int l =0; l<N;++l)
{
for(int k =0;k<N;++k)
{
Av_count[l][k]+= static_cast<double>(photon_count_utah[0][0][l][k]);
}
}
count1++;
}
}
//--------------Fill LED ON -------------//
for(int i =0; i<N;++i)
{
for(int j =0;j<N;++j)
{
LED_average->SetBinContent(i+1,j+1,Av_count[i][j]/count1);
//Poisson->Fill(Back_sum[i][j]/count1);
}
}
////////////////////////////////LED-BKG///////////
for(int i =0;i <N;++i)
{
for(int j =0;j<N;++j)
{
LED_BACK[i][j] = (Av_count[i][j]/count1) - (Back_sum[i][j]/count);
if(LED_BACK[i][j]<0) LED_BACK[i][j] =0;
}
}
//--------------- Fill LED-Background-------//
for(int i =0; i<N;++i)
{
for(int j =0;j<N;++j)
{
LED_Back->SetBinContent(i+1,j+1,LED_BACK[i][j]);
//Poisson->Fill(Back_sum[i][j]/count1);
}
}
//-----storing relative calibration in a vector -----//
const int size = 2304;
Double_t dummy = 0;
double rel[size] ={0};
for(int i=0;i<size;++i)
{
input>>dummy;
out<<dummy<<endl;
}
/////////////////////////////////////////////////////
int counter =0;
for(int i =0; i<N;++i)
{
for(int j =0;j<N;++j)
{
Re_LED[i][j]=LED_BACK[i][j]*rel[counter];
re_led->SetBinContent(i+1,j+1,Re_LED[i][j]);;
counter++;
}
}
//-----sum Relative led signal-------//
double sum_led=0;
for(int i=22;i<35;++i)
for(int j=19;j<34;++j)
{
sum_led+=Re_LED[i][j];
}
cout<<sum_led<<endl;
//Back_ground->Draw("colz");
//LED_average->Draw("colz");
//LED_Back->Draw("colz");
re_led->Draw("colz");
TBox* box = new TBox(0.1,0.1, 0.8,0.8);
box->SetFillStyle(0);
box->SetLineColor(2);
box->SetLineWidth(2);
box->Draw("I");
//////////
}