SiStripHitResolution/macros/Resolutions.cc

0001 using namespace ROOT;
0002 using ROOT::RDF::RNode;
0003 using floats = ROOT::VecOps::RVec<float>;
0004 using ints = ROOT::VecOps::RVec<int>;
0005 using bools = ROOT::VecOps::RVec<bool>;
0006 using chars = ROOT::VecOps::RVec<UChar_t>;
0007 using doubles = ROOT::VecOps::RVec<double>;
0008
0009 vector<float> HitResolutionVector;
0010 vector<float> DoubleDifferenceVector;
0011 vector<float> HitDXVector;
0012 vector<float> TrackDXVector;
0013 vector<float> TrackDXEVector;
0014
0015 std::string InputFileString;
0016 std::string HitResoFileName;
0017 std::string GaussianFitsFileName;
0018
0019 void ResolutionsCalculator(const string& region, const int& Unit_Int, const int& UL) {
0020   std::string CutFlowReportString;
0021   std::string DoubleDiffString;
0022   std::string HitDXString;
0023   std::string TrackDXString;
0024   std::string TrackDXEString;
0025   std::string ClusterW1String = "clusterW1";
0026   std::string ClusterW2String = "clusterW2";
0027
0028   switch (UL) {
0029     case 0:
0030       switch (Unit_Int) {
0031         case 0:
0032           GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO.root";
0033           HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO.txt";
0034           CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO.txt";
0035           DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch";
0036           HitDXString = "hitDX_OverPitch";
0037           TrackDXString = "trackDX_OverPitch";
0038           TrackDXEString = "trackDXE_OverPitch";
0039           break;
0040
0041         case 1:
0042           GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO.root";
0043           HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO.txt";
0044           CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO.txt";
0045           DoubleDiffString = "hitDX-trackDX";
0046           HitDXString = "hitDX";
0047           TrackDXString = "trackDX";
0048           TrackDXEString = "trackDXE";
0049           break;
0050
0051         default:
0052           std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl;
0053           break;
0054       }
0055
0056       InputFileString = "hitresol_ALCARECO.root";
0057       break;
0058
0059     case 1:
0060       switch (Unit_Int) {
0061         case 0:
0062           GaussianFitsFileName = "GaussianFits_PitchUnits_ALCARECO_UL.root";
0063           HitResoFileName = "HitResolutionValues_PitchUnits_ALCARECO_UL.txt";
0064           CutFlowReportString = "CutFlowReport_" + region + "_PitchUnits_ALCARECO_UL.txt";
0065           DoubleDiffString = "hitDX_OverPitch-trackDX_OverPitch";
0066           HitDXString = "hitDX_OverPitch";
0067           TrackDXString = "trackDX_OverPitch";
0068           TrackDXEString = "trackDXE_OverPitch";
0069           break;
0070
0071         case 1:
0072           GaussianFitsFileName = "GaussianFits_Centimetres_ALCARECO_UL.root";
0073           HitResoFileName = "HitResolutionValues_Centimetres_ALCARECO_UL.txt";
0074           CutFlowReportString = "CutFlowReport_" + region + "_Centimetres_ALCARECO_UL.txt";
0075           DoubleDiffString = "hitDX-trackDX";
0076           HitDXString = "hitDX";
0077           TrackDXString = "trackDX";
0078           TrackDXEString = "trackDXE";
0079           break;
0080
0081         default:
0082           std::cout << "ERROR: UnitInt must be 0 or 1." << std::endl;
0083           break;
0084       }
0085
0086       InputFileString = "hitresol_ALCARECO_UL.root";
0087       break;
0088     default:
0089       std::cout << "The UL input parameter must be set to 0 (for ALCARECO) or 1 (for UL ALCARECO)." << std::endl;
0090       break;
0091   }
0092
0093   //opening the root file
0094   ROOT::RDataFrame d("anResol/reso", InputFileString);
0095
0096   int RegionInt = 0;
0097
0098   if (region == "TIB_L1") {
0099     RegionInt = 1;
0100   } else if (region == "TIB_L2") {
0101     RegionInt = 2;
0102   } else if (region == "TIB_L3") {
0103     RegionInt = 3;
0104   } else if (region == "TIB_L4") {
0105     RegionInt = 4;
0106   } else if (region == "Side_TID") {
0107     RegionInt = 5;
0108   } else if (region == "Wheel_TID") {
0109     RegionInt = 6;
0110   } else if (region == "Ring_TID") {
0111     RegionInt = 7;
0112   } else if (region == "TOB_L1") {
0113     RegionInt = 8;
0114   } else if (region == "TOB_L2") {
0115     RegionInt = 9;
0116   } else if (region == "TOB_L3") {
0117     RegionInt = 10;
0118   } else if (region == "TOB_L4") {
0119     RegionInt = 11;
0120   } else if (region == "TOB_L5") {
0121     RegionInt = 12;
0122   } else if (region == "TOB_L6") {
0123     RegionInt = 13;
0124   } else if (region == "Side_TEC") {
0125     RegionInt = 14;
0126   } else if (region == "Wheel_TEC") {
0127     RegionInt = 15;
0128   } else if (region == "Ring_TEC") {
0129     RegionInt = 16;
0130   } else if (region == "TIB_All") {
0131     RegionInt = 17;
0132   } else if (region == "TOB_All") {
0133     RegionInt = 18;
0134   } else if (region == "TID_All") {
0135     RegionInt = 19;
0136   } else if (region == "TEC_All") {
0137     RegionInt = 20;
0138   } else if (region == "Pixel_Barrel") {
0139     RegionInt = 21;
0140   } else if (region == "Pixel_EndcapDisk") {
0141     RegionInt = 22;
0142   } else {
0143     std::cout << "Error: The tracker region " << region
0144               << " was chosen. Please choose a region out of: TIB L1, TIB L2, TIB L3, TIB L4, Side TID, Wheel TID, "
0145                  "Ring TID, TOB L1, TOB L2, TOB L3, TOB L4, TOB L5, TOB L6, Side TEC, Wheel TEC or Ring TEC."
0146               << std::endl;
0147     return 0;
0148   }
0149
0150   //Lambda function to filter the detID for different layers
0151   auto SubDet_Function{[&RegionInt](const int& detID1_input, const int& detID2_input) {
0152     bool OutputBool = 0;
0153
0154     switch (RegionInt) {
0155       case 1: {
0156         OutputBool = (((detID1_input >> 25) & 0x7) == 3) && ((detID1_input >> 14) & 0x7) == 1 &&
0157                      (((detID2_input >> 25) & 0x7) == 3) && ((detID2_input >> 14) & 0x7) == 1;  //TIB L1
0158         break;
0159       }
0160
0161       case 2: {
0162         OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 2) &&
0163                      (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 2);  //TIB L2
0164         break;
0165       }
0166
0167       case 3: {
0168         OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 3) &&
0169                      (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 3);  //TIB L3
0170         break;
0171       }
0172
0173       case 4: {
0174         OutputBool = (((detID1_input >> 25) & 0x7) == 3) && (((detID1_input >> 14) & 0x7) == 4) &&
0175                      (((detID2_input >> 25) & 0x7) == 3) && (((detID2_input >> 14) & 0x7) == 4);  //TIB L4
0176         break;
0177       }
0178
0179       case 5: {
0180         OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x3) == 1)) ||
0181                      ((((detID1_input >> 13) & 0x3) == 2) &&
0182                       (((detID2_input >> 13) & 0x3) == 2));  //TID Side (1 -> TID-, 2 -> TID+)
0183
0184         break;
0185       }
0186
0187       case 6: {
0188         OutputBool = (((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2);  //TID Wheel
0189
0190         break;
0191       }
0192
0193       case 7: {
0194         OutputBool = ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2));  //TID Ring
0195
0196         break;
0197       }
0198
0199       case 8: {
0200         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 1) &&
0201                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 1);  //TOB L1
0202         break;
0203       }
0204
0205       case 9: {
0206         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 2) &&
0207                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 2);  //TOB L2
0208         break;
0209       }
0210
0211       case 10: {
0212         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 3) &&
0213                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 3);  //TOB L3
0214         break;
0215       }
0216
0217       case 11: {
0218         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 4) &&
0219                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 4);  //TOB L4
0220         break;
0221       }
0222
0223       case 12: {
0224         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 5) &&
0225                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 5);  //TOB L5
0226         break;
0227       }
0228
0229       case 13: {
0230         OutputBool = (((detID1_input >> 25) & 0x7) == 5) && (((detID1_input >> 14) & 0x7) == 6) &&
0231                      (((detID2_input >> 25) & 0x7) == 5) && (((detID2_input >> 14) & 0x7) == 6);  //TOB L6
0232         break;
0233       }
0234
0235       case 14: {
0236         OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) ||
0237                      ((((detID1_input >> 18) & 0x3) == 2) &&
0238                       (((detID2_input >> 18) & 0x3) == 2));  //Side TEC (1 -> back, 2 -> front)
0239         break;
0240       }
0241
0242       case 15: {
0243         OutputBool = (((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4);  //Wheel TEC
0244         break;
0245       }
0246
0247       case 16: {
0248         OutputBool = (((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3);  //Ring TEC
0249
0250         break;
0251       }
0252
0253       case 17: {
0254         OutputBool = ((((detID1_input >> 25) & 0x7) == 3) && (((detID2_input >> 25) & 0x7) == 3));  //All TIB
0255
0256         break;
0257       }
0258
0259       case 18: {
0260         OutputBool = ((((detID1_input >> 25) & 0x7) == 5) && (((detID2_input >> 25) & 0x7) == 5));  //All TOB
0261
0262         break;
0263       }
0264
0265       case 19: {
0266         OutputBool = ((((detID1_input >> 13) & 0x3) == 1) && (((detID2_input >> 13) & 0x7) == 1)) ||
0267                      ((((detID1_input >> 13) & 0x3) == 2) && (((detID2_input >> 13) & 0x7) == 2)) ||
0268                      ((((detID1_input >> 11) & 0x3) == 2) && (((detID2_input >> 11) & 0x3) == 2)) ||
0269                      ((((detID1_input >> 9) & 0x3) == 2) && (((detID2_input >> 9) & 0x3) == 2)) ||
0270                      ((((detID1_input >> 7) & 0x3) == 1) && (((detID2_input >> 7) & 0x3) == 1)) ||
0271                      ((((detID1_input >> 7) & 0x3) == 2) && (((detID2_input >> 7) & 0x3) == 2)) ||
0272                      ((((detID1_input >> 2) & 0x1F) == 5) && (((detID2_input >> 2) & 0x1F) == 5)) ||
0273                      ((((detID1_input >> 0) & 0x3) == 0) && (((detID2_input >> 0) & 0x3) == 0)) ||
0274                      ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) ||
0275                      ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2));  //All TID
0276
0277         break;
0278       }
0279
0280       case 20: {
0281         OutputBool = ((((detID1_input >> 18) & 0x3) == 1) && (((detID2_input >> 18) & 0x3) == 1)) ||
0282                      ((((detID1_input >> 18) & 0x3) == 2) && (((detID2_input >> 18) & 0x3) == 2)) ||
0283                      ((((detID1_input >> 14) & 0xF) == 4) && (((detID2_input >> 14) & 0xF) == 4)) ||
0284                      ((((detID1_input >> 12) & 0x3) == 1) && (((detID2_input >> 12) & 0x3) == 1)) ||
0285                      ((((detID1_input >> 12) & 0x3) == 2) && (((detID2_input >> 12) & 0x3) == 2)) ||
0286                      ((((detID1_input >> 8) & 0xF) == 4) && (((detID2_input >> 8) & 0xF) == 4)) ||
0287                      ((((detID1_input >> 5) & 0x7) == 3) && (((detID2_input >> 5) & 0x7) == 3)) ||
0288                      ((((detID1_input >> 2) & 0x7) == 3) && (((detID2_input >> 2) & 0x7) == 3)) ||
0289                      ((((detID1_input >> 0) & 0x3) == 1) && (((detID2_input >> 0) & 0x3) == 1)) ||
0290                      ((((detID1_input >> 0) & 0x3) == 2) && (((detID2_input >> 0) & 0x3) == 2)) ||
0291                      ((((detID1_input >> 0) & 0x3) == 3) && (((detID2_input >> 0) & 0x3) == 3));  //All TEC
0292
0293         break;
0294       }
0295
0296       case 21: {
0297         OutputBool =
0298             (((detID1_input >> 20) & 0xF) == 4) && (((detID2_input >> 20) & 0xF) == 4);  //pixel barrel (phase 1)
0299         break;
0300       }
0301
0302       case 22: {
0303         OutputBool =
0304             (((detID1_input >> 18) & 0xF) == 4) && (((detID2_input >> 18) & 0xF) == 4);  //pixel endcap disk (phase 1)
0305         break;
0306       }
0307     }
0308
0309     return OutputBool;
0310   }};
0311
0312   //Function for expressing the hit resolution in either micrometres or pitch units.
0313   auto Pitch_Function{[&Unit_Int](const float& pitch, const float& input) {
0314     float InputOverPitch = input / pitch;
0315     return InputOverPitch;
0316   }};
0317
0318   //Defining columns needed for the unit conversion into pitch units, and applying the filter for the subdetector
0319   auto dataframe = d.Define("hitDX_OverPitch", Pitch_Function, {"pitch1", "hitDX"})
0320                        .Define("trackDX_OverPitch", Pitch_Function, {"pitch1", "trackDX"})
0321                        .Define("trackDXE_OverPitch", Pitch_Function, {"pitch1", "trackDXE"})
0322                        .Filter(SubDet_Function, {"detID1", "detID2"}, "Subdetector filter");
0323
0324   //Implementing selection criteria that were not implemented in HitResol.cc
0325   auto PairPathCriteriaFunction{[&RegionInt](const float& pairPath_input) {
0326     if ((RegionInt > 0 && RegionInt < 5) || (RegionInt > 7 || RegionInt < 13) || (RegionInt == 17) ||
0327         (RegionInt == 18)) {
0328       return abs(pairPath_input) < 7;
0329     }  //for TIB and TOB
0330     else if (RegionInt == 21 || RegionInt == 22) {
0331       return abs(pairPath_input) < 2;
0332     }  //for pixels
0333     else {
0334       return abs(pairPath_input) < 20;
0335     }  //for everything else (max value is 15cm so this will return all events anyway)
0336   }};
0337
0338   auto MomentaFunction{[&RegionInt](const float& momentum_input) {
0339     if (RegionInt == 21 || RegionInt == 22) {
0340       return momentum_input > 5;
0341     }  //pixels
0342     else {
0343       return momentum_input > 15;
0344     }  //strips
0345   }};
0346
0347   auto dataframe_filtered =
0348       dataframe.Filter(PairPathCriteriaFunction, {"pairPath"}, "Pair path criterion filter")
0349           .Filter(MomentaFunction, {"momentum"}, "Momentum criterion filter")
0350           .Filter("trackChi2 > 0.001", "chi2 criterion filter")
0351           .Filter("numHits > 6", "numHits filter")
0352           .Filter("trackDXE < 0.0025", "trackDXE filter")
0353           .Filter("(clusterW1 == clusterW2) && clusterW1 <= 4 && clusterW2 <= 4", "cluster filter");
0354
0355   //Creating histograms for the difference between the two hit positions, the difference between the two predicted positions and for the double difference
0356   //hitDX = the difference in the hit positions for the pair
0357   //trackDX =  the difference in the track positions for the pair
0358
0359   auto HistoName_DoubleDiff = "DoubleDifference_" + region;
0360   auto HistoName_HitDX = "HitDX_" + region;
0361   auto HistoName_TrackDX = "TrackDX_" + region;
0362   auto HistoName_TrackDXE = "TrackDXE_" + region;
0363   auto HistoName_ClusterW1 = "ClusterW1_" + region;
0364   auto HistoName_ClusterW2 = "ClusterW2_" + region;
0365
0366   auto h_DoubleDifference =
0367       dataframe_filtered.Define(HistoName_DoubleDiff, DoubleDiffString)
0368           .Histo1D({HistoName_DoubleDiff.c_str(), HistoName_DoubleDiff.c_str(), 40, -0.5, 0.5}, HistoName_DoubleDiff);
0369   auto h_hitDX = dataframe_filtered.Define(HistoName_HitDX, HitDXString).Histo1D(HistoName_HitDX);
0370   auto h_trackDX = dataframe_filtered.Define(HistoName_TrackDX, TrackDXString).Histo1D(HistoName_TrackDX);
0371   auto h_trackDXE = dataframe_filtered.Define(HistoName_TrackDXE, TrackDXEString).Histo1D(HistoName_TrackDXE);
0372
0373   auto h_clusterW1 = dataframe_filtered.Define(HistoName_ClusterW1, ClusterW1String).Histo1D(HistoName_ClusterW1);
0374   auto h_clusterW2 = dataframe_filtered.Define(HistoName_ClusterW2, ClusterW2String).Histo1D(HistoName_ClusterW2);
0375
0376   //Applying gaussian fits, taking the resolutions and squaring them
0377   h_DoubleDifference->Fit("gaus");
0378
0379   auto double_diff_StdDev = h_DoubleDifference->GetStdDev();
0380   auto hitDX_StdDev = h_hitDX->GetStdDev();
0381   auto trackDX_StdDev = h_trackDX->GetStdDev();
0382   auto trackDXE_Mean = h_trackDXE->GetMean();
0383
0384   auto sigma2_MeasMinusPred = pow(double_diff_StdDev, 2);
0385   auto sigma2_Meas = pow(hitDX_StdDev, 2);
0386   auto sigma2_Pred = pow(trackDX_StdDev, 2);
0387   auto sigma2_PredError = pow(trackDXE_Mean, 2);
0388
0389   DoubleDifferenceVector.push_back(sigma2_MeasMinusPred);
0390   HitDXVector.push_back(sigma2_Meas);
0391   TrackDXVector.push_back(sigma2_Pred);
0392   TrackDXEVector.push_back(sigma2_PredError);
0393
0394   //Saving the histograms with gaussian fits applied to an output root file
0395   TFile* output = new TFile(GaussianFitsFileName.c_str(), "UPDATE");
0396
0397   h_DoubleDifference->Write();
0398   h_hitDX->Write();
0399   h_trackDX->Write();
0400   h_trackDXE->Write();
0401   h_clusterW1->Write();
0402   h_clusterW2->Write();
0403
0404   output->Close();
0405
0406   //Calculating the hit resolution;
0407   auto numerator = sigma2_MeasMinusPred - sigma2_PredError;
0408
0409   auto HitResolution = sqrt(numerator / 2);
0410   HitResolutionVector.push_back(HitResolution);
0411
0412   //Printing the resolution
0413   std::cout << '\n' << std::endl;
0414   std::cout << "The hit resolution for tracker region " << region << " is: " << HitResolution << std::endl;
0415   std::cout << '\n' << std::endl;
0416
0417   //Cut flow report
0418   auto allCutsReport = d.Report();
0419   std::ofstream CutFlowReport;
0420
0421   CutFlowReport.open(CutFlowReportString.c_str());
0422
0423   for (auto&& cutInfo : allCutsReport) {
0424     CutFlowReport << cutInfo.GetName() << '\t' << cutInfo.GetAll() << '\t' << cutInfo.GetPass() << '\t'
0425                   << cutInfo.GetEff() << " %" << std::endl;
0426   }
0427 }
0428
0429 void Resolutions() {
0430   int UnitInteger = 0;
0431   int ULInteger = 0;
0432
0433   vector<std::string> LayerNames = {"TIB_L1",   "TIB_L2",   "TIB_L3",       "TIB_L4",          "Side_TID", "Wheel_TID",
0434                                     "Ring_TID", "TOB_L1",   "TOB_L2",       "TOB_L3",          "TOB_L4",   "TOB_L5",
0435                                     "TOB_L6",   "Side_TEC", "Wheel_TEC",    "Ring_TEC",        "TIB_All",  "TOB_All",
0436                                     "TID_All",  "TEC_All",  "Pixel_Barrel", "Pixel_EndcapDisk"};
0437
0438   for (int i = 0; i < LayerNames.size(); i++) {
0439     ResolutionsCalculator(LayerNames.at(i), UnitInteger, ULInteger);
0440   }
0441
0442   std::ofstream HitResoTextFile;
0443   HitResoTextFile.open(HitResoFileName);
0444
0445   auto Width = 28;
0446
0447   HitResoTextFile << std::right << "Layer " << std::setw(Width) << " Resolution " << std::setw(Width)
0448                   << " sigma2_HitDX " << std::setw(Width) << " sigma2_trackDX " << std::setw(Width)
0449                   << " sigma2_trackDXE " << std::setw(Width) << " sigma2_DoubleDifference " << std::endl;
0450
0451   for (int i = 0; i < HitResolutionVector.size(); i++) {
0452     HitResoTextFile << std::right << LayerNames.at(i) << std::setw(Width) << HitResolutionVector.at(i)
0453                     << std::setw(Width) << HitDXVector.at(i) << std::setw(Width) << TrackDXVector.at(i)
0454                     << std::setw(Width) << TrackDXEVector.at(i) << std::setw(Width) << DoubleDifferenceVector.at(i)
0455                     << std::endl;
0456   }
0457
0458   system(
0459       "mkdir HitResolutionValues; mkdir GaussianFits; mkdir CutFlowReports; mv CutFlowReport_* CutFlowReports/; mv "
0460       "HitResolutionValues_* HitResolutionValues/; mv GaussianFits_* GaussianFits/;");
0461 }