Skip to content
Snippets Groups Projects
Commit e7e0a4fa authored by Theodore Efremov's avatar Theodore Efremov :hibiscus:
Browse files

[AlPhaPha] Deleted redundant macro

parent 49646eb7
No related branches found
No related tags found
No related merge requests found
#include "TICPhysics.h"
#include <TChain.h>
#include <regex>
TChain* chain;
/////////////////////////////////////////////////////////////////
void FillHistoCharge( const char *Path){
chain = new TChain("PhysicsTree");
chain->Add(Path);
TICPhysics* IC = new TICPhysics() ;
chain->SetBranchStatus("IC", true);
chain->SetBranchAddress("IC", &IC);
double FF_IC_X;
chain->SetBranchStatus("FF_IC_X", true);
chain->SetBranchAddress("FF_IC_X", &FF_IC_X);
TH1F *hChargeAll = new TH1F("hcharge_all","hcharge_all",2000,20,80);
int Nentries = chain->GetEntries();
auto start = std::chrono::high_resolution_clock::now();
for (int e = 0; e < Nentries; e++) {
if (e % 100000 == 0 && e > 0 ) {
auto now = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = now - start;
double avgTimePerIteration = elapsed.count() / e;
double timeLeft = avgTimePerIteration * (Nentries - e);
std::cout << "********** Estimated time left: " << int(timeLeft) << " seconds **********" << "\r" << flush;
}
chain->GetEntry(e);
if(FF_IC_X>-530 ){
hChargeAll->Fill(IC->Chio_Z);
}
}
std::string strPath(Path);
std::regex pattern("Run(\\d{3})"); // Regex to match Run followed by exactly 3 digits
std::smatch matches;
if (std::regex_search(strPath, matches, pattern)) {
std::string number = matches[1]; // Extract the number
std::string fileName = "Output/histo_charge_" + number + ".root"; // Generate the filename
// Create the TFile
TFile* file = new TFile(fileName.c_str(), "RECREATE");
if (file->IsOpen()) {
std::cout << "File created: " << fileName << std::endl;
} else {
std::cout << "Failed to create the file." << std::endl;
}
hChargeAll->Write();
// Don't forget to close the file
file->Close();
} else {
std::cout << "Pattern not found!" << std::endl;
}
}
#include <TCanvas.h>
#include <TF1.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TH1.h>
using namespace std;
void Fit(TH1F* hcharge, string Energy);
TFile* ofile;
int rebin=1;
int RunNumber = 247 ;
//////////////////////////////////////////////////////
void FitCharge(){
//TFile* ifile = new TFile(Form("Output/histo_charge_%d.root",RunNumber),"read");
TFile* ifile = new TFile(Form("Output/histo_charge_merged.root"),"read");
ofile = new TFile("Output/histo_charge_fitted.root","recreate");
TH1F* hall = (TH1F*)ifile->FindObjectAny("hcharge_all");
Fit(hall,"AllEnergy");
ofile->Close();
}
//////////////////////////////////////////////////////
void Fit(TH1F* hcharge,string Energy){
TGraphErrors* gsigma = new TGraphErrors();
TGraphErrors* gsigmaZ = new TGraphErrors();
hcharge->Rebin(rebin);
hcharge->Draw();
int Zmin = 31;
int Zmax = 66;
int Z = Zmin;
int NumberOfZ = Zmax - Zmin;
Double_t para[3*NumberOfZ];
vector<TF1*> g(NumberOfZ);
vector<double> Integral(NumberOfZ);
vector<double> Integral_err(NumberOfZ);
double total_integral = 0;
vector<double> Yield(NumberOfZ);
for(int i=0; i<NumberOfZ; i++){
g[i] = new TF1(Form("g%i",i),"gaus",Z-0.3,Z+0.6);
if(i==0){
g[i]->SetParameter(0,100);
g[i]->SetParameter(1,Zmin);
g[i]->SetParameter(2,0.3);
}
else if(i>0){
g[i]->SetParameter(0,g[i-1]->GetParameter(0));
g[i]->SetParameter(1,Z);
g[i]->SetParameter(2,g[i-1]->GetParameter(2));
}
g[i]->SetParLimits(1,Z-0.2,Z+0.2);
g[i]->SetParLimits(2,0.22,0.4);
hcharge->Fit(g[i],"qBR");
//g[i]->Draw("lsame");
g[i]->GetParameters(&para[3*i]);
Z++;
}
TString total_func = "gaus(0)";
for(int i=1; i<NumberOfZ; i++){
TString gaus_func = Form("gaus(%i)",3*i);
total_func += "+" + gaus_func;
}
TF1* total = new TF1("total",total_func,Zmin-0.3, Zmin+NumberOfZ -0.7);
total->SetParameters(para);
Z = Zmin;
for(int i=0; i<NumberOfZ; i++){
total->SetParameter(3*i,g[i]->GetParameter(0));
total->SetParameter(3*i+1,g[i]->GetParameter(1));
total->SetParameter(3*i+2,g[i]->GetParameter(2));
total->SetParLimits(3*i+1,Z-0.2,Z+0.2);
total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.5,g[i]->GetParameter(2)*1.3);
Z++;
}
total->SetLineColor(4);
total->SetNpx(1000);
hcharge->Fit(total,"RBq");
Z = Zmin;
double MeanSigmaZ = 0 ;
//TF1* one_gauss = new TF1("one_gaus","gaus");
for(int i=0; i<NumberOfZ; i++){
double Amplitude = total->GetParameter(3*i);
double mean = total->GetParameter(3*i+1);
double sigma = total->GetParameter(3*i+2);
double Amplitude_err = total->GetParError(3*i);
double mean_err = total->GetParError(3*i+1);
double sigma_err = total->GetParError(3*i+2);
/*double Amplitude = g[i]->GetParameter(0);
double mean = g[i]->GetParameter(1);
double sigma = g[i]->GetParameter(2);
double Amplitude_err = g[i]->GetParError(0);
double mean_err = g[i]->GetParError(1);
double sigma_err = g[i]->GetParError(2);*/
//one_gauss->SetParameter(0,Amplitude);
//one_gauss->SetParameter(1,mean);
//one_gauss->SetParameter(2,sigma);
//one_gauss->SetParError(0,Amplitude_err);
//one_gauss->SetParError(1,mean_err);
//one_gauss->SetParError(2,sigma_err);
cout << "************************" << endl;
cout << "Z= " << Z << endl;
//cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
//cout << "Mean= " << mean << endl;
cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;
cout << "Sigma/Z = " << 2.35*sigma /Z*100. << endl;
//if (2.35*sigma/Z*100 > 1 && 2.35*sigma/Z*100 < 3){MeanSigmaZ += 2.35*sigma/Z*100;}
if (i==0 || i==NumberOfZ){continue;}
MeanSigmaZ += 2.35*sigma/Z*100;
gsigma->SetPoint(i-1,Z,sigma);
gsigma->SetPointError(i-1,0,sigma_err);
gsigmaZ->SetPoint(i-1,Z,2.35*sigma/Z*100);
Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
//Integral_err[i] = sqrt(Integral[i]);
total_integral += Integral[i];
Z++;
}
MeanSigmaZ = MeanSigmaZ/NumberOfZ;
cout << "********************" << endl << "Mean Sigma/Z " << MeanSigmaZ << endl;
double histo_integral = hcharge->Integral();
hcharge->Write();
total->Write();
string output_filename = "Charge_Yield_" + Energy + ".dat";
ofstream ofile;
ofile.open(output_filename.c_str());
// Yield calculation //
Z = Zmin;
for(int i=0; i<NumberOfZ; i++){
//Yield[i] = Integral[i]/total_integral*200;
Yield[i] = Integral[i]/histo_integral*200;
//double yield_err = Integral_err[i]/total_integral*200;
double yield_err = Integral_err[i]/histo_integral*200;
ofile << Z << " " << Yield[i] << " " << yield_err << endl;
Z++;
}
ofile.close();
gsigma->SetMarkerStyle(8);
gsigma->Write();
TCanvas* csig = new TCanvas("csig","csig",800,800);
csig->cd();
gsigma->Draw("ap");
TCanvas* csigZ = new TCanvas("csigZ","csigZ",800,800);
gsigmaZ->SetMarkerStyle(8);
gsigmaZ->GetYaxis()->SetTitle("2.35*#sigma/Z*100");
gsigmaZ->GetXaxis()->SetTitle("Z");
csigZ->cd();
gsigmaZ->Draw("ap");
total->Draw("lsame");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment