SpaceResolutionCalibration Class Reference

Class for Space resolution calibration. More...

#include <SpaceResolutionCalibration.h>

Public Member Functions
	SpaceResolutionCalibration ()
	Constructor.
virtual	~SpaceResolutionCalibration ()
	Destructor.
virtual void	setDebug (bool debug=false)
	Debug or not.
virtual void	setUseDB (bool useDB=false)
	Use database or text mode.
virtual void	setMinimumNDF (double minndf)
	minimum NDF required for track
virtual void	setMinimumPval (double minPval)
	Minimum Pval required.
virtual void	setBinWidth (double bw)
	Bin width of each slide.
virtual void	BField (bool bfield)
	Work with B field or not;.
virtual void	inputFileNames (std::string inputname)
	Input root file names, results of collector module.
virtual void	setStoreHisto (bool storeHist=false)
	Store histograms during the calibration or not.
virtual void	ProfileFileNames (std::string profileFileName)
	File name describe theta/alpha bin, if don't want to use default from input sigma.
virtual void	useProfileFromInputSigma (bool useProfileFromInputSigma)
	use sigma bin profile form input sigma or new one from input file
virtual void	setSigmaFileName (std::string name)
	Output sigma file name, for text mode.
void	execute ()
	execute all, make the interface the same as CAF

Protected Member Functions
virtual bool	calibrate ()
	Run algo on data.
virtual void	createHisto ()
	create histogram
virtual void	readProfile ()
	read sigma bining (alpha, theta bining)
virtual void	readSigma ()
	read sigma from previous calibration, (input sigma)
virtual void	readSigmaFromDB ()
	read sigma from DB
virtual void	readSigmaFromText ()
	read sigma from text file
virtual void	storeHisto ()
	store histogram
virtual void	write ()
	save calibration, in text file or db

Private Member Functions
double	getUpperBoundaryForFit (TGraphErrors *graph)
	search max point at boundary region

Private Attributes
double	m_ndfmin = 5
	Minimum NDF.
double	m_Pvalmin = 0.
	Minimum Prob(chi2) of track.
double	m_binWidth = 0.05
	width of each bin, unit cm
bool	m_debug = false
	Debug or not.
bool	m_storeHisto = false
	Store histogram or not.
bool	m_useDB = false
	use db or text mode
bool	m_useProfileFromInputSigma = true
	Use binning from old sigma or new one form input.
bool	m_BField = true
	Work with BField, fit range and initial parameters is different incase B and noB.
double	sigma_old [56][2][18][7][8]
	old sigma parameters.
double	sigma_new [56][2][18][7][8]
	new sigma parameters.
TGraphErrors *	gfit [56][2][18][7]
	sigma*sigma graph for fit
TGraphErrors *	gr [56][2][18][7]
	sigma graph.
TH2F *	hist_b [56][2][Max_nalpha][Max_ntheta]
	2D histogram of biased residual
TH2F *	hist_u [56][2][Max_nalpha][Max_ntheta]
	2D histogram of unbiased residual
TH1F *	hu_m [56][2][Max_nalpha][Max_ntheta]
	mean histogram biased residual
TH1F *	hu_s [56][2][Max_nalpha][Max_ntheta]
	sigma histogram of biased residual
TH1F *	hb_m [56][2][Max_nalpha][Max_ntheta]
	mean histogram of unbiased residual
TH1F *	hb_s [56][2][Max_nalpha][Max_ntheta]
	sigma histogram of ubiased residual
int	m_fitflag [56][2][Max_nalpha][Max_ntheta] = {{{{0}}}}
	Fit flag; 1:OK ; 0:error.
std::string	m_outputSigmaFileName = "sigma_new.dat"
	Output sigma file name.
std::string	m_inputRootFileNames = "rootfile/output*"
	Input root file names.
std::string	m_ProfileFileName = "sigma_profile"
	Profile file name.
DBObjPtr< CDCSpaceResols > *	m_sResolFromDB
	Database for sigma.
std::string	m_sigmafile = "cdc/data/sigma.dat"
	Sigma file name, for text mode.
int	m_nalpha
	number of alpha bins
int	m_ntheta
	number of theta bins
double	l_alpha [18]
	Lower boundaries of alpha bins.
double	u_alpha [18]
	Upper boundaries of alpha bins.
double	ialpha [18]
	represented alphas of alpha bins.
double	l_theta [7]
	Lower boundaries of theta bins.
double	u_theta [7]
	Upper boundaries of theta bins.
double	itheta [7]
	represented alphas of theta bins.
int	nalpha_old
	number of alpha bins from input
int	ntheta_old
	number of theta bins from input
double	l_alpha_old [18]
	Lower boundaries of alpha bins from input.
double	u_alpha_old [18]
	Upper boundaries of alpha bins from input.
double	ialpha_old [18]
	represented alphas of alpha bins from input.
double	l_theta_old [7]
	Lower boundaries of theta bins from input.
double	u_theta_old [7]
	Upper boundaries of theta bins from input.
double	itheta_old [7]
	represented alphas of theta bins from input.
unsigned short	m_sigmaParamMode_old
	sigma mode from input.

Static Private Attributes
static const int	Max_nalpha = 18
	Maximum alpha bin.
static const int	Max_ntheta = 7
	maximum theta bin
static const unsigned short	Max_np = 40
	Maximum number of point =1/binwidth.

Detailed Description

Class for Space resolution calibration.

Definition at line 25 of file SpaceResolutionCalibration.h.

Constructor & Destructor Documentation

SpaceResolutionCalibration()

SpaceResolutionCalibration

(

)

Constructor.

Definition at line 32 of file SpaceResolutionCalibration.cc.

{
  /*Space resolution calibration*/
}

~SpaceResolutionCalibration()

virtual ~SpaceResolutionCalibration ( )

inlinevirtual

Destructor.

Definition at line 31 of file SpaceResolutionCalibration.h.

{}

