I want to reconstruct the energy from the gaussian fit but instead i am getting the error.
terminate called after throwing an instance of 'std::out_of_range'
what(): vector::_M_range_check: __n (which is 53861) >= this->size() (which is 53861)
Can someone please identify the mistake.
Below is the problematic part of the code
TF1 *fun[points];
for (int i = 0; i < points; i++) {
for (int k = 0; k < N123[i].size(); k++) {
double alpha = p0A + p1A*N123[i][k] + p2A*N123[i][k]*N123[i][k];
double beta = p0B + p1B*N123[i][k] + p2B*N123[i][k]*N123[i][k];
double gamma = p0C + p1C*N123[i][k] + p2C*N123[i][k]*N123[i][k];
Erec[i].push_back(alpha * N1[i][k] + beta * N2[i][k] + gamma * N3[i][k]);
double sum = pow(N1[i][k]*AlphaErrors[i], 2) + pow(N2[i][k]*BetaErrors[i], 2) + pow(N3[i][k]*GammaErrors[i], 2) +
N1[i][k]*Alpha[i]*Alpha[i] + N2[i][k]*Beta[i]*Beta[i] + N3[i][k]*Gamma[i]*Gamma[i];
Error_Erec[i].push_back(pow(sum, 0.5));
energyrec[i]->Fill(Erec[i].back());
}
TString fname = Form("f_%d_GeV", i);
fun[i] = new TF1(fname, "gaus", 0, 100);
fun[i]->SetParameters(1, energyrec[i]->GetMean(), energyrec[i]->GetStdDev());
}
TString root_fileName = "SDHCAL_Erec_3cm_t3_12MIP.root";
TFile *fout = new TFile(root_fileName, "RECREATE");
fout->cd();
TTree *outTree = new TTree("SDHCAL_Erec", "SDHCAL_Erec");
Int_t k = 0; Double_t en[points] = {0.}, ener[points] = {0.};
outTree->Branch("evt", &k);
for(int i = 0; i < points; i++){
int j = i + 1;
TString branchName = Form("%dGeV", j);
TString branchName2 = Form("%dGeV_err", j);
outTree->Branch(branchName, &en[i]);
outTree->Branch(branchName2, &ener[i]);
}
for (int k = 0; k < Erec[0].size(); k++) {
for (int i = 0; i < points; i++) {
en[i] = Erec[i].at(k);
ener[i] = Error_Erec[i].at(k);
}
outTree->Fill();
}
fout->cd();
fout->Write();
fout->Close();
Double_t meanErec[points], sigmaErec[points], err_sigmaErec[points], err_meanErec[points];
Double_t resolution[points], err_resolution[points];
TCanvas* c5 = new TCanvas();
c5->SetTickx();
c5->SetTicky();
for (int i = 0; i < points; i++) {
TFitResultPtr r = energyrec[i]->Fit(fun[i], "SQ");
meanErec[i] = r->Parameter(1);
err_meanErec[i] = r->ParError(1);
sigmaErec[i] = r->Parameter(2);
err_sigmaErec[i] = r->ParError(2);
resolution[i] = sigmaErec[i] / meanErec[i];
err_resolution[i] = resolution[i] * sqrt(
pow(err_sigmaErec[i] / sigmaErec[i], 2) +
pow(err_meanErec[i] / meanErec[i], 2));
}
TCanvas* c6 = new TCanvas("c6", "Energy Resolution", 800, 600);
c6->SetTickx();
c6->SetTicky();
TGraphErrors* res = new TGraphErrors(points, energy, resolution, errorz, err_resolution);
res->SetMarkerStyle(22);
res->SetMarkerColor(kRed + 1);
res->SetTitle("");
res->GetXaxis()->SetTitle("E_{#pi} [GeV]");
res->GetYaxis()->SetTitle("#sigma(E)/E");
res->Draw("AP");
TF1* fit_resolution = new TF1("fit_resolution", "[0]/sqrt(x) + [1]", 1., 80.);
fit_resolution->SetParName(0, "a");
fit_resolution->SetParName(1, "b");
res->Fit(fit_resolution, "R");