.. _program_listing_file_Root_BasicEventSelection.cxx: Program Listing for File BasicEventSelection.cxx ================================================ |exhale_lsh| :ref:`Return to documentation for file ` (``Root/BasicEventSelection.cxx``) .. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS .. code-block:: cpp // EL include(s): #include #include #include "EventLoop/OutputStream.h" // EDM include(s): #include "xAODEventInfo/EventInfo.h" #include "xAODTracking/VertexContainer.h" #include "xAODCutFlow/CutBookkeeper.h" #include "xAODCutFlow/CutBookkeeperContainer.h" // package include(s): #include #include #include "PATInterfaces/CorrectionCode.h" //#include "AsgMessaging/StatusCode.h" // tools #include "GoodRunsLists/GoodRunsListSelectionTool.h" #include "PileupReweighting/PileupReweightingTool.h" #include "TrigConfxAOD/xAODConfigTool.h" //#include "PMGTools/PMGSherpa22VJetsWeightTool.h" // For reading metadata #include "AsgTools/AsgMetadataTool.h" #include "xAODMetaData/FileMetaData.h" // ROOT include(s): #include "TFile.h" #include "TTree.h" #include "TTreeFormula.h" #include "TSystem.h" #include "xAODCore/tools/IOStats.h" #include "xAODCore/tools/ReadStats.h" // this is needed to distribute the algorithm to the workers ClassImp(BasicEventSelection) BasicEventSelection :: BasicEventSelection () : Algorithm("BasicEventSelection") { } EL::StatusCode BasicEventSelection :: setupJob (EL::Job& job) { // Here you put code that sets up the job on the submission object // so that it is ready to work with your algorithm, e.g. you can // request the D3PDReader service or add output files. Any code you // put here could instead also go into the submission script. The // sole advantage of putting it here is that it gets automatically // activated/deactivated when you add/remove the algorithm from your // job, which may or may not be of value to you. // ANA_MSG_INFO( "Calling setupJob"); job.useXAOD(); // let's initialize the algorithm to use the xAODRootAccess package xAOD::Init("BasicEventSelection").ignore(); // call before opening first file EL::OutputStream outForCFlow(m_cutFlowStreamName); if(!job.outputHas(m_cutFlowStreamName) ){ job.outputAdd ( outForCFlow ); } EL::OutputStream outForMetadata(m_metaDataStreamName); if(!job.outputHas(m_metaDataStreamName) ){ job.outputAdd ( outForMetadata ); } EL::OutputStream outForDuplicates(m_duplicatesStreamName); if(!job.outputHas(m_duplicatesStreamName) ){ job.outputAdd ( outForDuplicates ); } return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: histInitialize () { // Here you do everything that needs to be done at the very // beginning on each worker node, e.g. create histograms and output // trees. This method gets called before any input files are // connected. ANA_MSG_INFO( "Calling histInitialize"); ANA_CHECK( xAH::Algorithm::algInitialize()); // write the metadata hist to this file so algos downstream can pick up the pointer TFile *fileMD = wk()->getOutputFile (m_metaDataStreamName); fileMD->cd(); // event counts from meta data if ( !m_histEventCount ) { m_histEventCount = new TH1D("MetaData_EventCount", "MetaData_EventCount", 6, 0.5, 6.5); m_histEventCount -> GetXaxis() -> SetBinLabel(1, "nEvents initial"); m_histEventCount -> GetXaxis() -> SetBinLabel(2, "nEvents selected"); m_histEventCount -> GetXaxis() -> SetBinLabel(3, "sumOfWeights initial"); m_histEventCount -> GetXaxis() -> SetBinLabel(4, "sumOfWeights selected"); m_histEventCount -> GetXaxis() -> SetBinLabel(5, "sumOfWeightsSquared initial"); m_histEventCount -> GetXaxis() -> SetBinLabel(6, "sumOfWeightsSquared selected"); } TFile *fileSumW = wk()->getOutputFile (m_metaDataStreamName); fileSumW->cd(); // event counts from meta data if ( !m_histSumW ) { m_histSumW = new TH1D("MetaData_SumW", "MetaData_SumW", 1, -0.5, 0.5); m_histSumW->SetCanExtend(TH1::kAllAxes); } ANA_MSG_INFO( "Creating histograms"); // write the cutflows to this file so algos downstream can pick up the pointer // TFile *fileCF = wk()->getOutputFile (m_cutFlowStreamName); fileCF->cd(); // Note: the following code is needed for anyone developing/running in ROOT 6.04.10+ // Bin extension is not done anymore via TH1::SetBit(TH1::kCanRebin), but with TH1::SetCanExtend(TH1::kAllAxes) //initialise event cutflow, which will be picked ALSO by the algos downstream where an event selection is applied (or at least can be applied) // // use 1,1,2 so Fill(bin) and GetBinContent(bin) refer to the same bin // m_cutflowHist = new TH1D("cutflow", "cutflow", 1, 1, 2); m_cutflowHist->SetCanExtend(TH1::kAllAxes); // use 1,1,2 so Fill(bin) and GetBinContent(bin) refer to the same bin // m_cutflowHistW = new TH1D("cutflow_weighted", "cutflow_weighted", 1, 1, 2); m_cutflowHistW->SetCanExtend(TH1::kAllAxes); // initialise object cutflows, which will be picked by the object selector algos downstream and filled. // m_el_cutflowHist_1 = new TH1D("cutflow_electrons_1", "cutflow_electrons_1", 1, 1, 2); m_el_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_el_cutflowHist_2 = new TH1D("cutflow_electrons_2", "cutflow_electrons_2", 1, 1, 2); m_el_cutflowHist_2->SetCanExtend(TH1::kAllAxes); m_mu_cutflowHist_1 = new TH1D("cutflow_muons_1", "cutflow_muons_1", 1, 1, 2); m_mu_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_mu_cutflowHist_2 = new TH1D("cutflow_muons_2", "cutflow_muons_2", 1, 1, 2); m_mu_cutflowHist_2->SetCanExtend(TH1::kAllAxes); m_ph_cutflowHist_1 = new TH1D("cutflow_photons_1", "cutflow_photons_1", 1, 1, 2); m_ph_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_tau_cutflowHist_1 = new TH1D("cutflow_taus_1", "cutflow_taus_1", 1, 1, 2); m_tau_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_tau_cutflowHist_2 = new TH1D("cutflow_taus_2", "cutflow_taus_2", 1, 1, 2); m_tau_cutflowHist_2->SetCanExtend(TH1::kAllAxes); m_jet_cutflowHist_1 = new TH1D("cutflow_jets_1", "cutflow_jets_1", 1, 1, 2); m_jet_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_trk_cutflowHist_1 = new TH1D("cutflow_trks_1", "cutflow_trks_1", 1, 1, 2); m_trk_cutflowHist_1->SetCanExtend(TH1::kAllAxes); m_truth_cutflowHist_1 = new TH1D("cutflow_truths_1", "cutflow_truths_1", 1, 1, 2); m_truth_cutflowHist_1->SetCanExtend(TH1::kAllAxes); // start labelling the bins for the event cutflow // m_cutflow_all = m_cutflowHist->GetXaxis()->FindBin("all"); m_cutflowHistW->GetXaxis()->FindBin("all"); m_cutflow_init = m_cutflowHist->GetXaxis()->FindBin("init"); m_cutflowHistW->GetXaxis()->FindBin("init"); // extra cutflows // create histograms when m_doRunByRunCutflows is set to true // later check if !isMC() for filling so that these only get // filled when running on data if (m_doRunByRunCutflows) { m_runByrun_beforeCuts = new TH1D("runByrun_cutflow_beforeCuts", "Run-by-Run cutflow before cuts", 1, 1, 2); m_runByrun_beforeCuts->SetCanExtend(TH1::kAllAxes); m_runByrun_afterCuts = new TH1D("runByrun_cutflow_afterCuts", "Run-by-Run cutflow after cuts", 1, 1, 2); m_runByrun_afterCuts->SetCanExtend(TH1::kAllAxes); } ANA_MSG_INFO( "Finished creating histograms"); return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: fileExecute () { // Here you do everything that needs to be done exactly once for every // single file, e.g. collect a list of all lumi-blocks processed ANA_MSG_INFO( "Calling fileExecute"); // get TEvent and TStore - must be done here b/c we need to retrieve CutBookkeepers container from TEvent! // m_event = wk()->xaodEvent(); m_store = wk()->xaodStore(); // get the MetaData tree once a new file is opened, with // TTree* MetaData = dynamic_cast( wk()->inputFile()->Get("MetaData") ); if ( !MetaData ) { ANA_MSG_ERROR( "MetaData tree not found! Exiting."); return EL::StatusCode::FAILURE; } MetaData->LoadTree(0); //--------------------------- // Meta data - CutBookkepers //--------------------------- // // Metadata for intial N (weighted) events are used to correctly normalise MC // if running on a MC DAOD which had some skimming applied at the derivation stage //check if file is from a DxAOD bool m_isDerivation = !MetaData->GetBranch("StreamAOD"); if ( m_useMetaData ) { // Check for potential file corruption // // If there are some Incomplete CBK, throw a WARNING, // unless ALL of them have inputStream == "unknownStream" // const xAOD::CutBookkeeperContainer* incompleteCBC(nullptr); if ( !m_event->retrieveMetaInput(incompleteCBC, "IncompleteCutBookkeepers").isSuccess() ) { ANA_MSG_ERROR("Failed to retrieve IncompleteCutBookkeepers from MetaData! Exiting."); return EL::StatusCode::FAILURE; } bool allFromUnknownStream(true); if ( incompleteCBC->size() != 0 ) { std::string stream(""); for ( auto cbk : *incompleteCBC ) { ANA_MSG_INFO("Incomplete cbk name: " << cbk->name() << " - stream: " << cbk->inputStream()); if ( cbk->inputStream() != "unknownStream" ) { allFromUnknownStream = false; stream = cbk->inputStream(); break; } } if ( !allFromUnknownStream ) { ANA_MSG_WARNING("Found incomplete CBK from stream: " << stream << ". This is not necessarily a sign of file corruption (incomplete CBK appear when 'maxevents' is set in the AOD jo, for instance), but you may still want to check input file for potential corruption..." ); } } // Now, let's find the actual information // const xAOD::CutBookkeeperContainer* completeCBC(nullptr); if ( !m_event->retrieveMetaInput(completeCBC, "CutBookkeepers").isSuccess() ) { ANA_MSG_ERROR("Failed to retrieve CutBookkeepers from MetaData! Exiting."); return EL::StatusCode::FAILURE; } // Find the smallest cycle number, the original first processing step/cycle int minCycle(10000); for ( auto cbk : *completeCBC ) { if ( !( cbk->name().empty() ) && ( minCycle > cbk->cycle() ) ){ minCycle = cbk->cycle(); } } // Now, let's actually find the right one that contains all the needed info... const xAOD::CutBookkeeper* allEventsCBK(nullptr); const xAOD::CutBookkeeper* DxAODEventsCBK(nullptr); if ( m_isDerivation ) { if(m_derivationName != ""){ ANA_MSG_INFO("Override auto config to look at DAOD made by Derivation Algorithm: " << m_derivationName); }else{ ANA_MSG_INFO("Will autoconfig to look at DAOD made by Derivation Algorithm."); } } int maxCycle(-1); for ( const xAOD::CutBookkeeper* cbk: *completeCBC ) { // Find initial book keeper ANA_MSG_INFO("Complete cbk name: " << cbk->name() << " - stream: " << cbk->inputStream() ); if( cbk->cycle() > maxCycle && cbk->name() == "AllExecutedEvents" && cbk->inputStream() == "StreamAOD" ) { allEventsCBK = cbk; maxCycle = cbk->cycle(); } // Find derivation book keeper if ( m_isDerivation ) { if(m_derivationName != "") { if ( cbk->name() == m_derivationName ) DxAODEventsCBK = cbk; } else if( cbk->name().find("Kernel") != std::string::npos ) { ANA_MSG_INFO("Auto config found DAOD made by Derivation Algorithm: " << cbk->name()); DxAODEventsCBK = cbk; } } // is derivation // Find and record AOD-level sumW information for all alternate weights // The nominal AllExecutedEvents will be filled later, due to posibility of multiple CBK entries if(cbk->name().substr(0,37) == "AllExecutedEvents_NonNominalMCWeight_" && cbk->name().length()!=17 && cbk->inputStream() == "StreamAOD") { // Extract weight index from name int32_t idx=std::stoi(cbk->name().substr(37)); // Fill histogram, making sure that there is space // Note will fill bin index = idx+1 while(idx>=m_histSumW->GetNbinsX()) m_histSumW->LabelsInflate("X"); m_histSumW->Fill(idx, cbk->sumOfEventWeights()); } } if(allEventsCBK == nullptr) { ANA_MSG_WARNING("No allEventsCBK found (this is expected for DataScouting, otherwise not). Event numbers set to 0."); m_MD_initialNevents = 0; m_MD_initialSumW = 0; m_MD_initialSumWSquared = 0; } else { m_MD_initialNevents = allEventsCBK->nAcceptedEvents(); m_MD_initialSumW = allEventsCBK->sumOfEventWeights(); m_MD_initialSumWSquared = allEventsCBK->sumOfEventWeightsSquared(); } if ( m_isDerivation && !DxAODEventsCBK ) { ANA_MSG_ERROR( "No CutBookkeeper corresponds to the selected Derivation Framework algorithm name. Check it with your DF experts! Aborting."); return EL::StatusCode::FAILURE; } m_MD_finalNevents = ( m_isDerivation ) ? DxAODEventsCBK->nAcceptedEvents() : m_MD_initialNevents; m_MD_finalSumW = ( m_isDerivation ) ? DxAODEventsCBK->sumOfEventWeights() : m_MD_initialSumW; m_MD_finalSumWSquared = ( m_isDerivation ) ? DxAODEventsCBK->sumOfEventWeightsSquared() : m_MD_initialSumWSquared; // Write metadata event bookkeepers to histogram // ANA_MSG_INFO( "Meta data from this file:"); ANA_MSG_INFO( "Initial events = " << static_cast(m_MD_initialNevents) ); ANA_MSG_INFO( "Selected events = " << static_cast(m_MD_finalNevents) ); ANA_MSG_INFO( "Initial sum of weights = " << m_MD_initialSumW); ANA_MSG_INFO( "Selected sum of weights = " << m_MD_finalSumW); ANA_MSG_INFO( "Initial sum of weights squared = " << m_MD_initialSumWSquared); ANA_MSG_INFO( "Selected sum of weights squared = " << m_MD_finalSumWSquared); m_cutflowHist ->Fill(m_cutflow_all, m_MD_initialNevents); m_cutflowHistW->Fill(m_cutflow_all, m_MD_initialSumW); m_histSumW->Fill(0., m_MD_initialSumW); m_histEventCount -> Fill(1, m_MD_initialNevents); m_histEventCount -> Fill(2, m_MD_finalNevents); m_histEventCount -> Fill(3, m_MD_initialSumW); m_histEventCount -> Fill(4, m_MD_finalSumW); m_histEventCount -> Fill(5, m_MD_initialSumWSquared); m_histEventCount -> Fill(6, m_MD_finalSumWSquared); } return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: changeInput (bool /*firstFile*/) { // Here you do everything you need to do when we change input files, // e.g. resetting branch addresses on trees. If you are using // D3PDReader or a similar service this method is not needed. return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: initialize () { // Here you do everything that you need to do after the first input // file has been connected and before the first event is processed, // e.g. create additional histograms based on which variables are // available in the input files. You can also create all of your // histograms and trees in here, but be aware that this method // doesn't get called if no events are processed. So any objects // you create here won't be available in the output if you have no // input events. ANA_MSG_INFO( "Initializing BasicEventSelection... "); // if truth level make sure parameters are set properly if( m_truthLevelOnly ) { ANA_MSG_INFO( "Truth only! Turn off trigger stuff"); m_triggerSelection = ""; m_extraTriggerSelection = ""; m_applyTriggerCut = m_storeTrigDecisions = m_storePassL1 = m_storePassHLT = m_storeTrigKeys = false; ANA_MSG_INFO( "Truth only! Turn off GRL"); m_applyGRLCut = false; ANA_MSG_INFO( "Truth only! Turn off Pile-up Reweight"); m_doPUreweighting = false; m_doPUreweightingSys = false; } // get TEvent and TStore // m_event = wk()->xaodEvent(); m_store = wk()->xaodStore(); const xAOD::EventInfo* eventInfo(nullptr); ANA_CHECK( HelperFunctions::retrieve(eventInfo, m_eventInfoContainerName, m_event, m_store, msg()) ); ANA_MSG_DEBUG( "Is MC? " << isMC() ); //Protection in case GRL does not apply to this run // if ( m_applyGRLCut ) { std::string runNumString = std::to_string(eventInfo->runNumber()); if ( m_GRLExcludeList.find( runNumString ) != std::string::npos ) { ANA_MSG_INFO( "RunNumber is in GRLExclusion list, setting applyGRL to false"); m_applyGRLCut = false; } } m_cleanPowheg = false; if ( eventInfo->runNumber() == 426005 ) { // Powheg+Pythia J5 m_cleanPowheg = true; ANA_MSG_INFO( "This is J5 Powheg - cleaning that nasty huge weight event"); } // https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/CentralMC15ProductionList#Sherpa_v2_2_0_V_jets_NJet_reweig //m_reweightSherpa22 = false; //if( isMC() && // ( (eventInfo->mcChannelNumber() >= 363331 && eventInfo->mcChannelNumber() <= 363483 ) || // (eventInfo->mcChannelNumber() >= 363102 && eventInfo->mcChannelNumber() <= 363122 ) || // (eventInfo->mcChannelNumber() >= 363361 && eventInfo->mcChannelNumber() <= 363435 ) ) ){ // ANA_MSG_INFO( "This is Sherpa 2.2 dataset and should be reweighted. An extra weight will be saved to EventInfo called \"weight_Sherpa22\"."); // m_reweightSherpa22 = true; // //Choose Jet Truth container, WZ has more information and is favored by the tool // std::string pmg_TruthJetContainer = ""; // if( m_event->contains("AntiKt4TruthWZJets") ){ // pmg_TruthJetContainer = "AntiKt4TruthWZJets"; // } else if( m_event->contains("AntiKt4TruthJets") ){ // pmg_TruthJetContainer = "AntiKt4TruthJets"; // } else { // ANA_MSG_WARNING( "No Truth Jet Container found for Sherpa 22 reweighting, weight_Sherpa22 will not be set."); // m_reweightSherpa22 = false; // } // //Initialize Tool // if( m_reweightSherpa22 ){ // ANA_CHECK( m_reweightSherpa22_tool_handle.setProperty( "TruthJetContainer", pmg_TruthJetContainer )); // ANA_CHECK( m_reweightSherpa22_tool_handle.setProperty( "OutputLevel", msg().level() )); // ANA_CHECK( m_reweightSherpa22_tool_handle.retrieve()); // ANA_MSG_DEBUG("Retrieved tool: " << m_reweightSherpa22_tool_handle); // } //}//if isMC and a Sherpa 2.2 sample ANA_MSG_INFO( "Setting up histograms based on isMC()"); if ( ( !isMC() && m_checkDuplicatesData ) || ( isMC() && m_checkDuplicatesMC ) ) { m_cutflow_duplicates = m_cutflowHist->GetXaxis()->FindBin("Duplicates"); m_cutflowHistW->GetXaxis()->FindBin("Duplicates"); } if ( !isMC() ) { if ( m_applyGRLCut ) { m_cutflow_grl = m_cutflowHist->GetXaxis()->FindBin("GRL"); m_cutflowHistW->GetXaxis()->FindBin("GRL"); } m_cutflow_lar = m_cutflowHist->GetXaxis()->FindBin("LAr"); m_cutflowHistW->GetXaxis()->FindBin("LAr"); m_cutflow_tile = m_cutflowHist->GetXaxis()->FindBin("tile"); m_cutflowHistW->GetXaxis()->FindBin("tile"); m_cutflow_SCT = m_cutflowHist->GetXaxis()->FindBin("SCT"); m_cutflowHistW->GetXaxis()->FindBin("SCT"); m_cutflow_core = m_cutflowHist->GetXaxis()->FindBin("core"); m_cutflowHistW->GetXaxis()->FindBin("core"); } if ( m_applyJetCleaningEventFlag ) { m_cutflow_jetcleaning = m_cutflowHist->GetXaxis()->FindBin("JetCleaning"); m_cutflowHistW->GetXaxis()->FindBin("JetCleaning"); } if ( !isMC() ) { if ( m_applyIsBadBatmanFlag ) { m_cutflow_isbadbatman = m_cutflowHist->GetXaxis()->FindBin("IsBadBatman"); m_cutflowHistW->GetXaxis()->FindBin("IsBadBatman"); } } m_cutflow_npv = m_cutflowHist->GetXaxis()->FindBin("NPV"); m_cutflowHistW->GetXaxis()->FindBin("NPV"); if ( !m_triggerSelection.empty() && m_applyTriggerCut ) { m_cutflow_trigger = m_cutflowHist->GetXaxis()->FindBin("Trigger"); m_cutflowHistW->GetXaxis()->FindBin("Trigger"); } ANA_MSG_INFO( "Histograms set up!"); // ------------------------------------------------------------------------------------------------- // Create TTree for bookeeeping duplicated events // TFile *fileDPL = wk()->getOutputFile (m_duplicatesStreamName); fileDPL->cd(); m_duplicatesTree = new TTree("duplicates","Info on duplicated events"); m_duplicatesTree->Branch("runNumber", &m_duplRunNumber, "runNumber/I"); m_duplicatesTree->Branch("eventNumber", &m_duplEventNumber, "eventNumber/L"); // ------------------------------------------------------------------------------------------------- ANA_MSG_INFO( "Setting Up Tools"); // 1. // initialize the GoodRunsListSelectionTool // if(m_applyGRLCut){ std::vector vecStringGRL; std::string grl; std::istringstream ss(m_GRLxml); while ( std::getline(ss, grl, ',') ) { std::string file = PathResolverFindCalibFile(grl); ANA_MSG_DEBUG("Found GRL: " << file); vecStringGRL.push_back(file); } ANA_CHECK( m_grl_handle.setProperty("GoodRunsListVec", vecStringGRL)); ANA_CHECK( m_grl_handle.setProperty("PassThrough", false)); ANA_CHECK( m_grl_handle.setProperty("OutputLevel", msg().level())); ANA_CHECK( m_grl_handle.retrieve()); ANA_MSG_DEBUG("Retrieved tool: " << m_grl_handle); } // 2. // initialize the Trig::TrigDecisionTool // if( !m_triggerSelection.empty() || !m_extraTriggerSelection.empty() || m_applyTriggerCut || m_storeTrigDecisions || m_storePassL1 || m_storePassHLT || m_storeTrigKeys ) { ANA_CHECK( m_trigConfTool_handle.setProperty("OutputLevel", msg().level())); ANA_CHECK( m_trigConfTool_handle.retrieve()); ANA_MSG_DEBUG("Retrieved tool: " << m_trigConfTool_handle); ANA_CHECK( m_trigDecTool_handle.setProperty( "ConfigTool", m_trigConfTool_handle )); ANA_CHECK( m_trigDecTool_handle.setProperty( "TrigDecisionKey", "xTrigDecision" )); ANA_CHECK( m_trigDecTool_handle.setProperty( "OutputLevel", msg().level() )); if ( m_useRun3navigation ) { ANA_CHECK( m_trigDecTool_handle.setProperty( "NavigationFormat", "TrigComposite") ); ANA_CHECK( m_trigDecTool_handle.setProperty( "HLTSummary", m_HLTSummary) ); } else { ANA_CHECK( m_trigDecTool_handle.setProperty( "NavigationFormat", "TriggerElement") ); } ANA_CHECK( m_trigDecTool_handle.retrieve()); ANA_MSG_DEBUG("Retrieved tool: " << m_trigDecTool_handle); // parse extra triggers list and split by comma std::string token; std::istringstream ss(m_extraTriggerSelection); while (std::getline(ss, token, ',')) { m_extraTriggerSelectionList.push_back(token); } }//end trigger configuration // 3. // initialize the CP::PileupReweightingTool // if ( m_doPUreweighting ) { std::vector PRWFiles; std::vector lumiCalcFiles; std::string tmp_lumiCalcFileNames = m_lumiCalcFileNames; std::string tmp_PRWFileNames = m_PRWFileNames; // Parse all comma seperated files // while ( tmp_PRWFileNames.size() > 0) { size_t pos = tmp_PRWFileNames.find_first_of(','); if ( pos == std::string::npos ) { pos = tmp_PRWFileNames.size(); PRWFiles.push_back(tmp_PRWFileNames.substr(0, pos)); tmp_PRWFileNames.erase(0, pos); } else { PRWFiles.push_back(tmp_PRWFileNames.substr(0, pos)); tmp_PRWFileNames.erase(0, pos+1); } } while ( tmp_lumiCalcFileNames.size() > 0) { size_t pos = tmp_lumiCalcFileNames.find_first_of(','); if ( pos == std::string::npos ) { pos = tmp_lumiCalcFileNames.size(); lumiCalcFiles.push_back(tmp_lumiCalcFileNames.substr(0, pos)); tmp_lumiCalcFileNames.erase(0, pos); } else { lumiCalcFiles.push_back(tmp_lumiCalcFileNames.substr(0, pos)); tmp_lumiCalcFileNames.erase(0, pos+1); } } // Find trigger specific lumicalc files if ( !isMC() ) { for ( const std::string &lumiCalcFile : lumiCalcFiles) { size_t pos = lumiCalcFile.find_first_of(':'); if ( pos != std::string::npos ) { m_triggerUnprescaleList.push_back(lumiCalcFile.substr(pos + 1)); } } if ( !m_triggerUnprescaleList.empty() ) { ANA_MSG_INFO("*** Trigger chains used for data unprescaling:"); for ( const std::string &trigger : m_triggerUnprescaleList ) { ANA_MSG_INFO( "\t " << trigger ); } } } if(m_autoconfigPRW) { ANA_CHECK( autoconfigurePileupRWTool() ); } else { if ( isMC() ) { ANA_MSG_INFO( "Adding Pileup files for CP::PileupReweightingTool:"); for( unsigned int i=0; i < PRWFiles.size(); ++i){ ANA_MSG_INFO( "\t- " << PRWFiles.at(i).c_str()); } ANA_CHECK( m_pileup_tool_handle.setProperty("ConfigFiles", PRWFiles)); } ANA_MSG_INFO( "Adding LumiCalc files for CP::PileupReweightingTool:"); for( unsigned int i=0; i < lumiCalcFiles.size(); ++i){ ANA_MSG_INFO( "\t- " << lumiCalcFiles.at(i).c_str()); } ANA_CHECK( m_pileup_tool_handle.setProperty("LumiCalcFiles", lumiCalcFiles)); } ANA_CHECK( m_pileup_tool_handle.setProperty("UsePeriodConfig", m_periodConfig) ); ANA_CHECK( m_pileup_tool_handle.setProperty("OutputLevel", msg().level() )); if ( !m_triggerUnprescaleList.empty() ) { // We need to make an instance of ITrigDecisionTool: asg::AnaToolHandle iTrigDecTool_handle {"Trig::TrigDecisionTool/TrigDecisionTool"}; if ( !iTrigDecTool_handle.isUserConfigured() ) { ANA_MSG_FATAL("A configured " << iTrigDecTool_handle.typeAndName() << " must have been previously created!"); return EL::StatusCode::FAILURE; } ANA_CHECK( iTrigDecTool_handle.retrieve() ); ANA_CHECK( m_pileup_tool_handle.setProperty("TrigDecisionTool", iTrigDecTool_handle )); } ANA_CHECK( m_pileup_tool_handle.retrieve()); ANA_MSG_DEBUG("Retrieved tool: " << m_pileup_tool_handle); } // As a check, let's see the number of events in our file (long long int) // ANA_MSG_INFO( "Number of events in file = " << m_event->getEntries()); // Initialize counter for number of entries m_eventCounter = 0; ANA_MSG_INFO( "BasicEventSelection succesfully initialized!"); return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: execute () { // Here you do everything that needs to be done on every single // events, e.g. read input variables, apply cuts, and fill // histograms and trees. This is where most of your actual analysis // code will go. ANA_MSG_DEBUG( "Basic Event Selection"); // Print every 1000 entries, so we know where we are: // if ( (m_eventCounter % 1000) == 0 ) { ANA_MSG_INFO( "Entry number = " << m_eventCounter); ANA_MSG_VERBOSE( "Store Content:"); if(msgLvl(MSG::VERBOSE)) m_store->print(); ANA_MSG_VERBOSE( "End Content"); } //----------------------------------------- // Print triggers used for first entry only // and fill the trigger expression for // unprescaling data //----------------------------------------- std::string TriggerExpression = ""; if ( !m_triggerSelection.empty() ) { if (m_eventCounter == 0) { if (m_eventCounter == 0) ANA_MSG_INFO( "*** Triggers used (in OR) are:\n"); auto printingTriggerChainGroup = m_trigDecTool_handle->getChainGroup(m_triggerSelection); std::vector triggersUsed = printingTriggerChainGroup->getListOfTriggers(); for ( unsigned int iTrigger = 0; iTrigger < triggersUsed.size(); ++iTrigger ) { if (m_eventCounter == 0) ANA_MSG_INFO("\t" << triggersUsed.at(iTrigger).c_str()); TriggerExpression.append(triggersUsed.at(iTrigger).c_str()); if ( iTrigger != triggersUsed.size() - 1) TriggerExpression.append(" || "); } } } if ( m_eventCounter == 0 && !m_extraTriggerSelection.empty() ) { ANA_MSG_INFO( "*** Extra Trigger Info Saved are :\n"); for ( const std::string &trigName : m_extraTriggerSelectionList ) { ANA_MSG_INFO("\t" << trigName.c_str()); ANA_MSG_INFO("\tEvaluates to:"); auto printingTriggerChainGroup = m_trigDecTool_handle->getChainGroup(trigName); std::vector triggersUsed = printingTriggerChainGroup->getListOfTriggers(); for ( unsigned int iTrigger = 0; iTrigger < triggersUsed.size(); ++iTrigger ) { ANA_MSG_INFO("\t- "<< triggersUsed.at(iTrigger).c_str()); } } } ++m_eventCounter; //------------------ // Grab event //------------------ const xAOD::EventInfo* eventInfo(nullptr); ANA_CHECK( HelperFunctions::retrieve(eventInfo, m_eventInfoContainerName, m_event, m_store, msg()) ); //------------------------------------------------------------------------------------------ // Declare an 'eventInfo' decorator with the MC event weight //------------------------------------------------------------------------------------------ static SG::AuxElement::Decorator< float > mcEvtWeightDecor("mcEventWeight"); static SG::AuxElement::Accessor< float > mcEvtWeightAcc("mcEventWeight"); float mcEvtWeight(1.0); // Check if need to create xAH event weight // if ( !mcEvtWeightDecor.isAvailable(*eventInfo) ) { if ( isMC() ) { const std::vector< float > weights = eventInfo->mcEventWeights(); // The weight (and systs) of all the MC events used in the simulation if ( weights.size() > 0 ) mcEvtWeight = weights[0]; //for ( auto& it : weights ) { ANA_MSG_INFO( "event weight: %2f.", it ); } // kill the powheg event with a huge weight if( m_cleanPowheg ) { if( eventInfo->eventNumber() == 1652845 ) { ANA_MSG_INFO("Dropping huge weight event. Weight should be 352220000"); ANA_MSG_INFO("WEIGHT : " << mcEvtWeight); wk()->skipEvent(); return EL::StatusCode::SUCCESS; // go to next event } } } // Decorate event with the *total* MC event weight // mcEvtWeightDecor(*eventInfo) = mcEvtWeight; } else { mcEvtWeight = mcEvtWeightAcc(*eventInfo); } //------------------------------------------------------------------------------------------ // Fill cutflows for run-by-run checks on data before cuts //------------------------------------------------------------------------------------------ uint32_t runNumberForCutflow = (uint32_t) eventInfo->runNumber(); if (m_doRunByRunCutflows && !isMC()) { m_runByrun_beforeCuts->Fill(TString(std::to_string(runNumberForCutflow)), 1.); } //------------------------------------------------------------------------------------------ // Declare an 'eventInfo' decorator with the Sherpa 2.2 reweight to multijet truth // https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/CentralMC15ProductionList#Sherpa_v2_2_0_V_jets_NJet_reweig //------------------------------------------------------------------------------------------ //if ( m_reweightSherpa22 ){ // static SG::AuxElement::Decorator< float > weight_Sherpa22Decor("weight_Sherpa22"); // // Check if weight needs to be added // if ( !weight_Sherpa22Decor.isAvailable(*eventInfo) ) { // float weight_Sherpa22 = -999.; // weight_Sherpa22 = m_reweightSherpa22_tool_handle->getWeight(); // weight_Sherpa22Decor( *eventInfo ) = weight_Sherpa22; // ANA_MSG_DEBUG("Setting Sherpa 2.2 reweight to " << weight_Sherpa22); // } // If not already decorated //} // if m_reweightSherpa22 //------------------------------------------------------------------------------------------ // Fill initial bin of cutflow //------------------------------------------------------------------------------------------ if( !m_useMetaData ) { m_cutflowHist ->Fill( m_cutflow_all, 1 ); m_cutflowHistW->Fill( m_cutflow_all, mcEvtWeight); m_histEventCount -> Fill(1, 1); m_histEventCount -> Fill(2, 1); m_histEventCount -> Fill(3, mcEvtWeight); m_histEventCount -> Fill(4, mcEvtWeight); m_histEventCount -> Fill(5, mcEvtWeight*mcEvtWeight); m_histEventCount -> Fill(6, mcEvtWeight*mcEvtWeight); } m_cutflowHist ->Fill( m_cutflow_init, 1 ); m_cutflowHistW->Fill( m_cutflow_init, mcEvtWeight); //-------------------------------------------------------------------------------------------------------- // Check current event is not a duplicate // This is done by checking against the std::set of filled for all previous events //-------------------------------------------------------------------------------------------------------- if ( ( !isMC() && m_checkDuplicatesData ) || ( isMC() && m_checkDuplicatesMC ) ) { std::pair thispair = std::make_pair(eventInfo->runNumber(),eventInfo->eventNumber()); if ( m_RunNr_VS_EvtNr.find(thispair) != m_RunNr_VS_EvtNr.end() ) { ANA_MSG_WARNING("Found duplicated event! runNumber = " << static_cast(eventInfo->runNumber()) << ", eventNumber = " << static_cast(eventInfo->eventNumber()) << ". Skipping this event"); // Bookkeep info in duplicates TTree // m_duplRunNumber = eventInfo->runNumber(); m_duplEventNumber = eventInfo->eventNumber(); m_duplicatesTree->Fill(); wk()->skipEvent(); return EL::StatusCode::SUCCESS; // go to next event } m_RunNr_VS_EvtNr.insert(thispair); m_cutflowHist ->Fill( m_cutflow_duplicates, 1 ); m_cutflowHistW->Fill( m_cutflow_duplicates, mcEvtWeight); } //------------------------------------------------------------------------------------------ // Update Pile-Up Reweighting //------------------------------------------------------------------------------------------ if ( m_doPUreweighting ) { m_pileup_tool_handle->applySystematicVariation(CP::SystematicSet()).ignore(); ANA_CHECK(m_pileup_tool_handle->apply( *eventInfo )); // NB: this call automatically decorates eventInfo with: // 1.) the PU weight ("PileupWeight") // 2.) the corrected mu ("corrected_averageInteractionsPerCrossing") // 3.) the random run number ("RandomRunNumber") // 4.) the random lumiblock number ("RandomLumiBlockNumber") // static SG::AuxElement::Decorator< float > correctedAvgMu("corrected_averageInteractionsPerCrossing"); static SG::AuxElement::Decorator< float > correctedAndScaledAvgMu("correctedScaled_averageInteractionsPerCrossing"); static SG::AuxElement::Decorator< float > correctedMu("corrected_actualInteractionsPerCrossing"); static SG::AuxElement::Decorator< float > correctedAndScaledMu("correctedScaled_actualInteractionsPerCrossing"); correctedAndScaledAvgMu( *eventInfo ) = m_pileup_tool_handle->getCorrectedAverageInteractionsPerCrossing( *eventInfo, true ); correctedMu( *eventInfo ) = m_pileup_tool_handle->getCorrectedActualInteractionsPerCrossing( *eventInfo ); correctedAndScaledMu( *eventInfo ) = m_pileup_tool_handle->getCorrectedActualInteractionsPerCrossing( *eventInfo, true ); if ( isMC() && m_doPUreweightingSys ) { CP::SystematicSet tmpSet;tmpSet.insert(CP::SystematicVariation("PRW_DATASF",1)); m_pileup_tool_handle->applySystematicVariation( tmpSet ).ignore(); eventInfo->auxdecor< float >( "PileupWeight_UP" )= m_pileup_tool_handle->getCombinedWeight( *eventInfo ); tmpSet.clear();tmpSet.insert(CP::SystematicVariation("PRW_DATASF",-1)); m_pileup_tool_handle->applySystematicVariation( tmpSet ).ignore(); eventInfo->auxdecor< float >( "PileupWeight_DOWN")= m_pileup_tool_handle->getCombinedWeight( *eventInfo ); } } if ( m_actualMuMin > 0 ) { // apply minimum pile-up cut if ( eventInfo->actualInteractionsPerCrossing() < m_actualMuMin ) { // veto event wk()->skipEvent(); return EL::StatusCode::SUCCESS; } } if ( m_actualMuMax > 0 ) { // apply maximum pile-up cut if ( eventInfo->actualInteractionsPerCrossing() > m_actualMuMax ) { // veto event wk()->skipEvent(); return EL::StatusCode::SUCCESS; } } //------------------------------------------------------ // If data, check if event passes GRL and event cleaning //------------------------------------------------------ if ( !isMC() ) { // Get the streams that the event was put in if (m_checkStreams){ const std::vector< xAOD::EventInfo::StreamTag > streams = eventInfo->streamTags(); for ( auto& stream : streams ) { ANA_MSG_DEBUG( "event has fired stream: " << stream.name() ); } } // GRL if ( m_applyGRLCut ) { if ( !m_grl_handle->passRunLB( *eventInfo ) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; // go to next event } m_cutflowHist ->Fill( m_cutflow_grl, 1 ); m_cutflowHistW->Fill( m_cutflow_grl, mcEvtWeight); } //------------------------------------------------------------ // Apply event cleaning to remove events due to // problematic regions of the detector, and incomplete events. // Apply to data. //------------------------------------------------------------ if ( m_applyEventCleaningCut && (eventInfo->errorState(xAOD::EventInfo::LAr)==xAOD::EventInfo::Error ) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } m_cutflowHist ->Fill( m_cutflow_lar, 1 ); m_cutflowHistW->Fill( m_cutflow_lar, mcEvtWeight); if ( m_applyEventCleaningCut && (eventInfo->errorState(xAOD::EventInfo::Tile)==xAOD::EventInfo::Error ) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } m_cutflowHist ->Fill( m_cutflow_tile, 1 ); m_cutflowHistW->Fill( m_cutflow_tile, mcEvtWeight); if ( m_applyEventCleaningCut && (eventInfo->errorState(xAOD::EventInfo::SCT)==xAOD::EventInfo::Error) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } m_cutflowHist ->Fill( m_cutflow_SCT, 1 ); m_cutflowHistW->Fill( m_cutflow_SCT, mcEvtWeight); if ( m_applyCoreFlagsCut && (eventInfo->isEventFlagBitSet(xAOD::EventInfo::Core, 18) ) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } m_cutflowHist ->Fill( m_cutflow_core, 1 ); m_cutflowHistW->Fill( m_cutflow_core, mcEvtWeight); } // more info: https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/HowToCleanJets2017 if ( m_applyJetCleaningEventFlag && eventInfo->isAvailable("DFCommonJets_eventClean_LooseBad") ) { if(eventInfo->auxdataConst("DFCommonJets_eventClean_LooseBad")<1) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } } m_cutflowHist ->Fill( m_cutflow_jetcleaning, 1 ); m_cutflowHistW->Fill( m_cutflow_jetcleaning, mcEvtWeight); // n.b. this cut should only be applied in 2015+16 data, and not to MC! // details here: https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/HowToCleanJets2017#IsBadBatMan_Event_Flag_and_EMEC if ( m_applyIsBadBatmanFlag && eventInfo->isAvailable("DFCommonJets_isBadBatman") && !isMC() ) { if(eventInfo->auxdataConst("DFCommonJets_isBadBatman")>0) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } } m_cutflowHist ->Fill( m_cutflow_isbadbatman, 1 ); m_cutflowHistW->Fill( m_cutflow_isbadbatman, mcEvtWeight); //----------------------------- // Primary Vertex 'quality' cut //----------------------------- const xAOD::VertexContainer* vertices(nullptr); if ( !m_truthLevelOnly && m_applyPrimaryVertexCut ) { ANA_CHECK( HelperFunctions::retrieve(vertices, m_vertexContainerName, m_event, m_store, msg()) ); if ( !HelperFunctions::passPrimaryVertexSelection( vertices, m_PVNTrack ) ) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } } m_cutflowHist ->Fill( m_cutflow_npv, 1 ); m_cutflowHistW->Fill( m_cutflow_npv, mcEvtWeight); //--------------------- // Trigger decision cut //--------------------- if ( !m_triggerSelection.empty() || m_storeTrigDecisions ) { auto triggerChainGroup = m_trigDecTool_handle->getChainGroup(m_triggerSelection); if ( m_applyTriggerCut ) { // additional DEBUG logging to validate conditional logic ANA_MSG_DEBUG("Applying trigger cut corresponding to chain group " << m_triggerSelection); ANA_MSG_DEBUG("Is Trigger-Level Analysis (TLA) data = " << int(m_isTLAData)); ANA_MSG_DEBUG("Trigger chain group is passed = " << int(m_isTLAData ? triggerChainGroup->isPassed(TrigDefs::requireDecision) : triggerChainGroup->isPassed())); // different behaviour for isPassed depending on whether you are running on TLA data or not // if running on TLA data, we only store the HLT part of the trigger decision i.e. the L1 part // will always be "false", so we need to use TrigDefs::requireDecision to limit the decision // to being satisfied by the HLT leg(s) of the trigger chain // TODO: check performance of this method when using trigger chains with the SAME HLT leg but different L1 seed // e.g. HLT_j20_pf_ftf_L1J100 vs. HLT_j20_pf_ftf_L1HT190-J15s5pETA21 if ( (m_isTLAData && !triggerChainGroup->isPassed(TrigDefs::requireDecision)) || (!m_isTLAData && !triggerChainGroup->isPassed()) ) { // if (!triggerChainGroup->isPassed(TrigDefs::requireDecision)) { wk()->skipEvent(); return EL::StatusCode::SUCCESS; } m_cutflowHist ->Fill( m_cutflow_trigger, 1 ); m_cutflowHistW->Fill( m_cutflow_trigger, mcEvtWeight); } // save passed triggers in eventInfo // if ( m_storeTrigDecisions ) { std::vector passedTriggers; std::vector disabledTriggers; std::vector triggerPrescales; std::vector triggerPrescalesLumi; std::vector isPassedBitsNames; std::vector isPassedBits; // Save info for the triggers used to skim events // for ( auto &trigName : triggerChainGroup->getListOfTriggers() ) { auto trigChain = m_trigDecTool_handle->getChainGroup( trigName ); if ( (m_isTLAData && trigChain->isPassed(TrigDefs::requireDecision)) || (!m_isTLAData && trigChain->isPassed()) ) { passedTriggers.push_back( trigName ); triggerPrescales.push_back( trigChain->getPrescale() ); bool doLumiPrescale = std::find(m_triggerUnprescaleList.begin(), m_triggerUnprescaleList.end(), trigName) != m_triggerUnprescaleList.end(); if ( doLumiPrescale ) { triggerPrescalesLumi.push_back( m_pileup_tool_handle->getDataWeight( *eventInfo, trigName, true ) ); } else { triggerPrescalesLumi.push_back( -1 ); } } isPassedBitsNames.push_back( trigName ); isPassedBits .push_back( m_trigDecTool_handle->isPassedBits(trigName) ); if(trigChain->getPrescale()<1) disabledTriggers.push_back( trigName ); } // Save info for extra triggers // if ( !m_extraTriggerSelection.empty() ) { for ( const std::string &trigName : m_extraTriggerSelectionList ) { auto trigChain = m_trigDecTool_handle->getChainGroup( trigName ); if ( (m_isTLAData && trigChain->isPassed(TrigDefs::requireDecision)) || (!m_isTLAData && trigChain->isPassed()) ) { passedTriggers.push_back( trigName ); triggerPrescales.push_back( trigChain->getPrescale() ); bool doLumiPrescale = true; for ( const std::string &trigPart : trigChain->getListOfTriggers() ) { if (std::find(m_triggerUnprescaleList.begin(), m_triggerUnprescaleList.end(), trigPart) == m_triggerUnprescaleList.end()) doLumiPrescale = false; } if ( doLumiPrescale ) { triggerPrescalesLumi.push_back( m_pileup_tool_handle->getDataWeight( *eventInfo, trigName, true ) ); } else { triggerPrescalesLumi.push_back( -1 ); } } isPassedBitsNames.push_back( trigName ); isPassedBits .push_back( m_trigDecTool_handle->isPassedBits(trigName) ); if(trigChain->getPrescale()<1) disabledTriggers.push_back( trigName ); } } static SG::AuxElement::Decorator< std::vector< std::string > > dec_passedTriggers("passedTriggers"); dec_passedTriggers ( *eventInfo ) = passedTriggers; static SG::AuxElement::Decorator< std::vector< std::string > > dec_disabledTriggers("disabledTriggers"); dec_disabledTriggers( *eventInfo ) = disabledTriggers; static SG::AuxElement::Decorator< std::vector< float > > dec_triggerPrescales("triggerPrescales"); dec_triggerPrescales( *eventInfo ) = triggerPrescales; static SG::AuxElement::Decorator< std::vector< float > > dec_triggerPrescalesLumi("triggerPrescalesLumi"); dec_triggerPrescalesLumi( *eventInfo ) = triggerPrescalesLumi; static SG::AuxElement::Decorator< std::vector< unsigned int > > dec_isPassedBits("isPassedBits"); dec_isPassedBits( *eventInfo ) = isPassedBits; static SG::AuxElement::Decorator< std::vector< std::string > > dec_isPassedBitsNames("isPassedBitsNames"); dec_isPassedBitsNames( *eventInfo ) = isPassedBitsNames; } if ( m_storePrescaleWeight ) { static SG::AuxElement::Decorator< float > weight_prescale("weight_prescale"); weight_prescale(*eventInfo) = triggerChainGroup->getPrescale(); } if ( m_storePassL1 ) { static SG::AuxElement::Decorator< int > passL1("passL1"); passL1(*eventInfo) = ( m_triggerSelection.find("L1_") != std::string::npos ) ? (int)m_trigDecTool_handle->isPassed(m_triggerSelection.c_str()) : -1; } if ( m_storePassHLT ) { static SG::AuxElement::Decorator< int > passHLT("passHLT"); if (!m_isTLAData) { passHLT(*eventInfo) = ( m_triggerSelection.find("HLT_") != std::string::npos ) ? (int)m_trigDecTool_handle->isPassed(m_triggerSelection.c_str()) : -1; } else { passHLT(*eventInfo) = ( m_triggerSelection.find("HLT_") != std::string::npos ) ? (int)m_trigDecTool_handle->isPassed(m_triggerSelection.c_str(), TrigDefs::requireDecision) : -1; } } } // if giving a specific list of triggers to look at if ( m_storeTrigKeys ) { static SG::AuxElement::Decorator< int > masterKey("masterKey"); masterKey(*eventInfo) = m_trigConfTool_handle->masterKey(); static SG::AuxElement::Decorator< int > L1PSKey("L1PSKey"); L1PSKey(*eventInfo) = m_trigConfTool_handle->lvl1PrescaleKey(); static SG::AuxElement::Decorator< int > HLTPSKey("HLTPSKey"); HLTPSKey(*eventInfo) = m_trigConfTool_handle->hltPrescaleKey(); } // Calculate distance to previous empty BCID and previous unpaired BCID, and save as decorations if( m_calcBCIDInfo && !isMC() && m_trigConfTool_handle.isInitialized() ){ //Distance to previous empty BCID for (int i = eventInfo->bcid() - 1; i >= 0; i--){ //get the bunch group pattern for bunch crossing i uint16_t bgPattern = m_trigConfTool_handle->bunchGroupSet()->bgPattern()[i]; bool isEmpty = (bgPattern >> 3) & 0x1; if (isEmpty){ static SG::AuxElement::Decorator< int > DistEmptyBCID("DistEmptyBCID"); DistEmptyBCID(*eventInfo) = eventInfo->bcid()-i; break; } }//for each bcid //Distance to previous unpaired crossing for (int i = eventInfo->bcid() - 1; i >= 0; i--){ //get the bunch group pattern for bunch crossing i uint16_t bgPattern = m_trigConfTool_handle->bunchGroupSet()->bgPattern()[i]; bool isUnpaired = !((bgPattern >> 1) & 0x1); if (isUnpaired){ static SG::AuxElement::Decorator< int > DistLastUnpairedBCID("DistLastUnpairedBCID"); DistLastUnpairedBCID(*eventInfo) = eventInfo->bcid()-i; break; } }//for each bcid //Distance to next unpaired crossing for (int i = eventInfo->bcid() + 1; i <= 3654; i++){ //get the bunch group pattern for bunch crossing i uint16_t bgPattern = m_trigConfTool_handle->bunchGroupSet()->bgPattern()[i]; bool isUnpaired = !((bgPattern >> 1) & 0x1); if (isUnpaired){ static SG::AuxElement::Decorator< int > DistNextUnpairedBCID("DistNextUnpairedBCID"); DistNextUnpairedBCID(*eventInfo) = i-eventInfo->bcid(); break; } }// for each bcid }//if data //------------------------------------------------------------------------------------------ // Fill cutflows for run-by-run checks on data after cuts //------------------------------------------------------------------------------------------ if (m_doRunByRunCutflows && !isMC()) { m_runByrun_afterCuts->Fill(TString(std::to_string(runNumberForCutflow)), 1.); } return EL::StatusCode::SUCCESS; } // "Borrowed" from SUSYTools // https://gitlab.cern.ch/atlas/athena/blob/3be30397de7c6cfdc15de38f532fdb4b9f338297/PhysicsAnalysis/SUSYPhys/SUSYTools/Root/SUSYObjDef_xAOD.cxx#L700 // Should track official https://gitlab.cern.ch/atlas/athena/-/blob/d09ebd5c011b654159e25ec3e205c09ac1bc64de/PhysicsAnalysis/AnalysisCommon/PileupReweighting/python/AutoconfigurePRW.py StatusCode BasicEventSelection::autoconfigurePileupRWTool() { // Don't do this if we aren't supposed to if (! (isMC() && m_autoconfigPRW )) return StatusCode::SUCCESS; const xAOD::EventInfo* eventInfo = 0; ANA_CHECK( m_event->retrieve( eventInfo, "EventInfo" ) ); // Determine simulation flavour std::string SimulationFlavour = isFastSim() ? ( isAF3() ? "AF3" : "AFII" ) : "FS"; // Extract campaign automatically from Run Number std::string mcCampaignMD = ""; uint32_t runNum = eventInfo->runNumber(); switch(runNum) { case 284500 : mcCampaignMD="mc20a"; break; case 300000 : mcCampaignMD="mc20d"; break; case 310000 : mcCampaignMD="mc20e"; break; case 410000 : mcCampaignMD="mc23a"; break; case 450000 : mcCampaignMD="mc23d"; break; case 470000 : mcCampaignMD="mc23e"; break; default : ANA_MSG_ERROR( "Could not determine mc campaign from run number! Impossible to autoconfigure PRW. Aborting." ); return StatusCode::FAILURE; break; } std::string mcCampaignMD_v2 = ""; const xAOD::FileMetaData* fmd(nullptr); if ( !m_event->retrieveMetaInput(fmd, "FileMetaData").isSuccess() || !fmd->value(xAOD::FileMetaData::mcCampaign, mcCampaignMD_v2) ) { ANA_MSG_WARNING("Failed to retrieve FileMetaData from MetaData! Using MC campaign from run number. PLEASE DOUBLE-CHECK this is the correct campaign for your samples!"); } else { fmd->value(xAOD::FileMetaData::mcCampaign, mcCampaignMD_v2); if(mcCampaignMD!=mcCampaignMD_v2){ std::string MetadataAndRunConflict(""); MetadataAndRunConflict += "autoconfigurePileupRWTool(): access to FileMetaData indicates a " + mcCampaignMD_v2; MetadataAndRunConflict += " sample, but the run number indiciates " +mcCampaignMD; MetadataAndRunConflict += ". Prioritizing the value from FileMetaData. This could occur if you are using an MC campaign with outdated pile-up reweighting. PLEASE DOUBLE-CHECK your samples!"; ANA_MSG_WARNING( MetadataAndRunConflict ); mcCampaignMD=mcCampaignMD_v2; } } ANA_MSG_INFO( "Determined MC campaign to be " << mcCampaignMD); // Extract campaign from user configuration std::string tmp_mcCampaign = m_mcCampaign; std::vector mcCampaignList; while ( tmp_mcCampaign.size() > 0) { size_t pos = tmp_mcCampaign.find_first_of(','); if ( pos == std::string::npos ) { pos = tmp_mcCampaign.size(); mcCampaignList.push_back(tmp_mcCampaign.substr(0, pos)); tmp_mcCampaign.erase(0, pos); } else { mcCampaignList.push_back(tmp_mcCampaign.substr(0, pos)); tmp_mcCampaign.erase(0, pos+1); } } // Sanity checks bool mc2XX_GoodFromProperty = !mcCampaignList.empty(); bool mc2XX_GoodFromMetadata = false; for(const auto& mcCampaignP : mcCampaignList) mc2XX_GoodFromProperty &= ( mcCampaignP == "mc20a" || mcCampaignP == "mc20d" || mcCampaignP == "mc20e" || mcCampaignP == "mc23a" || mcCampaignP == "mc23c" || mcCampaignP == "mc23d" || mcCampaignP == "mc23e"); if( mcCampaignMD == "mc20a" || mcCampaignMD == "mc20d" || mcCampaignMD == "mc20e" || mcCampaignMD == "mc23a" || mcCampaignMD == "mc23c" || mcCampaignMD == "mc23d" || mcCampaignMD == "mc23e") mc2XX_GoodFromMetadata = true; if( !mc2XX_GoodFromMetadata && !mc2XX_GoodFromProperty ) { // :: std::string MetadataAndPropertyBAD(""); MetadataAndPropertyBAD += "autoconfigurePileupRWTool(): access to FileMetaData failed, but don't panic. You can try to manually set the 'mcCampaign' BasicEventSelection property to "; MetadataAndPropertyBAD += "'mc20a', 'mc20c', 'mc20d', 'mc20e', 'mc20f', 'mc23a', 'mc23c', 'mc23d', or 'mc23e' and restart your job. If you set it to any other string, you will still incur in this error."; ANA_MSG_ERROR( MetadataAndPropertyBAD ); return StatusCode::FAILURE; // :: } if ( mc2XX_GoodFromProperty && mc2XX_GoodFromMetadata) { bool MDinP=false; for(const auto& mcCampaignP : mcCampaignList) MDinP |= (mcCampaignMD==mcCampaignP); if( !MDinP ) { // :: std::string MetadataAndPropertyConflict(""); MetadataAndPropertyConflict += "autoconfigurePileupRWTool(): access to FileMetaData indicates a " + mcCampaignMD; MetadataAndPropertyConflict += " sample, but the 'mcCampaign' property passed to BasicEventSelection is set to '" +m_mcCampaign; MetadataAndPropertyConflict += "'. Prioritizing the value set by user: PLEASE DOUBLE-CHECK the value you set the 'mcCampaign' property to!"; ANA_MSG_WARNING( MetadataAndPropertyConflict ); // :: } else { // :: std::string NoMetadataButPropertyOK(""); NoMetadataButPropertyOK += "autoconfigurePileupRWTool(): access to FileMetaData succeeded, but the 'mcCampaign' property is passed to BasicEventSelection as '"; NoMetadataButPropertyOK += m_mcCampaign; NoMetadataButPropertyOK += "'. Autoconfiguring PRW accordingly."; ANA_MSG_WARNING( NoMetadataButPropertyOK ); // :: } } // :: // Retrieve the input file if(!mc2XX_GoodFromProperty) { mcCampaignList.clear(); mcCampaignList.push_back(mcCampaignMD); } ANA_MSG_INFO( "Setting MC campgains for CP::PileupReweightingTool:"); for(const auto& mcCampaign : mcCampaignList) ANA_MSG_INFO( "\t" << mcCampaign.c_str() ); // float dsid = -999; dsid = eventInfo->mcChannelNumber(); int DSID_INT = (int) dsid; std::vector prwConfigFiles; for(const auto& mcCampaign : mcCampaignList) { std::string prwConfigFile; // If requested set the PRW file to common PRW file of the processed MC campaign if (m_useCommonPRWFiles) { if (mcCampaign == "mc20a") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC20a);} else if (mcCampaign == "mc20d") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC20d);} else if (mcCampaign == "mc20e") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC20e);} else if (mcCampaign == "mc23a") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC23a);} else if (mcCampaign == "mc23c") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC23c);} else if (mcCampaign == "mc23d") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC23d);} else if (mcCampaign == "mc23e") {prwConfigFile = PathResolverFindCalibFile(m_commonPRWFileMC23e);} else { ANA_MSG_ERROR("autoconfigurePileupRWTool(): no common PRW file known for MC campaign: " << mcCampaign); return StatusCode::FAILURE; } } else { prwConfigFile = PathResolverFindCalibFile("dev/PileupReweighting/share/DSID" + std::to_string(DSID_INT/1000) +"xxx/pileup_" + mcCampaign + "_dsid" + std::to_string(DSID_INT) + "_" + SimulationFlavour + ".root"); } TFile testF(prwConfigFile.data(),"read"); if(testF.IsZombie()) { ANA_MSG_ERROR("autoconfigurePileupRWTool(): Missing PRW config file for DSID " << std::to_string(DSID_INT) << " in campaign " << mcCampaign); return StatusCode::FAILURE; } else prwConfigFiles.push_back( prwConfigFile ); } // Add actualMu config files for(const auto& mcCampaign : mcCampaignList) { if( !m_prwActualMu2016File.empty() && mcCampaign == "mc20a" ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2016File)); if( !m_prwActualMu2017File.empty() && (mcCampaign == "mc20c" || mcCampaign=="mc20d") ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2017File)); if( !m_prwActualMu2018File.empty() && (mcCampaign == "mc20e" || mcCampaign=="mc20f") ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2018File)); if( !m_prwActualMu2022File.empty() && mcCampaign == "mc23a" ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2022File)); if( !m_prwActualMu2023File.empty() && (mcCampaign == "mc23c" || mcCampaign=="mc23d") ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2023File)); if( !m_prwActualMu2024File.empty() && (mcCampaign == "mc23e") ) prwConfigFiles.push_back(PathResolverFindCalibFile(m_prwActualMu2024File)); } // also need to handle lumicalc files: only use 2015+2016 with mc20a // and only use 2017 with mc20d and 2018 data with mc20e // according to instructions on https://twiki.cern.ch/twiki/bin/view/AtlasProtected/ExtendedPileupReweighting#Tool_Properties // Parse lumicalc file names std::vector allLumiCalcFiles; std::string tmp_lumiCalcFileNames = m_lumiCalcFileNames; while ( tmp_lumiCalcFileNames.size() > 0) { size_t pos = tmp_lumiCalcFileNames.find_first_of(','); if ( pos == std::string::npos ) { pos = tmp_lumiCalcFileNames.size(); allLumiCalcFiles.push_back(tmp_lumiCalcFileNames.substr(0, pos)); tmp_lumiCalcFileNames.erase(0, pos); } else { allLumiCalcFiles.push_back(tmp_lumiCalcFileNames.substr(0, pos)); tmp_lumiCalcFileNames.erase(0, pos+1); } } std::vector lumiCalcFiles; for(const auto& mcCampaign : mcCampaignList) { for(const auto& filename : allLumiCalcFiles) { // looking for things of format "stuff/data15_13TeV/stuff" etc size_t pos = filename.find("data"); std::string year = filename.substr(pos+4, 2); if (mcCampaign == "mc20a") { if (year == "15" || year == "16") { lumiCalcFiles.push_back(filename); } } else if (mcCampaign == "mc20d") { if (year == "17") { lumiCalcFiles.push_back(filename); } } else if (mcCampaign == "mc20e") { if (year == "18") { lumiCalcFiles.push_back(filename); } } else if (mcCampaign == "mc23a") { if (year == "22") { lumiCalcFiles.push_back(filename); } } else if (mcCampaign == "mc23c" || mcCampaign == "mc23d") { if (year == "23") { lumiCalcFiles.push_back(filename); } } else if (mcCampaign == "mc23e") { if (year == "24") { lumiCalcFiles.push_back(filename); } } else { ANA_MSG_ERROR( "No lumicalc file is suitable for your mc campaign!" ); } } } // Set everything and report on it. ANA_MSG_INFO( "Adding Pileup files for CP::PileupReweightingTool:"); for( unsigned int i=0; i < prwConfigFiles.size(); ++i) { ANA_MSG_INFO("\t- " << prwConfigFiles.at(i).c_str()); } ANA_CHECK( m_pileup_tool_handle.setProperty("ConfigFiles", prwConfigFiles)); ANA_MSG_INFO( "Adding LumiCalc files for CP::PileupReweightingTool:"); for( unsigned int i=0; i < lumiCalcFiles.size(); ++i) { ANA_MSG_INFO("\t- " << lumiCalcFiles.at(i).c_str()); } ANA_CHECK( m_pileup_tool_handle.setProperty("LumiCalcFiles", lumiCalcFiles)); // Return gracefully return StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: postExecute () { // Here you do everything that needs to be done after the main event // processing. This is typically very rare, particularly in user // code. It is mainly used in implementing the NTupleSvc. return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: finalize () { // This method is the mirror image of initialize(), meaning it gets // called after the last event has been processed on the worker node // and allows you to finish up any objects you created in // initialize() before they are written to disk. This is actually // fairly rare, since this happens separately for each worker node. // Most of the time you want to do your post-processing on the // submission node after all your histogram outputs have been // merged. This is different from histFinalize() in that it only // gets called on worker nodes that processed input events. ANA_MSG_INFO( "Number of processed events \t= " << m_eventCounter); m_RunNr_VS_EvtNr.clear(); if ( m_trigDecTool_handle.isInitialized() ){ if (asg::ToolStore::contains("ToolSvc.TrigDecisionTool") ){ m_trigDecTool_handle->finalize(); asg::ToolStore::remove("ToolSvc.TrigDecisionTool").ignore(); } } //after execution loop if(m_printBranchList){ xAOD::IOStats::instance().stats().printSmartSlimmingBranchList(); } return EL::StatusCode::SUCCESS; } EL::StatusCode BasicEventSelection :: histFinalize () { // This method is the mirror image of histInitialize(), meaning it // gets called after the last event has been processed on the worker // node and allows you to finish up any objects you created in // histInitialize() before they are written to disk. This is // actually fairly rare, since this happens separately for each // worker node. Most of the time you want to do your // post-processing on the submission node after all your histogram // outputs have been merged. This is different from finalize() in // that it gets called on all worker nodes regardless of whether // they processed input events. ANA_CHECK( xAH::Algorithm::algFinalize()); return EL::StatusCode::SUCCESS; }