Member Function Documentation

BField()

virtual void BField ( bool bfield )

inlinevirtual

Work with B field or not;.

Definition at line 43 of file SpaceResolutionCalibration.h.

{m_BField = bfield;}

calibrate()

bool calibrate ( )

protectedvirtual

Run algo on data.

Upper limit of fitting.

Definition at line 286 of file SpaceResolutionCalibration.cc.

{
  gROOT->SetBatch(1);
  gErrorIgnoreLevel = 3001;
  createHisto();
  B2INFO("Start to calibrate");
  gSystem->Exec("mkdir -p Sigma_Fit_Err"); //create a folder to store error histo
 
  TF1* func = new TF1("func", "[0]/(x*x + [1])+[2]* x+[3]+[4]*exp([5]*(x-[6])*(x-[6]))", 0, 1.);
  TH1F* hprob = new TH1F("h1", "", 20, 0, 1);
  double upFit; 
  double intp6;
 
  for (int i = 0; i < 56; ++i) {
    for (int lr = 0; lr < 2; ++lr) {
      for (int al = 0; al < m_nalpha; ++al) {
        for (int th = 0; th < m_ntheta; ++th) {
          if (m_fitflag[i][lr][al][th] != -1) { /*if graph exist, do fitting*/
 
 
            upFit = getUpperBoundaryForFit(gfit[i][lr][al][th]);
            intp6 = upFit + 0.2;
            B2DEBUG(199, "xmax for fitting: " << upFit);
            //          func->SetLineColor(1 + lr + al * 2 + th * 3);
            func->SetParameters(5E-6, 0.007, 1E-4, 1E-5, 0.00008, -30, intp6);
            func->SetParLimits(0, 1E-7, 1E-4);
            func->SetParLimits(1, 0.00001, 0.02);
            func->SetParLimits(2, 1E-6, 0.0005);
            func->SetParLimits(3, 1E-8, 0.0005);
            func->SetParLimits(4, 0., 0.001);
            func->SetParLimits(5, -40, 0.);
            func->SetParLimits(6, intp6 - 0.5, intp6 + 0.3);
            B2DEBUG(199, "FITTING for layer: " << i << "lr: " << lr << " ial" << al << " ith:" << th);
            B2DEBUG(199, "Fit flag before fit:" << m_fitflag[i][lr][al][th]);
            //    if(!(gfit[i][lr][al][th]->isValid())) continue;
            //      m_fitflag[i][lr][al][th] = 0;
            for (int j = 0; j < 10; j++) {
 
              B2DEBUG(199, "loop: " << j);
              B2DEBUG(199, "Int p6: " << intp6);
              B2DEBUG(199, "Number of Point: " << gfit[i][lr][al][th]->GetN());
              Int_t stat = gfit[i][lr][al][th]->Fit("func", "MQE", "", 0.05, upFit);
              B2DEBUG(199, "stat of fit" << stat);
              std::string Fit_status = gMinuit->fCstatu.Data();
              B2DEBUG(199, "FIT STATUS: " << Fit_status);
              //      stat=gfit[i]->Fit(Form("ffit[%d]",i),"M "+Q,"",0.0,cellsize(i)+0.05+j*0.005);
              if (Fit_status == "OK" ||
                  Fit_status == "SUCCESSFUL" ||
                  Fit_status == "CALL LIMIT" ||
                  Fit_status == "PROBLEMS") {
                if (fabs(func->Eval(0.3)) > 0.00035 || func->Eval(0.3) < 0) {
                  func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
                  //    func->SetParameters(defaultparsmall);
                  m_fitflag[i][lr][al][th] = 0;
                } else {
                  B2DEBUG(199, "Prob of fit: " << func->GetProb());
                  m_fitflag[i][lr][al][th] = 1;
                  break;
                }
              } else {
                m_fitflag[i][lr][al][th] = 0;
                func->SetParameters(5E-6, 0.007, 1E-4, 1E-7, 0.0007, -30, intp6 + 0.05 * j);
                //        upFit += 0.025;
                if (j == 9) {
                  TCanvas* c1 =  new TCanvas("c1", "", 600, 600);
                  gfit[i][lr][al][th]->Draw();
                  c1->SaveAs(Form("Sigma_Fit_Err/%d_%d_al%d_th%d_%s.png", i, lr, al, th, Fit_status.c_str()));
                  B2WARNING("Fit error: " << i << " " << lr << " " << al << " " << th);
                }
              }
            }
            if (m_fitflag[i][lr][al][th] == 1) {
              B2DEBUG(199, "ProbFit: Lay_lr_al_th: " << i << " " << lr << " " << al << " " << th << func->GetProb());
              hprob->Fill(func->GetProb());
              func->GetParameters(sigma_new[i][lr][al][th]);
            }
          }
        }
      }
    }
  }
  storeHisto();
  write();
 
  return true;
}

createHisto()

void createHisto ( )

protectedvirtual

create histogram

Definition at line 37 of file SpaceResolutionCalibration.cc.

