/******************************************************************************* * Description: * * ============ * * HighPtElectrons : an electroweak analysis module dealing with * * high pT electrons. * * * * The "fillHistograms" entry point is where you should add code to * * process event data; define histograms & ntuples in "beginJob" * * * * * * Date | Author | Comments : * * ---------+--------------+---------------------------------------------- * * 7.7.99 Dave Waters Original implementation. Simply fills a * * few histograms with electromagnetic cluster * * and matching track information, and fills * * a histogram with E/p. * * * * 26.7.01 Giulia Manca Left original version untouched and added many * * new variables in the electron Ntuple. * * The tree has different branches which can be * * filled or not according to the correspondent * * switch.For more explanation please contact me at * * manca@fnal.gov or (630) 840 2147 * ******************************************************************************/ //----------------------- // This Class's Header -- //----------------------- #include "HighPtElectrons.hh" //------------- // C Headers -- //------------- #include #include //------------------------------- // Collaborating Class Headers -- //------------------------------- #include "AbsEnv/AbsEnv.hh" #include "Edm/Id.hh" #include "CalorGeometry/TowerGeometry.hh" #include "CalorGeometry/CalConstants.hh" #include "CalorObjects/Pes2dClusterColl.hh" #include "Calor/CprWireCollectionMaker.hh" #include "Electron/emobj_alg.hh" #include "ElectronObjects/CdfEmObject.hh" #include "ElectronObjects/CdfEmObjectColl.hh" #include "ElectronObjects/EmCluster.hh" #include "HepTuple/HepHist1D.h" #include "HepTuple/HepHist2D.h" #include "HepTuple/HepNtuple.h" #include "MetObjects/CdfMet.hh" #include "RRL3/GrowResultList.hh" #include "SimulationObjects/HEPG_StorableBank.hh" // DSW 23.08.2002. TrackingHL/Regions REMOVED ! // #include "TrackingHL/Regions/RegionalSelectionCriteria.hh" #include "TrackingObjects/CT_Track/CT_TrackSet.hh" #include "TrackingObjects/SiData/SiClusterSet.hh" #include "TrackingObjects/Storable/CdfTrackColl.hh" #include "TrackingObjects/Storable/CdfTrackView.hh" // #include "TrackingObjects/Storable/StorableSiClusterData.hh" #include "TrackingObjects/Tracks/CdfTrack.hh" #include "TrackingObjects/Tracks/track_cut.hh" #include "TrackingUtils/GeneratorLevel/GenPrimaryVertex.hh" #include "TrackingUtils/GeneratorLevel/GenPrimaryVertexSet.hh" #include "TrackingUtils/GeneratorLevel/ParentedParticle.hh" #include "TrackingUtils/SelectionCriteria/SelectionCriteria.hh" #include "TrackingUtils/Trybos/TRYGenPrimaryVertexSet.hh" #include "Trybos/TRY_Record_Iter_Same.hh" //#include "Edm/ConstHandle.hh" //#include "Edm/Handle.hh" #include "CLHEP/Geometry/Point3D.h" #include "Calor/CprWireCollectionMaker.hh" #include #include "ShowerMax/StripCollectionMaker.hh" #include "SimulationObjects/OBSP_StorableBank.hh" #include "SimulationObjects/OBSV_StorableBank.hh" #include "StripChamberGeometry/SimpleCprPathFinder.hh" #include "StripChamberGeometry/CprIntersection.hh" #include "Trajectory/Helix.hh" #include "TriggerObjects/TFRD_StorableBank.hh" #include "TriggerObjects/TL2D_StorableBank.hh" #include "TriggerObjects/TL3D_StorableBank.hh" static const char rcsid[] = " "; //---------------- // Constructors -- //---------------- HighPtElectrons::HighPtElectrons(const char* const theName, const char* const theDescription ) : HepHistModule( theName, theDescription ), _debug("debug",this,false), _fill_OBS("fill_OBS",this,false), _fill_HEPG("fill_HEPG",this,false), _fill_ELEC("fill_ELEC",this,true), _fill_TRIG("fill_TRIG",this,false), _fill_MET("fill_MET",this,true), _useTrkZVert("useTrkZVert",this,false), _ClusterHadEmCut("ClusterHadEmCut",this,9999.0), _EM_EnergyScaleFactor("EM_EnergyScaleFactor",this,1.0), _EM_ClusterPtCut("EM_ClusterPtCut",this,0.0), _EM_ClusterAbsEtaRange("EM_ClusterAbsEtaRange",this,7.0), _MatchingTrackPtCut("MatchingTrackPtCut",this,0.0), _regionalHistos("RegionalHistos",this,false), _makeHistos("MakeHistos",this,false) { // ---------------------------------------------- // Debug flag : // ---------------------------------------------- std::ostringstream debugDesc; debugDesc << "\t Print out debug information ? (default = false)"; _debug.addDescription(debugDesc.str()); commands()->append(&_debug); std::ostringstream histoDesc; histoDesc << "\t Make some histograms (default = false)"; _makeHistos.addDescription(histoDesc.str()); commands()->append(&_makeHistos); std::ostringstream fill_OBSDesc; fill_OBSDesc << "\t Fill the OBS Banks information (default = false) " ; _fill_OBS.addDescription(fill_OBSDesc.str()); commands()->append(&_fill_OBS); std::ostringstream fill_HEPGDesc; fill_HEPGDesc << "\t Fill the HEPG Banks information (default = false)" ; _fill_HEPG.addDescription(fill_HEPGDesc.str()); commands()->append(&_fill_HEPG); std::ostringstream fill_ELECDesc; fill_ELECDesc << "\t Fill the Electron Block information (default = true)" ; _fill_ELEC.addDescription(fill_ELECDesc.str()); commands()->append(&_fill_ELEC); std::ostringstream fill_TRIGDesc; fill_TRIGDesc << "\t Fill the Trigger Block information (default = false)"; _fill_TRIG.addDescription(fill_TRIGDesc.str()); commands()->append(&_fill_TRIG); std::ostringstream fill_METDesc; fill_METDesc << "\t Fill the MET Block information (default = true) " ; _fill_MET.addDescription(fill_METDesc.str()); commands()->append(&_fill_MET); std::ostringstream useTrkZVertDesc; useTrkZVertDesc << "\t Set the Zvertex from the track to be used for the calculation of all the Em quantities (default = false) " ; _useTrkZVert.addDescription(useTrkZVertDesc.str()); commands()->append(&_useTrkZVert); // ---------------------------------------------- // Cuts Menu : // ---------------------------------------------- _CutsMenu.initialize("CutsMenu",this); _CutsMenu.initTitle("CutsMenu"); commands()->append(&_CutsMenu); std::ostringstream ClusterHadEmCutDesc; ClusterHadEmCutDesc << "\t Maximum Had/Em for EM cluster candidates (default = " << _ClusterHadEmCut.value() << ")."; _ClusterHadEmCut.addDescription(ClusterHadEmCutDesc.str()); _CutsMenu.commands()->append(&_ClusterHadEmCut); std::ostringstream EM_EnergyScaleFactorDesc; EM_EnergyScaleFactorDesc << "\t Scale factor for EM cluster energies (default = " << _EM_EnergyScaleFactor.value() << ")."; _EM_EnergyScaleFactor.addDescription(EM_EnergyScaleFactorDesc.str()); _CutsMenu.commands()->append(&_EM_EnergyScaleFactor); std::ostringstream EM_ClusterPtCutDesc; EM_ClusterPtCutDesc << "\t Minimum Pt of Electromagnetic Clusters (default = " << _EM_ClusterPtCut.value() << " GeV)."; _EM_ClusterPtCut.addDescription(EM_ClusterPtCutDesc.str()); _CutsMenu.commands()->append(&_EM_ClusterPtCut); std::ostringstream EM_ClusterAbsEtaRangeDesc; EM_ClusterAbsEtaRangeDesc << "\t Maximum |eta| of Electromagnetic Clusters (default = " << _EM_ClusterAbsEtaRange.value() << ")."; _EM_ClusterAbsEtaRange.addDescription(EM_ClusterAbsEtaRangeDesc.str()); _CutsMenu.commands()->append(&_EM_ClusterAbsEtaRange); std::ostringstream MatchingTrackPtCutDesc; MatchingTrackPtCutDesc << "\t Minimum Pt of Matching Tracks (default = " << _MatchingTrackPtCut.value() << " GeV)."; _MatchingTrackPtCut.addDescription(MatchingTrackPtCutDesc.str()); _CutsMenu.commands()->append(&_MatchingTrackPtCut); commands()->append(&_regionalHistos); } //-------------- // Destructor -- //-------------- HighPtElectrons::~HighPtElectrons( ) { } //-------------- // Operations -- //-------------- AppResult HighPtElectrons::bookHistograms( ) { HepFileManager* manager = fileManager( ); assert(0 != manager); if(_makeHistos.value()) { _ClusterEtot = &manager->hist1D("Cluster Etot (GeV)", 200 , 0.0, 100.0, 100); assert(0 != _ClusterEtot); _ClusterE_Em = &manager->hist1D("Cluster E_Em (GeV)", 200 , 0.0, 100.0, 101); assert(0 != _ClusterE_Em); _ClusterHadEm = &manager->hist1D("Cluster HadEm (GeV)", 200 , 0.0, 1.0, 102); assert(0 != _ClusterHadEm); _ClusterNumTow = &manager->hist1D("Cluster num. tow.", 100 , 0.0, 100.0, 190); assert(0 != _ClusterNumTow); _ClusterEtaPhiArea = &manager->hist1D("Cluster eta-phi area", 100 , 0.0, 1.0, 191); assert(0 != _ClusterEtaPhiArea); _EM_ClusterPt = &manager->hist1D("EM cluster pT (GeV)", 100 , 0.0, 100.0, 200); assert(0 != _EM_ClusterPt); _EM_ClusterEta = &manager->hist1D("EM cluster eta", 100 , -3.0, 3.0, 201); assert(0 != _EM_ClusterEta); _EM_ClusterPhi = &manager->hist1D("EM cluster phi", 100 , 0.0, 2*M_PI, 202); assert(0 != _EM_ClusterPhi); _MatchingTrackPt = &manager->hist1D("Matching Track pT (GeV)", 100, 0.0, 100.0, 300); assert(0 != _MatchingTrackPt); _MatchingTrackEta = &manager->hist1D("Matching Track eta", 100, -3.0, 3.0, 301); assert(0 != _MatchingTrackEta); _MatchingTrackPhi = &manager->hist1D("Matching Track phi", 100, 0.0, 2*M_PI, 302); assert(0 != _MatchingTrackPhi); _EoverP = &manager->hist1D("E over p for cluster, matching track", 200, 0.0, 4.0, 400); assert(0 != _EoverP); _TransverseMass = &manager->hist1D("Transverse Mass (GeV)", 150, 0.0, 150.0, 500); assert(0 != _TransverseMass); _InvariantMass = &manager->hist1D("Invariant Mass (2 EM Clusters, GeV)", 25, 0.0, 125.0, 600 ); assert(0 != _InvariantMass); // _nClustersVsEta = &manager->hist2D("No. Clusters vs. Eta)",40,1.0,3.0,25,0.,25., 1000); _nCesCls= &manager->hist1D("No.Shwr Max Clusters",100, 0 ,25 , 100); assert(0 != _nCesCls); _nCesDefWireCls= &manager->hist1D("No. Ces Default Wire Clusters",100, 0 ,25 , 101); assert(0 != _nCesDefWireCls); _nCesDefStripCls= &manager->hist1D("No. Ces Default Strip Clusters",100, 0 ,25 , 102); assert(0 != _nCesDefStripCls); _nMatchTrks= &manager->hist1D("No. CdfEmObj Matching Tracks",100, 0 ,25 , 104); assert(0 != _nMatchTrks); _nCprDefCls= &manager->hist1D("No. Cpr Default Clusters",100, 0 ,25 , 105); assert(0 != _nCprDefCls); } // -------------------------------------------------------------------------------- // Build up the blocks in the column-wise ntuple : _ElectronNtuple = &manager->ntuple("ElectronNtuple",2000); _ElectronNtuple->columnAt("GLB::nEvt",&_nEvt,0); _ElectronNtuple->columnAt("GLB::nRun",&_nRun,0); if(_fill_MET.value()) { _ElectronNtuple->columnAt("MET::MetEt",&_MetEt,0.); _ElectronNtuple->columnAt("MET::MetPhi",&_MetPhi,0.); } if(_fill_TRIG.value()) { //_MAXTRG = 128; _ElectronNtuple->columnAt("TRIG::nTrig",&_nTrig,0); _ElectronNtuple->span("TRIG::nTrig",0,_MAXTRG); _ElectronNtuple->columnArray("TRIG::L1Bit",0.,0); _ElectronNtuple->dimension("TRIG::L1Bit",_MAXTRG); _ElectronNtuple->index("TRIG::L1Bit","nTrig"); _ElectronNtuple->columnArray("TRIG::L2Bit",0.,0); _ElectronNtuple->dimension("TRIG::L2Bit",_MAXTRG); _ElectronNtuple->index("TRIG::L2Bit","nTrig"); _ElectronNtuple->columnArray("TRIG::L3Bit",0.,0); _ElectronNtuple->dimension("TRIG::L3Bit",_MAXTRG); _ElectronNtuple->index("TRIG::L3Bit","nTrig"); } if(_fill_ELEC.value()) { //_MAXELE = 100; _ElectronNtuple->columnAt("ELEC::nElec",&_nElec,0); _ElectronNtuple->span("ELEC::nElec",0,_MAXELE); _ElectronNtuple->columnArray("ELEC::ElecReg",0.,0); _ElectronNtuple->dimension("ELEC::ElecReg",_MAXELE); _ElectronNtuple->index("ELEC::ElecReg","nElec"); _ElectronNtuple->columnArray("ELEC::ElecEtot",0.,0.); _ElectronNtuple->dimension("ELEC::ElecEtot",_MAXELE); _ElectronNtuple->index("ELEC::ElecEtot","nElec"); _ElectronNtuple->columnArray("ELEC::ElecEmEt",0.,0.); _ElectronNtuple->dimension("ELEC::ElecEmEt",_MAXELE); _ElectronNtuple->index("ELEC::ElecEmEt","nElec"); _ElectronNtuple->columnArray("ELEC::ElecHadEt",0.,0.); _ElectronNtuple->dimension("ELEC::ElecHadEt",_MAXELE); _ElectronNtuple->index("ELEC::ElecHadEt","nElec"); _ElectronNtuple->columnArray("ELEC::ElecZCentroid",0.,0.); _ElectronNtuple->dimension("ELEC::ElecZCentroid",_MAXELE); _ElectronNtuple->index("ELEC::ElecZCentroid","nElec"); _ElectronNtuple->columnArray("ELEC::ElecHadEm",0.,0.); _ElectronNtuple->dimension("ELEC::ElecHadEm",_MAXELE); _ElectronNtuple->index("ELEC::ElecHadEm","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrigHadEm",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrigHadEm",_MAXELE); _ElectronNtuple->index("ELEC::ElecTrigHadEm","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoTot4",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoTot4",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoTot4","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoTot7",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoTot7",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoTot7","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoEm4",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoEm4",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoEm4","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoEm7",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoEm7",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoEm7","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoHad4",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoHad4",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoHad4","nElec"); _ElectronNtuple->columnArray("ELEC::ElecIsoHad7",0.,0.); _ElectronNtuple->dimension("ELEC::ElecIsoHad7",_MAXELE); _ElectronNtuple->index("ELEC::ElecIsoHad7","nElec"); _ElectronNtuple->columnArray("ELEC::ElecLShwr",0.,0.); _ElectronNtuple->dimension("ELEC::ElecLShwr",_MAXELE); _ElectronNtuple->index("ELEC::ElecLShwr","nElec"); _ElectronNtuple->columnArray("ELEC::ElecLShwr2",0.,0.); _ElectronNtuple->dimension("ELEC::ElecLShwr2",_MAXELE); _ElectronNtuple->index("ELEC::ElecLShwr2","nElec"); _ElectronNtuple->columnArray("ELEC::ElecPx",0.,0.); _ElectronNtuple->dimension("ELEC::ElecPx",_MAXELE); _ElectronNtuple->index("ELEC::ElecPx","nElec"); _ElectronNtuple->columnArray("ELEC::ElecPy",0.,0.); _ElectronNtuple->dimension("ELEC::ElecPy",_MAXELE); _ElectronNtuple->index("ELEC::ElecPy","nElec"); _ElectronNtuple->columnArray("ELEC::ElecPz",0.,0.); _ElectronNtuple->dimension("ELEC::ElecPz",_MAXELE); _ElectronNtuple->index("ELEC::ElecPz","nElec"); _ElectronNtuple->columnArray("ELEC::ElecEta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecEta",_MAXELE); _ElectronNtuple->index("ELEC::ElecEta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecPhi",0.,0.); _ElectronNtuple->dimension("ELEC::ElecPhi",_MAXELE); _ElectronNtuple->index("ELEC::ElecPhi","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTheta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTheta",_MAXELE); _ElectronNtuple->index("ELEC::ElecTheta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkPt",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkPt",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkPt","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkPhi",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkPhi",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkPhi","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkEta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkEta",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkEta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkEta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkEta",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkEta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkD0",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkD0",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkD0","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkZ0",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkZ0",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkZ0","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkCot",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkCot",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkCot","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkCurv",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkCurv",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkCurv","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkAxHits",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkAxHits",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkAxHits","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkStHits",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkStHits",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkStHits","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkAlg",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkAlg",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkAlg","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkId",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkId",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkId","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkXAtCes",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkXAtCes",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkXAtCes","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkZAtCes",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkZAtCes",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkZAtCes","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkEtaAtCes",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkEtaAtCes",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkEtaAtCes","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkPhiAtCes",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkPhiAtCes",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkPhiAtCes","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkSideOfCes",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkSideOfCes",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkSideOfCes","nElec"); _ElectronNtuple->columnArray("ELEC::ElecTrkCesWedge",0.,0.); _ElectronNtuple->dimension("ELEC::ElecTrkCesWedge",_MAXELE,_MAXTRK); _ElectronNtuple->index("ELEC::ElecTrkCesWedge","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkId",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkId",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkId","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkIso",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkIso",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkIso","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkPt",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkPt",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkPt","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkPhi",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkPhi",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkPhi","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkEta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkEta",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkEta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkTheta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkTheta",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkTheta","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkD0",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkD0",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkD0","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkZ0",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkZ0",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkZ0","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkX",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkX",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkX","nElec"); _ElectronNtuple->columnArray("ELEC::ElecMaxTrkZ",0.,0.); _ElectronNtuple->dimension("ELEC::ElecMaxTrkZ",_MAXELE); _ElectronNtuple->index("ELEC::ElecMaxTrkZ","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesX",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesX",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesX","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesXChiW",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesXChiW",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesXChiW","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesDeltaX",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesDeltaX",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesDeltaX","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesZ",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesZ",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesZ","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesZChiS",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesZChiS",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesZChiS","nElec"); _ElectronNtuple->columnArray("ELEC::ElecBestCesDeltaZ",0.,0.); _ElectronNtuple->dimension("ELEC::ElecBestCesDeltaZ",_MAXELE); _ElectronNtuple->index("ELEC::ElecBestCesDeltaZ","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesView",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesView",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesView","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesLocCoord",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesLocCoord",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesLocCoord","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesGlbCoord",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesGlbCoord",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesGlbCoord","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesEnergy",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesEnergy",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesEnergy","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesChi2",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesChi2",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesChi2","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesDistTrkCls",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesDistTrkCls",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesDistTrkCls","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCesAssTrkId",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCesAssTrkId",_MAXELE,_MAXCES); _ElectronNtuple->index("ELEC::ElecCesAssTrkId","nElec"); _ElectronNtuple->columnArray("ELEC::ElecNCesDefCls",0.,0); _ElectronNtuple->dimension("ELEC::ElecNCesDefCls",_MAXELE); _ElectronNtuple->index("ELEC::ElecNCesDefCls","nElec"); _ElectronNtuple->columnArray("ELEC::ElecNCesWirCls",0.,0); _ElectronNtuple->dimension("ELEC::ElecNCesWirCls",_MAXELE); _ElectronNtuple->index("ELEC::ElecNCesWirCls","nElec"); _ElectronNtuple->columnArray("ELEC::ElecNCesStrCls",0.,0); _ElectronNtuple->dimension("ELEC::ElecNCesStrCls",_MAXELE); _ElectronNtuple->index("ELEC::ElecNCesStrCls","nElec"); _ElectronNtuple->columnArray("ELEC::ElecNMatchTrks",0.,0); _ElectronNtuple->dimension("ELEC::ElecNMatchTrks",_MAXELE); _ElectronNtuple->index("ELEC::ElecNMatchTrks","nElec"); _ElectronNtuple->columnArray("ELEC::ElecNCprCls",0.,0.); _ElectronNtuple->dimension("ELEC::ElecNCprCls",_MAXELE); _ElectronNtuple->index("ELEC::ElecNCprCls","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCprCoord",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCprCoord",_MAXELE); _ElectronNtuple->index("ELEC::ElecCprCoord","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCprCharge",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCprCharge",_MAXELE); _ElectronNtuple->index("ELEC::ElecCprCharge","nElec"); _ElectronNtuple->columnArray("ELEC::ElecCprDelta",0.,0.); _ElectronNtuple->dimension("ELEC::ElecCprDelta",_MAXELE); _ElectronNtuple->index("ELEC::ElecCprDelta","nElec"); } if(_fill_OBS.value()) { //_MAXOBS = 250; _ElectronNtuple->columnAt("OBS::nElecOBS",&_nElecOBS,0); _ElectronNtuple->span("OBS::nElecOBS",0,_MAXOBS); _ElectronNtuple->columnArray("OBS::ElecOBSP_Px",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Px",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Px","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Py",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Py",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Py","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Pz",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Pz",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Pz","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Ptot",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Ptot",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Ptot","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Pt",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Pt",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Pt","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Theta",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Theta",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Theta","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Eta",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Eta",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Eta","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Phi",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Phi",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Phi","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_Curv",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_Curv",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_Curv","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprX",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprX",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprX","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_globCprX",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_globCprX",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_globCprX","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprY",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprY",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprY","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprZ",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprZ",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprZ","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprPhi",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprPhi",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprPhi","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprTheta",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprTheta",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprTheta","nElecOBS"); _ElectronNtuple->columnArray("OBS::ElecOBSP_cprEta",0.,0.); _ElectronNtuple->dimension("OBS::ElecOBSP_cprEta",_MAXVER,_MAXOBS); _ElectronNtuple->index("OBS::ElecOBSP_cprEta","nElecOBS"); //OBSV informations //_MAXVER = 100; _ElectronNtuple->columnAt("OBSV::nVert",&_nVert,0); _ElectronNtuple->span("OBSV::nVert",0,_MAXVER); _ElectronNtuple->columnArray("OBSV::ElecOBSV_Vx",0.,0.); _ElectronNtuple->dimension("OBSV::ElecOBSV_Vx",_MAXVER); _ElectronNtuple->index("OBSV::ElecOBSV_Vx","nVert"); _ElectronNtuple->columnArray("OBSV::ElecOBSV_Vy",0.,0.); _ElectronNtuple->dimension("OBSV::ElecOBSV_Vy",_MAXVER); _ElectronNtuple->index("OBSV::ElecOBSV_Vy","nVert"); _ElectronNtuple->columnArray("OBSV::ElecOBSV_Vz",0.,0.); _ElectronNtuple->dimension("OBSV::ElecOBSV_Vz",_MAXVER); _ElectronNtuple->index("OBSV::ElecOBSV_Vz","nVert"); _ElectronNtuple->columnArray("OBSV::ElecOBSV_nDau",0.,0.); _ElectronNtuple->dimension("OBSV::ElecOBSV_nDau",_MAXVER); _ElectronNtuple->index("OBSV::ElecOBSV_nDau","nVert"); } if(_fill_HEPG.value()) { //_MAXHEP = 1000; _ElectronNtuple->columnAt("HEPG::nHepgEvt",&_nHepgEvt,0); _ElectronNtuple->columnAt("HEPG::nHepg",&_nHepg,0); _ElectronNtuple->columnAt("HEPG::nHepgPrim",&_nHepgPrim,0); _ElectronNtuple->columnAt("HEPG::nHepgDispl",&_nHepgDispl,0); _ElectronNtuple->span("HEPG::nHepg",0,_MAXHEP); _ElectronNtuple->columnArray("HEPG::Hepg_M",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_M",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_M","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_ind",0,0); _ElectronNtuple->dimension("HEPG::Hepg_ind",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_ind","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_type",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_type",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_type","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_ist",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_ist",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_ist","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_moId",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_moId",_MAXHEP,_MAXMOT); _ElectronNtuple->index("HEPG::Hepg_moId","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_mo1",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_mo1",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_mo1","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_mo2",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_mo2",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_mo2","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_da1",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_da1",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_da1","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_da2",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_da2",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_da2","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Vx",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Vx",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Vx","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Vy",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Vy",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Vy","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Vz",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Vz",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Vz","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_E",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_E",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_E","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Px",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Px",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Px","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Py",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Py",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Py","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Pz",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Pz",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Pz","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Pt",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Pt",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Pt","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Phi",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Phi",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Phi","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Theta",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Theta",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Theta","nHepg"); _ElectronNtuple->columnArray("HEPG::Hepg_Eta",0.,0.); _ElectronNtuple->dimension("HEPG::Hepg_Eta",_MAXHEP); _ElectronNtuple->index("HEPG::Hepg_Eta","nHepg"); } // -------------------------------------------------------------------------------- return AppResult::OK; } AppResult HighPtElectrons::fillHistograms(AbsEvent* anEvent) { _nEvt = AbsEnv::instance()->trigNumber(); _nRun = AbsEnv::instance()->runNumber(); _mcflag = AbsEnv::instance( )->monteFlag(); if(_fill_ELEC.value()) { _nElec = 0; } if(_fill_OBS.value()) { _nElecOBS = 0; _nVert = 0; } if(_fill_HEPG.value()) { _nHepg = 0; } if(_fill_TRIG.value()) { _nTrig = 128; } _ElectronNtuple->clearDataBlock("GLB"); if(_fill_TRIG.value()) _ElectronNtuple->clearDataBlock("TRIG"); if(_fill_ELEC.value()) _ElectronNtuple->clearDataBlock("ELEC"); if(_fill_MET.value()) _ElectronNtuple->clearDataBlock("MET"); if(_fill_OBS.value()) _ElectronNtuple->clearDataBlock("OBS"); if(_fill_HEPG.value()) _ElectronNtuple->clearDataBlock("HEPG"); // -------------------------------------------------------------------------------- // TowerGeometry utility class for investigating tower coordinate boundaries etc. // -------------------------------------------------------------------------------- TowerGeometry towerGeometry; // -------------------------------------------------------------------------------- // Find the calorimeter measured missing Et : // -------------------------------------------------------------------------------- // this will return the first Met object it finds CdfMet_ch mets; CdfMet::Error metStatus = CdfMet::find(mets); if(_fill_MET.value()) { if(metStatus == CdfMet::ERROR){ if(_debug.value()) ERRLOG(ELerror,"CdfMetExampleModule: ") << "Error in CdfMetExampleModule: could not find CdfMet" << endmsg; } else{ _MetEt = mets->met(); _MetPhi = mets->phi(); } } CdfTrackColl_ch trackSet; StorableObject::SelectByDescription description("GlobalCOT_Tracking"); CdfTrackColl::Error error = CdfTrackColl::find(trackSet,description); if (error == CdfTrackColl::OK) { if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : CdfTrackColl for global COT tracking found with size = " << trackSet->contents().size() << std::endl; int iCOT_Track = 0; for (CdfTrackColl::const_iterator trackIter = trackSet->contents().begin(); trackIter != trackSet->contents().end(); ++trackIter) { iCOT_Track++; std::cout << " HighPtElectrons::fillHistograms : COT track " << iCOT_Track << " has ID = " << (*trackIter)->id() << " and pT = " << (*trackIter)->pt() << std::endl; } } } // -------------------------------------------------------------------------------- // Get the Pes2dClusterColl for this event : // -------------------------------------------------------------------------------- Pes2dClusterColl::const_handle pes2dClusterColl; Pes2dClusterColl::Error pes2dClusterCollResult = Pes2dClusterColl::find(pes2dClusterColl); if (pes2dClusterCollResult == Pes2dClusterColl::ERROR) { if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : Pes2dClusterColl does not exist" << std::endl;} _nPes2dCluster = 0; } else { _nPes2dCluster = pes2dClusterColl->contents().size(); if (_debug.value()) { std::cout << " Pes2dClusterColl size = " << pes2dClusterColl->contents().size() << std::endl; std::cout << " Loop through and print out the contents of Pes2dClusterColl : " << std::endl; } for (Pes2dClusterColl::const_iterator pes2dClusterIter = pes2dClusterColl->contents().begin(); pes2dClusterIter != pes2dClusterColl->contents().end(); ++pes2dClusterIter) { if (_debug.value()) std::cout << " New Pes2dCluster : " << std::endl; (*pes2dClusterIter).print(std::cout); } } // --------------------------------------------------------- // Find the electrons in the event and fill the ntuples : // --------------------------------------------------------- if (_fill_ELEC.value()) { CdfEmObjectColl::const_handle emObjectColl; CdfEmObjectColl::Error result = CdfEmObjectColl::find(emObjectColl); if (!result) { if (_debug.value()) std::cout << " HighPtElectrons::fillHistograms : CdfEmObject collection does not exist" << std::endl; } else { if (_debug.value()) std::cout << " HighPtElectrons::fillHistograms : number of CdfEmObjects = " << emObjectColl->contents().size() << std::endl; _nElec = emObjectColl->contents().size(); if (_debug.value()) std::cout << "HighPtElectrons : allocating _Elec arrays" << std::endl; for( int i = -1; i < _MAXELE;i++) { _ElecEtot[i] = 9999.0; _ElecZCentroid[i] = 9999.0; _ElecHadEm[i] = 9999.0; _ElecTrigHadEm[i] = 9999.0; _ElecLShwr[i] = 9999.0; _ElecLShwr2[i] = 9999.0; //non corrected quantities _ElecEmEt[i] = 9999.0; _ElecPx[i] = 9999.0; _ElecPy[i] = 9999.0; _ElecPz[i] = 9999.0; _ElecHadEt[i] = 9999.0; _ElecIsoTot4[i] = 9999.0; _ElecIsoTot7[i] = 9999.0; _ElecIsoEm4[i] = 9999.0; _ElecIsoEm7[i] = 9999.0; _ElecIsoHad4[i] = 9999.0; _ElecIsoHad7[i] = 9999.0; _ElecEta[i] = 9999.0; _ElecPhi[i] = 9999.0; _ElecTheta[i] = 9999.0; _ElecMaxTrkIso[i] = 9999.0; _ElecMaxTrkId[i] = 9999.0; _ElecMaxTrkPt[i] = 9999.0; _ElecMaxTrkPhi[i] = 9999.0; _ElecMaxTrkEta[i] = 9999.0; _ElecMaxTrkTheta[i]= 9999.0; _ElecMaxTrkD0[i] = 9999.0; _ElecMaxTrkZ0[i] = 9999.0; _ElecMaxTrkX[i] = 9999.0; _ElecMaxTrkZ[i] = 9999.0; _ElecBestCesX[i] = 9999.0; _ElecBestCesXChiW[i] = 9999.0; _ElecBestCesDeltaX[i] = 9999.0; _ElecBestCesZ[i] = 9999.0; _ElecBestCesZChiS[i] = 9999.0; _ElecBestCesDeltaZ[i] = 9999.0; _ElecNMatchTrks[i] = 9999; _ElecNCesDefCls[i] = 9999; _ElecNCesWirCls[i] = 9999; _ElecNCesStrCls[i] = 9999; _ElecNCprCls[i] = 9999.0; _ElecCprCoord[i] = 9999.0; _ElecCprCharge[i] = 9999.0; _ElecCprDelta[i] = 9999.0; } // Iterate over the collection of high level objects : int iElec = -1; for(CdfEmObjectColl::const_iterator emIter = emObjectColl->contents().begin();emIter != emObjectColl->contents().end(); ++emIter) { if (_debug.value())std::cout << " HighPtElectrons::fillHistograms : examine new CdfEmObject " << std::endl; iElec++; if(iElec <_MAXELE) { const CdfEmObject& thisCdfEmObject = *emIter; // Get the cluster pointed to by the CdfEmObject : const EmCluster* emClust = thisCdfEmObject.getEmCluster(); //Do we want the ZVertex from the track ?? double ZVert = 0.0; if(_useTrkZVert.value() && thisCdfEmObject.maxPtTrack().is_nonnull()) { // if(thisCdfEmObject.matchingTracks().contents().size() > 0) { //setting the vertex //cluster->vertexZ(); ZVert = thisCdfEmObject.maxPtTrack()->z0(); } //} emClust->setVertexZ(ZVert); _ElecReg[iElec] = emClust->region(); _ElecEtot[iElec] = emClust->totalEnergy(); _ElecEmEt[iElec] = emClust->emEt(); _ElecHadEt[iElec] = emClust->hadEt(); _ElecPx[iElec] = emClust->fourMomentum().px(); _ElecPy[iElec] = emClust->fourMomentum().py(); _ElecPz[iElec] = emClust->fourMomentum().pz(); double px = emClust->fourMomentum().px(); double py = emClust->fourMomentum().py(); double pz = emClust->fourMomentum().pz(); double pt = sqrt(px*px + py*py); // double phi = 9999.0 ; // if(py>0) phi = acos(px/pt); // if(py<0) phi = 2*M_PI-acos(px/pt); // if(py>0&&px==0) phi = 0.5*M_PI; // if(py<0&&px==0) phi = -0.5*M_PI; // if(py==0&&px==0)phi = 0.0; _ElecPhi[iElec] = emClust->emPhi(); double theta = acos(pz/sqrt(px*px + py*py + pz*pz)); double eta = log(1/tan(0.5*theta)); _ElecEta[iElec] = emClust->emEtaEvent(); std::cout <<"Just a check: etaEv, etaDet at z=0 are"<emEtaEvent()<<" , "<emEtaDetector()<fourMomentum().theta(); _ElecZCentroid[iElec]= emClust->zCentroid(); _ElecHadEm[iElec] = emClust->hadEm(); _ElecTrigHadEm[iElec]= emClust->hadEmTriggerTower(); _ElecIsoTot4[iElec] = emClust->totalIsolationEt4(); _ElecIsoTot7[iElec] = emClust->totalIsolationEt7(); _ElecIsoEm4[iElec] = emClust->emIsolationEt4(); _ElecIsoEm7[iElec] = emClust->emIsolationEt7(); _ElecIsoHad4[iElec] = emClust->hadIsolationEt4(); _ElecIsoHad7[iElec] = emClust->hadIsolationEt7(); _ElecLShwr[iElec] = thisCdfEmObject.lshr(); _ElecLShwr2[iElec] = thisCdfEmObject.lshr(2); if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : theta, eta = " << _ElecTheta[iElec] << " , " << _ElecEta[iElec] << std::endl; } if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : electron has " << thisCdfEmObject.matchingTracks().contents().size() << " matching tracks." << std::endl; } int nMatchingTracks = thisCdfEmObject.matchingTracks().contents().size(); _ElecNMatchTrks[iElec] = nMatchingTracks ; int itrk = -1; for(CdfTrackView::const_iterator iter = thisCdfEmObject.matchingTracks().contents().begin(); iter != thisCdfEmObject.matchingTracks().contents().end() ; iter++) { if ((int)*iter == NULL) continue; const CdfTrack & trkIter = **iter; itrk++; if(itrk < _MAXTRK) { _ElecTrkPt[iElec][itrk] = trkIter.pt(); _ElecTrkPhi[iElec][itrk]= trkIter.phi0(); _ElecTrkEta[iElec][itrk]= trkIter.pseudoRapidity(); _ElecTrkD0[iElec][itrk] = trkIter.d0(); _ElecTrkZ0[iElec][itrk] = trkIter.z0(); _ElecTrkCot[iElec][itrk] = trkIter.lambda(); _ElecTrkCurv[iElec][itrk] = trkIter.curvature(); _ElecTrkAxHits[iElec][itrk] = trkIter.numCTHitsAx(); _ElecTrkStHits[iElec][itrk] = trkIter.numCTHitsSt(); _ElecTrkAlg[iElec][itrk] = float(trkIter.algorithm().value()); _ElecTrkId[iElec][itrk] = 0;//trkIter.object_id().value(); //trying to extrapolate the track to the CES radius to get the coordinate; StripCollectionMaker strip_coll; CdfTrack_clnk cdfTrack=(*iter); strip_coll.extrapolate_track (cdfTrack,_mcflag ); // Extrapolation _ElecTrkXAtCes[iElec][itrk] = strip_coll.get_Xlocal(); _ElecTrkZAtCes[iElec][itrk] = strip_coll.get_Zlocal(); _ElecTrkEtaAtCes[iElec][itrk]= strip_coll.Eta_global(); _ElecTrkPhiAtCes[iElec][itrk] = strip_coll.Phi_global(); _ElecTrkSideOfCes[iElec][itrk]= strip_coll.getSide();//barrel of the cluster,0=west, 1 = east _ElecTrkCesWedge[iElec][itrk] = strip_coll.getModule();//phi wedge of the cluster } } if(_makeHistos.value()) _nMatchTrks->accumulate(nMatchingTracks); //if (thisCdfEmObject.matchingTracks().contents().size() > 0) if ( thisCdfEmObject.maxPtTrack().is_nonnull() ) { if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : highest matching track pT = " << thisCdfEmObject.maxPtTrack()->pt() << std::endl; } _ElecMaxTrkId[iElec] = 0;//thisCdfEmObject.maxPtTrack().object_id().value(); _ElecMaxTrkPt[iElec] = thisCdfEmObject.maxPtTrack()->pt(); _ElecMaxTrkPhi[iElec] = thisCdfEmObject.maxPtTrack()->phi0(); //this is actually the phi0 of the input track to SimpleExtrapolatedTrack _ElecMaxTrkEta[iElec] = thisCdfEmObject.maxPtTrack()->pseudoRapidity(); _ElecMaxTrkTheta[iElec] = thisCdfEmObject.maxPtTrack()->theta(); _ElecMaxTrkZ0[iElec] = thisCdfEmObject.maxPtTrack()->z0(); _ElecMaxTrkD0[iElec] = thisCdfEmObject.maxPtTrack()->d0(); _ElecMaxTrkX[iElec] = thisCdfEmObject.maxPtTrackLocalCoord().x(); _ElecMaxTrkZ[iElec] = thisCdfEmObject.maxPtTrackLocalCoord().z(); // Track Isolation Selection float coneSize = 0.4; float total_pt = 0.0; CdfTrackView_h view_hndl; CdfTrackView::Error trackStatus = CdfTrackView::defTracks(view_hndl); if (trackStatus == CdfTrackView::OK) { CdfTrackView* coneTracks = new CdfTrackView(thisCdfEmObject.matchingTracks().contents().begin(), thisCdfEmObject.matchingTracks().contents().end(), track_cut::InConeAround(thisCdfEmObject.maxPtTrack()->momentum(),coneSize) ); for (CdfTrackView::iterator coneIterator = coneTracks->contents().begin(); coneIterator != coneTracks->contents().end(); coneIterator++) { const CdfTrack* compare_track = coneIterator->operator->(); total_pt += compare_track->pt(); } } float TrackPt = thisCdfEmObject.maxPtTrack()->pt(); if (total_pt > 0.0&&TrackPt!=0) _ElecMaxTrkIso[iElec] = (total_pt - TrackPt)/total_pt; // } ValueVectorLink strip_cluster; ValueVectorLink wire_cluster; double XCes = 9999.0; double ChiW = 9999.0; double DeltaX = 9999.0; wire_cluster = thisCdfEmObject.bestMatchingCesCluster(CesCluster::WIRE); if (wire_cluster.is_nonnull()) { XCes = wire_cluster->fitted_position(); //ChiW = wire_cluster->chi2_strip(); ChiW = scaleCesChisq(wire_cluster->chi2_strip() , emClust->totalEnergy()); if( thisCdfEmObject.maxPtTrack().is_nonnull() )DeltaX = thisCdfEmObject.maxPtTrackLocalCoord().x() - wire_cluster->fitted_position(); } _ElecBestCesDeltaX[iElec] = DeltaX; _ElecBestCesX[iElec] = XCes; _ElecBestCesXChiW[iElec] = ChiW; double ZCes = 9999.0; double ZCesGlob = 9999.0; double ChiS = 9999.0; double DeltaZ = 9999.0; strip_cluster = thisCdfEmObject.bestMatchingCesCluster(CesCluster::STRIP); if (strip_cluster.is_nonnull()) { //HepPoint3D globalPos = strip_cluster->global_position(); //ZCesGlob = globalPos.z(); ZCes = strip_cluster->fitted_position(); // ZCesGlob = strip_cluster->globalCoord();this should work in the next version int side = strip_cluster->side(); if(side==0)ZCesGlob=-ZCes; if(side==1)ZCesGlob=ZCes; ChiS = scaleCesChisq(strip_cluster->chi2_strip(), emClust->totalEnergy()); if( thisCdfEmObject.maxPtTrack().is_nonnull() )DeltaZ = thisCdfEmObject.maxPtTrackLocalCoord().z() - ZCesGlob; } _ElecBestCesDeltaZ[iElec] = DeltaZ; _ElecBestCesZ[iElec] = ZCesGlob; _ElecBestCesZChiS[iElec] = ChiS; if (_debug.value()) { std::cout << "Xces= "< _MAXCES) continue; _ElecCesView[iElec][ices] = 9999.0; _ElecCesLocCoord[iElec][ices] = 9999.0; _ElecCesGlbCoord[iElec][ices] = 9999.0; _ElecCesEnergy[iElec][ices] = 9999.0; _ElecCesChi2[iElec][ices] = 9999.0; _ElecCesDistTrkCls[iElec][ices] = 9999.0; _ElecCesAssTrkId[iElec][ices] = 9999.0; if ((*ciIndxLinksCes)->view() == CesCluster::WIRE) { nCesWire++ ; _ElecCesView[iElec][ices] = 0; _ElecCesLocCoord[iElec][ices] = (*ciIndxLinksCes)->localCoord();//fitted_position(); //Z for strips,X for wires _ElecCesGlbCoord[iElec][ices] = (*ciIndxLinksCes)->globalCoord();//global_position(); //Z for strips,Phi for wires _ElecCesEnergy[iElec][ices] = (*ciIndxLinksCes)->raw_energy();//fitted_energy() is always equal to zero... _ElecCesChi2[iElec][ices] = (*ciIndxLinksCes)->chi2(); //const CdfTrack & theTrack() const; if((*ciIndxLinksCes)->theTrack().is_nonnull()) { _ElecCesAssTrkId[iElec][ices] = 0;// (*ciIndxLinksCes)->theTrack().object_id().value(); _ElecCesDistTrkCls[iElec][ices] = (*ciIndxLinksCes)->dist_track_cluster(); } else{ _ElecCesAssTrkId[iElec][ices] = 9999.0; _ElecCesDistTrkCls[iElec][ices] = 9999.0; } } if ((*ciIndxLinksCes)->view() == CesCluster::STRIP) { nCesStrip++ ; _ElecCesView[iElec][ices] = 1; int side = (*ciIndxLinksCes)->side();//0 = west,1 = east; if(side==0) _ElecCesGlbCoord[iElec][ices] = -( (*ciIndxLinksCes)->global_position()); else _ElecCesGlbCoord[iElec][ices] = (*ciIndxLinksCes)->global_position(); _ElecCesLocCoord[iElec][ices] = (*ciIndxLinksCes)->fitted_position(); _ElecCesEnergy[iElec][ices] = (*ciIndxLinksCes)->fitted_energy(); _ElecCesChi2[iElec][ices] = (*ciIndxLinksCes)->chi2(); if((*ciIndxLinksCes)->theTrack().is_nonnull()) { _ElecCesDistTrkCls[iElec][ices] = (*ciIndxLinksCes)->dist_track_cluster(); _ElecCesAssTrkId[iElec][ices] = 0;//(*ciIndxLinksCes)->theTrack().object_id().value(); } else{ _ElecCesAssTrkId[iElec][ices] = 9999.0; _ElecCesDistTrkCls[iElec][ices] = 9999.0; } } } nShwrMaxClusters = thisCdfEmObject.matchingCesClusters("DEFAULT").contents().size(); if(_makeHistos.value()) { nDefWireCesCls = nCesWire; nDefStripCesCls = nCesStrip; _nCesDefWireCls->accumulate(nDefWireCesCls); _nCesDefStripCls->accumulate(nDefStripCesCls); } } else if (emClust->region() == PEM) { nShwrMaxClusters = thisCdfEmObject.matchingPes2dClusters().contents().size(); } if(_makeHistos.value())_nCesCls->accumulate(nShwrMaxClusters); _ElecNCesDefCls[iElec] = nShwrMaxClusters ; _ElecNCesWirCls[iElec] = nCesWire ; _ElecNCesStrCls[iElec] = nCesStrip ; //CPR informations if(_makeHistos.value()) _nCprDefCls->accumulate(thisCdfEmObject.matchingCprClusters("DEFAULT").contents().size()); _ElecNCprCls[iElec] = thisCdfEmObject.matchingCprClusters("DEFAULT").contents().size(); if (_debug.value()) { std::cout << "n cpr clusters = " <::const_iterator ciIndxLinksCpr = thisCdfEmObject.matchingCprClusters( "DEFAULT" ).contents().begin(); std::vector::const_iterator lastIndxLink = thisCdfEmObject.matchingCprClusters( "DEFAULT" ).contents().end(); for ( ; ciIndxLinksCpr != lastIndxLink; ++ciIndxLinksCpr) { cprCoord = (*ciIndxLinksCpr)->localCoord(); cprCharge= (*ciIndxLinksCpr)->energy(); // extrapolate a track : CprWireCollectionMaker wireCol; if( thisCdfEmObject.maxPtTrack().is_nonnull() ) { wireCol.extrapolate_track(thisCdfEmObject.maxPtTrack()); float Xloc = wireCol.get_Xlocal(); //float Zloc = wireCol.get_Zlocal(); diff = Xloc - (*ciIndxLinksCpr)->localCoord(); if (_debug.value()) std::cout << "diff= "< M_PI) deltaPhi = 2*M_PI - deltaPhi; etaPhiArea = etaPhiArea + deltaEta*deltaPhi; if (towerGeometry.etaMin(thisTowerKey) < etaMin) etaMin = towerGeometry.etaMin(thisTowerKey); if (towerGeometry.etaMax(thisTowerKey) > etaMax) etaMax = towerGeometry.etaMax(thisTowerKey); } if (_debug.value()) { std::cout << " HighPtElectrons::fillHistograms : numTowers = " << numTowers << std::endl; std::cout << " HighPtElectrons::fillHistograms : eta size from towers = " << etaMax-etaMin << std::endl; std::cout << " HighPtElectrons::fillHistograms : eta RMS from cluster = " << emClust->emEtaRMS() << std::endl; std::cout << " HighPtElectrons::fillHistograms : eta-phi area from towers = " << etaPhiArea << std::endl; std::cout << " HighPtElectrons::fillHistograms : eta-phi area from RMS = " << M_PI*emClust->emEtaRMS()*emClust->emPhiRMS() << std::endl; } // Fill certain histograms for all clusters : if(_makeHistos.value()){ _ClusterEtot->accumulate(emClust->totalEnergy()); _ClusterE_Em->accumulate(emClust->emEnergy()); _ClusterHadEm->accumulate(emClust->hadEm()); _ClusterNumTow->accumulate(numTowers); _ClusterEtaPhiArea->accumulate(etaPhiArea); } if (emClust->emEt()*_EM_EnergyScaleFactor.value() > _EM_ClusterPtCut.value() && fabs(emClust->emEtaEvent()) < _EM_ClusterAbsEtaRange.value() && // Cut on HAD/EM : (emClust->hadEm() < _ClusterHadEmCut.value())) { // Fill cluster Pt histogram : if(_makeHistos.value()){ _EM_ClusterPt->accumulate(emClust->emEt()*_EM_EnergyScaleFactor.value()); _EM_ClusterEta->accumulate(emClust->emEtaEvent()); _EM_ClusterPhi->accumulate(emClust->emPhi()); } // A pointer to the highest pT matching track, set to NULL : const CdfTrack* highestPtMatchingTrack = 0; // Keep track of the matching track pT values : double highestPt = 0.0; const CdfTrackView& matchedTracks = thisCdfEmObject.matchingTracks(); if (_debug.value()) std::cout << " HighPtElectrons::fillHistograms : size of matching track set = " << matchedTracks.contents().size() << std::endl; for (CdfTrackView::const_iterator trkIter = matchedTracks.contents().begin(); trkIter != matchedTracks.contents().end(); ++trkIter) { double matchingTrackPt = (*trkIter)->pt(); if (_debug.value()) std::cout << " HighPtElectrons::fillHistograms : matching track pT = " << matchingTrackPt << std::endl; // Only consider the track if its pT is higher than previous matching tracks : if ( matchingTrackPt > _MatchingTrackPtCut.value() && matchingTrackPt > highestPt ) { highestPtMatchingTrack = (*trkIter).operator->(); highestPt = matchingTrackPt; } } // Check to see if a high pT matching track was found : if (highestPtMatchingTrack) { // Fill matching track histograms : if(_makeHistos.value()){ _MatchingTrackPt->accumulate(highestPtMatchingTrack->pt()); _MatchingTrackEta->accumulate(log(tan((atan(highestPtMatchingTrack->lambda())+M_PI/2.0)/2.0))); _MatchingTrackPhi->accumulate(highestPtMatchingTrack->phi0()); _EoverP->accumulate(emClust->emEt()*_EM_EnergyScaleFactor.value()/highestPtMatchingTrack->pt()); } // This is a candidate electron, so calculate a transverse mass : if(!metStatus == CdfMet::ERROR){ HepVector electron(2), neutrino(2); electron[0] = emClust->fourMomentum()[0]; electron[1] = emClust->fourMomentum()[1]; neutrino[0] = -mets->exSum(); neutrino[1] = -mets->eySum(); double mT = transverseMass(electron, neutrino); // Fill transverse mass histogram : if(_makeHistos.value())_TransverseMass->accumulate(mT); } // Store this emClust in a local array to form a multi-electron invariant mass distribution : } // ---------------------------------------------------------------------------------- // Find a matching track for this cluster using a separate track matching algorithm : // CT_Track* nearestTrack = findMatchingTrackForCluster(*emClust); // if (nearestTrack&&_makeHistos.value()) // { // std::cout << " DW matching track : pT = " << nearestTrack->pt() << std::endl; // _MatchingTrackPt->accumulate(nearestTrack->pt()); // _EoverP->accumulate(emClust->calEtEm()/nearestTrack->pt()); // } // Fill an invariant mass histogram in the case of exactly two clusters : if (emObjectColl->contents().size() == 2) { CdfEmObjectColl::const_iterator EmObjectIter = emObjectColl->contents().begin(); const EmCluster* pEmCluster1 = (*EmObjectIter).getEmCluster(); ++EmObjectIter; const EmCluster* pEmCluster2 = (*EmObjectIter).getEmCluster(); HepLorentzVector fourMom1 = pEmCluster1->fourMomentum(); HepLorentzVector fourMom2 = pEmCluster2->fourMomentum(); HepLorentzVector fourMomSum = fourMom1 + fourMom2; double iMass = fourMomSum.mag(); if(_makeHistos.value()) _InvariantMass->accumulate( iMass ); } } } } if (_regionalHistos.value()) { // Make some regional histograms. Intended at this point for single particle events only. // First, find some generator level information : const ParentedParticle* particleOfInterest = NULL; TRYGenPrimaryVertexSet myGenSet(*anEvent); myGenSet.initializeFromObs(); for (GenPrimaryVertexSet::ConstIterator vtx=myGenSet.begin();vtx!=myGenSet.end();vtx++) { for (GenPrimaryVertex::ConstParticleIterator part=(*vtx).begin();part!=(*vtx).end();part++) { bool electron = (labs((*part).type()) == 11); const SelectionCriterion &Prompt = IsPrompt(); if (electron && Prompt(*part)) { particleOfInterest = &(*part); } } } if (particleOfInterest) { double genEta = particleOfInterest->momentum().pseudoRapidity(); // Now loop through the regions : StorableObject::SelectByClassName class_name("GrowResultList"); for (EventRecord::ConstIterator GrowResultListIter(anEvent,class_name); GrowResultListIter != anEvent->end(); GrowResultListIter++) { GrowResultList_ch GrowResultListCH(GrowResultListIter); // DSW 23.08.2002. TrackingHL/Regions REMOVED ! // // Want to only consider those regions for which the seed matches the electron. Technique is // // similar to that applied in benchmarking code : // const SelectionCriterion& matchesSeed = DoesMatchSeed(GrowResultListCH.operator*(), // // MC-COT matching chi-squared cut (irrelevant) : // 0.0, // // MC-Plug matching DCA cut : // 5.0, // // Verbose level : // 0); // if (matchesSeed(*particleOfInterest)) // { // // Find the cluster data : // // GrowResultList::typeIterator siClusterDataIter(GrowResultListCH); // // if (siClusterDataIter.is_valid()&&_makeHistos.value()) // // { // // _nClustersVsEta->accumulate(genEta,(*siClusterDataIter)->clusterSet()->nEntries()); // // } // } } } } } } //Trigger information if(_fill_TRIG.value()) { for (int i = 0; i < 128 ; i++) { _L1Bit[i] = 9999; _L2Bit[i] = 9999; _L3Bit[i] = 9999; } // Get L1 and L2 trigger decision bits long bitmask[32]; long tl2d_l2trig[4],tfrd_l1trig[2]; for(int i = 0; i < 32; i++) { bitmask[i] = 1 << i; //std::cout << std::hex << i << "," << bitmask[i] << std::endl; } EventRecord::ConstIterator TFRD_iter( anEvent, "TFRD_StorableBank" ); if ( TFRD_iter.is_valid() ) { ConstHandle tfrd( TFRD_iter ); // Get the L1 decision from FRED. for(int iword = 0; iword < 2; iword++) { tfrd_l1trig[iword] = tfrd->getL1PrescaleTrigDecision(iword*32); } if(_debug.value()) std::cout << " Looking at L1 bits from FRED:" ; for (int iword = 0 ; iword < 2; iword++){ for(int ibit = 0; ibit < 32; ibit++) { if (ibit%16 ==0) if(_debug.value()) std::cout << std::endl << (ibit+(iword)*32) << ": "; if(_debug.value()) std::cout << int(1&&( tfrd_l1trig[iword]&bitmask[ibit] )) << " "; _L1Bit[ibit+(iword)*32] = int(1&&( tfrd_l1trig[iword]&bitmask[ibit] )); } } } EventRecord::ConstIterator TL2D_iter( anEvent, "TL2D_StorableBank" ); if ( TL2D_iter.is_valid() ) { ConstHandle tl2d( TL2D_iter ); // Get the L2 prescaled decision bits for(int iword = 0; iword < 4; iword++) { tl2d_l2trig[iword] = tl2d->getL2TriggerBit(iword*32); } if(_debug.value()) std::cout << std::endl << " Looking at the L2 prescaled bits:"; for (int iword = 0 ; iword < 4; iword++){ for(int ibit = 0; ibit < 32; ibit++) { // prescaled decision [127:0], words 1-4 if (ibit%16 ==0) if(_debug.value()) std::cout << std::endl << (ibit+iword*32) << ": "; if(_debug.value()) std::cout << int(1&&( tl2d_l2trig[iword]&bitmask[ibit] )) << " "; _L2Bit[ibit+iword*32] = int(1&&( tl2d_l2trig[iword]&bitmask[ibit] )); } } } EventRecord::ConstIterator TL3D_iter( anEvent, "TL3D_StorableBank" ); if ( TL3D_iter.is_valid() ) { ConstHandle tl3d( TL3D_iter ); if(_debug.value()) std::cout << endl << "looking at L3 bits:" << std::endl; for ( int ibit = 0; (ibit < tl3d->nPaths()) ; ibit++ ) { if(_debug.value()) std::cout << "Mod #" << std::dec << ibit << " " << tl3d->didPathPass(ibit) << std::endl; _L3Bit[ibit] = tl3d->didPathPass(ibit); } } } if(_fill_OBS.value()) { if (_debug.value()) std::cout << "HighPtElectrons : allocating _ElecOBSV arrays" << std::endl; for( int i = 0; i < _MAXVER; i++) { _ElecOBSV_Vx[i] = 9999.0; _ElecOBSV_Vy[i] = 9999.0; _ElecOBSV_Vz[i] = 9999.0; _ElecOBSV_nDau[i] = 9999.0; // OBSP informations from OBSP Bank if (_debug.value()) std::cout << "HighPtElectrons : allocating _ElecOBSP arrays" << std::endl; for( int j = 0; j < _MAXOBS; j++) { _ElecOBSP_Px[i][j]= 9999.0; _ElecOBSP_Py[i][j]= 9999.0; _ElecOBSP_Pz[i][j]= 9999.0; _ElecOBSP_Ptot[i][j]= 9999.0; _ElecOBSP_Pt[i][j] = 9999.0; _ElecOBSP_Theta[i][j]= 9999.0; _ElecOBSP_Eta[i][j]= 9999.0; _ElecOBSP_Phi[i][j]= 9999.0; _ElecOBSP_Curv[i][j]= 9999.0; _ElecOBSP_globCprX[i][j]= 9999.0; _ElecOBSP_cprX[i][j]= 9999.0; _ElecOBSP_cprY[i][j]= 9999.0; _ElecOBSP_cprZ[i][j]= 9999.0; _ElecOBSP_cprPhi[i][j]= 9999.0; _ElecOBSP_cprTheta[i][j]= 9999.0; _ElecOBSP_cprEta[i][j]= 9999.0; } } EventRecord::ConstIterator obsvIter(anEvent,"OBSV_StorableBank"); EventRecord::ConstIterator obspIter(anEvent,"OBSP_StorableBank"); while(obsvIter.is_valid()) { ConstHandle obsvCH(obsvIter); const OBSV_StorableBank& obsv = obsvCH.operator*(); int iVERT=-1; for(OBSV_StorableBank::vertexIter vx(obsv);vx.is_valid();++vx) { ++iVERT; if(iVERT<_MAXVER) { _ElecOBSV_Vx[iVERT] = obsv.X(vx); _ElecOBSV_Vy[iVERT] = obsv.Y(vx); _ElecOBSV_Vz[iVERT] = obsv.Z(vx); _ElecOBSV_nDau[iVERT] = obsv.NDaughters(vx); int iOBSP = 0; while(obspIter.is_valid()) { ConstHandle obspCH(obspIter); const OBSP_StorableBank& obsp = obspCH.operator*(); for (OBSP_StorableBank::particleIter i(obsp); i.is_valid(); ++i) { if(iOBSP < _MAXOBS ) { double ptot = sqrt(obsp.Px(i)*obsp.Px(i) + obsp.Py(i)*obsp.Py(i)+ obsp.Pz(i)*obsp.Pz(i)); double pt = sqrt(obsp.Px(i)*obsp.Px(i) + obsp.Py(i)*obsp.Py(i)); //double cot = obsp.Pz(i)/pt; double theta = 9999.0; double cot = 9999.0; if(ptot!=0) { theta = acos(obsp.Pz(i)/ptot); if(abs(theta)>_myMin) cot = 1/tan(theta); } if(abs(tan(0.5*theta))>= _myMin)_ElecOBSP_Eta[iVERT][iOBSP] = log(1/tan(0.5*theta)); double px=obsp.Px(i); double py=obsp.Py(i); double pz=obsp.Pz(i); double phi0; if(py>0) phi0 = acos(px/pt); if(py<0) phi0 = 2*M_PI-acos(px/pt); if(py>0&&px==0) phi0 = 0.5*M_PI; if(py<0&&px==0) phi0 = -0.5*M_PI; if(py==0&&px==0)phi0 = 0.0; double curv = ((0.001489)*(1.4116))/pt; // B (in Tesla) = 1.4116 Tesla double z0=obsv.Z(vx); double d0=0.0; Helix myhelix(cot,curv,z0,d0,phi0); // if (_debug.value()) std::cout << "helix = "<< myhelix.getParameters() << std::endl; SimpleCprPathFinder finder; const CprIntersection *intersect; if(abs(phi0)>_myMin&&abs(phi0)<_myPi&&abs(theta)>_myMin&&abs(theta)<_myPi) intersect = finder.newIntersection(myhelix); else intersect = NULL; float z_local = 9999.0; float x_local = 9999.0; float y_local = 9999.0; float x_global = 9999.0; float theta_global = 9999.0; float eta_global = 9999.0; float phi_global = 9999.0; if(intersect) { HepPoint3D globalPos = intersect->globalPosition(); x_local = intersect->localX(); z_local = intersect->localZ(); y_local = globalPos.y(); x_global = globalPos.x(); if (_debug.value()) { std::cout << "Checking CPR: is GlobPos.z " < chHEPG(hepgIter); const HEPG_StorableBank& hepg = chHEPG.operator*(); _nHepg = min(_MAXHEP,hepg.n_particles()); for (HEPG_StorableBank::particleIter partIter(hepg); partIter.is_valid(); ++partIter) { int iHEPG = partIter.index(); float hepg_type = hepg.idhep(partIter); //fill up this array only for electrons if(abs(hepg_type)==11) { int mo_index = hepg.jmo1hep(partIter); for(int i = 0;i< _MAXMOT;i++) { double mo_type = hepg.idhep(mo_index-1); _Hepg_moId[iHEPG][i]= mo_type; //doing only in case the particle is repeated ONCE; if(mo_type== hepg.idhep( hepg.jmo1hep(mo_index-1)-1))mo_type= hepg.idhep(mo_index-1); if(_debug.value()) { std::cout << "Particle #" <capture("GLB::nRun",_nRun); _ElectronNtuple->captureBlock("GLB"); if(_fill_TRIG.value()&&_nTrig!=0) { _ElectronNtuple->capture("TRIG::nTrig",_nTrig); _ElectronNtuple->capture("TRIG::L1Bit",&(_L1Bit[0])); _ElectronNtuple->capture("TRIG::L2Bit",&(_L2Bit[0])); _ElectronNtuple->capture("TRIG::L3Bit",&(_L3Bit[0])); _ElectronNtuple->captureBlock("TRIG"); } if(_fill_ELEC.value()&&_nElec!=0) { _ElectronNtuple->capture("ELEC::nElec",_nElec); _ElectronNtuple->capture("ELEC::ElecReg",&(_ElecReg[0])); _ElectronNtuple->capture("ELEC::ElecEtot",&(_ElecEtot[0])); _ElectronNtuple->capture("ELEC::ElecEmEt",&(_ElecEmEt[0])); _ElectronNtuple->capture("ELEC::ElecHadEt",&(_ElecHadEt[0])); _ElectronNtuple->capture("ELEC::ElecHadEm",&(_ElecHadEm[0])); _ElectronNtuple->capture("ELEC::ElecTrigHadEm",&(_ElecTrigHadEm[0])); _ElectronNtuple->capture("ELEC::ElecZCentroid",&(_ElecZCentroid[0])); _ElectronNtuple->capture("ELEC::ElecIsoTot4",&(_ElecIsoTot4[0])); _ElectronNtuple->capture("ELEC::ElecIsoTot7",&(_ElecIsoTot7[0])); _ElectronNtuple->capture("ELEC::ElecIsoEm4",&(_ElecIsoEm4[0])); _ElectronNtuple->capture("ELEC::ElecIsoEm7",&(_ElecIsoEm7[0])); _ElectronNtuple->capture("ELEC::ElecIsoHad4",&(_ElecIsoHad4[0])); _ElectronNtuple->capture("ELEC::ElecIsoHad7",&(_ElecIsoHad7[0])); _ElectronNtuple->capture("ELEC::ElecLShwr",&(_ElecLShwr[0])); _ElectronNtuple->capture("ELEC::ElecLShwr2",&(_ElecLShwr2[0])); _ElectronNtuple->capture("ELEC::ElecPx",&(_ElecPx[0])); _ElectronNtuple->capture("ELEC::ElecPy",&(_ElecPy[0])); _ElectronNtuple->capture("ELEC::ElecPz",&(_ElecPz[0])); _ElectronNtuple->capture("ELEC::ElecEta",&(_ElecEta[0])); _ElectronNtuple->capture("ELEC::ElecPhi",&(_ElecPhi[0])); _ElectronNtuple->capture("ELEC::ElecTheta",&(_ElecTheta[0])); _ElectronNtuple->capture("ELEC::ElecTrkPt",&(_ElecTrkPt[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkPhi",&(_ElecTrkPhi[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkEta",&(_ElecTrkEta[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkD0",&(_ElecTrkD0[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkZ0",&(_ElecTrkZ0[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkCot",&(_ElecTrkCot[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkCurv",&(_ElecTrkCurv[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkAxHits",&(_ElecTrkAxHits[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkStHits",&(_ElecTrkStHits[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkAlg",&(_ElecTrkAlg[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkId",&(_ElecTrkId[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkXAtCes",&(_ElecTrkXAtCes[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkZAtCes",&(_ElecTrkZAtCes[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkEtaAtCes",&(_ElecTrkEtaAtCes[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkPhiAtCes",&(_ElecTrkPhiAtCes[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkSideOfCes",&(_ElecTrkSideOfCes[0][0])); _ElectronNtuple->capture("ELEC::ElecTrkCesWedge",&(_ElecTrkCesWedge[0][0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkId",&(_ElecMaxTrkId[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkIso",&(_ElecMaxTrkIso[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkPt",&(_ElecMaxTrkPt[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkPhi",&(_ElecMaxTrkPhi[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkEta",&(_ElecMaxTrkEta[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkTheta",&(_ElecMaxTrkTheta[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkZ0",&(_ElecMaxTrkZ0[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkD0",&(_ElecMaxTrkD0[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkX",&(_ElecMaxTrkX[0])); _ElectronNtuple->capture("ELEC::ElecMaxTrkZ",&(_ElecMaxTrkZ[0])); _ElectronNtuple->capture("ELEC::ElecBestCesX",&(_ElecBestCesX[0])); _ElectronNtuple->capture("ELEC::ElecBestCesXChiW",&(_ElecBestCesXChiW[0])); _ElectronNtuple->capture("ELEC::ElecBestCesDeltaX",&(_ElecBestCesDeltaX[0])); _ElectronNtuple->capture("ELEC::ElecBestCesZ",&(_ElecBestCesZ[0])); _ElectronNtuple->capture("ELEC::ElecBestCesZChiS",&(_ElecBestCesZChiS[0])); _ElectronNtuple->capture("ELEC::ElecBestCesDeltaZ",&(_ElecBestCesDeltaZ[0])); _ElectronNtuple->capture("ELEC::ElecNCesDefCls",&(_ElecNCesDefCls[0])); _ElectronNtuple->capture("ELEC::ElecNCesWirCls",&(_ElecNCesWirCls[0])); _ElectronNtuple->capture("ELEC::ElecNCesStrCls",&(_ElecNCesStrCls[0])); _ElectronNtuple->capture("ELEC::ElecNMatchTrks",&(_ElecNMatchTrks[0])); _ElectronNtuple->capture("ELEC::ElecCesView",&(_ElecCesView[0][0])); _ElectronNtuple->capture("ELEC::ElecCesLocCoord",&(_ElecCesLocCoord[0][0])); _ElectronNtuple->capture("ELEC::ElecCesGlbCoord",&(_ElecCesGlbCoord[0][0])); _ElectronNtuple->capture("ELEC::ElecCesEnergy",&(_ElecCesEnergy[0][0])); _ElectronNtuple->capture("ELEC::ElecCesChi2",&(_ElecCesChi2[0][0])); _ElectronNtuple->capture("ELEC::ElecCesDistTrkCls",&(_ElecCesDistTrkCls[0][0])); _ElectronNtuple->capture("ELEC::ElecCesAssTrkId",&(_ElecCesAssTrkId[0][0])); _ElectronNtuple->capture("ELEC::ElecNCprCls", &(_ElecNCprCls[0])); _ElectronNtuple->capture("ELEC::ElecCprCoord", &(_ElecCprCoord[0])); _ElectronNtuple->capture("ELEC::ElecCprCharge",&(_ElecCprCharge[0])); _ElectronNtuple->capture("ELEC::ElecCprDelta",&(_ElecCprDelta[0])); _ElectronNtuple->captureBlock("ELEC"); } if(_fill_MET.value()&&_MetEt) { _ElectronNtuple->capture("MET::MetEt",_MetEt); _ElectronNtuple->capture("MET::MetPhi",_MetPhi); _ElectronNtuple->captureBlock("MET"); } if(_fill_OBS.value()) { _ElectronNtuple->capture("OBS::nElecOBS",_nElecOBS); _ElectronNtuple->capture("OBS::ElecOBSP_Px",&(_ElecOBSP_Px[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Py",&(_ElecOBSP_Py[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Pz",&(_ElecOBSP_Pz[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Ptot",&(_ElecOBSP_Ptot[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Pt",&(_ElecOBSP_Pt[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Theta",&(_ElecOBSP_Theta[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Eta",&(_ElecOBSP_Eta[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Phi",&(_ElecOBSP_Phi[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_Curv",&(_ElecOBSP_Curv[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_globCprX",&(_ElecOBSP_globCprX[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprX",&(_ElecOBSP_cprX[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprY",&(_ElecOBSP_cprY[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprZ",&(_ElecOBSP_cprZ[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprPhi",&(_ElecOBSP_cprPhi[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprTheta",&(_ElecOBSP_cprTheta[0][0])); _ElectronNtuple->capture("OBS::ElecOBSP_cprEta",&(_ElecOBSP_cprEta[0][0])); _ElectronNtuple->captureBlock("OBS"); _ElectronNtuple->capture("OBSV::nVert",_nVert); _ElectronNtuple->capture("OBSV::ElecOBSV_Vx",&(_ElecOBSV_Vx[0])); _ElectronNtuple->capture("OBSV::ElecOBSV_Vy",&(_ElecOBSV_Vy[0])); _ElectronNtuple->capture("OBSV::ElecOBSV_Vz",&(_ElecOBSV_Vz[0])); _ElectronNtuple->capture("OBSV::ElecOBSV_nDau",&(_ElecOBSV_nDau[0])); _ElectronNtuple->captureBlock("OBSV"); } if(_fill_HEPG.value()) { _ElectronNtuple->capture("HEPG::nHepgEvt",_nHepgEvt); _ElectronNtuple->capture("HEPG::nHepgPrim",_nHepgPrim); _ElectronNtuple->capture("HEPG::nHepgDispl",_nHepgDispl); _ElectronNtuple->capture("HEPG::nHepg",_nHepg); _ElectronNtuple->capture("HEPG::Hepg_ind",&(_Hepg_ind[0])); _ElectronNtuple->capture("HEPG::Hepg_type",&(_Hepg_type[0])); _ElectronNtuple->capture("HEPG::Hepg_ist",&(_Hepg_ist[0])); _ElectronNtuple->capture("HEPG::Hepg_mo1",&(_Hepg_mo1[0])); _ElectronNtuple->capture("HEPG::Hepg_mo2",&(_Hepg_mo2[0])); _ElectronNtuple->capture("HEPG::Hepg_moId",&(_Hepg_moId[0][0])); _ElectronNtuple->capture("HEPG::Hepg_da1",&(_Hepg_da1[0])); _ElectronNtuple->capture("HEPG::Hepg_da2",&(_Hepg_da2[0])); _ElectronNtuple->capture("HEPG::Hepg_M",&(_Hepg_M[0])); _ElectronNtuple->capture("HEPG::Hepg_Vx",&(_Hepg_Vx[0])); _ElectronNtuple->capture("HEPG::Hepg_Vy",&(_Hepg_Vy[0])); _ElectronNtuple->capture("HEPG::Hepg_Vz",&(_Hepg_Vz[0])); _ElectronNtuple->capture("HEPG::Hepg_E",&(_Hepg_E[0])); _ElectronNtuple->capture("HEPG::Hepg_Px",&(_Hepg_Px[0])); _ElectronNtuple->capture("HEPG::Hepg_Py",&(_Hepg_Py[0])); _ElectronNtuple->capture("HEPG::Hepg_Pz",&(_Hepg_Pz[0])); _ElectronNtuple->capture("HEPG::Hepg_Pt",&(_Hepg_Pt[0])); _ElectronNtuple->capture("HEPG::Hepg_Phi",&(_Hepg_Phi[0])); _ElectronNtuple->capture("HEPG::Hepg_Theta",&(_Hepg_Theta[0])); _ElectronNtuple->capture("HEPG::Hepg_Eta",&(_Hepg_Eta[0])); _ElectronNtuple->captureBlock("HEPG"); } _ElectronNtuple->storeCapturedData(); _ElectronNtuple->clearData(); return AppResult::OK; } AppModule* HighPtElectrons::clone(const char* cloneName) { return new HighPtElectrons(cloneName,"this module is a clone HighPtElectrons module"); } const char * HighPtElectrons::rcsId( ) const { return rcsid; } //this function gives you the id of the mother of the particle with index pIndex // double HighPtElectrons::getMother(const HEPG_StorableBank & _hepg, int pIndex) { // //index of the mother of pIndex; // int jmo1 = _hepg.jmo1hep(pIndex); // //id of pIndex; // int jtype = _hepg.idhep(pIndex-1); // if(jmo1==-999||jtype==-999||jtype==0||jmo1==0||jmo1<=2) return 0; // if(_hepg.idhep(jmo1-1)!=jtype) return jtype;// _hepg.idhep(jmo1-1); // else return getMother(_hepg,jmo1); // } const CdfTrack* HighPtElectrons::findMatchingTrackForCluster(EmCluster& emClust) { // Find a matching COT track for the given CdfEmObject. This is currently based on eta-phi matching : double etaPhiSepMin = 0.5; const CdfTrack* nearestTrack = 0; CdfTrackColl_ch trackSet; StorableObject::SelectByDescription description("GlobalCOT_Tracking"); CdfTrackColl::Error error = CdfTrackColl::find(trackSet,description); if (error == CdfTrackColl::OK) { for (CdfTrackColl::const_iterator trackIter = trackSet->contents().begin(); trackIter != trackSet->contents().end(); ++trackIter) if ((*trackIter)->pt() > _MatchingTrackPtCut.value()) { { double etaTrk = log(tan((atan((*trackIter)->lambda())+M_PI/2.0)/2.0)); double etaCls = emClust.emEtaEvent(); double phiSep = (*trackIter)->phi0()-emClust.emPhi(); while (phiSep > M_PI) phiSep= -2.0*M_PI+phiSep; while (phiSep < -M_PI) phiSep= 2.0*M_PI+phiSep; double etaPhiSep = sqrt(pow((etaTrk-etaCls),2)+pow(phiSep,2)); if (etaPhiSep < etaPhiSepMin) { nearestTrack = (*trackIter); etaPhiSepMin = etaPhiSep; } } } } return nearestTrack; }