//
// pulseFit.C
// Obtaining the integral of the pulse shape (between the roots)
// representing the amount of collected charge
//
// Created by Christian Elsasser on 31.05.13.
// University of Zurich, elsasser@cern.ch
// Copyright only for commercial use
//
// Compile with:
// g++ -o pulseFit pulseFit.C `root-config --glibs --cflags` -lRooFit -lRooFitCore -lEG
//
// General include
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <stdlib.h>
#include <limits>
#include <string>
#include <time.h>
// ROOT include
#include "../../include/incBasicROOT.h"
#include "../../include/incFuncROOT.h"
#include "../../include/incDrawROOT.h"
#include "../../include/incGraphROOT.h"
#include "../../include/incIOROOT.h"
#include "../../include/basicROOTTools.h"
#include "../../include/lhcbstyle.C"
#include "deplTool.h"
// -d detector (TT/IT)
// -l layer (TTaX,TTaU,TTbV,TTbX, T[1-3]{X1,U,V,X2})
// -s read-out sector number
// -f fill number
// -c Voltage step [0-10]
// -r consider only hit in a distance closer than r mm from the beam axis
// -v number of summed up strips
// -ns Do not save the value in the root file
// -ss Perform a temporary save
int pulseFit(int argc, char* argv[]);
int main(int argc, char* argv[]){
TStyle* style = lhcbStyle();
int argc_copy = argc;
TApplication* theApp = new TApplication("Analysis", &argc, argv);
argc = theApp->Argc();
argv = theApp->Argv();
bool runBatch = false;
for (int i = 1; i<argc; i++) {
if (strcmp(argv[i],"-b")==0) {
runBatch = true;
gROOT->SetBatch();
}
}
int exit = pulseFit(argc,argv);
Printf("End");
if (!runBatch) {
theApp->Run(!runBatch);
}
return exit;
}
// Executable method
int pulseFit(int argc, char* argv[]){
//TVectorD v_res(11);
TVectorD v_res(12); //MIchele
gStyle->SetPadLeftMargin(0.20);
gStyle->SetTitleOffset(1.5,"Y");
gStyle->SetTitleOffset(1.0,"X");
const char* diskVar = std::getenv ("DISK");
int fill = 2083;
TString det("TT");
TString lay("TTaU");
int sector = 2634;
int caliStep = 0;
int val = 3;
bool notSave = false;
bool temSave = false;
double radius = -10.0;
for (int i=1; i<argc; i++) {
if (strcmp(argv[i],"-d")==0 && i+1<argc) {
det = argv[++i];
}else if (strcmp(argv[i],"-l")==0 && i+1<argc) {
lay = argv[++i];
}else if (strcmp(argv[i],"-s")==0 && i+1<argc) {
sector = std::atoi(argv[++i]);
}else if (strcmp(argv[i],"-f")==0 && i+1<argc) {
fill = std::atoi(argv[++i]);
}else if (strcmp(argv[i],"-v")==0 && i+1<argc) {
val = std::atoi(argv[++i]);
}else if (strcmp(argv[i],"-c")==0 && i+1<argc) {
caliStep = std::atoi(argv[++i]);
}else if (strcmp(argv[i],"-r")==0 && i+1<argc) {
radius = std::atof(argv[++i]);
}else if (strcmp(argv[i],"-ns")==0) {
notSave = true;
}else if (strcmp(argv[i],"-ss")==0) {
temSave = true;
}
}
int voltStep = caliStep;
double volt = 0.0;
// Vectors to store the ADC and time values obtained in the landau fit
std::vector<double> v_time_val;
std::vector<double> v_time_err;
std::vector<double> v_adc_val;
std::vector<double> v_adc_err;
double time = STTool::GetTime(caliStep, fill);
// Get the correct voltage value
if (det.EqualTo("TT")) {
volt = STTool::GetTTVoltage(caliStep*6);
}else{
volt = STTool::GetITVoltage(caliStep*6);
}
// Name of the vector where the result is stored
TString vn_result(Form("v_pulse%d_val%d",caliStep,val));
// Directory and file name for the plots
TString dn_graph(Form("%s/data/ST/graphics/%s/%s/%d/%d",
diskVar,det.Data(),lay.Data(),sector,fill));
TString fn_graph_p(Form("%s/pulse_%d_val%d",
dn_graph.Data(),caliStep,val));
// Dealing with the radius
if (radius<0.0) {
}else{
vn_result +=Form("_r%d",(int)radius);
fn_graph_p +=Form("_r%d",(int)radius);
}
fn_graph_p += ".pdf";
Printf("=========================================");
Info("pulseFit","Analyze: ");
Printf(" Sector: %d",sector);
if (radius>0.0) {
Printf(" Radius: %4.2f mm",radius);
}
Printf(" Fill: %d",fill);
Printf(" Detector: %s",det.Data());
Printf(" Layer: %s",lay.Data());
Printf(" Voltage: %4.2f V",volt);
// Directory and file name for the input data
TString dn_data(Form("%s/data/ST/Aging"
,diskVar));
TString fn_data(Form("%s/CCEScan.root",
dn_data.Data()));
TString fn_dummy(Form("%s/check",dn_data.Data()));
if (temSave) {
fn_data= TString(Form("%s/temp_%s_%s_%d_%d_r%d.root",
dn_data.Data(),det.Data(),lay.Data(),
sector,fill,(int)radius));
}
TFile* f_data = TFile::Open(fn_data.Data());
if (!f_data) {
Error("pulseFit","Data file does not exist!");
Printf("=========================================\n");
return EXIT_FAILURE;
}
// Loading the values from the Landau distribution fits and the corresponding time
Info("pulseFit","Loading data:");
TString dn_input(Form("%s/%s/%d/%d",
det.Data(),
lay.Data(),
sector,
fill));
for (int i=6*caliStep; i<6*caliStep+6; i++) {
TString vn_input(Form("%s/v_%d_val%d",dn_input.Data(),i,val));
if (radius>0.0) {
vn_input +=Form("_r%d",(int)radius);
}
TVectorD* v_input = (TVectorD*)f_data->Get(vn_input.Data());
if (!v_input) {
Warning("pulseFit","Vector %s not available",vn_input.Data());
continue;
}
Printf(" %8.2f ns: %8.4f+/-%8.4f",STTool::GetTime(i,fill),(*v_input)(0),(*v_input)(1));
if (i==6*caliStep) {
v_res(6) = (*v_input)(0);
v_res(7) = (*v_input)(1);
}
// Excluding steps bad data quality
if (i==0 && (fill==2797 || fill==3108)) {
Info("pulseFit","Bad calibration step! Skip it!");
continue;
}
v_adc_val.push_back((*v_input)(0));
v_adc_err.push_back((*v_input)(1));
v_time_val.push_back(STTool::GetTime(i,fill));
v_time_err.push_back(0.0);
}
if (v_adc_val.size()<3) {
Error("pulseFit","Too few (%d) data points",v_adc_val.size());
Printf("=========================================\n");
return EXIT_FAILURE;
}
// Transform STL vectors into ROOT vectors
TVectorD vr_adc_val = getROOTVector(v_adc_val);
TVectorD vr_adc_err = getROOTVector(v_adc_err);
TVectorD vr_time_val = getROOTVector(v_time_val);
TVectorD vr_time_err = getROOTVector(v_time_err);
TGraphErrors* ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err);
// Fit the pulse shape
TF1* fu_pulse = STTool::halfGauss;
ge_pulse->Fit(fu_pulse,"M0Q","0");
double chi2ndf = fu_pulse->GetChisquare()/fu_pulse->GetNDF();
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());
// Make PDG-like correction of the uncertainties
//Commented out by Michele - 8/31 2017
//if (chi2ndf>1.0) {
// Info("pulseFit","scaling");
// vr_adc_err*=TMath::Sqrt(chi2ndf);
//}
//Michele
TMultiGraph* mg_pulse = new TMultiGraph("mg_pulse","Pulse multigraph");
ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err);
ge_pulse->Fit(fu_pulse,"M0Q","0");
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());
//Evaluation of the systematic error on the integral of the pulse shape using a polinomial
//Define the interval were most likely the zeros should be
double low_limit = -50;
double high_limit = 50;
//POLY3
cout << "================ poly3 ====================" << endl<<endl;
TF1 *poly3 = new TF1("poly3","[0] + [1]*x +[2]*x**2 + [3]*x**3 ",-40.,40.);
poly3->SetLineColor(1);
ge_pulse->Fit("poly3","M");
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly3->GetChisquare(),poly3->GetNDF());
double x1_poly3 = poly3->GetX(0,low_limit,0); //GetX will return x* such that f(x*)=0 only if this exists, will (probably) return the argmin(f(x)) in the interval otherwis
double x2_poly3 = poly3->GetX(0,0, high_limit);
double z1_poly3 = -99;
double z2_poly3 = -99;
bool b_poly3 = false;
double integral_poly3_val = -999;
double integral_poly3_err = 0;
if( TMath::Abs(poly3->Eval(x1_poly3)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly3->Eval(x1_poly3)<<endl;}
else{
z1_poly3 = x1_poly3;
cout <<"Zero found! z1: "<< z1_poly3<<", y1: "<<poly3->Eval(x1_poly3)<<endl;
if( TMath::Abs(poly3->Eval(x2_poly3)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly3->Eval(x2_poly3)<<endl;}
else{
z2_poly3 = x2_poly3;
cout << "Zero found! z2: "<<z2_poly3<<", y2: "<<poly3->Eval(x2_poly3)<<endl;
b_poly3 = true;
}
}
//This has to be done before any other fit is performed
if(b_poly3 == true){
cout << "Poly3 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
integral_poly3_val = poly3->Integral(z1_poly3,z2_poly3)/20.;
integral_poly3_err = poly3->IntegralError(z1_poly3,z2_poly3)/20.;
cout << "The integral between the two zeros is "<< integral_poly3_val <<"+/-"<<integral_poly3_err << endl;
}
//POLY4
cout << "============ poly4 ==============" << endl<<endl;
TF1 *poly4 = new TF1("poly4","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 ",-40.,40.);
poly4->SetLineColor(4);
ge_pulse->Fit("poly4","M0Q","0");
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly4->GetChisquare(),poly4->GetNDF());
double x1_poly4 = poly4->GetX(0,low_limit,0);
double x2_poly4 = poly4->GetX(0,0, high_limit);
double z1_poly4 = -99;
double z2_poly4 = -99;
bool b_poly4 = false;
double integral_poly4_val = -999;
double integral_poly4_err = 0;
if( TMath::Abs(poly4->Eval(x1_poly4)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly4->Eval(x1_poly4)<<endl;}
else{
z1_poly4 = x1_poly4;
cout <<"Zero found! z1: "<< z1_poly4<<", y1: "<<poly4->Eval(x1_poly4)<<endl;
if( TMath::Abs(poly4->Eval(x2_poly4)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly4->Eval(x2_poly4)<<endl;}
else{
z2_poly4 = x2_poly4;
cout << "Zero found! z2: "<<z2_poly4<<", y2: "<<poly4->Eval(x2_poly4)<<endl;
b_poly4 = true;
}
}
//This has to be done before any other fit is performed
if(b_poly4 == true){
cout << "Poly4 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
integral_poly4_val = poly4->Integral(z1_poly4,z2_poly4)/20.;
integral_poly4_err = poly4->IntegralError(z1_poly4,z2_poly4)/20.;
cout << "The integral between the two zeros is "<< integral_poly4_val <<"+/-"<<integral_poly4_err << endl;
}
//POLY5
cout << "================ poly5 =================" << endl<<endl;
TF1 *poly5 = new TF1("poly5","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 + [5]*x**5",-40.,40.);
poly5->SetLineColor(6);
ge_pulse->Fit("poly5","M");
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly5->GetChisquare(),poly5->GetNDF());
double x1_poly5 = poly5->GetX(0,low_limit,0);
double x2_poly5 = poly5->GetX(0,0, high_limit);
double z1_poly5 = -99;
double z2_poly5 = -99;
bool b_poly5 = false;
double integral_poly5_val = -999;
double integral_poly5_err = 0;
if( TMath::Abs(poly5->Eval(x1_poly5)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly5->Eval(x1_poly5)<<endl;}
else{
z1_poly5 = x1_poly5;
cout <<"Zero found! z1: "<< z1_poly5<<", y1: "<<poly5->Eval(x1_poly5)<<endl;
if( TMath::Abs(poly5->Eval(x2_poly5)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly5->Eval(x2_poly5)<<endl;}
else{
z2_poly5 = x2_poly5;
cout << "Zero found! z2: "<<z2_poly5<<", y2: "<<poly5->Eval(x2_poly5)<<endl;
b_poly5 = true;
}
}
//This has to be done before any other fit is performed
if(b_poly5 == true){
cout << "Poly5 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
integral_poly5_val = poly5->Integral(z1_poly5,z2_poly5)/20.;
integral_poly5_err = poly5->IntegralError(z1_poly5,z2_poly5)/20.;
cout << "The integral between the two zeros is "<< integral_poly5_val <<"+/-"<<integral_poly5_err << endl;
}
//MIchele - End
/*
TF1* fu_pulse2 = STTool::halfGauss2;
ge_pulse->Fit(fu_pulse2,"MQ","0");
Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse2->GetChisquare(),fu_pulse2->GetNDF());
*/
TVirtualFitter *fit_volt = TVirtualFitter::GetFitter();
Double_t* cov_volt = fit_volt->GetCovarianceMatrix();
double maxPos = 3-TMath::Sqrt(3);
double maxFac = TMath::Exp(-maxPos)*(maxPos*maxPos/2.0-maxPos*maxPos*maxPos/6.0);
double max_val = maxFac*fu_pulse->GetParameter(0);
double max_err = maxFac*fu_pulse->GetParError(0);
double tau_val = fu_pulse->GetParameter(2);
double tau_err = fu_pulse->GetParError(2);
double t0_val = fu_pulse->GetParameter(1);
double t0_err = fu_pulse->GetParError(1);
// Calculate integrale; factor 1/20 is included such that the maximal value of the pulse and the integral have
// a similar numerical value
double int_val = 9.0/2.0*TMath::Exp(-3.0)*fu_pulse->GetParameter(0)*tau_val/20.0;
double int_err = int_val*TMath::Sqrt(TMath::Power(fu_pulse->GetParError(0)/fu_pulse->GetParameter(0),2)+
TMath::Power(tau_err/tau_val,2)+
2.0*cov_volt[2]/(tau_val*fu_pulse->GetParameter(0))
);
//Michele 3/9/2017
//Estimate the systematics: # if a polynomial with zeros is found take as sys the semidifference
// between its integral and the integral obtained from the fitted theoretical model
// # if no polynomial behaved properly set the systematic error to -999
double int_sys = -999.0;
double b_int_sys = true;
if(b_poly3 == true){
int_sys = TMath::Abs(integral_poly3_val - int_val)/2.;
}else{
if(b_poly4 == true){
int_sys = TMath::Abs(integral_poly4_val - int_val)/2.;
}else{
if(b_poly5 == true){
int_sys = TMath::Abs(integral_poly5_val - int_val)/2.;
}else{
b_int_sys = false;
}
}
}
//Michele - end
mg_pulse->Add(ge_pulse,"PE");
Info("pulseFit","pulse Height: %8.4f+/-%8.4f",max_val,max_err);
Info("pulseFit"," tau: %8.4f+/-%8.4f ns",tau_val,tau_err);
//Michele
if(b_int_sys == true){
Info("pulseFit"," charge eq: %8.4f +/- %8.4f(stat) +/- %8.4f (sys)",int_val,int_err, int_sys);
}else{
Info("pulseFit"," charge eq: %8.4f +/- %8.4f(stat) +/- No sys eval.",int_val,int_err);
}
//
// Plot pulse
TCanvas* c_pul = new TCanvas("c_pul","Pulse Canvas",800,600);
mg_pulse->Draw("A");
fu_pulse->Draw("Same");
//fu_pulse2->Draw("Same");
poly3->Draw("Same");
poly4->Draw("Same");
poly5->Draw("Same");
mg_pulse->Draw("");
mg_pulse->GetXaxis()->SetLimits(-30.0,30.0);
mg_pulse->GetXaxis()->SetTitle("#delta#it{t} [ns]");
mg_pulse->SetMaximum( 1.2*max_val);
mg_pulse->SetMinimum(0.0);
mg_pulse->GetYaxis()->SetTitle("MPV [ADC value]");
TLine* line = new TLine(-30.0,0.0,30.0,0.0);
//line->Draw();
mg_pulse->GetHistogram()->Draw("axissame");
// 0,1: Maximal signal height
v_res(0) = max_val;
v_res(1) = max_err;
//v_res(1) = max_sys;
// 2,3: tau (width of pulse)
v_res(2) = tau_val;
v_res(3) = tau_err;
// 4,5: t0 (first root)b
v_res(4) = t0_val;
v_res(5) = t0_err;
// 6,7: covariance (tau,A0) (t0,A0) (t0,tau)
v_res(8) = cov_volt[1];
v_res(9) = cov_volt[2];
v_res(10)= cov_volt[5];
v_res(11) = int_sys;
// Write the data and the plots
while (!gSystem->Exec(Form("ls %s &> /dev/null",fn_dummy.Data()))) {
Info("landauFit","Waiting until file is closed");
gSystem->Sleep(1000);
}
gSystem->Exec(Form("touch %s &> /dev/null",fn_dummy.Data()));
f_data = new TFile(fn_data.Data(),"UPDATE");
TString dn_result = TString::Format("%s/%s/%d/%d",
det.Data(),
lay.Data(),
sector,
fill);
// Get right directory
bool isThere = f_data->GetListOfKeys()->Contains(det.Data());
if (!isThere) {
f_data->mkdir(det.Data());
}
isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s",
det.Data(),
lay.Data()));
if (!isThere) {
f_data->mkdir(Form("%s/%s",
det.Data(),
lay.Data()));
}
isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d",
det.Data(),
lay.Data(),
sector));
if (!isThere) {
f_data->mkdir(Form("%s/%s/%d",
det.Data(),
lay.Data(),
sector));
}
isThere = f_data->GetListOfKeys()->Contains(Form("%s/%s/%d/%d",
det.Data(),
lay.Data(),
sector,
fill));
if (!isThere) {
f_data->mkdir(Form("%s/%s/%d/%d",
det.Data(),
lay.Data(),
sector,
fill));
}
f_data->cd(Form("%s/%s/%d/%d",
det.Data(),
lay.Data(),
sector,
fill));
if (!notSave) {
gDirectory->Delete(Form("%s;*",
vn_result.Data()));
v_res.Write(Form("%s",
vn_result.Data()));
Info("pulseFit","Result written to %s",fn_data.Data());
}
f_data->Close();
gSystem->Exec(Form("rm -f %s &> /dev/null",fn_dummy.Data()));
if (!notSave) {
// Save picture
gSystem->Exec(Form("mkdir -p %s",dn_graph.Data()));
gSystem->Exec(Form("rm %s &> /dev/null",fn_graph_p.Data()));
c_pul->SaveAs(fn_graph_p.Data());
}
Printf("=========================================\n");
return EXIT_SUCCESS;
}