{
 
  B2INFO("createHisto");
  readSigma();
  readProfile();
  const int m_np = floor(1 / m_binWidth);
 
  TChain* tree = new TChain("tree");
  tree->Add(m_inputRootFileNames.c_str());
  B2INFO(" Open file name: " << m_inputRootFileNames.c_str());
  if (!tree->GetBranch("ndf")) {
    B2FATAL("input data do not exits, please check!");
    gSystem->Exec("echo rootfile do not exits or something wrong >> error");
    return;
  }
 
  int lay;
  double w;
  double x_u;
  double x_b;
  double x_mea;
  double Pval;
  double alpha;
  double theta;
  double ndf;
  double absRes_u;
  double absRes_b;
  tree->SetBranchAddress("lay", &lay);
  tree->SetBranchAddress("ndf", &ndf);
  tree->SetBranchAddress("Pval", &Pval);
  tree->SetBranchAddress("x_u", &x_u);
  tree->SetBranchAddress("x_b", &x_b);
  tree->SetBranchAddress("x_mea", &x_mea);
  tree->SetBranchAddress("weight", &w);
  tree->SetBranchAddress("alpha", &alpha);
  tree->SetBranchAddress("theta", &theta);
 
  /* Disable unused branch */
  std::vector<TString> list_vars = {"lay", "ndf", "Pval", "x_u", "x_b", "x_mea", "weight", "alpha", "theta"};
  tree->SetBranchStatus("*", 0);
 
  for (TString brname : list_vars) {
    tree->SetBranchStatus(brname, 1);
  }
 
 
  vector<double> yu;
  vector <double> yb;
  for (int i = 0; i < 50; ++i) {
    yb.push_back(-0.07 + i * (0.14 / 50));
  }
  for (int i = 0; i < 50; ++i) {
    yu.push_back(-0.08 + i * (0.16 / 50));
  }
 
  vector<double> xbin;
  xbin.push_back(0.);
  xbin.push_back(0.02);
  for (int i = 1; i < m_np; ++i) {
    xbin.push_back(i * m_binWidth);
  }
 
  for (int il = 0; il < 56; ++il) {
    for (int lr = 0; lr < 2; ++lr) {
      for (int al = 0; al < m_nalpha; ++al) {
        for (int th = 0; th < m_ntheta; ++th) {
          hist_b[il][lr][al][th] = new TH2F(Form("hb_%d_%d_%d_%d", il, lr, al, th),
                                            Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, ialpha[al], itheta[th]),
                                            xbin.size() - 1, &xbin.at(0), yb.size() - 1, &yb.at(0));
          hist_u[il][lr][al][th] = new TH2F(Form("hu_%d_%d_%d_%d", il, lr, al, th),
                                            Form("lay_%d_lr%d_al_%3.0f_th_%3.0f;Drift Length [cm];#DeltaX", il, lr, ialpha[al], itheta[th]),
                                            xbin.size() - 1, &xbin.at(0), yu.size() - 1, &yu.at(0));
        }
      }
    }
  }
 
 
  const int nEntries = tree->GetEntries();
  B2INFO("Number of entries: " << nEntries);
  int ith = -99;
  int ial = -99;
  int ilr = -99;
  for (int i = 0; i < nEntries; ++i) {
    tree->GetEntry(i);
    //cut
    if (std::fabs(x_b) < 0.02 || std::fabs(x_u) < 0.02) continue;
    if (Pval < m_Pvalmin) continue;
    if (ndf < m_ndfmin) continue;
    for (int k = 0; k < m_nalpha; ++k) {
      if (alpha < u_alpha[k]) {
        ial = k;
        break;
      }
    }
 
    for (int j = 0; j < m_ntheta; ++j) {
      if (theta < u_theta[j]) {
        ith = j;
        break;
      }
    }
 
    ilr = x_u > 0 ? 1 : 0;
 
    if (ial == -99 || ith == -99 || ilr == -99) {
      TString command = Form("Error in alpha=%3.2f and theta = %3.2f>> error", alpha, theta);
      //      gSystem->Exec(command);
      B2FATAL("ERROR" << command);
    }
    absRes_u = fabs(x_mea) - fabs(x_u);
    absRes_b = fabs(x_mea) - fabs(x_b);
 
    hist_u[lay][ilr][ial][ith]->Fill(fabs(x_u), absRes_u, w);
    hist_b[lay][ilr][ial][ith]->Fill(fabs(x_b), absRes_b, w);
 
  }
 
  B2INFO("Finish reading data");
 
  TF1* gb = new TF1("gb", "gaus", -0.05, 0.05);
  TF1* gu = new TF1("gu", "gaus", -0.06, 0.06);
  TF1* g0b = new TF1("g0b", "gaus", -0.015, 0.07);
  TF1* g0u = new TF1("g0u", "gaus", -0.015, 0.08);
  g0b->SetParLimits(1, -0.01, 0.004);
  g0u->SetParLimits(1, -0.01, 0.004);
 
  std::vector<double> sigma;
  std::vector<double> dsigma;
  std::vector<double> s2;
  std::vector<double> ds2;
  std::vector<double> xl;
  std::vector<double> dxl;
  std::vector<double> dxl0;
 
  const int ib1 = int(0.1 / m_binWidth) + 1;
  int firstbin = 1;
  int minEntry = 10;
  for (int il = 0; il < 56; ++il) {
    for (int lr = 0; lr < 2; ++lr) {
      for (int al = 0; al < m_nalpha; ++al) {
        for (int th = 0; th < m_ntheta; ++th) {
          //fit half gaus for first range near sense wire
          B2DEBUG(199, "layer-lr-al-th " << il << " - " << lr << " - " << al << " - " << th);
          if (hist_b[il][lr][al][th]->GetEntries() < 20000) {
            m_fitflag[il][lr][al][th] = -1;
            continue;
          }
          B2DEBUG(199, "Nentries: " << hist_b[il][lr][al][th]->GetEntries());
          hist_b[il][lr][al][th]->SetDirectory(0);
          hist_u[il][lr][al][th]->SetDirectory(0);
 
          hist_b[il][lr][al][th]->FitSlicesY(g0b, firstbin, ib1, minEntry);
 
          TH1F* hm1 = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th));
          TH1F* hs1 = (TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th));
          if (!hm1 || !hs1) {
            m_fitflag[il][lr][al][th] = -1;
            continue;
          }
 
          hb_m[il][lr][al][th] = (TH1F*)hm1->Clone(Form("hb_%d_%d_%d_%d_m", il, lr, al, th));
          hb_s[il][lr][al][th] = (TH1F*)hs1->Clone(Form("hb_%d_%d_%d_%d_s", il, lr, al, th));//sigma
 
          hb_m[il][lr][al][th]->SetDirectory(0);
          hb_s[il][lr][al][th]->SetDirectory(0);
 
 
          //slice other bin with full gaus func
          hist_b[il][lr][al][th]->FitSlicesY(gb, ib1 + 1, m_np, minEntry);
          hb_m[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_1", il, lr, al, th))); //mean
          hb_s[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hb_%d_%d_%d_%d_2", il, lr, al, th))); //sigma
          B2DEBUG(199, "entries (2nd): " << hb_s[il][lr][al][th]->GetEntries());
 
          //fit half gaus for first range near sense wire
          hist_u[il][lr][al][th]->FitSlicesY(g0u, firstbin, ib1, minEntry);
          hu_m[il][lr][al][th] = (TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th))->Clone(Form("hu_%d_%d_%d_%d_m", il, lr, al,
                                 th));//mean
          hu_s[il][lr][al][th] = (TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th))->Clone(Form("hu_%d_%d_%d_%d_s", il, lr, al,
                                 th));//sigma
          hu_m[il][lr][al][th]->SetDirectory(0);
          hu_s[il][lr][al][th]->SetDirectory(0);
 
          //slice other bin with full gaus func
          hist_u[il][lr][al][th]->FitSlicesY(gu, ib1 + 1, m_np, minEntry);
          hu_m[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_1", il, lr, al, th))); //mean
          hu_s[il][lr][al][th]->Add((TH1F*)gDirectory->Get(Form("hu_%d_%d_%d_%d_2", il, lr, al, th))); //sigma
          if (!hu_s[il][lr][al][th] || !hb_s[il][lr][al][th]) {
            B2WARNING("slice histo do not found");
            m_fitflag[il][lr][al][th] = -1;
            continue;
          }
          //clean up container before adding new values
          xl.clear(); dxl.clear(); dxl0.clear(); sigma.clear(); dsigma.clear(); s2.clear(); ds2.clear();
          for (Int_t j = 1; j < hu_s[il][lr][al][th]->GetNbinsX(); j++) {
            if (hu_s[il][lr][al][th]->GetBinContent(j) == 0) continue;
            if (hb_s[il][lr][al][th]->GetBinContent(j) == 0) continue;
            double sb = hb_s[il][lr][al][th]->GetBinContent(j);
            double su = hu_s[il][lr][al][th]->GetBinContent(j);
 
            double dsb = hb_s[il][lr][al][th]->GetBinError(j);
            double dsu = hu_s[il][lr][al][th]->GetBinError(j);
            double XL = hb_s[il][lr][al][th]->GetXaxis()->GetBinCenter(j);
            double dXL = (hb_s[il][lr][al][th]->GetXaxis()->GetBinWidth(j)) / 2;
            double  s_int = std::sqrt(sb * su);
            double  ds_int = 0.5 * s_int * (dsb / sb + dsu / su);
            if (ds_int > 0.02) continue;
            xl.push_back(XL);
            dxl.push_back(dXL);
            dxl0.push_back(0.);
            sigma.push_back(s_int);
            dsigma.push_back(ds_int);
            s2.push_back(s_int * s_int);
            ds2.push_back(2 * s_int * ds_int);
          }
 
          if (xl.size() < 7 || xl.size() > Max_np) {
            m_fitflag[il][lr][al][th] = -1;
            B2WARNING("number of element might out of range"); continue;
          }
 
          //Intrinsic resolution
          B2DEBUG(199, "Create Histo for layer-lr: " << il << " " << lr);
          gr[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &sigma.at(0), &dxl.at(0), &dsigma.at(0));
          gr[il][lr][al][th]->SetMarkerSize(0.5);
          gr[il][lr][al][th]->SetMarkerStyle(8);
          //          gr[il][lr][al][th]->SetMarkerColor(kBlack);
          //gr[il][lr][al][th]->SetLineColor(1 + lr + al * 2 + th * 3);
          gr[il][lr][al][th]->SetTitle(Form("Layer_%d_lr%d | #alpha = %3.0f | #theta = %3.0f", il, lr, ialpha[al], itheta[th]));
          gr[il][lr][al][th]->SetName(Form("lay%d_lr%d_al%d_th%d", il, lr, al, th));
 
          //s2 for fitting
          gfit[il][lr][al][th] = new TGraphErrors(xl.size(), &xl.at(0), &s2.at(0), &dxl0.at(0), &ds2.at(0));
          gfit[il][lr][al][th]->SetMarkerSize(0.5);
          gfit[il][lr][al][th]->SetMarkerStyle(8);
          //          gfit[il][lr][al][th]->SetMarkerColor(1 + lr + al * 2 + th * 3);
          //gfit[il][lr][al][th]->SetLineColor(1 + lr + al * 2 + th * 3);
          gfit[il][lr][al][th]->SetTitle(Form("L%d-lr%d | #alpha = %3.0f | #theta = %3.0f ", il, lr, ialpha[al], itheta[th]));
          gfit[il][lr][al][th]->SetName(Form("sigma2_lay%d_lr%d_al%d_th%d", il, lr, al, th));
 
 
          gDirectory->Delete("hu_%d_%d_%d_%d_0");
        }
      }
    }
  }
}

