// // voltageFit_spline.C // Extraction of the depletion voltage based on a spline fit // // Created by Christian Elsasser on 31.05.13. // University of Zurich, elsasser@cern.ch // Copyright only for commercial use // // 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 "TExMap.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 // -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 voltageFit(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 = voltageFit(argc,argv); Printf("End"); if (!runBatch) { theApp->Run(!runBatch); } return exit; } // Executable method int voltageFit(int argc, char* argv[]){ TVectorD v_res(5); gStyle->SetPadLeftMargin(0.20); gStyle->SetTitleOffset(1.5,"Y"); gStyle->SetTitleOffset(1.0,"X"); const char* diskVar = std::getenv ("DISK"); const char* homeVar = std::getenv ("CCEHOME"); int fill = 2083; TString det("TT"); TString lay("TTaU"); int sector = 2634; int val = 3; bool notSave = false; bool temSave = false; bool use_height = false; bool average_ratio = false; double radius = -10.0; // Reading command line arguments for (int i=1; i<argc; i++) { if (strcmp(argv[i],"-d")==0 && i+1<argc) { // Detector det = argv[++i]; }else if (strcmp(argv[i],"-l")==0 && i+1<argc) { // Layer lay = argv[++i]; }else if (strcmp(argv[i],"-s")==0 && i+1<argc) { // Sector sector = std::atoi(argv[++i]); }else if (strcmp(argv[i],"-f")==0 && i+1<argc) { // Fill fill = std::atoi(argv[++i]); }else if (strcmp(argv[i],"-v")==0 && i+1<argc) { // Number of strips val = std::atoi(argv[++i]); }else if (strcmp(argv[i],"-r")==0 && i+1<argc) { // Radius radius = std::atof(argv[++i]); }else if (strcmp(argv[i],"-use_height")==0) { // use height of pulse shape use_height = true; }else if (strcmp(argv[i],"-average_ratio")==0) { // use the average ratio found for each layer and val average_ratio = true; }else if (strcmp(argv[i],"-ns")==0) { // Do not save notSave = true; }else if (strcmp(argv[i],"-ss")==0) { // Save it in a temporary file temSave = true; } } std::vector<double> v_volt_val; std::vector<double> v_volt_err; std::vector<double> v_adc_val; std::vector<double> v_adc_err; // Vector name where the result is stored TString vn_result(Form("v_volt_val%d",val)); // 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/volt_val%d", dn_graph.Data(),val)); if (radius<0.0) { }else{ vn_result +=Form("_r%d",(int)radius); fn_graph_p +=Form("_r%d",(int)radius); } if(use_height){ vn_result +=Form("_height"); fn_graph_p +=Form("_height"); }else{ vn_result +=Form("_int"); fn_graph_p +=Form("_int"); } if(average_ratio){ vn_result +=Form("_avgratio"); fn_graph_p +=Form("_avgratio"); }else{ vn_result +=Form("_indratio"); fn_graph_p +=Form("_indratio"); } fn_graph_p += ".pdf"; Printf("========================================="); Info("voltageFit_spline","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()); 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())); // TString dn_ratio(Form("%s/cmtuser/Vetra_v13r2/ST/STAging/data/db" // ,homeVar)); TString dn_ratio(Form("%s/data/db", homeVar)); TString fn_ratio(Form("%s/ratio", dn_ratio.Data())); if(use_height){ fn_ratio +=Form("_height.root"); }else{ fn_ratio +=Form("_int.root"); } TString fn_syst(Form("%s/systVolt.root", dn_ratio.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("voltageFit_spline","Data file does not exist!"); Printf("=========================================\n"); return EXIT_FAILURE; } TFile* f_ratio = TFile::Open(fn_ratio.Data()); if (!f_ratio) { Error("voltageFit_spline","Ratio data file does not exist!"); Printf("=========================================\n"); return EXIT_FAILURE; } // Load the data Info("voltageFit","Loading data:"); TString dn_input(Form("%s/%s/%d/%d", det.Data(), lay.Data(), sector, fill)); for (int i=0; i<11; i++) { TString vn_input(Form("%s/v_pulse%d_val%d",dn_input.Data(),i,val)); if (radius>0.0) { vn_input +=Form("_r%d",(int)radius); } TVectorD* vp_input = (TVectorD*)f_data->Get(vn_input.Data()); if (!vp_input) { Warning("voltageFit_spline","Vector %s not available",vn_input.Data()); continue; } double volt = 0.0; // Extract the voltage value if (det.EqualTo("TT")) { volt = STTool::GetTTVoltage(i*6); }else{ volt = STTool::GetITVoltage(i*6); } if (det.EqualTo("IT") && (volt>250.0 || volt<30.0)){ continue; } if (det.EqualTo("TT") && (volt<80.0)){ continue; } // Get the input values TVectorD v_input = (*vp_input); double height_val = v_input(0); double height_err = v_input(1); double pulse_quality = v_input(11); if(not use_height){ if (true || det.EqualTo("IT")) { height_val /= 0.1306; // transform back to the amplitude height_err /= 0.1306; // transform back to the amplitude height_val *= v_input(2)*9.0/2.0*TMath::Exp(-3.0)/(20.0); height_err = height_val*TMath::Sqrt(TMath::Power(height_err/height_val,2)+ TMath::Power(v_input(3)/v_input(2),2)+ 2*v_input(8)/(height_val*v_input(2))); } } Printf(" %8.2f V: %8.4f+/-%8.4f",volt,height_val,height_err); if (pulse_quality == -1){ cout << "Ignoring this point for the voltageFit..."<<endl; }else{ // Load the input values in to vector v_adc_val.push_back(height_val); v_adc_err.push_back(height_err); v_volt_val.push_back(volt); v_volt_err.push_back(0.0); } } if (v_adc_val.size()<3) { Error("voltageFit_spline","Too few (%d) data points",v_adc_val.size()); Printf("=========================================\n"); return EXIT_FAILURE; } // Add zero point v_adc_val.push_back(0.0); v_adc_err.push_back(0.5); v_volt_val.push_back(-1e-9); v_volt_err.push_back(0.0); // Get the corresponding root vectors TVectorD vr_adc_val = getROOTVector(v_adc_val); TVectorD vr_adc_err = getROOTVector(v_adc_err); TVectorD vr_volt_val = getROOTVector(v_volt_val); TVectorD vr_volt_err = getROOTVector(v_volt_err); TGraphErrors* ge_volt = new TGraphErrors(vr_volt_val,vr_adc_val,vr_volt_err,vr_adc_err); TMultiGraph* mg_volt = new TMultiGraph("mg_volt","Voltage multigraph"); // Extract the ratio values for the sampling point TVectorD* vp_ratio = (TVectorD*)f_ratio->Get(Form("v_ratio_%s_val%d",lay.Data(),val)); TVectorD* vp_ratio_err = (TVectorD*)f_ratio->Get(Form("v_ratio_err_%s_val%d",lay.Data(),val)); TVectorD* vp_secRatio = (TVectorD*)f_ratio->Get(Form("v_sector_%s_val%d",lay.Data(),val)); if (!vp_ratio || !vp_secRatio || !vp_ratio_err) { Error("voltageFit_spline","Ratio vectors not available"); Printf("=========================================\n"); return EXIT_FAILURE; } TVectorD v_ratio = *vp_ratio; TVectorD v_ratio_err = *vp_ratio_err; TVectorD v_secRatio = *vp_secRatio; TF1* fu_volt = STTool::Constfunc; double ratio_val = STTool::ITratio_val; double ratio_err = STTool::ITratio_err; if (det.EqualTo("TT")) { /* ratio_val = STTool::TTratio_val; ratio_err = STTool::TTratio_err; */ if(val == 3){ ratio_val = STTool::TTratio_val_3; ratio_err = STTool::TTratio_err_3; } else if(val == 5){ ratio_val = STTool::TTratio_val_5; ratio_err = STTool::TTratio_err_5; } else if(val == 7){ ratio_val = STTool::TTratio_val_7; ratio_err = STTool::TTratio_err_7; } } if (not average_ratio){ //if (false && det.EqualTo("TT")) { // Using ratio per sector //This is excluded by default //if (true) { int ratioIdx = getClosestIndex<double>(v_secRatio,(double)sector); if (!(ratioIdx<0) && (int)v_secRatio(ratioIdx)==sector && true) { Info("voltageFit_spline","Using specific ratio for sector %d",sector); ratio_val = v_ratio(ratioIdx); ratio_err = v_ratio_err(ratioIdx); Printf(" Individual ratio: %f",ratio_val); } } Printf(" Ratio: %f",ratio_val); // Define fit range double fitMax = 510.0; double fitMin = 280.0; if (det.EqualTo("IT")) { fitMax = 280.0; fitMin = 160.0; } ge_volt->Fit(fu_volt,"M0","0",fitMin,fitMax); ge_volt->Fit(fu_volt,"M0","0",fitMin,fitMax); double chi2ndf = fu_volt->GetChisquare()/fu_volt->GetNDF(); // Adding 0/0 and saturation point to the data set v_adc_val.insert(v_adc_val.begin(),fu_volt->GetParameter("A_{0}")); v_adc_err.insert(v_adc_err.begin(),fu_volt->GetParError(fu_volt->GetParNumber("A_{0}"))); v_volt_val.insert(v_volt_val.begin(),fitMax); v_volt_err.insert(v_volt_err.begin(),0.0); // Reverse vector to be compatible with TSpline std::reverse(v_adc_val.begin(), v_adc_val.end()); std::reverse(v_adc_err.begin(), v_adc_err.end()); std::reverse(v_volt_val.begin(), v_volt_val.end()); std::reverse(v_volt_err.begin(), v_volt_err.end()); // Get the new corresponding root vectors TVectorD vr_adc_val_corr = getROOTVector(v_adc_val); TVectorD vr_adc_err_corr = getROOTVector(v_adc_err); TVectorD vr_volt_val_corr = getROOTVector(v_volt_val); TVectorD vr_volt_err_corr = getROOTVector(v_volt_err); // Scale them Info("voltageFit_spline","#chi^{2}/ndf: %6.2f/%d",fu_volt->GetChisquare(),fu_volt->GetNDF()); /*Michele if (chi2ndf>0.00) { Info("voltageFit_spline","scaling"); vr_adc_err_corr*=TMath::Sqrt(chi2ndf); } //Michele*/ // Get new graph ge_volt = new TGraphErrors(vr_volt_val_corr,vr_adc_val_corr,vr_volt_err_corr,vr_adc_err_corr); // Create Spline with zero first and second derivate TSpline5* s5_volt = new TSpline5("s5_volt",ge_volt,"b2 e1 e2",0.0,0.0,0.0,0.0); s5_volt->SetLineColor(kRed); // Get the value for the depletion voltage double tarADC_val = ratio_val; double tarADC_sys = ratio_err; tarADC_val *= fu_volt->GetParameter("A_{0}"); tarADC_sys *= fu_volt->GetParameter("A_{0}"); // Get depletion voltage double Vdepl_val = STTool::GetX(s5_volt,tarADC_val,10.0,400.0); // Get the uncertainty on the sampling point // with 1/distance^2 weighted uncertainty on the ADC values of the two closest data points int iFirst = -1; int iSecon = -1; // Excluding first and last point as they have been artifically added for(int i=1; i<vr_adc_err_corr.GetNoElements()-1; i++){ if (iFirst==-1) { iFirst = i; }else{ if (TMath::Abs(vr_volt_val_corr(i)-Vdepl_val)<TMath::Abs(vr_volt_val_corr(iFirst)-Vdepl_val)) { iFirst = i; continue; // Successful -> do not fill second } } if (iSecon==-1) { iSecon = i; }else{ if (TMath::Abs(vr_volt_val_corr(i)-Vdepl_val)<TMath::Abs(vr_volt_val_corr(iSecon)-Vdepl_val)) { iSecon = i; } } } double tarADC_err = vr_adc_err_corr(iFirst)/TMath::Power(vr_volt_val_corr(iFirst)-Vdepl_val,2)+vr_adc_err_corr(iSecon)/TMath::Power(vr_volt_val_corr(iSecon)-Vdepl_val,2); tarADC_err /= (1.0/TMath::Power(vr_volt_val_corr(iFirst)-Vdepl_val,2)+1.0/TMath::Power(vr_volt_val_corr(iSecon)-Vdepl_val,2)); Printf("tarADC_err: %f",tarADC_err); // Translate the uncertainty to the uncertainty on the depletion voltage double Vdepl_err_high = STTool::GetX(s5_volt,tarADC_val+tarADC_err,10.0,400.0); double Vdepl_err_low = STTool::GetX(s5_volt,tarADC_val-tarADC_err,10.0,400.0); double Vdepl_sys_high = STTool::GetX(s5_volt,tarADC_val+tarADC_sys,10.0,400.0); double Vdepl_sys_low = STTool::GetX(s5_volt,tarADC_val-tarADC_sys,10.0,400.0); double Vdepl_err = TMath::Abs((Vdepl_err_high-Vdepl_err_low)/2.0); double Vdepl_sys = TMath::Abs((Vdepl_sys_high-Vdepl_sys_low)/2.0); // Alternative determination of the systematic uncertainty by taking the // difference between the determined depletion voltage and the crossing point of the horizontal // top and the line extrapolated from (0,0) double Vdepl_sys1 = Vdepl_val-fu_volt->GetParameter("A_{0}")/s5_volt->Derivative(0.0); // Alternative determination of the systematic uncertainty by taking the // difference between the determined depletion voltage and the extracted value from the sigmoid function TF1* fu_voltAlt = STTool::sigmoid4; ge_volt->Fit(fu_voltAlt,"M0","0"); double Vdepl_sys2 = Vdepl_val-fu_voltAlt->GetX(tarADC_val,10.0,400.0);; // Draw the result mg_volt->Add(ge_volt,"pe"); TCanvas* c_vol = new TCanvas("c_vol","Voltage Canvas",800,600); double maxX = 500.0; if (det.EqualTo("IT")) { maxX = 250.0; } double maxY = 75.0; // 50.0; mg_volt->Draw("A"); mg_volt->GetXaxis()->SetLimits(0.0,maxX); mg_volt->GetXaxis()->SetTitle("#it{V}_{bias} [V]"); mg_volt->SetMaximum(maxY); mg_volt->SetMinimum(0.0); mg_volt->GetYaxis()->SetTitle("Charge Equivalent [ADC value]"); Int_t err_n = 2 * (Int_t)maxY + 3; Double_t err_x[err_n], err_y[err_n]; for (int i=0; i<err_n-2; i++) { err_x[i] = Vdepl_val + (i-maxY)/maxY*TMath::Sqrt(Vdepl_err*Vdepl_err+Vdepl_sys*Vdepl_sys); err_y[i] = s5_volt->Eval(err_x[i]); } err_x[err_n-2] = err_x[err_n-3]; err_x[err_n-1] = err_x[0]; err_y[err_n-2] = 0.0; err_y[err_n-1] = 0.0; TPolyLine *err_poly = new TPolyLine(err_n,err_x,err_y); err_poly->SetFillColor(kGray); err_poly->SetLineColor(kGray); err_poly->SetLineWidth(0); err_poly->Draw("f"); for (int i=0; i<err_n-2; i++) { err_x[i] = Vdepl_val + (i-maxY)/maxY*Vdepl_err; err_y[i] = s5_volt->Eval(err_x[i]); } err_x[err_n-2] = err_x[err_n-3]; err_x[err_n-1] = err_x[0]; err_y[err_n-2] = 0.0; err_y[err_n-1] = 0.0; err_poly = new TPolyLine(err_n,err_x,err_y); err_poly->SetFillColor(kGray+1); err_poly->SetLineColor(kGray+1); err_poly->SetLineWidth(0); err_poly->Draw("f"); TLine* val_line = new TLine(Vdepl_val,0.0,Vdepl_val,s5_volt->Eval(Vdepl_val)); val_line->DrawClone(); val_line = new TLine(Vdepl_val,0.0,Vdepl_val,maxY); val_line->SetLineStyle(2); val_line->DrawClone(); Printf("Drawing"); TF1* fu_lin = new TF1("f_lin","[0]*x",-100,800); fu_lin->SetParameter(0,s5_volt->Derivative(0.0)); fu_lin->SetLineStyle(2); //fu_lin->DrawCopy("Same"); fu_volt->SetLineColor(kBlack); fu_volt->SetLineStyle(2); fu_volt->DrawCopy("Same"); s5_volt->Draw("Same"); mg_volt->Draw(""); mg_volt->GetHistogram()->Draw("axissame"); Info("voltageFit","V_{depl} = (%8.4f+/-%8.4f(stat.)+/-%8.4f(syst.)) V",Vdepl_val,Vdepl_err,Vdepl_sys); // Load map /* TFile* f_map = new TFile(fn_syst.Data(),"UPDATE"); TExMap* m_lin = (TExMap*)f_map->Get(Form("m_lin_%s",lay.Data())); // Deviations from linear approach TExMap* m_sig = (TExMap*)f_map->Get(Form("m_sig_%s",lay.Data())); // Deviations from sigmoid approach if(!m_lin){ Info("voltageFit_spline","Linear map does not exist! Creating it..."); m_lin = new TExMap(1000); } if(!m_sig){ Info("voltageFit_spline","Sigmoid map does not exist! Creating it..."); m_sig = new TExMap(1000); } (*m_lin)((Long64_t)(sector+10000*fill)) = *(reinterpret_cast<Long64_t*>(&Vdepl_sys1)); (*m_sig)((Long64_t)(sector+10000*fill)) = *(reinterpret_cast<Long64_t*>(&Vdepl_sys2)); Info("voltageFit_spline","Deviation to linear approach: %8.4f V",Vdepl_sys1); Info("voltageFit_spline","Deviation to sigmoid approach: %8.4f V",Vdepl_sys2); f_map->Delete(Form("m_lin_%s;*",lay.Data())); f_map->Delete(Form("m_sig_%s",lay.Data())); m_lin->Write(Form("m_lin_%s",lay.Data())); m_sig->Write(Form("m_sig_%s",lay.Data())); f_map->Close(); */ // Save the result v_res(0) = Vdepl_val; v_res(1) = Vdepl_err; v_res(2) = Vdepl_sys; v_res(3) = fu_volt->GetParameter("A_{0}"); v_res(4) = fu_volt->GetParError(fu_volt->GetParNumber("A_{0}")); while (!gSystem->Exec(Form("ls %s &> /dev/null",fn_dummy.Data()))) { Info("voltageFit_spline","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("voltageFit_spline","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_vol->SaveAs(fn_graph_p.Data()); } Printf("=========================================\n"); return EXIT_SUCCESS; }