CMS 3D CMS Logo

Resolutions.cc
Go to the documentation of this file.
1 using namespace ROOT;
2 using ROOT::RDF::RNode;
3 using floats = ROOT::VecOps::RVec<float>;
4 using ints = ROOT::VecOps::RVec<int>;
5 using bools = ROOT::VecOps::RVec<bool>;
6 using chars = ROOT::VecOps::RVec<UChar_t>;
7 using doubles = ROOT::VecOps::RVec<double>;
8 
9 vector<float> HitResolutionVector;
10 vector<float> DoubleDifferenceVector;
11 vector<float> HitDXVector;
12 vector<float> TrackDXVector;
13 vector<float> TrackDXEVector;
14 
18 
19 void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& UL) {
20  std::string CutFlowReportString;
21  std::string DoubleDiffString;
22  std::string HitDXString;
23  std::string TrackDXString;
24  std::string TrackDXEString;
25  std::string ClusterW1String = "clusterW1";
26  std::string ClusterW2String = "clusterW2";
27 
28  switch (UL) {
29  case 0:
30  switch (Unit_Int) {
31  case 0:
32  GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO.root";
33  HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO.txt";
34  CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO.txt";
35  DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch";
36  HitDXString = "hitDX_OverPitch";
37  TrackDXString = "trackDX_OverPitch";
38  TrackDXEString = "trackDXE_OverPitch";
39  break;
40 
41  case 1:
42  GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO.root";
43  HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO.txt";
44  CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO.txt";
45  DoubleDiffString = "hitDX-trackDX";
46  HitDXString = "hitDX";
47  TrackDXString = "trackDX";
48  TrackDXEString = "trackDXE";
49  break;
50 
51  default:
52  std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl;
53  break;
54  }
55 
56  InputFileString = "hitresol_ALCARECO.root";
57  break;
58 
59  case 1:
60  switch (Unit_Int) {
61  case 0:
62  GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO_UL.root";
63  HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO_UL.txt";
64  CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO_UL.txt";
65  DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch";
66  HitDXString = "hitDX_OverPitch";
67  TrackDXString = "trackDX_OverPitch";
68  TrackDXEString = "trackDXE_OverPitch";
69  break;
70 
71  case 1:
72  GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO_UL.root";
73  HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO_UL.txt";
74  CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO_UL.txt";
75  DoubleDiffString = "hitDX-trackDX";
76  HitDXString = "hitDX";
77  TrackDXString = "trackDX";
78  TrackDXEString = "trackDXE";
79  break;
80 
81  default:
82  std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl;
83  break;
84  }
85 
86  InputFileString = "hitresol_ALCARECO_UL.root";
87  break;
88  default:
89  std::cout << "The UL input parameter must be set to 0 (for ALCARECO) or 1 (for UL ALCARECO)." << std::endl;
90  break;
91  }
92 
93  //opening the root file
94  ROOT::RDataFrame d("anResol/reso", InputFileString);
95 
96  int RegionInt = 0;
97 
98  if (region == "TIB_L1") {
99  RegionInt = 1;
100  } else if (region == "TIB_L2") {
101  RegionInt = 2;
102  } else if (region == "TIB_L3") {
103  RegionInt = 3;
104  } else if (region == "TIB_L4") {
105  RegionInt = 4;
106  } else if (region == "Side_TID") {
107  RegionInt = 5;
108  } else if (region == "Wheel_TID") {
109  RegionInt = 6;
110  } else if (region == "Ring_TID") {
111  RegionInt = 7;
112  } else if (region == "TOB_L1") {
113  RegionInt = 8;
114  } else if (region == "TOB_L2") {
115  RegionInt = 9;
116  } else if (region == "TOB_L3") {
117  RegionInt = 10;
118  } else if (region == "TOB_L4") {
119  RegionInt = 11;
120  } else if (region == "TOB_L5") {
121  RegionInt = 12;
122  } else if (region == "TOB_L6") {
123  RegionInt = 13;
124  } else if (region == "Side_TEC") {
125  RegionInt = 14;
126  } else if (region == "Wheel_TEC") {
127  RegionInt = 15;
128  } else if (region == "Ring_TEC") {
129  RegionInt = 16;
130  } else if (region == "TIB_All") {
131  RegionInt = 17;
132  } else if (region == "TOB_All") {
133  RegionInt = 18;
134  } else if (region == "TID_All") {
135  RegionInt = 19;
136  } else if (region == "TEC_All") {
137  RegionInt = 20;
138  } else if (region == "Pixel_Barrel") {
139  RegionInt = 21;
140  } else if (region == "Pixel_EndcapDisk") {
141  RegionInt = 22;
142  } else {
143  std::cout << "Error: The tracker region " << region
144  << " was chosen. Please choose a region out of: TIB L1, TIB L2, TIB L3, TIB L4, Side TID, Wheel TID, "
145  "Ring TID, TOB L1, TOB L2, TOB L3, TOB L4, TOB L5, TOB L6, Side TEC, Wheel TEC or Ring TEC."
146  << std::endl;
147  return 0;
148  }
149 
150  //Lambda function to filter the detID for different layers
151  auto SubDet_Function{[&RegionInt](const int& detID1_input, const int& detID2_input) {
152  bool OutputBool = 0;
153 
154  switch (RegionInt) {
155  case 1: {
156  OutputBool = (((detID1_input >> 25) & 0x7) == 3) && ((detID1_input >> 14) & 0x7) == 1 &&
157  (((detID2_input >> 25) & 0x7) == 3) && ((detID2_input >> 14) & 0x7) == 1; //TIB L1
158  break;
159  }
160 
161  case 2: {
162  OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 2) &&
163  (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 2); //TIB L2
164  break;
165  }
166 
167  case 3: {
168  OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 3) &&
169  (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 3); //TIB L3
170  break;
171  }
172 
173  case 4: {
174  OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 4) &&
175  (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 4); //TIB L4
176  break;
177  }
178 
179  case 5: {
180  OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x3) == 1)) ||
181  ((((detID1_input >> 13) & 0x3) == 2) &&
182  (((detID2_input >> 13) & 0x3) == 2)); //TID Side (1 -> TID-, 2 -> TID+)
183 
184  break;
185  }
186 
187  case 6: {
188  OutputBool = (((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2); //TID Wheel
189 
190  break;
191  }
192 
193  case 7: {
194  OutputBool = ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2)); //TID Ring
195 
196  break;
197  }
198 
199  case 8: {
200  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 1) &&
201  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 1); //TOB L1
202  break;
203  }
204 
205  case 9: {
206  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 2) &&
207  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 2); //TOB L2
208  break;
209  }
210 
211  case 10: {
212  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 3) &&
213  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 3); //TOB L3
214  break;
215  }
216 
217  case 11: {
218  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 4) &&
219  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 4); //TOB L4
220  break;
221  }
222 
223  case 12: {
224  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 5) &&
225  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 5); //TOB L5
226  break;
227  }
228 
229  case 13: {
230  OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 6) &&
231  (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 6); //TOB L6
232  break;
233  }
234 
235  case 14: {
236  OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) ||
237  ((((detID1_input >> 18) & 0x3) == 2) &&
238  (((detID2_input >> 18) & 0x3) == 2)); //Side TEC (1 -> back, 2 -> front)
239  break;
240  }
241 
242  case 15: {
243  OutputBool = (((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4); //Wheel TEC
244  break;
245  }
246 
247  case 16: {
248  OutputBool = (((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3); //Ring TEC
249 
250  break;
251  }
252 
253  case 17: {
254  OutputBool = ((((detID1_input >> 25) & 0x7) == 3) && (((detID2_input >> 25) & 0x7) == 3)); //All TIB
255 
256  break;
257  }
258 
259  case 18: {
260  OutputBool = ((((detID1_input >> 25) & 0x7) == 5) && (((detID2_input >> 25) & 0x7) == 5)); //All TOB
261 
262  break;
263  }
264 
265  case 19: {
266  OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x7) == 1)) ||
267  ((((detID1_input >> 13) & 0x3) == 2) && (((detID2_input >> 13) & 0x7) == 2)) ||
268  ((((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2)) ||
269  ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2)) ||
270  ((((detID1_input >> 7) & 0x3) == 1) && (((detID2_input >> 7) & 0x3) == 1)) ||
271  ((((detID1_input >> 7) & 0x3) == 2) && (((detID2_input >> 7) & 0x3) == 2)) ||
272  ((((detID1_input >> 2) & 0x1F) == 5) && (((detID2_input >> 2) & 0x1F) == 5)) ||
273  ((((detID1_input >> 0) & 0x3) == 0) && (((detID2_input >> 0) & 0x3) == 0)) ||
274  ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) ||
275  ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2)); //All TID
276 
277  break;
278  }
279 
280  case 20: {
281  OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) ||
282  ((((detID1_input >> 18) & 0x3) == 2) && (((detID2_input >> 18) & 0x3) == 2)) ||
283  ((((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4)) ||
284  ((((detID1_input >> 12) & 0x3) == 1) && (((detID2_input >> 12) & 0x3) == 1)) ||
285  ((((detID1_input >> 12) & 0x3) == 2) && (((detID2_input >> 12) & 0x3) == 2)) ||
286  ((((detID1_input >> 8) & 0xF) == 4) && (((detID2_input >> 8) & 0xF) == 4)) ||
287  ((((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3)) ||
288  ((((detID1_input >> 2) & 0x7) == 3) && (((detID2_input >> 2) & 0x7) == 3)) ||
289  ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) ||
290  ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2)) ||
291  ((((detID1_input >> 0) & 0x3) == 3) && (((detID2_input >> 0) & 0x3) == 3)); //All TEC
292 
293  break;
294  }
295 
296  case 21: {
297  OutputBool =
298  (((detID1_input >> 20) & 0xF) == 4) && (((detID2_input >> 20) & 0xF) == 4); //pixel barrel (phase 1)
299  break;
300  }
301 
302  case 22: {
303  OutputBool =
304  (((detID1_input >> 18) & 0xF) == 4) && (((detID2_input >> 18) & 0xF) == 4); //pixel endcap disk (phase 1)
305  break;
306  }
307  }
308 
309  return OutputBool;
310  }};
311 
312  //Function for expressing the hit resolution in either micrometres or pitch units.
313  auto Pitch_Function{[&Unit_Int](const float& pitch, const float& input) {
314  float InputOverPitch = input / pitch;
315  return InputOverPitch;
316  }};
317 
318  //Defining columns needed for the unit conversion into pitch units, and applying the filter for the subdetector
319  auto dataframe = d.Define("hitDX_OverPitch", Pitch_Function, {"pitch1", "hitDX"})
320  .Define("trackDX_OverPitch", Pitch_Function, {"pitch1", "trackDX"})
321  .Define("trackDXE_OverPitch", Pitch_Function, {"pitch1", "trackDXE"})
322  .Filter(SubDet_Function, {"detID1", "detID2"}, "Subdetector filter");
323 
324  //Implementing selection criteria that were not implemented in HitResol.cc
325  auto PairPathCriteriaFunction{[&RegionInt](const float& pairPath_input) {
326  if ((RegionInt > 0 && RegionInt < 5) || (RegionInt > 7 || RegionInt < 13) || (RegionInt == 17) ||
327  (RegionInt == 18)) {
328  return abs(pairPath_input) < 7;
329  } //for TIB and TOB
330  else if (RegionInt == 21 || RegionInt == 22) {
331  return abs(pairPath_input) < 2;
332  } //for pixels
333  else {
334  return abs(pairPath_input) < 20;
335  } //for everything else (max value is 15cm so this will return all events anyway)
336  }};
337 
338  auto MomentaFunction{[&RegionInt](const float& momentum_input) {
339  if (RegionInt == 21 || RegionInt == 22) {
340  return momentum_input > 5;
341  } //pixels
342  else {
343  return momentum_input > 15;
344  } //strips
345  }};
346 
347  auto dataframe_filtered =
348  dataframe.Filter(PairPathCriteriaFunction, {"pairPath"}, "Pair path criterion filter")
349  .Filter(MomentaFunction, {"momentum"}, "Momentum criterion filter")
350  .Filter("trackChi2 > 0.001", "chi2 criterion filter")
351  .Filter("numHits > 6", "numHits filter")
352  .Filter("trackDXE < 0.0025", "trackDXE filter")
353  .Filter("(clusterW1 == clusterW2) && clusterW1 <= 4 && clusterW2 <= 4", "cluster filter");
354 
355  //Creating histograms for the difference between the two hit positions, the difference between the two predicted positions and for the double difference
356  //hitDX = the difference in the hit positions for the pair
357  //trackDX = the difference in the track positions for the pair
358 
359  auto HistoName_DoubleDiff = "DoubleDifference_" + region;
360  auto HistoName_HitDX = "HitDX_" + region;
361  auto HistoName_TrackDX = "TrackDX_" + region;
362  auto HistoName_TrackDXE = "TrackDXE_" + region;
363  auto HistoName_ClusterW1 = "ClusterW1_" + region;
364  auto HistoName_ClusterW2 = "ClusterW2_" + region;
365 
366  auto h_DoubleDifference =
367  dataframe_filtered.Define(HistoName_DoubleDiff, DoubleDiffString)
368  .Histo1D({HistoName_DoubleDiff.c_str(), HistoName_DoubleDiff.c_str(), 40, -0.5, 0.5}, HistoName_DoubleDiff);
369  auto h_hitDX = dataframe_filtered.Define(HistoName_HitDX, HitDXString).Histo1D(HistoName_HitDX);
370  auto h_trackDX = dataframe_filtered.Define(HistoName_TrackDX, TrackDXString).Histo1D(HistoName_TrackDX);
371  auto h_trackDXE = dataframe_filtered.Define(HistoName_TrackDXE, TrackDXEString).Histo1D(HistoName_TrackDXE);
372 
373  auto h_clusterW1 = dataframe_filtered.Define(HistoName_ClusterW1, ClusterW1String).Histo1D(HistoName_ClusterW1);
374  auto h_clusterW2 = dataframe_filtered.Define(HistoName_ClusterW2, ClusterW2String).Histo1D(HistoName_ClusterW2);
375 
376  //Applying gaussian fits, taking the resolutions and squaring them
377  h_DoubleDifference->Fit("gaus");
378 
379  auto double_diff_StdDev = h_DoubleDifference->GetStdDev();
380  auto hitDX_StdDev = h_hitDX->GetStdDev();
381  auto trackDX_StdDev = h_trackDX->GetStdDev();
382  auto trackDXE_Mean = h_trackDXE->GetMean();
383 
384  auto sigma2_MeasMinusPred = pow(double_diff_StdDev, 2);
385  auto sigma2_Meas = pow(hitDX_StdDev, 2);
386  auto sigma2_Pred = pow(trackDX_StdDev, 2);
387  auto sigma2_PredError = pow(trackDXE_Mean, 2);
388 
389  DoubleDifferenceVector.push_back(sigma2_MeasMinusPred);
390  HitDXVector.push_back(sigma2_Meas);
391  TrackDXVector.push_back(sigma2_Pred);
392  TrackDXEVector.push_back(sigma2_PredError);
393 
394  //Saving the histograms with gaussian fits applied to an output root file
395  TFile* output = new TFile(GaussianFitsFileName.c_str(), "UPDATE");
396 
397  h_DoubleDifference->Write();
398  h_hitDX->Write();
399  h_trackDX->Write();
400  h_trackDXE->Write();
401  h_clusterW1->Write();
402  h_clusterW2->Write();
403 
404  output->Close();
405 
406  //Calculating the hit resolution;
407  auto numerator = sigma2_MeasMinusPred - sigma2_PredError;
408 
409  auto HitResolution = sqrt(numerator / 2);
410  HitResolutionVector.push_back(HitResolution);
411 
412  //Printing the resolution
413  std::cout << '\n' << std::endl;
414  std::cout << "The hit resolution for tracker region " << region << " is: " << HitResolution << std::endl;
415  std::cout << '\n' << std::endl;
416 
417  //Cut flow report
418  auto allCutsReport = d.Report();
419  std::ofstream CutFlowReport;
420 
421  CutFlowReport.open(CutFlowReportString.c_str());
422 
423  for (auto&& cutInfo : allCutsReport) {
424  CutFlowReport << cutInfo.GetName() << '\t' << cutInfo.GetAll() << '\t' << cutInfo.GetPass() << '\t'
425  << cutInfo.GetEff() << " %" << std::endl;
426  }
427 }
428 
429 void Resolutions() {
430  int UnitInteger = 0;
431  int ULInteger = 0;
432 
433  vector<std::string> LayerNames = {"TIB_L1", "TIB_L2", "TIB_L3", "TIB_L4", "Side_TID", "Wheel_TID",
434  "Ring_TID", "TOB_L1", "TOB_L2", "TOB_L3", "TOB_L4", "TOB_L5",
435  "TOB_L6", "Side_TEC", "Wheel_TEC", "Ring_TEC", "TIB_All", "TOB_All",
436  "TID_All", "TEC_All", "Pixel_Barrel", "Pixel_EndcapDisk"};
437 
438  for (int i = 0; i < LayerNames.size(); i++) {
439  ResolutionsCalculator(LayerNames.at(i), UnitInteger, ULInteger);
440  }
441 
442  std::ofstream HitResoTextFile;
443  HitResoTextFile.open(HitResoFileName);
444 
445  auto Width = 28;
446 
447  HitResoTextFile << std::right << "Layer " << std::setw(Width) << " Resolution " << std::setw(Width)
448  << " sigma2_HitDX " << std::setw(Width) << " sigma2_trackDX " << std::setw(Width)
449  << " sigma2_trackDXE " << std::setw(Width) << " sigma2_DoubleDifference " << std::endl;
450 
451  for (int i = 0; i < HitResolutionVector.size(); i++) {
452  HitResoTextFile << std::right << LayerNames.at(i) << std::setw(Width) << HitResolutionVector.at(i)
453  << std::setw(Width) << HitDXVector.at(i) << std::setw(Width) << TrackDXVector.at(i)
454  << std::setw(Width) << TrackDXEVector.at(i) << std::setw(Width) << DoubleDifferenceVector.at(i)
455  << std::endl;
456  }
457 
458  system(
459  "mkdir HitResolutionValues; mkdir GaussianFits; mkdir CutFlowReports; mv CutFlowReport_* CutFlowReports/; mv "
460  "HitResolutionValues_* HitResolutionValues/; mv GaussianFits_* GaussianFits/;");
461 }
constexpr int pow(int x)
Definition: conifer.h:24
ROOT::VecOps::RVec< double > doubles
Definition: Resolutions.cc:7
static std::string const input
Definition: EdmProvDump.cc:50
std::string HitResoFileName
Definition: Resolutions.cc:16
ROOT::VecOps::RVec< float > floats
Definition: Resolutions.cc:3
ROOT::VecOps::RVec< UChar_t > chars
Definition: Resolutions.cc:6
void ResolutionsCalculator(const string &region, const int &Unit_Int, const int &UL)
Definition: Resolutions.cc:19
std::string GaussianFitsFileName
Definition: Resolutions.cc:17
T sqrt(T t)
Definition: SSEVec.h:19
dd4hep::Filter Filter
vector< float > TrackDXVector
Definition: Resolutions.cc:12
vector< float > HitResolutionVector
Definition: Resolutions.cc:9
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void Resolutions()
Definition: Resolutions.cc:429
vector< float > DoubleDifferenceVector
Definition: Resolutions.cc:10
std::string InputFileString
Definition: Resolutions.cc:15
ROOT::VecOps::RVec< bool > bools
Definition: Resolutions.cc:5
d
Definition: ztail.py:151
vector< float > HitDXVector
Definition: Resolutions.cc:11
vector< float > TrackDXEVector
Definition: Resolutions.cc:13
ROOT::VecOps::RVec< int > ints
Definition: Resolutions.cc:4