execute()

void execute ( )

inline

execute all, make the interface the same as CAF

Definition at line 65 of file SpaceResolutionCalibration.h.

      {
        calibrate();
      }

getUpperBoundaryForFit()

double getUpperBoundaryForFit ( TGraphErrors * graph )

inlineprivate

search max point at boundary region

Definition at line 151 of file SpaceResolutionCalibration.h.

      {
        double ymax = 0;
        double xmax = 0;
        int imax = 0;
        double x, y;
        int unCount = floor(0.05 / m_binWidth);
        int N = graph->GetN();
        int Nstart = floor(0.5 * (N - unCount));
        int Nend = N - unCount;
        for (int i  = Nstart; i < Nend; ++i) {
          graph->GetPoint(i, x, y);
          if (graph->GetErrorY(i) > 0.06E-3) continue;
          if (y > ymax) {
            xmax = x; ymax = y;
            imax = i;
          }
        }
        if (imax <= Nstart) {
          graph->GetPoint(Nend, x, y);
          xmax = x;
        }
        return xmax;
      }

inputFileNames()

virtual void inputFileNames ( std::string inputname )

inlinevirtual

Input root file names, results of collector module.

Definition at line 45 of file SpaceResolutionCalibration.h.

      {
        m_inputRootFileNames.assign(inputname);
      }

