diff --git a/Trigger/TrigT1/TrigT1NSW/share/NSWL1.py b/Trigger/TrigT1/TrigT1NSW/share/NSWL1.py index 43a234d75a563c9c221ab9da66ba78d769cfbb51..753b5e96e3d22210db30161c31a1e3762d066468 100644 --- a/Trigger/TrigT1/TrigT1NSW/share/NSWL1.py +++ b/Trigger/TrigT1/TrigT1NSW/share/NSWL1.py @@ -1,53 +1,64 @@ -# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration +import glob +import os -####################################### -# some basic settings -####################################### +################################################### +# Some basic settings and input file submission +################################################### MessageSvc.defaultLimit=100 MessageSvc.useColors = True MessageSvc.Format = "% F%30W%S%7W%R%T %0W%M" -from AthenaCommon.AthenaCommonFlags import athenaCommonFlags +from AthenaCommon.AthenaCommonFlags import athenaCommonFlags, jobproperties +athenaCommonFlags.AllowIgnoreConfigError=False #This job will stop if an include fails athenaCommonFlags.EvtMax = -1 athenaCommonFlags.SkipEvents = 0 from AthenaCommon.AppMgr import ServiceMgr as svcMgr import AthenaPoolCnvSvc.ReadAthenaPool # needed to use EventSelector -svcMgr.EventSelector.InputCollections=["input.rdo.pool.root"] -####################################### - - -####################################### -# initialize the geometry -####################################### -from AthenaCommon.GlobalFlags import jobproperties -jobproperties.Global.DetDescrVersion="ATLAS-R3S-2021-01-00-00" +#to run on mutliple files at once please use -c "customInput='/some/path/*pattern*.root'" +if 'customInput' not in locals() or 'customInput' not in globals(): + print("customInput not defined: the job will be launched on the GRID") + svcMgr.EventSelector.InputCollections=jobproperties.AthenaCommonFlags.FilesInput.get_Value() +else: + if(os.path.isdir(customInput)): customInput+="/*.root" + svcMgr.EventSelector.InputCollections=glob.glob(customInput) + +################################################### +# Detector flags: needed to setup geometry tools +################################################### +from AthenaCommon.DetFlags import DetFlags +DetFlags.detdescr.Muon_setOn() +DetFlags.sTGC_setOff() +DetFlags.Micromegas_setOn() +DetFlags.digitize.Micromegas_setOn() +DetFlags.digitize.sTGC_setOff() +DetFlags.Truth_setOn() +DetFlags.Print() + +################################################### +# Initialize ATLAS geometry +################################################### +jobproperties.Global.DetDescrVersion="ATLAS-R3S-2021-01-00-02" from AtlasGeoModel import SetGeometryVersion from AtlasGeoModel import GeoModelInit -from GeoModelSvc.GeoModelSvcConf import GeoModelSvc -GeoModelSvc = GeoModelSvc() +from AtlasGeoModel import MuonGM +from AtlasGeoModel import Agdd2Geo + from AtlasGeoModel.MuonGMJobProperties import MuonGeometryFlags -# initialize the MuonIdHelperService from MuonIdHelpers.MuonIdHelpersConf import Muon__MuonIdHelperSvc svcMgr += Muon__MuonIdHelperSvc("MuonIdHelperSvc",HasCSC=MuonGeometryFlags.hasCSC(), HasSTgc=MuonGeometryFlags.hasSTGC(), HasMM=MuonGeometryFlags.hasMM()) -# create the MuonDetectorTool (which creates the MuonDetectorManager needed by PadTdsOfflineTool) -from MuonGeoModel.MuonGeoModelConf import MuonDetectorTool -GeoModelSvc.DetectorTools += [ MuonDetectorTool(HasCSC=MuonGeometryFlags.hasCSC(), HasSTgc=MuonGeometryFlags.hasSTGC(), HasMM=MuonGeometryFlags.hasMM()) ] -####################################### - - -####################################### -# now the trigger related things -####################################### +################################################### +# Initialize the trigger related things +################################################### include('TrigT1NSW/TrigT1NSW_jobOptions.py') #Switch on and off trigger simulaton components sTGC / MicroMegas -#October 2019 : MM not working so keep it False until fixed - -topSequence.NSWL1Simulation.DosTGC=True +topSequence.NSWL1Simulation.DosTGC=False topSequence.NSWL1Simulation.UseLookup=False #use lookup table for the pad trigger -topSequence.NSWL1Simulation.DoMM=False +topSequence.NSWL1Simulation.DoMM=True +topSequence.NSWL1Simulation.DoMMDiamonds=True topSequence.NSWL1Simulation.NSWTrigRDOContainerName="NSWTRGRDO" topSequence.NSWL1Simulation.StripSegmentTool.rIndexScheme=0 @@ -76,12 +87,12 @@ topSequence.NSWL1Simulation.StripSegmentTool.OutputLevel=INFO #S.I 2019 : Below code handles enabling/disabling of ntuples for the Tools according to the master flag (NSWL1Simulation) #in principle we wouldnt need the tuning here but ntuple making and the trigger code is quite tangled so this is just a workaround for now -#ntuple code must be totally stripped off from trigger Tools. One way of doing is to create a separate tool and implement methods that takes -# std::vector<std::shared_ptr<PadData>> pads; +#ntuple code must be totally stripped off from trigger Tools. One way of doing is to create a separate tool and implement methods that takes +# std::vector<std::shared_ptr<PadData>> pads; # std::vector<std::unique_ptr<PadTrigger>> padTriggers; # std::vector<std::unique_ptr<StripData>> strips; # std::vector< std::unique_ptr<StripClusterData> > clusters; -# as arguments (see NSWL1Simulation.cxx) +# as arguments (see NSWL1Simulation.cxx) if topSequence.NSWL1Simulation.DoNtuple: @@ -99,7 +110,7 @@ if topSequence.NSWL1Simulation.DoNtuple: ServiceMgr += NTupleSvc() ServiceMgr.THistSvc.Output += [ "NSWL1Simulation DATAFILE='NSWL1Simulation.root' OPT='RECREATE'" ] - print ServiceMgr + print(ServiceMgr) else:#to avoid any possible crash. If DoNtuple is set to true for a tool but false for NSWL1Simulation the code will crash # ideally that should have been already handled in C++ side diff --git a/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.cxx b/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.cxx index 5db7bcfa5cdff3a3757e6638963a2c52057bea97..99e2e1731aa8bb885c5b2912aadd8f05197f92a9 100644 --- a/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.cxx +++ b/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.cxx @@ -25,8 +25,8 @@ namespace NSWL1 { m_pad_trigger("NSWL1::PadTriggerLogicOfflineTool",this), m_pad_trigger_lookup("NSWL1::PadTriggerLookupTool",this), m_strip_tds("NSWL1::StripTdsOfflineTool",this), - m_strip_cluster("NSWL1::StripClusterTool",this), - m_strip_segment("NSWL1::StripSegmentTool",this), + //m_strip_cluster("NSWL1::StripClusterTool",this), + //m_strip_segment("NSWL1::StripSegmentTool",this), TODO: this line makes the code crash in initialization... please, sTGC friends, fix it!!! m_mmstrip_tds("NSWL1::MMStripTdsOfflineTool",this), m_mmtrigger("NSWL1::MMTriggerTool",this), m_tree(nullptr), @@ -36,24 +36,24 @@ namespace NSWL1 { // Property setting general behaviour: declareProperty( "DoOffline", m_doOffline = false, "Steers the offline emulation of the LVL1 logic" ); - declareProperty( "UseLookup", m_useLookup = false, "Toggle Lookup mode on and off default is the otf(old) mode" ); - declareProperty( "DoNtuple", m_doNtuple = false, "Create an ntuple for data analysis" ); - declareProperty( "DoMM", m_doMM = false, "Run data analysis for MM" ); - declareProperty( "DosTGC", m_dosTGC = true, "Run data analysis for sTGCs" ); + declareProperty( "UseLookup", m_useLookup = false, "Toggle Lookup mode on and off default is the otf(old) mode" ); + declareProperty( "DoNtuple", m_doNtuple = true, "Create an ntuple for data analysis" ); + declareProperty( "DoMM", m_doMM = true, "Run data analysis for MM" ); + declareProperty( "DoMMDiamonds", m_doMMDiamonds = false, "Run data analysis for MM using Diamond Roads algorithm" ); + declareProperty( "DosTGC", m_dosTGC = false, "Run data analysis for sTGCs" ); // declare monitoring tools - declareProperty( "AthenaMonTools", m_monitors, "List of monitoring tools to be run with this instance, if incorrect then tool is silently skipped."); - - declareProperty( "PadTdsTool", m_pad_tds, "Tool that simulates the functionalities of the PAD TDS"); + declareProperty( "AthenaMonTools", m_monitors, "List of monitoring tools to be run with this instance, if incorrect then tool is silently skipped."); + declareProperty( "PadTdsTool", m_pad_tds, "Tool that simulates the functionalities of the PAD TDS"); //PadTriggerTool : in principle can be totally wuiped out. necesary for ntuples currently. Once you isolate ntuple making code and the trigger code you can abandon this method. Things ae still tangled a bit somewhow so keep it just in case - declareProperty( "PadTriggerTool", m_pad_trigger, "Tool that simulates the pad trigger logic"); - declareProperty( "PadTriggerLookupTool", m_pad_trigger_lookup, "Tool that is used to lookup pad trigger patterns per execute against the same LUT as in trigger FPGA"); - declareProperty( "StripTdsTool", m_strip_tds, "Tool that simulates the functionalities of the Strip TDS"); - declareProperty( "StripClusterTool",m_strip_cluster, "Tool that simulates the Strip Clustering"); - declareProperty( "StripSegmentTool",m_strip_segment, "Tool that simulates the Segment finding"); - declareProperty( "MMStripTdsTool", m_mmstrip_tds, "Tool that simulates the functionalities of the MM STRIP TDS"); - declareProperty( "MMTriggerTool", m_mmtrigger, "Tool that simulates the MM Trigger"); - declareProperty("NSWTrigRDOContainerName", m_trigRdoContainer = "NSWTRGRDO"," Give a name to NSW trigger rdo container"); + declareProperty( "PadTriggerTool", m_pad_trigger, "Tool that simulates the pad trigger logic"); + declareProperty( "PadTriggerLookupTool", m_pad_trigger_lookup, "Tool that is used to lookup pad trigger patterns per execute against the same LUT as in trigger FPGA"); + declareProperty( "StripTdsTool", m_strip_tds, "Tool that simulates the functionalities of the Strip TDS"); + declareProperty( "StripClusterTool", m_strip_cluster, "Tool that simulates the Strip Clustering"); + declareProperty( "StripSegmentTool", m_strip_segment, "Tool that simulates the Segment finding"); + declareProperty( "MMStripTdsTool", m_mmstrip_tds, "Tool that simulates the functionalities of the MM STRIP TDS"); + declareProperty( "MMTriggerTool", m_mmtrigger, "Tool that simulates the MM Trigger"); + declareProperty( "NSWTrigRDOContainerName", m_trigRdoContainer = "NSWTRGRDO"," Give a name to NSW trigger rdo container"); } @@ -80,7 +80,6 @@ namespace NSWL1 { } // retrieving the private tools implementing the simulation - if(m_dosTGC){ ATH_CHECK(m_pad_tds.retrieve()); if(m_useLookup){ @@ -90,13 +89,15 @@ namespace NSWL1 { ATH_CHECK(m_pad_trigger.retrieve()); } ATH_CHECK(m_strip_tds.retrieve()); - ATH_CHECK(m_strip_cluster.retrieve()); - ATH_CHECK(m_strip_segment.retrieve()); + //ATH_CHECK(m_strip_cluster.retrieve()); + //ATH_CHECK(m_strip_segment.retrieve()); } if(m_doMM ){ ATH_CHECK(m_mmtrigger.retrieve()); + if(m_doMMDiamonds) ATH_CHECK( m_mmtrigger->initDiamondAlgorithm() ); } + // Connect to Monitoring Service ATH_CHECK(m_monitors.retrieve()); return StatusCode::SUCCESS; @@ -105,7 +106,6 @@ namespace NSWL1 { StatusCode NSWL1Simulation::start() { ATH_MSG_INFO("start " << name() ); - for ( auto& mon : m_monitors ) { ATH_CHECK(mon->bookHists()); } @@ -137,8 +137,8 @@ namespace NSWL1 { } ATH_CHECK( m_strip_tds->gather_strip_data(strips,padTriggers) ); - ATH_CHECK( m_strip_cluster->cluster_strip_data(strips,clusters) ); - ATH_CHECK( m_strip_segment->find_segments(clusters,trgContainer) ); + //ATH_CHECK( m_strip_cluster->cluster_strip_data(strips,clusters) ); + //ATH_CHECK( m_strip_segment->find_segments(clusters,trgContainer) ); auto rdohandle = SG::makeHandle( m_trigRdoContainer ); ATH_CHECK( rdohandle.record( std::move(trgContainer))); @@ -146,7 +146,7 @@ namespace NSWL1 { //retrive the MM Strip hit data if(m_doMM){ - ATH_CHECK( m_mmtrigger->runTrigger() ); + ATH_CHECK( m_mmtrigger->runTrigger(m_doMMDiamonds) ); } for ( auto& mon : m_monitors) { ATH_CHECK(mon->fillHists()); @@ -162,6 +162,7 @@ namespace NSWL1 { for ( auto& mon : m_monitors ) { ATH_CHECK(mon->finalHists()); } + if(m_doMM) ATH_CHECK( m_mmtrigger->finalizeDiamondAlgorithm(m_doMMDiamonds) ); return StatusCode::SUCCESS; } diff --git a/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.h b/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.h index a5ddeb28f75076b8380c1f4b123ee07aa4078756..119a2dc270beaf43d9fbbb058ca712234da6a7eb 100644 --- a/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.h +++ b/Trigger/TrigT1/TrigT1NSW/src/NSWL1Simulation.h @@ -42,7 +42,7 @@ namespace NSWL1 { * or the detailed hardware simulation. * * @authors Alessandro Di Mattia <dimattia@cern.ch>, Geraldine Conti <geraldine.conti@cern.ch> - * + * Major updates for Release 22 processing: Francesco Giuseppe Gravili <francesco.giuseppe.gravili@cern.ch> * */ @@ -74,7 +74,7 @@ namespace NSWL1 { ToolHandle < IStripClusterTool > m_strip_cluster; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink ToolHandle < IStripSegmentTool > m_strip_segment; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink ToolHandle < IMMStripTdsTool > m_mmstrip_tds; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink - ToolHandle < IMMTriggerTool > m_mmtrigger; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink + ToolHandle < IMMTriggerTool > m_mmtrigger; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink // put analysis variables here TTree* m_tree; //!< analysis ntuple @@ -87,6 +87,7 @@ namespace NSWL1 { bool m_useLookup; bool m_doNtuple; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink bool m_doMM; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink + bool m_doMMDiamonds; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink bool m_dosTGC; //!< property, see @link NSWL1Simulation::NSWL1Simulation @endlink diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/CMakeLists.txt b/Trigger/TrigT1/TrigT1NSWSimTools/CMakeLists.txt index 2b00061e32c04208d9a53ea0c63c85a8f2bfc8e7..f19a4b86d19b8095029f0e69421d5a86e839cef7 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/CMakeLists.txt +++ b/Trigger/TrigT1/TrigT1NSWSimTools/CMakeLists.txt @@ -13,7 +13,7 @@ atlas_add_library( TrigT1NSWSimToolsLib src/*.cxx PUBLIC_HEADERS TrigT1NSWSimTools INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} - LINK_LIBRARIES ${Boost_LIBRARIES} ${ROOT_LIBRARIES} AthenaBaseComps AthenaKernel GaudiKernel GeoPrimitives Identifier MuonDigitContainer MuonIdHelpersLib MuonRDO MuonReadoutGeometry MuonSimData MuonSimEvent RegSelLUT TrkSurfaces + LINK_LIBRARIES ${Boost_LIBRARIES} ${ROOT_LIBRARIES} AthenaBaseComps AthenaKernel FourMomUtils GaudiKernel GeoPrimitives Identifier MuonDigitContainer MuonIdHelpersLib MuonRDO MuonReadoutGeometry MuonSimData MuonSimEvent RegSelLUT TrkSurfaces PRIVATE_LINK_LIBRARIES ${CLHEP_LIBRARIES} AGDDKernel AtlasHepMCLib EventInfo GeneratorObjects MuonAGDDDescription MuonRegionSelectorLib PathResolver StoreGateLib TrackRecordLib ) # Component(s) in the package: diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/IMMTriggerTool.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/IMMTriggerTool.h index e88ebb9618d7fbe9c57a91c71d01fe360781df10..dc45f7af1a01f77de1bbdcf51b36c350c2dd4fca 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/IMMTriggerTool.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/IMMTriggerTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #ifndef IMMTRIGGERTOOL_H @@ -23,8 +23,9 @@ namespace NSWL1 { public: virtual ~IMMTriggerTool() {} - virtual StatusCode runTrigger() = 0; - + virtual StatusCode runTrigger(const bool do_MMDiamonds) = 0; + virtual StatusCode initDiamondAlgorithm() = 0; + virtual StatusCode finalizeDiamondAlgorithm(const bool do_MMDiamonds) = 0; static const InterfaceID& interfaceID() { static const InterfaceID IID_IMMTriggerTool("NSWL1::IMMTriggerTool", 1 ,0); diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMLoadVariables.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMLoadVariables.h index bc0c99e3f1f0d5c497781c45a9ae92cb90507e74..af03f93404497b66f7ea0ad1c36585cf5d07e53d 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMLoadVariables.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMLoadVariables.h @@ -18,7 +18,6 @@ class MmDigit; class StoreGateSvc; class MMT_Parameters; class MmDigitContainer; -class TVector3; namespace MuonGM { class MuonDetectorManager; @@ -28,18 +27,20 @@ namespace MuonGM { public: - MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetectorManager* detManager, const MmIdHelper* idhelper, MMT_Parameters *par); + MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetectorManager* detManager, const MmIdHelper* idhelper, std::map<std::string,MMT_Parameters*> pars); - StatusCode getMMDigitsInfo(std::vector<digitWrapper>& entries, std::map<hitData_key,hitData_entry>& Hits_Data_Set_Time, std::map<int,evInf_entry>& Event_Info); + StatusCode getMMDigitsInfo(std::map<std::pair<int,unsigned int>,std::vector<digitWrapper> >& entries, + std::map<std::pair<int,unsigned int>,std::map<hitData_key,hitData_entry> >& Hits_Data_Set_Time, + std::map<std::pair<int,unsigned int>,evInf_entry>& Event_Info); //Import_Athena..._.m stuff double phi_shift(double athena_phi,const std::string& wedgeType, int stationPhi) const; int Get_VMM_chip(int strip) const; //*** Not Finished... Rough - int strip_number(int station,int plane,int spos)const; + int strip_number(int station, int plane, int spos, const MMT_Parameters *par)const; int Get_Strip_ID(double X,double Y,int plane) const; bool Mimic_VMM_Chip_Deadtime(hitData_entry& candy); - void xxuv_to_uvxx(TVector3& hit,int plane)const; - void hit_rot_stereo_fwd(TVector3& hit)const;//x to v, u to x - void hit_rot_stereo_bck(TVector3& hit)const;//x to u, v to x + void xxuv_to_uvxx(ROOT::Math::XYZVector& hit, int plane, const MMT_Parameters *par)const; + void hit_rot_stereo_fwd(ROOT::Math::XYZVector& hit, const MMT_Parameters *par)const;//x to v, u to x + void hit_rot_stereo_bck(ROOT::Math::XYZVector& hit, const MMT_Parameters *par)const;//x to u, v to x struct histogramVariables{ std::vector<std::string> *m_NSWMM_dig_stationName; @@ -142,7 +143,7 @@ namespace MuonGM { const MuonGM::MuonDetectorManager* m_detManager; //!< MuonDetectorManager const MmIdHelper* m_MmIdHelper; //!< MM offline Id helper StoreGateSvc* m_evtStore; - MMT_Parameters* m_par{}; + std::map<std::string, MMT_Parameters*> m_pars{}; bool m_striphack{}; std::string getWedgeType(const MmDigitContainer *nsw_MmDigitContainer); diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMStripTdsOfflineTool.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMStripTdsOfflineTool.h index edb468731a08bdbe1d130adb3255eb59458915b1..2a79648977dfa381bb1266b6f1cc5e9144c0e545 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMStripTdsOfflineTool.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMStripTdsOfflineTool.h @@ -94,9 +94,9 @@ namespace NSWL1 { int Get_Strip_ID(double X,double Y,int plane) const; //x <---> u/v switches - void xxuv_to_uvxx(TVector3& hit,int plane)const; - void hit_rot_stereo_fwd(TVector3& hit)const;//x to v, u to x - void hit_rot_stereo_bck(TVector3& hit)const;//x to u, v to x + void xxuv_to_uvxx(ROOT::Math::XYZVector& hit,int plane)const; + void hit_rot_stereo_fwd(ROOT::Math::XYZVector& hit)const;//x to v, u to x + void hit_rot_stereo_bck(ROOT::Math::XYZVector& hit)const;//x to u, v to x //MMT_Loader stuff end diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Diamond.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Diamond.h new file mode 100644 index 0000000000000000000000000000000000000000..2c7082833c7d1b75c167a1981750242ca6266193 --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Diamond.h @@ -0,0 +1,101 @@ +/* + * Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + */ + +#ifndef MMT_DIAMOND_H +#define MMT_DIAMOND_H + +#include "AthenaBaseComps/AthMessaging.h" +#include "MMT_Road.h" +#include <vector> + +struct slope_t { + slope_t(int ev=-1, int bc=-1, unsigned int tC=999, unsigned int rC=999, int iX=-1, int iU=-1, int iV=-1, unsigned int uvb=999, unsigned int xb=999, + unsigned int uvm=999, unsigned int xm=999, int age=-1, double mxl=999., double xavg=999., double uavg=999., double vavg=999., + double th=999., double eta=999., double dth=999., double phi=999., double phiS=999.); + unsigned int event; + int BC; + unsigned int totalCount; + unsigned int realCount; + int iRoad; + int iRoadu; + int iRoadv; + unsigned int uvbkg; + unsigned int xbkg; + unsigned int uvmuon; + unsigned int xmuon; + int age; + double mxl; + double xavg; + double uavg; + double vavg; + double theta; + double eta; + double dtheta; + double phi; + double phiShf; +}; + +struct diamond_t { + unsigned int wedgeCounter; + char sector; + int phi; + std::vector<MMT_Road*> ev_roads; + std::vector<slope_t> slopes; + std::vector<MMT_Hit*> ev_hits; +}; + +class MMT_Diamond : public AthMessaging { + public: + MMT_Diamond(const MuonGM::MuonDetectorManager* detManager); + + void clearEvent(); + void createRoads_fillHits(const unsigned int iterator, std::vector<hitData_entry> &hitDatas, const MuonGM::MuonDetectorManager* detManager, MMT_Parameters *par, const int phi); + void findDiamonds(const unsigned int iterator, const double &sm_bc, const int &event); + double phiShift(const int &n, const double &phi, const char &side); + std::vector<diamond_t> getDiamondVector() const { return m_diamonds; } + diamond_t getDiamond(const unsigned int iterator) const { return m_diamonds.at(iterator); } + std::vector<double> getHitSlopes() const { return m_hitslopes; } + std::vector<MMT_Hit*> getHitVector(const unsigned int iterator) const { return m_diamonds.at(iterator).ev_hits; } + std::vector<slope_t> getSlopeVector(const unsigned int iterator) const { return m_diamonds.at(iterator).slopes; } + unsigned int getDiamondSize() const { return m_diamonds.size(); } + int getUVfactor() const { return m_uvfactor; } + void printHits(const unsigned int iterator); + void resetSlopes(); + void setTrapezoidalShape(bool flag) { m_trapflag = flag; } + void setUVfactor(int factor) { m_uvfactor = factor; } + bool isTrapezoidalShape() const { return m_trapflag; } + + int getStationPhi() const { return m_phi; } + int getRoadSize() const { return m_roadSize; } + void setRoadSize(int size) { m_roadSize = size; } + int getRoadSizeUpX() const { return m_roadSizeUpX; } + void setRoadSizeUpX(int sizeUp) { m_roadSizeUpX = sizeUp; } + int getRoadSizeDownX() const { return m_roadSizeDownX; } + void setRoadSizeDownX(int sizeDown) { m_roadSizeDownX = sizeDown; } + int getRoadSizeUpUV() const { return m_roadSizeUpUV; } + void setRoadSizeUpUV(int sizeUpUV) { m_roadSizeUpUV = sizeUpUV; } + int getRoadSizeDownUV() const { return m_roadSizeDownUV; } + void setRoadSizeDownUV(int sizeDownUV) { m_roadSizeDownUV = sizeDownUV; } + char getSector() const { return m_sector; } + unsigned int getXthreshold() const { return m_xthr; } + void setXthreshold(int threshold) { m_xthr = threshold; } + bool getUV() const { return m_uvflag; } + void setUV(bool flag) { m_uvflag = flag; } + unsigned int getUVthreshold() const { return m_uvthr; } + void setUVthreshold(int threshold) { m_uvthr = threshold; } + + private: + const MuonGM::MuonDetectorManager* m_detManager{}; //!< MuonDetectorManager + bool m_trapflag; + int m_uvfactor; + bool m_uvflag; + int m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV; + int m_xthr, m_uvthr; + int m_phi; + char m_sector; + + std::vector<diamond_t> m_diamonds; + std::vector<double> m_hitslopes; +}; +#endif diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Finder.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Finder.h index ead0f1119c097774cfa028b8269d06ae6d53ec09..764936c7934167e21aae74c2e9caeab47cf424e8 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Finder.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Finder.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #ifndef MMT_FINDER_H diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Fitter.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Fitter.h index 3f6fb557f0b3173a835dfc5d6fdcdb5843973f0c..1ba687f7b3dede7335a5ea39c03afceaf8c44128 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Fitter.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Fitter.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #ifndef MMT_FITTER_H @@ -27,13 +27,13 @@ class MMT_Fitter : AthMessaging { //functions translated int Filter_UV(std::vector<Hit>& track) const; - float32fixed<2> Get_Global_Slope(const std::vector<Hit>& track, const std::string& type) const; - float32fixed<2> Get_Local_Slope(const std::vector<Hit>& Track,double theta=-999.,double phi=-999.) const; - ROI Get_ROI(float32fixed<2> M_x,float32fixed<2> M_u,float32fixed<2> M_v,const std::vector<Hit>&track) const; + double Get_Global_Slope(const std::vector<Hit>& track, const std::string& type) const; + double Get_Local_Slope(const std::vector<Hit>& Track,double theta=-999.,double phi=-999.) const; + ROI Get_ROI(double M_x,double M_u,double M_v,const std::vector<Hit>&track) const; double phi_correct_factor(const std::vector<Hit>&track)const; - float32fixed<2> Get_Delta_Theta(float32fixed<2> M_local,float32fixed<2> M_global) const; - float32fixed<2> Get_Delta_Theta_division(float32fixed<2> M_local,float32fixed<2> M_global, float32fixed<4> a=1.) const; - int Rough_ROI_temp(float32fixed<4> theta, float32fixed<4> phi) const; + double Get_Delta_Theta(double M_local,double M_global) const; + double Get_Delta_Theta_division(double M_local,double M_global, double a=1.) const; + int Rough_ROI_temp(double theta, double phi) const; //sim hit code stuff int track_to_index(const std::vector<Hit>&track)const; @@ -45,12 +45,12 @@ class MMT_Fitter : AthMessaging { double ideal_zbar(const std::vector<Hit>& Track)const; //translated from Table_Generators.m - float32fixed<2> LG_lgr(int ilgr, double a, int number_LG_regions, float32fixed<2> min, float32fixed<2> max) const; - float32fixed<2> mult_factor_lgr(int ilgr, double a, int number_LG_regions, float32fixed<2> min, float32fixed<2> max) const; - float32fixed<4> Slope_Components_ROI_val(int jy, int ix, int thetaphi) const; - float32fixed<4> Slope_Components_ROI_theta(int jy, int ix) const; - float32fixed<4> Slope_Components_ROI_phi(int jy, int ix) const; - float32fixed<2> DT_Factors_val(int i, int j) const; + double LG_lgr(int ilgr, double a, int number_LG_regions, double min, double max) const; + double mult_factor_lgr(int ilgr, double a, int number_LG_regions, double min, double max) const; + double Slope_Components_ROI_val(int jy, int ix, int thetaphi) const; + double Slope_Components_ROI_theta(int jy, int ix) const; + double Slope_Components_ROI_phi(int jy, int ix) const; + double DT_Factors_val(int i, int j) const; private: @@ -59,7 +59,7 @@ class MMT_Fitter : AthMessaging { //Fitter components int m_number_LG_regions,m_n_fit; - float32fixed<2> m_LG_min,m_LG_max; + double m_LG_min,m_LG_max; std::vector<int> q_planes(char type) const;//return position of what planes are where }; diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Hit.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Hit.h new file mode 100644 index 0000000000000000000000000000000000000000..86e0057b8c0fadef3860acb4ad58eaa9729057ce --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Hit.h @@ -0,0 +1,79 @@ +/* + * Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + */ + +#ifndef MMT_HIT_H +#define MMT_HIT_H + +#include "AthenaBaseComps/AthMessaging.h" +#include "MMT_struct.h" + +namespace MuonGM { + class MuonDetectorManager; +} + +class MMT_Hit : public AthMessaging { + public: + MMT_Hit(char wedge, hitData_entry entry, const MuonGM::MuonDetectorManager* detManager); + MMT_Hit(const MMT_Hit &hit); + MMT_Hit& operator=(const MMT_Hit&); + + int getART() const { return m_ART_ASIC; } + int getAge() const { return m_age; } + int getBC() const { return m_BC_time; } + int getChannel() const { return m_strip; } + int getGasGap() const { return m_gasgap; } + std::string getModule() const { return m_module; } + int getMultiplet() const { return m_multiplet; } + int getPlane() const { return m_plane; } + char getSector() const { return m_sector; } + double getRZSlope() const { return m_RZslope; } + double getYZSlope() const { return m_YZslope; } + int getVMM() const { return m_VMM_chip; } + int getMMFE8() const { return m_MMFE_VMM; } + std::string getStationName() const { return m_station_name; } + int getStationEta() const { return m_station_eta; } + int getStationPhi() const { return m_station_phi; } + double getR() const { return m_R; } + double getRi() const { return m_Ri; } + double getX() const { return m_localX; } + double getY() const { return m_Y; } + double getZ() const { return m_Z; } + bool isNoise() const { return m_isNoise; } + bool isX() const; + bool isU() const; + bool isV() const; + void printHit(); + void setAge(int age) { m_age = age; } + void setAsNoise() { m_isNoise = true; } + void setBC(int bc) { m_BC_time = bc; } + void setHitProperties(const Hit &hit); + void setRZSlope(double slope) { m_RZslope = slope; } + void setYZSlope(double slope) { m_YZslope = slope; } + void setY(double y) { m_Y = y; } + void setZ(double z) { m_Z = z; } + void updateHitProperties(const MMT_Parameters *par); + bool verifyHit(); + + private: + char m_sector; + std::string m_module, m_station_name; + int m_VMM_chip; + int m_MMFE_VMM; + int m_ART_ASIC; + int m_plane; + int m_station_eta; + int m_station_phi; + int m_multiplet; + int m_gasgap; + int m_strip; + double m_localX; + double m_RZslope, m_YZslope; + int m_BC_time, m_age; + double m_Y, m_Z, m_R, m_Ri; + bool m_isNoise; + + const MuonGM::MuonDetectorManager* m_detManager{}; //!< MuonDetectorManager + const MuonGM::MuonDetectorManager* getDetManager() { return m_detManager; } +}; +#endif diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Road.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Road.h new file mode 100644 index 0000000000000000000000000000000000000000..635b9f778414888c42aed64e6c635139f4e10277 --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_Road.h @@ -0,0 +1,94 @@ +/* + * Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration + */ + +#ifndef MMT_ROAD_H +#define MMT_ROAD_H + +#include "AthenaBaseComps/AthMessaging.h" +#include "MMT_Hit.h" +#include <numeric> + +namespace MuonGM { + class MuonDetectorManager; +} + +struct micromegas_t { + int roadSize{-1}; + int nstrip_up_XX{-1}; + int nstrip_dn_XX{-1}; + int nstrip_up_UV{-1}; + int nstrip_dn_UV{-1}; + unsigned int strips{0}; + unsigned int nMissedTopEta{0}; + unsigned int nMissedBottomEta{0}; + unsigned int nMissedTopStereo{0}; + unsigned int nMissedBottomStereo{0}; + double pitch{0.}; + double dimensions_top{0.}; + double dimensions_bottom{0.}; + double dimensions_height{0.}; + double activeArea_top{0.}; + double activeArea_bottom{0.}; + double activeArea_height{0.}; + double innerRadiusEta1{0.}; + double innerRadiusEta2{0.}; + std::vector<double> stereoAngles{}; + std::vector<ROOT::Math::XYZVector> planeCoordinates{}; +}; + +class MMT_Road : public AthMessaging { + public: + MMT_Road(const char sector, const MuonGM::MuonDetectorManager* detManager, const micromegas_t mm, int xthr, int uvthr, int iroadx, int iroadu = -1, int iroadv = -1); + ~MMT_Road(); + + void addHits(std::vector<MMT_Hit*> &hits); + double avgSofX(); + double avgSofUV(const int uv1, const int uv2); + double avgZofUV(const int uv1, const int uv2); + bool checkCoincidences(const int &bcwind); + bool containsNeighbors(const MMT_Hit &hit); + unsigned int countHits() const { return m_road_hits.size(); } + unsigned int countRealHits(); + unsigned int countUVHits(bool flag); + unsigned int countXHits(bool flag); + bool horizontalCheck(); + void incrementAge(const int &bcwind); + double getB() const { return m_B; } + int getRoadSize() const { return m_roadSize; } + int getRoadSizeUpX() const { return m_roadSizeUpX; } + int getRoadSizeDownX() const { return m_roadSizeDownX; } + int getRoadSizeUpUV() const { return m_roadSizeUpUV; } + int getRoadSizeDownUV() const { return m_roadSizeDownUV; } + std::vector<MMT_Hit> getHitVector() const { return m_road_hits; } + char getSector() const { return m_sector; } + int getXthreshold() const { return m_xthr; } + int getUVthreshold() const { return m_uvthr; } + int iRoad() const { return m_iroad; } + int iRoadx() const { return m_iroadx; } + int iRoadu() const { return m_iroadu; } + int iRoadv() const { return m_iroadv; } + bool matureCheck(const int &bcwind); + double mxl(); + void reset(); + bool stereoCheck(); + micromegas_t getMM() const { return m_micromegas; } + + private: + const MuonGM::MuonDetectorManager* m_detManager{}; //!< MuonDetectorManager + const MuonGM::MuonDetectorManager* GetDetManager() { return m_detManager; } + double m_B; + int m_iroad; + int m_iroadx; + int m_iroadu; + int m_iroadv; + micromegas_t m_micromega; + char m_sector; + int m_xthr, m_uvthr; + int m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV; + bool m_trig; + + micromegas_t m_micromegas; + std::vector<MMT_Hit> m_road_hits; +}; +#endif diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_struct.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_struct.h index c1303670be43a91034d45a935acba90d328ff06f..8605a8d949f1ebe066339332fdacfebe0b37ff40 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_struct.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMT_struct.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #ifndef MM_STRUCT_H @@ -14,7 +14,9 @@ #include "AthenaBaseComps/AthMessaging.h" -#include "TLorentzVector.h" +#include <Math/Vector2D.h> +#include <Math/Vector3D.h> +#include <Math/Vector4D.h> #include "TMath.h" //flags @@ -56,7 +58,7 @@ public: scale = static_cast<int>( ::truncf(std::log((std::pow(2., static_cast<int>(nBits)) - 1.)/absVal)*(1./std::log(2.))) ); } //=== return input value with fixed point precision - return ::roundf(m_fixp_content * std::pow(2., scale)) / std::pow(2., scale); + return std::roundf(m_fixp_content * std::pow(2., scale)) / std::pow(2., scale); } @@ -115,7 +117,7 @@ struct std_align{ //so let's do show the correspondence between par_cor_mis and this new std_align, struct (we only do the back quadruplet) //(--,-dy,dz;dt,--,dp)=(s,z,t;gamma,beta,alpha) //our internal stuff will still match the (x,y,z) coordinate scheme used in the algorithm; we will just translate results into the standard coordinates for plots, etc. - std_align(int qcm=0,const TVector3& trans=TVector3(),const TVector3& ang=TVector3()); + std_align(int qcm=0,const ROOT::Math::XYZVector& trans=ROOT::Math::XYZVector(),const ROOT::Math::XYZVector& ang=ROOT::Math::XYZVector()); std::string par_title(int par_num,bool small_unit=false) const; std::string par_name(int par_num) const; std::string par_title_val(int par_num)const; @@ -129,7 +131,7 @@ struct std_align{ //members int type;//corresponds to qcm; 0 (nominal case: no misal/correction), 1 (misalignment), 2 (correction), 3 (sim_correction; do corrections based on simulation results) - TVector3 translate,rotate; + ROOT::Math::XYZVector translate,rotate; }; @@ -217,8 +219,6 @@ class MMT_Parameters : public AthMessaging { par_par param_par() const; double y_from_eta_wedge(double eta,int plane)const; double eta_wedge_from_y(double y,int plane)const; - int ybin(float32fixed<18> y,int plane=0)const; - int ybin(float32fixed<yzdex> y,int plane=0)const; int ybin(double y,int plane=0)const; //fill the tables @@ -230,6 +230,7 @@ class MMT_Parameters : public AthMessaging { void fill_crep_table(const std::string&dir,const std::string&tag); void fill_yzmod(); void index_key_test(); + char getSector() const { return sector; } //eta-phi stuff int eta_bin(double theta) const; @@ -255,16 +256,16 @@ class MMT_Parameters : public AthMessaging { std::string bool_to_hit_str(const std::vector<bool>&track)const; //table - std::map<std::vector<int>,std::pair<float32fixed<2>,float32fixed<2> > > AB_k_local; - std::vector<std::vector<std::vector<float32fixed<zbardex> > > > Ak_local_slim;//[x_hit combo][ybin][case #] - std::vector<std::vector<std::vector<float32fixed<bkdex> > > > Bk_local_slim;//[x_hit combo][ybin][case #] - std::vector<std::vector<std::vector<float32fixed<4> > > > Slope_to_ROI; - std::vector<std::vector<float32fixed<2> > > DT_Factors; + std::map<std::vector<int>,std::pair<double,double > > AB_k_local; + std::vector<std::vector<std::vector<double > > > Ak_local_slim;//[x_hit combo][ybin][case #] + std::vector<std::vector<std::vector<double > > > Bk_local_slim;//[x_hit combo][ybin][case #] + std::vector<std::vector<std::vector<double > > > Slope_to_ROI; + std::vector<std::vector<double > > DT_Factors; //theta, phi, hit code, theta/phi/dtheta //hit code: binary stuff.... //old hit code...%mis X, %mis UV: 2-4 X, 1-4 UV possible fractions: 0,1/2,1/3,2/3,1/4,(3/4,not possible with only one misaligned multiplet), 1: 0,3,4,6,8,12 std::vector<std::vector<std::vector<std::vector<float> > > >crep_table; - std::vector<std::vector<std::vector<float32fixed<yzdex> > > >ymod,zmod; + std::vector<std::vector<std::vector<double > > >ymod,zmod; //a toggle bool fill0; @@ -274,12 +275,13 @@ class MMT_Parameters : public AthMessaging { std::vector<double>m_etabins,m_phibins; //currently configurable parameters bool diag,dlm_new; - float32fixed<2> h; + double h; int CT_x,CT_uv; - float32fixed<2> uv_error; + double uv_error; double dtheta_cut; std::string setup; bool islarge,genbg; + char sector; double chargeThreshold; //new, standardized, misalignment and correction information std_align misal,correct; @@ -289,42 +291,42 @@ class MMT_Parameters : public AthMessaging { bool misalign,val_tbl; //dimensions - float32fixed<18> w1, w2, w3, h1, h2, h3, H, Hnom, L, wedge_opening_angle; - float32fixed<4> strip_width; - float32fixed<4> stereo_degree; + double w1, w2, w3, h1, h2, h3, H, Hnom, L, wedge_opening_angle; + double strip_width; + double stereo_degree; double stereo_strip_separation_top; double stereo_strip_separation_bottom; - std::vector<float32fixed<18> > z_nominal; - std::vector<std::vector<float32fixed<18> > > z_large;//[y bin][plane] - std::vector<std::vector<float32fixed<18> > > ybases;//by stationEta--saved from file, hardcoded, alternative is equally spaced, in MMT_Loader::Get_Strip_Id - float32fixed<2> m_x_min,m_x_max,m_y_min,m_y_max,h_mx, h_my; + std::vector<double > z_nominal; + std::vector<std::vector<double > > z_large;//[y bin][plane] + std::vector<std::vector<double > > ybases;//by stationEta--saved from file, hardcoded, alternative is equally spaced, in MMT_Loader::Get_Strip_Id + double m_x_min,m_x_max,m_y_min,m_y_max,h_mx, h_my; int n_x,n_y; - float32fixed<3> slope_min, slope_max; - float32fixed<2> x_error; + double slope_min, slope_max; + double x_error; int CT, CT_u, CT_v; - float32fixed<4> minimum_large_theta, maximum_large_theta; - float32fixed<4> minimum_large_phi, maximum_large_phi; + double minimum_large_theta, maximum_large_theta; + double minimum_large_phi, maximum_large_phi; int n_theta_rois, n_phi_rois, BC_window; - float32fixed<18> mid_plane_large_X, mid_plane_large, mid_plane_large_UV; - float32fixed<4> vertical_strip_width_UV; + double mid_plane_large_X, mid_plane_large, mid_plane_large_UV; + double vertical_strip_width_UV; }; struct mm_digit_entry{ - int multiplet,gas_gap; - double global_time,time; - TVector2 local_pos, stripl_pos; - TVector3 stripg_pos; - double charge; + int multiplet{0},gas_gap{0}; + double global_time{0.},time{0.}; + ROOT::Math::XYVector local_pos{0.,0.}, stripl_pos{0.,0.}; + ROOT::Math::XYZVector stripg_pos{0.,0.,0.}; + double charge{0.}; }; struct evInf_entry{ evInf_entry(int event=0,int pdg=0,double e=0,double p=0,double ieta=0,double peta=0,double eeta=0,double iphi=0,double pphi=0,double ephi=0, - double ithe=0,double pthe=0,double ethe=0,double dth=0,int trn=0,int mun=0,const TVector3& tex=TVector3(), + double ithe=0,double pthe=0,double ethe=0,double dth=0,int trn=0,int mun=0,const ROOT::Math::XYZVector& tex=ROOT::Math::XYZVector(), int troi=0,int antev=0,int postv=0,int nxh=0,int nuvh=0,int nxbg=0,int nuvbg=0,double adt=0,int difx=0, int difuv=0,bool cut=false,bool bad=false); void print() const; @@ -332,7 +334,7 @@ struct evInf_entry{ int athena_event,pdg_id; double E,pt,eta_ip,eta_pos,eta_ent,phi_ip,phi_pos,phi_ent,theta_ip,theta_pos,theta_ent,dtheta; int truth_n,mu_n; - TVector3 vertex; + ROOT::Math::XYZVector vertex; int truth_roi; int N_hits_preVMM,N_hits_postVMM/*<--true*/; int N_X_hits,N_UV_hits;//signal pre vmm; the X and UV hits that count in efficiency denominator @@ -368,18 +370,18 @@ struct hitData_key{ }; struct evFit_entry{ - evFit_entry(int event=0,float32fixed<4> fthe=0.,float32fixed<4> fphi=0., - float32fixed<2> fdth=0.,int roi=-1,int xhit=false,int uvhit=false, - int bgx=false,int bguv=false,float32fixed<2> dth_nd=0.,int hc=0,int tph=0, + evFit_entry(int event=0,double fthe=0.,double fphi=0., + double fdth=0.,int roi=-1,int xhit=false,int uvhit=false, + int bgx=false,int bguv=false,double dth_nd=0.,int hc=0,int tph=0, int bgph=0); void print() const; int athena_event; - float32fixed<4> fit_theta,fit_phi; - float32fixed<2> fit_dtheta; + double fit_theta,fit_phi; + double fit_dtheta; int fit_roi,X_hits_in_fit,UV_hits_in_fit,bg_X_fit,bg_UV_fit; - float32fixed<2> dtheta_nodiv; + double dtheta_nodiv; int hcode,truth_planes_hit,bg_planes_hit; std::vector<hitData_key> fit_hit_keys; @@ -387,15 +389,15 @@ struct evFit_entry{ }; struct evAna_entry{ - int athena_event; - double theta_err,phi_err,dtheta_err; - bool qualified_event, fit, good_fit; - int X_hits_in_fit,UV_hits_in_fit; - int bg_X_fit,bg_UV_fit; + int athena_event{0}; + double theta_err{0.},phi_err{0.},dtheta_err{0.}; + bool qualified_event{false}, fit{false}, good_fit{false}; + int X_hits_in_fit{0},UV_hits_in_fit{0}; + int bg_X_fit{0},bg_UV_fit{0}; }; struct hitData_info{ - hitData_info(int plane,int station_eta,int strip,MMT_Parameters *par,const TVector3& tru,double tpos,double ppos); + hitData_info(int plane,int station_eta,int strip,MMT_Parameters *par,const ROOT::Math::XYZVector& tru,double tpos,double ppos); hitData_info(int the_pl=0,double the_y=0,double the_z=-999); double mis_dy(int pl,MMT_Parameters *m_par,double tpos,double ppos)const; std::string hdr()const; @@ -405,8 +407,8 @@ struct hitData_info{ //members int plane; - float32fixed<yzdex> y,z;//actual values divided by store_const() to make fixed point calculations doable--all this stuff is dimensionless in the end, so it's okay. - float32fixed<2> slope; + double y,z;//actual values divided by store_const() to make fixed point calculations doable--all this stuff is dimensionless in the end, so it's okay. + double slope; }; @@ -423,29 +425,29 @@ struct Hit{ }; struct hitData_entry{ - hitData_entry(int ev=0, double gt=0, double q=0, int vmm=0, int pl=0, int st=0, int est=0, double tr_the=0, double tru_phi=0, - bool q_tbg=0, int bct=0, double time=0,const TVector3& tru=TVector3(), const TVector3& rec=TVector3(), + hitData_entry(int ev=0, double gt=0, double q=0, int vmm=0, int mmfe=0, int pl=0, int st=0, int est=0, int phi=0, int mult=0, int gg=0, double locX=0, double tr_the=0, double tru_phi=0, + bool q_tbg=0, int bct=0, double time=0,const ROOT::Math::XYZVector& tru=ROOT::Math::XYZVector(), const ROOT::Math::XYZVector& rec=ROOT::Math::XYZVector(), double fit_the=0, double fit_phi=0, double fit_dth=0, double tru_dth=0,// double tru_thl=0, double tru_thg=0, double mxg=0, double mug=0, double mvg=0, double mxl=0, double the_mx=0, double the_my=0, int the_roi=0); Hit entry_hit(MMT_Parameters *m_par)const; hitData_key entry_key() const; hitData_info entry_info(MMT_Parameters *m_par)const; - void fit_fill(float32fixed<4> fthe,float32fixed<4> fphi, float32fixed<2> fdth, float32fixed<2> mxg=0., float32fixed<2> mug=0., float32fixed<2> mvg=0., float32fixed<2> mxl=0., float32fixed<2> m_x=0., float32fixed<2> m_y=0., int king=0); + void fit_fill(double fthe,double fphi, double fdth, double mxg=0., double mug=0., double mvg=0., double mxl=0., double m_x=0., double m_y=0., int king=0); void print() const; int event; double gtime,charge; - int VMM_chip,plane,strip,station_eta; - double tru_theta_ip,tru_phi_ip; + int VMM_chip,MMFE_VMM,plane,strip,station_eta,station_phi,multiplet,gasgap; + double localX,tru_theta_ip,tru_phi_ip; bool truth_nbg;//truth (i.e. not bg) if true, int BC_time;/*fit theta, phi, dtheta originally here*/ double time; - TVector3 truth,recon; - float32fixed<4> fit_theta,fit_phi; - float32fixed<2> fit_dtheta; + ROOT::Math::XYZVector truth,recon; + double fit_theta,fit_phi; + double fit_dtheta; double tru_dtheta;//,tru_theta_local,tru_theta_global; - float32fixed<2> M_x_global,M_u_global,M_v_global,M_x_local,mx,my; + double M_x_global,M_u_global,M_v_global,M_x_local,mx,my; int roi; @@ -469,50 +471,52 @@ struct ROI{ ROI(double the_theta, double the_phi, double the_m_x, double the_m_y, int the_roi); //the members: - float32fixed<4> theta,phi; - float32fixed<2> m_x,m_y; + double theta,phi; + double m_x,m_y; int roi; }; struct athena_header{ //make a well-behaved constructor - athena_header(const TLorentzVector& par=TLorentzVector(), int tpn=0, double etp=0, double ete=0, double php=0, double phe=0, int mun=0, const TVector3& ver=TVector3()); + athena_header(const ROOT::Math::PtEtaPhiEVector& par=ROOT::Math::PtEtaPhiEVector(), int tpn=0, double etp=0, double ete=0, double php=0, double phe=0, int mun=0, const ROOT::Math::XYZVector& ver=ROOT::Math::XYZVector()); //the members: - TLorentzVector the_part; + ROOT::Math::PtEtaPhiEVector the_part; int trupart_n; double etapos,etaent,phipos,phient; int muent_n; - TVector3 vertex; + ROOT::Math::XYZVector vertex; }; struct digitWrapper{ digitWrapper(const MmDigit* digit=0, + const std::string& stationName=std::string(), double tmpGTime=0, - const TVector3& truthLPos=TVector3(), - const TVector3& stripLPos=TVector3(), - const TVector3& stripGPos=TVector3() + const ROOT::Math::XYZVector& truthLPos=ROOT::Math::XYZVector(), + const ROOT::Math::XYZVector& stripLPos=ROOT::Math::XYZVector(), + const ROOT::Math::XYZVector& stripGPos=ROOT::Math::XYZVector() ); const MmDigit* digit; + std::string stName; double gTime; - TVector3 truth_lpos;//4,5 - TVector3 strip_lpos; - TVector3 strip_gpos;//6-11 + ROOT::Math::XYZVector truth_lpos;//4,5 + ROOT::Math::XYZVector strip_lpos; + ROOT::Math::XYZVector strip_gpos;//6-11 - inline Identifier id(){ return digit->identify(); }; + inline Identifier id() const { return digit->identify(); }; }; struct track_address{//cartesian_hit_rot entry - track_address(int bct=0,bool big=false,int wed=0,int pl=-1,int sh=0,const TVector3& chr=TVector3()); + track_address(int bct=0,bool big=false,int wed=0,int pl=-1,int sh=0,const ROOT::Math::XYZVector& chr=ROOT::Math::XYZVector()); int BC; bool islarge; int wedge,plane,strip_hit; - TVector3 cart_hit; + ROOT::Math::XYZVector cart_hit; }; #endif diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMTriggerTool.h b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMTriggerTool.h index 4ade632b0c32a33832564cb572cf2bc37ff2da34..6a841f0fe9b49bbabd824f5f0e449a73dffc3aaf 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMTriggerTool.h +++ b/Trigger/TrigT1/TrigT1NSWSimTools/TrigT1NSWSimTools/MMTriggerTool.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #ifndef MMTRIGGERTOOL_H @@ -11,8 +11,7 @@ //local includes #include "TrigT1NSWSimTools/IMMTriggerTool.h" - - +#include "MMT_Diamond.h" //forward declarations class IIncidentSvc; @@ -45,6 +44,8 @@ namespace NSWL1 { MMT_Parameters *m_par_large; MMT_Parameters *m_par_small; + MMT_Diamond *diamond; + //MMT_Loader stuff end MMTriggerTool(const std::string& type, @@ -57,7 +58,9 @@ namespace NSWL1 { virtual void handle (const Incident& inc); - StatusCode runTrigger(); + StatusCode initDiamondAlgorithm(); + StatusCode runTrigger(const bool do_MMDiamonds); + StatusCode finalizeDiamondAlgorithm(const bool do_MMDiamonds); private: @@ -72,6 +75,7 @@ namespace NSWL1 { void clear_ntuple_variables(); //!< clear the variables used in the analysis ntuple void fillNtuple(const MMLoadVariables& loadedVariables); + void storeEventProperties(const unsigned int iterator); //!< Store event variable for Diamond Road algorithm // properties: container and service names StringProperty m_MmDigitContainer; //!< property, see @link MMStripTdsOfflineTool::MMStripTdsOfflineTool @endlink @@ -79,6 +83,31 @@ namespace NSWL1 { BooleanProperty m_doNtuple; //!< property, see @link MMStripTdsOfflineTool::MMStripTdsOfflineTool @endlink TTree* m_tree; //!< ntuple for analysis + std::vector<unsigned int>* m_trigger_diamond_ntrig; + std::vector<int>* m_trigger_diamond_bc; + std::vector<char>* m_trigger_diamond_sector; + std::vector<int>* m_trigger_diamond_stationPhi; + std::vector<unsigned int>* m_trigger_diamond_totalCount; + std::vector<unsigned int>* m_trigger_diamond_realCount; + std::vector<int>* m_trigger_diamond_iX; + std::vector<int>* m_trigger_diamond_iU; + std::vector<int>* m_trigger_diamond_iV; + std::vector<unsigned int>* m_trigger_diamond_XbkgCount; + std::vector<unsigned int>* m_trigger_diamond_UVbkgCount; + std::vector<unsigned int>* m_trigger_diamond_XmuonCount; + std::vector<unsigned int>* m_trigger_diamond_UVmuonCount; + std::vector<double>* m_trigger_diamond_age; + std::vector<double>* m_trigger_diamond_Xavg; + std::vector<double>* m_trigger_diamond_Uavg; + std::vector<double>* m_trigger_diamond_Vavg; + std::vector<double>* m_trigger_diamond_mxl; + std::vector<double>* m_trigger_diamond_theta; + std::vector<double>* m_trigger_diamond_eta; + std::vector<double>* m_trigger_diamond_dtheta; + std::vector<double>* m_trigger_diamond_phi; + std::vector<double>* m_trigger_diamond_phiShf; + + std::vector<double>* m_trigger_RZslopes; std::vector<double>* m_trigger_fitThe; std::vector<double>* m_trigger_fitPhi; std::vector<double>* m_trigger_fitDth; @@ -95,6 +124,9 @@ namespace NSWL1 { std::vector<double>* m_trigger_large_fitDth; std::vector<double>* m_trigger_large_trueEtaRange; std::vector<double>* m_trigger_large_truePtRange; + std::vector<double>* m_trigger_large_trueThe; + std::vector<double>* m_trigger_large_truePhi; + std::vector<double>* m_trigger_large_trueDth; std::vector<double>* m_trigger_large_fitEtaRange; std::vector<double>* m_trigger_large_fitPtRange; std::vector<double>* m_trigger_large_resThe; @@ -106,6 +138,9 @@ namespace NSWL1 { std::vector<double>* m_trigger_small_fitDth; std::vector<double>* m_trigger_small_trueEtaRange; std::vector<double>* m_trigger_small_truePtRange; + std::vector<double>* m_trigger_small_trueThe; + std::vector<double>* m_trigger_small_truePhi; + std::vector<double>* m_trigger_small_trueDth; std::vector<double>* m_trigger_small_fitEtaRange; std::vector<double>* m_trigger_small_fitPtRange; std::vector<double>* m_trigger_small_resThe; diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMLoadVariables.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMLoadVariables.cxx index e0be7ccd2cc44744eb0f4dcce70e9efe3e234e2d..af71007d3a77d0ab768ccfa8ce48f8e5a378922b 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMLoadVariables.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMLoadVariables.cxx @@ -17,26 +17,30 @@ #include "EventInfo/EventID.h" #include "StoreGate/StoreGateSvc.h" #include "MuonIdHelpers/MmIdHelper.h" -#include "AthenaBaseComps/AthCheckMacros.h" #include "AthenaKernel/getMessageSvc.h" +#include "AthenaBaseComps/AthAlgorithm.h" -#include "TVector3.h" +#include "Math/Vector4D.h" #include <cmath> #include <stdexcept> +#include "FourMomUtils/xAODP4Helpers.h" + using std::map; using std::vector; using std::string; -MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetectorManager* detManager, const MmIdHelper* idhelper, MMT_Parameters *par): +MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetectorManager* detManager, const MmIdHelper* idhelper, std::map<std::string,MMT_Parameters*> pars): AthMessaging(Athena::getMessageSvc(), "MMLoadVariables") { - m_par = par; + m_pars = pars; m_evtStore = evtStore; m_detManager = detManager; m_MmIdHelper = idhelper; } - StatusCode MMLoadVariables::getMMDigitsInfo(vector<digitWrapper>& entries, map<hitData_key,hitData_entry>& Hits_Data_Set_Time, map<int,evInf_entry>& Event_Info){ + StatusCode MMLoadVariables::getMMDigitsInfo(map<std::pair<int,unsigned int>,std::vector<digitWrapper> >& entries, + map<std::pair<int,unsigned int>,map<hitData_key,hitData_entry> >& Hits_Data_Set_Time, + map<std::pair<int,unsigned int>,evInf_entry>& Event_Info) { //*******Following MuonPRD code to access all the variables********** histogramVariables fillVars; @@ -57,20 +61,21 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec const MmDigitContainer *nsw_MmDigitContainer = nullptr; ATH_CHECK( m_evtStore->retrieve(nsw_MmDigitContainer,"MM_DIGITS") ); - std::string wedgeType = getWedgeType(nsw_MmDigitContainer); - - float phiEntry = 0; - float phiPosition = 0; - float etaEntry = 0; - float etaPosition = 0; + std::vector<ROOT::Math::PtEtaPhiEVector> truthParticles, truthParticles_ent, truthParticles_pos; + std::vector<int> pdg; + std::vector<ROOT::Math::XYZVector> vertex; + float phiEntry_tmp = 0; + float phiPosition_tmp = 0; + float etaEntry_tmp = 0; + float etaPosition_tmp = 0; + int pdg_tmp = 0; + ROOT::Math::XYZVector vertex_tmp(0.,0.,0.); - TLorentzVector thePart, theInfo; - TVector3 vertex; - int pdg=0; + ROOT::Math::PtEtaPhiEVector thePart, theInfo; auto MuEntry_Particle_n = trackRecordCollection->size(); - int j=0; //# iteration of particle entries + int j=0; // iteration of particle entries - for(const auto it : *truthContainer) { //first loop in MMT_loader::load_event + for(const auto it : *truthContainer) { const HepMC::GenEvent *subEvent = it; #ifdef HEPMC3 for(auto particle : subEvent->particles()) { @@ -80,23 +85,23 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec const HepMC::GenParticle *particle = pit; #endif const HepMC::FourVector momentum = particle->momentum(); - int k=trackRecordCollection->size(); //number of mu entries - if(HepMC::barcode(particle) == 10001 && std::abs(particle->pdg_id())==13){ - thePart.SetPtEtaPhiE(momentum.perp(),momentum.eta(),momentum.phi(),momentum.e()); + int k=trackRecordCollection->size(); //number of mu + if(particle->barcode() < 1e06 && std::abs(particle->pdg_id())==13){ + thePart.SetCoordinates(momentum.perp(),momentum.eta(),momentum.phi(),momentum.e()); for(const auto & mit : *trackRecordCollection ) { - if(k>0&&j<k){ - const CLHEP::Hep3Vector mumomentum = mit.GetMomentum(); - const CLHEP::Hep3Vector muposition = mit.GetPosition(); - pdg=HepMC::barcode(particle); - phiEntry = mumomentum.getPhi(); - etaEntry = mumomentum.getEta(); - phiPosition = muposition.getPhi(); - etaPosition = muposition.getEta(); + const CLHEP::Hep3Vector mumomentum = mit.GetMomentum(); + const CLHEP::Hep3Vector muposition = mit.GetPosition(); + if(k>0 && j<k && particle->barcode()==mit.GetBarCode()) { + pdg_tmp = HepMC::barcode(particle); + phiEntry_tmp = mumomentum.getPhi(); + etaEntry_tmp = mumomentum.getEta(); + phiPosition_tmp = muposition.getPhi(); + etaPosition_tmp = muposition.getEta(); } }//muentry loop - int l=0; + int l=0; #ifdef HEPMC3 - for(auto vertex1 : subEvent->vertices()) { + for(auto vertex1 : subEvent->vertices()) { #else HepMC::ConstGenEventVertexRange vertex_range = subEvent->vertex_range(); for(const auto vit : vertex_range) { @@ -104,34 +109,51 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec const HepMC::GenVertex *vertex1 = vit; #endif const HepMC::FourVector& position = vertex1->position(); - vertex=TVector3(position.x(),position.y(),position.z()); + vertex_tmp.SetXYZ(position.x(),position.y(),position.z()); l++; }//end vertex loop } j++; + + if(thePart.Pt() > 0. && particle->barcode() < 1e06){ + bool addIt = true; + for(unsigned int ipart=0; ipart < truthParticles.size(); ipart++){ + if( std::abs(thePart.Pt()-truthParticles[ipart].Pt()) < 0.001 || + std::abs(thePart.Eta()-truthParticles[ipart].Eta()) < 0.001 || + std::abs(xAOD::P4Helpers::deltaPhi(thePart.Phi(), truthParticles[ipart].Phi())) < 0.001 || + std::abs(thePart.E()-truthParticles[ipart].E()) < 0.001 ) addIt = false; + } + if(addIt){ + truthParticles.push_back(thePart); + //new stuff + vertex.push_back(vertex_tmp); + pdg.push_back(pdg_tmp); + truthParticles_ent.push_back(ROOT::Math::PtEtaPhiEVector(momentum.perp(),etaEntry_tmp ,phiEntry_tmp ,momentum.e())); + truthParticles_pos.push_back(ROOT::Math::PtEtaPhiEVector(momentum.perp(),etaPosition_tmp,phiPosition_tmp,momentum.e())); + } + } + } //end particle loop - } //end container loop (should be only 1 container per event) + } //end truth container loop (should be only 1 container per event) const EventInfo* pevt = 0; ATH_CHECK( m_evtStore->retrieve(pevt) ); int event = pevt->event_ID()->event_number(); int TruthParticle_n = j; - evFit_entry fit;fit.athena_event=event;//-1; **May not be necessary - - //Truth information for theta, phi - double theta_pos=std::atan(std::exp(-etaPosition))*2; - double theta_ent=std::atan(std::exp(-etaEntry))*2; - double phi_pos=phiPosition; + evFit_entry fit; + fit.athena_event=event; //get hits container const MMSimHitCollection *nswContainer = nullptr; ATH_CHECK( m_evtStore->retrieve(nswContainer,"MicromegasSensitiveDetector") ); + unsigned int digit_particles = 0; for(auto digitCollectionIter : *nsw_MmDigitContainer) { // a digit collection is instanciated for each container, i.e. holds all digits of a multilayer const MmDigitCollection* digitCollection = digitCollectionIter; // loop on all digits inside a collection, i.e. multilayer int digit_count =0; + std::vector<digitWrapper> entries_tmp; for (const auto item:*digitCollection) { // get specific digit and identify it @@ -148,6 +170,7 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec int gas_gap = m_MmIdHelper->gasGap(id); int channel = m_MmIdHelper->channel(id); + if (stName.substr(0,2) != "MM") continue; int isSmall = (stName[2] == 'S'); const MuonGM::MMReadoutElement* rdoEl = m_detManager->getMMRElement_fromIdFields(isSmall, stationEta, stationPhi, multiplet ); @@ -166,21 +189,13 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec fillVars.NSWMM_dig_gas_gap.push_back(gas_gap); fillVars.NSWMM_dig_channel.push_back(channel); - //match to truth particle - TLorentzVector truthPart; - for(auto it1 : *truthContainer) { //Must be a more elegant way... should work for now though - for(auto particle1 : *it1) { - const HepMC::FourVector momentum1 = particle1->momentum(); - truthPart.SetPtEtaPhiE(momentum1.perp(),momentum1.eta(),momentum1.phi(),momentum1.e()); - }//end particle loop - }//end truth container loop (1 iteration) for matching - std::vector<double> localPosX; std::vector<double> localPosY; std::vector<double> globalPosX; std::vector<double> globalPosY; std::vector<double> globalPosZ; + int nstrip = 0; //counter of the number of firing strips for (const auto &i: stripPosition) { // take strip index form chip information int cr_strip = i; @@ -189,6 +204,7 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec globalPosX.push_back(0.); globalPosY.push_back(0.); globalPosZ.push_back(0.); + ++nstrip; Identifier cr_id = m_MmIdHelper->channelID(stationName, stationEta, stationPhi, multiplet, gas_gap, cr_strip, true, &isValid); if (!isValid) { @@ -201,22 +217,19 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec ATH_MSG_WARNING("MicroMegas digitization: failed to associate a valid local position for (chip response) strip n. " << cr_strip << "; associated positions will be set to 0.0."); } else { - localPosX.at(i) = cr_strip_pos.x(); - localPosY.at(i) = cr_strip_pos.y(); + localPosX[nstrip-1] = cr_strip_pos.x(); + localPosY[nstrip-1] = cr_strip_pos.y(); } // asking the detector element to transform this local to the global position Amg::Vector3D cr_strip_gpos(0., 0., 0.); rdoEl->surface(cr_id).localToGlobal(cr_strip_pos, Amg::Vector3D(0., 0., 0.), cr_strip_gpos); - globalPosX.at(i) = cr_strip_gpos[0]; - globalPosY.at(i) = cr_strip_gpos[1]; - globalPosZ.at(i) = cr_strip_gpos[2]; - - } - + globalPosX[nstrip-1] = cr_strip_gpos[0]; + globalPosY[nstrip-1] = cr_strip_gpos[1]; + globalPosZ[nstrip-1] = cr_strip_gpos[2]; + } }//end of strip position loop - //NTUPLE FILL DIGITS fillVars.NSWMM_dig_time.push_back(time); fillVars.NSWMM_dig_charge.push_back(charge); @@ -267,8 +280,6 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec int sim_layer = hitHelper->GetLayer(simId); int sim_side = hitHelper->GetSide(simId); - // Fill Ntuple with SimId data - //fillVars.NSWMM_sim_stationName .push_back(sim_stationName); if(digit_count){ fillVars.NSWMM_sim_stationEta .push_back(sim_stationEta); fillVars.NSWMM_sim_stationPhi .push_back(sim_stationPhi); @@ -279,204 +290,165 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec if(hit.depositEnergy()==0.) continue; // SimHits without energy loss are not recorded. if(digit_count==hit_count) { - entries.push_back( + entries_tmp.emplace_back( digitWrapper(digit, + stName, hit.globalTime(), - TVector3(-999, -99, -999),//Digits_MM_truth_localPosZ->at(i)), - TVector3(localPosX.at(indexOfFastestSignal), - localPosY.at(indexOfFastestSignal), - -999), //Digits_MM_stripLposZ->at(i).at(indexOfFastestSignal)), - TVector3(globalPosX.at(indexOfFastestSignal), - globalPosY.at(indexOfFastestSignal), - globalPosZ.at(indexOfFastestSignal) ) + ROOT::Math::XYZVector(-999, -999, -999), + ROOT::Math::XYZVector(localPosX[indexOfFastestSignal], + localPosY[indexOfFastestSignal], + -999), + ROOT::Math::XYZVector(globalPosX[indexOfFastestSignal], + globalPosY[indexOfFastestSignal], + globalPosZ[indexOfFastestSignal] ) ) - ); + ); } hit_count++; - }//end of hit cotainer loop */ - // } //end if sdo + }//end of hit container loop digit_count++; } //end iterator digit loop - } // end digit container loop (1 for event?) - - vector<digitWrapper> dummy; - vector<int> indices; - //Truth info for particle (originally primer) - int phiSt = 0; - if(entries.size() > 0) phiSt = m_MmIdHelper->stationPhi( entries.at(0).id() ); - evInf_entry particle_info(event,pdg,thePart.E(),thePart.Pt(),thePart.Eta(),etaPosition,etaEntry,phi_shift(thePart.Phi(),wedgeType,phiSt),phi_shift(phi_pos,wedgeType,phiSt),phi_shift(phiEntry,wedgeType,phiSt),thePart.Theta(),theta_pos,theta_ent,theta_ent-theta_pos,TruthParticle_n,MuEntry_Particle_n,vertex); - if(wedgeType == "Neither") particle_info.bad_wedge=true; - else particle_info.bad_wedge=false; - - vector<digitWrapper> origEntries = entries; - - for(unsigned int i=0; i<entries.size(); i++){ - if(entries.size() < 8) continue; - Identifier tmpID = entries.at(i).id(); - if ( m_MmIdHelper->multilayer( tmpID )==1 and m_MmIdHelper->gasGap(tmpID)==1) entries.at(i).gTime = origEntries.at(0).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==1 and m_MmIdHelper->gasGap(tmpID)==2) entries.at(i).gTime = origEntries.at(1).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==1 and m_MmIdHelper->gasGap(tmpID)==3) entries.at(i).gTime = origEntries.at(2).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==1 and m_MmIdHelper->gasGap(tmpID)==4) entries.at(i).gTime = origEntries.at(3).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==2 and m_MmIdHelper->gasGap(tmpID)==1) entries.at(i).gTime = origEntries.at(4).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==2 and m_MmIdHelper->gasGap(tmpID)==2) entries.at(i).gTime = origEntries.at(5).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==2 and m_MmIdHelper->gasGap(tmpID)==3) entries.at(i).gTime = origEntries.at(6).gTime; - else if( m_MmIdHelper->multilayer( tmpID )==2 and m_MmIdHelper->gasGap(tmpID)==4) entries.at(i).gTime = origEntries.at(7).gTime; - } - for(unsigned int i=0; i<entries.size(); i++){ - float min = 100000; - int minIndex=-999; - for(unsigned int j=0; j<entries.size(); j++){ - bool notmindex = true; - for (unsigned int k=0; k<indices.size(); k++){ - if(j==(unsigned int) indices[k]) notmindex=false; - } - if(notmindex){ - if(min > entries.at(j).gTime ){ - minIndex = j; - min = entries.at(minIndex).gTime ; - } + + if (!entries_tmp.empty()) { + std::vector<std::string> stNames; + for(const auto &dW : entries_tmp) stNames.push_back(dW.stName); + if(std::all_of(stNames.begin(), stNames.end(), [&] (const std::string & name) { return name == stNames[0]; })) { + entries[std::make_pair(event,digit_particles)]=entries_tmp; + digit_particles++; } + stNames.clear(); } - if(minIndex < 0) { minIndex = i; } - dummy.push_back(entries.at(minIndex)); - indices.push_back(minIndex); + } // end digit container loop + + for(unsigned int i=0; i<truthParticles.size(); i++) { + evInf_entry particle_info(event, pdg[i], + truthParticles[i].E(), truthParticles[i].Pt(), + truthParticles[i].Eta(), truthParticles_pos[i].Eta(), truthParticles_ent[i].Eta(), + truthParticles[i].Phi(), truthParticles_pos[i].Phi(), truthParticles_ent[i].Phi(), + truthParticles[i].Theta(), truthParticles_pos[i].Theta(), truthParticles_ent[i].Theta(), truthParticles_ent[i].Theta()-truthParticles_pos[i].Theta(), + TruthParticle_n,MuEntry_Particle_n,vertex[i]); + particle_info.NUV_bg_preVMM = 0; + Event_Info[std::make_pair(event,i)] = particle_info; } - entries = dummy; - int min_hits = 1,max_hits = 10000,nent=entries.size(); - - m_uvxxmod=(m_par->setup.compare("xxuvuvxx")==0); - //Number of hits cut - if(!particle_info.bad_wedge) particle_info.pass_cut=true; //default is false - if(nent<min_hits||nent>max_hits) particle_info.pass_cut=false; - - double tru_phi = -999; - if (entries.size() >0) tru_phi=phi_shift(thePart.Phi(),wedgeType,phiSt ); - double tru_theta=thePart.Theta(); + for (auto it=entries.begin(); it!=entries.end(); it++) { + std::sort(it->second.begin(), it->second.end(), [](const auto &dW1, const auto &dW2){ return (dW1.gTime < dW2.gTime); }); + } - //Hit information in Stephen's code... Starts getting a little weird. - map<hitData_key,hitData_entry> hit_info; //Originally "targaryen" - vector<hitData_key> keys; //Loop over entries, which has digitization info for each event - for(unsigned int ient=0; ient<entries.size(); ient++){ - digitWrapper thisSignal=entries.at(ient); - Identifier tmpID = thisSignal.id(); - int thisMultiplet = m_MmIdHelper->multilayer( tmpID ); - int thisGasGap = m_MmIdHelper->gasGap( tmpID ); - int thisTime = thisSignal.digit->stripTimeForTrigger().at(0); - int thisCharge = 2; //thisSignal.digit->stripChargeForTrigger().at(0); - int thisStripPosition = thisSignal.digit->stripPositionForTrigger().at(0); - int thisStationEta = m_MmIdHelper->stationEta( tmpID ); - //DLM_NEW plane assignments - //stated [3,2,1,0;7,6,5,4] - int thisPlane = (thisMultiplet-1)*4+thisGasGap-1; - - int strip = strip_number(thisStationEta, - thisPlane, - thisStripPosition); - int thisVMM = Get_VMM_chip(strip); - int BC_id = std::ceil( thisTime / 25. ); - TVector3 mazin_check( - thisSignal.strip_gpos.X(), - thisSignal.strip_gpos.Y(), - thisSignal.strip_gpos.Z() - ); - - TVector3 athena_tru( - thisSignal.strip_gpos.X(), - thisSignal.strip_gpos.Y()-thisSignal.truth_lpos.Y(), - thisSignal.strip_gpos.Z()); - - if(m_par->dlm_new){ - athena_tru.SetX(thisSignal.strip_gpos.X()-thisSignal.strip_lpos.X()); + unsigned int ient=0; + for (auto it=entries.begin(); it!=entries.end(); it++){ + + /* Identifying the wedge from digits: + * now the digit is associated with the corresponding station, so lambda function can be exploited to check if they are all the same + */ + double tru_phi = -999, tru_theta = -999; + std::pair<int, unsigned int> pair (event,ient); + auto tru_it = Event_Info.find(pair); + if (tru_it != Event_Info.end()) { + tru_phi = tru_it->second.phi_pos; + tru_theta = tru_it->second.theta_pos; } - TVector3 athena_rec(thisSignal.strip_gpos); - //now store some variables BLC initializes; we might trim this down for efficiency later - //the following line should be rather easy to one-to-one replace + std::string station = it->second[0].stName; + m_uvxxmod=(m_pars[station]->setup.compare("xxuvuvxx")==0); + + map<hitData_key,hitData_entry> hit_info; + vector<hitData_key> keys; + + //Now we need to loop on digits + for (const auto &dW : it->second) { + Identifier tmpID = dW.id(); + int thisMultiplet = m_MmIdHelper->multilayer( tmpID ); + int thisGasGap = m_MmIdHelper->gasGap( tmpID ); + int thisTime = dW.digit->stripTimeForTrigger()[0]; + int thisCharge = 2; //dW.digit->stripChargeForTrigger().at(0); + int thisStripPosition = dW.digit->stripPositionForTrigger()[0]; + double thisLocalPosX = dW.strip_lpos.X(); + int thisVMM = dW.digit->VMM_idForTrigger()[0]; + int thisMMFE_VMM = dW.digit->MMFE_VMM_idForTrigger()[0]; + int thisStationEta = m_MmIdHelper->stationEta( tmpID ); + int thisStationPhi = m_MmIdHelper->stationPhi( tmpID ); + int thisPlane = (thisMultiplet-1)*4+thisGasGap-1; + int BC_id = std::ceil( thisTime / 25. ); + ROOT::Math::XYZVector mazin_check( + dW.strip_gpos.X(), + dW.strip_gpos.Y(), + dW.strip_gpos.Z() + ); + + ROOT::Math::XYZVector athena_rec(dW.strip_gpos); + ROOT::Math::XYZVector recon(athena_rec.Y(),-athena_rec.X(),athena_rec.Z()); + + if(m_uvxxmod){ + xxuv_to_uvxx(recon,thisPlane,m_pars[station]); + } - TVector3 truth(athena_tru.Y(),-athena_tru.X(),athena_tru.Z()), recon(athena_rec.Y(),-athena_rec.X(),athena_rec.Z()); + //We're doing everything by the variable known as "athena_event" to reflect C++ vs MATLAB indexing + int btime=(event+1)*10+(BC_id-1); + int special_time = thisTime + (event+1)*100; + + hitData_entry hit_entry(event, + dW.gTime, + thisCharge, + thisVMM, + thisMMFE_VMM, + thisPlane, + thisStripPosition, + thisStationEta, + thisStationPhi, + thisMultiplet, + thisGasGap, + thisLocalPosX, + tru_theta, + tru_phi, + true, + btime, + special_time, + mazin_check, + mazin_check); + + hit_info[hit_entry.entry_key()]=hit_entry; + ATH_MSG_DEBUG("Filling hit_info slot: "); + if (msgLvl(MSG::DEBUG)){ + hit_entry.entry_key().print(); + } + keys.push_back(hit_entry.entry_key()); + }//end digit wrapper loop - if(m_uvxxmod){ - xxuv_to_uvxx(truth,thisPlane);xxuv_to_uvxx(recon,thisPlane); + if (tru_it != Event_Info.end()) { + tru_it->second.N_hits_preVMM=hit_info.size(); + tru_it->second.N_hits_postVMM=0; } - //we're doing everything by the variable known as "athena_event" to reflect C++ vs MATLAB indexing - int btime=(event+1)*10+(BC_id-1); - particle_info.NUV_bg_preVMM = 0; //thisSignal.gtime; - int special_time = thisTime + (event+1)*100; - - hitData_entry hit_entry(event, - thisSignal.gTime, - thisCharge, - thisVMM, - thisPlane, - thisStripPosition, - thisStationEta, - tru_theta, - tru_phi, - true, - btime, - special_time, - mazin_check, - mazin_check); - - hit_info[hit_entry.entry_key()]=hit_entry; - ATH_MSG_DEBUG("Filling hit_info slot: "); - if (msgLvl(MSG::DEBUG)){ - hit_entry.entry_key().print(); - } - keys.push_back(hit_entry.entry_key()); //may be only used when "incoherent background" is generated (not included for now) - - }//end entries loop - - particle_info.N_hits_preVMM=hit_info.size(); - particle_info.N_hits_postVMM=0; - // unsigned int ir=0; - - //might want to move these somewhere smarter in future - m_VMM_Deadtime = 100; - m_numVMM_PerPlane = 1000; - m_VMM_ChipStatus=vector<vector<bool> >(m_numVMM_PerPlane,vector<bool>(8,true)); - m_VMM_ChipLastHitTime=vector<vector<int> >(m_numVMM_PerPlane,vector<int>(8,0)); - - //*** FIGURE OUT WHAT TO DO WITH THE TIES--IS MIMIC VMM BUSTED? DO WE PLACE THE HIT REQUIREMENTS HERE? (PROBABLY) - int xhit=0,uvhit=0,strip_X_tot=0,strip_UV_tot=0; - vector<bool>plane_hit(m_par->setup.size(),false); - - for(map<hitData_key,hitData_entry>::iterator it=hit_info.begin(); it!=hit_info.end(); ++it){ - int plane=it->second.plane; - plane_hit[plane]=true; - particle_info.N_hits_postVMM++; - Hits_Data_Set_Time[it->first]=it->second; - if(m_par->setup.substr(plane,1).compare("x")==0){ - ATH_MSG_DEBUG("ADD X STRIP VALUE "<<it->second.strip); - strip_X_tot+=it->second.strip; - } - else{ - strip_UV_tot+=it->second.strip; + //might want to move these somewhere smarter in future + m_VMM_Deadtime = 100; + m_numVMM_PerPlane = 1000; + m_VMM_ChipStatus=vector<vector<bool> >(m_numVMM_PerPlane,vector<bool>(8,true)); + m_VMM_ChipLastHitTime=vector<vector<int> >(m_numVMM_PerPlane,vector<int>(8,0)); + + int xhit=0,uvhit=0; + vector<bool>plane_hit(m_pars[station]->setup.size(),false); + + for(map<hitData_key,hitData_entry>::iterator it=hit_info.begin(); it!=hit_info.end(); ++it){ + int plane=it->second.plane; + plane_hit[plane]=true; + if (tru_it != Event_Info.end()) tru_it->second.N_hits_postVMM++; } - }//end map iterator loop + Hits_Data_Set_Time[std::make_pair(event,ient)] = hit_info; - for(unsigned int ipl=0;ipl<plane_hit.size();ipl++){ - if(plane_hit[ipl]){ - if(m_par->setup.substr(ipl,1)=="x") xhit++; - else if(m_par->setup.substr(ipl,1)=="u"||m_par->setup.substr(ipl,1)=="v") uvhit++; + for(unsigned int ipl=0;ipl<plane_hit.size();ipl++){ + if(plane_hit[ipl]){ + if(m_pars[station]->setup.substr(ipl,1)=="x") xhit++; + else if(m_pars[station]->setup.substr(ipl,1)=="u" || m_pars[station]->setup.substr(ipl,1)=="v") uvhit++; + } } + + histVars = fillVars; + ient++; } - - particle_info.N_X_hits=xhit; - particle_info.N_UV_hits=uvhit; - //X and UV hits minumum cut - if(xhit<m_par->CT_x) particle_info.pass_cut=false;//return; - if(uvhit<m_par->CT_uv) particle_info.pass_cut=false;//return; - - //Moved the removing of some entries to the end of ~Trigger - Event_Info[event]=particle_info; - histVars = fillVars; - - return StatusCode::SUCCESS; - } + return StatusCode::SUCCESS; + } double MMLoadVariables::phi_shift(double athena_phi,const std::string& wedgeType, int stationPhi) const{ float n = 2*(stationPhi-1); @@ -490,16 +462,16 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec else return athena_phi; } - void MMLoadVariables::xxuv_to_uvxx(TVector3& hit,const int plane) const{ + void MMLoadVariables::xxuv_to_uvxx(ROOT::Math::XYZVector& hit,const int plane, const MMT_Parameters *par) const{ if(plane<4)return; - else if(plane==4)hit_rot_stereo_bck(hit);//x to u - else if(plane==5)hit_rot_stereo_fwd(hit);//x to v - else if(plane==6)hit_rot_stereo_fwd(hit);//u to x - else if(plane==7)hit_rot_stereo_bck(hit);//v to x + else if(plane==4)hit_rot_stereo_bck(hit, par);//x to u + else if(plane==5)hit_rot_stereo_fwd(hit, par);//x to v + else if(plane==6)hit_rot_stereo_fwd(hit, par);//u to x + else if(plane==7)hit_rot_stereo_bck(hit, par);//v to x } - void MMLoadVariables::hit_rot_stereo_fwd(TVector3& hit)const{ - double degree=TMath::DegToRad()*(m_par->stereo_degree.getFixed()); + void MMLoadVariables::hit_rot_stereo_fwd(ROOT::Math::XYZVector& hit, const MMT_Parameters *par)const{ + double degree=TMath::DegToRad()*(par->stereo_degree); if(m_striphack) hit.SetY(hit.Y()*cos(degree)); else{ double xnew=hit.X()*std::cos(degree)+hit.Y()*std::sin(degree),ynew=-hit.X()*std::sin(degree)+hit.Y()*std::cos(degree); @@ -507,8 +479,8 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec } } - void MMLoadVariables::hit_rot_stereo_bck(TVector3& hit)const{ - double degree=-TMath::DegToRad()*(m_par->stereo_degree.getFixed()); + void MMLoadVariables::hit_rot_stereo_bck(ROOT::Math::XYZVector& hit, const MMT_Parameters *par)const{ + double degree=-TMath::DegToRad()*(par->stereo_degree); if(m_striphack) hit.SetY(hit.Y()*std::cos(degree)); else{ double xnew=hit.X()*std::cos(degree)+hit.Y()*std::sin(degree),ynew=-hit.X()*std::sin(degree)+hit.Y()*std::cos(degree); @@ -516,35 +488,6 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec } } - - int - MMLoadVariables::Get_Strip_ID(double X,double Y,int plane) const{ //athena_strip_id,module_y_center,plane) - if(Y==-9999) return -1; - string setup(m_par->setup); - double strip_width=m_par->strip_width.getFixed(), degree=TMath::DegToRad()*(m_par->stereo_degree.getFixed());//,vertical_strip_width_UV = strip_width/cos(degree); - double y_hit=Y; - int setl=setup.length(); - if(plane>=setl||plane<0){ - ATH_MSG_FATAL("Pick a plane in [0,"<<setup.length()<<"] not "<<plane); - throw std::runtime_error("MMLoadVariables::Get_Strip_ID: invalid plane"); - } - string xuv=setup.substr(plane,1); - if(xuv=="u"){//||xuv=="v"){ - if(m_striphack)return ceil(Y*cos(degree)/strip_width); - y_hit = X*sin(degree)+Y*cos(degree); - } - else if(xuv=="v"){ - if(m_striphack)return ceil(Y*cos(degree)/strip_width); - y_hit = -X*sin(degree)+Y*cos(degree); - } - else if(xuv!="x"){ - ATH_MSG_FATAL("Invalid plane option " << xuv ); - throw std::runtime_error("MMLoadVariables::Get_Strip_ID: invalid plane"); - } - double strip_hit = ceil(y_hit*1./strip_width); - return strip_hit; - } - int MMLoadVariables::Get_VMM_chip(int strip) const{ //Not Finished... Rough int strips_per_VMM = 64; @@ -552,12 +495,12 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec } int - MMLoadVariables::strip_number(int station,int plane,int spos)const{ - if (station<=0||station>m_par->n_stations_eta) { + MMLoadVariables::strip_number(int station, int plane, int spos, const MMT_Parameters *par)const{ + if (station<=0||station>par->n_stations_eta) { int base_strip = 0; return base_strip; } - if (plane<0||plane>(int)m_par->setup.size()) { + if (plane<0||plane>(int)par->setup.size()) { int base_strip = 0; return base_strip; @@ -566,33 +509,3 @@ MMLoadVariables::MMLoadVariables(StoreGateSvc* evtStore, const MuonGM::MuonDetec int base_strip=spos; return base_strip; } - - std::string - MMLoadVariables::getWedgeType(const MmDigitContainer *nsw_MmDigitContainer){ - std::vector<bool> isLargeWedge; - //Digit loop to match to truth - for(auto digitCollectionIter : *nsw_MmDigitContainer) { - - const MmDigitCollection* digitCollection = digitCollectionIter; - for (unsigned int item=0; item<digitCollection->size(); item++) { - - const MmDigit* digit = digitCollection->at(item); - Identifier id = digit->identify(); - - std::string stName = m_MmIdHelper->stationNameString(m_MmIdHelper->stationName(id)); - const string& sname(stName); - if (sname.compare("MML")==0)isLargeWedge.push_back(true); - else isLargeWedge.push_back(false); - } - } - bool allLarge = true; - bool allSmall = true; - for(unsigned int i=0; i<isLargeWedge.size(); i++){ - if (isLargeWedge.at(i)) allSmall = false; - else allLarge = false; - } - std::string wedgeType = "Neither"; - if (allLarge) wedgeType = "Large"; - if (allSmall) wedgeType = "Small"; - return wedgeType; - } diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMStripTdsOfflineTool.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMStripTdsOfflineTool.cxx index ae616fb99c57f47e1915c9632826f14d51760515..4f3289e4b8fec55fa4ae42b2bf6492ba6c2d89ad 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMStripTdsOfflineTool.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMStripTdsOfflineTool.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ // Athena/Gaudi includes @@ -16,7 +16,6 @@ #include "EventInfo/EventID.h" // ROOT includes -#include "TVector3.h" #include "TMath.h" const bool striphack=true; @@ -563,7 +562,7 @@ namespace NSWL1 { return athena_phi+(athena_phi>=0?-1:1)*TMath::Pi(); return (-1.*athena_phi+(athena_phi>=0?1.5*TMath::Pi():-0.5*TMath::Pi())); } - void MMStripTdsOfflineTool::xxuv_to_uvxx(TVector3& hit,int plane)const{ + void MMStripTdsOfflineTool::xxuv_to_uvxx(ROOT::Math::XYZVector& hit,int plane)const{ if(plane<4)return; else if(plane==4)hit_rot_stereo_bck(hit);//x to u else if(plane==5)hit_rot_stereo_fwd(hit);//x to v @@ -571,8 +570,8 @@ namespace NSWL1 { else if(plane==7)hit_rot_stereo_bck(hit);//v to x } - void MMStripTdsOfflineTool::hit_rot_stereo_fwd(TVector3& hit)const{ - double degree=TMath::DegToRad()*(m_par->stereo_degree.getFixed()); + void MMStripTdsOfflineTool::hit_rot_stereo_fwd(ROOT::Math::XYZVector& hit)const{ + double degree=TMath::DegToRad()*(m_par->stereo_degree); if(striphack) hit.SetY(hit.Y()*cos(degree)); else{ double xnew=hit.X()*cos(degree)+hit.Y()*sin(degree),ynew=-hit.X()*sin(degree)+hit.Y()*cos(degree); @@ -580,8 +579,8 @@ namespace NSWL1 { } } - void MMStripTdsOfflineTool::hit_rot_stereo_bck(TVector3& hit)const{ - double degree=-TMath::DegToRad()*(m_par->stereo_degree.getFixed()); + void MMStripTdsOfflineTool::hit_rot_stereo_bck(ROOT::Math::XYZVector& hit)const{ + double degree=-TMath::DegToRad()*(m_par->stereo_degree); if(striphack) hit.SetY(hit.Y()*cos(degree)); else{ double xnew=hit.X()*cos(degree)+hit.Y()*sin(degree),ynew=-hit.X()*sin(degree)+hit.Y()*cos(degree); @@ -606,7 +605,7 @@ namespace NSWL1 { int MMStripTdsOfflineTool::Get_Strip_ID(double X,double Y,int plane) const{ //athena_strip_id,module_y_center,plane) if(Y==-9999) return -1; std::string setup(m_par->setup); - double strip_width=m_par->strip_width.getFixed(), degree=TMath::DegToRad()*(m_par->stereo_degree.getFixed());//,vertical_strip_width_UV = strip_width/cos(degree); + double strip_width=m_par->strip_width, degree=TMath::DegToRad()*(m_par->stereo_degree);//,vertical_strip_width_UV = strip_width/cos(degree); double y_hit=Y; int setl=setup.length(); if(plane>=setl||plane<0){ @@ -648,7 +647,7 @@ namespace NSWL1 { } bool do_auto=false; //if true do strip # (ceil(Y/strip_width); what's currently fed into the algorithm) calculation based on evenly spaced eta assumption of stations - double H=m_par->H.getFixed()/*,h=m_par->h1,z=m_par->z_nominal[plane],z0=m_par->z_nominal.front()*/,ybase=m_par->ybases[plane][station-1].getFixed(); + double H=m_par->H/*,h=m_par->h1,z=m_par->z_nominal[plane],z0=m_par->z_nominal.front()*/,ybase=m_par->ybases[plane][station-1]; if(do_auto){ //-log(tan(0.5(atan(y/z))))=eta //this is the even y spacing @@ -660,7 +659,7 @@ namespace NSWL1 { ybase=z*tan(2*atan(exp(-1.*this_eta))); */ } - double width=m_par->strip_width.getFixed(); + double width=m_par->strip_width; std::string plane_char=m_par->setup.substr(plane,1); int base_strip=ceil(ybase/width)+spos; return base_strip; diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Diamond.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Diamond.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c3ba01156d0b7cb41d4971200daf5f892630b4f4 --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Diamond.cxx @@ -0,0 +1,306 @@ +#include "TrigT1NSWSimTools/MMT_Diamond.h" +#include "AthenaKernel/getMessageSvc.h" +#include "MuonAGDDDescription/MMDetectorDescription.h" +#include "MuonAGDDDescription/MMDetectorHelper.h" + +MMT_Diamond::MMT_Diamond(const MuonGM::MuonDetectorManager* detManager): AthMessaging(Athena::getMessageSvc(), "MMT_Diamond") { + m_detManager = detManager; +} + +void MMT_Diamond::clearEvent() { + if (!m_diamonds.empty()) { + for (auto &diam : m_diamonds) { + if (!diam.ev_roads.empty()) { + for (auto &road : diam.ev_roads) delete road; + diam.ev_roads.clear(); + } + if (!diam.ev_hits.empty()) { + for (auto &hit : diam.ev_hits) delete hit; + diam.ev_hits.clear(); + } + diam.slopes.clear(); + } + m_diamonds.clear(); + } +} + +void MMT_Diamond::createRoads_fillHits(const unsigned int iterator, std::vector<hitData_entry> &hitDatas, const MuonGM::MuonDetectorManager* detManager, MMT_Parameters *par, const int phi) { + ATH_MSG_DEBUG("createRoads_fillHits: Feeding hitDatas Start"); + + diamond_t entry; + entry.wedgeCounter = iterator; + entry.sector = par->getSector(); + entry.phi = phi; + + micromegas_t micromegas; + MMDetectorHelper aHelper; + MMDetectorDescription* mm = aHelper.Get_MMDetector(par->getSector(), 1, 5, 1, 'A'); + MMReadoutParameters roP = mm->GetReadoutParameters(); + MMDetectorDescription* mm2 = aHelper.Get_MMDetector(par->getSector(), 2, 5, 1, 'A'); + MMReadoutParameters roP_2 = mm2->GetReadoutParameters(); + + micromegas.roadSize = this->getRoadSize(); + micromegas.nstrip_up_XX = this->getRoadSizeUpX(); + micromegas.nstrip_dn_XX = this->getRoadSizeDownX(); + micromegas.nstrip_up_UV = this->getRoadSizeUpUV(); + micromegas.nstrip_dn_UV = this->getRoadSizeDownUV(); + micromegas.strips = roP.tStrips; + micromegas.pitch = roP.stripPitch; + micromegas.nMissedTopEta = roP.nMissedTopEta; + micromegas.nMissedBottomEta = roP.nMissedBottomEta; + micromegas.nMissedTopStereo = roP.nMissedTopStereo; + micromegas.nMissedBottomStereo = roP.nMissedBottomStereo; + + micromegas.dimensions_top = mm->lWidth(); + micromegas.dimensions_bottom = mm->sWidth(); + micromegas.dimensions_height = mm->Length(); + + micromegas.activeArea_top = roP.activeTopLength; + micromegas.activeArea_bottom = roP.activeBottomLength; + micromegas.activeArea_height = roP.activeH; + micromegas.innerRadiusEta1 = roP.distanceFromZAxis; + micromegas.innerRadiusEta2 = roP_2.distanceFromZAxis; + micromegas.stereoAngles = roP.stereoAngle; + + std::string sector = (par->getSector() == 'L') ? "MML" : "MMS"; + /* + * The following for-loop merges all plane global coordinates in one single shot: + * X & Y are constant in all layers, for a given phi (Not used at the moment) + * Z (layer coordinates) changes only for Small and Large sectors --> USED and working! + */ + for (unsigned int phi = 0; phi < 8; phi++) { + Amg::Vector3D globalPos(0.0, 0.0, 0.0); + int multilayer = (phi < 4) ? 1 : 2; + int gasgap = ((phi+1)%4 == 0) ? 4 : (phi+1)%4; + /* + * Strip 100 (or 200) chosen as good compromise, considering offline strips. + * channelMin() function (https://acode-browser1.usatlas.bnl.gov/lxr/source/athena/MuonSpectrometer/MuonIdHelpers/MuonIdHelpers/MmIdHelper.h?v=21.3#0123) + * could be used, instead of hardcoding the strip id, but it returns (always?) as initialized in level ranges + */ + int strip = (phi > 1 && phi < 6) ? roP.nMissedBottomStereo+1 : roP.nMissedBottomEta+1; + Identifier strip_id = detManager->mmIdHelper()->channelID(sector, 1, phi+1, multilayer, gasgap, strip); + const MuonGM::MMReadoutElement* readout = detManager->getMMReadoutElement(strip_id); + + ROOT::Math::XYZVector coord(0.,0.,0.); + if (readout->stripGlobalPosition(strip_id, globalPos)) coord.SetXYZ(globalPos.x(), globalPos.y(), globalPos.z()); + else ATH_MSG_WARNING("Wedge " << sector << " phi: " << phi << " mult. " << multilayer << " gas " << gasgap << " | Unable to retrieve global positions"); + micromegas.planeCoordinates.push_back(coord); + } + + int nroad = 8192/this->getRoadSize(); + double B = (1./TMath::Tan(1.5/180.*TMath::Pi())); + int uvfactor = std::round( mm2->lWidth() / (B * 0.4 * 2.)/this->getRoadSize() ); // mm2 pointer is used because the full wedge has to be considered, i.e. S(L/M)2 + this->setUVfactor(uvfactor); + + for (int ihds = 0; ihds < (int)hitDatas.size(); ihds++) { + MMT_Hit *myhit = new MMT_Hit(par->getSector(), hitDatas[ihds], detManager); + myhit->updateHitProperties(par); + if (myhit->verifyHit()) { + m_hitslopes.push_back(myhit->getRZSlope()); + entry.ev_hits.push_back(myhit); + } else delete myhit; + } + + std::vector<MMT_Road*> temp_roads; + for (int i = 0; i < nroad; i++) { + + MMT_Road* myroad = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i); + temp_roads.push_back(myroad); + + int nuv = (this->getUV()) ? this->getUVfactor() : 0; + for (int uv = 1; uv <= nuv; uv++) { + if (i-uv < 0) continue; + + MMT_Road* myroad_0 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i+uv, i-uv); + temp_roads.push_back(myroad_0); + + MMT_Road* myroad_1 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i-uv, i+uv); + temp_roads.push_back(myroad_1); + + MMT_Road* myroad_2 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i+uv-1, i-uv); + temp_roads.push_back(myroad_2); + + MMT_Road* myroad_3 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i-uv, i+uv-1); + temp_roads.push_back(myroad_3); + + MMT_Road* myroad_4 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i-uv+1, i+uv); + temp_roads.push_back(myroad_4); + + MMT_Road* myroad_5 = new MMT_Road(par->getSector(), detManager, micromegas, this->getXthreshold(), this->getUVthreshold(), i, i+uv, i-uv+1); + temp_roads.push_back(myroad_5); + } + } + + for (const auto &road : temp_roads) { + if (this->isTrapezoidalShape()) { + if (road->iRoadu() < 0. && road->iRoadv() < 0.) continue; + /* + * Here some condition(s) about the trapezoidal shape should be added + */ + entry.ev_roads.push_back(road); + } else entry.ev_roads.push_back(road); + } + + temp_roads.clear(); + m_diamonds.push_back(entry); + + ATH_MSG_DEBUG("CreateRoadsAndFillHits: Feeding hitDatas Ended"); +} + +void MMT_Diamond::printHits(const unsigned int iterator) { + if (iterator < m_diamonds.size() && !m_diamonds[iterator].ev_hits.empty()) { + for (const auto &hit : m_diamonds[iterator].ev_hits) hit->printHit(); + } else { + ATH_MSG_DEBUG("Hit vector is empty!"); + } +} + +void MMT_Diamond::findDiamonds(const unsigned int iterator, const double &sm_bc, const int &event) { + auto t0 = std::chrono::high_resolution_clock::now(); + int ntrig = 0; + int bc_start = 999999; + int bc_end = -1; + int bc_wind = 4; // fixed time window (in bunch crossings) during which the algorithm collects ART hits + unsigned int ibc = 0; + + m_diamonds[iterator].slopes.clear(); + + // Comparison with lambda function (easier to implement) + std::sort(m_diamonds[iterator].ev_hits.begin(), m_diamonds[iterator].ev_hits.end(), [](MMT_Hit *h1, MMT_Hit *h2){ return h1->getAge() < h2->getAge(); }); + bc_start = m_diamonds[iterator].ev_hits.front()->getAge() - bc_wind*2; + bc_end = m_diamonds[iterator].ev_hits.back()->getAge() + bc_wind*2; + ATH_MSG_DEBUG("Window Start: " << bc_start << " - Window End: " << bc_end); + + for (const auto &road : m_diamonds[iterator].ev_roads) road->reset(); + + std::vector<MMT_Hit*> hits_now = {}; + std::vector<int> vmm_same = {}; + std::vector< std::pair<int, int> > addc_same = {}; + std::vector<int> to_erase = {}; + int n_vmm = std::ceil(8192/64.); + int n_addc = std::ceil(8192/2048.); + + // each road makes independent triggers, evaluated on each BC + for (int bc = m_diamonds[iterator].ev_hits.front()->getBC(); bc < bc_end; bc++) { + // Cleaning stuff + hits_now.clear(); + + for (unsigned int j = ibc; j < m_diamonds[iterator].ev_hits.size(); j++) { + if (m_diamonds[iterator].ev_hits[j]->getAge() == bc) hits_now.push_back(m_diamonds[iterator].ev_hits[j]); + else if (m_diamonds[iterator].ev_hits[j]->getAge() > bc) { + ibc = j; + break; + } + } + + // Implement VMM Art signal + for (unsigned int ib = 1; ib < 9; ib++) { + for (int j = 1; j < n_vmm; j++) { + vmm_same.clear(); + for (unsigned int k = 0; k < hits_now.size(); k++) { + if ((unsigned int)(hits_now[k]->getPlane()) != ib) continue; + if (hits_now[k]->getVMM() == j) vmm_same.push_back(k); + } + /* + * Is this really necessary with signal only? It kills most of the hits + * + * if (vmm_same.size() > 1) { + * int the_chosen_one = ran->Integer((int)(vmm_same.size())); + * for (int k = vmm_same.size()-1; k > -1; k--) if (k != the_chosen_one) hits_now.erase(hits_now.begin()+vmm_same[k]); + * } + */ + } + + // Implement ADDC-like handling + for (int ia = 1; ia < n_addc; ia++) { // From 1 to 4 + addc_same.clear(); + for (unsigned int k = 0; k < hits_now.size(); k++) { + if ((unsigned int)(hits_now[k]->getPlane()) == ib && hits_now[k]->getART() == ia) addc_same.push_back( std::make_pair(k, hits_now[k]->getChannel()) ); + } + + if (addc_same.size() <= 8) continue; + else std::cout << "Going further in ADDC loop" << std::endl; + + // priority encode the hits by channel number; remember hits 8+ + to_erase.clear(); + + std::sort(addc_same.begin(), addc_same.end(), [](std::pair<int, int> p1, std::pair<int, int> p2) { return p1.second < p2.second; }); + for (unsigned int it = 8; it < addc_same.size(); it++) to_erase.push_back(addc_same[it].first); + + // reverse and erase + std::sort(to_erase.rbegin(), to_erase.rend()); + for (auto k : to_erase) { + if (hits_now[k]->isNoise() == false) continue; + else hits_now.erase(hits_now.begin() + k); + } + } + } + + for (auto &road : m_diamonds[iterator].ev_roads) { + road->incrementAge(bc_wind); + if (!hits_now.empty()) road->addHits(hits_now); + + if (road->checkCoincidences(bc_wind) && bc >= (sm_bc - 1)) { + ntrig++; + + ATH_MSG_DEBUG("------------------------------------------------------------------"); + ATH_MSG_DEBUG("Coincidence FOUND @BC: " << bc); + ATH_MSG_DEBUG("Road (i, u, v, count): (" << road->iRoad() << ", " << road->iRoadu() << ", " << road->iRoadv() << ", " << road->countHits() << ")"); + ATH_MSG_DEBUG("------------------------------------------------------------------"); + for (const auto &hit: road->getHitVector()) { + ATH_MSG_DEBUG("Hit (board, BC, strip, eta): (" << hit.getPlane() << ", " << hit.getBC() << ", " << hit.getChannel() << ", " << hit.getStationEta() << ")"); + } + + slope_t slope; + slope.event = event; + slope.BC = bc; + slope.totalCount = road->countHits(); + slope.realCount = road->countRealHits(); + slope.iRoad = road->iRoad(); + slope.iRoadu = road->iRoadu(); + slope.iRoadv = road->iRoadv(); + slope.uvbkg = road->countUVHits(true); // the bool in the following 4 functions refers to background/noise hits + slope.xbkg = road->countXHits(true); + slope.uvmuon = road->countUVHits(false); + slope.xmuon = road->countXHits(false); + slope.age = bc - sm_bc; + slope.mxl = road->mxl(); + slope.xavg = road->avgSofX(); // defined as my in ATL-COM-UPGRADE-2015-033 + slope.uavg = road->avgSofUV(2,4); + slope.vavg = road->avgSofUV(3,5); + double mx = (slope.uavg-slope.vavg)/(2.*TMath::Tan(0.02618)); // The stereo angle is fixed and can be hardcoded + double theta = TMath::ATan(TMath::Sqrt(TMath::Power(mx,2) + TMath::Power(slope.xavg,2))); + slope.theta = (slope.xavg > 0.) ? theta : TMath::Pi() - theta; + slope.eta = -1.*TMath::Log(TMath::Tan(slope.theta/2.)); + slope.dtheta = (slope.mxl - slope.xavg)/(1. + slope.mxl*slope.xavg); + int wedge = (this->getDiamond(iterator).sector == 'L') ? 0 : 1; + int n = 2*(this->getDiamond(iterator).phi -1) + wedge; + char side = (slope.xavg > 0.) ? 'A' : 'C'; + double phi = TMath::ATan(mx/slope.xavg), phiShifted = this->phiShift(n, phi, side); + slope.phi = phi; + slope.phiShf = phiShifted; + + m_diamonds[iterator].slopes.push_back(slope); + } + } + } + auto t1 = std::chrono::high_resolution_clock::now(); + ATH_MSG_DEBUG("Processing roads took " << std::chrono::duration_cast<std::chrono::milliseconds>(t1 - t0).count() << " ms"); +} + +double MMT_Diamond::phiShift(const int &n, const double &phi, const char &side) { + double Phi = (side == 'A') ? phi : -phi; + float shift = (n > 8) ? (16-n)*TMath::Pi()/8. : n*TMath::Pi()/8.; + if (n < 8) return (Phi + shift); + else if (n == 8) return (Phi + ((Phi > 0.) ? -1. : 1.)*shift); + else return (Phi - shift); +} + +void MMT_Diamond::resetSlopes() { + if (!m_hitslopes.empty()) m_hitslopes.clear(); +} + +slope_t::slope_t(int ev, int bc, unsigned int tC, unsigned int rC, int iX, int iU, int iV, unsigned int uvb, unsigned int xb, unsigned int uvm, unsigned int xm, + int age, double mxl, double xavg, double uavg, double vavg, double th, double eta, double dth, double phi, double phiS) : + event(ev), BC(bc), totalCount(tC), realCount(rC), iRoad(iX), iRoadu(iU), iRoadv(iV), uvbkg(uvb), xbkg(xb), uvmuon(uvm), xmuon(xm), + age(age), mxl(mxl), xavg(xavg), uavg(uavg), vavg(vavg), theta(th), eta(eta), dtheta(dth), phi(phi), phiShf(phiS) {} diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Finder.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Finder.cxx index d923ed329bd7f0a45ed0e8d056f9f4c65d7872bf..e91416e5237b4805ce51709e36e28b1161bccfb5 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Finder.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Finder.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "AthenaKernel/getMessageSvc.h" @@ -18,7 +18,7 @@ MMT_Finder::MMT_Finder(MMT_Parameters *par, int nUVRoads) : m_par = par; m_nUVRoads = nUVRoads; - m_nRoads = ceil( ( ( m_par->slope_max - m_par->slope_min ) / m_par->h.getFixed() ).getFixed() ); //initialization, can use floats + m_nRoads = ceil( ( ( m_par->slope_max - m_par->slope_min ) / m_par->h ) ); //initialization, can use floats if(m_nUVRoads>1){ m_nRoads *= m_nUVRoads; // This should probably be configurable and dynamic based on the geometry of the chamber @@ -26,7 +26,7 @@ MMT_Finder::MMT_Finder(MMT_Parameters *par, int nUVRoads) : int nplanes=m_par->setup.size(); - ATH_MSG_DEBUG( "MMT_Find::finder entries " << m_nRoads << " " << m_par->slope_max.getFixed() << " " << m_par->slope_min.getFixed() << " " << m_par->h.getFixed() ); + ATH_MSG_DEBUG( "MMT_Find::finder entries " << m_nRoads << " " << m_par->slope_max << " " << m_par->slope_min << " " << m_par->h << " " << nplanes ); m_gateFlags = vector<vector<double> >(m_nRoads,(vector<double>(2,0)));// sloperoad, m_finder = vector<vector<finder_entry> >(m_nRoads,(vector<finder_entry>(nplanes,finder_entry()))); //[strip,slope,hit_index]; @@ -42,28 +42,31 @@ void MMT_Finder::fillHitBuffer( map< pair<int,int> , finder_entry > & hitBuffer, // This function takes in the Hit object and places it into the hit buffer hitBuffer, putting it in any relevant (road,plane) //Get initial parameters: tolerance, step size (h), slope of hit - float32fixed<3> tol; - float32fixed<3> h=m_par->h.getFixed(); + double tol; + double h=m_par->h; //Conver hit to slope here - float32fixed<3> slope=hit.info.slope.getFixed(); - ATH_MSG_DEBUG("SLOPE " << hit.info.slope.getFixed()); + double slope=hit.info.slope; + ATH_MSG_DEBUG("SLOPE " << hit.info.slope); //Plane and key info of the hit int plane=hit.info.plane; string plane_type=m_par->setup.substr(plane,1); - if(plane_type=="x") tol=m_par->x_error.getFixed(); - else if(plane_type=="u"||plane_type=="v") tol=m_par->uv_error.getFixed(); - else return; //if it's an unsupported plane option, don't fill + if(plane_type=="x") tol=m_par->x_error; + else if(plane_type=="u"||plane_type=="v") tol=m_par->uv_error; + else { + ATH_MSG_WARNING("WARNING: unsupported plane option!"); + return; + } //if it's an unsupported plane option, don't fill //---slope road boundaries based on hit_slope +/- tolerance---; if min or max is out of bounds, put it at the limit - float32fixed<3> s_min = slope - tol, s_max = slope + tol; + double s_min = slope - tol, s_max = slope + tol; - int road_min = round( ( (s_min - m_par->slope_min)/h ).getFixed() ); - int road_max = round( ( (s_max - m_par->slope_min)/h ).getFixed() ); + int road_min = round( ( (s_min - m_par->slope_min)/h ) ); + int road_max = round( ( (s_max - m_par->slope_min)/h ) ); if( road_min < 0 ) road_min = 0 ; if( road_max >= (m_nRoads/m_nUVRoads) ){ road_max = (m_nRoads/m_nUVRoads) - 1 ; } diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Fitter.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Fitter.cxx index f4ef8b40b834e5b1b13d72d1f9d28cad71f8ff19..d2dd972d7a464b961cbd2e3013838028adc65281 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Fitter.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Fitter.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "AthenaKernel/getMessageSvc.h" @@ -34,79 +34,86 @@ evFit_entry MMT_Fitter::fit_event(int event, vector<Hit>& track, vector<hitData_ vector<int> upl=m_par->q_planes("u"); vector<int> vpl=m_par->q_planes("v"); //---- Calc global slopes and local X slope ----- - float32fixed<2> M_x_global = Get_Global_Slope(track,"x"); - float32fixed<2> M_u_global = (check>=10?float32fixed<2>(-999.):Get_Global_Slope(track,"u")); - float32fixed<2> M_v_global = (check%10==1?float32fixed<2>(-999.):Get_Global_Slope(track,"v")); + double M_x_global = Get_Global_Slope(track,"x"); + double M_u_global = (check>=10?(-999.):Get_Global_Slope(track,"u")); + double M_v_global = (check%10==1?(-999.):Get_Global_Slope(track,"v")); //---- Calc delta theta ---------- //----- Calc ROI ---------- ROI ROI = Get_ROI(M_x_global,M_u_global,M_v_global,track); - mxmy.push_back(pair<double,double>(ROI.m_x.getFixed(),ROI.m_y.getFixed())); - mu = M_u_global.getFixed(); - mv = M_v_global.getFixed(); + mxmy.push_back(pair<double,double>(ROI.m_x,ROI.m_y)); + mu = M_u_global; + mv = M_v_global; if(ROI.theta==-999){ for(unsigned int i=0;i<track.size();i++) track[i].print(); - ATH_MSG_WARNING("SOMETHING IS OFF! fit_event\n"); - throw std::runtime_error("MMT_Fitter::fit_event: invalid ROI.theta"); - } - static_assert(std::is_trivially_copyable<float32fixed<2>>::value); - static_assert(std::is_trivially_destructible<float32fixed<2>>::value); - float32fixed<2> M_x_local = Get_Local_Slope(track,ROI.theta.getFixed(),ROI.phi.getFixed()),Delta_Theta_division = Get_Delta_Theta_division(M_x_local,M_x_global,1.), Delta_Theta = Get_Delta_Theta(M_x_local,M_x_global), dtheta_idl=Get_Delta_Theta_division(ideal_local_slope(track),M_x_global); - - mxl = M_x_local.getFixed(); - - if(dtheta_idl-Delta_Theta_division.fabs()>2.e-3)m_par->fill0=true; - ATH_MSG_DEBUG("Mxg="<<M_x_global.getFixed()<<",Mug="<<M_u_global.getFixed()<<",Mvg="<<M_v_global.getFixed()<<",Mxl="<<M_x_local.getFixed()<<",dth="<<Delta_Theta.getFixed()); + ATH_MSG_WARNING("SOMETHING IS OFF! fit_event, ROI has theta=-999 (bad reconstruction) \n"); + evFit_entry empty(event,-999,-999,-999,ROI.roi,-999,-999,-999,-999,-999,track_to_index(track)); + return empty; + //throw std::runtime_error("MMT_Fitter::fit_event: invalid ROI.theta"); + } + static_assert(std::is_trivially_copyable<double>::value); + static_assert(std::is_trivially_destructible<double>::value); + double M_x_local = Get_Local_Slope(track, ROI.theta, ROI.phi); + double Delta_Theta_division = Get_Delta_Theta_division(M_x_local, M_x_global, 1.); + double Delta_Theta = Get_Delta_Theta(M_x_local, M_x_global); + double dtheta_idl = Get_Delta_Theta_division(ideal_local_slope(track), M_x_global); + + mxl = M_x_local; + + if(dtheta_idl-std::abs(Delta_Theta_division) > 2.e-3) m_par->fill0 = true; + ATH_MSG_DEBUG("Mxg=" << M_x_global << ", Mug=" << M_u_global << ", Mvg=" << M_v_global << ", Mxl=" << M_x_local << ", dth=" <<Delta_Theta); //@@@@@@@@ Begin Info Storage for Later Analysis @@@@@@@@@@@@@@@@@@@@@@@ - vector<bool> planes_hit_tr(8,false),planes_hit_bg(8,false); + vector<bool> planes_hit_tr(8,false), planes_hit_bg(8,false); for(unsigned int ihit=0; ihit<track.size(); ihit++){ - int hitData_pos=find_hitData(hitDatas,track[ihit].key); - if(hitData_pos==-1) continue; + int hitData_pos = find_hitData(hitDatas,track[ihit].key); + if(hitData_pos == -1) continue; if(hitDatas[hitData_pos].truth_nbg) planes_hit_tr[track[ihit].info.plane]=true; else planes_hit_bg[track[ihit].info.plane]=true; } - int n_xpl_tr=0,n_xpl_bg=0,n_uvpl_tr=0,n_uvpl_bg=0; + int n_xpl_tr=0, n_xpl_bg=0, n_uvpl_tr=0, n_uvpl_bg=0; for(int xp=0; xp<(int)xpl.size(); xp++){ - int plane=xpl[xp]; - n_xpl_tr+=planes_hit_tr[plane];n_xpl_bg+=planes_hit_bg[plane]; + int plane = xpl[xp]; + n_xpl_tr += planes_hit_tr[plane]; + n_xpl_bg += planes_hit_bg[plane]; } if(check<10){ for(unsigned int up=0; up<upl.size(); up++){ - int plane=upl[up]; - n_uvpl_tr+=planes_hit_tr[plane]; - n_uvpl_bg+=planes_hit_bg[plane]; + int plane = upl[up]; + n_uvpl_tr += planes_hit_tr[plane]; + n_uvpl_bg += planes_hit_bg[plane]; } } if(check%10==0){ for(unsigned int vp=0; vp<vpl.size(); vp++){ - int plane=vpl[vp]; - n_uvpl_tr+=planes_hit_tr[plane]; - n_uvpl_bg+=planes_hit_bg[plane]; + int plane = vpl[vp]; + n_uvpl_tr += planes_hit_tr[plane]; + n_uvpl_bg += planes_hit_bg[plane]; } } //FINISH ME!!!!!!! - float32fixed<4> candtheta=ROI.theta,candphi=ROI.phi; + double candtheta=ROI.theta,candphi=ROI.phi; /* I think this bit appears redundant but could end up being stupid; the fitter shouldn't care about CT stuff (beyond having min num hits to do fit, which the -999 or w/e is responsible for taking care of) bool xfail=(n_xpl_tr+n_xpl_bg<m_par->CT_x),uvfail= (n_uvpl_tr+n_uvpl_bg<m_par->CT_uv); msg(MSG::DEBUG) << n_xpl_tr+n_xpl_bg<<" x hits"); if(xfail)candtheta=-999.; if(uvfail)candphi=-999.; */ - bool fitkill=(ROI.theta==-999 || Delta_Theta==-999||Delta_Theta==-4);// ||xfail||uvfail); - ATH_MSG_DEBUG("HIT CODE: "<<track_to_index(track)); + bool fitkill = (ROI.theta==-999 || Delta_Theta==-999||Delta_Theta==-4);// ||xfail||uvfail); + ATH_MSG_DEBUG("HIT CODE: " << track_to_index(track)); evFit_entry aemon(event,candtheta,candphi,Delta_Theta_division,ROI.roi,n_xpl_tr,n_uvpl_tr,n_xpl_bg,n_uvpl_bg,Delta_Theta,track_to_index(track)); if(fitkill) return aemon; - int nplanes=m_par->setup.size(); + + int nplanes = m_par->setup.size(); for(int plane=0; plane<nplanes; plane++){ - if(track[plane].info.slope==-999) continue; //&& Delta_Theta_division~=-999 - int hitData_pos=find_hitData(hitDatas,track[plane].key); - if(hitData_pos==-1) continue; - did_fit=true; + if(track[plane].info.slope == -999) continue; //&& Delta_Theta_division~=-999 + int hitData_pos = find_hitData(hitDatas, track[plane].key); + if(hitData_pos == -1) continue; + did_fit = true; hitDatas[hitData_pos].fit_fill(ROI.theta,ROI.phi,Delta_Theta,M_x_global,M_u_global,M_v_global,M_x_local,ROI.m_x,ROI.m_y,ROI.roi); aemon.fit_hit_keys.push_back(track[plane].key); - if(hitDatas[hitData_pos].truth_nbg) aemon.truth_planes_hit+=pow(10,nplanes-plane-1); - else aemon.bg_planes_hit+=pow(10,nplanes-plane-1); + if(hitDatas[hitData_pos].truth_nbg) aemon.truth_planes_hit += std::pow(10,nplanes-plane-1); + else aemon.bg_planes_hit += std::pow(10,nplanes-plane-1); } if(did_fit) nfit++; return aemon; @@ -114,7 +121,7 @@ evFit_entry MMT_Fitter::fit_event(int event, vector<Hit>& track, vector<hitData_ int MMT_Fitter::find_hitData(const vector<hitData_entry>& hitDatas, const hitData_key& key) const{ for(unsigned int i=0;i<hitDatas.size();i++){ - if(hitDatas[i].BC_time==key.BC_time&&hitDatas[i].time==key.time&&hitDatas[i].gtime==key.gtime&&hitDatas[i].VMM_chip==key.VMM_chip) return i; + if(hitDatas[i].BC_time==key.BC_time && hitDatas[i].time==key.time && hitDatas[i].gtime==key.gtime && hitDatas[i].VMM_chip==key.VMM_chip) return i; } return -1; } @@ -122,35 +129,36 @@ int MMT_Fitter::find_hitData(const vector<hitData_entry>& hitDatas, const hitDat int MMT_Fitter::Filter_UV(vector<Hit>& track) const{ return 0; - float32fixed<2> h=m_par->h, tolerance = h;//*2; //Can be optimized... - vector<int> u_planes=m_par->q_planes("u"),v_planes=m_par->q_planes("v"); - vector<Hit> u_hits=q_hits("u",track),v_hits=q_hits("v",track); - bool pass_u=!u_hits.empty(),pass_v=!v_hits.empty(); + double h=m_par->h, tolerance = h;//*2; //Can be optimized... + vector<int> u_planes = m_par->q_planes("u"), v_planes = m_par->q_planes("v"); + vector<Hit> u_hits = q_hits("u",track), v_hits = q_hits("v",track); + bool pass_u =! u_hits.empty(), pass_v =! v_hits.empty(); + //if the difference in slope between the first and last u/v planes is too great don't pass, set track hits to zero if(pass_u){ - if((float32fixed<2>(u_hits.front().info.slope-u_hits.back().info.slope)).fabs()>tolerance) pass_u=false; + if(std::abs( double(u_hits.front().info.slope-u_hits.back().info.slope)) > tolerance) pass_u=false; } if(pass_v){ - if((float32fixed<2>(v_hits.front().info.slope-v_hits.back().info.slope)).fabs()>tolerance) pass_v=false; + if(std::abs( double(v_hits.front().info.slope-v_hits.back().info.slope)) > tolerance) pass_v=false; } - int return_val=0;//of form (bool ubad, bool vbad), so 10 would be bad u, good v + int return_val = 0; //of form (bool ubad, bool vbad), so 10 would be bad u, good v if(!pass_u){ - for(unsigned int iup=0;iup<u_planes.size();iup++) track[u_planes[iup]].info.slope=-999; - return_val+=10; + for(unsigned int iup=0; iup<u_planes.size(); iup++) track[u_planes[iup]].info.slope = -999; + return_val += 10; } if(!pass_v){ - for(unsigned int ivp=0;ivp<v_planes.size();ivp++) track[v_planes[ivp]].info.slope=-999; - return_val+=1; + for(unsigned int ivp=0; ivp<v_planes.size(); ivp++) track[v_planes[ivp]].info.slope = -999; + return_val += 1; } return return_val; } -float32fixed<2> MMT_Fitter::Get_Global_Slope (const vector<Hit>& track, const string& type) const +double MMT_Fitter::Get_Global_Slope (const vector<Hit>& track, const string& type) const { vector<Hit> qhits = q_hits(type,track); - float32fixed<2> sum = 0.; + double sum = 0.; if(qhits.size()==0) return -999; - float32fixed<2> nhitdiv= 1./qhits.size(); + double nhitdiv= 1./qhits.size(); for(int ihit=0;ihit<(int)qhits.size();ihit++){ sum+=(qhits[ihit].info.slope*nhitdiv); } @@ -158,44 +166,55 @@ float32fixed<2> MMT_Fitter::Get_Global_Slope (const vector<Hit>& track, const st } //CHANGE! -float32fixed<2> MMT_Fitter::Get_Local_Slope (const vector<Hit>& Track,double theta,double phi) const +double MMT_Fitter::Get_Local_Slope (const vector<Hit>& Track,double theta,double phi) const { - vector<int> x_planes=m_par->q_planes("x"),ybin_hits(x_planes.size(),-1); + vector<int> x_planes = m_par->q_planes("x"), ybin_hits(x_planes.size(),-1); int nxp=x_planes.size(); for(int ipl=0; ipl<nxp; ipl++){ - ybin_hits[ipl]=((Track[x_planes[ipl]].info.slope==-999||Track[x_planes[ipl]].info.slope==-4)?-1:m_par->ybin(Track[x_planes[ipl]].info.y)); + ybin_hits[ipl]=(((Track[x_planes[ipl]].info.slope==-999) || (Track[x_planes[ipl]].info.slope==-4)) ? -1 : m_par->ybin(Track[x_planes[ipl]].info.y,x_planes[ipl])); } bool hit=false; - float32fixed<yzdex> yzsum=0; + double yzsum=0; + double ysum=0;//added to try alternative calculus + double zsum=0;//added to try alternative calculus float mxlf=0; int xdex=-1; int ybin=-1; int which=-1; - m_par->key_to_indices(ybin_hits,xdex,ybin,which); - if(xdex<0||ybin<0||which<0) return -999; - float32fixed<zbardex> zbar=m_par->Ak_local_slim[xdex][ybin][which]; - float32fixed<bkdex>bk=m_par->Bk_local_slim[xdex][ybin][which]; - ATH_MSG_DEBUG("zbar is "<<zbar.getFixed()<<", and bk is "<<bk.getFixed()); - int ebin=m_par->eta_bin(theta),pbin=m_par->phi_bin(phi); + m_par->key_to_indices(ybin_hits, xdex, ybin, which); + if(xdex<0 || ybin<0 || which<0) return -999; + + double zbar=m_par->Ak_local_slim[xdex][ybin][which]; + double bk=m_par->Bk_local_slim[xdex][ybin][which]; + ATH_MSG_DEBUG("zbar is " << zbar << ", and bk is " << bk); + int ebin = m_par->eta_bin(theta); + int pbin = m_par->phi_bin(phi); for(int ipl=0; ipl<nxp; ipl++){ - float32fixed<yzdex> z=Track[x_planes[ipl]].info.z,y=Track[x_planes[ipl]].info.y; - if(ebin!=-1){ - z=z+m_par->zmod[ebin][pbin][ipl]; - y=y+m_par->ymod[ebin][pbin][ipl]; + double z = Track[x_planes[ipl]].info.z; + double y = Track[x_planes[ipl]].info.y; + double zfirst = Track[x_planes[0]].info.z; + double yfirst = Track[x_planes[0]].info.y; + if(ebin != -1){ + z = z+m_par->zmod[ebin][pbin][ipl]; + y = y+m_par->ymod[ebin][pbin][ipl]; } - if(Track[x_planes[ipl]].info.slope==-999||Track[x_planes[ipl]].info.slope==-4) continue; + if(Track[x_planes[ipl]].info.slope == -999 || Track[x_planes[ipl]].info.slope == -4) continue; hit=true; yzsum+=y*(z*zbar-1.); - mxlf += bk.getFixed()*y.getFixed()*(z.getFixed()*zbar.getFixed()-1.); - } - float32fixed<2> mxl=float32fixed<2>(bk.getFixed()*yzsum.getFixed()); - if(!hit) {return float32fixed<2>(999);} + mxlf += bk*y*(z*zbar-1.); + //alternative calculus + ysum += y-yfirst; + zsum += z-zfirst; + } + //double mxl=double(bk*yzsum); //old calculus + double mxl=(ysum/zsum); //new calculus + if(!hit) {return double(999);} return mxl; } int MMT_Fitter::track_to_index(const vector<Hit>&track)const{ vector<bool>hits(m_par->setup.size(),false); - for(int ihit=0;ihit<(int)track.size();ihit++)hits[track[ihit].info.plane]=(hits[track[ihit].info.plane]?true:track[ihit].info.slope>-2.); + for(int ihit=0; ihit<(int)track.size(); ihit++) hits[track[ihit].info.plane] = (hits[track[ihit].info.plane] ? true : track[ihit].info.slope>-2.); return m_par->bool_to_index(hits); } @@ -203,17 +222,18 @@ double MMT_Fitter::ideal_local_slope(const vector<Hit>& Track)const{ vector<vector<double> > z_hit; for(unsigned int i = 0; i<m_par->z_large.size(); i++){ vector<double> temp; - for(unsigned int j = 0; j<m_par->z_large[i].size(); j++) - temp.push_back(m_par->z_large[i][j].getFixed()); + for(unsigned int j = 0; j<m_par->z_large[i].size(); j++) temp.push_back(m_par->z_large[i][j]); z_hit.push_back(temp); } vector<int> x_planes=m_par->q_planes("x"); int nxp=x_planes.size(); bool hit=false; - double sum_xy=0,sum_y=0; - double ak_idl=ideal_ak(Track),bk_idl=ak_idl*ideal_zbar(Track); + double sum_xy=0, sum_y=0; + double ak_idl = ideal_ak(Track); + double bk_idl = ak_idl*ideal_zbar(Track); for(int ipl=0; ipl<nxp; ipl++){ - double y=Track[x_planes[ipl]].info.y.getFixed(),z=Track[x_planes[ipl]].info.z.getFixed(); + double y = Track[x_planes[ipl]].info.y; + double z = Track[x_planes[ipl]].info.z; if(y==-999) continue; hit=true; sum_xy += ak_idl*z*y; @@ -226,8 +246,11 @@ double MMT_Fitter::ideal_local_slope(const vector<Hit>& Track)const{ double MMT_Fitter::ideal_z(const Hit& hit)const{ int plane=hit.info.plane; - double tilt=(plane<4?m_par->correct.rotate.X():0),dz=(plane<4?m_par->correct.translate.Z():0), - nominal=m_par->z_nominal[plane].getFixed(),y=hit.info.y.getFixed()-m_par->ybases[plane].front().getFixed(),z=nominal+dz+y*tan(tilt); + double tilt = (plane<4 ? m_par->correct.rotate.X() : 0); + double dz = (plane<4 ? m_par->correct.translate.Z() : 0); + double nominal = m_par->z_nominal[plane]; + double y = hit.info.y - m_par->ybases[plane].front(); + double z = nominal+dz+y*std::tan(tilt); return z; } @@ -241,25 +264,27 @@ double MMT_Fitter::ideal_ak(const vector<Hit>& Track)const{ hits++;//there's a hit sum_x += addme; sum_xx += addme*addme; - ATH_MSG_DEBUG("z["<<ip<<"]="<<addme<<", sum_x="<<sum_x<<", sum_xx="<<sum_xx<<", hits="<<hits); + ATH_MSG_DEBUG("z[" << ip << "]=" << addme << ", sum_x=" << sum_x << ", sum_xx=" << sum_xx << ", hits=" << hits); } double diff = hits*sum_xx-sum_x*sum_x; return hits/diff; } double MMT_Fitter::ideal_zbar(const vector<Hit>& Track)const{ - vector<int> x_planes=m_par->q_planes("x");//this tells us which planes are x planes - double ztot=0;int nhit=0; - for(unsigned int ip=0;ip<x_planes.size();ip++){ - if(Track[x_planes[ip]].info.slope==-999)continue; - nhit++;ztot+=ideal_z(Track[x_planes[ip]]); + vector<int> x_planes = m_par->q_planes("x");//this tells us which planes are x planes + double ztot=0; + int nhit=0; + for(unsigned int ip=0; ip<x_planes.size(); ip++) { + if(Track[x_planes[ip]].info.slope == -999) continue; + nhit++; + ztot += ideal_z(Track[x_planes[ip]]); } return ztot/nhit; } -float32fixed<2> MMT_Fitter::Get_Delta_Theta(float32fixed<2> M_local, float32fixed<2> M_global) const{ +double MMT_Fitter::Get_Delta_Theta(double M_local, double M_global) const{ int region=-1; - float32fixed<2> LG = M_local * M_global; + double LG = M_local * M_global; for(int j=0;j<m_number_LG_regions;j++){ //m_number_LG_regions if(LG <= DT_Factors_val(j,0)){ region = j; @@ -270,15 +295,15 @@ float32fixed<2> MMT_Fitter::Get_Delta_Theta(float32fixed<2> M_local, float32fixe return DT_Factors_val(region,1)*(M_local - M_global); } -float32fixed<2> MMT_Fitter::DT_Factors_val(int i, int j) const{ +double MMT_Fitter::DT_Factors_val(int i, int j) const{ if(m_par->val_tbl){ return m_par->DT_Factors[i][j]; } - if(j<0||j>1){ + if(j<0 || j>1){ ATH_MSG_WARNING("DT_Factors only has two entries on the second index (for LG and mult_factor); you inputed an index of " << j ); throw std::runtime_error("MMT_Fitter::DT_Factors_val: invalid index"); } - if(i<0||i>=m_number_LG_regions){ + if(i<0 || i>=m_number_LG_regions){ ATH_MSG_WARNING("There are " << m_number_LG_regions << " in DT_Factors(_val); you inputed an index of " << i ); throw std::runtime_error("MMT_Fitter::DT_Factors_val: invalid index"); } @@ -287,21 +312,21 @@ float32fixed<2> MMT_Fitter::DT_Factors_val(int i, int j) const{ return LG_lgr(i,a,m_number_LG_regions,m_LG_min,m_LG_max); } -float32fixed<2> MMT_Fitter::LG_lgr(int ilgr, double a, int number_LG_regions, float32fixed<2> min, float32fixed<2> max) const{ +double MMT_Fitter::LG_lgr(int ilgr, double a, int number_LG_regions, double min, double max) const{ a+=0; - return min+float32fixed<2>(ilgr/number_LG_regions)*(max-min); + return min+double(ilgr/number_LG_regions)*(max-min); } -float32fixed<2> MMT_Fitter::mult_factor_lgr(int ilgr, double a, int number_LG_regions, float32fixed<2> min, float32fixed<2> max) const{ - return float32fixed<2>(1.) / float32fixed<2>( ( LG_lgr(ilgr,a,number_LG_regions,min,max) / a + a ) ); +double MMT_Fitter::mult_factor_lgr(int ilgr, double a, int number_LG_regions, double min, double max) const{ + return double(1.) / double( ( LG_lgr(ilgr,a,number_LG_regions,min,max) / a + a ) ); } -float32fixed<2> MMT_Fitter::Get_Delta_Theta_division(float32fixed<2> M_local, float32fixed<2> M_global, float32fixed<4> a) const{ +double MMT_Fitter::Get_Delta_Theta_division(double M_local, double M_global, double a) const{ //we could use 2 bits for the numerator and 3 for the denominator, but then //fixed_point doesn't know how to do the algebra. Assume we know how to do //this division (we don't, efficiently, thus the method Get_Delta_Theta - return float32fixed<2>( (M_local - M_global).getFixed() / ( (M_local*M_global).getFixed() / a.getFixed() + a.getFixed() )); + return double( (M_local - M_global) / ( (M_local*M_global) / a + a )); } vector<Hit> MMT_Fitter::q_hits(const string& type,const vector<Hit>& track) const{ @@ -320,95 +345,106 @@ vector<Hit> MMT_Fitter::q_hits(const string& type,const vector<Hit>& track) cons } //change this to take u and/or v out of the roi calculation -ROI MMT_Fitter::Get_ROI(float32fixed<2> M_x,float32fixed<2> M_u,float32fixed<2> M_v,const vector<Hit>&track) const{ +ROI MMT_Fitter::Get_ROI(double M_x, double M_u, double M_v, const vector<Hit>&track) const{ //M_* are all global slopes - ATH_MSG_DEBUG("\nGet_ROI("<<M_x.getFixed()<<","<<M_u.getFixed()<<","<<M_v.getFixed()<<") "); + ATH_MSG_DEBUG("\nGet_ROI(" << M_x << "," << M_u << "," << M_v << ") "); //--- calc constants ------ - float32fixed<2> b=TMath::DegToRad()*(m_par->stereo_degree.getFixed()); - float32fixed<7> A=1./tan(b.getFixed()); - float32fixed<7> B=1./tan(b.getFixed()); + double b=TMath::DegToRad()*(m_par->stereo_degree); + double A=1./std::tan(b); + double B=1./std::tan(b); //--- slope conversion equations ---- - float32fixed<2> my = M_x; - float32fixed<2> mxu = ( A*M_u.getFixed() - B*my.getFixed() ).getFixed(); - float32fixed<2> mxv = ( B*my.getFixed() - A*M_v.getFixed() ).getFixed(); + double my = M_x; + double mxu = ( A*M_u - B*my ); + double mxv = ( B*my - A*M_v ); //--- which slopes are truly present ---- //Note that bad slopes are not necessarily 0 as I often use -999 to denote something missing //we have -999 for M_u or M_v to denote that it didn't pass filtering - int nu=1,nv=1; - if(M_u<0||M_u==float32fixed<2>(-999)){ - mxu = 0;nu=0; + int nu=1, nv=1; + if(M_u<0 || M_u==-999){ + mxu=0; + nu=0; } - if(M_v<0||M_v==float32fixed<2>(-999)){ - mxv=0;nv=0; + if(M_v<0 || M_v==-999){ + mxv=0; + nv=0; } - if(nu==0&&nv==0) return ROI(-999,-999,-999,-999,-999); + if(nu==0 && nv==0) return ROI(-999,-999,-999,-999,-999); //--- average of 2 mx slope values ----if both u and v were bad, give it a -999 value to know not to use m_x //*** check to see if U and V are necessary for fit - float32fixed<2> mx = (nu+nv==0?0:(mxv+mxu)/(nu+nv)); - if(m_par->correct.translate.X()!=0&&m_par->correct.type==2){ - mx+=phi_correct_factor(track)*m_par->correct.translate.X()/m_par->z_nominal[3].getFixed(); + double mx = (nu+nv==0 ? 0 : (mxv+mxu)/(nu+nv)); + if(m_par->correct.translate.X()!=0 && m_par->correct.type==2){ + mx += phi_correct_factor(track)*m_par->correct.translate.X()/m_par->z_nominal[3]; } - ATH_MSG_DEBUG("(b,A,B,my,mxu,mxv,mx)=("<<b.getFixed()<<","<<A.getFixed()<<","<<B.getFixed()<<","<<my.getFixed()<<","<<mxu.getFixed()<<","<<mxv.getFixed()<<","<<mx.getFixed()<<")\n"); + ATH_MSG_DEBUG("(b,A,B,my,mxu,mxv,mx)=(" << b << "," << A << "," << B << "," << my << "," << mxu << "," << mxv << "," << mx << ")\n"); //Get mx and my in parameterized values - int a_x = round((mx.getFixed()-m_par->m_x_min.getFixed())/m_par->h_mx.getFixed()), a_y = round((my.getFixed()-m_par->m_y_min.getFixed())/m_par->h_my.getFixed()); + int a_x = round((mx-m_par->m_x_min)/m_par->h_mx); + int a_y = round((my-m_par->m_y_min)/m_par->h_my); // Generally, this offers a reality check or cut. The only reason a slope // should be "out of bounds" is because it represents a weird UV combination // -- ie. highly background influenced if(a_y>m_par->n_y || a_y<0){ - ATH_MSG_DEBUG( "y slope (theta) out of bounds in Get_ROI....(a_x,a_y,m_par->n_x,m_par->n_y)=("<<a_x<<","<<a_y<<","<<m_par->n_x<<","<<m_par->n_y<<")"); + ATH_MSG_DEBUG( "y slope (theta) out of bounds in Get_ROI....(a_x,a_y,m_par->n_x,m_par->n_y)=("<< a_x << "," << a_y << "," << m_par->n_x << "," << m_par->n_y << ")"); return ROI(-999,-999,-999,-999,-999); } if(a_x>m_par->n_x || a_x<0){ - ATH_MSG_DEBUG( "x slope (phi) out of bounds in Get_ROI....(a_x,a_y,m_par->n_x,m_par->n_y)=("<<a_x<<","<<a_y<<","<<m_par->n_x<<","<<m_par->n_y<<")"); + ATH_MSG_DEBUG( "x slope (phi) out of bounds in Get_ROI....(a_x,a_y,m_par->n_x,m_par->n_y)=(" << a_x << "," << a_y << "," << m_par->n_x << "," <<m_par->n_y << ")"); return ROI(-999,-999,-999,-999,-999); } ATH_MSG_DEBUG("fv_angles...(a_x,a_y)=("<<a_x<<","<<a_y<<")"); double phicor=0.; - if(m_par->correct.rotate.Z()!=0&&m_par->correct.type==2){ + if(m_par->correct.rotate.Z()!=0 && m_par->correct.type==2){ phicor=-0.2*m_par->correct.rotate.Z(); } - float32fixed<4> fv_theta=Slope_Components_ROI_theta(a_y,a_x), fv_phi=(mx.getFixed()==0?-999:Slope_Components_ROI_phi(a_y,a_x).getFixed()+phicor); - ATH_MSG_DEBUG("fv_theta="<<fv_theta.getFixed()<<", fv_phi="<<fv_phi.getFixed()); + double fv_theta = Slope_Components_ROI_theta(a_y,a_x); + double fv_phi = (mx==0) ? -999 : Slope_Components_ROI_phi(a_y,a_x)+phicor; + ATH_MSG_DEBUG("fv_theta=" << fv_theta << ", fv_phi=" << fv_phi); //--- More hardware realistic approach but need fine tuning ---- int roi = Rough_ROI_temp(fv_theta,fv_phi); //--- current "roi" which is not an actual roi but an approx phi and theta - return ROI(fv_theta.getFixed(),fv_phi.getFixed(),mx.getFixed(),my.getFixed(),roi); + return ROI(fv_theta,fv_phi,mx,my,roi); } double MMT_Fitter::phi_correct_factor(const vector<Hit>&track)const{ - if((m_par->correct.rotate.Z()==0&&m_par->correct.translate.X()==0)||m_par->correct.type!=2)return 0.; - int nxmis=0,nx=0,numis=0,nu=0,nvmis=0,nv=0; - double xpart=0.5,upart=0.5,vpart=0.5; + if((m_par->correct.rotate.Z()==0 && m_par->correct.translate.X()==0) || m_par->correct.type!=2) return 0.; + int nxmis=0, nx=0, numis=0, nu=0, nvmis=0, nv=0; + double xpart=0.5, upart=0.5, vpart=0.5; string set=m_par->setup; for(int ihit=0;ihit<(int)track.size();ihit++){ int n_pln=track[ihit].info.plane; bool ismis=n_pln<4; char pln=set[n_pln]; - if(pln=='x'||pln=='X'){nx++;nxmis+=ismis;} - else if(pln=='u'||pln=='U'){nu++;numis+=ismis;} - else if(pln=='v'||pln=='V'){nv++;nvmis+=ismis;} + if(pln=='x'||pln=='X'){ + nx++; + nxmis+=ismis; + } else if(pln=='u'||pln=='U') { + nu++; + numis+=ismis; + } else if(pln=='v'||pln=='V') { + nv++; + nvmis+=ismis; + } } - if(nu==0&&nv==0)return 0.; - if(nu==0)upart=0.; + if(nu==0 && nv==0)return 0.; + if(nu==0) upart=0.; else if(nv==0)vpart=0.; else xpart=0.; return xpart*1.*nxmis/nx+upart*1.*numis/nu+vpart*1.*nvmis/nv; } -float32fixed<4> MMT_Fitter::Slope_Components_ROI_val(int jy, int ix, int thetaphi) const{ +double MMT_Fitter::Slope_Components_ROI_val(int jy, int ix, int thetaphi) const{ if(m_par->val_tbl){ return m_par->Slope_to_ROI[jy][ix][thetaphi]; } - if(thetaphi<0||thetaphi>1){ + if(thetaphi<0 || thetaphi>1){ ATH_MSG_WARNING("Slope_Components_ROI only has two entries on the third index (for theta and phi); you inputed an index of " << thetaphi); throw std::runtime_error("MMT_Fitter::Slope_Components_ROI_val: invalid number of entries"); } @@ -416,67 +452,72 @@ float32fixed<4> MMT_Fitter::Slope_Components_ROI_val(int jy, int ix, int thetaph return Slope_Components_ROI_phi(jy,ix); } -float32fixed<4> MMT_Fitter::Slope_Components_ROI_theta(int jy, int ix) const{ +double MMT_Fitter::Slope_Components_ROI_theta(int jy, int ix) const{ //get some parameter information - if(jy<0||jy>=m_par->n_y){ + if(jy<0 || jy>=m_par->n_y){ ATH_MSG_WARNING("You picked a y slope road index of " << jy << " in Slope_Components_ROI_theta; there are only " << m_par->n_y << " of these.\n"); - if(jy>=m_par->n_y)jy=m_par->n_y-1; + if(jy >= m_par->n_y) jy=m_par->n_y-1; else jy=0; } if(ix<0||ix>=m_par->n_x){ ATH_MSG_WARNING("You picked an x slope road index of " << ix << " in Slope_Components_ROI_theta; there are only " << m_par->n_x << " of these.\n"); - if(ix>=m_par->n_x)ix=m_par->n_x-1; + if(ix >= m_par->n_x) ix=m_par->n_x-1; else ix=0; } - int xdex=ix,ydex=jy+1; - if(xdex==0)xdex++; - float32fixed<2> mx=m_par->m_x_min+m_par->h_mx*xdex, my=m_par->m_y_min+m_par->h_my*ydex; - float32fixed<4> theta=atan(sqrt( (mx*mx+my*my).getFixed() )); + int xdex=ix, ydex=jy+1; + if(xdex==0) xdex++; + double mx = m_par->m_x_min+m_par->h_mx*xdex; + double my = m_par->m_y_min+m_par->h_my*ydex; + double theta = std::atan(sqrt( (mx*mx+my*my) )); if(theta<m_par->minimum_large_theta || theta>m_par->maximum_large_theta){ theta=0; } return theta; } -float32fixed<4> MMT_Fitter::Slope_Components_ROI_phi(int jy, int ix) const{ - if(jy<0||jy>=m_par->n_y){ +double MMT_Fitter::Slope_Components_ROI_phi(int jy, int ix) const{ + if(jy<0 || jy>=m_par->n_y){ ATH_MSG_WARNING("You picked a y slope road index of " << jy << " in Slope_Components_ROI_phi; there are only " << m_par->n_y << " of these.\n"); - if(jy>=m_par->n_y)jy=m_par->n_y-1; + if(jy >= m_par->n_y) jy=m_par->n_y-1; else jy=0; } - if(ix<0||ix>=m_par->n_x){ + if(ix<0 || ix>=m_par->n_x){ ATH_MSG_WARNING("You picked an x slope road index of " << ix << " in Slope_Components_ROI_phi; there are only " << m_par->n_x << " of these.\n"); //right now we're assuming these are cases just on the edges and so put the values to the okay limits - if(ix>=m_par->n_x)ix=m_par->n_x-1; + if(ix >= m_par->n_x) ix=m_par->n_x-1; else ix=0; } - int xdex=ix,ydex=jy+1; - float32fixed<2> mx=m_par->m_x_min+m_par->h_mx*xdex, my=m_par->m_y_min+m_par->h_my*ydex; - ATH_MSG_DEBUG( "m_par->m_x_min+m_par->h_mx*xdex="<<m_par->m_x_min.getFixed()<<"+"<<m_par->h_mx.getFixed()<<"*"<<xdex<<"="<<mx.getFixed()<<", "); - ATH_MSG_DEBUG( "m_par->m_y_min+m_par->h_my*ydex="<<m_par->m_y_min.getFixed()<<"+"<<m_par->h_my.getFixed()<<"*"<<ydex<<"="<<my.getFixed()<<", "); - float32fixed<4> phi(atan2(mx.getFixed(),my.getFixed()));//the definition is flipped from what you'd normally think - msg(MSG::DEBUG) << "for a phi of "<<phi.getFixed()<< endmsg; - if(phi<m_par->minimum_large_phi || phi>m_par->maximum_large_phi){ - msg(MSG::DEBUG) << "Chucking phi of " << phi.getFixed()<<" which registers as not in ["<<m_par->minimum_large_phi.getFixed()<<","<<m_par->maximum_large_phi.getFixed()<<"]"<< endmsg; + int xdex=ix, ydex=jy+1; + double mx=m_par->m_x_min+m_par->h_mx*xdex; + double my=m_par->m_y_min+m_par->h_my*ydex; + ATH_MSG_DEBUG("m_par->m_x_min+m_par->h_mx*xdex=" << m_par->m_x_min << "+" << m_par->h_mx << "*" << xdex << "=" << mx << ", "); + ATH_MSG_DEBUG("m_par->m_y_min+m_par->h_my*ydex=" << m_par->m_y_min << "+" << m_par->h_my << "*" << ydex << "=" << my << ", "); + + double phi = std::atan2(mx,my);//the definition is flipped from what you'd normally think + ATH_MSG_DEBUG("for a phi of " << phi); + if(phi < m_par->minimum_large_phi || phi > m_par->maximum_large_phi){ + ATH_MSG_DEBUG("Chucking phi of " << phi <<" which registers as not in [" << m_par->minimum_large_phi << "," << m_par->maximum_large_phi << "]"); phi=999; } return phi; } -int MMT_Fitter::Rough_ROI_temp(float32fixed<4> theta, float32fixed<4> phi) const{ +int MMT_Fitter::Rough_ROI_temp(double theta, double phi) const{ //temporary function to identify areas of the wedge. - float32fixed<4> minimum_large_theta=m_par->minimum_large_theta, maximum_large_theta=m_par->maximum_large_theta,minimum_large_phi=m_par->minimum_large_phi, maximum_large_phi=m_par->maximum_large_phi; + double minimum_large_theta = m_par->minimum_large_theta; + double maximum_large_theta = m_par->maximum_large_theta; + double minimum_large_phi = m_par->minimum_large_phi; + double maximum_large_phi = m_par->maximum_large_phi; int n_theta_rois=32, n_phi_rois=16;//*** ASK BLC WHAT THESE VALUES OUGHT TO BE! - float32fixed<4> h_theta = (maximum_large_theta - minimum_large_theta)/n_theta_rois; - float32fixed<4> h_phi = (maximum_large_phi - minimum_large_phi)/n_phi_rois; + double h_theta = (maximum_large_theta - minimum_large_theta)/n_theta_rois; + double h_phi = (maximum_large_phi - minimum_large_phi)/n_phi_rois; //how is this done in the FPGA? No such division, for sure - double roi_t = ceil( ( (theta - minimum_large_theta)/h_theta).getFixed()); - double roi_p = ceil( ( (phi - minimum_large_phi)/h_phi).getFixed()); + double roi_t = std::ceil( ( (theta - minimum_large_theta)/h_theta)); + double roi_p = std::ceil( ( (phi - minimum_large_phi)/h_phi)); if(theta<minimum_large_theta || theta>maximum_large_theta) roi_t = 0; if(phi<minimum_large_phi || phi>maximum_large_phi) roi_p = 0; int ret_val=roi_t * 1000 + roi_p; return ret_val; } - diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Hit.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Hit.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0ed24ee33d8650cef66d6eefb5512ea416152159 --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Hit.cxx @@ -0,0 +1,172 @@ +#include "TrigT1NSWSimTools/MMT_Hit.h" +#include "AthenaKernel/getMessageSvc.h" +#include "MuonAGDDDescription/MMDetectorDescription.h" +#include "MuonAGDDDescription/MMDetectorHelper.h" +#include "MuonReadoutGeometry/MuonChannelDesign.h" +#include "MuonReadoutGeometry/MuonDetectorManager.h" +#include "MuonReadoutGeometry/MMReadoutElement.h" + +MMT_Hit::MMT_Hit(char wedge, hitData_entry entry, const MuonGM::MuonDetectorManager* detManager) : AthMessaging(Athena::getMessageSvc(), "MMT_Hit") { + m_sector = wedge; + + std::string module(1, wedge); + module += (TMath::Abs(entry.station_eta) == 1) ? "M1" : "M2"; + m_module = module; + + m_station_name = "MM"; + m_station_name += wedge; + m_VMM_chip = entry.VMM_chip; + m_MMFE_VMM = entry.MMFE_VMM; + m_ART_ASIC = std::ceil(1.*entry.MMFE_VMM/2); + m_station_eta = entry.station_eta; + m_station_phi = entry.station_phi; + m_multiplet = entry.multiplet; + m_gasgap = entry.gasgap; + m_plane = entry.plane; + m_strip = entry.strip; + m_localX = entry.localX; + m_BC_time = entry.BC_time; + m_age = entry.BC_time; + m_Y = -1.; + m_Z = -1.; + m_R = -1.; + m_Ri = -1.; + m_isNoise = false; + m_detManager = detManager; +} + +MMT_Hit::MMT_Hit(const MMT_Hit &hit) : AthMessaging(Athena::getMessageSvc(), "MMT_Hit") { + m_sector = hit.m_sector; + m_module = hit.m_module; + m_station_name = hit.m_station_name; + m_VMM_chip = hit.m_VMM_chip; + m_MMFE_VMM = hit.m_MMFE_VMM; + m_ART_ASIC = hit.m_ART_ASIC; + m_station_eta = hit.m_station_eta; + m_station_phi = hit.m_station_phi; + m_multiplet = hit.m_multiplet; + m_gasgap = hit.m_gasgap; + m_plane = hit.m_plane; + m_strip = hit.m_strip; + m_localX = hit.m_localX; + m_BC_time = hit.m_BC_time; + m_age = hit.m_age; + m_Y = hit.m_Y; + m_Z = hit.m_Z; + m_R = hit.m_R; + m_Ri = hit.m_Ri; + m_RZslope = hit.m_RZslope; + m_YZslope = hit.m_YZslope; + m_isNoise = hit.m_isNoise; + m_detManager = NULL; +} + +MMT_Hit& MMT_Hit::operator=(const MMT_Hit& hit) { + m_sector = hit.m_sector; + m_module = hit.m_module; + m_station_name = hit.m_station_name; + m_VMM_chip = hit.m_VMM_chip; + m_MMFE_VMM = hit.m_MMFE_VMM; + m_ART_ASIC = hit.m_ART_ASIC; + m_station_eta = hit.m_station_eta; + m_station_phi = hit.m_station_phi; + m_multiplet = hit.m_multiplet; + m_gasgap = hit.m_gasgap; + m_plane = hit.m_plane; + m_strip = hit.m_strip; + m_localX = hit.m_localX; + m_BC_time = hit.m_BC_time; + m_age = hit.m_age; + m_Y = hit.m_Y; + m_Z = hit.m_Z; + m_R = hit.m_R; + m_Ri = hit.m_Ri; + m_RZslope = hit.m_RZslope; + m_YZslope = hit.m_YZslope; + m_isNoise = hit.m_isNoise; + m_detManager = NULL; + + return *this; +} + +bool MMT_Hit::isX() const { + int id = this->getPlane(); + return (id == 0 || id == 1 || id == 6 || id == 7) ? true : false; +} + +bool MMT_Hit::isU() const { + int id = this->getPlane(); + return (id == 2 || id == 4) ? true : false; +} + +bool MMT_Hit::isV() const { + int id = this->getPlane(); + return (id == 3 || id == 5) ? true : false; +} + +void MMT_Hit::printHit() { + ATH_MSG_DEBUG("****************** HIT PROPERTIES ******************\n"<< + "\t\t\t\t *** wedge: " << this->getStationName() << "\n"<< + "\t\t\t\t *** VMM: " << this->getVMM() << "\n"<< + "\t\t\t\t *** MMFE: " << this->getMMFE8() << "\n"<< + "\t\t\t\t *** ART_ASIC: " << this->getART() << "\n"<< + "\t\t\t\t *** plane: " << this->getPlane() << "\n"<< + "\t\t\t\t *** st. Eta: " << this->getStationEta() << "\n"<< + "\t\t\t\t *** st. Phi: " << this->getStationPhi() << "\n"<< + "\t\t\t\t *** Multip.: " << this->getMultiplet() << "\n"<< + "\t\t\t\t *** Gas Gap: " << this->getGasGap() << "\n"<< + "\t\t\t\t *** strip: " << this->getChannel() << "\n"<< + "\t\t\t\t *** slopeRZ: " << this->getRZSlope() << "\n"<< + "\t\t\t\t *** BC: " << this->getBC() << "\n"<< + "\t\t\t\t *** X: " << this->getX() << "\n"<< + "\t\t\t\t *** Y: " << this->getY() << "\n"<< + "\t\t\t\t *** Z: " << this->getZ() << "\n"<< + "\t\t\t\t *** R: " << this->getR() << "\n"<< + "\t\t\t\t ****************************************************"); +} + +void MMT_Hit::setHitProperties(const Hit &hit) { + m_Y = hit.info.y; + m_Z = hit.info.z; +} + +void MMT_Hit::updateHitProperties(const MMT_Parameters *par) { + Identifier strip_id = this->getDetManager()->mmIdHelper()->channelID(this->getStationName(), this->getStationEta(), this->getStationPhi(), this->getMultiplet(), this->getGasGap(), this->getChannel()); + const MuonGM::MMReadoutElement* readout = this->getDetManager()->getMMReadoutElement(strip_id); + Amg::Vector3D globalPos(0.0, 0.0, 0.0); + readout->stripGlobalPosition(strip_id, globalPos); + + MMDetectorHelper aHelper; + char side = (globalPos.z() > 0.) ? 'A' : 'C'; + MMDetectorDescription* mm = aHelper.Get_MMDetector(this->getSector(), TMath::Abs(this->getStationEta()), this->getStationPhi(), this->getMultiplet(), side); + MMReadoutParameters roP = mm->GetReadoutParameters(); + + double R = roP.distanceFromZAxis + this->getChannel()*roP.stripPitch - roP.stripPitch/2.; + m_R = R; + m_Z = globalPos.z(); + m_Ri = this->getChannel()*roP.stripPitch; + m_RZslope = R / this->getZ(); + m_Z = this->getZ(); + ATH_MSG_DEBUG("Module: " << this->getModule() << + " HIT --> R_0: " << roP.distanceFromZAxis << + " + Ri " << this->getChannel()*roP.stripPitch << " corresponding to channel " << this->getChannel() << + " ----- Z: " << this->getZ() << ", Plane: " << this->getPlane() << ", eta " << this->getStationEta() << " -- BC: " << this->getBC() << + " RZslope: " << this->getRZSlope()); + + int eta = TMath::Abs(this->getStationEta())-1; + double base = par->ybases[this->getPlane()][eta]; + double Y = base + this->getChannel()*roP.stripPitch - roP.stripPitch/2.; + m_Y = Y; + m_YZslope = Y / this->getZ(); +} + +bool MMT_Hit::verifyHit() { + /* + * Put here all Hit checks, probably redundant if digitization is ok + */ + if (this->getBC() < 0.) return false; + else if (isinf(this->getRZSlope())) return false; + else if (this->getChannel() < 1 || this->getChannel() > 8192) return false; + else return true; +} + diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Road.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Road.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a89be54c507823150b40acb2e9e67a3dc044124a --- /dev/null +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_Road.cxx @@ -0,0 +1,223 @@ +#include "TrigT1NSWSimTools/MMT_Road.h" +#include "AthenaKernel/getMessageSvc.h" + +MMT_Road::MMT_Road(const char sector, const MuonGM::MuonDetectorManager* detManager, const micromegas_t mm, int xthr, int uvthr, int iroadx, int iroadu, int iroadv): AthMessaging(Athena::getMessageSvc(), "MMT_Road") { + m_sector = sector; + m_iroad = iroadx; + m_iroadx = iroadx; + m_iroadu = (iroadu != -1) ? iroadu : iroadx; + m_iroadv = (iroadv != -1) ? iroadv : iroadx; + m_trig = false; + m_xthr = xthr; + m_uvthr = uvthr; + + if (iroadu == -1 && iroadv != -1) ATH_MSG_WARNING("iroadu = -1 but iroadv ain't"); + if (iroadu != -1 && iroadv == -1) ATH_MSG_WARNING("iroadv = -1 but iroadu ain't"); + m_detManager = detManager; + m_micromegas = mm; + m_roadSize = mm.roadSize; + m_roadSizeUpX = mm.nstrip_up_XX; + m_roadSizeDownX = mm.nstrip_dn_XX; + m_roadSizeUpUV = mm.nstrip_up_UV; + m_roadSizeDownUV = mm.nstrip_dn_UV; + m_B = (1./TMath::Tan(1.5/180.*TMath::Pi())); +} + +MMT_Road::~MMT_Road() { + this->reset(); +} + +void MMT_Road::addHits(std::vector<MMT_Hit*> &hits) { + for (auto hit_i : hits) { + int bo = hit_i->getPlane(); + bool has_hit = false; + if( this->containsNeighbors(*hit_i) ) { + for (const auto &hit_j : m_road_hits) { + if (hit_j.getPlane() == bo) { + has_hit = true; + break; + } + } + if (hit_i->isNoise() == false) { + int erase_me = -1; + for (unsigned int j = 0; j < m_road_hits.size(); j++) { + if (m_road_hits[j].getPlane() == bo && m_road_hits[j].isNoise()) { + erase_me = j; + has_hit = false; + break; + } + } + if (erase_me > -1) m_road_hits.erase(m_road_hits.begin() + erase_me); + } + + if (has_hit) continue; + m_road_hits.push_back(*hit_i); + m_road_hits.back().setAge(0); + } + } +} + +bool MMT_Road::containsNeighbors(const MMT_Hit &hit) { + int iroad = 0; + if (hit.isX()) iroad = this->iRoadx(); + else if (hit.isU()) iroad = this->iRoadu(); + else if (hit.isV()) iroad = this->iRoadv(); + else { + std::cerr << "ERROR! Unknown hit type: " << hit.getPlane() << std::endl; + return false; + } + + double pitch = this->getMM().pitch; + double R = (TMath::Abs(hit.getStationEta()) == 1) ? this->getMM().innerRadiusEta1 : this->getMM().innerRadiusEta2; + double Z = hit.getZ(); + + double index = std::round((TMath::Abs(hit.getRZSlope())-0.1)/5e-04); // 0.0005 is approx. the step in slope achievable with a road size of 8 strips + double roundedSlope = 0.1 + index*((0.6 - 0.1)/1000.); + double shift = roundedSlope*(this->getMM().planeCoordinates[hit.getPlane()].Z() - this->getMM().planeCoordinates[0].Z()); + + double olow = 0., ohigh = 0.; + if (hit.isX()) { + olow = this->getRoadSizeDownX()*pitch; + ohigh = this->getRoadSizeUpX()*pitch; + } else { + olow = this->getRoadSizeDownUV()*pitch; + ohigh = this->getRoadSizeUpUV()*pitch; + } + + double slow = (R + (this->getRoadSize()*iroad )*pitch + shift + pitch/2. - olow)/Z; + double shigh = (R + (this->getRoadSize()*(iroad+1))*pitch + shift + pitch/2. + ohigh)/Z; + + double slope = hit.getRZSlope(); + if (this->getSector() != hit.getSector()) { + ATH_MSG_DEBUG("Mismatch between road (" << this->getSector() << ") and hit (" << hit.getSector() << ") sectors"); + return false; + } + if (slope > 0.) return (slope >= slow && slope < shigh); + else return (slope >= shigh && slope < slow); +} + +double MMT_Road::avgSofX() { + std::vector<double> sl; + for (const auto &hit : m_road_hits) { + int bo = hit.getPlane(); + if (bo < 2 || bo > 5) sl.push_back(hit.getRZSlope()); + } + double avg_x = std::accumulate(sl.begin(), sl.end(), 0.0)/(double)sl.size(); + return avg_x; +} + +double MMT_Road::avgSofUV(const int uv1, const int uv2) { + std::vector<double> sl; + for (const auto &hit : m_road_hits) { + int bo = hit.getPlane(); + if (bo == uv1 || bo == uv2) sl.push_back(hit.getRZSlope()); + } + double avg_uv = std::accumulate(sl.begin(), sl.end(), 0.0)/(double)sl.size(); + return avg_uv; +} + +double MMT_Road::avgZofUV(const int uv1, const int uv2) { + std::vector<double> zs; + for (const auto &hit : m_road_hits) { + int bo = hit.getPlane(); + if (bo == uv1 || bo == uv2) zs.push_back(hit.getZ()); + } + double avg_z = std::accumulate(zs.begin(), zs.end(), 0.0)/(double)zs.size(); + return avg_z; +} + +bool MMT_Road::checkCoincidences(const int &bcwind) { + bool passHorizontalCheck = this->horizontalCheck(); + bool passStereoCheck = this->stereoCheck(); + bool passMatureCheck = this->matureCheck(bcwind); + return (passHorizontalCheck && passStereoCheck && passMatureCheck) ? true : false; +} + +unsigned int MMT_Road::countRealHits() { + int nreal = 0; + for (const auto &hit : m_road_hits) { + if (hit.isNoise() == false) nreal++; + } + return nreal; +} + +unsigned int MMT_Road::countUVHits(bool flag) { + unsigned int nuv = 0; + for (const auto &hit : m_road_hits) { + if (hit.getPlane() == 2 || hit.getPlane() == 4) { + if (hit.isNoise() == flag) nuv++; + } + if (hit.getPlane() == 3 || hit.getPlane() == 5) { + if (hit.isNoise() == flag) nuv++; + } + } + return nuv; +} + +unsigned int MMT_Road::countXHits(bool flag) { + unsigned int nx = 0; + for (const auto &hit : m_road_hits) { + if (hit.getPlane() < 2 || hit.getPlane() > 5) { + if (hit.isNoise() == flag) nx++; + } + } + return nx; +} + +bool MMT_Road::horizontalCheck() { + int nx1 = 0, nx2 = 0; + for (const auto &hit : m_road_hits) { + if (hit.getPlane() >-1 && hit.getPlane() < 2) nx1++; + if (hit.getPlane() > 5 && hit.getPlane() < 8) nx2++; + } + if (nx1 > 0 && nx2 > 0 && (nx1+nx2) >= this->getXthreshold()) return true; + return false; +} + +void MMT_Road::incrementAge(const int &bcwind) { + std::vector<unsigned int> old_ihits; + for (unsigned int j = 0; j < m_road_hits.size(); j++) { + m_road_hits[j].setAge(m_road_hits[j].getAge() +1); + if (m_road_hits[j].getAge() > (bcwind-1)) old_ihits.push_back(j); + } + for (int j = old_ihits.size()-1; j > -1; j--) m_road_hits.erase(m_road_hits.begin()+j); +} + +bool MMT_Road::matureCheck(const int &bcwind) { + for (const auto &hit : m_road_hits) { + if (hit.getAge() == (bcwind - 1)) return true; + } + return false; +} + +double MMT_Road::mxl() { + std::vector<double> ys, zs; + for (const auto &hit : m_road_hits) { + int bo = hit.getPlane(); + if (bo < 2 || bo > 5) { + ys.push_back(hit.getR()); + zs.push_back(hit.getZ()); + } + } + double mxl = 0; + double avg_z = std::accumulate(zs.begin(), zs.end(), 0.0)/(double)zs.size(); + double sum_sq_z = std::inner_product(zs.begin(), zs.end(), zs.begin(), 0.0); + for (unsigned int i = 0; i < ys.size(); i++) mxl += ys[i]*( (zs[i]-avg_z) / (sum_sq_z - zs.size()*TMath::Power(avg_z,2)) ); + + return mxl; +} + +void MMT_Road::reset() { + if (!m_road_hits.empty()) m_road_hits.clear(); +} + +bool MMT_Road::stereoCheck() { + int nu = 0, nv = 0; + for (const auto &hit : m_road_hits) { + if (hit.getPlane() == 2 || hit.getPlane() == 4) nu++; + if (hit.getPlane() == 3 || hit.getPlane() == 5) nv++; + } + if (this->getUVthreshold() == 0) return true; + if (nu > 0 && nv > 0 && (nu+nv) >= this->getUVthreshold()) return true; + return false; +} diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_struct.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_struct.cxx index 84525fbc3802ccfa9f6467f2eab562bd65a02fc3..3a3578b397f005e0aadff6df31b8550c488c9ad0 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_struct.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMT_struct.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "TrigT1NSWSimTools/MMT_struct.h" @@ -25,9 +25,10 @@ using std::setw; const float MMTStructConst = 8192.; -std_align::std_align(int qcm,const TVector3& trans,const TVector3& ang):type(qcm),translate(trans),rotate(ang){ +std_align::std_align(int qcm,const ROOT::Math::XYZVector& trans,const ROOT::Math::XYZVector& ang):type(qcm),translate(trans),rotate(ang){ if(type==0){ - translate=TVector3();rotate=TVector3(); + translate=ROOT::Math::XYZVector(); + rotate=ROOT::Math::XYZVector(); } } @@ -43,7 +44,7 @@ string std_align::par_title(int par_num,bool small_unit)const{ throw std::runtime_error("MMT_struct: Invalid parameter index"); } if(par_num>2){ - if(small_unit)par_title+=" [mrad]"; + if(small_unit) par_title+=" [mrad]"; else par_title+=" [rad]"; } return par_title; @@ -64,50 +65,55 @@ string std_align::par_name(int par_num)const{ } string std_align::par_title_val(int par_num)const{ - assert(par_num>=0&&par_num<parmax()); - ostringstream par_title;bool do_it=false; + assert(par_num>=0 && par_num<parmax()); + ostringstream par_title; + bool do_it=false; for(int i=0;i<6;i++){ - if(par_val(i)!=0)do_it=true; + if(par_val(i)!=0) do_it=true; } - if(!do_it)return "nominal"; - if(par_num==0) par_title<<"#Deltas="<<par_val(0)<<" mm"; - else if(par_num==1) par_title<<"#Deltaz="<<par_val(1)<<" mm"; - else if(par_num==2) par_title<<"#Deltat="<<par_val(2)<<" mm"; - else if(par_num==3) par_title<<"#gamma_{s}="<<par_val(3)<<" rad"; - else if(par_num==4) par_title<<"#beta_{z}="<<par_val(4)<<" rad"; - else if(par_num==5) par_title<<"#alpha_{t}="<<par_val(5)<<" rad"; + if(!do_it) return "nominal"; + if(par_num==0) par_title << "#Deltas=" << par_val(0) << " mm"; + else if(par_num==1) par_title << "#Deltaz=" << par_val(1)<< " mm"; + else if(par_num==2) par_title << "#Deltat=" << par_val(2)<< " mm"; + else if(par_num==3) par_title << "#gamma_{s}=" << par_val(3) << " rad"; + else if(par_num==4) par_title << "#beta_{z}=" << par_val(4) << " rad"; + else if(par_num==5) par_title << "#alpha_{t}=" << par_val(5) << " rad"; return par_title.str(); } string std_align::par_name_val(int par_num)const{ - assert(par_num>=0&&par_num<parmax()); + assert(par_num>=0 && par_num<parmax()); ostringstream par_name; - if(par_num==0) par_name<<"dts"; - else if(par_num==1) par_name<<"dtz"; - else if(par_num==2) par_name<<"dtt"; - else if(par_num==3) par_name<<"drs"; - else if(par_num==4) par_name<<"drz"; - else if(par_num==5) par_name<<"drt"; + if(par_num==0) par_name << "dts"; + else if(par_num==1) par_name << "dtz"; + else if(par_num==2) par_name << "dtt"; + else if(par_num==3) par_name << "drs"; + else if(par_num==4) par_name << "drz"; + else if(par_num==5) par_name << "drt"; par_name<<par_val(par_num); return par_name.str(); } string std_align::print()const{ if(type==0) return string("_NOM");//nominal - string gordon;int maxpar=0; + string gordon; + int maxpar=0; if(type==1){ - gordon+="_MIS";maxpar=parmax(); + gordon+="_MIS"; + maxpar=parmax(); } else if(type==2){ - gordon+="_COR";maxpar=parmax(); + gordon+="_COR"; + maxpar=parmax(); } else if(type==3){ - gordon+="_3CR";maxpar=parmax(); + gordon+="_3CR"; + maxpar=parmax(); } for(int imal=0; imal<maxpar;imal++){ ostringstream addme; if(par_val(imal)!=0){ - addme<<"_"<<par_name(imal)<<par_val(imal); + addme << "_" << par_name(imal) << par_val(imal); } gordon+=addme.str(); } @@ -115,7 +121,7 @@ string std_align::print()const{ } double std_align::par_val(int par_num) const{ - assert(par_num>=0&&par_num<parmax()); + assert(par_num>=0 && par_num<parmax()); double par_val=-999.; if(par_num==0) par_val=translate.X(); else if(par_num==1) par_val=translate.Y(); @@ -127,8 +133,8 @@ double std_align::par_val(int par_num) const{ } void std_align::set_val(int par_num,double par_val){ - assert(par_num>=0&&par_num<parmax()); - if(log10(abs(par_val))<-7)par_val=0; + assert(par_num>=0 && par_num<parmax()); + if(std::log10(std::abs(par_val))<-7) par_val=0; if(par_num==0) translate.SetX(par_val); else if(par_num==1) translate.SetY(par_val); else if(par_num==2) translate.SetZ(par_val); @@ -144,10 +150,10 @@ string std_align::detail()const{ else if(type==2) logan+="COR"; for(int imal=0; imal<parmax();imal++){ ostringstream addme; - if(par_val(imal)!=0)addme<<", "<<par_title(imal)<<"="<<par_val(imal); + if(par_val(imal)!=0) addme << ", " << par_title(imal) << "=" << par_val(imal); logan+=addme.str(); } - if(logan.compare("")==0)logan+=" nominal"; + if(logan.compare("")==0) logan+=" nominal"; return logan; } @@ -162,7 +168,7 @@ gcm_key::gcm_key(int the_pt,int the_ct,int the_mis,int the_corr,int the_eta,int :pt(the_pt),ct(the_ct),mis(the_mis),correct(the_corr),eta(the_eta),qt(the_qt),bgcode(the_bg){} bool gcm_key::operator==(const gcm_key& rhs) const{ - if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis==rhs.mis&&this->correct==rhs.correct&&this->eta==rhs.eta&&this->qt==rhs.qt&&this->bgcode==rhs.bgcode) return true; + if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis==rhs.mis && this->correct==rhs.correct && this->eta==rhs.eta && this->qt==rhs.qt && this->bgcode==rhs.bgcode) return true; return false; } @@ -171,12 +177,12 @@ bool gcm_key::operator!=(const gcm_key& rhs) const{ } bool gcm_key::operator<(const gcm_key& rhs) const{ if(this->pt<rhs.pt) return true; - else if(this->pt==rhs.pt&&this->ct<rhs.ct) return true; - else if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis<rhs.mis) return true; - else if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis==rhs.mis&&this->correct<rhs.correct) return true; - else if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis==rhs.mis&&this->correct==rhs.correct&&this->eta<rhs.eta) return true; - else if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis==rhs.mis&&this->correct==rhs.correct&&this->eta==rhs.eta&&this->qt<rhs.qt) return true; - else if(this->pt==rhs.pt&&this->ct==rhs.ct&&this->mis==rhs.mis&&this->correct==rhs.correct&&this->eta==rhs.eta&&this->qt==rhs.qt&&this->bgcode<rhs.bgcode) return true; + else if(this->pt==rhs.pt && this->ct<rhs.ct) return true; + else if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis<rhs.mis) return true; + else if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis==rhs.mis && this->correct<rhs.correct) return true; + else if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis==rhs.mis && this->correct==rhs.correct && this->eta<rhs.eta) return true; + else if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis==rhs.mis && this->correct==rhs.correct && this->eta==rhs.eta && this->qt<rhs.qt) return true; + else if(this->pt==rhs.pt && this->ct==rhs.ct && this->mis==rhs.mis && this->correct==rhs.correct && this->eta==rhs.eta && this->qt==rhs.qt && this->bgcode<rhs.bgcode) return true; return false; } bool gcm_key::operator<=(const gcm_key& rhs) const{ @@ -190,30 +196,31 @@ bool gcm_key::operator>=(const gcm_key& rhs) const{ } string gcm_key::print()const{ - ostringstream ricola;ricola<<"(p"<<pt<<",C"<<ct<<",m"<<mis<<",c"<<correct<<",e"<<eta<<",q"<<qt<<",b"<<bgcode<<")"; + ostringstream ricola; + ricola << "(p" << pt << ",C" << ct << ",m" << mis << ",c" << correct << ",e" << eta << ",q" << qt << ",b" << bgcode << ")"; return ricola.str(); } void gcm_key::set_var(int var,int val){ - assert(var>=0&&var<varmax()); - if(var==0)pt=val; - else if(var==1)ct=val; - else if(var==2)mis=val; - else if(var==3)correct=val; - else if(var==4)eta=val; - else if(var==5)qt=val; - else if(var==6)bgcode=val; + assert(var>=0 && var<varmax()); + if(var==0) pt=val; + else if(var==1) ct=val; + else if(var==2) mis=val; + else if(var==3) correct=val; + else if(var==4) eta=val; + else if(var==5) qt=val; + else if(var==6) bgcode=val; } int gcm_key::get_var(int var)const{ - assert(var>=0&&var<varmax()); - if(var==0)return pt; - else if(var==1)return ct; - else if(var==2)return mis; - else if(var==3)return correct; - else if(var==4)return eta; - else if(var==5)return qt; - else if(var==6)return bgcode; + assert(var>=0 && var<varmax()); + if(var==0) return pt; + else if(var==1) return ct; + else if(var==2) return mis; + else if(var==3) return correct; + else if(var==4) return eta; + else if(var==5) return qt; + else if(var==6) return bgcode; return -999; } @@ -223,17 +230,18 @@ par_par::par_par(double the_h,int xct,int uvct,double uver,const string& set,boo string par_par::print_pars(const vector<int>&hide) const{ vector<bool>sho(gcm_key().varmax(),true); for(int i=0;i<(int)hide.size();i++){ - if(hide[i]<0||hide[i]>=gcm_key().varmax())continue; - sho[hide[i]]=false;sho[2]=false; + if(hide[i]<0 || hide[i]>=gcm_key().varmax()) continue; + sho[hide[i]]=false; + sho[2]=false; } ostringstream mars; - if(sho[6])mars<<"_h"<<h; - if(sho[1])mars<<"_ctx"<<ctx<<"_ctuv"<<ctuv; - if(sho[6])mars<<"_uverr"<<uverr; - mars<<"_set"<<setup<<"_ql"<<islarge<<"_qdlm"<<q_dlm; - if(sho[6])mars<<"_qbg"<<genbg; - if(sho[5])mars<<"_qt"<<qt; - if(sho[2])mars<<misal.print()<<corr.print(); + if(sho[6]) mars << "_h" << h; + if(sho[1]) mars << "_ctx" << ctx << "_ctuv" << ctuv; + if(sho[6]) mars << "_uverr" << uverr; + mars << "_set" << setup << "_ql" << islarge << "_qdlm" << q_dlm; + if(sho[6]) mars << "_qbg" << genbg; + if(sho[5]) mars << "_qt" << qt; + if(sho[2]) mars << misal.print() << corr.print(); return mars.str(); } @@ -243,9 +251,9 @@ string par_par::detail()const{ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM::MuonDetectorManager* detManager) : AthMessaging(Athena::getMessageSvc(), "MMT_Parameters") { - if(inputParams.misal.is_nominal())inputParams.misal.type=0; + if(inputParams.misal.is_nominal()) inputParams.misal.type=0; //can still do sim_corrections for nominal - if(inputParams.corr.is_nominal()&&inputParams.corr.type==2)inputParams.corr.type=0; + if(inputParams.corr.is_nominal() && inputParams.corr.type==2) inputParams.corr.type=0; n_etabins=2; n_etabins=10; n_phibins=10; @@ -253,10 +261,10 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM diag=true; fill0=false; //stuff pulled from the inputParams - h=float32fixed<2>(inputParams.h); + h=inputParams.h; CT_x=inputParams.ctx; CT_uv=inputParams.ctuv; - uv_error=float32fixed<2>(inputParams.uverr); + uv_error=inputParams.uverr; setup=(inputParams.setup); islarge=inputParams.islarge; dlm_new=inputParams.q_dlm; @@ -265,15 +273,15 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM misal=inputParams.misal; misalign=(misal.type==1); correct=inputParams.corr; - ybins=8; n_stations_eta=(dlm_new?2:4); + ybins=8; + n_stations_eta=(dlm_new ? 2 : 4); val_tbl=inputParams.fill_val; - // Get the modules for each multiplet, the sector and side arguement (third and fifth argument, respectively) shouldn't matter // since the algorithm computes locally char side = 'A'; - char sector = wedgeSize; + sector = wedgeSize; MMDetectorHelper aHelper; - MMDetectorDescription* mm_top_mult1 = aHelper.Get_MMDetector(sector, 2, 5, 1, side); + MMDetectorDescription* mm_top_mult1 = aHelper.Get_MMDetector(sector, 2, 5, 1, side); // Parameters: sector, eta, phi, layer/multiplet MMDetectorDescription* mm_top_mult2 = aHelper.Get_MMDetector(sector, 2, 5, 2, side); MMDetectorDescription* mm_bottom_mult1 = aHelper.Get_MMDetector(sector, 1, 5, 1, side); MMDetectorDescription* mm_bottom_mult2 = aHelper.Get_MMDetector(sector, 1, 5, 2, side); @@ -289,26 +297,67 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM if(!islarge){ ATH_MSG_WARNING("We haven't configured the small wedge parameters yet! Go bother the developer...kindly...this will likely crash now!\n"); - return; + //return; TODO: is this still the case? } //y = vertical distance from beamline //z = into wedge along beamline //x = horizontal distance from beam looking down + + if (sector == 'L') { + ATH_MSG_INFO("LM1 \n" << + "\t\t\t\t Total Strips: " << roParam_bottom_mult1.tStrips << " with pitch: " << roParam_bottom_mult1.stripPitch << "\n" << + "\t\t\t\t KO strips TopEta: " << roParam_bottom_mult1.nMissedTopEta << " - BottomEta: " << roParam_bottom_mult1.nMissedBottomEta << "\n" << + "\t\t\t\t KO strips TopStereo: " << roParam_bottom_mult1.nMissedTopStereo << " - BottomStereo: " << roParam_bottom_mult1.nMissedBottomStereo << "\n" << + "\t\t\t\t (Top, Bottom, Length): (" << mm_bottom_mult1->lWidth() << ", " << mm_bottom_mult1->sWidth() << ", " << mm_bottom_mult1->Length() << ")\n" << + "\t\t\t\t Zpos - Distance from ZAxis: " << roParam_bottom_mult1.zpos << " - " << roParam_bottom_mult1.distanceFromZAxis << "\n" << + "\t\t\t\t dlStereoTop/Bottom: " << roParam_bottom_mult1.dlStereoTop << " " << roParam_bottom_mult1.dlStereoBottom << "\n" << + "\t\t\t\t Active area --> (Bottom, Top, Height) : (" << roParam_bottom_mult1.activeBottomLength << ", " << roParam_bottom_mult1.activeTopLength << ", " << roParam_bottom_mult1.activeH << ")"); + } else if (sector == 'S') { + ATH_MSG_INFO("SM1 \n" << + "\t\t\t\t KO strips TopEta: " << roParam_bottom_mult1.nMissedTopEta << " - BottomEta: " << roParam_bottom_mult1.nMissedBottomEta << "\n" << + "\t\t\t\t KO strips TopStereo: " << roParam_bottom_mult1.nMissedTopStereo << " - BottomStereo: " << roParam_bottom_mult1.nMissedBottomStereo << "\n" << + "\t\t\t\t Total Strips: " << roParam_bottom_mult1.tStrips << " with pitch: " << roParam_bottom_mult1.stripPitch << "\n" << + "\t\t\t\t (Top, Bottom, Length): (" << mm_bottom_mult1->lWidth() << ", " << mm_bottom_mult1->sWidth() << ", " << mm_bottom_mult1->Length() << ")\n" << + "\t\t\t\t Zpos - Distance from ZAxis: " << roParam_bottom_mult1.zpos << " - " << roParam_bottom_mult1.distanceFromZAxis << "\n" << + "\t\t\t\t dlStereoTop/Bottom: " << roParam_bottom_mult1.dlStereoTop << " " << roParam_bottom_mult1.dlStereoBottom << "\n" << + "\t\t\t\t Active area --> (Bottom, Top, Height) : (" << roParam_bottom_mult1.activeBottomLength << ", " << roParam_bottom_mult1.activeTopLength << ", " << roParam_bottom_mult1.activeH << ")"); + } + + if (sector == 'L') { + ATH_MSG_INFO("LM2 \n" << + "\t\t\t\t Total Strips: " << roParam_top_mult1.tStrips << " with pitch: " << roParam_top_mult1.stripPitch << "\n" << + "\t\t\t\t KO strips TopEta: " << roParam_top_mult1.nMissedTopEta << " - BottomEta: " << roParam_top_mult1.nMissedBottomEta << "\n" << + "\t\t\t\t KO strips TopStereo: " << roParam_top_mult1.nMissedTopStereo << " - BottomStereo: " << roParam_top_mult1.nMissedBottomStereo << "\n" << + "\t\t\t\t (Top, Bottom, Length): (" << mm_top_mult1->lWidth() << ", " << mm_top_mult1->sWidth() << ", " << mm_top_mult1->Length() << ")\n" << + "\t\t\t\t Zpos - Distance from ZAxis: " << roParam_top_mult1.zpos << " - " << roParam_top_mult1.distanceFromZAxis << "\n" << + "\t\t\t\t dlStereoTop/Bottom: " << roParam_top_mult1.dlStereoTop << " " << roParam_top_mult1.dlStereoBottom << "\n" << + "\t\t\t\t Active area --> (Top, Bottom, Height) : (" << roParam_top_mult1.activeTopLength << ", " << roParam_top_mult1.activeBottomLength << ", " << roParam_top_mult1.activeH << ")"); + } else if (sector == 'S') { + ATH_MSG_INFO("SM2 " << + "\t\t\t\t KO strips TopEta: " << roParam_top_mult1.nMissedTopEta << " - BottomEta: " << roParam_top_mult1.nMissedBottomEta << "\n" << + "\t\t\t\t KO strips TopStereo: " << roParam_top_mult1.nMissedTopStereo << " - BottomStereo: " << roParam_top_mult1.nMissedBottomStereo << "\n" << + "\t\t\t\t (Top, Bottom, Length): (" << mm_top_mult1->lWidth() << ", " << mm_top_mult1->sWidth() << ", " << mm_top_mult1->Length() << ")\n" << + "\t\t\t\t Zpos - Distance from ZAxis: " << roParam_top_mult1.zpos << " - " << roParam_top_mult1.distanceFromZAxis << "\n" << + "\t\t\t\t dlStereoTop/Bottom: " << roParam_top_mult1.dlStereoTop << " " << roParam_top_mult1.dlStereoBottom << "\n" << + "\t\t\t\t Active area --> (Top, Bottom, Height) : (" << roParam_top_mult1.activeTopLength << ", " << roParam_top_mult1.activeBottomLength << ", " << roParam_top_mult1.activeH << ")"); + for (const auto &angle : roParam_top_mult1.stereoAngle) ATH_MSG_INFO("Stereo angle: " << angle); + } + //////////// Define the large wedge ///////////////// - w1=float32fixed<18>(mm_top_mult1->lWidth()); //top - w2=float32fixed<18>(mm_top_mult1->lWidth()*1./(cos(roParam_top_mult1.stereoAngle.at(3)))); //determined by 33deg angle //shelves part - w3=float32fixed<18>(mm_bottom_mult1->sWidth());//582.3); //bottom - h1=float32fixed<18>(roParam_top_mult1.roLength+roParam_bottom_mult1.roLength+5.0); //how tall wedge is at w1, ie total height + w1=double(mm_top_mult1->lWidth()); //top + w2=double(mm_top_mult1->lWidth()*1./(std::cos(roParam_top_mult1.stereoAngle[3]))); //determined by 33deg angle //shelves part + w3=double(mm_bottom_mult1->sWidth());//582.3); //bottom + h1=double(roParam_top_mult1.roLength+roParam_bottom_mult1.roLength+5.0); //how tall wedge is at w1, ie total height - wedge_opening_angle = float32fixed<18>(33.); //degree + wedge_opening_angle = 33.; //degree - z_large=vector<vector<float32fixed<18> > >(ybins,vector<float32fixed<18> >(setup.length(),float32fixed<18>(0.))); + z_large=vector<vector<double > >(ybins,vector<double >(setup.length(),0.)); // retrieve the z-position of the planes z_nominal.clear(); int eta = 1; for (auto pos: MM_firststrip_positions(detManager, wedgeString, eta)){ - z_nominal.push_back(float32fixed<18>(pos.z())); + z_nominal.push_back(double(pos.z())); } if(z_nominal.size() != setup.size()){ @@ -317,37 +366,38 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM throw std::runtime_error("MMT_Parameters: Invalid number of planes"); } - mid_plane_large=float32fixed<18>(0.); - for(unsigned int iz=0; iz<z_nominal.size(); iz++) mid_plane_large+=z_nominal[iz].getFixed(); - mid_plane_large/=z_nominal.size(); + mid_plane_large=0.; + for(unsigned int iz=0; iz<z_nominal.size(); iz++) mid_plane_large += z_nominal[iz]; + mid_plane_large /= z_nominal.size(); - vector<int> xp=q_planes("x"),up=q_planes("u"),vp=q_planes("v"); - mid_plane_large_X=float32fixed<18>(0.); - for(unsigned int ix=0;ix<xp.size();ix++) mid_plane_large_X+=z_nominal[xp[ix]].getFixed(); + vector<int> xp=q_planes("x"), up=q_planes("u"), vp=q_planes("v"); + mid_plane_large_X=0.; + for(unsigned int ix=0;ix<xp.size();ix++) mid_plane_large_X += z_nominal[xp[ix]]; mid_plane_large_X /= 1.*xp.size(); - mid_plane_large_UV=float32fixed<18>(0.); - for(unsigned int iu=0;iu<up.size();iu++) mid_plane_large_UV+=z_nominal[up[iu]].getFixed(); - for(unsigned int iv=0;iv<vp.size();iv++) mid_plane_large_UV+=z_nominal[vp[iv]].getFixed(); + mid_plane_large_UV=0.; + for(unsigned int iu=0;iu<up.size();iu++) mid_plane_large_UV+=z_nominal[up[iu]]; + for(unsigned int iv=0;iv<vp.size();iv++) mid_plane_large_UV+=z_nominal[vp[iv]]; mid_plane_large_UV /= 1.*(up.size()+vp.size()); - H=float32fixed<18>(roParam_bottom_mult1.distanceFromZAxis);//982.); //bottom of wedge to the beamline - Hnom=float32fixed<18>(roParam_bottom_mult1.distanceFromZAxis);//982.); //bottom of wedge to the beamline + H=double(roParam_bottom_mult1.distanceFromZAxis);//982.); //bottom of wedge to the beamline + Hnom=double(roParam_bottom_mult1.distanceFromZAxis);//982.); //bottom of wedge to the beamline + - strip_width = float32fixed<4>(roParam_top_mult1.stripPitch); // 0.5; - stereo_degree = float32fixed<4>(TMath::RadToDeg()*roParam_top_mult1.stereoAngle.at(2)); //0.75 //3 in degrees! - float32fixed<2> degree=roParam_top_mult1.stereoAngle.at(2); - vertical_strip_width_UV = strip_width.getFixed()/cos(degree.getFixed()); - ybases=vector<vector<float32fixed<18> > >(setup.size(),vector<float32fixed<18> >(n_stations_eta,float32fixed<18>(0.))); + strip_width = double(roParam_top_mult1.stripPitch); // 0.5; + stereo_degree = double(TMath::RadToDeg()*roParam_top_mult1.stereoAngle[2]); //0.75 //3 in degrees! + double degree=roParam_top_mult1.stereoAngle[2]; + vertical_strip_width_UV = strip_width/std::cos(degree); + ybases=vector<vector<double > >(setup.size(), vector<double >(n_stations_eta,0.)); //gposx for X1; gpos-lpos X2,UV1,UV2 - vector<float32fixed<18> > xeta,ueta,veta; + vector<double > xeta,ueta,veta; // MM_firststrip_positions returns the position for phi sector 1. // => for the small sectors, rotate this by -1/16 of a rotation to make our lives easier. // in this coordinate basis, x is up/down the wedge (radial), and y is left/right (phi). - float cos_rotation = cos(-2*TMath::Pi() / 16.0); - float sin_rotation = sin(-2*TMath::Pi() / 16.0); + float cos_rotation = std::cos(-2*TMath::Pi() / 16.0); + float sin_rotation = std::sin(-2*TMath::Pi() / 16.0); float x_rotated = 0.0; float y_rotated = 0.0; @@ -370,13 +420,13 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM y_rotated = pos.y(); } - if (is_u(layer)) st_angle = -1*abs(stereo_degree.getFixed()); - else if (is_v(layer)) st_angle = abs(stereo_degree.getFixed()); + if (is_u(layer)) st_angle = -1*abs(stereo_degree); + else if (is_v(layer)) st_angle = abs(stereo_degree); else st_angle = 0; // walk from the center of the strip to the position of the strip at the center of the wedge. // NB: for X-planes, this is simply the center of the strip: tan(0) = 0. - radial_pos = abs(x_rotated - y_rotated*tan(st_angle * TMath::Pi()/180.0)); + radial_pos = abs(x_rotated - y_rotated*std::tan(st_angle * TMath::Pi()/180.0)); if (is_x(layer) && eta==1) radial_pos_xx_1 = radial_pos; if (is_x(layer) && eta==2) radial_pos_xx_2 = radial_pos; @@ -388,8 +438,8 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM } // store radial positions as fixed point - std::vector< float32fixed<18> > radial_pos_xx = {float32fixed<18>(radial_pos_xx_1), float32fixed<18>(radial_pos_xx_2)}; - std::vector< float32fixed<18> > radial_pos_uv = {float32fixed<18>(radial_pos_uv_1), float32fixed<18>(radial_pos_uv_2)}; + std::vector< double > radial_pos_xx = {double(radial_pos_xx_1), double(radial_pos_xx_2)}; + std::vector< double > radial_pos_uv = {double(radial_pos_uv_1), double(radial_pos_uv_2)}; ybases[0] = radial_pos_xx; ybases[1] = radial_pos_xx; ybases[2] = radial_pos_uv; @@ -399,58 +449,65 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM ybases[6] = radial_pos_xx; ybases[7] = radial_pos_xx; - //now put in the positions at `ly spaced points for a y dependent z bool hack=false; - z_large=vector<vector<float32fixed<18> > >(ybins,z_nominal); - double pitch_f=1.*sin(correct.rotate.X())*h1.getFixed()/ybins,pitch_b=0,bumper_up=0.0; - if(hack){double factor=-1;pitch_f*=factor;pitch_b*=factor;} - ATH_MSG_DEBUG("Specs: correct.rotate.X()="<<correct.rotate.X()<<",correct.translate.Z()="<<correct.translate.Z()<<",pitch_f="<<pitch_f); + z_large=vector<vector<double > >(ybins,z_nominal); + double pitch_f=1.*std::sin(correct.rotate.X())*h1/ybins, pitch_b=0, bumper_up=0.0; + if(hack) { + double factor=-1; + pitch_f*=factor; + pitch_b*=factor; + } + ATH_MSG_DEBUG("Specs: correct.rotate.X()=" << correct.rotate.X() << ",correct.translate.Z()=" << correct.translate.Z() << ",pitch_f=" << pitch_f); for(int iy=0;iy<ybins;iy++){ //z axis in misalignment line points in the opposite direction....TOO LATE! (right handed coordinate system the way we expect for x and y makes this happen) - double over_f=pitch_f*(iy+bumper_up)-correct.translate.Z(),over_b=0;//pitch_b*(iy+bumper_up)+correct.dzb; - ATH_MSG_DEBUG("iy="<<iy<<"over_f="<<over_f<<",over_b="<<over_b); + double over_f=pitch_f*(iy+bumper_up)-correct.translate.Z(); + double over_b=0;//pitch_b*(iy+bumper_up)+correct.dzb; + ATH_MSG_DEBUG("iy=" << iy << "over_f=" << over_f << ",over_b=" << over_b); for(int jp=0;jp<8;jp++){ - ATH_MSG_DEBUG("z_large["<<iy<<"]["<<jp<<"]"<<z_large[iy][jp].getFixed()<<"--->"); - if(jp<4)z_large[iy][jp]+=over_f; + ATH_MSG_DEBUG("z_large[" << iy << "][" << jp << "]" << z_large[iy][jp] << "--->"); + if(jp<4) z_large[iy][jp]+=over_f; else z_large[iy][jp]+=over_b; - ATH_MSG_DEBUG(z_large[iy][jp].getFixed()); + ATH_MSG_DEBUG(z_large[iy][jp]); } } //////// TABLE GENERATORS /////////////// //size of cartesian steps - h_mx = float32fixed<2>(0.0001); // 0.005; //0.001; - h_my = float32fixed<2>(0.0001); //0.005; //0.001; + h_mx = 0.0001; // 0.005; //0.001; + h_my = 0.0001; //0.005; //0.001; m_y_max = ( (Hnom+h1)/z_nominal.front() ); m_y_min = ( H/z_nominal.back() ); m_x_max = (w2/2.)/z_nominal.front(); m_x_min = (w2/-2.)/z_nominal.back(); //-2; - n_x = ceil((m_x_max - m_x_min).getFixed()/h_mx.getFixed()); - n_y = ceil((m_y_max - m_y_min).getFixed()/h_my.getFixed()); + n_x = std::ceil((m_x_max - m_x_min)/h_mx); + n_y = std::ceil((m_y_max - m_y_min)/h_my); ///////////////////////////////////////// ////////// for cut applications ///////////////// double tol = 0; //0.02; //BLC had some interesting bounds...let's do the ones that make sense to me - minimum_large_theta = float32fixed<4>(atan(H.getFixed()/z_nominal.back().getFixed())+tol); - maximum_large_theta = float32fixed<4>(atan(sqrt(pow( (Hnom+h1).getFixed(),2)+0.25*pow(w1.getFixed(),2))/z_nominal.back().getFixed())-tol); - minimum_large_phi = float32fixed<4>(-TMath::DegToRad()*0.5*wedge_opening_angle.getFixed())+tol; maximum_large_phi = float32fixed<4>(TMath::DegToRad()*(0.5*wedge_opening_angle.getFixed())-tol); + minimum_large_theta = std::atan(H/z_nominal.back()) + tol; + maximum_large_theta = std::atan(std::sqrt(std::pow( (Hnom+h1),2) + 0.25*std::pow(w1,2))/z_nominal.back()) - tol; + minimum_large_phi = -TMath::DegToRad()*0.5*wedge_opening_angle + tol; + maximum_large_phi = TMath::DegToRad()*(0.5*wedge_opening_angle) - tol; /////////////////////////////////////////////////// - double phiseg= ((maximum_large_phi-minimum_large_phi)/n_phibins).getFixed(); + double phiseg= ((maximum_large_phi-minimum_large_phi)/n_phibins); m_phibins.clear(); - for(int i=0;i<=n_phibins;i++) m_phibins.push_back( (minimum_large_phi+phiseg*i).getFixed() ); - double etalo=-log(tan(0.5*maximum_large_theta.getFixed())),etahi=-log(tan(0.5*minimum_large_theta.getFixed()));bool custom=false; + for(int i=0;i<=n_phibins;i++) m_phibins.push_back( (minimum_large_phi+phiseg*i) ); + double etalo=-std::log(std::tan(0.5*maximum_large_theta)); + double etahi=-std::log(std::tan(0.5*minimum_large_theta)); + bool custom=false; if(custom){ //these custom eta bins have the top station as one bin (since that has normal shape for the phi resolutions) //and separates the old part from gap to 2.5 (where the simulation used to be limited) and then one bin from 2.5 to max (new stuff) m_etabins.clear(); m_etabins.push_back(etalo); - double etamid=asinh( (z_nominal.front()/ybases.front()[1] ).getFixed() ); + double etamid=std::asinh( (z_nominal.front()/ybases.front()[1] ) ); double intermid=(2.5-etamid)/8.; - for(int i=0;i<8;i++)m_etabins.push_back(etamid+intermid*i); + for(int i=0;i<8;i++) m_etabins.push_back(etamid+intermid*i); m_etabins.push_back(2.5); m_etabins.push_back(etahi); n_etabins=m_etabins.size()-1; @@ -458,12 +515,12 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM else if(n_etabins==(int)xeta.size()){ m_etabins.push_back(etalo); for(int i=n_etabins-1;i>=0;i--){ - m_etabins.push_back(asinh( (z_nominal.front()/xeta[i]).getFixed() )); + m_etabins.push_back(std::asinh( (z_nominal.front()/xeta[i]) )); } } else{ double increment=(etahi-etalo)/n_etabins*1.; - for(int i=0;i<=n_etabins;i++)m_etabins.push_back(etalo+increment*i); + for(int i=0;i<=n_etabins;i++) m_etabins.push_back(etalo+increment*i); } /////// Rough ROI //////////// @@ -472,11 +529,11 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM ////////////////////////////// /////// Front Filter ///////////// - double theta_max=maximum_large_theta.getFixed(); - double theta_min=minimum_large_theta.getFixed(); + double theta_max=maximum_large_theta; + double theta_min=minimum_large_theta; - slope_max = float32fixed<3>(tan(theta_max)); - slope_min = float32fixed<3>(tan(theta_min)); + slope_max = double(std::tan(theta_max)); + slope_min = double(std::tan(theta_min)); CT_u = 0; CT_v = 0; @@ -496,7 +553,7 @@ MMT_Parameters::MMT_Parameters(par_par inputParams, char wedgeSize, const MuonGM Delta_theta_optimization_LG(); //simulation-based correction fill...if available - if(correct.type==3)fill_crep_table(inputParams.pcrep_dir,inputParams.tag); + if(correct.type==3) fill_crep_table(inputParams.pcrep_dir,inputParams.tag); fill_yzmod(); } } @@ -507,7 +564,8 @@ std::vector<Amg::Vector3D> MMT_Parameters::MM_firststrip_positions(const MuonGM: positions.clear(); Identifier strip_id; - int phi = 1, strip = 1; + int phi = 1; + int strip = 200; for (unsigned int mult = 1; mult <= n_multiplets; mult++) { for (unsigned int layer = 1; layer <= n_layers; layer++) { @@ -535,58 +593,44 @@ int MMT_Parameters::is_v(int plane){ return (std::find(planes_v.begin(), planes_ vector<int> MMT_Parameters::q_planes(const string& type) const{ if(type.length()!=1) throw std::runtime_error("MMT_Parameters::q_planes: Invalid type"); - if(type.compare("x")!=0&&type.compare("u")!=0&&type.compare("v")!=0){ + if(type.compare("x")!=0 && type.compare("u")!=0 && type.compare("v")!=0){ ATH_MSG_WARNING("Unsupported plane type " << type << " in q_planes...aborting....\n"); throw std::runtime_error("MMT_Parameters::q_planes: Unsupported plane type"); } vector<int> q_planes; for(unsigned int ip=0;ip<setup.size();ip++){ - if(setup.compare(ip,1,type)==0){ q_planes.push_back(ip); - } + if(setup.compare(ip,1,type)==0) q_planes.push_back(ip); } return q_planes; } par_par MMT_Parameters::param_par() const{ - return par_par(h.getFixed(),CT_x,CT_uv,uv_error.getFixed(),setup,islarge,dlm_new,genbg,chargeThreshold,misal,correct,val_tbl); + return par_par(h,CT_x,CT_uv,uv_error,setup,islarge,dlm_new,genbg,chargeThreshold,misal,correct,val_tbl); } double MMT_Parameters::y_from_eta_wedge(double eta,int plane)const{ //assumes wedge geometry--average x^2, is 1/3 y^2*tan^2(stereo_degree), for eta/y correspondence - double z=z_nominal[plane].getFixed(),zeta=TMath::DegToRad()*(0.5*stereo_degree.getFixed()); - return z*tan(2*atan(exp(-1.*eta)))/sqrt(1+tan(zeta)*tan(zeta)/3.); + double z=z_nominal[plane]; + double zeta=TMath::DegToRad()*(0.5*stereo_degree); + return z*std::tan(2*std::atan(std::exp(-1.*eta)))/std::sqrt(1+std::tan(zeta)*std::tan(zeta)/3.); } double MMT_Parameters::eta_wedge_from_y(double y,int plane)const{ - double z=z_nominal[plane].getFixed(),zeta=TMath::DegToRad()*(0.5*stereo_degree.getFixed()); - return -1.*log(tan(0.5*atan(y/z*sqrt(1+tan(zeta)*tan(zeta)/3.)))); + double z=z_nominal[plane]; + double zeta=TMath::DegToRad()*(0.5*stereo_degree); + return -1.*std::log(std::tan(0.5*std::atan(y/z*std::sqrt(1+std::tan(zeta)*std::tan(zeta)/3.)))); } int MMT_Parameters::ybin(double y,int plane)const{ - double base=ybases[plane].front().getFixed(),seg_len=h1.getFixed()/ybins; + double base=ybases[plane].front(), seg_len=h1/ybins; //if it's over, just keep it...for now--if it's below keep it, too: perhaps necessary for larger slope road sizes //update: DON'T keep the below ones....terrible things for drs int the_bin=ybases.front().size()-1; for(int i=0;i<=ybins;i++){ - if(y<base+seg_len*i)return i-1; + if(y<(base+seg_len*i)) return i-1; } return the_bin; - - double eta_min_y=eta_wedge_from_y(base+h1.getFixed(),plane),eta_max_y=eta_wedge_from_y(base,plane),eta_y=eta_wedge_from_y(y,plane); - double segment_length=(eta_max_y-eta_min_y)/ybins; - for(int i=ybins;i>=0;i--){ - if(eta_y>(eta_max_y-segment_length*i)){ - return i-1; - } - } -} - -int MMT_Parameters::ybin(float32fixed<18> y,int plane)const{ - return ybin(y.getFixed()*MMTStructConst,plane); -} -int MMT_Parameters::ybin(float32fixed<yzdex> y,int plane)const{ - return ybin(y.getFixed()*MMTStructConst,plane); } //make the local slope ab for all ybins @@ -595,11 +639,12 @@ int MMT_Parameters::ybin(float32fixed<yzdex> y,int plane)const{ void MMT_Parameters::fill_full_Ak_Bk(){ AB_k_local.clear(); - vector<int> x_planes=q_planes("x"); int nx=x_planes.size(); + vector<int> x_planes=q_planes("x"); + int nx=x_planes.size(); vector<int> variable_key(nx,-1);//a grid with all possible do{ - AB_k_local[variable_key]=ak_bk_hit_bins(variable_key); - }while(!toggle_key(variable_key,-1,ybins-1)); + AB_k_local[variable_key] = ak_bk_hit_bins(variable_key); + } while(!toggle_key(variable_key,-1,ybins-1)); } //The problem is this: we don't a priori know how many X planes there will be (configured at run time) @@ -609,11 +654,11 @@ void MMT_Parameters::fill_full_Ak_Bk(){ //We start at the rightmost entry and check values until we get a value != highest possible value and then // increment this up, resetting all values to the right of it to lo_value. bool MMT_Parameters::toggle_key(vector<int>& key,int lo_value,int hi_value)const{ - for(int idex=(int)(key.size()-1);idex>=0;idex--){ + for(int idex=(int)(key.size()-1); idex>=0; idex--){ if(key[idex]==hi_value) continue; else{ key[idex]++; - for(unsigned int reset_dex=idex+1;reset_dex<key.size();reset_dex++) key[reset_dex]=lo_value; + for(unsigned int reset_dex=idex+1; reset_dex<key.size(); reset_dex++) key[reset_dex]=lo_value; return false; } } @@ -621,8 +666,9 @@ bool MMT_Parameters::toggle_key(vector<int>& key,int lo_value,int hi_value)const } int MMT_Parameters::eta_bin(double theta) const{ - int ebin=-999; double eta=-log(tan(0.5*theta)); - if(theta==-999||theta==-16)return -1; + int ebin=-999; + double eta=-std::log(std::tan(0.5*theta)); + if(theta==-999 || theta==-16) return -1; for(int i=0;i<=n_etabins;i++){ if(eta<m_etabins[i]){ ebin=i-1; @@ -637,15 +683,16 @@ int MMT_Parameters::eta_bin(double theta) const{ } string MMT_Parameters::eta_str(int eta) const{ - ostringstream what; what<<"_eta"; - if(eta>=0&&eta<n_etabins) what<<m_etabins[eta]<<"_"<<m_etabins[eta+1]; + ostringstream what; + what<<"_eta"; + if(eta>=0 && eta<n_etabins) what << m_etabins[eta] << "_" << m_etabins[eta+1]; else what << "all"; return what.str(); } int MMT_Parameters::phi_bin(double phi) const{ int pbin=-999; - if(phi==-999||phi==-16)return n_phibins*0.5; + if(phi==-999 || phi==-16) return n_phibins*0.5; for(int i=0;i<=n_phibins;i++){ if(phi<m_phibins[i]){ pbin=i-1; @@ -660,8 +707,9 @@ int MMT_Parameters::phi_bin(double phi) const{ } string MMT_Parameters::phi_str(int phi) const{ - ostringstream what; what<<"_phi"; - if(phi>=0&&phi<n_phibins) what<<m_phibins[phi]<<"_"<<m_phibins[phi+1]; + ostringstream what; + what<<"_phi"; + if(phi>=0 && phi<n_phibins) what << m_phibins[phi] << "_" << m_phibins[phi+1]; else what << "all"; return what.str(); } @@ -673,10 +721,11 @@ string MMT_Parameters::phi_str(int phi) const{ pair<double,double> MMT_Parameters::ak_bk_hit_bins(const vector<int>& hits)const{ vector<int> x_planes=q_planes("x"); assert(hits.size()==x_planes.size()); - int nhits=0; double sum_x=0,sum_xx=0; + int nhits=0; + double sum_x=0,sum_xx=0; for(int ih=0;ih<(int)(hits.size());ih++){ - if(hits[ih]<0||hits[ih]>=ybins)continue; - double addme=z_large[hits[ih]][x_planes[ih]].getFixed()/MMTStructConst; + if(hits[ih]<0 || hits[ih]>=ybins) continue; + double addme=z_large[hits[ih]][x_planes[ih]]/MMTStructConst; sum_x += addme; sum_xx += addme*addme; nhits++; @@ -686,22 +735,22 @@ pair<double,double> MMT_Parameters::ak_bk_hit_bins(const vector<int>& hits)const } void MMT_Parameters::fill_slims(){ - Ak_local_slim=vector<vector<vector<float32fixed<zbardex> > > >(11,vector<vector<float32fixed<zbardex> > >(ybins,vector<float32fixed<zbardex> >())); - Bk_local_slim=vector<vector<vector<float32fixed<bkdex> > > >(11,vector<vector<float32fixed<bkdex> > >(ybins,vector<float32fixed<bkdex> >())); + Ak_local_slim=vector<vector<vector<double > > >(11,vector<vector<double > >(ybins,vector<double >())); + Bk_local_slim=vector<vector<vector<double > > >(11,vector<vector<double > >(ybins,vector<double >())); for(int xdex=0;xdex<11;xdex++){ vector<bool>xhits(lcl_int_to_xhits(xdex)); int ntru=0; - for(int k=0;k<(int)xhits.size();k++)ntru+=xhits[k]; - ATH_MSG_DEBUG("xdex="<<xdex<<" ("<<ntru<<" hits)"); - for(int ybin=0;ybin<ybins;ybin++){ - ATH_MSG_DEBUG("\tybin="<<ybin ); - if(ybin==ybins-1)ntru=1; - for(int which=0;which<ntru;which++){ - ATH_MSG_DEBUG("\t\twhich="<<which); + for(int k=0; k<(int)xhits.size(); k++) ntru+=xhits[k]; + ATH_MSG_DEBUG("xdex=" << xdex << " (" << ntru << " hits)"); + for(int ybin=0; ybin<ybins; ybin++){ + ATH_MSG_DEBUG("\tybin=" << ybin ); + if(ybin==ybins-1) ntru=1; + for(int which=0; which<ntru; which++){ + ATH_MSG_DEBUG("\t\twhich=" << which); vector<int>indices(indices_to_key(xdex,ybin,which)); for(int i=0;i<(int)indices.size();i++) ATH_MSG_DEBUG(indices[i]); pair<double,double>howdy(ak_bk_hit_bins(indices)); - ATH_MSG_DEBUG(setprecision(8)<<"---akbk.first="<<howdy.first<<", second="<<howdy.second); + ATH_MSG_DEBUG(setprecision(8) << "---akbk.first=" << howdy.first << ", second=" << howdy.second); Ak_local_slim[xdex][ybin].push_back(howdy.first); Bk_local_slim[xdex][ybin].push_back(howdy.second); } @@ -725,14 +774,17 @@ void MMT_Parameters::fill_crep_table(const string&dir,const string&tag){ for(int j=0;j<n_phibins;j++){ string pstr=phi_str(i); for(int k=0;k<nk;k++){ - string title;double the,phi,dth; + string title; + double the,phi,dth; crep>>title; if(title!=estr+pstr+index_to_hit_str(k)){ - ATH_MSG_WARNING("Something's wrong with your simulation-based correct read-in...you want entries for "<<estr+pstr+index_to_hit_str(k)<<", but you got "<<title<<" in "<<crep_nom.str()); - throw std::runtime_error("MMT_Parameters::fill_crep_table: Invalid configuration"); + ATH_MSG_WARNING("Something's wrong with your simulation-based correct read-in...you want entries for " << estr+pstr+index_to_hit_str(k) << ", but you got "<< title << " in " << crep_nom.str()); + throw std::runtime_error("MMT_Parameters::fill_crep_table: Invalid configuration"); } - crep>>the>>phi>>dth; - crep_table[i][j][k][0]=the*fudge_factor;crep_table[i][j][k][1]=phi;crep_table[i][j][k][2]=dth*fudge_factor; + crep >> the >> phi >> dth; + crep_table[i][j][k][0]=the*fudge_factor; + crep_table[i][j][k][1]=phi; + crep_table[i][j][k][2]=dth*fudge_factor; } } } @@ -742,22 +794,22 @@ void MMT_Parameters::fill_crep_table(const string&dir,const string&tag){ void MMT_Parameters::fill_yzmod(){ vector<int> x_planes=q_planes("x");//this tells us which planes are x planes int nxp=x_planes.size(); - ymod=vector<vector<vector<float32fixed<yzdex> > > >(n_etabins,vector<vector<float32fixed<yzdex> > >(n_phibins,vector<float32fixed<yzdex> >(nxp,float32fixed<yzdex>(0.)))); - zmod=vector<vector<vector<float32fixed<yzdex> > > >(n_etabins,vector<vector<float32fixed<yzdex> > >(n_phibins,vector<float32fixed<yzdex> >(nxp,float32fixed<yzdex>(0.)))); - double pbump=0.5*(maximum_large_phi.getFixed()-minimum_large_phi.getFixed())/(n_phibins); + ymod=vector<vector<vector<double > > >(n_etabins,vector<vector<double > >(n_phibins,vector<double >(nxp,double(0.)))); + zmod=vector<vector<vector<double > > >(n_etabins,vector<vector<double > >(n_phibins,vector<double >(nxp,double(0.)))); + double pbump=0.5*(maximum_large_phi-minimum_large_phi)/(n_phibins); for(int et=0;et<n_etabins;et++){ - double theta=2*atan(exp(-1.*m_etabins[et])); + double theta=2*std::atan(std::exp(-1.*m_etabins[et])); for(int ph=0;ph<n_phibins;ph++){ double phi=m_phibins[ph]+pbump; for(int pl=0;pl<nxp;pl++){ - if(correct.type!=2||x_planes[pl]>3)continue;//if we don't correct or if it's the back multiplet, just leave these as zero + if(correct.type!=2 || x_planes[pl]>3) continue;//if we don't correct or if it's the back multiplet, just leave these as zero double yadd=0.; double zadd=0.; - double zflt=z_nominal[x_planes[pl]].getFixed(); - double x=zflt*tan(theta)*sin(phi); - double yflt=zflt*tan(theta)*cos(phi); - double yup=yflt-ybases[x_planes[pl]][0].getFixed(); + double zflt=z_nominal[x_planes[pl]]; + double x=zflt*std::tan(theta)*std::sin(phi); + double yflt=zflt*std::tan(theta)*std::cos(phi); + double yup=yflt-ybases[x_planes[pl]][0]; double alpha=correct.rotate.Z(); double beta=correct.rotate.Y(); @@ -767,14 +819,15 @@ void MMT_Parameters::fill_yzmod(){ //beta angle--rotation around the y axis; the x coordinate or phi taken care of in correction to slope road limits in x in constructor //z is changed, and so y must be scaled (should it?): THIS DOESN'T WORK.... // zadd-=1.*x*sin(beta)/MMTStructConst;yadd+=yup*zadd/zflt/MMTStructConst; - zadd-=tan(beta)*tan(theta)*sin(phi)*zflt/MMTStructConst; - yadd-=tan(beta)*tan(theta)*sin(phi)*yflt/MMTStructConst; + zadd-=std::tan(beta)*std::tan(theta)*std::sin(phi)*zflt/MMTStructConst; + yadd-=std::tan(beta)*std::tan(theta)*std::sin(phi)*yflt/MMTStructConst; //alpha angle--rotation around the z axis; z unchanged, y the expected rotation (final rotation taken care of at start) - yadd+=((cos(alpha)-1.)*yup+x*sin(alpha))/MMTStructConst; + yadd+=((std::cos(alpha)-1.)*yup+x*std::sin(alpha))/MMTStructConst; //add the entries - ymod[et][ph][pl]=float32fixed<yzdex>(yadd);zmod[et][ph][pl]=float32fixed<yzdex>(zadd); + ymod[et][ph][pl]=double(yadd); + zmod[et][ph][pl]=double(zadd); } } } @@ -784,15 +837,16 @@ void MMT_Parameters::index_key_test(){ for(int xdex=0;xdex<11;xdex++){ vector<bool>xhits(lcl_int_to_xhits(xdex)); int ntru=0; - for(int k=0;k<(int)xhits.size();k++)ntru+=xhits[k]; - for(int ybin=0;ybin<ybins;ybin++){ - if(ybin==ybins-1)ntru=1; - for(int which=0;which<ntru;which++){ - ATH_MSG_DEBUG("Let's start with ["<<xdex<<"]["<<ybin<<"]["<<which<<"] makes "); + for(int k=0; k<(int)xhits.size(); k++) ntru+=xhits[k]; + for(int ybin=0; ybin<ybins; ybin++){ + if(ybin==ybins-1) ntru=1; + for(int which=0; which<ntru; which++){ + ATH_MSG_DEBUG("Let's start with [" << xdex << "][" << ybin << "][" << which << "] makes "); vector<int>key(indices_to_key(xdex,ybin,which)); - for(int i=0;i<(int)key.size();i++) ATH_MSG_DEBUG(key[i]<<" "); - int x,y,w; key_to_indices(key,x,y,w); - ATH_MSG_DEBUG("...and back to ["<<x<<"]["<<y<<"]["<<w<<"] "); + for(int i=0; i<(int)key.size(); i++) ATH_MSG_DEBUG(key[i] << " "); + int x,y,w; + key_to_indices(key,x,y,w); + ATH_MSG_DEBUG("...and back to [" << x << "][" << y << "][" << w << "] "); } } } @@ -802,26 +856,29 @@ void MMT_Parameters::index_key_test(){ void MMT_Parameters::key_to_indices(const vector<int>& key,int& xdex,int& ybin,int& which)const{ vector<bool> boogah;vector<int> was_hit; for(unsigned int i=0;i<key.size();i++){ - boogah.push_back(key[i]>=0&&key[i]<ybins); - if(boogah.back())was_hit.push_back(key[i]); + boogah.push_back(key[i]>=0 && key[i]<ybins); + if(boogah.back()) was_hit.push_back(key[i]); } - bool even=true;int dev_pos=0,dev_bin=-1; + bool even=true; + int dev_pos=0, dev_bin=-1; for(unsigned int i=0;i<was_hit.size();i++){ if(was_hit[i]!=was_hit.front()){ if(even){ dev_pos=i; if(was_hit[i]!=was_hit.front()+1){ - dev_pos=-2;break; + dev_pos=-2; + break; } else dev_bin=was_hit[i]; } else if(was_hit[i]!=dev_bin){ - dev_pos=-2;break; + dev_pos=-2; + break; } even=false; } } - if(dev_pos<0||dev_pos>=(int)boogah.size())which=0;//if the track isn't upward going, make the single bin assumption + if(dev_pos<0 || dev_pos>=(int)boogah.size()) which=0;//if the track isn't upward going, make the single bin assumption else which=dev_pos; xdex=xhits_to_lcl_int(boogah); if(was_hit.size()==0) was_hit.push_back(-999); @@ -831,11 +888,12 @@ void MMT_Parameters::key_to_indices(const vector<int>& key,int& xdex,int& ybin,i vector<int> MMT_Parameters::indices_to_key(int xdex,int ybin,int which)const{ vector<bool> xhits=lcl_int_to_xhits(xdex); vector<int> key(xhits.size(),-1); - int count=0;int raise_it=(which==0?xhits.size():which); + int count=0; + int raise_it=(which==0?xhits.size():which); for(unsigned int ix=0;ix<xhits.size();ix++){ if(xhits[ix]){ key[ix]=ybin; - if(count>=raise_it)key[ix]++; + if(count>=raise_it) key[ix]++; count++; } } @@ -850,17 +908,17 @@ void MMT_Parameters::Local_Slope_A_B(){ void MMT_Parameters::Slope_Components_ROI(){ //these are mm/mm i.e. not units of strip # - Slope_to_ROI=vector<vector<vector<float32fixed<4> > > >(n_y,vector<vector<float32fixed<4> > >(n_x,vector<float32fixed<4> >(2,float32fixed<4>(0.)))); - for(int ix=0;ix<n_x;ix++){ - for(int jy=0;jy<n_y;jy++){ - int xdex=ix,ydex=jy+1; - if(xdex==0)xdex++; + Slope_to_ROI=vector<vector<vector<double > > >(n_y,vector<vector<double > >(n_x,vector<double >(2,double(0.)))); + for(int ix=0; ix<n_x; ix++){ + for(int jy=0; jy<n_y; jy++){ + int xdex=ix, ydex=jy+1; + if(xdex==0) xdex++; - double M_x = ( m_x_min + h_mx * xdex ).getFixed(); - double M_y = ( m_y_min + h_my * ydex ).getFixed(); + double M_x = ( m_x_min + h_mx * xdex ); + double M_y = ( m_y_min + h_my * ydex ); - double theta = atan(sqrt(M_x*M_x+M_y*M_y)); - double phi=atan2(M_y,M_x); + double theta = std::atan(std::sqrt(M_x*M_x+M_y*M_y)); + double phi= std::atan2(M_y,M_x); if(minimum_large_phi>phi || maximum_large_phi<phi) phi=999; if(minimum_large_theta>theta || maximum_large_theta<theta) theta=0; Slope_to_ROI[jy][ix][0] = theta; @@ -879,9 +937,9 @@ void MMT_Parameters::Delta_theta_optimization_LG(){ double LG_region_width = (LG_max - LG_min)/number_LG_regions;//LG_min=0 means e no neg slopes for(int ilgr=0; ilgr<number_LG_regions; ilgr++){ - vector<float32fixed<2> > tmpVector; + vector<double > tmpVector; tmpVector.push_back(LG_min+LG_region_width *ilgr);//LG - tmpVector.push_back( float32fixed<2>(1.) / (tmpVector.front()/a + a));//mult_factor + tmpVector.push_back( double(1.) / (tmpVector.front()/a + a));//mult_factor DT_Factors.push_back(tmpVector); } ATH_MSG_DEBUG( "return Delta_theta_optimization_LG" ); @@ -889,40 +947,40 @@ void MMT_Parameters::Delta_theta_optimization_LG(){ int MMT_Parameters::xhits_to_lcl_int(const vector<bool>& xhits) const{ if(xhits.size()!=4){ - ATH_MSG_WARNING("There should be 4 xplanes, only "<<xhits.size()<<" in the xhit vector given to local_slope_index()\n"); + ATH_MSG_WARNING("There should be 4 xplanes, only " << xhits.size() << " in the xhit vector given to local_slope_index()\n"); throw std::runtime_error("MMT_Parameters::xhits_to_lcl_int: Invalid number of planes"); } - if(xhits[0]&& xhits[1]&& xhits[2]&& xhits[3]) return 0; - else if( xhits[0]&& xhits[1]&& xhits[2]&&!xhits[3]) return 1; - else if( xhits[0]&& xhits[1]&&!xhits[2]&& xhits[3]) return 2; - else if( xhits[0]&&!xhits[1]&& xhits[2]&& xhits[3]) return 3; - else if(!xhits[0]&& xhits[1]&& xhits[2]&& xhits[3]) return 4; - else if( xhits[0]&& xhits[1]&&!xhits[2]&&!xhits[3]) return 5; - else if( xhits[0]&&!xhits[1]&& xhits[2]&&!xhits[3]) return 6; - else if( xhits[0]&&!xhits[1]&&!xhits[2]&& xhits[3]) return 7; - else if(!xhits[0]&& xhits[1]&& xhits[2]&&!xhits[3]) return 8; - else if(!xhits[0]&& xhits[1]&&!xhits[2]&& xhits[3]) return 9; - else if(!xhits[0]&&!xhits[1]&& xhits[2]&& xhits[3]) return 10; + if(xhits[0] && xhits[1] && xhits[2] && xhits[3]) return 0; + else if( xhits[0] && xhits[1] && xhits[2] && !xhits[3]) return 1; + else if( xhits[0] && xhits[1] && !xhits[2] && xhits[3]) return 2; + else if( xhits[0] && !xhits[1] && xhits[2] && xhits[3]) return 3; + else if(!xhits[0] && xhits[1] && xhits[2] && xhits[3]) return 4; + else if( xhits[0] && xhits[1] && !xhits[2] && !xhits[3]) return 5; + else if( xhits[0] && !xhits[1] && xhits[2] && !xhits[3]) return 6; + else if( xhits[0] && !xhits[1] && !xhits[2] && xhits[3]) return 7; + else if(!xhits[0] && xhits[1] && xhits[2] && !xhits[3]) return 8; + else if(!xhits[0] && xhits[1] && !xhits[2] && xhits[3]) return 9; + else if(!xhits[0] && !xhits[1] && xhits[2] && xhits[3]) return 10; return -999; } int MMT_Parameters::nsimmax_1d()const{ - return pow(2.,setup.size()); + return std::pow(2.,setup.size()); } int MMT_Parameters::bool_to_index(const vector<bool>&track)const{ assert(track.size()==setup.size()); int index=0; - for(int plane=0;plane<(int)track.size();plane++){ - if(!track[plane])continue; - index = index | int(pow(2.,plane)); + for(int plane=0; plane<(int)track.size(); plane++){ + if(!track[plane]) continue; + index = index | int(std::pow(2.,plane)); } return index; } vector<bool> MMT_Parameters::index_to_bool(int index)const{ vector<bool>code(setup.size(),false); - for(int i=0;i<(int)code.size();i++)code[i]=(index&int(pow(2.,i))); + for(int i=0; i<(int)code.size(); i++) code[i]=(index&int(std::pow(2.,i))); return code; } @@ -931,34 +989,36 @@ string MMT_Parameters::index_to_hit_str(int index)const{ } string MMT_Parameters::bool_to_hit_str(const vector<bool>&track)const{ - ostringstream ok;ok<<"_ht"; - for(int i=0;i<(int)track.size();i++)ok<<track[i]; + ostringstream ok; + ok<<"_ht"; + for(int i=0; i<(int)track.size(); i++) ok << track[i]; return ok.str(); } vector<bool> MMT_Parameters::lcl_int_to_xhits(int lcl_int)const{ vector<bool> xhits(4,true); - if(lcl_int<0||lcl_int>10){ - ATH_MSG_WARNING("Wherefore dost thou chooseth the hits of planes of X for thy lcl_int "<<lcl_int<<"? 'Tis not in [0,10]!"); + if(lcl_int<0 || lcl_int>10){ + ATH_MSG_WARNING("Wherefore dost thou chooseth the hits of planes of X for thy lcl_int " << lcl_int << "? 'Tis not in [0,10]!"); throw std::runtime_error("MMT_Parameters::lcl_int_to_xhits: invalid value for lcl_int"); } if(lcl_int==0) return xhits; - if(lcl_int==1||lcl_int==5||lcl_int==6||lcl_int==8)xhits[3]=false; - if(lcl_int==2||lcl_int==5||lcl_int==7||lcl_int==9)xhits[2]=false; - if(lcl_int==3||lcl_int==6||lcl_int==7||lcl_int==10)xhits[1]=false; - if(lcl_int==4||lcl_int==8||lcl_int==9||lcl_int==10)xhits[0]=false; + if(lcl_int==1 || lcl_int==5 || lcl_int==6 || lcl_int==8) xhits[3]=false; + if(lcl_int==2 || lcl_int==5 || lcl_int==7 || lcl_int==9) xhits[2]=false; + if(lcl_int==3 || lcl_int==6 || lcl_int==7 || lcl_int==10) xhits[1]=false; + if(lcl_int==4 || lcl_int==8 || lcl_int==9 || lcl_int==10) xhits[0]=false; return xhits; } string MMT_Parameters::lcl_int_to_xhit_str(int lcl_int)const{ vector<bool>hits=lcl_int_to_xhits(lcl_int); - ostringstream out;out<<"_xht"; - for(int i=0;i<(int)hits.size();i++)out<<hits[i]; + ostringstream out; + out<<"_xht"; + for(int i=0;i<(int)hits.size();i++) out<<hits[i]; return out.str(); } evInf_entry::evInf_entry(int event,int pdg,double e,double p,double ieta,double peta,double eeta,double iphi,double pphi,double ephi,double ithe,double pthe,double ethe,double dth, - int trn,int mun,const TVector3& tex,int troi,int antev,int postv,int nxh,int nuvh,int nxbg,int nuvbg,double adt,int difx, + int trn,int mun,const ROOT::Math::XYZVector& tex,int troi,int antev,int postv,int nxh,int nuvh,int nxbg,int nuvbg,double adt,int difx, int difuv,bool cut,bool bad) :athena_event(event),pdg_id(pdg),E(e),pt(p),eta_ip(ieta),eta_pos(peta),eta_ent(eeta),phi_ip(iphi),phi_pos(pphi),phi_ent(ephi),theta_ip(ithe),theta_pos(pthe),theta_ent(ethe), dtheta(dth),truth_n(trn),mu_n(mun),vertex(tex),truth_roi(troi),N_hits_preVMM(antev),N_hits_postVMM(postv),N_X_hits(nxh),N_UV_hits(nuvh),NX_bg_preVMM(nxbg), @@ -967,7 +1027,7 @@ evInf_entry::evInf_entry(int event,int pdg,double e,double p,double ieta,double void evInf_entry::print() const{ } -evFit_entry::evFit_entry(int event,float32fixed<4> fthe,float32fixed<4> fphi,float32fixed<2> fdth,int roi,int xhit,int uvhit,int bgx,int bguv,float32fixed<2> dth_nd,int hc,int tph,int bgph) +evFit_entry::evFit_entry(int event,double fthe,double fphi,double fdth,int roi,int xhit,int uvhit,int bgx,int bguv,double dth_nd,int hc,int tph,int bgph) :athena_event(event),fit_theta(fthe),fit_phi(fphi),fit_dtheta(fdth),fit_roi(roi),X_hits_in_fit(xhit),UV_hits_in_fit(uvhit),bg_X_fit(bgx), bg_UV_fit(bguv),dtheta_nodiv(dth_nd),hcode(hc),truth_planes_hit(tph),bg_planes_hit(bgph) {} @@ -977,7 +1037,7 @@ void evFit_entry::print()const{ hitData_key::hitData_key(int bct, double t, double gt, int vmm, int ev):BC_time(bct),time(t),gtime(gt),VMM_chip(vmm),event(ev) {} bool hitData_key::operator==(const hitData_key& rhs) const{ - if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime==rhs.gtime&&this->VMM_chip==rhs.VMM_chip&&this->event==rhs.event) return true; + if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime==rhs.gtime && this->VMM_chip==rhs.VMM_chip && this->event==rhs.event) return true; return false; } @@ -986,19 +1046,19 @@ bool hitData_key::operator!=(const hitData_key& rhs) const{ } bool hitData_key::operator<(const hitData_key& rhs) const{ if(this->BC_time<rhs.BC_time) return true; - else if(this->BC_time==rhs.BC_time&&this->time<rhs.time) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime<rhs.gtime) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime==rhs.gtime&&this->VMM_chip<rhs.VMM_chip) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime==rhs.gtime&&this->VMM_chip==rhs.VMM_chip&&this->event<rhs.event) return true; + else if(this->BC_time==rhs.BC_time && this->time<rhs.time) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime<rhs.gtime) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime==rhs.gtime && this->VMM_chip<rhs.VMM_chip) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime==rhs.gtime && this->VMM_chip==rhs.VMM_chip && this->event<rhs.event) return true; return false; } bool hitData_key::operator>(const hitData_key& rhs) const{ if(this->BC_time>rhs.BC_time) return true; - else if(this->BC_time==rhs.BC_time&&this->time>rhs.time) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime>rhs.gtime) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime==rhs.gtime&&this->VMM_chip>rhs.VMM_chip) return true; - else if(this->BC_time==rhs.BC_time&&this->time==rhs.time&&this->gtime==rhs.gtime&&this->VMM_chip==rhs.VMM_chip&&this->event>rhs.event) return true; + else if(this->BC_time==rhs.BC_time && this->time>rhs.time) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime>rhs.gtime) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime==rhs.gtime && this->VMM_chip>rhs.VMM_chip) return true; + else if(this->BC_time==rhs.BC_time && this->time==rhs.time && this->gtime==rhs.gtime && this->VMM_chip==rhs.VMM_chip && this->event>rhs.event) return true; return false; } @@ -1012,39 +1072,41 @@ bool hitData_key::operator>=(const hitData_key& rhs) const{ string hitData_key::hdr()const{ ostringstream out; - out<<setw(9)<<"BC_t"<<setw(9)<<"t"<<setw(9)<<"g_t"<<setw(9)<<"VMM"<<setw(9)<<"event"; + out << setw(12) << "BC_t" << setw(12) << "t" << setw(12) << "g_t" << setw(12) << "VMM" << setw(12) << "event"; return out.str(); } string hitData_key::str()const{ ostringstream out; - out<<setw(9)<<this->BC_time<<setw(9)<<this->time<<setw(9)<<this->gtime<<setw(9)<<this->VMM_chip<<setw(9)<<this->event; + out << setw(12) << this->BC_time << setw(12) << this->time << setw(12) << this->gtime << setw(12) << this->VMM_chip << setw(12) << this->event; return out.str(); } void hitData_key::print()const{ } -hitData_info::hitData_info(int pl,int station_eta,int strip,MMT_Parameters *par,const TVector3&tru,double tpos,double ppos):plane(pl){ +hitData_info::hitData_info(int pl,int station_eta,int strip,MMT_Parameters *par,const ROOT::Math::XYZVector &tru,double tpos,double ppos):plane(pl){ (void) tru; //The idea here is to calculate/assign a y and a z to a given hit based on its pl/station/strip, the geometry of the detector (in par), and misalignment based on position. //We start by assigning the plane dependent strip width (the stereo planes come in skew and so get divided by cos(stereo_angle) char schar=par->setup[plane]; bool horizontal=(schar=='x'||schar=='X'); - double swidth=par->strip_width.getFixed(); - double base=par->ybases[plane][station_eta-1].getFixed(); - if(!horizontal)swidth/=cos(TMath::DegToRad()*(par->stereo_degree.getFixed())); + double swidth=par->strip_width; + int eta = TMath::Abs(station_eta) -1; + double base=par->ybases[plane][eta]; + if(!horizontal) swidth/=std::cos(TMath::DegToRad()*(par->stereo_degree)); //Next, we initialize some constants--y will eventually be calculated as y=base+scale*(width*strip+delta_y(misalignment,correction)). //The correction portion of delta_y is done in in the calculations on the ybases object in par //yup, or "y up" is the portion of y above the base of the plane (i.e. in the detector) - double delta_y=mis_dy(pl,par,tpos,ppos),yup_scale=1.; + double delta_y=mis_dy(pl,par,tpos,ppos), yup_scale=1.; double yflt=-999; yflt=base+(strip*swidth+delta_y)*yup_scale; //Note that ybin should be calculated on the basis of a misaligned y, not a corrected one (to get the position on the wedge just right), //but the difference for this part of the calculation should be negligible (an effect of <~ misal/(bin length) ~ (5 mm/(3600 mm/15)) ~ %. int bin=par->ybin(yflt,plane); + if(bin==-1) bin=0; //We have here both the "perfect" and binned z calculations/lookups; the commented out one is there for reference - double zflt=par->z_large[bin][plane].getFixed(); + double zflt=par->z_large[bin][plane]; // if(par->correct.type==2&&plane<4&&par->correct.rotate.X()!=0){ // // double fltyup=(base+(strip*swidth+delta_y)-planebase); // double angle=par->correct.rotate.X(); @@ -1053,69 +1115,82 @@ hitData_info::hitData_info(int pl,int station_eta,int strip,MMT_Parameters *par, // // } // y = yflt / MMTStructConst; // z = zflt / MMTStructConst; - + y = yflt; + z = zflt; slope = (zflt!=0.) ? yflt / zflt : 0.; } double hitData_info::mis_dy(int plane,MMT_Parameters *par,double tpos,double ppos)const{ - if(par->misal.type!=1||plane>3)return 0.; - double zplane=par->z_nominal[plane].getFixed(); - double base=par->ybases[plane].front().getFixed(); + if(par->misal.type!=1 || plane>3) return 0.; + double zplane=par->z_nominal[plane]; + double base=par->ybases[plane].front(); // double s_z0=zplane; - double s_x0=zplane*tan(tpos)*sin(ppos); - double s_y0=zplane*tan(tpos)*cos(ppos);//initial position + double s_x0=zplane*std::tan(tpos)*std::sin(ppos); + double s_y0=zplane*std::tan(tpos)*std::cos(ppos);//initial position - double hats_z0=cos(tpos),hats_x0=sin(tpos)*sin(ppos),hats_y0=sin(tpos)*cos(ppos);//muon track unit vector + double hats_z0=std::cos(tpos); + double hats_x0=std::sin(tpos)*std::sin(ppos); + double hats_y0=std::sin(tpos)*std::cos(ppos);//muon track unit vector double zeta_y0=s_y0-base;//height in y in the wedge local coordinates--this is what we have to compare in the end - double alpha=par->misal.rotate.Z(),beta=par->misal.rotate.Y(),gamma=par->misal.rotate.X();//rotation angles - double ds=par->misal.translate.X(),dz=par->misal.translate.Y(),dt=-1.*par->misal.translate.Z();//t comes in -z - double O_bxf=ds,O_byf=base+dz,O_bzf=zplane+dt;//position of bottom of the wedge in global coordinates; subtract this from the final y position (s_yf) for final zeta (comparison to add to y position) - double yhat_x=-1.*sin(alpha)*cos(beta),yhat_y=(cos(alpha)*cos(gamma)-sin(alpha)*sin(beta)*sin(gamma)),yhat_z=(cos(alpha)*sin(gamma)+sin(alpha)*sin(beta)*cos(gamma));//new y direction after rotations; horizontal case + double alpha=par->misal.rotate.Z(); + double beta=par->misal.rotate.Y(); + double gamma=par->misal.rotate.X();//rotation angles + double ds=par->misal.translate.X(); + double dz=par->misal.translate.Y(); + double dt=-1.*par->misal.translate.Z();//t comes in -z + double O_bxf=ds; + double O_byf=base+dz; + double O_bzf=zplane+dt;//position of bottom of the wedge in global coordinates; subtract this from the final y position (s_yf) for final zeta (comparison to add to y position) + double yhat_x=-1.*std::sin(alpha)*std::cos(beta); + double yhat_y=(std::cos(alpha)*std::cos(gamma)-std::sin(alpha)*std::sin(beta)*std::sin(gamma)); + double yhat_z=(std::cos(alpha)*std::sin(gamma)+std::sin(alpha)*std::sin(beta)*std::cos(gamma));//new y direction after rotations; horizontal case char schar=par->setup[plane]; if(!(schar=='x'||schar=='X')){ //if we're in a stereo plane, calculate different coefficients. - double omega=TMath::DegToRad()*(par->stereo_degree.getFixed()); + double omega=TMath::DegToRad()*(par->stereo_degree); double pm=(schar=='u'||schar=='U'?1.:-1.); - yhat_x=pm*cos(alpha)*cos(beta)*sin(omega)-sin(alpha)*cos(beta)*cos(omega); - yhat_y=pm*sin(omega)*(sin(alpha)*cos(gamma)+cos(alpha)*sin(beta)*sin(gamma))+cos(omega)*(cos(alpha)*cos(gamma)-sin(alpha)*sin(beta)*sin(gamma)); - yhat_z=pm*sin(omega)*(sin(alpha)*sin(gamma)-cos(alpha)*sin(beta)*cos(gamma))+cos(omega)*(cos(alpha)*sin(gamma)+sin(alpha)*sin(beta)*cos(gamma)); - zeta_y0=pm*sin(omega)*s_x0+cos(omega)*(s_y0-base); + yhat_x=pm*std::cos(alpha)*std::cos(beta)*std::sin(omega)-std::sin(alpha)*std::cos(beta)*std::cos(omega); + yhat_y=pm*std::sin(omega)*(std::sin(alpha)*std::cos(gamma)+std::cos(alpha)*std::sin(beta)*std::sin(gamma))+std::cos(omega)*(std::cos(alpha)*std::cos(gamma)-std::sin(alpha)*std::sin(beta)*std::sin(gamma)); + yhat_z=pm*std::sin(omega)*(std::sin(alpha)*std::sin(gamma)-std::cos(alpha)*std::sin(beta)*std::cos(gamma))+std::cos(omega)*(std::cos(alpha)*std::sin(gamma)+std::sin(alpha)*std::sin(beta)*std::cos(gamma)); + zeta_y0=pm*std::sin(omega)*s_x0+std::cos(omega)*(s_y0-base); } - double kprime=(sin(beta)*O_bxf-cos(beta)*sin(gamma)*O_byf+cos(beta)*cos(gamma)*O_bzf)/(sin(beta)*hats_x0-cos(beta)*sin(gamma)*hats_y0+cos(beta)*cos(gamma)*hats_z0); - double zeta_xf=kprime*hats_x0-O_bxf,zeta_yf=kprime*hats_y0-O_byf,zeta_zf=kprime*hats_z0-O_bzf; + double kprime=(std::sin(beta)*O_bxf-std::cos(beta)*std::sin(gamma)*O_byf+std::cos(beta)*std::cos(gamma)*O_bzf)/(std::sin(beta)*hats_x0-std::cos(beta)*std::sin(gamma)*hats_y0+std::cos(beta)*std::cos(gamma)*hats_z0); + double zeta_xf=kprime*hats_x0-O_bxf; + double zeta_yf=kprime*hats_y0-O_byf; + double zeta_zf=kprime*hats_z0-O_bzf; double zetayf_yhatf=zeta_xf*yhat_x+zeta_yf*yhat_y+zeta_zf*yhat_z; return zetayf_yhatf-zeta_y0; } hitData_info::hitData_info(int the_pl,double the_y,double the_z) : plane(the_pl),y(the_y),z(the_z){ - if(the_z==0||the_z==-999)slope=-999; + if(the_z==0 || the_z==-999) slope=-999; else slope=the_y/the_z; } string hitData_info::hdr()const{ ostringstream out; - out<<setw(9)<<"plane"<<setw(9)<<"y"<<setw(9)<<"z"<<setw(9)<<"slope"; + out << setw(12) << "plane" << setw(12) << "y" << setw(12) << "z" << setw(12) << "slope"; return out.str(); } string hitData_info::str()const{ ostringstream out; - out<<setprecision(4)<<setw(9)<<plane<<setw(9)<<y.getFixed()<<setw(9)<<z.getFixed()<<setw(9)<<slope.getFixed(); + out << setprecision(4) << setw(12) << plane << setw(12) << y << setw(12) << z << setw(12) << slope; return out.str(); } void hitData_info::print()const{ } bool hitData_info::operator==(const hitData_info& rhs) const{ - if(this->plane==rhs.plane&&this->y==rhs.y&&this->z==rhs.z&&this->slope==rhs.slope)return true; + if(this->plane==rhs.plane && this->y==rhs.y && this->z==rhs.z && this->slope==rhs.slope)return true; return false; } Hit::Hit(const hitData_key& k,const hitData_info&i):key(k),info(i) {} bool Hit::operator==(const Hit& rhs) const{ - if(this->key==rhs.key&&this->info==rhs.info)return true; + if(this->key==rhs.key && this->info==rhs.info)return true; return false; } @@ -1126,12 +1201,11 @@ void Hit::print_track(const vector<Hit>& track) const{ void Hit::print() const{ } -hitData_entry::hitData_entry(int ev, double gt, double q, int vmm, int pl, int st, int est, double tr_the, double tru_phi, - bool q_tbg, int bct, double t, const TVector3& tru, const TVector3& rec, +hitData_entry::hitData_entry(int ev, double gt, double q, int vmm, int mmfe, int pl, int st, int est, int phi, int mult, int gg, double locX, double tr_the, double tru_phi, + bool q_tbg, int bct, double t, const ROOT::Math::XYZVector& tru, const ROOT::Math::XYZVector& rec, double fit_the, double fit_ph, double fit_dth, double tru_dth,// double tru_thl, double tru_thg, double mxg, double mug, double mvg, double mxl, double the_mx, double the_my, int the_roi): - event(ev),gtime(gt),charge(q),VMM_chip(vmm),plane(pl),strip(st),station_eta(est),tru_theta_ip(tr_the),tru_phi_ip(tru_phi),truth_nbg(q_tbg), - BC_time(bct),time(t),truth(tru),recon(rec),fit_theta(fit_the),fit_phi(fit_ph),fit_dtheta(fit_dth),tru_dtheta(tru_dth), + event(ev),gtime(gt),charge(q),VMM_chip(vmm),MMFE_VMM(mmfe),plane(pl),strip(st),station_eta(est),station_phi(phi),multiplet(mult),gasgap(gg),localX(locX),tru_theta_ip(tr_the),tru_phi_ip(tru_phi),truth_nbg(q_tbg),BC_time(bct),time(t),truth(tru),recon(rec),fit_theta(fit_the),fit_phi(fit_ph),fit_dtheta(fit_dth),tru_dtheta(tru_dth), /*tru_theta_local(tru_thl),tru_theta_global(tru_thg),*/M_x_global(mxg),M_u_global(mug),M_v_global(mvg),M_x_local(mxl),mx(the_mx),my(the_my),roi(the_roi) {} Hit hitData_entry::entry_hit(MMT_Parameters *par)const{ @@ -1143,11 +1217,19 @@ hitData_key hitData_entry::entry_key() const{ hitData_info hitData_entry::entry_info(MMT_Parameters *par)const{ hitData_info spade(plane,station_eta,strip,par,recon,tru_theta_ip,tru_phi_ip);//truth or recon? doesn't matter too much--it's for misalignment -// spade.y=recon.Y(); return spade; } -void hitData_entry::fit_fill(float32fixed<4> fthe,float32fixed<4> fphi, float32fixed<2> fdth, float32fixed<2> mxg, float32fixed<2> mug, float32fixed<2> mvg, float32fixed<2> mxl, float32fixed<2> M_x, float32fixed<2> M_y, int king){ - this->fit_theta=fthe; this->fit_phi=fphi; this->fit_dtheta=fdth; this->M_x_global=mxg; this->M_u_global=mug; this->M_v_global=mvg; this->M_x_local=mxl; this->mx=M_x; this->my=M_y; this->roi=king; +void hitData_entry::fit_fill(double fthe,double fphi, double fdth, double mxg, double mug, double mvg, double mxl, double M_x, double M_y, int king){ + this->fit_theta=fthe; + this->fit_phi=fphi; + this->fit_dtheta=fdth; + this->M_x_global=mxg; + this->M_u_global=mug; + this->M_v_global=mvg; + this->M_x_local=mxl; + this->mx=M_x; + this->my=M_y; + this->roi=king; } void hitData_entry::print() const{ @@ -1171,7 +1253,7 @@ finder_entry::finder_entry(bool the_is_hit, int the_clock,const Hit& the_k): bool finder_entry::operator==(const finder_entry& merp) const{ - if(merp.is_hit==this->is_hit&&merp.clock==this->clock&&this->hit==merp.hit) return true; + if(merp.is_hit==this->is_hit && merp.clock==this->clock && this->hit==merp.hit) return true; else return false; } bool finder_entry::operator!=(const finder_entry& merp) const{ @@ -1182,20 +1264,22 @@ ROI::ROI(double the_theta, double the_phi, double the_m_x, double the_m_y, int t theta(the_theta), phi(the_phi), m_x(the_m_x), m_y(the_m_y), roi(the_roi) {} -athena_header::athena_header(const TLorentzVector& par, int tpn, double etp, double ete, double php, double phe, int mun, const TVector3& ver): +athena_header::athena_header(const ROOT::Math::PtEtaPhiEVector& par, int tpn, double etp, double ete, double php, double phe, int mun, const ROOT::Math::XYZVector& ver): the_part(par),trupart_n(tpn),etapos(etp),etaent(ete),phipos(php),phient(phe),muent_n(mun),vertex(ver) {} digitWrapper::digitWrapper(const MmDigit* digit, + const std::string &stationName, double tmpGTime, - const TVector3& truthLPos, - const TVector3& stripLPos, - const TVector3& stripGPos + const ROOT::Math::XYZVector& truthLPos, + const ROOT::Math::XYZVector& stripLPos, + const ROOT::Math::XYZVector& stripGPos ): digit(digit), + stName(stationName), gTime(tmpGTime), truth_lpos(truthLPos), strip_lpos(stripLPos), strip_gpos(stripGPos){} -track_address::track_address(int bct,bool big,int wed,int pl,int sh,const TVector3& chr): +track_address::track_address(int bct,bool big,int wed,int pl,int sh,const ROOT::Math::XYZVector& chr): BC(bct),islarge(big),wedge(wed),plane(pl),strip_hit(sh),cart_hit(chr) {} diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerTool.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerTool.cxx index 2bfd7b2a662a8da71df9feafc381c7348299b204..c19b9565be82f701a3db38d8f3c06f6edcc7f945 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerTool.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerTool.cxx @@ -9,6 +9,7 @@ // Athena/Gaudi includes #include "GaudiKernel/ITHistSvc.h" #include "GaudiKernel/IIncidentSvc.h" +#include "MuonIdHelpers/MmIdHelper.h" //Muon software includes #include "MuonDigitContainer/MmDigit.h" @@ -45,6 +46,7 @@ namespace NSWL1 { declareProperty("MM_DigitContainerName", m_MmDigitContainer = "MM_DIGITS", "the name of the MM digit container"); declareProperty("DoNtuple", m_doNtuple = true, "input the MMStripTds branches into the analysis ntuple"); + diamond = nullptr; } MMTriggerTool::~MMTriggerTool() { @@ -117,278 +119,329 @@ namespace NSWL1 { return StatusCode::SUCCESS; } + StatusCode MMTriggerTool::initDiamondAlgorithm() { + diamond = new MMT_Diamond(m_detManager); + diamond->setTrapezoidalShape(true); + diamond->setXthreshold(2); + diamond->setUV(true); + diamond->setUVthreshold(2); + diamond->setRoadSize(8); + diamond->setRoadSizeUpX(4); + diamond->setRoadSizeDownX(0); + diamond->setRoadSizeUpUV(4); + diamond->setRoadSizeDownUV(0); - StatusCode MMTriggerTool::runTrigger() { - - //Retrieve the current run number and event number - const EventInfo* pevt = 0; - ATH_CHECK( evtStore()->retrieve(pevt) ); - int event = pevt->event_ID()->event_number(); - - ////////////////////////////////////////////////////////////// - // // - // Load Variables From Containers into our Data Structures // - // // - ////////////////////////////////////////////////////////////// - - map<hitData_key,hitData_entry> Hits_Data_Set_Time; - map<int,evInf_entry> Event_Info; - - const MmDigitContainer *nsw_MmDigitContainer = nullptr; - ATH_CHECK( evtStore()->retrieve(nsw_MmDigitContainer,"MM_DIGITS") ); - - std::string wedgeType = getWedgeType(nsw_MmDigitContainer); - if(wedgeType=="Large") m_par = m_par_large; - if(wedgeType=="Small") m_par = m_par_small; - if(wedgeType=="Neither") { - ATH_MSG_INFO( "SMALL AND LARGE!! Event did (NOT) pass " ); - return StatusCode::SUCCESS; - } - - MMLoadVariables load = MMLoadVariables(&(*(evtStore())), m_detManager, m_MmIdHelper, m_par); - - std::vector<digitWrapper> entries; - - ATH_CHECK( load.getMMDigitsInfo(entries, Hits_Data_Set_Time, Event_Info) ); - this->fillNtuple(load); - - //Originally boom, this is the saved "particle_info" (originally primer) - evInf_entry truth_info(Event_Info.find(pevt->event_ID()->event_number())->second); - - bool pass_cuts = truth_info.pass_cut; - double trueta = truth_info.eta_ip; - double trupt = truth_info.pt; - - - double tent=truth_info.theta_ent; - double tpos=truth_info.theta_pos; - double ppos=truth_info.phi_pos; - double dt=truth_info.dtheta; - m_trigger_trueThe->push_back(tent); - m_trigger_truePhi->push_back(ppos); - m_trigger_trueDth->push_back(dt); - - //from MMT_Loader >>>> If entry matches find(event) adds element to vector - std::vector<hitData_entry> hitDatas(event_hitDatas(event,Hits_Data_Set_Time)); - //Only consider fits if they satisfy CT and fall in wedge - if(pass_cuts){ - //Make sure hit info is not empy - if(!hitDatas.empty()){ - - ////////////////////////////////////////////////////////////// - // // - // Finder Applied Here // - // // - ////////////////////////////////////////////////////////////// - - //Initialization of the finder: defines all the roads - MMT_Finder find = MMT_Finder(m_par, 1); - - ATH_MSG_DEBUG( "Number of Roads Configured " << find.get_roads() ); - - //Convert hits to slopes and fill the buffer - map<pair<int,int>,finder_entry> hitBuffer; - for(int ihds=0; ihds<(int)hitDatas.size(); ihds++){ - find.fillHitBuffer( hitBuffer, // Map (road,plane) -> Finder entry - hitDatas[ihds].entry_hit(m_par) ); // Hit object - - hitData_info hitInfo = hitDatas[ihds].entry_hit(m_par).info; + return StatusCode::SUCCESS; + } - m_trigger_VMM->push_back(hitDatas[ihds].VMM_chip); - m_trigger_plane->push_back(hitDatas[ihds].plane); - m_trigger_station->push_back(hitDatas[ihds].station_eta); - m_trigger_strip->push_back(hitDatas[ihds].strip); - m_trigger_slope->push_back(hitInfo.slope.getFixed()); + StatusCode MMTriggerTool::runTrigger(const bool do_MMDiamonds) { + //Retrieve the current run number and event number + const EventInfo* pevt = 0; + ATH_CHECK( evtStore()->retrieve(pevt) ); + int event = pevt->event_ID()->event_number(); + ATH_MSG_INFO("********************************************************* EVENT NUMBER = " << event); + + ////////////////////////////////////////////////////////////// + // // + // Load Variables From Containers into our Data Structures // + // // + ////////////////////////////////////////////////////////////// + + const MmDigitContainer *nsw_MmDigitContainer = nullptr; + ATH_CHECK( evtStore()->retrieve(nsw_MmDigitContainer,"MM_DIGITS") ); + + std::map<std::string, MMT_Parameters*> pars; + pars["MML"] = m_par_large; + pars["MMS"] = m_par_small; + MMLoadVariables load = MMLoadVariables(&(*(evtStore())), m_detManager, m_MmIdHelper, pars); + + std::map<std::pair<int, unsigned int>,std::vector<digitWrapper> > entries; + std::map<std::pair<int, unsigned int>,map<hitData_key,hitData_entry> > Hits_Data_Set_Time; + std::map<std::pair<int, unsigned int>,evInf_entry> Event_Info; + ATH_CHECK( load.getMMDigitsInfo(entries, Hits_Data_Set_Time, Event_Info) ); + this->fillNtuple(load); + + unsigned int particles = entries.rbegin()->first.second +1, nskip=0; + + if (entries.empty()) { + ATH_MSG_ERROR("No entries available, something is going wrong"); + return StatusCode::FAILURE; + } - } - if(hitDatas.size()==8){ + for (unsigned int i=0; i<particles; i++) { + double trueta = -999., trupt = -999., dt = -999., tpos = -999., ppos = -999.; + // We need to extract truth info, if available + std::pair<int, unsigned int> pair_event (event,i); + if (!Event_Info.empty()) { + auto tru_it = Event_Info.find(pair_event); + if (tru_it != Event_Info.end()) { + evInf_entry truth_info(tru_it->second); + trueta = truth_info.eta_ip; m_trigger_trueEtaRange->push_back(trueta); - m_trigger_truePtRange->push_back(trupt); - if(wedgeType=="Large") { - m_trigger_large_trueEtaRange->push_back(trueta); - m_trigger_large_truePtRange->push_back(trupt); - } - if(wedgeType=="Small") { - m_trigger_small_trueEtaRange->push_back(trueta); - m_trigger_small_truePtRange->push_back(trupt); - } - - } - - ////////////////////////////////////////////////////////////// - // // - // Fitter Applied Here // - // // - ////////////////////////////////////////////////////////////// - - MMT_Fitter fit = MMT_Fitter(m_par); - - //First loop over the roads and planes and apply the fitter - int fits_occupied=0; - const int nfit_max=1; //MOVE THIS EVENTUALLY - int nRoads = find.get_roads(); - vector<evFit_entry> road_fits = vector<evFit_entry>(nRoads,evFit_entry()); - - //Variables saved for Alex T. for hardware validation - double mxl; - double fillmxl=-999; - double muGlobal; - double mvGlobal; - vector<pair<double,double> > mxmy; - - for(int iRoad=0; iRoad<nRoads; iRoad++){ - - vector<bool> plane_is_hit; - vector<Hit> track; - - //Check if there are hits in the buffer - find.checkBufferForHits( plane_is_hit, // Empty, To be filled by function. - track, // Empty, To be filled by function. - iRoad, // roadID - hitBuffer // All hits. Map ( (road,plane) -> finder_entry ) - ); - - //Look for coincidences - int road_num=find.Coincidence_Gate(plane_is_hit); - - if(road_num>0){ - - if(fits_occupied>=nfit_max) break; - - //Perform the fit -> calculate local, global X, UV slopes -> calculate ROI and TriggerTool signal (theta, phi, deltaTheta) - evFit_entry candidate=fit.fit_event(event,track,hitDatas,fits_occupied,mxmy,mxl,mvGlobal,muGlobal); - - ATH_MSG_DEBUG( "THETA " << candidate.fit_theta.getFixed() << " PHI " << candidate.fit_phi.getFixed() << " DTH " << candidate.fit_dtheta.getFixed() ); - road_fits[iRoad]=candidate; - fillmxl = mxl; - fits_occupied++; - - } - - road_fits[iRoad].hcode=road_num; - - } //end roads + trupt = truth_info.pt; + m_trigger_truePtRange->push_back(trupt); - ////////////////////////////////////////////////////////////// - // // - // Pass the ROI as Signal // - // // - ////////////////////////////////////////////////////////////// + tpos=truth_info.theta_pos; + m_trigger_trueThe->push_back(truth_info.theta_ip); + ppos=truth_info.phi_pos; + m_trigger_truePhi->push_back(truth_info.phi_ip); - if(road_fits.size()==0 and hitDatas.size()==8 ) { - ATH_MSG_DEBUG( "TruthRF0 " << tpos << " " << ppos << " " << dt << " " << trueta ); + dt=truth_info.dtheta; + m_trigger_trueDth->push_back(dt); + } else ATH_MSG_WARNING( "Extra reco particle with no truth candidate available" ); + } + // Now let's switch to reco hits: firstly, extracting the station name we're working on... + std::string station = "-"; + auto event_it = entries.find(pair_event); + station = event_it->second[0].stName; // Station name is taken from the first digit! In MMLoadVariables there's a check to ensure all digits belong to the same station + + // Secondly, extracting the Phi of the station we're working on... + int stationPhi = -999; + digitWrapper dW = event_it->second[0]; + Identifier tmpID = dW.id(); + stationPhi = m_MmIdHelper->stationPhi(tmpID); + + // Finally, let's start with hits + auto reco_it = Hits_Data_Set_Time.find(pair_event); + if (reco_it != Hits_Data_Set_Time.end()) { + if (!reco_it->second.empty()) { + std::vector<hitData_entry> hitDatas; + for (auto hit_it = reco_it->second.begin(); hit_it != reco_it->second.end(); hit_it++) hitDatas.push_back(hit_it->second); + if (do_MMDiamonds) { + /* + * Filling hits for each event: a new class, MMT_Hit, is called in + * order to use both algorithms witghout interferences + */ + diamond->createRoads_fillHits(i-nskip, hitDatas, m_detManager, pars[station], stationPhi); + double smallest_bc = 999999.; + for(int ihds=0; ihds<(int)hitDatas.size(); ihds++) { + if (hitDatas[ihds].BC_time < 0.) continue; + else if (hitDatas[ihds].BC_time < smallest_bc) smallest_bc = hitDatas[ihds].BC_time; + + // The PrintHits function below gives identical results to the following one: hitDatas[ihds].print(); + m_trigger_VMM->push_back(hitDatas[ihds].VMM_chip); + m_trigger_plane->push_back(hitDatas[ihds].plane); + m_trigger_station->push_back(hitDatas[ihds].station_eta); + m_trigger_strip->push_back(hitDatas[ihds].strip); + } + diamond->printHits(i-nskip); + std::vector<double> slopes = diamond->getHitSlopes(); + for (const auto &s : slopes) m_trigger_RZslopes->push_back(s); + diamond->resetSlopes(); + slopes.clear(); + /* + * Here we create roads with all MMT_Hit collected before (if any), then we save the results + */ + if (diamond->getHitVector(i-nskip).size() >= (diamond->getXthreshold()+diamond->getUVthreshold())) { + diamond->findDiamonds(i-nskip, smallest_bc, event); + storeEventProperties(i-nskip); + } + } else { + ////////////////////////////////////////////////////////////// + // // + // Finder Applied Here // + // // + ////////////////////////////////////////////////////////////// + + //Initialization of the finder: defines all the roads + MMT_Finder find = MMT_Finder(pars[station], 1); + ATH_MSG_DEBUG( "Number of Roads Configured " << find.get_roads() ); + + std::map<pair<int,int>,finder_entry> hitBuffer; + for (auto hit_it = reco_it->second.begin(); hit_it != reco_it->second.end(); hit_it++) { + find.fillHitBuffer( hitBuffer, hit_it->second.entry_hit(pars[station]) ); // Hit object, Map (road,plane) -> Finder entry + + hitData_info hitInfo = hit_it->second.entry_hit(pars[station]).info; + m_trigger_VMM->push_back(hit_it->second.VMM_chip); + m_trigger_plane->push_back(hit_it->second.plane); + m_trigger_station->push_back(hit_it->second.station_eta); + m_trigger_strip->push_back(hit_it->second.strip); + m_trigger_slope->push_back(hitInfo.slope); + } + if (reco_it->second.size() > 7) { + m_trigger_trueEtaRange->push_back(trueta); + m_trigger_truePtRange->push_back(trupt); + if (station == "MML") { + m_trigger_large_trueEtaRange->push_back(trueta); + m_trigger_large_truePtRange->push_back(trupt); + } + if (station == "MMS") { + m_trigger_small_trueEtaRange->push_back(trueta); + m_trigger_small_truePtRange->push_back(trupt); + } + } + + ////////////////////////////////////////////////////////////// + // // + // Fitter Applied Here // + // // + ////////////////////////////////////////////////////////////// + + MMT_Fitter fit = MMT_Fitter(pars[station]); + + //First loop over the roads and planes and apply the fitter + int fits_occupied = 0; + const int nfit_max = 1; //MOVE THIS EVENTUALLY + int nRoads = find.get_roads(); + + vector<evFit_entry> road_fits = vector<evFit_entry>(nRoads,evFit_entry()); + + //Variables saved for Alex T. for hardware validation + double mxl; + double fillmxl = -999; + double muGlobal; + double mvGlobal; + vector<pair<double,double> > mxmy; + + for (int iRoad = 0; iRoad < nRoads; iRoad++) { + vector<bool> plane_is_hit; + vector<Hit> track; + + //Check if there are hits in the buffer + find.checkBufferForHits( plane_is_hit, // Empty, To be filled by function. + track, // Empty, To be filled by function. + iRoad, // roadID + hitBuffer // All hits. Map ( (road,plane) -> finder_entry ) + ); + + //Look for coincidences + int road_num = find.Coincidence_Gate(plane_is_hit); + if (road_num > 0) { + if (fits_occupied >= nfit_max) break; + + //Perform the fit -> calculate local, global X, UV slopes -> calculate ROI and TriggerTool signal (theta, phi, deltaTheta) + evFit_entry candidate = fit.fit_event(event,track,hitDatas,fits_occupied,mxmy,mxl,mvGlobal,muGlobal); + + ATH_MSG_DEBUG( "THETA " << candidate.fit_theta << " PHI " << candidate.fit_phi << " DTH " << candidate.fit_dtheta ); + road_fits[iRoad] = candidate; + fillmxl = mxl; + fits_occupied++; + } + road_fits[iRoad].hcode = road_num; + } //end roads + + ////////////////////////////////////////////////////////////// + // // + // Pass the ROI as Signal // + // // + ////////////////////////////////////////////////////////////// + + if (road_fits.size() == 0 and hitDatas.size() == 8) ATH_MSG_DEBUG( "TruthRF0 " << tpos << " " << ppos << " " << dt << " " << trueta ); + + for (unsigned int i = 0; i < road_fits.size(); i++) { + if (road_fits[i].fit_roi == 0 and hitDatas.size() == 8) { + ATH_MSG_DEBUG( "TruthROI0 " << tpos << " " << ppos << " " << dt << " " << trueta ); + } + if (road_fits[i].fit_roi > 0) { + //For the future: how do we want these to pass on as the signal? Some new data structure? + double fitTheta = road_fits[i].fit_theta; + double fitPhi = road_fits[i].fit_phi; + double fitDeltaTheta = road_fits[i].fit_dtheta; + + //need a correction for the fitted phi, taken from phi_shift in MMLoadVariables.cxx (now it's local) + int wedge = 0; + if (station == "MML") wedge = 0; + else if (station == "MMS") wedge = 1; + float n = 2*(stationPhi-1) + wedge; + float shift = n * M_PI/8.; + if(n > 8) shift = (16.-n)*M_PI/8.; + if(n < 8) fitPhi = (fitPhi + shift); + else if(n == 8) fitPhi = (fitPhi + (fitPhi >= 0 ? -1 : 1)*shift); + else if(n > 8) fitPhi = (fitPhi - shift); + + double fitEta = -1. * TMath::Log(TMath::Tan(fitTheta/2.)); //VALE: trueta was filled!!!! + + ATH_MSG_DEBUG( "Truth " << tpos << " " << ppos << " " << dt ); + ATH_MSG_DEBUG( "FIT!! " << fitTheta << " " << fitPhi << " " << fitDeltaTheta ); + m_trigger_fitThe->push_back(fitTheta); + m_trigger_fitPhi->push_back(fitPhi); + m_trigger_fitDth->push_back(fitDeltaTheta); + + m_trigger_resThe->push_back(fitTheta-tpos); + m_trigger_resPhi->push_back(fitPhi-ppos); + m_trigger_resDth->push_back(fitDeltaTheta-dt); + + m_trigger_mx->push_back(mxmy.front().first); + m_trigger_my->push_back(mxmy.front().second); + m_trigger_mxl->push_back(fillmxl); + + m_trigger_mu->push_back(muGlobal); + m_trigger_mv->push_back(mvGlobal); + + m_trigger_fitEtaRange->push_back(fitEta); + m_trigger_fitPtRange->push_back(trupt); + if (station == "MML") { + m_trigger_large_fitEtaRange->push_back(fitEta); + m_trigger_large_fitPtRange->push_back(trupt); + } + if (station == "MMS") { + m_trigger_small_fitEtaRange->push_back(fitEta); + m_trigger_small_fitPtRange->push_back(trupt); + } + } + } + } // if-else Diamond clause + hitDatas.clear(); + } else { + ATH_MSG_WARNING( "Available hits are less than X+UV thresholds, skipping" ); + nskip++; } - for(unsigned int i=0; i<road_fits.size(); i++){ - if(road_fits[i].fit_roi==0 and hitDatas.size()==8) { - ATH_MSG_DEBUG( "TruthROI0 " << tpos << " " << ppos << " " << dt << " " << trueta ); - } - if(road_fits[i].fit_roi>0){ - //For the future: how do we want these to pass on as the signal? Some new data structure? - double fitTheta = road_fits[i].fit_theta.getFixed(); - double fitPhi = road_fits[i].fit_phi.getFixed(); - double fitDeltaTheta = road_fits[i].fit_dtheta.getFixed(); - - - ATH_MSG_DEBUG( "Truth " << tpos << " " << ppos << " " << dt ); - ATH_MSG_DEBUG( "FIT!! " << fitTheta << " " << fitPhi << " " << fitDeltaTheta ); - m_trigger_fitThe->push_back(fitTheta); - m_trigger_fitPhi->push_back(fitPhi); - m_trigger_fitDth->push_back(fitDeltaTheta); - - m_trigger_resThe->push_back(fitTheta-tpos); - m_trigger_resPhi->push_back(fitPhi-ppos); - m_trigger_resDth->push_back(fitDeltaTheta-dt); - - m_trigger_mx->push_back(mxmy.front().first); - m_trigger_my->push_back(mxmy.front().second); - m_trigger_mxl->push_back(fillmxl); - - m_trigger_mu->push_back(muGlobal); - m_trigger_mv->push_back(mvGlobal); - - m_trigger_fitEtaRange->push_back(trueta); - m_trigger_fitPtRange->push_back(trupt); - if(wedgeType=="Large") { - m_trigger_large_fitEtaRange->push_back(trueta); - m_trigger_large_fitPtRange->push_back(trupt); - } - if(wedgeType=="Small") { - m_trigger_small_fitEtaRange->push_back(trueta); - m_trigger_small_fitPtRange->push_back(trupt); - } - - }//fit roi > 0 - } // end road_fits - }//end if hitDataS EMPTY - }//end if PASS_CUTS - - //clear pointers, filled hit info - - Event_Info.erase(Event_Info.find(event)); - vector<hitData_key> kill_keys(event_hitData_keys(event,Hits_Data_Set_Time)); - return StatusCode::SUCCESS; - } - - //Function that find the hits information and hits keys that get stored throughout the run. - //The data structures are defined in MMT_struct - - vector<hitData_key> MMTriggerTool::event_hitData_keys(int find_event, map<hitData_key,hitData_entry>& Hits_Data_Set_Time) const{ - vector<hitData_key> ravel; - int fnd_entries=0; - for(map<hitData_key,hitData_entry>::const_iterator entry=Hits_Data_Set_Time.begin(); entry!=Hits_Data_Set_Time.end(); ++entry){ - if(entry->second.event==find_event){ - ravel.push_back(entry->first); - fnd_entries++; + } else { + ATH_MSG_WARNING( "Empty hit map, skipping" ); + nskip++; } - else if(fnd_entries>0) break;//we use the fact that maps store things according to the strict weak ordering of the key's comparison operator - } - return ravel; + } // Main particle loop + entries.clear(); + Hits_Data_Set_Time.clear(); + Event_Info.clear(); + if(do_MMDiamonds) diamond->clearEvent(); + + return StatusCode::SUCCESS; } - vector<hitData_entry> MMTriggerTool::event_hitDatas(int find_event, map<hitData_key,hitData_entry>& Hits_Data_Set_Time) const{ - vector<hitData_entry> bolero; - int fnd_entries=0; - for(map<hitData_key,hitData_entry>::const_iterator entry=Hits_Data_Set_Time.begin(); entry!=Hits_Data_Set_Time.end(); ++entry){ - if(entry->second.event==find_event){ - bolero.push_back(entry->second); - fnd_entries++; + void MMTriggerTool::storeEventProperties(const unsigned int iterator) { + if (diamond->getDiamond(iterator).wedgeCounter == iterator) { + if (diamond->getSlopeVector(iterator).empty()) return; + m_trigger_diamond_ntrig->push_back(diamond->getSlopeVector(iterator).size()); + + for (const auto &slope : diamond->getSlopeVector(iterator)) { + m_trigger_diamond_stationPhi->push_back(diamond->getDiamond(iterator).phi); + m_trigger_diamond_sector->push_back(diamond->getDiamond(iterator).sector); + m_trigger_diamond_bc->push_back(slope.BC); + m_trigger_diamond_totalCount->push_back(slope.totalCount); + m_trigger_diamond_realCount->push_back(slope.realCount); + m_trigger_diamond_XbkgCount->push_back(slope.xbkg); + m_trigger_diamond_UVbkgCount->push_back(slope.uvbkg); + m_trigger_diamond_XmuonCount->push_back(slope.xmuon); + m_trigger_diamond_UVmuonCount->push_back(slope.uvmuon); + m_trigger_diamond_iX->push_back(slope.iRoad); + m_trigger_diamond_iU->push_back(slope.iRoadu); + m_trigger_diamond_iV->push_back(slope.iRoadv); + m_trigger_diamond_age->push_back(slope.age); + m_trigger_diamond_Xavg->push_back(slope.xavg); + m_trigger_diamond_Uavg->push_back(slope.uavg); + m_trigger_diamond_Vavg->push_back(slope.vavg); + m_trigger_diamond_mxl->push_back(slope.mxl); + m_trigger_diamond_theta->push_back(slope.theta); + m_trigger_diamond_eta->push_back(slope.eta); + m_trigger_diamond_dtheta->push_back(slope.dtheta); + m_trigger_diamond_phi->push_back(slope.phi); + m_trigger_diamond_phiShf->push_back(slope.phiShf); } - else if(fnd_entries>0) break;//we use the fact that maps store things according to the strict weak ordering of the key's comparison operator } - return bolero; + else ATH_MSG_FATAL( "Iterators don't match! Cannot retrieve corresponding variables" ); } - std::string MMTriggerTool::getWedgeType(const MmDigitContainer *nsw_MmDigitContainer){ - - std::vector<bool> isLargeWedge; - //Digit loop to match to truth - for(auto dit : *nsw_MmDigitContainer) { - - const MmDigitCollection* coll = dit; - for (unsigned int item=0; item<coll->size(); item++) { + StatusCode MMTriggerTool::finalizeDiamondAlgorithm(const bool do_MMDiamonds) { + /* + * Place here all the cleaning stuff: pointer deletion etc... + */ + if(do_MMDiamonds) delete diamond; + delete m_par_large; + delete m_par_small; - const MmDigit* digit = coll->at(item); - Identifier Id = digit->identify(); - - std::string stName = m_MmIdHelper->stationNameString(m_MmIdHelper->stationName(Id)); - const string& sname(stName); - if (sname.compare("MML")==0)isLargeWedge.push_back(true); - else isLargeWedge.push_back(false); - } - } - - bool allLarge = true; - bool allSmall = true; - for(unsigned int i=0; i<isLargeWedge.size(); i++){ - if (isLargeWedge.at(i)) allSmall = false; - else allLarge = false; - } - std::string wedgeType = "Neither"; - if (allLarge) wedgeType = "Large"; - if (allSmall) wedgeType = "Small"; - return wedgeType; + return StatusCode::SUCCESS; } }//end namespace - diff --git a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerToolTree.cxx b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerToolTree.cxx index 2088839fcd3364811e9c4e21daf3bda77ec2f90d..97d6b2e21f4b41d00a95bb40dd204cd0a9e1dd04 100644 --- a/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerToolTree.cxx +++ b/Trigger/TrigT1/TrigT1NSWSimTools/src/MMTriggerToolTree.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ // Athena/Gaudi includes @@ -25,7 +25,31 @@ namespace NSWL1 { } StatusCode MMTriggerTool::book_branches() { - + m_trigger_diamond_ntrig = new std::vector<unsigned int>(); + m_trigger_diamond_bc = new std::vector<int>(); + m_trigger_diamond_sector = new std::vector<char>(); + m_trigger_diamond_stationPhi = new std::vector<int>(); + m_trigger_diamond_totalCount = new std::vector<unsigned int>(); + m_trigger_diamond_realCount = new std::vector<unsigned int>(); + m_trigger_diamond_XbkgCount = new std::vector<unsigned int>(); + m_trigger_diamond_UVbkgCount = new std::vector<unsigned int>(); + m_trigger_diamond_XmuonCount = new std::vector<unsigned int>(); + m_trigger_diamond_UVmuonCount = new std::vector<unsigned int>(); + m_trigger_diamond_iX = new std::vector<int>(); + m_trigger_diamond_iU = new std::vector<int>(); + m_trigger_diamond_iV = new std::vector<int>(); + m_trigger_diamond_age = new std::vector<double>(); + m_trigger_diamond_Xavg = new std::vector<double>(); + m_trigger_diamond_Uavg = new std::vector<double>(); + m_trigger_diamond_Vavg = new std::vector<double>(); + m_trigger_diamond_mxl = new std::vector<double>(); + m_trigger_diamond_theta = new std::vector<double>(); + m_trigger_diamond_eta = new std::vector<double>(); + m_trigger_diamond_dtheta = new std::vector<double>(); + m_trigger_diamond_phi = new std::vector<double>(); + m_trigger_diamond_phiShf = new std::vector<double>(); + + m_trigger_RZslopes = new std::vector<double>(); m_trigger_fitThe = new std::vector<double>(); m_trigger_fitPhi = new std::vector<double>(); m_trigger_fitDth = new std::vector<double>(); @@ -34,7 +58,7 @@ namespace NSWL1 { m_trigger_fitEtaRange = new std::vector<double>(); m_trigger_fitPtRange = new std::vector<double>(); m_trigger_resThe = new std::vector<double>(); - m_trigger_resPhi = new std::vector<double>(); + m_trigger_resPhi = new std::vector<double>(); m_trigger_resDth = new std::vector<double>(); m_trigger_large_fitThe = new std::vector<double>(); @@ -42,10 +66,13 @@ namespace NSWL1 { m_trigger_large_fitDth = new std::vector<double>(); m_trigger_large_trueEtaRange = new std::vector<double>(); m_trigger_large_truePtRange = new std::vector<double>(); + m_trigger_large_trueThe = new std::vector<double>(); + m_trigger_large_truePhi = new std::vector<double>(); + m_trigger_large_trueDth = new std::vector<double>(); m_trigger_large_fitEtaRange = new std::vector<double>(); m_trigger_large_fitPtRange = new std::vector<double>(); m_trigger_large_resThe = new std::vector<double>(); - m_trigger_large_resPhi = new std::vector<double>(); + m_trigger_large_resPhi = new std::vector<double>(); m_trigger_large_resDth = new std::vector<double>(); m_trigger_small_fitThe = new std::vector<double>(); @@ -53,10 +80,13 @@ namespace NSWL1 { m_trigger_small_fitDth = new std::vector<double>(); m_trigger_small_trueEtaRange = new std::vector<double>(); m_trigger_small_truePtRange = new std::vector<double>(); + m_trigger_small_trueThe = new std::vector<double>(); + m_trigger_small_truePhi = new std::vector<double>(); + m_trigger_small_trueDth = new std::vector<double>(); m_trigger_small_fitEtaRange = new std::vector<double>(); m_trigger_small_fitPtRange = new std::vector<double>(); m_trigger_small_resThe = new std::vector<double>(); - m_trigger_small_resPhi = new std::vector<double>(); + m_trigger_small_resPhi = new std::vector<double>(); m_trigger_small_resDth = new std::vector<double>(); m_trigger_VMM = new std::vector<int>(); @@ -161,6 +191,31 @@ namespace NSWL1 { std::string ToolName = name().substr( name().find("::")+2,std::string::npos ); const char* n = ToolName.c_str(); + m_tree->Branch(TString::Format("%s_trigger_diamond_bc",n).Data(), &m_trigger_diamond_bc); + m_tree->Branch(TString::Format("%s_trigger_diamond_ntrig",n).Data(), &m_trigger_diamond_ntrig); + m_tree->Branch(TString::Format("%s_trigger_diamond_sector",n).Data(), &m_trigger_diamond_sector); + m_tree->Branch(TString::Format("%s_trigger_diamond_stationPhi",n).Data(), &m_trigger_diamond_stationPhi); + m_tree->Branch(TString::Format("%s_trigger_diamond_totalCount",n).Data(), &m_trigger_diamond_totalCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_realCount",n).Data(), &m_trigger_diamond_realCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_XbkgCount",n).Data(), &m_trigger_diamond_XbkgCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_UVbkgCount",n).Data(), &m_trigger_diamond_UVbkgCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_XmuonCount",n).Data(), &m_trigger_diamond_XmuonCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_UVmuonCount",n).Data(), &m_trigger_diamond_UVmuonCount); + m_tree->Branch(TString::Format("%s_trigger_diamond_iX",n).Data(), &m_trigger_diamond_iX); + m_tree->Branch(TString::Format("%s_trigger_diamond_iU",n).Data(), &m_trigger_diamond_iU); + m_tree->Branch(TString::Format("%s_trigger_diamond_iV",n).Data(), &m_trigger_diamond_iV); + m_tree->Branch(TString::Format("%s_trigger_diamond_age",n).Data(), &m_trigger_diamond_age); + m_tree->Branch(TString::Format("%s_trigger_diamond_Xavg",n).Data(), &m_trigger_diamond_Xavg); + m_tree->Branch(TString::Format("%s_trigger_diamond_Uavg",n).Data(), &m_trigger_diamond_Uavg); + m_tree->Branch(TString::Format("%s_trigger_diamond_Vavg",n).Data(), &m_trigger_diamond_Vavg); + m_tree->Branch(TString::Format("%s_trigger_diamond_mxl",n).Data(), &m_trigger_diamond_mxl); + m_tree->Branch(TString::Format("%s_trigger_diamond_theta",n).Data(), &m_trigger_diamond_theta); + m_tree->Branch(TString::Format("%s_trigger_diamond_eta",n).Data(), &m_trigger_diamond_eta); + m_tree->Branch(TString::Format("%s_trigger_diamond_dtheta",n).Data(), &m_trigger_diamond_dtheta); + m_tree->Branch(TString::Format("%s_trigger_diamond_phi",n).Data(), &m_trigger_diamond_phi); + m_tree->Branch(TString::Format("%s_trigger_diamond_phiShf",n).Data(), &m_trigger_diamond_phiShf); + + m_tree->Branch(TString::Format("%s_trigger_RZslopes",n).Data(),&m_trigger_RZslopes); m_tree->Branch(TString::Format("%s_trigger_fitThe",n).Data(),&m_trigger_fitThe); m_tree->Branch(TString::Format("%s_trigger_fitPhi",n).Data(), &m_trigger_fitPhi); m_tree->Branch(TString::Format("%s_trigger_fitDth",n).Data(), &m_trigger_fitDth); @@ -171,12 +226,14 @@ namespace NSWL1 { m_tree->Branch(TString::Format("%s_trigger_resThe",n).Data(), &m_trigger_resThe); m_tree->Branch(TString::Format("%s_trigger_resPhi",n).Data(), &m_trigger_resPhi); m_tree->Branch(TString::Format("%s_trigger_resDth",n).Data(), &m_trigger_resDth); - m_tree->Branch(TString::Format("%s_trigger_large_fitThe",n).Data(),&m_trigger_large_fitThe); m_tree->Branch(TString::Format("%s_trigger_large_fitPhi",n).Data(), &m_trigger_large_fitPhi); m_tree->Branch(TString::Format("%s_trigger_large_fitDth",n).Data(), &m_trigger_large_fitDth); m_tree->Branch(TString::Format("%s_trigger_large_trueEtaRange",n).Data(), &m_trigger_large_trueEtaRange); m_tree->Branch(TString::Format("%s_trigger_large_truePtRange",n).Data(), &m_trigger_large_truePtRange); + m_tree->Branch(TString::Format("%s_trigger_large_trueThe",n).Data(), &m_trigger_large_trueThe); + m_tree->Branch(TString::Format("%s_trigger_large_truePhi",n).Data(), &m_trigger_large_truePhi); + m_tree->Branch(TString::Format("%s_trigger_large_trueDth",n).Data(), &m_trigger_large_trueDth); m_tree->Branch(TString::Format("%s_trigger_large_fitEtaRange",n).Data(), &m_trigger_large_fitEtaRange); m_tree->Branch(TString::Format("%s_trigger_large_fitPtRange",n).Data(), &m_trigger_large_fitPtRange); m_tree->Branch(TString::Format("%s_trigger_large_resThe",n).Data(), &m_trigger_large_resThe); @@ -188,6 +245,9 @@ namespace NSWL1 { m_tree->Branch(TString::Format("%s_trigger_small_fitDth",n).Data(), &m_trigger_small_fitDth); m_tree->Branch(TString::Format("%s_trigger_small_trueEtaRange",n).Data(), &m_trigger_small_trueEtaRange); m_tree->Branch(TString::Format("%s_trigger_small_truePtRange",n).Data(), &m_trigger_small_truePtRange); + m_tree->Branch(TString::Format("%s_trigger_small_trueThe",n).Data(), &m_trigger_small_trueThe); + m_tree->Branch(TString::Format("%s_trigger_small_truePhi",n).Data(), &m_trigger_small_truePhi); + m_tree->Branch(TString::Format("%s_trigger_small_trueDth",n).Data(), &m_trigger_small_trueDth); m_tree->Branch(TString::Format("%s_trigger_small_fitEtaRange",n).Data(), &m_trigger_small_fitEtaRange); m_tree->Branch(TString::Format("%s_trigger_small_fitPtRange",n).Data(), &m_trigger_small_fitPtRange); m_tree->Branch(TString::Format("%s_trigger_small_resThe",n).Data(), &m_trigger_small_resThe); @@ -301,6 +361,31 @@ namespace NSWL1 { //clear the ntuple variables if(m_tree==0) return; + m_trigger_diamond_bc->clear(); + m_trigger_diamond_ntrig->clear(); + m_trigger_diamond_sector->clear(); + m_trigger_diamond_stationPhi->clear(); + m_trigger_diamond_totalCount->clear(); + m_trigger_diamond_realCount->clear(); + m_trigger_diamond_XbkgCount->clear(); + m_trigger_diamond_UVbkgCount->clear(); + m_trigger_diamond_XmuonCount->clear(); + m_trigger_diamond_UVmuonCount->clear(); + m_trigger_diamond_iX->clear(); + m_trigger_diamond_iU->clear(); + m_trigger_diamond_iV->clear(); + m_trigger_diamond_age->clear(); + m_trigger_diamond_Xavg->clear(); + m_trigger_diamond_Uavg->clear(); + m_trigger_diamond_Vavg->clear(); + m_trigger_diamond_mxl->clear(); + m_trigger_diamond_theta->clear(); + m_trigger_diamond_eta->clear(); + m_trigger_diamond_dtheta->clear(); + m_trigger_diamond_phi->clear(); + m_trigger_diamond_phiShf->clear(); + + m_trigger_RZslopes->clear(); m_trigger_fitThe->clear(); m_trigger_fitPhi->clear(); m_trigger_fitDth->clear(); @@ -317,6 +402,9 @@ namespace NSWL1 { m_trigger_large_fitDth->clear(); m_trigger_large_trueEtaRange->clear(); m_trigger_large_truePtRange->clear(); + m_trigger_large_trueThe->clear(); + m_trigger_large_truePhi->clear(); + m_trigger_large_trueDth->clear(); m_trigger_large_fitEtaRange->clear(); m_trigger_large_fitPtRange->clear(); m_trigger_large_resThe->clear(); @@ -328,6 +416,9 @@ namespace NSWL1 { m_trigger_small_fitDth->clear(); m_trigger_small_trueEtaRange->clear(); m_trigger_small_truePtRange->clear(); + m_trigger_small_trueThe->clear(); + m_trigger_small_truePhi->clear(); + m_trigger_small_trueDth->clear(); m_trigger_small_fitEtaRange->clear(); m_trigger_small_fitPtRange->clear(); m_trigger_small_resThe->clear();