Merge branch '21.2-truthClassifierInMSTT' into '21.2'

21.2 MC truth classifier in MuonSelectorToolsTester

See merge request atlas/athena!38341

PhysicsAnalysis/MuonID/MuonSelectorTools/CMakeLists.txt

@@ -33,7 +33,7 @@ if( XAOD_STANDALONE )
    atlas_add_executable( MuonSelectorToolsTester
       util/MuonSelectorToolsTester.cxx
       INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
-      LINK_LIBRARIES ${ROOT_LIBRARIES} xAODRootAccess xAODEventInfo xAODMuon
+      LINK_LIBRARIES ${ROOT_LIBRARIES} xAODRootAccess xAODEventInfo xAODMuon MCTruthClassifierLib
       xAODTracking xAODCore MuonSelectorToolsLib )
 endif()
 

PhysicsAnalysis/MuonID/MuonSelectorTools/util/MuonSelectorToolsTester.cxx

@@ -38,6 +38,10 @@
 #include "xAODCore/tools/IOStats.h"
 #include "xAODCore/tools/ReadStats.h"
 
+// Truth classification
+#include "MCTruthClassifier/MCTruthClassifier.h"
+#include "MCTruthClassifier/MCTruthClassifierDefs.h"
+
 /// Example of how to run the MuonSelectorTools package to obtain information from muons
 
 int main( int argc, char* argv[] ) {
@@ -143,6 +147,7 @@ int main( int argc, char* argv[] ) {
    for (int wp = 0; wp < Nwp; wp++)
      selectedMuonsNotBad[wp] = 0;
 
+
    //Obtain summary information also split by muon type
    const int Ntype = 5;
 
@@ -163,7 +168,6 @@ int main( int argc, char* argv[] ) {
        typeNames[type] = "unknown";
    }
 
-
    //Muon counters for each type
    int allMuonsType[Ntype];
    for (int type = 0; type < Ntype; type++)
@@ -206,6 +210,35 @@ int main( int argc, char* argv[] ) {
 
 
 
+   //Obtain summary information also split by muon truth type
+   const int NtruthType = 5;
+
+   std::string truthTypeNames[NtruthType] = {"Prompt", "Non-isolated", "Hadron", "Background", "Other"};
+
+   //Muon counters for each truth type
+   int allMuonsTruthType[NtruthType];
+   for (int truthType = 0; truthType < NtruthType; truthType++)
+     allMuonsTruthType[truthType] = 0;
+
+   //Muon counters for muons of each type passing each working point
+   int selectedMuonsTruthType[NtruthType][Nwp];
+   for (int truthType = 0; truthType < NtruthType; truthType++)
+     for (int wp = 0; wp < Nwp; wp++)
+       selectedMuonsTruthType[truthType][wp] = 0;
+
+   int selectedMuonsTruthTypeNotBad[NtruthType][Nwp];
+   for (int truthType = 0; truthType < NtruthType; truthType++)
+     for (int wp = 0; wp < Nwp; wp++)
+       selectedMuonsTruthTypeNotBad[truthType][wp] = 0;
+
+
+
+   //Truth classifier
+   bool isMC = false;
+   MCTruthClassifier truthClassifier("truthClassifier");
+   CHECK( truthClassifier.initialize() );
+
+
    // Loop over the events:
    for( Long64_t entry = 0; entry < entries; ++entry ) {
 
@@ -265,20 +298,39 @@ int main( int argc, char* argv[] ) {
 	      static_cast< int >( muCounter ) );
 
 
+	//Check truth origin
+	isMC = ei->eventType(xAOD::EventInfo::IS_SIMULATION);
+	std::pair<MCTruthPartClassifier::ParticleType,MCTruthPartClassifier::ParticleOrigin> truthClassification = truthClassifier.particleTruthClassifier(*mu_itr);
+
+	int truthType;
+	if (truthClassification.first == MCTruthPartClassifier::IsoMuon)
+	  truthType = 0;
+	else if (truthClassification.first == MCTruthPartClassifier::NonIsoMuon)
+	  truthType = 1;
+	else if (truthClassification.first == MCTruthPartClassifier::Hadron)
+	  truthType = 2;
+	else if (truthClassification.first == MCTruthPartClassifier::BkgMuon)
+	  truthType = 3;
+	else
+	  truthType = 4;
+
+	allMuonsTruthType[truthType]++;
+
 	if((*mu_itr)->charge() > 0)
 	  nPositive++;
 	else
 	  nNegative++;
 
-
         passesIDRequirements = selectorTools[0]->passedIDCuts(**mu_itr);
         passesPreselectionCuts = selectorTools[0]->passedMuonCuts(**mu_itr);
 	my_quality = selectorTools[0]->getQuality(**mu_itr);
 
 	//Print some general information about the muon
-	Info( APP_NAME, "Muon pT [GeV]:       %g ", std::abs((*mu_itr)->pt())/1000.);
-	Info( APP_NAME, "Muon eta, phi:       %g, %g ", (*mu_itr)->eta(),(*mu_itr)->phi());
-	Info( APP_NAME, "Muon muonType:       %d (%s)", (*mu_itr)->muonType(), typeNames[(*mu_itr)->muonType()].c_str());
+	if (isMC)
+	  Info( APP_NAME, "Muon truthType:       %d (%s)", truthClassification.first, truthTypeNames[truthType].c_str());
+	Info( APP_NAME, "Muon pT [GeV]:        %g ", std::abs((*mu_itr)->pt())/1000.);
+	Info( APP_NAME, "Muon eta, phi:        %g, %g ", (*mu_itr)->eta(),(*mu_itr)->phi());
+	Info( APP_NAME, "Muon muonType:        %d (%s)", (*mu_itr)->muonType(), typeNames[(*mu_itr)->muonType()].c_str());
 
 	Info( APP_NAME, "Muon quality (from tool, from xAOD):      %d, %d", my_quality, (*mu_itr)->quality());
 	Info( APP_NAME, "Muon passes cuts (ID hits, preselection): %d, %d", passesIDRequirements, passesPreselectionCuts);
@@ -333,6 +385,7 @@ int main( int argc, char* argv[] ) {
 	    selectedMuons[wp]++;
 	    selectedMuonsEvent[wp]++;
 	    selectedMuonsType[(*mu_itr)->muonType()][wp]++;
+	    selectedMuonsTruthType[truthType][wp]++;
 	    selectedMuonsEta[etaIndex][wp]++;
 	    selectionResults += "pass     ";
 
@@ -340,6 +393,7 @@ int main( int argc, char* argv[] ) {
 	      selectedMuonsNotBad[wp]++;
 	      selectedMuonsEventNotBad[wp]++;
 	      selectedMuonsTypeNotBad[(*mu_itr)->muonType()][wp]++;
+	      selectedMuonsTruthTypeNotBad[truthType][wp]++;
 	      selectedMuonsEtaNotBad[etaIndex][wp]++;
 	    }
 	  }
@@ -466,6 +520,43 @@ int main( int argc, char* argv[] ) {
    }
    Info(APP_NAME, "---------------------------------------------------------------------------------------");
 
+
+
+   //Make table of selected muons by truth type and working point
+   if (isMC) {
+     Info(APP_NAME, "Selected muons by truth classification and working point (numbers in parenthesis include bad muon veto):");
+     Info(APP_NAME, "---------------------------------------------------------------------------------------");
+     for (int l = 0; l < Nwp+2; l++) {
+
+       std::string line = "";
+       if (l == 0) { //line with truth classification labels
+	 line += "              ";
+	 for (int truthType = 0; truthType < NtruthType; truthType++)
+	   line += truthTypeNames[truthType] + "        ";
+       }
+       else if (l == 1) { //line for all muons inclusive
+	 line += "All muons:      ";
+	 for (int truthType = 0; truthType < NtruthType; truthType++) {
+	   std::stringstream ss;
+	   ss << std::left << std::setw(16) << std::to_string(allMuonsTruthType[truthType]);
+	   line += ss.str();
+	 }
+       }
+       else { //lines for each of the working points
+	 int wp = l - 2;
+	 line += WPnames[wp] + ":" + padding[wp] + "     ";
+	 for (int truthType = 0; truthType < NtruthType; truthType++) {
+	   std::stringstream ss;
+	   ss << std::left << std::setw(16) << (std::to_string(selectedMuonsTruthType[truthType][wp]) + " (" + std::to_string(selectedMuonsTruthTypeNotBad[truthType][wp]) + ")");
+	   line += ss.str();
+	 }
+       }
+
+       Info(APP_NAME, "%s", line.c_str());
+     }
+     Info(APP_NAME, "---------------------------------------------------------------------------------------");
+   }
+
    // Needed for Smart Slimming
    xAOD::IOStats::instance().stats().printSmartSlimmingBranchList();