ProfileFileNames()

virtual void ProfileFileNames ( std::string profileFileName )

inlinevirtual

File name describe theta/alpha bin, if don't want to use default from input sigma.

Definition at line 52 of file SpaceResolutionCalibration.h.

      {
        m_ProfileFileName.assign(profileFileName);
      }

readProfile()

void readProfile ( )

protectedvirtual

read sigma bining (alpha, theta bining)

Definition at line 612 of file SpaceResolutionCalibration.cc.

{
  B2INFO("readProfile");
  /*Read profile for xt*/
  if (m_useProfileFromInputSigma) {
    B2INFO("use Sigma bining from input Sigma");
    m_nalpha = nalpha_old;
    m_ntheta = ntheta_old;
    B2INFO("Number of alpha bins: " << m_nalpha);
    for (int i = 0; i < m_nalpha; ++i) {
      l_alpha[i] = l_alpha_old[i];    u_alpha[i] = u_alpha_old[i];     ialpha[i] = ialpha_old[i];
      B2INFO("" << i << " | " << l_alpha[i] << " " << u_alpha[i] << " " << ialpha[i]);
    }
    B2INFO("Number of theta bins: " << m_ntheta);
    for (int i = 0; i < m_ntheta; ++i) {
      l_theta[i] = l_theta_old[i];    u_theta[i] = u_theta_old[i];     itheta[i] = itheta_old[i];
      B2INFO("" << i << " |" << l_theta[i] << " " << u_theta[i] << " " << itheta[i]);
    }
  } else {
    B2INFO("use Sigma bining from profile file");
    ifstream proxt(m_ProfileFileName.c_str());
    if (!proxt) {
      B2FATAL("file not found: " << m_ProfileFileName);
    }
    double dumy1, dumy2, dumy3;
    proxt >> m_nalpha;
    B2INFO("Number of alpha bins: " << m_nalpha);
    if (m_nalpha > Max_nalpha) {B2FATAL("number of  alpha bin excess limit; please increase uplimit: " << m_nalpha << " > " << Max_nalpha);}
    for (int i = 0; i < m_nalpha; ++i) {
      proxt >> dumy1 >> dumy2 >> dumy3;
      l_alpha[i] = dumy1;    u_alpha[i] = dumy2;     ialpha[i] = dumy3;
      B2INFO("" << i << " | " << l_alpha[i] << " " << u_alpha[i] << " " << ialpha[i]);
    }
    proxt >> m_ntheta;
    B2INFO("Number of theta bins: " << m_ntheta);
    if (m_ntheta > Max_ntheta) {B2FATAL("number of  theta bin excess limit; please increase uplimit: " << m_ntheta << " > " << Max_ntheta);}
    for (int i = 0; i < m_ntheta; ++i) {
      proxt >> dumy1 >> dumy2 >> dumy3;
      l_theta[i] = dumy1;    u_theta[i] = dumy2;     itheta[i] = dumy3;
      B2INFO("" << i << " |" << l_theta[i] << " " << u_theta[i] << " " << itheta[i]);
    }
  }
  B2INFO("Finish assign sigma bining");
}

readSigma()

void readSigma ( )

protectedvirtual

read sigma from previous calibration, (input sigma)

Definition at line 471 of file SpaceResolutionCalibration.cc.

{
  B2INFO("readSigma");
  if (m_useDB) {
    B2INFO("reading sigma from DB");
    m_sResolFromDB = new DBObjPtr<CDCSpaceResols>;
    readSigmaFromDB();
    B2INFO("Number of theta bins from input sigma: " << ntheta_old);
  } else {
    B2INFO("Read sigma from text");
    readSigmaFromText();
    B2INFO("number of alpha bins from input sigma: " << nalpha_old);
  }
}

readSigmaFromDB()

void readSigmaFromDB ( )

protectedvirtual

read sigma from DB

< angle bin info.

Definition at line 570 of file SpaceResolutionCalibration.cc.

{
  typedef std::array<float, 3> array3; 
  //  std::cout <<"setSResol called" << std::endl;
  nalpha_old = (*m_sResolFromDB)->getNoOfAlphaBins();
  double rad2deg = 180 / M_PI;
  for (unsigned short i = 0; i < nalpha_old; ++i) {
    //    m_alphaPoints[i] = (*dbXT_old)->getAlphaPoint(i);
    array3 alpha = (*m_sResolFromDB)->getAlphaBin(i);
    l_alpha_old[i] = alpha[0] * rad2deg;
    u_alpha_old[i] = alpha[1] * rad2deg;
    ialpha_old[i] = alpha[2] * rad2deg;
    //    std::cout << m_alphaPoints[i]*180./M_PI << std::endl;
  }
 
  ntheta_old = (*m_sResolFromDB)->getNoOfThetaBins();
  B2INFO("Ntheta: " << ntheta_old);
  for (unsigned short i = 0; i < ntheta_old; ++i) {
    //    m_thetaPoints[i] = (*dbXT_old).getThetaPoint(i);
    array3 theta = (*m_sResolFromDB)->getThetaBin(i);
    l_theta_old[i] = theta[0] * rad2deg;
    u_theta_old[i] = theta[1] * rad2deg;
    itheta_old[i] = theta[2] * rad2deg;
  }
  m_sigmaParamMode_old = (*m_sResolFromDB)->getSigmaParamMode();
 
  for (unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
    for (unsigned short iLR = 0; iLR < 2; ++iLR) {
      for (unsigned short iA = 0; iA < nalpha_old; ++iA) {
        for (unsigned short iT = 0; iT < ntheta_old; ++iT) {
          const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
          unsigned short np = params.size();
          //    std::cout <<"np4sigma= " << np << std::endl;
          for (unsigned short i = 0; i < np; ++i) {
            sigma_old[iCL][iLR][iA][iT][i] = params[i];
          }
        }
      }
    }
  }
 
}

