groot/macros/DNL_Correction.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 // PC = PhiCut (efficiency/track attempt before Hough)      2018-12-05 Vlad
0013
0014 TH2D* BH_5Residual_vsPhi_Corrected;
0015 TH1D* ProjectionY;
0016
0017 // Histograming ALL Radii through an array
0018 char name[100];
0019 TH2D* BH_Residual_vsPhi_Corrected_[16];
0020 TH1D* ProjectionY_[16];
0021  TH2D* BH_Residual_vsPhi_PC_Corrected_[16];
0022  TH1D* ProjectionY_PC_[16];
0023
0024
0025 void DNL_Correction(int REBIN=1)
0026 {
0027
0028   BH_5Residual_vsPhi->Rebin2D(1, REBIN);
0029
0030   BH_5Residual_vsPhi_Corrected = (TH2D*)BH_5Residual_vsPhi->Clone("BH_5Residual_vsPhi_Corrected");
0031   BH_5Residual_vsPhi_Corrected->Reset();
0032
0033   for (int i=1; i<15; i++) // 0 & 15 currently aren't "new'ed" in FillBlobHist.C
0034     {
0035
0036       BH_Residual_vsPhi_1-> Rebin2D(1, REBIN);
0037       BH_Residual_vsPhi_2-> Rebin2D(1, REBIN);
0038       BH_Residual_vsPhi_3-> Rebin2D(1, REBIN);
0039       BH_Residual_vsPhi_4-> Rebin2D(1, REBIN);
0040       BH_Residual_vsPhi_5-> Rebin2D(1, REBIN);
0041       BH_Residual_vsPhi_6-> Rebin2D(1, REBIN);
0042       BH_Residual_vsPhi_7-> Rebin2D(1, REBIN);
0043       BH_Residual_vsPhi_8-> Rebin2D(1, REBIN);
0044       BH_Residual_vsPhi_9-> Rebin2D(1, REBIN);
0045       BH_Residual_vsPhi_10-> Rebin2D(1, REBIN);
0046       BH_Residual_vsPhi_11-> Rebin2D(1, REBIN);
0047       BH_Residual_vsPhi_12-> Rebin2D(1, REBIN);
0048       BH_Residual_vsPhi_13-> Rebin2D(1, REBIN);
0049       BH_Residual_vsPhi_14-> Rebin2D(1, REBIN);
0050
0051       //      BH_Residual_vsPhi_PC_[i]->Rebin2D(1, REBIN);
0052
0053
0054       // sprintf(name, "BH_Residual_vsPhi_PC_%d", i);
0055       // name->Rebin2D(1, REBIN);
0056
0057       //BH_Residual_vsPhi_[i] = new TH2D(Form(BH_Residual_vsPhi_%i, i));
0058
0059       sprintf(name, "BH_Residual_vsPhi_Corrected_%d", i);
0060       BH_Residual_vsPhi_Corrected_[i] = (TH2D*)Form("BH_Residual_vsPhi_%d",i)->Clone(name);
0061       BH_Residual_vsPhi_Corrected_[i]->Reset();
0062
0063       sprintf(name, "BH_Residual_vsPhi_PC_Corrected_%d", i);
0064       BH_Residual_vsPhi_PC_Corrected_[i] = (TH2D*)BH_Residual_vsPhi_PC_[i]->Clone(name);
0065       BH_Residual_vsPhi_PC_Corrected_[i]->Reset();
0066     }
0067
0068   for (int i=1; i<BH_5Residual_vsPhi->GetNbinsX()+1; i++)
0069     {
0070       double correction = BH_5Residual_vsPhi->ProfileX()->GetBinContent(i); // This ultimately gives a y-value
0071       double xValue     = BH_5Residual_vsPhi->GetXaxis()->GetBinCenter(i);
0072
0073       for (int j=1; j< BH_5Residual_vsPhi->GetNbinsY()+1; j++)
0074     {
0075       int getBin     = BH_5Residual_vsPhi->GetBin(i,j);
0076       double content = BH_5Residual_vsPhi->GetBinContent(getBin);
0077       double yValue  = BH_5Residual_vsPhi->GetYaxis()->GetBinCenter(j);
0078       BH_5Residual_vsPhi_Corrected->Fill(xValue, yValue-correction, content); // "content" in this 2D histogram is the weight
0079     }
0080     }
0081
0082
0083   for (int k=1; k<15; k++)
0084     {
0085       for (int i=1; i<BH_Residual_vsPhi_[k]->GetNbinsX()+1; i++)
0086     {
0087       double correction = BH_Residual_vsPhi_[k]->ProfileX()->GetBinContent(i); // This ultimately gives a y-value
0088       double xValue     = BH_Residual_vsPhi_[k]->GetXaxis()->GetBinCenter(i);
0089
0090       double correction_PC = BH_Residual_vsPhi_PC_[k]->ProfileX()->GetBinContent(i);
0091       double xValue_PC     = BH_Residual_vsPhi_PC_[k]->GetXaxis()->GetBinCenter(i);
0092
0093       for (int j=1; j< BH_Residual_vsPhi_[k]->GetNbinsY()+1; j++)
0094         {
0095           int getBin     = BH_Residual_vsPhi_[k]->GetBin(i,j);
0096           double content = BH_Residual_vsPhi_[k]->GetBinContent(getBin);
0097           double yValue  = BH_Residual_vsPhi_[k]->GetYaxis()->GetBinCenter(j);
0098           BH_Residual_vsPhi_Corrected_[k]->Fill(xValue, yValue-correction, content); // "content" in this 2D histogram is the weight
0099
0100           int getBin_PC     = BH_Residual_vsPhi_PC_[k]->GetBin(i,j);
0101           double content_PC = BH_Residual_vsPhi_PC_[k]->GetBinContent(getBin_PC);
0102           double yValue_PC  = BH_Residual_vsPhi_PC_[k]->GetYaxis()->GetBinCenter(j);
0103           BH_Residual_vsPhi_PC_Corrected_[k]->Fill(xValue_PC, yValue_PC-correction_PC, content_PC); // "content" in this 2D histogram is the weight
0104         }
0105     }
0106     }
0107
0108   ProjectionY = BH_5Residual_vsPhi_Corrected->ProjectionY(); // Sums (cause that's what histograms do) all x-values along that particular y-value
0109   double min = ProjectionY->GetBinCenter(1); // first bin is minimum (in Gauss tail)
0110   double max = ProjectionY->GetBinCenter(ProjectionY->GetNbinsX());  // last bin is maximum (in Gauss tail)
0111
0112   double min[16], max[16], min_PC[16], max_PC[16];
0113   for (int i=1; i<15; i++)
0114     {
0115       ProjectionY_[i] = BH_Residual_vsPhi_Corrected_[i]->ProjectionY(); // Sums (cause that's what histograms do) all x-values along that particular y-value
0116       min[i] = ProjectionY_[i]->GetBinCenter(1);
0117       max[i] = ProjectionY_[i]->GetBinCenter(ProjectionY_[i]->GetNbinsX());
0118
0119       ProjectionY_PC_[i] = BH_Residual_vsPhi_PC_Corrected_[i]->ProjectionY();
0120       min_PC[i] = ProjectionY_PC_[i]->GetBinCenter(1);
0121       max_PC[i] = ProjectionY_PC_[i]->GetBinCenter(ProjectionY_PC_[i]->GetNbinsX());
0122     }
0123
0124   TF1*   GAUSS = new TF1("GAUSS", "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min, max);
0125   double amp   = ProjectionY->GetBinContent(ProjectionY->GetMaximumBin());
0126   double mean  = ProjectionY->GetMean();
0127   double sigma = ProjectionY->GetRMS();
0128
0129
0130   double amp[16], mean[16], sigma[16];
0131   double amp_PC[16], mean_PC[16], sigma_PC[16];
0132   for (int i=1; i<15; i++)
0133     {
0134       sprintf(name, "GAUSS_%d", i);
0135       TF1* GAUSS_[i] = new TF1(name, "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min[i], max[i]);
0136       amp[i]   = ProjectionY_[i]->GetBinContent(ProjectionY_[i]->GetMaximumBin());
0137       mean[i]  = ProjectionY_[i]->GetMean();
0138       sigma[i] = ProjectionY_[i]->GetRMS();
0139
0140       sprintf(name, "GAUSS_PC_%d", i);
0141       TF1* GAUSS_PC_[i] = new TF1(name, "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min_PC[i], max_PC[i]);
0142       amp_PC[i]   = ProjectionY_PC_[i]->GetBinContent(ProjectionY_PC_[i]->GetMaximumBin());
0143       mean_PC[i]  = ProjectionY_PC_[i]->GetMean();
0144       sigma_PC[i] = ProjectionY_PC_[i]->GetRMS();
0145     }
0146
0147
0148   for (int j=0; j<3; j++)
0149     {
0150       GAUSS->SetParameter(0,amp);
0151       GAUSS->SetParameter(1,mean);
0152       GAUSS->SetParameter(2,sigma);
0153
0154       ProjectionY->Fit(GAUSS, "", "", min, max);
0155
0156       amp   = ProjectionY->GetFunction("GAUSS")->GetParameter(0);
0157       mean  = ProjectionY->GetFunction("GAUSS")->GetParameter(1);
0158       sigma = ProjectionY->GetFunction("GAUSS")->GetParameter(2);
0159
0160       min = mean - 2.0*sigma;
0161       max = mean + 2.0*sigma;
0162     }
0163
0164
0165   for (int i=1; i<15; i++)
0166     {
0167       for (int j=0; j<3; j++)
0168     {
0169       GAUSS_[i]->SetParameter(0,amp[i]);      GAUSS_PC_[i]->SetParameter(0,amp_PC[i]);
0170       GAUSS_[i]->SetParameter(1,mean[i]);     GAUSS_PC_[i]->SetParameter(1,mean_PC[i]);
0171       GAUSS_[i]->SetParameter(2,sigma[i]);    GAUSS_PC_[i]->SetParameter(2,sigma_PC[i]);
0172
0173       ProjectionY_[i]->Fit(GAUSS_[i], "", "", min[i], max[i]);
0174       ProjectionY_PC_[i]->Fit(GAUSS_PC_[i], "", "", min_PC[i], max_PC[i]);
0175
0176       sprintf(name, "GAUSS_%d", i);
0177       amp[i]   = ProjectionY_[i]->GetFunction(name)->GetParameter(0);
0178       mean[i]  = ProjectionY_[i]->GetFunction(name)->GetParameter(1);
0179       sigma[i] = ProjectionY_[i]->GetFunction(name)->GetParameter(2);
0180
0181       sprintf(name, "GAUSS_PC_%d", i);
0182       amp_PC[i]   = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(0);
0183       mean_PC[i]  = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(1);
0184       sigma_PC[i] = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(2);
0185
0186       min[i] = mean[i] - 2.0*sigma[i];    min_PC[i] = mean_PC[i] - 2.0*sigma_PC[i];
0187       max[i] = mean[i] + 2.0*sigma[i];    max_PC[i] = mean_PC[i] + 2.0*sigma_PC[i];
0188     }
0189     }
0190
0191 }