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>;
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";
44 CutFlowReportString =
"CutFlowReport_" +
region +
"_Centimetres_ALCARECO.txt";
45 DoubleDiffString =
"hitDX-trackDX";
46 HitDXString =
"hitDX";
47 TrackDXString =
"trackDX";
48 TrackDXEString =
"trackDXE";
52 std::cout <<
"ERROR: UnitInt must be 0 or 1." << std::endl;
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";
74 CutFlowReportString =
"CutFlowReport_" +
region +
"_Centimetres_ALCARECO_UL.txt";
75 DoubleDiffString =
"hitDX-trackDX";
76 HitDXString =
"hitDX";
77 TrackDXString =
"trackDX";
78 TrackDXEString =
"trackDXE";
82 std::cout <<
"ERROR: UnitInt must be 0 or 1." << std::endl;
89 std::cout <<
"The UL input parameter must be set to 0 (for ALCARECO) or 1 (for UL ALCARECO)." << std::endl;
100 }
else if (
region ==
"TIB_L2") {
102 }
else if (
region ==
"TIB_L3") {
104 }
else if (
region ==
"TIB_L4") {
106 }
else if (
region ==
"Side_TID") {
108 }
else if (
region ==
"Wheel_TID") {
110 }
else if (
region ==
"Ring_TID") {
112 }
else if (
region ==
"TOB_L1") {
114 }
else if (
region ==
"TOB_L2") {
116 }
else if (
region ==
"TOB_L3") {
118 }
else if (
region ==
"TOB_L4") {
120 }
else if (
region ==
"TOB_L5") {
122 }
else if (
region ==
"TOB_L6") {
124 }
else if (
region ==
"Side_TEC") {
126 }
else if (
region ==
"Wheel_TEC") {
128 }
else if (
region ==
"Ring_TEC") {
130 }
else if (
region ==
"TIB_All") {
132 }
else if (
region ==
"TOB_All") {
134 }
else if (
region ==
"TID_All") {
136 }
else if (
region ==
"TEC_All") {
138 }
else if (
region ==
"Pixel_Barrel") {
140 }
else if (
region ==
"Pixel_EndcapDisk") {
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." 151 auto SubDet_Function{[&RegionInt](
const int& detID1_input,
const int& detID2_input) {
156 OutputBool = (((detID1_input >> 25) & 0x7) == 3) && ((detID1_input >> 14) & 0x7) == 1 &&
157 (((detID2_input >> 25) & 0x7) == 3) && ((detID2_input >> 14) & 0x7) == 1;
162 OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 2) &&
163 (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 2);
168 OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 3) &&
169 (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 3);
174 OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 4) &&
175 (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 4);
180 OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x3) == 1)) ||
181 ((((detID1_input >> 13) & 0x3) == 2) &&
182 (((detID2_input >> 13) & 0x3) == 2));
188 OutputBool = (((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2);
194 OutputBool = ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2));
200 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 1) &&
201 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 1);
206 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 2) &&
207 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 2);
212 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 3) &&
213 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 3);
218 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 4) &&
219 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 4);
224 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 5) &&
225 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 5);
230 OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 6) &&
231 (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 6);
236 OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) ||
237 ((((detID1_input >> 18) & 0x3) == 2) &&
238 (((detID2_input >> 18) & 0x3) == 2));
243 OutputBool = (((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4);
248 OutputBool = (((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3);
254 OutputBool = ((((detID1_input >> 25) & 0x7) == 3) && (((detID2_input >> 25) & 0x7) == 3));
260 OutputBool = ((((detID1_input >> 25) & 0x7) == 5) && (((detID2_input >> 25) & 0x7) == 5));
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));
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));
298 (((detID1_input >> 20) & 0xF) == 4) && (((detID2_input >> 20) & 0xF) == 4);
304 (((detID1_input >> 18) & 0xF) == 4) && (((detID2_input >> 18) & 0xF) == 4);
313 auto Pitch_Function{[&Unit_Int](
const float& pitch,
const float&
input) {
314 float InputOverPitch =
input / pitch;
315 return InputOverPitch;
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");
325 auto PairPathCriteriaFunction{[&RegionInt](
const float& pairPath_input) {
326 if ((RegionInt > 0 && RegionInt < 5) || (RegionInt > 7 || RegionInt < 13) || (RegionInt == 17) ||
328 return abs(pairPath_input) < 7;
330 else if (RegionInt == 21 || RegionInt == 22) {
331 return abs(pairPath_input) < 2;
334 return abs(pairPath_input) < 20;
338 auto MomentaFunction{[&RegionInt](
const float& momentum_input) {
339 if (RegionInt == 21 || RegionInt == 22) {
340 return momentum_input > 5;
343 return momentum_input > 15;
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");
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;
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);
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);
377 h_DoubleDifference->Fit(
"gaus");
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();
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);
397 h_DoubleDifference->Write();
401 h_clusterW1->Write();
402 h_clusterW2->Write();
407 auto numerator = sigma2_MeasMinusPred - sigma2_PredError;
414 std::cout <<
"The hit resolution for tracker region " <<
region <<
" is: " << HitResolution << std::endl;
418 auto allCutsReport =
d.Report();
419 std::ofstream CutFlowReport;
421 CutFlowReport.open(CutFlowReportString.c_str());
423 for (
auto&& cutInfo : allCutsReport) {
424 CutFlowReport << cutInfo.GetName() <<
'\t' << cutInfo.GetAll() <<
'\t' << cutInfo.GetPass() <<
'\t' 425 << cutInfo.GetEff() <<
" %" << std::endl;
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"};
438 for (
int i = 0;
i < LayerNames.size();
i++) {
442 std::ofstream HitResoTextFile;
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;
459 "mkdir HitResolutionValues; mkdir GaussianFits; mkdir CutFlowReports; mv CutFlowReport_* CutFlowReports/; mv " 460 "HitResolutionValues_* HitResolutionValues/; mv GaussianFits_* GaussianFits/;");
ROOT::VecOps::RVec< double > doubles
static std::string const input
std::string HitResoFileName
ROOT::VecOps::RVec< float > floats
ROOT::VecOps::RVec< UChar_t > chars
void ResolutionsCalculator(const string ®ion, const int &Unit_Int, const int &UL)
std::string GaussianFitsFileName
vector< float > TrackDXVector
vector< float > HitResolutionVector
Abs< T >::type abs(const T &t)
vector< float > DoubleDifferenceVector
std::string InputFileString
ROOT::VecOps::RVec< bool > bools
vector< float > HitDXVector
vector< float > TrackDXEVector
ROOT::VecOps::RVec< int > ints