readSigmaFromText()

void readSigmaFromText ( )

protectedvirtual

read sigma from text file

Definition at line 486 of file SpaceResolutionCalibration.cc.

{
  ifstream ifs;
  std::string fileName1 = "/data/cdc" + m_sigmafile;
  std::string fileName = FileSystem::findFile(fileName1);
  if (fileName == "") {
    fileName = FileSystem::findFile(m_sigmafile);
  }
  if (fileName == "") {
    cout << "CDCGeometryPar: " << fileName << " not exist!" << endl;
  } else {
    cout << "CDCGeometryPar: " << fileName << " exists." << endl;
    ifs.open(fileName);
    if (!ifs) cout << "CDCGeometryPar: cannot open " << fileName << " !" << endl;
  }
 
  //read alpha bin info.
  if (ifs >> nalpha_old) {
    if (nalpha_old == 0 || nalpha_old > 18) cout << "Fail to read alpha bins !" << endl;
  } else {
    cout << "Fail to read alpha bins !" << endl; return;
  }
  double alpha0, alpha1, alpha2;
  for (unsigned short i = 0; i < nalpha_old; ++i) {
    ifs >> alpha0 >> alpha1 >> alpha2;
    l_alpha_old[i] = alpha0;
    u_alpha_old[i] = alpha1;
    ialpha_old[i] = alpha2;
  }
 
  //read theta bin info.
  if (ifs >> ntheta_old) {
    if (ntheta_old == 0 || ntheta_old > 7) cout << "CDCGeometryPar: fail to read theta bins !" << endl;
  } else {
    cout << "CDCGeometryPar: fail to read theta bins !" << endl;
  }
  double theta0, theta1, theta2;
  for (unsigned short i = 0; i < ntheta_old; ++i) {
    ifs >> theta0 >> theta1 >> theta2;
    l_theta_old[i] = theta0;
    u_theta_old[i] = theta1;
    itheta_old[i] = theta2;
  }
 
  unsigned short np = 0;
  unsigned short iCL, iLR;
  double theta, alpha;
 
  ifs >> m_sigmaParamMode_old >> np;
  double sigma[8]; // cppcheck-suppress constVariable
  //  if (m_sigmaParamMode < 0 || m_sigmaParamMode > 1) cout<<"CDCGeometryPar: invalid sigma-parameterization mode read !"<<endl;
  //if (m_sigmaParamMode == 1) cout<<"CDCGeometryPar: sigma-parameterization mode=1 not ready yet"<<endl;
  //  if (np <= 0 || np > nSigmaParams) cout<<"CDCGeometryPar: no. of sigma-params. outside limits !"<<endl;
  const double epsi = 0.1;
  while (ifs >> iCL) {
    ifs >> theta >> alpha >> iLR;
    for (int i = 0; i < np; ++i) {
      ifs >> sigma[i];
    }
 
    int ith = -99;
    for (unsigned short i = 0; i < ntheta_old; ++i) {
      if (fabs(theta - itheta_old[i]) < epsi) {
        ith = i;
        break;
      }
    }
    if (ith < 0) cout << "CDCGeometryPar: thetas in sigma.dat are inconsistent !" << endl;
 
    int ial = -99;
    for (unsigned short i = 0; i < nalpha_old; ++i) {
      if (fabs(alpha - ialpha_old[i]) < epsi) {
        ial = i;
        break;
      }
    }
    if (ial < 0) cout << "CDCGeometryPar: alphas in sigma.dat are inconsistent !" << endl;
 
    for (int i = 0; i < np; ++i) {
      sigma_old[iCL][iLR][ial][ith][i] = sigma[i];
    }
  }  //end of while loop
  ifs.close();
}

setBinWidth()

virtual void setBinWidth ( double bw )

inlinevirtual

Bin width of each slide.

Definition at line 41 of file SpaceResolutionCalibration.h.

{m_binWidth = bw;}

setDebug()

virtual void setDebug ( bool debug = false )

inlinevirtual

Debug or not.

Definition at line 33 of file SpaceResolutionCalibration.h.

{m_debug = debug; }

setMinimumNDF()

virtual void setMinimumNDF ( double minndf )

inlinevirtual

minimum NDF required for track

Definition at line 37 of file SpaceResolutionCalibration.h.

{m_ndfmin = minndf;}

setMinimumPval()

virtual void setMinimumPval ( double minPval )

inlinevirtual

Minimum Pval required.

Definition at line 39 of file SpaceResolutionCalibration.h.

{m_Pvalmin = minPval;}

setSigmaFileName()

virtual void setSigmaFileName ( std::string name )

inlinevirtual

Output sigma file name, for text mode.

Definition at line 62 of file SpaceResolutionCalibration.h.

{m_sigmafile.assign(name);}

setStoreHisto()

virtual void setStoreHisto ( bool storeHist = false )

inlinevirtual

Store histograms during the calibration or not.

Definition at line 50 of file SpaceResolutionCalibration.h.

{m_storeHisto = storeHist;}

setUseDB()

virtual void setUseDB ( bool useDB = false )

inlinevirtual

Use database or text mode.

Definition at line 35 of file SpaceResolutionCalibration.h.

{m_useDB = useDB; }

storeHisto()

void storeHisto ( )

protectedvirtual

store histogram

Definition at line 372 of file SpaceResolutionCalibration.cc.

