Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DNL_Correction.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DNL_Correction.C
1 /* This is a decent routine to fix differential non-linearity (DNL)
2  and more rapidly get at the actual resolution. We assume the following:
3  Since the FTBF tracks are "nearly" horizontal, the "phi" coordinate causes
4  all of sets of triple-layers to smoothly vary across the DNL parameter.
5  In this case, the residual will have a deterministic dependence on phi.
6  We'll use a profile histogram from the #Residual_vsPhi result to extract a
7  correction function. This will be turned into a "shift". Then we will copy
8  the #Residual_vsPhi into a new histogram while applying the appropriate shifts.
9 
10  2018-10-30 TKH & Vlad */
11 
12 // PC = PhiCut (efficiency/track attempt before Hough) 2018-12-05 Vlad
13 
14 TH2D* BH_5Residual_vsPhi_Corrected;
15 TH1D* ProjectionY;
16 
17 // Histograming ALL Radii through an array
18 char name[100];
19 TH2D* BH_Residual_vsPhi_Corrected_[16];
20 TH1D* ProjectionY_[16];
21  TH2D* BH_Residual_vsPhi_PC_Corrected_[16];
22  TH1D* ProjectionY_PC_[16];
23 
24 
25 void DNL_Correction(int REBIN=1)
26 {
27 
28  BH_5Residual_vsPhi->Rebin2D(1, REBIN);
29 
30  BH_5Residual_vsPhi_Corrected = (TH2D*)BH_5Residual_vsPhi->Clone("BH_5Residual_vsPhi_Corrected");
31  BH_5Residual_vsPhi_Corrected->Reset();
32 
33  for (int i=1; i<15; i++) // 0 & 15 currently aren't "new'ed" in FillBlobHist.C
34  {
35 
36  BH_Residual_vsPhi_1-> Rebin2D(1, REBIN);
37  BH_Residual_vsPhi_2-> Rebin2D(1, REBIN);
38  BH_Residual_vsPhi_3-> Rebin2D(1, REBIN);
39  BH_Residual_vsPhi_4-> Rebin2D(1, REBIN);
40  BH_Residual_vsPhi_5-> Rebin2D(1, REBIN);
41  BH_Residual_vsPhi_6-> Rebin2D(1, REBIN);
42  BH_Residual_vsPhi_7-> Rebin2D(1, REBIN);
43  BH_Residual_vsPhi_8-> Rebin2D(1, REBIN);
44  BH_Residual_vsPhi_9-> Rebin2D(1, REBIN);
45  BH_Residual_vsPhi_10-> Rebin2D(1, REBIN);
46  BH_Residual_vsPhi_11-> Rebin2D(1, REBIN);
47  BH_Residual_vsPhi_12-> Rebin2D(1, REBIN);
48  BH_Residual_vsPhi_13-> Rebin2D(1, REBIN);
49  BH_Residual_vsPhi_14-> Rebin2D(1, REBIN);
50 
51  // BH_Residual_vsPhi_PC_[i]->Rebin2D(1, REBIN);
52 
53 
54  // sprintf(name, "BH_Residual_vsPhi_PC_%d", i);
55  // name->Rebin2D(1, REBIN);
56 
57  //BH_Residual_vsPhi_[i] = new TH2D(Form(BH_Residual_vsPhi_%i, i));
58 
59  sprintf(name, "BH_Residual_vsPhi_Corrected_%d", i);
60  BH_Residual_vsPhi_Corrected_[i] = (TH2D*)Form("BH_Residual_vsPhi_%d",i)->Clone(name);
61  BH_Residual_vsPhi_Corrected_[i]->Reset();
62 
63  sprintf(name, "BH_Residual_vsPhi_PC_Corrected_%d", i);
64  BH_Residual_vsPhi_PC_Corrected_[i] = (TH2D*)BH_Residual_vsPhi_PC_[i]->Clone(name);
65  BH_Residual_vsPhi_PC_Corrected_[i]->Reset();
66  }
67 
68  for (int i=1; i<BH_5Residual_vsPhi->GetNbinsX()+1; i++)
69  {
70  double correction = BH_5Residual_vsPhi->ProfileX()->GetBinContent(i); // This ultimately gives a y-value
71  double xValue = BH_5Residual_vsPhi->GetXaxis()->GetBinCenter(i);
72 
73  for (int j=1; j< BH_5Residual_vsPhi->GetNbinsY()+1; j++)
74  {
75  int getBin = BH_5Residual_vsPhi->GetBin(i,j);
76  double content = BH_5Residual_vsPhi->GetBinContent(getBin);
77  double yValue = BH_5Residual_vsPhi->GetYaxis()->GetBinCenter(j);
78  BH_5Residual_vsPhi_Corrected->Fill(xValue, yValue-correction, content); // "content" in this 2D histogram is the weight
79  }
80  }
81 
82 
83  for (int k=1; k<15; k++)
84  {
85  for (int i=1; i<BH_Residual_vsPhi_[k]->GetNbinsX()+1; i++)
86  {
87  double correction = BH_Residual_vsPhi_[k]->ProfileX()->GetBinContent(i); // This ultimately gives a y-value
88  double xValue = BH_Residual_vsPhi_[k]->GetXaxis()->GetBinCenter(i);
89 
90  double correction_PC = BH_Residual_vsPhi_PC_[k]->ProfileX()->GetBinContent(i);
91  double xValue_PC = BH_Residual_vsPhi_PC_[k]->GetXaxis()->GetBinCenter(i);
92 
93  for (int j=1; j< BH_Residual_vsPhi_[k]->GetNbinsY()+1; j++)
94  {
95  int getBin = BH_Residual_vsPhi_[k]->GetBin(i,j);
96  double content = BH_Residual_vsPhi_[k]->GetBinContent(getBin);
97  double yValue = BH_Residual_vsPhi_[k]->GetYaxis()->GetBinCenter(j);
98  BH_Residual_vsPhi_Corrected_[k]->Fill(xValue, yValue-correction, content); // "content" in this 2D histogram is the weight
99 
100  int getBin_PC = BH_Residual_vsPhi_PC_[k]->GetBin(i,j);
101  double content_PC = BH_Residual_vsPhi_PC_[k]->GetBinContent(getBin_PC);
102  double yValue_PC = BH_Residual_vsPhi_PC_[k]->GetYaxis()->GetBinCenter(j);
103  BH_Residual_vsPhi_PC_Corrected_[k]->Fill(xValue_PC, yValue_PC-correction_PC, content_PC); // "content" in this 2D histogram is the weight
104  }
105  }
106  }
107 
108  ProjectionY = BH_5Residual_vsPhi_Corrected->ProjectionY(); // Sums (cause that's what histograms do) all x-values along that particular y-value
109  double min = ProjectionY->GetBinCenter(1); // first bin is minimum (in Gauss tail)
110  double max = ProjectionY->GetBinCenter(ProjectionY->GetNbinsX()); // last bin is maximum (in Gauss tail)
111 
112  double min[16], max[16], min_PC[16], max_PC[16];
113  for (int i=1; i<15; i++)
114  {
115  ProjectionY_[i] = BH_Residual_vsPhi_Corrected_[i]->ProjectionY(); // Sums (cause that's what histograms do) all x-values along that particular y-value
116  min[i] = ProjectionY_[i]->GetBinCenter(1);
117  max[i] = ProjectionY_[i]->GetBinCenter(ProjectionY_[i]->GetNbinsX());
118 
119  ProjectionY_PC_[i] = BH_Residual_vsPhi_PC_Corrected_[i]->ProjectionY();
120  min_PC[i] = ProjectionY_PC_[i]->GetBinCenter(1);
121  max_PC[i] = ProjectionY_PC_[i]->GetBinCenter(ProjectionY_PC_[i]->GetNbinsX());
122  }
123 
124  TF1* GAUSS = new TF1("GAUSS", "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min, max);
125  double amp = ProjectionY->GetBinContent(ProjectionY->GetMaximumBin());
126  double mean = ProjectionY->GetMean();
127  double sigma = ProjectionY->GetRMS();
128 
129 
130  double amp[16], mean[16], sigma[16];
131  double amp_PC[16], mean_PC[16], sigma_PC[16];
132  for (int i=1; i<15; i++)
133  {
134  sprintf(name, "GAUSS_%d", i);
135  TF1* GAUSS_[i] = new TF1(name, "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min[i], max[i]);
136  amp[i] = ProjectionY_[i]->GetBinContent(ProjectionY_[i]->GetMaximumBin());
137  mean[i] = ProjectionY_[i]->GetMean();
138  sigma[i] = ProjectionY_[i]->GetRMS();
139 
140  sprintf(name, "GAUSS_PC_%d", i);
141  TF1* GAUSS_PC_[i] = new TF1(name, "[0]*exp(-(x-[1])*(x-[1])/(2.0*[2]*[2]))", min_PC[i], max_PC[i]);
142  amp_PC[i] = ProjectionY_PC_[i]->GetBinContent(ProjectionY_PC_[i]->GetMaximumBin());
143  mean_PC[i] = ProjectionY_PC_[i]->GetMean();
144  sigma_PC[i] = ProjectionY_PC_[i]->GetRMS();
145  }
146 
147 
148  for (int j=0; j<3; j++)
149  {
150  GAUSS->SetParameter(0,amp);
151  GAUSS->SetParameter(1,mean);
152  GAUSS->SetParameter(2,sigma);
153 
154  ProjectionY->Fit(GAUSS, "", "", min, max);
155 
156  amp = ProjectionY->GetFunction("GAUSS")->GetParameter(0);
157  mean = ProjectionY->GetFunction("GAUSS")->GetParameter(1);
158  sigma = ProjectionY->GetFunction("GAUSS")->GetParameter(2);
159 
160  min = mean - 2.0*sigma;
161  max = mean + 2.0*sigma;
162  }
163 
164 
165  for (int i=1; i<15; i++)
166  {
167  for (int j=0; j<3; j++)
168  {
169  GAUSS_[i]->SetParameter(0,amp[i]); GAUSS_PC_[i]->SetParameter(0,amp_PC[i]);
170  GAUSS_[i]->SetParameter(1,mean[i]); GAUSS_PC_[i]->SetParameter(1,mean_PC[i]);
171  GAUSS_[i]->SetParameter(2,sigma[i]); GAUSS_PC_[i]->SetParameter(2,sigma_PC[i]);
172 
173  ProjectionY_[i]->Fit(GAUSS_[i], "", "", min[i], max[i]);
174  ProjectionY_PC_[i]->Fit(GAUSS_PC_[i], "", "", min_PC[i], max_PC[i]);
175 
176  sprintf(name, "GAUSS_%d", i);
177  amp[i] = ProjectionY_[i]->GetFunction(name)->GetParameter(0);
178  mean[i] = ProjectionY_[i]->GetFunction(name)->GetParameter(1);
179  sigma[i] = ProjectionY_[i]->GetFunction(name)->GetParameter(2);
180 
181  sprintf(name, "GAUSS_PC_%d", i);
182  amp_PC[i] = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(0);
183  mean_PC[i] = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(1);
184  sigma_PC[i] = ProjectionY_PC_[i]->GetFunction(name)->GetParameter(2);
185 
186  min[i] = mean[i] - 2.0*sigma[i]; min_PC[i] = mean_PC[i] - 2.0*sigma_PC[i];
187  max[i] = mean[i] + 2.0*sigma[i]; max_PC[i] = mean_PC[i] + 2.0*sigma_PC[i];
188  }
189  }
190 
191 }