diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h index a5c1a56ac4cb9042ba599e27aaf1d54ca1a775a2..ae88412a4e3b3b6b83fb2b4f0862d2d3cb772c7f 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JGTowerReader.h @@ -53,7 +53,6 @@ class JGTowerReader: public ::AthAlgorithm { bool m_debugJetAlg; bool m_dumpTowerInfo; bool m_dumpSeedsEtaPhi; - bool m_makeRoundJetsPUsub; bool m_makeSquareJets; bool m_buildgBlockJets; @@ -98,22 +97,13 @@ class JGTowerReader: public ::AthAlgorithm { bool m_plotSeeds; bool m_saveSeeds; - float m_gJet_seed_size; - float m_gJet_max_r; float m_gJet_r; - float m_gJet_seed_tower_noise_multiplier; - float m_gJet_seed_total_noise_multiplier; - float m_gJet_seed_min_ET_MeV; - float m_gJet_jet_tower_noise_multiplier; - float m_gJet_jet_total_noise_multiplier; + float m_gJet_block_tower_noise_multiplier; + float m_gJet_block_min_ET_MeV; + float m_gJet_tower_noise_multiplier; float m_gJet_jet_min_ET_MeV; std::string m_noise_file; - - std::string m_jXERHO_correction_file; - float m_jXERHO_fixed_noise_cut; - float m_jXERHO_rho_up_threshold; - float m_jXERHO_min_noise_cut; //job options for gFEX MET algorithms bool m_useRMS; @@ -137,10 +127,6 @@ class JGTowerReader: public ::AthAlgorithm { virtual StatusCode BuildJetsFromMap(const xAOD::JGTowerContainer*jTs); std::vector jT_noise; - std::vector jTowerArea; - std::vector < std::vector < int > > jFEX_bins; - std::vector < std::vector < int > > jFEX_bins_core; - bool buildbins=false; std::vector jJet_thr; std::vector jJet_jet_thr; std::vector gT_noise; diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JetAlg.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JetAlg.h index 6ac8f608eb90c53796386de9cc2ff4c78a882851..16e677217a9943875e4dd31af4fe3b02a08a4d6a 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JetAlg.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/JetAlg.h @@ -60,9 +60,9 @@ class JetAlg{ static StatusCode SeedGrid(const xAOD::JGTowerContainer*towers, TString seedname, bool &m_dumpSeedsEtaPhi); static StatusCode SeedFinding(const xAOD::JGTowerContainer*towers, TString seedname, float seed_size, float range, std::vector noise, float seed_tower_noise_multiplier, float seed_total_noise_multiplier, float seed_min_ET_MeV); - static StatusCode BuildFatJet(const xAOD::JGTowerContainer towers, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float jet_total_noise_multiplier, float jet_min_ET_MeV, float rho); + static StatusCode BuildFatJet(const xAOD::JGTowerContainer towers, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float seed_threshold, float jet_min_ET_MeV, float rhoA, float rhoB, bool useNegTowers); static StatusCode BuildJet(const xAOD::JGTowerContainer*towers, TString seedname, TString jetname, float jet_size, std::vector noise, float jet_tower_noise_multiplier, float jet_total_noise_multiplier, float jet_min_ET_MeV, bool m_saveSeeds = false); - static StatusCode BuildRoundJet(const xAOD::JGTowerContainer*towers, TString seedname, TString jetname, float jet_size, std::vector noise, float jet_tower_noise_multiplier, float jet_total_noise_multiplier, float jet_min_ET_MeV, bool m_saveSeeds = false); - static StatusCode BuildgBlocksJets(const xAOD::JGTowerContainer* blocks, TString jetname, float rho); + static StatusCode BuildRoundJet(const xAOD::JGTowerContainer*towers, TString seedname, TString jetname, float jet_size, std::vector noise, float jet_tower_noise_multiplier, float jet_min_ET_MeV, bool m_saveSeeds = false); + static StatusCode BuildgBlocksJets(const xAOD::JGTowerContainer* blocks, TString jetname, float rhoA, float rhoB); }; #endif diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h index 33567280fc2d663d4748e365f8d49e5f6fc68997..6f8a17ed9741886a53983fb1f948d7f6626b558e 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h +++ b/Trigger/TrigT1/TrigT1CaloFexSim/TrigT1CaloFexSim/METAlg.h @@ -33,16 +33,13 @@ #include "Objects.h" class METAlg{ - private: - std::vector < std::vector < int > > m_jFEX_bins; - std::vector < std::vector < int > > m_jFEX_bins_core; - bool m_buildbins=false; - public: + + public: struct MET{ float ex; float ey; - float phi; + //float phi; float et; float rho = 0; float mht = 0; @@ -54,7 +51,7 @@ class METAlg{ float scalar_Et = 0; }; - static std::map> m_METMap; + static std::map> m_METMap; /** *@brief Calculate MET using a fixed 4 sigma noise cut */ @@ -63,12 +60,6 @@ class METAlg{ *@brief Calculates MET with pileup subtraction */ static StatusCode SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metname, bool useEtaBins, bool useRMS, bool useMedian, bool useNegTowers); - /** - *@brief Calculates MET with pileup subtraction in jFEX - */ - static std::vector check_in_bin (const float &eta, const float &phi, const std::vector < std::pair < float, float> > &eta_bins, const std::vector < std::pair < float, float> > &phi_bins, const float &phi_offset); - static StatusCode build_jFEX_bins( std::vector < std::vector < int > > &bins, std::vector < std::vector < int > > &bins_core, const xAOD::JGTowerContainer* towers ); - static StatusCode jXERHO(const xAOD::JGTowerContainer* towers, TString metName, const std::vector jTowerArea, const std::vector < std::vector < int > > jFEX_bins, const std::vector < std::vector < int > > jFEX_bins_core, float fixed_noise_cut, float rho_up_threshold, float min_noise_cut, xAOD::JGTowerContainer* towers_PUsub ); /** *@brief Calculates MET with Softkiller */ @@ -76,7 +67,7 @@ class METAlg{ /** *@brief Calculates MET with Jets without Jets */ - static StatusCode JwoJ_MET(const xAOD::JGTowerContainer* towers, const std::vector gBlocks, TString metname, float pTcone_cut, bool useEtaBins, bool useRho, bool useNegTowers); + static StatusCode JwoJ_MET(const xAOD::JGTowerContainer* towers, const std::vector gBlocks, TString metname, float pTcone_cut, bool useEtaBins, bool useRho, bool useNegTowers); /** *@brief Calculates MET using PUfit */ diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/python/TrigT1CaloFexSimConfig.py b/Trigger/TrigT1/TrigT1CaloFexSim/python/TrigT1CaloFexSimConfig.py index b2c43d0e0caf32447b11a47f1be55c751a1a4ec9..3c063fc134e90578a2e162dabce3f6efd9a396ea 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/python/TrigT1CaloFexSimConfig.py +++ b/Trigger/TrigT1/TrigT1CaloFexSim/python/TrigT1CaloFexSimConfig.py @@ -56,8 +56,6 @@ def createJGTowerReader( SuperCellType = "SCell", **kwargs ) : jJetRound_jet_total_noise_multiplier = 0.0, jJetRound_jet_min_ET_MeV = 5000, - makeRoundJetsPUsub=True, - makeRoundLargeRJets = True, jJetRound_LargeR_seed_size = 0.31, # seed square of side this. 0.3 for 3x3 towers jJetRound_LargeR_max_r = 0.26, # distance (in eta and phi, not a radius) within @@ -75,12 +73,12 @@ def createJGTowerReader( SuperCellType = "SCell", **kwargs ) : map_jet_total_noise_multiplier = 0.0, map_jet_min_ET_MeV = 500, - gJet_seed_size=0.2,#Not for gFEX - gJet_max_r=1.0, - gJet_r=1.0, - gJet_seed_tower_noise_multiplier = 0.0, - gJet_seed_min_ET_MeV = 500, + gJet_r=0.9, + gJet_block_tower_noise_multiplier = 0., + gJet_block_min_ET_MeV = 0., + gJet_tower_noise_multiplier = 0., gJet_jet_min_ET_MeV = 5000., + useNegTowers = True ) for arg in kwargs: towerReader.__setattr__(arg,kwargs[arg]) diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.cxx new file mode 100644 index 0000000000000000000000000000000000000000..aa4dbdc8c08c3321fc943652c32081a233b6c3c2 --- /dev/null +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.cxx @@ -0,0 +1,93 @@ +// TrigT1CaloFexSim includes +#include "JFexEleTau.h" + +//#include "xAODEventInfo/EventInfo.h" + + + + +JFexEleTau::JFexEleTau( const std::string& name, ISvcLocator* pSvcLocator ) : AthAnalysisAlgorithm( name, pSvcLocator ){ + + //declareProperty( "Property", m_nProperty = 0, "My Example Integer Property" ); //example property declaration + +} + + +JFexEleTau::~JFexEleTau() {} + + +StatusCode JFexEleTau::initialize() { + ATH_MSG_INFO ("Initializing " << name() << "..."); + // + //This is called once, before the start of the event loop + //Retrieves of tools you have configured in the joboptions go here + // + + //HERE IS AN EXAMPLE + //We will create a histogram and a ttree and register them to the histsvc + //Remember to configure the histsvc stream in the joboptions + // + //m_myHist = new TH1D("myHist","myHist",100,0,100); + //CHECK( histSvc()->regHist("/MYSTREAM/myHist", m_myHist) ); //registers histogram to output stream + //m_myTree = new TTree("myTree","myTree"); + //CHECK( histSvc()->regTree("/MYSTREAM/SubDirectory/myTree", m_myTree) ); //registers tree to output stream inside a sub-directory + + + return StatusCode::SUCCESS; +} + +StatusCode JFexEleTau::finalize() { + ATH_MSG_INFO ("Finalizing " << name() << "..."); + // + //Things that happen once at the end of the event loop go here + // + + + return StatusCode::SUCCESS; +} + +StatusCode JFexEleTau::execute() { + ATH_MSG_DEBUG ("Executing " << name() << "..."); + setFilterPassed(false); //optional: start with algorithm not passed + + + + // + //Your main analysis code goes here + //If you will use this algorithm to perform event skimming, you + //should ensure the setFilterPassed method is called + //If never called, the algorithm is assumed to have 'passed' by default + // + + + //HERE IS AN EXAMPLE + //const xAOD::EventInfo* ei = 0; + //CHECK( evtStore()->retrieve( ei , "EventInfo" ) ); + //ATH_MSG_INFO("eventNumber=" << ei->eventNumber() ); + //m_myHist->Fill( ei->averageInteractionsPerCrossing() ); //fill mu into histogram + + + setFilterPassed(true); //if got here, assume that means algorithm passed + return StatusCode::SUCCESS; +} + +StatusCode JFexEleTau::beginInputFile() { + // + //This method is called at the start of each input file, even if + //the input file contains no events. Accumulate metadata information here + // + + //example of retrieval of CutBookkeepers: (remember you will need to include the necessary header files and use statements in requirements file) + // const xAOD::CutBookkeeperContainer* bks = 0; + // CHECK( inputMetaStore()->retrieve(bks, "CutBookkeepers") ); + + //example of IOVMetaData retrieval (see https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/AthAnalysisBase#How_to_access_file_metadata_in_C) + //float beamEnergy(0); CHECK( retrieveMetadata("/TagInfo","beam_energy",beamEnergy) ); + //std::vector bunchPattern; CHECK( retrieveMetadata("/Digitiation/Parameters","BeamIntensityPattern",bunchPattern) ); + + + + return StatusCode::SUCCESS; +} + + diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.h b/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.h new file mode 100644 index 0000000000000000000000000000000000000000..06225f1e51db62b251d57952f7c91acb2412dd95 --- /dev/null +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JFexEleTau.h @@ -0,0 +1,49 @@ +#ifndef TRIGT1CALOFEXSIM_JFEXELETAU_H +#define TRIGT1CALOFEXSIM_JFEXELETAU_H 1 + +#include "AthAnalysisBaseComps/AthAnalysisAlgorithm.h" + +//Example ROOT Includes +//#include "TTree.h" +//#include "TH1D.h" + + + +class JFexEleTau: public ::AthAnalysisAlgorithm { + public: + JFexEleTau( const std::string& name, ISvcLocator* pSvcLocator ); + virtual ~JFexEleTau(); + + ///uncomment and implement methods as required + + //IS EXECUTED: + virtual StatusCode initialize(); //once, before any input is loaded + virtual StatusCode beginInputFile(); //start of each input file, only metadata loaded + //virtual StatusCode firstExecute(); //once, after first eventdata is loaded (not per file) + virtual StatusCode execute(); //per event + //virtual StatusCode endInputFile(); //end of each input file + //virtual StatusCode metaDataStop(); //when outputMetaStore is populated by MetaDataTools + virtual StatusCode finalize(); //once, after all events processed + + + ///Other useful methods provided by base class are: + ///evtStore() : ServiceHandle to main event data storegate + ///inputMetaStore() : ServiceHandle to input metadata storegate + ///outputMetaStore() : ServiceHandle to output metadata storegate + ///histSvc() : ServiceHandle to output ROOT service (writing TObjects) + ///currentFile() : TFile* to the currently open input file + ///retrieveMetadata(...): See twiki.cern.ch/twiki/bin/view/AtlasProtected/AthAnalysisBase#ReadingMetaDataInCpp + + + private: + + //Example algorithm property, see constructor for declaration: + //int m_nProperty = 0; + + //Example histogram, see initialize method for registration to output histSvc + //TH1D* m_myHist = 0; + //TTree* m_myTree = 0; + +}; + +#endif //> !TRIGT1CALOFEXSIM_JFEXELETAU_H diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTower.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTower.cxx index fcfa0b280e0e05581c091432b456689248e2b757..2d6726f40e71993d9e6fb1d02a4fb9ec62de7651 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTower.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTower.cxx @@ -42,11 +42,11 @@ bool withinRadius(float eta1, float eta2, float phi1, float phi2, float dR, bool //https://root.cern.ch/doc/master/TLorentzVector_8h_source.html#l00456 float dPhi = deltaPhi(phi1, phi2); float dEta = eta1-eta2; - float dR2 = dPhi*dPhi + dEta*dEta; + float test_dR = TMath::Sqrt(dPhi*dPhi + dEta*dEta); if(acceptEqual) - return dR2 <= dR*dR; + return test_dR <= dR; else - return dR2 < dR*dR; + return test_dR < dR; } float deltaPhi(float phi1,float phi2){ diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx index 0cbc928d3c4d0228d87f56c7071b0661f865425a..96e37145b04468b75d5f948be269213dbf1ffcfe 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JGTowerReader.cxx @@ -50,7 +50,6 @@ histSvc("THistSvc",name){ declareProperty("dumpTowerInfo", m_dumpTowerInfo=false); declareProperty("dumpSeedsEtaPhi", m_dumpSeedsEtaPhi=false); declareProperty("noise_file", m_noise_file="Run3L1CaloSimulation/Noise/noise_r10684.root"); - declareProperty("makeRoundJetsPUsub", m_makeRoundJetsPUsub=false); declareProperty("makeSquareJets", m_makeSquareJets = true); declareProperty("jJet_seed_size", m_jJet_seed_size=0.3); @@ -79,7 +78,7 @@ histSvc("THistSvc",name){ declareProperty("jJetRound_LargeR_max_r", m_jJetRound_LargeR_max_r=0.4); declareProperty("jJetRound_LargeR_r", m_jJetRound_LargeR_r=0.8); declareProperty("jJetRound_LargeR_seed_min_ET_MeV", m_jJetRound_LargeR_seed_min_ET_MeV = 2000.0); - declareProperty("jJetRound_LargeR_jet_min_ET_MeV", m_jJetRound_jet_min_ET_MeV = 2000.0); + declareProperty("jJetRound_LargeR_jet_min_ET_MeV", m_jJetRound_LargeR_jet_min_ET_MeV = 2000.0); declareProperty("makeJetsFromMap", m_makeJetsFromMap = false); declareProperty("towerMap", m_towerMap = ""); @@ -94,14 +93,10 @@ histSvc("THistSvc",name){ declareProperty("saveSeeds", m_saveSeeds = false); declareProperty("buildgBlockJets", m_buildgBlockJets=true); - declareProperty("gJet_seed_size", m_gJet_seed_size=0.2); - declareProperty("gJet_max_r", m_gJet_max_r=1.0); //gFEX constructs large radius jets declareProperty("gJet_r", m_gJet_r=1.0); - declareProperty("gJet_seed_tower_noise_multiplier", m_gJet_seed_tower_noise_multiplier = 1.0); - declareProperty("gJet_seed_total_noise_multiplier", m_gJet_seed_total_noise_multiplier = 1.0); - declareProperty("gJet_seed_min_ET_MeV", m_gJet_seed_min_ET_MeV = 2000.0); - declareProperty("gJet_jet_tower_noise_multiplier", m_gJet_jet_tower_noise_multiplier = 1.0); - declareProperty("gJet_jet_total_noise_multiplier", m_gJet_jet_total_noise_multiplier = 0.0); + declareProperty("gJet_block_tower_noise_multiplier", m_gJet_block_tower_noise_multiplier = 0.0); + declareProperty("gJet_block_min_ET_MeV", m_gJet_block_min_ET_MeV = 0.0); + declareProperty("gJet_tower_noise_multiplier", m_gJet_tower_noise_multiplier = 0.0); declareProperty("gJet_jet_min_ET_MeV", m_gJet_jet_min_ET_MeV = 2000.0); declareProperty("developerMET",m_developerMET=false); @@ -109,11 +104,6 @@ histSvc("THistSvc",name){ declareProperty("useMedian", m_useMedian=false); declareProperty("useNegTowers", m_useNegTowers=false); declareProperty("pTcone_cut", m_pTcone_cut=25); //cone threshold for Jets without Jets: declared in GeV - - declareProperty("jXERHO_correction_file" , m_jXERHO_correction_file="Run3L1CaloSimulation/Noise/jTowerCorrection.20200302.r11364.root"); //correction file for jXERHO - declareProperty("jXERHO_fixed_noise_cut" , m_jXERHO_fixed_noise_cut=0.0); - declareProperty("jXERHO_rho_up_threshold" , m_jXERHO_rho_up_threshold=1000.0); - declareProperty("jXERHO_min_noise_cut" , m_jXERHO_min_noise_cut=100.0); } @@ -163,31 +153,6 @@ StatusCode JGTowerReader::initialize() { gT_noise.push_back( gh_noise->GetBinContent(i+1) ); } } - - std::string fullPathTo_jXERHO_correction_file = PathResolverFindCalibFile(m_jXERHO_correction_file); - std::ifstream jXERHO_correction_exist(fullPathTo_jXERHO_correction_file.c_str()); - jTowerArea.clear(); - jTowerArea.resize(7712,1.0); - if(jXERHO_correction_exist){ - ATH_MSG_INFO ("jXERHO_correction_file is set with root file:" << fullPathTo_jXERHO_correction_file << "..."); - TFile *correctionfile = new TFile(fullPathTo_jXERHO_correction_file.c_str()); - TH1F *jTowerArea_final_hist = (TH1F*)correctionfile->Get("jTowerArea_final_hist"); - for (size_t i=0; i < 7712; i++) { - jTowerArea[i]=jTowerArea_final_hist->GetBinContent(i+1); - if (jTowerArea[i] == 0) { - ATH_MSG_DEBUG ("jTowerArea should never be zero, this concerns tower "<& std::vector neighbors = grid.neighbors(*seed, blockRows, blockCols); float seed_Et = seed->et(); - if(!useNegTowers) seed_Et = TMath::Abs(seed_Et); + if(!useNegTowers && seed->et() < 0)continue; + //if(!useNegTowers) seed_Et = TMath::Abs(seed_Et); double block_area(0.0); double sum_deta(0.0); double sum_dphi(0.0); @@ -401,7 +367,8 @@ StatusCode JGTowerReader::BuildBlocksFromTowers(std::vector& sum_deta += neighbor->deta(); sum_dphi += neighbor->dphi(); neighbor_pt = neighbor->et(); - if(!useNegTowers) neighbor_pt = TMath::Abs(neighbor_pt); + //if(!useNegTowers) neighbor_pt = TMath::Abs(neighbor_pt); + if(!useNegTowers && neighbor->et() < 0)continue; block_pt += neighbor_pt; } @@ -427,37 +394,7 @@ StatusCode JGTowerReader::JFexAlg(const xAOD::JGTowerContainer* jTs){ ATH_MSG_DEBUG("Found " << jTs->size() << " jTowers"); - ATH_MSG_DEBUG("JFexAlg: BuildMET"); - CHECK(METAlg::Baseline_MET(jTs, "jNOISECUT", jT_noise, m_useNegTowers)); - if (!buildbins) { - CHECK(METAlg::build_jFEX_bins(jFEX_bins, jFEX_bins_core, jTs )); - buildbins=true; - } - - xAOD::JGTowerContainer* jTs_PUsub = new xAOD::JGTowerContainer(); - xAOD::JGTowerAuxContainer* jTs_PUsubAux = new xAOD::JGTowerAuxContainer(); - jTs_PUsub->setStore(jTs_PUsubAux); - CHECK(METAlg::jXERHO(jTs , "jXERHO" , jTowerArea , jFEX_bins , jFEX_bins_core , m_jXERHO_fixed_noise_cut , m_jXERHO_rho_up_threshold , m_jXERHO_min_noise_cut , jTs_PUsub)); - ATH_MSG_DEBUG("JFexAlg: Done"); - if (m_makeRoundJetsPUsub) { - ATH_MSG_DEBUG("JFexAlg: JetsPUsub"); - if(m_makeRoundJets) { - ATH_MSG_DEBUG("JFexAlg: JetsPUsub, makeRoundJets"); - if( JetAlg::m_SeedMap.find("jRoundSeedsPUsub") == JetAlg::m_SeedMap.end() ) - CHECK( JetAlg::SeedGrid(jTs_PUsub, "jRoundSeedsPUsub", m_dumpSeedsEtaPhi) ); - - ATH_MSG_DEBUG("JFexAlg using JetsPUsub: SeedFinding with jJetSeeds; m_jJet_seed_size = " << m_jJet_seed_size << ", m_jJet_max_r = " << m_jJet_max_r); - CHECK( JetAlg::SeedFinding( jTs_PUsub, "jRoundSeedsPUsub", m_jJetRound_seed_size, m_jJetRound_max_r, jT_noise, - m_jJetRound_seed_tower_noise_multiplier, m_jJetRound_seed_total_noise_multiplier, - m_jJetRound_seed_min_ET_MeV) ); - - ATH_MSG_DEBUG("JFexAlg usinf JetsPUsub: BuildRoundJet"); - CHECK( JetAlg::BuildRoundJet(jTs_PUsub, "jRoundSeedsPUsub", "jRoundJetsPUsub", m_jJetRound_r, jT_noise, m_jJetRound_jet_tower_noise_multiplier, m_jJetRound_jet_total_noise_multiplier, m_jJetRound_jet_min_ET_MeV, m_saveSeeds) ); - } - delete jTs_PUsub; - delete jTs_PUsubAux; - } if(m_makeSquareJets) { // find all seeds @@ -487,18 +424,20 @@ StatusCode JGTowerReader::JFexAlg(const xAOD::JGTowerContainer* jTs){ m_jJetRound_seed_min_ET_MeV) ); ATH_MSG_DEBUG("JFexAlg: BuildRoundJet"); - CHECK( JetAlg::BuildRoundJet(jTs, "jRoundSeeds", "jRoundJets", m_jJetRound_r, jT_noise, m_jJetRound_jet_tower_noise_multiplier, m_jJetRound_jet_total_noise_multiplier, m_jJetRound_jet_min_ET_MeV, m_saveSeeds) ); + CHECK( JetAlg::BuildRoundJet(jTs, "jRoundSeeds", "jRoundJets", m_jJetRound_r, jT_noise, m_jJetRound_jet_tower_noise_multiplier, m_jJetRound_jet_min_ET_MeV, m_saveSeeds) ); } if(m_makeRoundLargeRJets) { + /* if( JetAlg::m_SeedMap.find("jRoundLargeRSeeds") == JetAlg::m_SeedMap.end() ) CHECK( JetAlg::SeedGrid(jTs, "jRoundLargeRSeeds", m_dumpSeedsEtaPhi) ); - - ATH_MSG_DEBUG("JFexAlg: SeedFinding with jJetLargeRSeeds; m_jJet_LargeR_seed_size = " << m_jJetRound_LargeR_seed_size << ", m_jJetRound_LargeR_max_r = " << m_jJetRound_LargeR_max_r); - CHECK( JetAlg::SeedFinding(jTs, "jRoundLargeRSeeds", m_jJetRound_LargeR_seed_size, m_jJetRound_LargeR_max_r, jT_noise, m_jJetRound_seed_tower_noise_multiplier, m_jJetRound_seed_total_noise_multiplier, m_jJetRound_LargeR_seed_min_ET_MeV) ); + */ + //ATH_MSG_DEBUG("JFexAlg: SeedFinding with jJetLargeRSeeds; m_jJet_LargeR_seed_size = " << m_jJetRound_LargeR_seed_size << ", m_jJetRound_LargeR_max_r = " << m_jJetRound_LargeR_max_r); + //changing seed from jRoundLargeRSeeds to jRoundSeeds + //CHECK( JetAlg::SeedFinding(jTs, "jRoundLargeRSeeds", m_jJetRound_LargeR_seed_size, m_jJetRound_LargeR_max_r, jT_noise, m_jJetRound_seed_tower_noise_multiplier, m_jJetRound_seed_total_noise_multiplier, m_jJetRound_LargeR_seed_min_ET_MeV) ); - ATH_MSG_DEBUG("JFexAlg: BuildRoundLargeRJet"); - CHECK( JetAlg::BuildRoundJet(jTs, "jRoundLargeRSeeds", "jRoundLargeRJets", m_jJetRound_LargeR_r, jT_noise, m_jJetRound_jet_tower_noise_multiplier, m_jJetRound_jet_total_noise_multiplier, m_jJetRound_LargeR_jet_min_ET_MeV, m_saveSeeds) ); + ATH_MSG_INFO("JFexAlg: BuildRoundLargeRJet"); + CHECK( JetAlg::BuildRoundJet(jTs, "jRoundSeeds", "jRoundLargeRJets", m_jJetRound_LargeR_r, jT_noise, m_jJetRound_jet_tower_noise_multiplier,m_jJetRound_LargeR_jet_min_ET_MeV, /*m_saveSeeds*/ false) ); } if(m_makeJetsFromMap) { @@ -506,6 +445,9 @@ StatusCode JGTowerReader::JFexAlg(const xAOD::JGTowerContainer* jTs){ CHECK( BuildJetsFromMap(jTs) ); } + ATH_MSG_DEBUG("JFexAlg: BuildMET"); + CHECK(METAlg::Baseline_MET(jTs, "jNOISECUT", jT_noise, m_useNegTowers)); + ATH_MSG_DEBUG("JFexAlg: Done"); return StatusCode::SUCCESS; } @@ -550,6 +492,7 @@ StatusCode JGTowerReader::GFexAlg(const xAOD::JGTowerContainer* gTs){ //identifiers are configured for now... totalEt=gTs->at(t)->et(); } + if(gt_em->sampling()!=0)continue; //only build 1 layer for the calo towers to avoid double counting; const float tEt = totalEt; const std::vector index(2, 0); @@ -605,12 +548,12 @@ StatusCode JGTowerReader::GFexAlg(const xAOD::JGTowerContainer* gTs){ xAOD::JGTowerAuxContainer* pu_subAux = new xAOD::JGTowerAuxContainer(); xAOD::JGTowerContainer* pu_sub = new xAOD::JGTowerContainer(); pu_sub->setStore(pu_subAux); - - float rhoA = METAlg::Rho_avg_etaRings(fpga_a, 1, false); - float rhoB = METAlg::Rho_avg_etaRings(fpga_b, 2, false); - float rhoC = METAlg::Rho_avg_etaRings(fpga_c, 3, false); - float rho_barrel = METAlg::Rho_avg(gCaloTowers, false); + float rhoA = METAlg::Rho_avg_etaRings(fpga_a, 1, m_useNegTowers); + float rhoB = METAlg::Rho_avg_etaRings(fpga_b, 2, m_useNegTowers); + float rhoC = METAlg::Rho_avg_etaRings(fpga_c, 3, m_useNegTowers); + + float rho_barrel = METAlg::Rho_avg(gCaloTowers, m_useNegTowers); (*acc_rhoA)(*RhoCont) = rhoA; (*acc_rhoB)(*RhoCont) = rhoB; @@ -699,16 +642,19 @@ StatusCode JGTowerReader::GFexAlg(const xAOD::JGTowerContainer* gTs){ } // jet algorithms + + //Prefer not to delete this yet, but it also is a little confusing to execute... since the seeding should not be run for gFEX. + /* if(JetAlg::m_SeedMap.find("gSeeds") == JetAlg::m_SeedMap.end()) CHECK(JetAlg::SeedGrid(gTs,"gSeeds",m_dumpSeedsEtaPhi)); CHECK(JetAlg::SeedFinding(gTs, "gSeeds", m_gJet_seed_size, m_gJet_max_r, gT_noise, m_gJet_seed_tower_noise_multiplier, m_gJet_seed_total_noise_multiplier, m_gJet_seed_min_ET_MeV)); - + */ // CHECK(JetAlg::SeedFinding(gTs,gSeeds,m_gJet_seed_size,m_gJet_max_r,gJet_thr)); // the diameter of seed, and its range to be local maximum // Careful to ensure the range set to be no tower double counted //CHECK(JetAlg::BuildJet(gTs,gSeeds,gL1Jets,m_gJet_r,gJet_thr)); //default gFex jets are cone jets wih radius of 1.0 - - CHECK(JetAlg::BuildFatJet(*gCaloTowers, "gL1Jets", m_gJet_r, gT_noise, m_gJet_jet_tower_noise_multiplier, m_gJet_jet_total_noise_multiplier, m_gJet_jet_min_ET_MeV, rho_barrel)); - if(m_buildgBlockJets) CHECK(JetAlg::BuildgBlocksJets(gBs, "gBlockJets",rho_barrel)); + //*gCaloTowers + CHECK(JetAlg::BuildFatJet(*gCaloTowers, "gL1Jets", m_gJet_r, gT_noise, m_gJet_tower_noise_multiplier,m_gJet_block_min_ET_MeV, m_gJet_jet_min_ET_MeV, rhoA, rhoB, m_useNegTowers)); + if(m_buildgBlockJets) CHECK(JetAlg::BuildgBlocksJets(gBs, "gBlockJets",rhoA,rhoB)); //gFEX MET algorithms std::vector noNoise; if(m_developerMET){ @@ -733,8 +679,8 @@ StatusCode JGTowerReader::GFexAlg(const xAOD::JGTowerContainer* gTs){ else{ CHECK(METAlg::Baseline_MET(pu_sub, "gXERHO", noNoise, m_useNegTowers)); CHECK(METAlg::Baseline_MET(gCaloTowers,"gXENOISECUT",gT_noise, m_useNegTowers)); - CHECK(METAlg::JwoJ_MET(pu_sub, puSub_gBlocks,"gXEJWOJRHO",m_pTcone_cut,false, false, m_useNegTowers)); - CHECK(METAlg::JwoJ_MET(gTs, gBlocks, "gXEJWOJ",m_pTcone_cut,false, true, /*m_useNegTowers*/ true));//by default, m_useNegTowers=false, but setting it = true here for consistency + CHECK(METAlg::JwoJ_MET(pu_sub, puSub_gBlocks,"gXEJWOJRHO",m_pTcone_cut,false, false, false));//do not use negative towers, so if the subtracted ET is < 0 you remove the tower + CHECK(METAlg::JwoJ_MET(gTs, gBlocks, "gXEJWOJ",m_pTcone_cut,/*bool useEtaBins*/ false, /*bool useRho*/ false, /*m_useNegTowers*/ m_useNegTowers)); CHECK(METAlg::Pufit_MET(gCaloTowers,"gXEPUFIT", m_useNegTowers) ); @@ -781,7 +727,7 @@ StatusCode JGTowerReader::ProcessObjects(){ } CHECK(evtStore()->record(JetCont,it->first.Data())); CHECK(evtStore()->record(JetContAux,Form("%sAux.",it->first.Data()))); - ATH_MSG_DEBUG("Recording JetRoIContainer with name " << it->first.Data() << " and size " << JetCont->size()); + ATH_MSG_INFO("Recording JetRoIContainer with name " << it->first.Data() << " and size " << JetCont->size()); } // Output Seeds @@ -815,7 +761,6 @@ StatusCode JGTowerReader::ProcessObjects(){ CHECK(evtStore()->record(METCont,Form("%s_MET",it->first.Data()))); CHECK(evtStore()->record(METContAux,Form("%s_METAux.",it->first.Data()))); ATH_MSG_DEBUG("Recording EnergySumRoI with name " << it->first.Data() << "_MET"); - ATH_MSG_DEBUG("Recording EnergySumRoI with name " << it->first.Data() << "_MET"); } diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/JetAlg.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/JetAlg.cxx index 9b08d34e9803573550c5d38cb473f11e543e6132..f3d33ec24a4a72b4d2f0145894286db14695fc65 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/JetAlg.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/JetAlg.cxx @@ -18,7 +18,7 @@ StatusCode JetAlg::SeedGrid(const xAOD::JGTowerContainer*towers, TString seedNam std::vector seed_candi_eta; unsigned t_size = towers->size(); - + //find t_maxi=(max eta of all towers) float t_maxi=-999; for(unsigned i=0; iat(i); std::vector SC_indices = tower->SCIndex(); - bool isTile = (fabs(tower->eta())<1.5 && tower->sampling()==1); + bool isTile = (fabs(tower->eta())<1.5 && tower->sampling()==1); if(SC_indices.size()==0 && !isTile) continue; // only want EM towers as centres of barrel seeds (overlap in position with hadronic ones). 0 = barrel EM. 1 = barrel had @@ -102,7 +102,7 @@ StatusCode JetAlg::SeedGrid(const xAOD::JGTowerContainer*towers, TString seedNam } if(m_dumpSeedsEtaPhi) { - ATH_REPORT_MESSAGE_WITH_CONTEXT(MSG::INFO,"JetAlg") << "Dumping seed eta phi"; + ATH_REPORT_MESSAGE_WITH_CONTEXT(MSG::DEBUG,"JetAlg") << "Dumping seed eta phi"; std::cout << "seed eta phi" << std::endl; std::cout << "i_eta" << "\t" << "i_phi" << "\t" << "eta" << "\t" << "phi" << std::endl; for(unsigned i=0; ieta.size(); i++){ @@ -330,8 +330,8 @@ StatusCode JetAlg::SeedFinding(const xAOD::JGTowerContainer*towers, TString seed return StatusCode::SUCCESS; } -StatusCode JetAlg::BuildFatJet(const xAOD::JGTowerContainer towers, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float /*jet_total_noise_multiplier*/, float jet_min_ET_MeV, float rho){ - +StatusCode JetAlg::BuildFatJet(const xAOD::JGTowerContainer towers, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float seed_block_threshold, float jet_min_ET_MeV, float rhoA, float rhoB, bool useNegTowers){ + //std::cout << jetname << " " << jet_r << " noise " << jet_tower_noise_multiplier << " seed " << seed_block_threshold << " jet min " << jet_min_ET_MeV << " rhoA,B" << rhoA <<"," << rhoB << std::endl; std::vector blocks; std::vector> js; @@ -340,6 +340,8 @@ StatusCode JetAlg::BuildFatJet(const xAOD::JGTowerContainer towers, TString jetn TowerObject::TowerGrid grid = TowerObject::TowerGrid(towers); for(const xAOD::JGTower* seed: towers){ + if(fabs(seed->eta())>2.4) continue; + if(!useNegTowers && seed->et() <0) continue; int seedIndex = std::find(towers.begin(), towers.end(), seed) - towers.begin(); std::vector neighbors = grid.neighbors(*seed, 3, 3); @@ -353,7 +355,6 @@ StatusCode JetAlg::BuildFatJet(const xAOD::JGTowerContainer towers, TString jetn const xAOD::JGTower* neighbor = towers.at(neighborIndex); block_area += neighbor->deta()*neighbor->dphi(); neighbor_pt = neighbor->et(); - block_pt += neighbor_pt; } @@ -368,27 +369,49 @@ StatusCode JetAlg::BuildFatJet(const xAOD::JGTowerContainer towers, TString jetn } std::sort(blocks.rbegin(), blocks.rend()); - float pt_cone_cut = 25*Gaudi::Units::GeV; - for(unsigned b = 0; b < blocks.size(); b++){ // const xAOD::JGTower* seed = towers.at(blocks[b].seedIndex()); float block_phi = blocks[b].Phi(); float block_eta = blocks[b].Eta(); - float pt_cone = blocks[b].Pt(); + float pt_block = blocks[b].Pt(); + + float rho = -999; + + if(block_eta<0)rho = rhoA; + if(block_eta>=0)rho = rhoB; float j_Et = 0; - if(pt_cone > pt_cone_cut){ - for(unsigned int t = 0; t < towers.size(); t++){ - const xAOD::JGTower* tower = towers.at(t); + bool print = false; + if(pt_block <= seed_block_threshold) continue; + //if(pt_block > 10e3)print = true; + if(print)std::cout << "pt block " << pt_block << " eta,phi " << block_eta << "," << block_phi << std::endl; + int NtowersInJet = 0; + for(unsigned int t = 0; t < towers.size(); t++){ + const xAOD::JGTower* tower = towers.at(t); + //std::cout << "tower ET: " << tower->et() << " noise " << noise.at(t) << std::endl; + if(!useNegTowers && tower->et()<0)continue; + if(jet_tower_noise_multiplier>0){ + //only apply noise cuts if the noise multiplier is non-zero if(fabs(tower->et()) < noise.at(t) * jet_tower_noise_multiplier) continue; - if(!withinRadius(block_eta, tower->eta(), block_phi, tower->phi(), jet_r, false) ) continue; - j_Et += tower->et(); - }//looping over all towers - j_Et -= 64*rho; - if(j_Et < jet_min_ET_MeV) continue; - std::shared_ptr j = std::make_shared(block_eta, block_phi, j_Et); - js.push_back(j); - }//cone above threshold + } + //std::cout << "tower ET: " << tower->et() << " noise " << noise.at(t) << std::endl; + //withinRadius(float eta1, float eta2, float phi1, float phi2, float dR, bool acceptEqual=false); + if(!withinRadius(block_eta, tower->eta(), block_phi, tower->phi(), jet_r, /*acceptEqual*/ false) ) continue; + float dphi = deltaPhi(block_phi,tower->phi() ); + float deta= block_eta-tower->eta(); + float dR2 = dphi*dphi+deta*deta; + if(print) std::cout << "tower ET: " << tower->et() << " " << tower->eta() << " , " << tower->phi() << std::endl; + if(print)std::cout << "dphi: " << dphi << " deta " << deta << " dR^2 " << dR2 << " dR " << TMath::Sqrt(dR2) << std::endl; + j_Et += tower->et(); + NtowersInJet++; + }//looping over all towers + if(print)std::cout << "N towers: " << NtowersInJet << " jet pT " << j_Et << " 69*Rho: " << 69*rho << " "; + if(rho>0)j_Et -= 69*rho; + if(print)std::cout << "jet pT - rho " << j_Et << std::endl; + if(j_Et < jet_min_ET_MeV) continue; + std::shared_ptr j = std::make_shared(block_eta, block_phi, j_Et); + js.push_back(j); + } m_JetMap[jetname] = js; @@ -460,8 +483,8 @@ StatusCode JetAlg::BuildJet(const xAOD::JGTowerContainer*towers, TString seednam } -StatusCode JetAlg::BuildRoundJet(const xAOD::JGTowerContainer*towers, TString seedname, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float jet_total_noise_multiplier, float jet_min_ET_MeV, bool m_saveSeeds){ - +StatusCode JetAlg::BuildRoundJet(const xAOD::JGTowerContainer*towers, TString seedname, TString jetname, float jet_r, std::vector noise, float jet_tower_noise_multiplier, float jet_min_ET_MeV, bool m_saveSeeds){ + std::vector> js; std::shared_ptr seeds = m_SeedMap[seedname]; if(m_JetMap.find(jetname)!=m_JetMap.end()) js = m_JetMap[jetname]; @@ -490,7 +513,6 @@ StatusCode JetAlg::BuildRoundJet(const xAOD::JGTowerContainer*towers, TString se if(!seeds->local_max.at(eta_ind).at(phi_ind)) continue; float j_et = 0; - float j_totalnoise = 0; ATH_REPORT_MESSAGE_WITH_CONTEXT(MSG::DEBUG,"JetAlg::BuildRoundJet") << "BuildRoundJet: found a local max seed at (eta,phi)=("< 0) { @@ -506,12 +528,10 @@ StatusCode JetAlg::BuildRoundJet(const xAOD::JGTowerContainer*towers, TString se if(applyNoise && fabs(tower->et()) < noise.at(t)*jet_tower_noise_multiplier) continue; if(!withinRadius(eta, tower->eta(), phi, tower->phi(), jet_r,1)) continue; j_et += tower->et(); - j_totalnoise += noise.at(t); ATH_REPORT_MESSAGE_WITH_CONTEXT(MSG::DEBUG,"JetAlg::BuildRoundJet") << " adding tower at (eta,phi)=("<eta()<<","<phi()<<")"; } ATH_REPORT_MESSAGE_WITH_CONTEXT(MSG::DEBUG,"JetAlg::BuildRoundJet") << " final jet has et = " << j_et; if(j_et j = std::make_shared(eta, phi, j_et); @@ -531,7 +551,7 @@ StatusCode JetAlg::BuildRoundJet(const xAOD::JGTowerContainer*towers, TString se return StatusCode::SUCCESS; } -StatusCode JetAlg::BuildgBlocksJets(const xAOD::JGTowerContainer* gBs, TString jetname, float rho){ +StatusCode JetAlg::BuildgBlocksJets(const xAOD::JGTowerContainer* gBs, TString jetname, float rhoA, float rhoB){ std::vector> gbJets; if(m_JetMap.find(jetname)!=m_JetMap.end()) gbJets = m_JetMap[jetname]; @@ -545,6 +565,7 @@ StatusCode JetAlg::BuildgBlocksJets(const xAOD::JGTowerContainer* gBs, TString j const xAOD::JGTower* block = gBs->at(b); const float eta = block->eta(); + if(fabs(eta)>2.4)continue; const float pt = block->et(); // fpga_a @@ -611,17 +632,24 @@ StatusCode JetAlg::BuildgBlocksJets(const xAOD::JGTowerContainer* gBs, TString j } } } - - rho *= 9; //assuming there are ~9 towers per jet + //rho *= 9 bc. assuming there are ~9 towers per jet for(unsigned int i=0; iat(i_max)->eta(); float phi_max = gBs->at(i_max)->phi(); + float rho = -999; + if(eta_max > 0) rho = rhoB; + else if (eta_max <0) rho = rhoA; + rho *=9; float et_max = gBs->at(i_max)->et() - rho; float eta_sec = gBs->at(i_sec)->eta(); + if(eta_sec >0) rho = rhoB; + else if(eta_sec<=0)rho = rhoA; + rho*=9; + float phi_sec = gBs->at(i_sec)->phi(); float et_sec = gBs->at(i_sec)->et() - rho; diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx index 9193e577f87c4a7729c5255aa0342e04c30c5d8d..da5c85223e802335748c4e226dd6d2f11ce36d19 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/METAlg.cxx @@ -17,146 +17,6 @@ #include "TrigT1CaloFexSim/Pufit.h" std::map> METAlg::m_METMap; -std::vector METAlg::check_in_bin (const float &eta, const float &phi, const std::vector < std::pair < float, float> > &eta_bins, const std::vector < std::pair < float, float> > &phi_bins, const float &phi_offset) { - int calculation_bin=-1; - std::vector < int > result; - for (const auto &eta_bin : eta_bins) { - for (const auto &phi_bin : phi_bins) { - calculation_bin++; - if (eta< eta_bin.first || eta > eta_bin.second) continue;; - float phi_t=(phi/TMath::Pi())+1; - if (phi_t phi_bin.second+phi_offset) continue; - result.push_back(calculation_bin); - } - } - return result; -} - -StatusCode METAlg::build_jFEX_bins( std::vector < std::vector < int > > &bins, std::vector < std::vector < int > > &bins_core, const xAOD::JGTowerContainer* towers ) { - std::vector *HAD1towers = new std::vector ; - std::vector *HAD2towers = new std::vector ; - std::vector *HAD3towers = new std::vector ; - std::vector *HADtowers = new std::vector ; - std::vector *EMtowers = new std::vector ; - std::vector *FCAL1 = new std::vector ; - std::vector *FCAL2 = new std::vector ; - std::vector *FCAL3 = new std::vector ; - std::vector *fullFCAL = new std::vector ; - std::vector < std::pair < float , float > > eta_bins = { - {-1.6 , 0.8} , - {-0.8 , 1.6} , - {-2.4 , 0.0} , - {0.0 , 2.4} , - {-4.9 , -0.8} , - {0.8 , 4.9} - }; - std::vector < std::pair < float , float > > eta_bins_core = { - {-0.8 , 0.0} , - {0.0 , 0.8} , - {-1.6 , -0.8} , - {0.8 , 1.6} , - {-4.9 , -1.6} , - {1.6 , 4.9} - }; - - // mutiple of pi , only positive - std::vector < std::pair < float , float > > phi_bins = { {0.0 , 1.0} , {1.0 , 2.0} }; - float phi_offset = +0.25; - for (unsigned int i=0; i <7744; i++){ - if ( i> 0 && i< 1696 ) EMtowers->push_back(i); - if ( i>=3392 && i< 5088 ) EMtowers->push_back(i); - if ( i>=6784 && i< 6816 ) EMtowers->push_back(i); - if ( i>=6848 && i< 6880 ) EMtowers->push_back(i); - if ( i>=6912 && i< 6976 ) EMtowers->push_back(i); - if ( i>=1696 && i< 3392 ) HADtowers->push_back(i); - if ( i>=5088 && i< 6784 ) HADtowers->push_back(i); - if ( i>=6816 && i< 6848 ) HADtowers->push_back(i); - if ( i>=6880 && i< 6912 ) HADtowers->push_back(i); - if ( i>=6976 && i< 7168 ) FCAL1->push_back(i); - if ( i>=7168 && i< 7296 ) FCAL2->push_back(i); - if ( i>=7296 && i< 7360 ) FCAL3->push_back(i); - if ( i>=7360 && i< 7552 ) FCAL1->push_back(i); - if ( i>=7552 && i< 7680 ) FCAL2->push_back(i); - if ( i>=7680 && i< 7744 ) FCAL3->push_back(i); - } - HAD1towers->clear(); - HAD2towers->clear(); - HAD3towers->clear(); - for (const auto &i : *HADtowers) { - const xAOD::JGTower* tower = towers->at(i); - float aeta = std::abs(tower->eta()); - if (aeta< 1.5) HAD1towers->push_back(i); - else if (aeta< 1.6) HAD2towers->push_back(i); - else HAD3towers->push_back(i); - } - - fullFCAL->reserve(fullFCAL->size() + distance(FCAL1->begin(), FCAL1->end())); - fullFCAL->insert(fullFCAL->end(), FCAL1->begin(), FCAL1->end()); - fullFCAL->reserve(fullFCAL->size() + distance(FCAL2->begin(), FCAL2->end())); - fullFCAL->insert(fullFCAL->end(), FCAL2->begin(), FCAL2->end()); - fullFCAL->reserve(fullFCAL->size() + distance(FCAL3->begin(), FCAL3->end())); - fullFCAL->insert(fullFCAL->end(), FCAL3->begin(), FCAL3->end()); - - std::vector < float > rho ( bins.size () , 0.0 ) ; - std::vector < std::vector * > regions; - regions = { HAD1towers, HAD2towers, HAD3towers, EMtowers, fullFCAL }; - - int totalnum=0; - for (const auto ®ion: regions) { - totalnum+=region->size(); - } - - int offset_steps=eta_bins.size()*phi_bins.size(); - int offset_core_steps=eta_bins_core.size()*phi_bins.size(); - - bins.clear(); - bins_core.clear(); - int offset=0; - int offset_core=0; - for (const auto ®ion: regions) { - bins.resize(offset_steps+bins.size()); - bins_core.resize(offset_core_steps+bins_core.size()); - for(const auto &i : *region){ - const xAOD::JGTower* tower = towers->at(i); - float eta = tower->eta(); - float phi = tower->phi(); - std::vector bin= check_in_bin(eta,phi,eta_bins,phi_bins, phi_offset); - std::vector bin_core= check_in_bin(eta,phi,eta_bins_core,phi_bins, phi_offset); - for (const auto j: bin) bins.at(offset+j).push_back(i); - for (const auto j: bin_core) bins_core.at(offset+j).push_back(i); - } - offset+=offset_steps; - offset_core+=offset_core_steps; - } - - int num_bins_core=0; - int total_bins_core=0; - for (const auto &bin: bins_core) { - if (bin.size() !=0 ) num_bins_core++; - total_bins_core+=bin.size(); - } - - int num_bins=0; - int total_bins=0; - for (const auto &bin: bins) { - if (bin.size() !=0 ) num_bins++; - total_bins+=bin.size(); - } - - for (size_t bin=0; bin< bins_core.size(); bin++) { - if (bins_core.at(bin).size() == 0) bins.at(bin).clear(); - } - - int num_bins_c=0; - int total_bins_c=0; - for (const auto &bin: bins) { - if (bin.size() !=0 ) num_bins_c++; - total_bins_c+=bin.size(); - } - return StatusCode::SUCCESS; -} - //------------------------------------------------------------------------------------------------ StatusCode METAlg::Baseline_MET(const xAOD::JGTowerContainer*towers, TString metName, std::vector noise, bool useNegTowers){ @@ -187,92 +47,6 @@ StatusCode METAlg::Baseline_MET(const xAOD::JGTowerContainer*towers, TString met return StatusCode::SUCCESS; } //------------------------------------------------------------------------------------------------ -StatusCode METAlg::jXERHO(const xAOD::JGTowerContainer* towers, TString metName, const std::vector jTowerArea, const std::vector < std::vector < int > > jFEX_bins, const std::vector < std::vector < int > > jFEX_bins_core, float fixed_noise_cut, float rho_up_threshold, float min_noise_cut, xAOD::JGTowerContainer* towers_PUsub ){ - std::vector < int > count ( jFEX_bins.size ( ) , -1 ) ; - std::vector < float > jT_Sum ( jFEX_bins.size ( ) , 0.0 ) ; - std::vector < float > rho ( jFEX_bins.size ( ) , 0.0 ) ; - - // rho calculation - for (unsigned int calculation_bin=0; calculation_binat(i); - float Et = tower->et(); - float Et_area = Et / jTowerArea[i]; - - if (count.at(calculation_bin) == -1) count[calculation_bin]=0; - - if ((fixed_noise_cut rho_tower_et; - rho_tower_et.resize(7744,0.0); - for (unsigned int calculation_bin=0; calculation_binat(i); - float phi=tower->phi(); - float Et = tower->et(); - float Et_sub=0.; - float Et_diff=Et; - Et_diff=(Et-rho.at(calculation_bin)*jTowerArea[i]); - if (Et_diff/jTowerArea[i]>min_noise_cut) { - Et_sub=Et_diff; - } - MET_sub_x -= Et_sub*cos(phi); - MET_sub_y -= Et_sub*sin(phi); - rho_tower_et[i]=Et_sub; - } - } - - towers_PUsub->clear(); - - for ( size_t i =0 ; i< 7744; i++) { - xAOD::JGTower* new_tower = new xAOD::JGTower(*(towers->at(i))); - towers_PUsub->push_back(new_tower); - const xAOD::JGTower* tower = towers->at(i); - float eta=tower->eta(); - float phi=tower->phi(); - float deta=tower->deta(); - float dphi=tower->dphi(); - const int t_ = i; - const int sampling = tower->sampling(); - new_tower->initialize(t_, eta, phi); - new_tower->setdEta(deta); - new_tower->setdPhi(dphi); - new_tower->setEt(rho_tower_et[i]); - const std::vector SCindex= tower->SCIndex(); - const std::vector TileIndex= tower->TileIndex(); - new_tower->setSCIndex(SCindex); - new_tower->setTileIndex(TileIndex); - new_tower->setSampling(sampling); - } - - float et_MET_sub = sqrt(MET_sub_x * MET_sub_x + MET_sub_y * MET_sub_y); - float phi_MET_sub=TMath::ACos(MET_sub_x/et_MET_sub); - if (MET_sub_y<0) { - phi_MET_sub = -phi_MET_sub; - } - - std::shared_ptr met = std::make_shared(); - met->ex=MET_sub_x; - met->ey=MET_sub_y; - met->phi = phi_MET_sub; - met->et = et_MET_sub;; - - - if(m_METMap.find(metName)==m_METMap.end()) m_METMap[metName] = met; - - return StatusCode::SUCCESS; -} - StatusCode METAlg::SubtractRho_MET(const xAOD::JGTowerContainer* towers, TString metName, bool useEtaBins, bool useRMS, bool useMedian, bool useNegTowers){ float EtMiss = 0; @@ -619,6 +393,7 @@ float METAlg::Rho_avg_etaRings(const xAOD::JGTowerContainer* towers, int fpga, b float Et = tower->et(); if(!useNegTowers && Et < 0) continue; + if(Et > et_max)continue; if(eta < 0 && eta >=-2.5){ if(eta > -2.4) area_a += 1.; @@ -632,8 +407,8 @@ float METAlg::Rho_avg_etaRings(const xAOD::JGTowerContainer* towers, int fpga, b if(fabs(eta) < 3.2) area_c += 1.; else area_c += 4.; } - - if(Et < et_max) rho+=Et; //sum ET together to compute total + rho+=Et; + //if(Et < et_max) rho+=Et; //sum ET together to compute total } //after computing total, divide by total area diff --git a/Trigger/TrigT1/TrigT1CaloFexSim/src/components/TrigT1CaloFexSim_entries.cxx b/Trigger/TrigT1/TrigT1CaloFexSim/src/components/TrigT1CaloFexSim_entries.cxx index 5a5d6464a1bd821ff50a86167d4b32ecdfae4ac6..1d828b932c2b39981f7334fb1cbbdc5bce443943 100644 --- a/Trigger/TrigT1/TrigT1CaloFexSim/src/components/TrigT1CaloFexSim_entries.cxx +++ b/Trigger/TrigT1/TrigT1CaloFexSim/src/components/TrigT1CaloFexSim_entries.cxx @@ -9,8 +9,13 @@ DECLARE_ALGORITHM_FACTORY( JGTowerReader ) #include "TrigT1CaloFexSim/JGTowerMaker.h" DECLARE_ALGORITHM_FACTORY( JGTowerMaker ) + +#include "../JFexEleTau.h" +DECLARE_ALGORITHM_FACTORY( JFexEleTau ) + DECLARE_FACTORY_ENTRIES( TrigT1CaloFexSim ) { + DECLARE_ALGORITHM( JFexEleTau ); DECLARE_ALGORITHM( JGTowerMaker ); DECLARE_ALGORITHM( JGTowerReader ); }