{
  B2INFO("storeHisto");
  TFile*  ff = new TFile("sigma_histo.root", "RECREATE");
  TDirectory* top = gDirectory;
  TDirectory* Direct[56];
 
  for (int il = 0; il < 56; ++il) {
    top->cd();
    Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
    Direct[il]->cd();
    for (int lr = 0; lr < 2; ++lr) {
      for (int al = 0; al < m_nalpha; ++al) {
        for (int th = 0; th < m_ntheta; ++th) {
          if (!gr[il][lr][al][th]) continue;
          if (!gfit[il][lr][al][th]) continue;
          if (m_fitflag[il][lr][al][th] == 1) {
            hist_b[il][lr][al][th]->Write();
            hist_u[il][lr][al][th]->Write();
            hb_m[il][lr][al][th]->Write();
            hb_s[il][lr][al][th]->Write();
            hu_m[il][lr][al][th]->Write();
            hu_s[il][lr][al][th]->Write();
            gr[il][lr][al][th]->Write();
            gfit[il][lr][al][th]->Write();
          }
        }
      }
    }
  }
  ff->Close();
  B2INFO("Finish store histogram");
}

useProfileFromInputSigma()

virtual void useProfileFromInputSigma ( bool useProfileFromInputSigma )

inlinevirtual

use sigma bin profile form input sigma or new one from input file

Definition at line 57 of file SpaceResolutionCalibration.h.

      {
        m_useProfileFromInputSigma = useProfileFromInputSigma;
      }

write()

void write ( )

protectedvirtual

save calibration, in text file or db

Definition at line 405 of file SpaceResolutionCalibration.cc.

{
 
  B2INFO("Exporting parameters...");
  int nfitted = 0;
  int nfailure = 0;
  /* Write the fit params*/
 
  ofstream ofs(m_outputSigmaFileName.c_str());
  ofs << m_nalpha << endl;
  for (int i = 0; i < m_nalpha; ++i) {
    ofs << std::setprecision(4) << l_alpha[i] << "   " << std::setprecision(4) << u_alpha[i] << "   " << std::setprecision(
          4) << ialpha[i] << endl;
  }
 
  ofs << m_ntheta << endl;
  for (int i = 0; i < m_ntheta; ++i) {
    ofs << std::setprecision(4) << l_theta[i] << "   " << std::setprecision(4) << u_theta[i] << "   " << std::setprecision(
          4) << itheta[i] << endl;
  }
 
  ofs << 0 << "  " << 7 << endl; //mode and number of params;
  for (int al = 0; al < m_nalpha; ++al) {
    for (int th = 0; th < m_ntheta; ++th) {
      for (int i = 0; i < 56; ++i) {
        for (int lr = 1; lr >= 0; --lr) {
          //          ffit[i][lr][al][th]->GetParameters(par);
          ofs << i << std::setw(4) << itheta[th] << std::setw(4) << ialpha[al] << std::setw(4) << lr << std::setw(15);
          if (m_fitflag[i][lr][al][th] == 1) {
            nfitted += 1;
            for (int p = 0; p < 7; ++p) {
              if (p != 6) { ofs << std::setprecision(7) << sigma_new[i][lr][al][th][p] << std::setw(15);}
              if (p == 6) { ofs << std::setprecision(7) << sigma_new[i][lr][al][th][p] << std::endl;}
            }
          } else {
            B2WARNING("Fitting error and old sigma will be used. (Layer " << i << ") (lr = " << lr << ") (al = " << al << ") (th = " << th <<
                      ")");
            nfailure += 1;
            int ial_old = 0;
            int ith_old = 0;
            for (int k = 0; k < nalpha_old; ++k) {
              if (ialpha[al] < u_alpha_old[k]) {ial_old = k; break;}
            }
            for (int j = 0; j < ntheta_old; ++j) {
              if (itheta[th] < u_theta_old[j]) {ith_old = j; break;}
            }
            for (int p = 0; p < 7; ++p) {
              if (p != 6) { ofs << std::setprecision(7) << sigma_old[i][lr][ial_old][ith_old][p] << std::setw(15);}
              if (p == 6) { ofs << std::setprecision(7) << sigma_old[i][lr][ial_old][ith_old][p] << std::endl;}
            }
          }
        }
      }
    }
  }
  ofs.close();
  B2RESULT("Number of histogram: " << 56 * 2 * m_nalpha * m_ntheta);
  B2RESULT("Histos successfully fitted: " << nfitted);
  B2RESULT("Histos fit failure: " << nfailure);
  if (m_useDB) {
    CDCDatabaseImporter import(0, 0, -1, -1);
    import.importSigma(m_outputSigmaFileName.c_str());
  }
 
 
}

Member Data Documentation

gfit

TGraphErrors* gfit[56][2][18][7]

private

sigma*sigma graph for fit

Definition at line 106 of file SpaceResolutionCalibration.h.

gr

TGraphErrors* gr[56][2][18][7]

private

sigma graph.

Definition at line 107 of file SpaceResolutionCalibration.h.

hb_m

TH1F* hb_m[56][2][Max_nalpha][Max_ntheta]

private

mean histogram of unbiased residual

Definition at line 112 of file SpaceResolutionCalibration.h.

hb_s

TH1F* hb_s[56][2][Max_nalpha][Max_ntheta]

private

sigma histogram of ubiased residual

Definition at line 113 of file SpaceResolutionCalibration.h.

hist_b

TH2F* hist_b[56][2][Max_nalpha][Max_ntheta]

private

2D histogram of biased residual

Definition at line 108 of file SpaceResolutionCalibration.h.

hist_u

TH2F* hist_u[56][2][Max_nalpha][Max_ntheta]

private

2D histogram of unbiased residual

Definition at line 109 of file SpaceResolutionCalibration.h.

hu_m

TH1F* hu_m[56][2][Max_nalpha][Max_ntheta]

private

mean histogram biased residual

Definition at line 110 of file SpaceResolutionCalibration.h.

hu_s

TH1F* hu_s[56][2][Max_nalpha][Max_ntheta]

private

sigma histogram of biased residual

Definition at line 111 of file SpaceResolutionCalibration.h.

ialpha

double ialpha[18]

private

represented alphas of alpha bins.

Definition at line 130 of file SpaceResolutionCalibration.h.

ialpha_old

double ialpha_old[18]

private

represented alphas of alpha bins from input.

Definition at line 140 of file SpaceResolutionCalibration.h.

itheta

double itheta[7]

private

represented alphas of theta bins.

Definition at line 133 of file SpaceResolutionCalibration.h.

itheta_old

double itheta_old[7]

private

represented alphas of theta bins from input.

Definition at line 143 of file SpaceResolutionCalibration.h.

l_alpha

double l_alpha[18]

private

Lower boundaries of alpha bins.

Definition at line 128 of file SpaceResolutionCalibration.h.

l_alpha_old

double l_alpha_old[18]

private

Lower boundaries of alpha bins from input.

Definition at line 138 of file SpaceResolutionCalibration.h.

l_theta

double l_theta[7]

private

Lower boundaries of theta bins.

Definition at line 131 of file SpaceResolutionCalibration.h.

l_theta_old

double l_theta_old[7]

private

Lower boundaries of theta bins from input.

Definition at line 141 of file SpaceResolutionCalibration.h.

m_BField

bool m_BField = true

private

Work with BField, fit range and initial parameters is different incase B and noB.

Definition at line 101 of file SpaceResolutionCalibration.h.

m_binWidth

double m_binWidth = 0.05

private

width of each bin, unit cm

Definition at line 95 of file SpaceResolutionCalibration.h.

m_debug

bool m_debug = false

private

Debug or not.

Definition at line 96 of file SpaceResolutionCalibration.h.

m_fitflag

int m_fitflag[56][2][Max_nalpha][Max_ntheta] = {{{{0}}}}

private

Fit flag; 1:OK ; 0:error.

Definition at line 114 of file SpaceResolutionCalibration.h.

{{{{0}}}} ; 

m_inputRootFileNames

std::string m_inputRootFileNames = "rootfile/output*"

private

Input root file names.

Definition at line 117 of file SpaceResolutionCalibration.h.

m_nalpha

int m_nalpha

private

number of alpha bins

Definition at line 126 of file SpaceResolutionCalibration.h.

m_ndfmin

double m_ndfmin = 5

private

Minimum NDF.

Definition at line 93 of file SpaceResolutionCalibration.h.

m_ntheta

int m_ntheta

private

number of theta bins

Definition at line 127 of file SpaceResolutionCalibration.h.

m_outputSigmaFileName

std::string m_outputSigmaFileName = "sigma_new.dat"

private

Output sigma file name.

Definition at line 116 of file SpaceResolutionCalibration.h.

m_ProfileFileName

std::string m_ProfileFileName = "sigma_profile"

private

Profile file name.

Definition at line 118 of file SpaceResolutionCalibration.h.

m_Pvalmin

double m_Pvalmin = 0.

private

Minimum Prob(chi2) of track.

Definition at line 94 of file SpaceResolutionCalibration.h.

m_sigmafile

std::string m_sigmafile = "cdc/data/sigma.dat"

private

Sigma file name, for text mode.

Definition at line 120 of file SpaceResolutionCalibration.h.

m_sigmaParamMode_old

unsigned short m_sigmaParamMode_old

private

sigma mode from input.

Definition at line 145 of file SpaceResolutionCalibration.h.

m_sResolFromDB

DBObjPtr<CDCSpaceResols>* m_sResolFromDB

private

Database for sigma.

Definition at line 119 of file SpaceResolutionCalibration.h.

m_storeHisto

bool m_storeHisto = false

private

Store histogram or not.

Definition at line 98 of file SpaceResolutionCalibration.h.

m_useDB

bool m_useDB = false

private

use db or text mode

Definition at line 99 of file SpaceResolutionCalibration.h.

m_useProfileFromInputSigma

bool m_useProfileFromInputSigma = true

private

Use binning from old sigma or new one form input.

Definition at line 100 of file SpaceResolutionCalibration.h.

Max_nalpha

const int Max_nalpha = 18

staticprivate

Maximum alpha bin.

Definition at line 89 of file SpaceResolutionCalibration.h.

Max_np

const unsigned short Max_np = 40

staticprivate

Maximum number of point =1/binwidth.

Definition at line 91 of file SpaceResolutionCalibration.h.

Max_ntheta

const int Max_ntheta = 7

staticprivate

maximum theta bin

Definition at line 90 of file SpaceResolutionCalibration.h.

nalpha_old

int nalpha_old

private

number of alpha bins from input

Definition at line 136 of file SpaceResolutionCalibration.h.

ntheta_old

int ntheta_old

private

number of theta bins from input

Definition at line 137 of file SpaceResolutionCalibration.h.

sigma_new

double sigma_new[56][2][18][7][8]

private

new sigma parameters.

Definition at line 104 of file SpaceResolutionCalibration.h.

sigma_old

double sigma_old[56][2][18][7][8]

private

old sigma parameters.

Definition at line 103 of file SpaceResolutionCalibration.h.

u_alpha

double u_alpha[18]

private

Upper boundaries of alpha bins.

Definition at line 129 of file SpaceResolutionCalibration.h.

u_alpha_old

double u_alpha_old[18]

private

Upper boundaries of alpha bins from input.

Definition at line 139 of file SpaceResolutionCalibration.h.

u_theta

double u_theta[7]

private

Upper boundaries of theta bins.

Definition at line 132 of file SpaceResolutionCalibration.h.

u_theta_old

double u_theta_old[7]

private

Upper boundaries of theta bins from input.

Definition at line 142 of file SpaceResolutionCalibration.h.

---

The documentation for this class was generated from the following files:

• cdc/calibration/include/SpaceResolutionCalibration.h
• cdc/calibration/src/SpaceResolutionCalibration.cc