groot/macros/DNL_CorrectionArray.C

0001 /*    This is a decent routine to fix differential non-linearity (DNL)
0002   and  more rapidly get at the actual resolution. We assume the following:
0003   Since the FTBF tracks are "nearly" horizontal, the "phi" coordinate causes
0004   all of sets of triple-layers to smoothly vary across the DNL parameter.
0005       In this case, the residual will have a deterministic dependence on phi.
0006   We'll use a profile histogram from the #Residual_vsPhi result to extract a
0007   correction function. This will be turned into a "shift". Then we will copy
0008   the #Residual_vsPhi into a new histogram while applying the appropriate shifts.
0009
0010                                                             2018-10-30 TKH & Vlad  */
0011
0012
0013 // PC = PhiCut (efficiency/track attempt before Hough)      2018-12-05 Vlad
0014
0015 // New pointers need to be created and assigned to the original histograms in the code
0016 char name[100];
0017 TH2D* BH_Residual_vsPhi_[16];
0018 TH2D* BH_Residual_vsPhi_PC_[16];
0019
0020 // Histograming ALL Radii through an array
0021 TH2D* BH_Residual_vsPhi_Corrected_[16];
0022 TH2D* BH_Residual_vsPhi_PC_Corrected_[16];
0023
0024 void TripleFit(TH2*);
0025
0026 void DNL_CorrectionArray(int REBIN=1)
0027 {
0028
0029   // Assign the array of pointers.
0030   for( int i=0; i<16; i++)
0031     {
0032       sprintf(name, "BH_Residual_vsPhi_%d", i);
0033       _file0->GetObject(name, BH_Residual_vsPhi_[i]);
0034
0035       sprintf(name, "BH_Residual_vsPhi_PC_%d", i);
0036       _file0->GetObject(name, BH_Residual_vsPhi_PC_[i]);
0037     }
0038
0039   // Clone these into the "correct ones"
0040   for (int i=1; i<15; i++) // 0 & 15 currently aren't "new'ed" in FillBlobHist.C
0041     {
0042       BH_Residual_vsPhi_[i]->Rebin2D(1, REBIN);
0043       BH_Residual_vsPhi_PC_[i]->Rebin2D(1, REBIN);
0044
0045       sprintf(name, "BH_Residual_vsPhi_Corrected_%d", i);
0046       BH_Residual_vsPhi_Corrected_[i] = (TH2D*)BH_Residual_vsPhi_[i]->Clone(name);
0047       BH_Residual_vsPhi_Corrected_[i]->Reset();
0048
0049       sprintf(name, "BH_Residual_vsPhi_PC_Corrected_%d", i);
0050       BH_Residual_vsPhi_PC_Corrected_[i] = (TH2D*)BH_Residual_vsPhi_PC_[i]->Clone(name);
0051       BH_Residual_vsPhi_PC_Corrected_[i]->Reset();
0052     }
0053
0054   // Fill the "corrected" histograms
0055   for (int k=1; k<15; k++)
0056     {
0057       cout << "Calculating corrections for layer " << k << endl;
0058       for (int i=1; i<BH_Residual_vsPhi_[k]->GetNbinsX()+1; i++)
0059     {
0060       double correction    = BH_Residual_vsPhi_[k]->ProfileX()->GetBinContent(i); // This ultimately gives a y-value
0061       double correction_PC = BH_Residual_vsPhi_PC_[k]->ProfileX()->GetBinContent(i);
0062
0063       double xValue    = BH_Residual_vsPhi_[k]->GetXaxis()->GetBinCenter(i);
0064       double xValue_PC = BH_Residual_vsPhi_PC_[k]->GetXaxis()->GetBinCenter(i);
0065
0066       for (int j=1; j< BH_Residual_vsPhi_[k]->GetNbinsY()+1; j++)
0067         {
0068           int getBin;
0069           double content, yValue;
0070
0071           getBin  = BH_Residual_vsPhi_[k]->GetBin(i,j);
0072           content = BH_Residual_vsPhi_[k]->GetBinContent(getBin);
0073           yValue  = BH_Residual_vsPhi_[k]->GetYaxis()->GetBinCenter(j);
0074           BH_Residual_vsPhi_Corrected_[k]->Fill(xValue, yValue-correction, content); // "content" in this 2D histogram is the weight
0075
0076           getBin  = BH_Residual_vsPhi_PC_[k]->GetBin(i,j);
0077           content = BH_Residual_vsPhi_PC_[k]->GetBinContent(getBin);
0078           yValue  = BH_Residual_vsPhi_PC_[k]->GetYaxis()->GetBinCenter(j);
0079           BH_Residual_vsPhi_PC_Corrected_[k]->Fill(xValue_PC, yValue-correction_PC, content);
0080         }
0081     }
0082     }
0083
0084   for (int i=1; i<15; i++)
0085     {
0086       TripleFit(BH_Residual_vsPhi_Corrected_[i]);
0087       TripleFit(BH_Residual_vsPhi_PC_Corrected_[i]);
0088     }
0089 }
0090
0091 TF1* GAUSS = 0;
0092
0093 void TripleFit(TH2* hist)
0094 {
0095
0096   if (hist->Integral()<1000) { return; }
0097
0098   TH1D* ProjectionY = hist->ProjectionY(); // Sums (cause that's what histograms do) all x-values along that particular y-value
0099   double min = ProjectionY->GetBinCenter(1); // first bin is minimum (in Gauss tail)
0100   double max = ProjectionY->GetBinCenter(ProjectionY->GetNbinsX());  // last bin is maximum (in Gauss tail)
0101
0102   double amp   = ProjectionY->GetBinContent(ProjectionY->GetMaximumBin());
0103   double mean  = ProjectionY->GetMean();
0104   double sigma = ProjectionY->GetRMS();
0105
0106   if (!GAUSS) { GAUSS = new TF1("GAUSS", "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min, max); }
0107
0108   for (int j=0; j<3; j++)
0109     {
0110       GAUSS->SetParameter(0,amp);
0111       GAUSS->SetParameter(1,mean);
0112       GAUSS->SetParameter(2,sigma);
0113
0114       ProjectionY->Fit(GAUSS, "Q0", "", min, max);
0115
0116       amp   = ProjectionY->GetFunction("GAUSS")->GetParameter(0);
0117       mean  = ProjectionY->GetFunction("GAUSS")->GetParameter(1);
0118       sigma = ProjectionY->GetFunction("GAUSS")->GetParameter(2);
0119
0120       min = mean - 2.0*sigma;
0121       max = mean + 2.0*sigma;
0122     }
0123
0124
0125   cout << "Name: " << hist->GetName();
0126   cout << "\t Amp:" << amp;
0127   cout << "\t Mean:" << mean;
0128   cout << "\t Sigma:" << sigma;
0129   cout << endl;
0130 }
0131