Commit 2dfcb4ec authored by Keerthi Nakkalil's avatar Keerthi Nakkalil
Browse files

A fine_toacut is implemented for fitting timewalk function. The code is still...

A fine_toacut is implemented for fitting timewalk function. The code is still not working. All the pixels somehow give 0 points for fitting.
parent e878adcd
void TreeEnergyWrite(){
int Ecol1;
double Ecol2,Ecol3,Ecol4;
string line_E;
TFile* Efile = new TFile("Energy_table.root","recreate");
TTree* ETree = new TTree("ETree","Energy tree");
ETree->Branch("Ecol1",&Ecol1,"Ecol1/I"); // ECol1 = Energy values = {200,300.....14000}
ETree->Branch("Ecol4",&Ecol4,"Ecol4/D"); // Voltage values
std::ifstream Energy_file("Energy_table.dat");
if(Energy_file.is_open()){
while(getline(Energy_file,line_E)){
istringstream ssE(line_E);
ssE >> Ecol1 >> Ecol2 >> Ecol3 >> Ecol4;
ETree->Fill();
}
}
ETree->Write();
Efile->Close();
}
void TreeTPWrite(){
int energy, col, row, pc;
double voltage,tot, tot_var, toa, toa_var, time_ref,rough_toacut=0,fine_toacut=0,mean=0,RMS=0;
string line_TP;
// std::ofstream skippixels;
// skippixels.open("toa_skip_pixels_1000_1190.txt",std::ofstream::app);
TFile* Efile = new TFile("Energy_table.root");
TTree* ETree = (TTree*)Efile->Get("ETree");
ETree->SetBranchAddress("Ecol1",&energy);
ETree->SetBranchAddress("Ecol4",&voltage);
TFile* TPfile = new TFile("toa_new.root","recreate");
TTree* TPtree = new TTree("TPtree","TP tree");
TPtree->Branch("col",&col,"col/I");
TPtree->Branch("row",&row,"row/I");
TPtree->Branch("pc",&pc,"pc/I");
TPtree->Branch("tot",&tot,"tot/D");
TPtree->Branch("tot_var",&tot_var,"tot_var/D");
TPtree->Branch("toa",&toa,"toa/D");
TPtree->Branch("toa_var",&toa_var,"toa_var/D");
TPtree->Branch("time_ref",&time_ref,"time_ref/D");
auto voltageBranch = TPtree->Branch("voltage",&voltage,"voltage/D");// make a new branch and connect it to the same variable as the branch address of the energy tree
auto Branch_coarse_toacut = TPtree->Branch("rough_toacut",&rough_toacut,"rough_toacut/D");
auto Branch_fine_toacut = TPtree->Branch("fine_toacut",&fine_toacut,"fine_toacut/D");
auto Branch_mean = TPtree->Branch("mean",&mean,"mean/D");
auto Branch_RMS = TPtree->Branch("RMS",&RMS,"RMS/D");
for(int i = 0; i < ETree->GetEntries(); i++){
ETree->GetEntry(i);
// std:: cout << energy << " " << voltage << std::endl;
std::string TPfilename = "tp_dat_" + std::to_string(energy) + ".dat";
//std:: cout << " The TPfilenames are : " << TPfilename.c_str() <<std::endl;
std::ifstream TPdatafile(TPfilename.c_str());
voltageBranch->Fill();
if(TPdatafile.is_open()){
std::cout << " The data file " << TPfilename.c_str() << " is opened ..." <<std::endl;
// reading through individial tp_dat_ files
while(getline(TPdatafile,line_TP)){
if(line_TP.rfind("#", 0) == 0) {
continue;
}
istringstream ssTP(line_TP);
ssTP >> col >> row >> pc >> tot >> tot_var >> toa >> toa_var >> time_ref;
TPtree->Fill();
} //While loop getline
}//if statement TPdatafile open
TPtree->Draw("toa",Form("voltage==(%f)",voltage),"goff");
mean = TMath::Mean(TPtree->GetSelectedRows(), TPtree->GetV1());
Branch_mean->Fill();
RMS = TMath::RMS(TPtree->GetSelectedRows(), TPtree->GetV1());
Branch_RMS->Fill();
// std::cout <<" mean " << mean << " RMS " << RMS <<std::endl;
/* for (int isigma=0; isigma < 15; isigma++){
std::cout << " cut at " << isigma << " times the RMS leaves : "
<< 100.0*TPtree->GetEntries(Form("toa<(%f) && voltage==(%f)",mean+isigma*RMS,voltage))/TPtree->GetEntries(Form("voltage==(%f)",voltage))
<< " percentage of entries " << std::endl;
}*/
rough_toacut = mean + 6*RMS;
// cout << " toa : " << toa << " " << "rough toa cut : " << rough_toacut <<std::endl;
Branch_coarse_toacut->Fill();
double entries = TPtree->GetEntries(Form("voltage==(%f) && pc>250",voltage));
double tot_entries = TPtree->GetEntries(Form("voltage==(%f)",voltage)) ;
double ratio = (entries)/(tot_entries);
// cout << " entries " << entries << " " << " tot_entries " << tot_entries << " ratio " << ratio << std::endl;
if(ratio < 0.019)
{
std::cout << "skipping the energy value" << std::endl;
}
else {
// c->SetLogy();
gStyle->SetOptFit(1111);
TF1 *fgaus = new TF1("fgaus", "gaus") ;
TPtree->Fit("fgaus", "toa>>hfit(1000)",Form("voltage==(%f)&& toa<(%f)&& toa>0 && pc>250",voltage,rough_toacut));
TPtree->Draw(Form("toa>>htoa_%d", i),Form("voltage==(%f)&& toa<(%f)&& toa>0 && pc>250",voltage,rough_toacut),"same goff");//coarse selection
// gPad->SaveAs(Form("toa>>htoa_coarsecut_%d.png", i));
auto gauss_mean = fgaus->GetParameter(1);
auto gauss_sigma = fgaus->GetParameter(2);
// cout << " gauss mean : " << gauss_mean << " " << "gauss sigma : " << gauss_sigma << std::endl;
fine_toacut = gauss_mean + 6*gauss_sigma;
Branch_fine_toacut->Fill();
// std::cout << " coarse toa cut : " << rough_toacut << " " << "fine toa cut : " << fine_toacut << std::endl;
}
}//For loop GetEntries
//Save only the new version of the tree
TPtree->Write("",TObject::kOverwrite);
TPfile->Write();//KN
TPfile->Close();
}
void timewalkfit(){
int col, row, pc ;
double tot, tot_var, toa, toa_var,voltage,fine_toacut,rough_toacut;
TFile* TPfile = new TFile("toa.root");
TTree* TPtree = (TTree*)TPfile->Get("TPtree");
TPtree->SetBranchAddress("col",&col);
TPtree->SetBranchAddress("row",&row);
TPtree->SetBranchAddress("pc",&pc);
TPtree->SetBranchAddress("tot",&tot);
TPtree->SetBranchAddress("tot_var",&tot_var);
TPtree->SetBranchAddress("toa",&toa);
TPtree->SetBranchAddress("toa_var",&toa_var);
TPtree->SetBranchAddress("voltage",&voltage);
TPtree->SetBranchAddress("fine_toacut",&fine_toacut);
TPtree->SetBranchAddress("rough_toacut",&rough_toacut);
for(int irow = 0; irow < 256; irow++){
for (int icol = 0; icol <256; icol++){
std::cout << icol << " " << irow << " " << std::endl;
auto N = TPtree->Draw("voltage:toa:sqrt(toa_var)",Form("col==(%i) && row==(%i) && toa>0 && toa<(%f)", icol, irow, fine_toacut),"goff");
cout << N <<std::endl;
std::ofstream fitpars;
if(N==0)
{
std::cout << "Dead pixel or no enough points for fitting. Skipping pixel " << icol << "_" << irow << std::endl;
std::ofstream skippixels;
skippixels.open("toa_skip_pixels_new_1000_1190.txt",std::ofstream::app);
skippixels << "Skipping the pixel " << icol << "_" << irow <<std::endl;
fitpars.open("toa_fit_parameters_new_1000_1190.txt",std::ofstream::app);
fitpars << icol << " " << irow << " " << "0.0"
<< " " << "0.0" << " "
<< "0.0" << " " << "0.0" << " " << "0.0" << std::endl;
}
else{
TGraphErrors *graph = new TGraphErrors(N,TPtree->GetV1(),TPtree->GetV2(),0,TPtree->GetV3());
// std::cout << icol << " " << irow << std::endl; //Correct values printed here..but text file still prints 255_255
auto fitmin = TMath::MinElement(graph->GetN(),graph->GetX());
// std::cout << "fitmin" << fitmin << std::endl;
auto fitmax = TMath::MaxElement(graph->GetN(),graph->GetX());
// std::cout << "fitmax " << fitmax << std::endl;
TF1 *fitfuncTOA = new TF1("fitfuncTOA","(([2]/(x-[0]))+ [1])",fitmin,fitmax);
fitfuncTOA->SetParNames("asymtote","offset","curvature");
int iymax=TMath::LocMax(graph->GetN(),graph->GetY());
// cout << "iymax " << iymax << std::endl;
auto ttoa = graph->GetPointX(iymax);
// cout << "ttoa" << " " << ttoa << std::endl;;
fitfuncTOA->SetParameter(0,ttoa);
auto xmax = TMath::MaxElement(graph->GetN(),graph->GetX());
auto dtoa = graph->Eval(xmax);
// cout << "dtoa" << " " << dtoa <<std::endl;
fitfuncTOA->SetParameter(1,dtoa);
auto testx = 0.5*fitmax;
auto ctoa_y = graph->Eval(testx);
auto ctoa = (ctoa_y - dtoa)*(testx - ttoa);
// cout << "ctoa" << " " << ctoa <<std::endl;
fitfuncTOA->SetParameter(2,ctoa);
// gStyle->SetOptFit(1111);
// fitfuncTOA->SetLineColor(kRed);
graph->Fit(fitfuncTOA,"R");
/* graph->GetXaxis()->SetTitle("Volatge[mV]");
graph->GetYaxis()->SetTitle ("time difference[ns]");
graph->SetTitle(" ");
graph->SetMarkerStyle(7);*/
// graph->Draw("AP");
fitpars.open("toa_fit_parameters_new_1000_1190.txt",std::ofstream::app);
fitpars << icol << " " << irow << " " << fitfuncTOA->GetParameter(2)
<< " " << fitfuncTOA->GetParameter(0) << " " << fitfuncTOA->GetParameter(1) << " " << (fitfuncTOA->GetChisquare()/fitfuncTOA->GetNDF()) << std::endl; // chi2/ndf included for the time-being!
}
}
}
TPfile->Close();
}
void toa_caliberation(){
TreeEnergyWrite();
TreeTPWrite();
timewalkfit();
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment