CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Member Functions | Private Member Functions | Private Attributes
CSCEfficiency Class Reference

#include <CSCEfficiency.h>

Inheritance diagram for CSCEfficiency:
edm::EDFilter edm::ProducerBase edm::ProductRegistryHelper

Classes

struct  ChamberHistos
 
struct  StationHistos
 

Public Member Functions

 CSCEfficiency (const edm::ParameterSet &pset)
 Constructor. More...
 
virtual ~CSCEfficiency ()
 Destructor. More...
 
- Public Member Functions inherited from edm::EDFilter
 EDFilter ()
 
virtual ~EDFilter ()
 
- Public Member Functions inherited from edm::ProducerBase
 ProducerBase ()
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
boost::function< void(const
BranchDescription &)> 
registrationCallback () const
 used by the fwk to register list of products More...
 
virtual ~ProducerBase ()
 

Private Member Functions

bool applyTrigger (edm::Handle< edm::TriggerResults > &hltR, const edm::TriggerNames &triggerNames)
 
virtual void beginJob ()
 
void chamberCandidates (int station, int ring, float phi, std::vector< int > &coupleOfChambers)
 
bool checkLocal (double yLocal, double yBoundary, int station, int ring)
 
void chooseDirection (CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition)
 
bool efficienciesPerChamber (CSCDetId &id, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
 
virtual void endJob ()
 
double extrapolate1D (double initPosition, double initDirection, double parameterOfTheLine)
 
void fillDigiInfo (edm::Handle< CSCALCTDigiCollection > &alcts, edm::Handle< CSCCLCTDigiCollection > &clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > &correlatedlcts, edm::Handle< CSCWireDigiCollection > &wires, edm::Handle< CSCStripDigiCollection > &strips, edm::Handle< edm::PSimHitContainer > &simhits, edm::Handle< CSCRecHit2DCollection > &rechits, edm::Handle< CSCSegmentCollection > &segments, edm::ESHandle< CSCGeometry > &cscGeom)
 
void fillLCT_info (edm::Handle< CSCALCTDigiCollection > &alcts, edm::Handle< CSCCLCTDigiCollection > &clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > &correlatedlcts)
 
void fillRechitsSegments_info (edm::Handle< CSCRecHit2DCollection > &rechits, edm::Handle< CSCSegmentCollection > &segments, edm::ESHandle< CSCGeometry > &cscGeom)
 
void fillSimhit_info (edm::Handle< edm::PSimHitContainer > &simHits)
 
void fillStrips_info (edm::Handle< CSCStripDigiCollection > &strips)
 
void fillWG_info (edm::Handle< CSCWireDigiCollection > &wires, edm::ESHandle< CSCGeometry > &cscGeom)
 
virtual bool filter (edm::Event &event, const edm::EventSetup &eventSetup)
 
FreeTrajectoryState getFromCLHEP (const CLHEP::Hep3Vector &p3, const CLHEP::Hep3Vector &r3, int charge, const AlgebraicSymMatrix66 &cov, const MagneticField *field)
 
void getFromFTS (const FreeTrajectoryState &fts, CLHEP::Hep3Vector &p3, CLHEP::Hep3Vector &r3, int &charge, AlgebraicSymMatrix66 &cov)
 
bool inSensitiveLocalRegion (double xLocal, double yLocal, int station, int ring)
 
void linearExtrapolation (GlobalPoint initialPosition, GlobalVector initialDirection, float zSurface, std::vector< float > &posZY)
 
double lineParameter (double initZPosition, double destZPosition, double initZDirection)
 
TrajectoryStateOnSurface propagate (FreeTrajectoryState &ftsStart, const BoundPlane &bp)
 
const Propagatorpropagator (std::string propagatorName) const
 
bool recHitSegment_Efficiencies (CSCDetId &cscDetId, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
 
bool recSimHitEfficiency (CSCDetId &id, FreeTrajectoryState &ftsChamber)
 
void returnTypes (CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing)
 
void ringCandidates (int station, float absEta, std::map< std::string, bool > &chamberTypes)
 
bool stripWire_Efficiencies (CSCDetId &cscDetId, FreeTrajectoryState &ftsChamber)
 

Private Attributes

edm::InputTag alctDigiTag_
 
TH1F * ALCTPerEvent
 
bool allALCT [2][4][4][NumCh]
 
bool allCLCT [2][4][4][NumCh]
 
bool allCorrLCT [2][4][4][NumCh]
 
std::vector< std::pair
< LocalPoint, bool > > 
allRechits [2][4][4][NumCh][6]
 
std::vector< std::pair
< LocalPoint, LocalVector > > 
allSegments [2][4][4][NumCh]
 
std::vector< std::pair
< LocalPoint, int > > 
allSimhits [2][4][4][NumCh][6]
 
std::vector< std::pair< int,
float > > 
allStrips [2][4][4][NumCh][6]
 
std::vector< std::pair
< std::pair< int, float >, int > > 
allWG [2][4][4][NumCh][6]
 
bool alongZ
 
bool andOr
 
bool applyIPangleCuts
 
struct CSCEfficiency::ChamberHistos ChHist [2][4][3][LastCh-FirstCh+1]
 
edm::InputTag clctDigiTag_
 
TH1F * CLCTPerEvent
 
edm::InputTag corrlctDigiTag_
 
TH1F * DataFlow
 
double distanceFromDeadZone
 
bool emptyChambers [2][4][4][NumCh]
 
bool getAbsoluteEfficiency
 
edm::InputTag hlTriggerResults_
 
bool isBeamdata
 
bool isData
 
bool isIPdata
 
double local_DX_DZ_Max
 
double local_DY_DZ_Max
 
double local_DY_DZ_Min
 
bool magField
 
double maxNormChi2
 
double maxP
 
double minP
 
unsigned int minTrackHits
 
std::vector< std::string > myTriggers
 
int nEventsAnalyzed
 
bool passTheEvent
 
std::vector< int > pointToTriggers
 
bool printalot
 
unsigned int printout_NEvents
 
edm::InputTag rechitDigiTag_
 
TH1F * recHitsPerEvent
 
std::string rootFileName
 
edm::InputTag segmentDigiTag_
 
TH1F * segmentsPerEvent
 
edm::InputTag simHitTag
 
struct CSCEfficiency::StationHistos StHist [2][4]
 
edm::InputTag stripDigiTag_
 
TFile * theFile
 
MuonServiceProxytheService
 
edm::InputTag tracksTag
 
TH1F * TriggersFired
 
bool useDigis
 
bool useTrigger
 
edm::InputTag wireDigiTag_
 

Additional Inherited Members

- Public Types inherited from edm::EDFilter
typedef EDFilter ModuleType
 
typedef WorkerT< EDFilterWorkerType
 
- Public Types inherited from edm::ProducerBase
typedef
ProductRegistryHelper::TypeLabelList 
TypeLabelList
 
- Static Public Member Functions inherited from edm::EDFilter
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
- Protected Member Functions inherited from edm::EDFilter
CurrentProcessingContext const * currentContext () const
 
- Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)
 

Detailed Description

Efficiency calculations Stoyan Stoynev, Northwestern University

Definition at line 116 of file CSCEfficiency.h.

Constructor & Destructor Documentation

CSCEfficiency::CSCEfficiency ( const edm::ParameterSet pset)

Constructor.

Definition at line 1639 of file CSCEfficiency.cc.

References ExpressReco_HICollisions_FallBack::andOr, FirstCh, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), alignBH_cfg::minP, MuonServiceProxy_cff::MuonServiceProxy, NumCh, and interactiveExample::theFile.

1639  {
1640 
1641  // const float Xmin = -70;
1642  //const float Xmax = 70;
1643  //const int nXbins = int(4.*(Xmax - Xmin));
1644  const float Ymin = -165;
1645  const float Ymax = 165;
1646  const int nYbins = int((Ymax - Ymin)/2);
1647  const float Layer_min = -0.5;
1648  const float Layer_max = 9.5;
1649  const int nLayer_bins = int(Layer_max - Layer_min);
1650  //
1651 
1652  //---- Get the input parameters
1653  printout_NEvents = pset.getUntrackedParameter<unsigned int>("printout_NEvents",0);
1654  rootFileName = pset.getUntrackedParameter<string>("rootFileName","cscHists.root");
1655 
1656  isData = pset.getUntrackedParameter<bool>("runOnData",true);//
1657  isIPdata = pset.getUntrackedParameter<bool>("IPdata",false);//
1658  isBeamdata = pset.getUntrackedParameter<bool>("Beamdata",false);//
1659  getAbsoluteEfficiency = pset.getUntrackedParameter<bool>("getAbsoluteEfficiency",true);//
1660  useDigis = pset.getUntrackedParameter<bool>("useDigis", true);//
1661  distanceFromDeadZone = pset.getUntrackedParameter<double>("distanceFromDeadZone", 10.);//
1662  minP = pset.getUntrackedParameter<double>("minP",20.);//
1663  maxP = pset.getUntrackedParameter<double>("maxP",100.);//
1664  maxNormChi2 = pset.getUntrackedParameter<double>("maxNormChi2", 3.);//
1665  minTrackHits = pset.getUntrackedParameter<unsigned int>("minTrackHits",10);//
1666 
1667  applyIPangleCuts = pset.getUntrackedParameter<bool>("applyIPangleCuts", false);//
1668  local_DY_DZ_Max = pset.getUntrackedParameter<double>("local_DY_DZ_Max",-0.1);//
1669  local_DY_DZ_Min = pset.getUntrackedParameter<double>("local_DY_DZ_Min",-0.8);//
1670  local_DX_DZ_Max = pset.getUntrackedParameter<double>("local_DX_DZ_Max",0.2);//
1671 
1672  alctDigiTag_ = pset.getParameter<edm::InputTag>("alctDigiTag") ;
1673  clctDigiTag_ = pset.getParameter<edm::InputTag>("clctDigiTag") ;
1674  corrlctDigiTag_ = pset.getParameter<edm::InputTag>("corrlctDigiTag") ;
1675  stripDigiTag_ = pset.getParameter<edm::InputTag>("stripDigiTag") ;
1676  wireDigiTag_ = pset.getParameter<edm::InputTag>("wireDigiTag") ;
1677  rechitDigiTag_ = pset.getParameter<edm::InputTag>("rechitDigiTag") ;
1678  segmentDigiTag_ = pset.getParameter<edm::InputTag>("segmentDigiTag") ;
1679  simHitTag = pset.getParameter<edm::InputTag>("simHitTag");
1680  tracksTag = pset.getParameter< edm::InputTag >("tracksTag");
1681 
1682  ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
1683  // maybe use the service for getting magnetic field, propagators, etc. ...
1684  theService = new MuonServiceProxy(serviceParameters);
1685 
1686  // Trigger
1687  useTrigger = pset.getUntrackedParameter<bool>("useTrigger", false);
1688  hlTriggerResults_ = pset.getParameter<edm::InputTag> ("HLTriggerResults");
1689  myTriggers = pset.getParameter<std::vector <std::string> >("myTriggers");
1690  andOr = pset.getUntrackedParameter<bool>("andOr");
1691  pointToTriggers.clear();
1692 
1693 
1694  //---- set counter to zero
1695  nEventsAnalyzed = 0;
1696  //---- set presence of magnetic field
1697  magField = true;
1698  //
1699  std::string Path = "AllChambers/";
1700  std::string FullName;
1701  //---- File with output histograms
1702  theFile = new TFile(rootFileName.c_str(), "RECREATE");
1703  theFile->cd();
1704  //---- Book histograms for the analysis
1705  char SpecName[50];
1706 
1707  sprintf(SpecName,"DataFlow");
1708  DataFlow =
1709  new TH1F(SpecName,"Data flow;condition number;entries",40,-0.5,39.5);
1710  //
1711  sprintf(SpecName,"TriggersFired");
1712  TriggersFired =
1713  new TH1F(SpecName,"Triggers fired;trigger number;entries",140,-0.5,139.5);
1714  //
1715  int Chan = 50;
1716  float minChan = -0.5;
1717  float maxChan = 49.5;
1718  //
1719  sprintf(SpecName,"ALCTPerEvent");
1720  ALCTPerEvent = new TH1F(SpecName,"ALCTs per event;N digis;entries",Chan,minChan,maxChan);
1721  //
1722  sprintf(SpecName,"CLCTPerEvent");
1723  CLCTPerEvent = new TH1F(SpecName,"CLCTs per event;N digis;entries",Chan,minChan,maxChan);
1724  //
1725  sprintf(SpecName,"recHitsPerEvent");
1726  recHitsPerEvent = new TH1F(SpecName,"RecHits per event;N digis;entries",150,-0.5,149.5);
1727  //
1728  sprintf(SpecName,"segmentsPerEvent");
1729  segmentsPerEvent = new TH1F(SpecName,"segments per event;N digis;entries",Chan,minChan,maxChan);
1730  //
1731  //---- Book groups of histograms (for any chamber)
1732 
1733  map<std::string,bool>::iterator iter;
1734  for(int ec = 0;ec<2;++ec){
1735  for(int st = 0;st<4;++st){
1736  theFile->cd();
1737  sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
1738  theFile->mkdir(SpecName);
1739  theFile->cd(SpecName);
1740 
1741  //
1742  sprintf(SpecName,"segmentChi2_ndf_St%d",st+1);
1743  StHist[ec][st].segmentChi2_ndf =
1744  new TH1F(SpecName,"Chi2/ndf of a segment;chi2/ndf;entries",100,0.,20.);
1745  //
1746  sprintf(SpecName,"hitsInSegment_St%d",st+1);
1747  StHist[ec][st].hitsInSegment =
1748  new TH1F(SpecName,"Number of hits in a segment;nHits;entries",7,-0.5,6.5);
1749  //
1750  Chan = 170;
1751  minChan = 0.85;
1752  maxChan = 2.55;
1753  //
1754  sprintf(SpecName,"AllSegments_eta_St%d",st+1);
1755  StHist[ec][st].AllSegments_eta =
1756  new TH1F(SpecName,"All segments in eta;eta;entries",Chan,minChan,maxChan);
1757  //
1758  sprintf(SpecName,"EfficientSegments_eta_St%d",st+1);
1759  StHist[ec][st].EfficientSegments_eta =
1760  new TH1F(SpecName,"Efficient segments in eta;eta;entries",Chan,minChan,maxChan);
1761  //
1762  sprintf(SpecName,"ResidualSegments_St%d",st+1);
1763  StHist[ec][st].ResidualSegments =
1764  new TH1F(SpecName,"Residual (segments);residual,cm;entries",75,0.,15.);
1765  //
1766  Chan = 200;
1767  minChan = -800.;
1768  maxChan = 800.;
1769  int Chan2 = 200;
1770  float minChan2 = -800.;
1771  float maxChan2 = 800.;
1772 
1773  sprintf(SpecName,"EfficientSegments_XY_St%d",st+1);
1774  StHist[ec][st].EfficientSegments_XY = new TH2F(SpecName,"Efficient segments in XY;X;Y",
1775  Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
1776  sprintf(SpecName,"InefficientSegments_XY_St%d",st+1);
1777  StHist[ec][st].InefficientSegments_XY = new TH2F(SpecName,"Inefficient segments in XY;X;Y",
1778  Chan,minChan,maxChan,Chan2,minChan2,maxChan2);
1779  //
1780  Chan = 80;
1781  minChan = 0;
1782  maxChan = 3.2;
1783  sprintf(SpecName,"EfficientALCT_momTheta_St%d",st+1);
1784  StHist[ec][st].EfficientALCT_momTheta = new TH1F(SpecName,"Efficient ALCT in theta (momentum);theta, rad;entries",
1785  Chan,minChan,maxChan);
1786  //
1787  sprintf(SpecName,"InefficientALCT_momTheta_St%d",st+1);
1788  StHist[ec][st].InefficientALCT_momTheta = new TH1F(SpecName,"Inefficient ALCT in theta (momentum);theta, rad;entries",
1789  Chan,minChan,maxChan);
1790  //
1791  Chan = 160;
1792  minChan = -3.2;
1793  maxChan = 3.2;
1794  sprintf(SpecName,"EfficientCLCT_momPhi_St%d",st+1);
1795  StHist[ec][st].EfficientCLCT_momPhi = new TH1F(SpecName,"Efficient CLCT in phi (momentum);phi, rad;entries",
1796  Chan,minChan,maxChan);
1797  //
1798  sprintf(SpecName,"InefficientCLCT_momPhi_St%d",st+1);
1799  StHist[ec][st].InefficientCLCT_momPhi = new TH1F(SpecName,"Inefficient CLCT in phi (momentum);phi, rad;entries",
1800  Chan,minChan,maxChan);
1801  //
1802  theFile->cd();
1803  for(int rg = 0;rg<3;++rg){
1804  if(0!=st && rg>1){
1805  continue;
1806  }
1807  else if(1==rg && 3==st){
1808  continue;
1809  }
1810  for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
1811  if(0!=st && 0==rg && iChamber >18){
1812  continue;
1813  }
1814  theFile->cd();
1815  sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
1816  theFile->mkdir(SpecName);
1817  theFile->cd(SpecName);
1818  //
1819 
1820  sprintf(SpecName,"EfficientRechits_inSegment_Ch%d",iChamber);
1821  ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment =
1822  new TH1F(SpecName,"Existing RecHit given a segment;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1823  //
1824  sprintf(SpecName,"InefficientSingleHits_Ch%d",iChamber);
1825  ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits =
1826  new TH1F(SpecName,"Single RecHits not in the segment;layers (1-6);entries ",nLayer_bins,Layer_min,Layer_max);
1827  //
1828  sprintf(SpecName,"AllSingleHits_Ch%d",iChamber);
1829  ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits =
1830  new TH1F(SpecName,"Single RecHits given a segment; layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1831  //
1832  sprintf(SpecName,"digiAppearanceCount_Ch%d",iChamber);
1833  ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount =
1834  new TH1F(SpecName,"Digi appearance (no-yes): segment(0,1), ALCT(2,3), CLCT(4,5), CorrLCT(6,7); digi type;entries",
1835  8,-0.5,7.5);
1836  //
1837  Chan = 100;
1838  minChan = -1.1;
1839  maxChan = 0.9;
1840  sprintf(SpecName,"EfficientALCT_dydz_Ch%d",iChamber);
1841  ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz =
1842  new TH1F(SpecName,"Efficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
1843  Chan, minChan, maxChan);
1844  //
1845  sprintf(SpecName,"InefficientALCT_dydz_Ch%d",iChamber);
1846  ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz =
1847  new TH1F(SpecName,"Inefficient ALCT; local dy/dz (ME 3 and 4 flipped);entries",
1848  Chan, minChan, maxChan);
1849  //
1850  Chan = 100;
1851  minChan = -1.;
1852  maxChan = 1.0;
1853  sprintf(SpecName,"EfficientCLCT_dxdz_Ch%d",iChamber);
1854  ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz =
1855  new TH1F(SpecName,"Efficient CLCT; local dxdz;entries",
1856  Chan, minChan, maxChan);
1857  //
1858  sprintf(SpecName,"InefficientCLCT_dxdz_Ch%d",iChamber);
1859  ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz =
1860  new TH1F(SpecName,"Inefficient CLCT; local dxdz;entries",
1861  Chan, minChan, maxChan);
1862  //
1863  sprintf(SpecName,"EfficientRechits_good_Ch%d",iChamber);
1864  ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good =
1865  new TH1F(SpecName,"Existing RecHit - sensitive area only;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1866  //
1867  sprintf(SpecName,"EfficientStrips_Ch%d",iChamber);
1868  ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips =
1869  new TH1F(SpecName,"Existing strip;layer (1-6); entries",nLayer_bins,Layer_min,Layer_max);
1870  //
1871  sprintf(SpecName,"EfficientWireGroups_Ch%d",iChamber);
1872  ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups =
1873  new TH1F(SpecName,"Existing WireGroups;layer (1-6); entries ",nLayer_bins,Layer_min,Layer_max);
1874  //
1875  sprintf(SpecName,"StripWiresCorrelations_Ch%d",iChamber);
1876  ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations =
1877  new TH1F(SpecName,"StripWire correlations;; entries ",5,0.5,5.5);
1878  //
1879  Chan = 80;
1880  minChan = 0;
1881  maxChan = 3.2;
1882  sprintf(SpecName,"NoWires_momTheta_Ch%d",iChamber);
1883  ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta =
1884  new TH1F(SpecName,"No wires (all strips present) - in theta (momentum);theta, rad;entries",
1885  Chan,minChan,maxChan);
1886  //
1887  Chan = 160;
1888  minChan = -3.2;
1889  maxChan = 3.2;
1890  sprintf(SpecName,"NoStrips_momPhi_Ch%d",iChamber);
1891  ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi =
1892  new TH1F(SpecName,"No strips (all wires present) - in phi (momentum);phi, rad;entries",
1893  Chan,minChan,maxChan);
1894  //
1895  for(int iLayer=0; iLayer<6;iLayer++){
1896  sprintf(SpecName,"Y_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
1897  ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment.push_back
1898  (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
1899  nYbins,Ymin, Ymax));
1900  //
1901  sprintf(SpecName,"Y_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
1902  ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment.push_back
1903  (new TH1F(SpecName,"Efficient (extrapolated from the segment) RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
1904  nYbins,Ymin, Ymax));
1905  //
1906  Chan = 200;
1907  minChan = -0.2;
1908  maxChan = 0.2;
1909  sprintf(SpecName,"Phi_InefficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
1910  ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment.push_back
1911  (new TH1F(SpecName,"Missing RecHit/layer in a segment (local system, whole chamber);Phi, rad; entries",
1912  Chan, minChan, maxChan));
1913  //
1914  sprintf(SpecName,"Phi_EfficientRecHits_inSegment_Ch%d_L%d",iChamber,iLayer);
1915  ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment.push_back
1916  (new TH1F(SpecName,"Efficient (extrapolated from the segment) in a segment (local system, whole chamber);Phi, rad; entries",
1917  Chan, minChan, maxChan));
1918 
1919  }
1920  //
1921  sprintf(SpecName,"Sim_Rechits_Ch%d",iChamber);
1922  ChHist[ec][st][rg][iChamber-FirstCh].SimRechits =
1923  new TH1F(SpecName,"Existing RecHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1924  //
1925  sprintf(SpecName,"Sim_Simhits_Ch%d",iChamber);
1926  ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits =
1927  new TH1F(SpecName,"Existing SimHit (Sim);layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1928  //
1929  /*
1930  sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber);
1931  ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each =
1932  new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1933  //
1934  sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber);
1935  ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each =
1936  new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1937  */
1938  theFile->cd();
1939  }
1940  }
1941  }
1942  }
1943 }
edm::InputTag corrlctDigiTag_
T getParameter(std::string const &) const
struct CSCEfficiency::StationHistos StHist[2][4]
T getUntrackedParameter(std::string const &, T const &) const
std::string rootFileName
edm::InputTag segmentDigiTag_
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
TH1F * CLCTPerEvent
std::vector< TH1F * > Y_InefficientRecHits_inSegment
double local_DX_DZ_Max
bool getAbsoluteEfficiency
std::vector< std::string > myTriggers
std::vector< TH1F * > Phi_InefficientRecHits_inSegment
double local_DY_DZ_Max
std::vector< TH1F * > Y_EfficientRecHits_inSegment
MuonServiceProxy * theService
Definition: Path.h:29
edm::InputTag rechitDigiTag_
edm::InputTag hlTriggerResults_
double distanceFromDeadZone
std::vector< TH1F * > Phi_EfficientRecHits_inSegment
double maxNormChi2
double local_DY_DZ_Min
#define FirstCh
edm::InputTag tracksTag
unsigned int printout_NEvents
TH1F * TriggersFired
edm::InputTag clctDigiTag_
TH1F * segmentsPerEvent
TH1F * recHitsPerEvent
edm::InputTag simHitTag
std::vector< int > pointToTriggers
TH1F * ALCTPerEvent
edm::InputTag stripDigiTag_
#define NumCh
edm::InputTag alctDigiTag_
edm::InputTag wireDigiTag_
unsigned int minTrackHits
CSCEfficiency::~CSCEfficiency ( )
virtual

Destructor.

Definition at line 1946 of file CSCEfficiency.cc.

References gather_cfg::cout, FirstCh, NumCh, and interactiveExample::theFile.

1946  {
1947  if (theService) delete theService;
1948  // Write the histos to a file
1949  theFile->cd();
1950  //
1951  char SpecName[20];
1952  std::vector<float> bins, Efficiency, EffError;
1953  std::vector<float> eff(2);
1954 
1955  //---- loop over chambers
1956  std::map <std::string, bool> chamberTypes;
1957  chamberTypes["ME11"] = false;
1958  chamberTypes["ME12"] = false;
1959  chamberTypes["ME13"] = false;
1960  chamberTypes["ME21"] = false;
1961  chamberTypes["ME22"] = false;
1962  chamberTypes["ME31"] = false;
1963  chamberTypes["ME32"] = false;
1964  chamberTypes["ME41"] = false;
1965 
1966  map<std::string,bool>::iterator iter;
1967  std::cout<<" Writing proper histogram structure (patience)..."<<std::endl;
1968  for(int ec = 0;ec<2;++ec){
1969  for(int st = 0;st<4;++st){
1970  sprintf(SpecName,"Stations__E%d_S%d",ec+1, st+1);
1971  theFile->cd(SpecName);
1972  StHist[ec][st].segmentChi2_ndf->Write();
1973  StHist[ec][st].hitsInSegment->Write();
1974  StHist[ec][st].AllSegments_eta->Write();
1975  StHist[ec][st].EfficientSegments_eta->Write();
1976  StHist[ec][st].ResidualSegments->Write();
1977  StHist[ec][st].EfficientSegments_XY->Write();
1978  StHist[ec][st].InefficientSegments_XY->Write();
1979  StHist[ec][st].EfficientALCT_momTheta->Write();
1980  StHist[ec][st].InefficientALCT_momTheta->Write();
1981  StHist[ec][st].EfficientCLCT_momPhi->Write();
1982  StHist[ec][st].InefficientCLCT_momPhi->Write();
1983  for(int rg = 0;rg<3;++rg){
1984  if(0!=st && rg>1){
1985  continue;
1986  }
1987  else if(1==rg && 3==st){
1988  continue;
1989  }
1990  for(int iChamber=FirstCh;iChamber<FirstCh+NumCh;iChamber++){
1991  if(0!=st && 0==rg && iChamber >18){
1992  continue;
1993  }
1994  sprintf(SpecName,"Chambers__E%d_S%d_R%d_Chamber_%d",ec+1, st+1, rg+1,iChamber);
1995  theFile->cd(SpecName);
1996 
1997  ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_inSegment->Write();
1998  ChHist[ec][st][rg][iChamber-FirstCh].AllSingleHits->Write();
1999  ChHist[ec][st][rg][iChamber-FirstCh].digiAppearanceCount->Write();
2000  ChHist[ec][st][rg][iChamber-FirstCh].EfficientALCT_dydz->Write();
2001  ChHist[ec][st][rg][iChamber-FirstCh].InefficientALCT_dydz->Write();
2002  ChHist[ec][st][rg][iChamber-FirstCh].EfficientCLCT_dxdz->Write();
2003  ChHist[ec][st][rg][iChamber-FirstCh].InefficientCLCT_dxdz->Write();
2004  ChHist[ec][st][rg][iChamber-FirstCh].InefficientSingleHits->Write();
2005  ChHist[ec][st][rg][iChamber-FirstCh].EfficientRechits_good->Write();
2006  ChHist[ec][st][rg][iChamber-FirstCh].EfficientStrips->Write();
2007  ChHist[ec][st][rg][iChamber-FirstCh].StripWiresCorrelations->Write();
2008  ChHist[ec][st][rg][iChamber-FirstCh].NoWires_momTheta->Write();
2009  ChHist[ec][st][rg][iChamber-FirstCh].NoStrips_momPhi->Write();
2010  ChHist[ec][st][rg][iChamber-FirstCh].EfficientWireGroups->Write();
2011  for(unsigned int iLayer = 0; iLayer< 6; iLayer++){
2012  ChHist[ec][st][rg][iChamber-FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write();
2013  ChHist[ec][st][rg][iChamber-FirstCh].Y_EfficientRecHits_inSegment[iLayer]->Write();
2014  ChHist[ec][st][rg][iChamber-FirstCh].Phi_InefficientRecHits_inSegment[iLayer]->Write();
2015  ChHist[ec][st][rg][iChamber-FirstCh].Phi_EfficientRecHits_inSegment[iLayer]->Write();
2016  }
2017  ChHist[ec][st][rg][iChamber-FirstCh].SimRechits->Write();
2018  ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits->Write();
2019  /*
2020  ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each->Write();
2021  ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each->Write();
2022  */
2023  //
2024  theFile->cd(SpecName);
2025  theFile->cd();
2026  }
2027  }
2028  }
2029  }
2030  //
2031  sprintf(SpecName,"AllChambers");
2032  theFile->mkdir(SpecName);
2033  theFile->cd(SpecName);
2034  DataFlow->Write();
2035  TriggersFired->Write();
2036  ALCTPerEvent->Write();
2037  CLCTPerEvent->Write();
2038  recHitsPerEvent->Write();
2039  segmentsPerEvent->Write();
2040  //
2041  theFile->cd(SpecName);
2042  //---- Close the file
2043  theFile->Close();
2044 }
struct CSCEfficiency::StationHistos StHist[2][4]
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
TH1F * CLCTPerEvent
std::vector< TH1F * > Y_InefficientRecHits_inSegment
std::vector< TH1F * > Phi_InefficientRecHits_inSegment
std::vector< TH1F * > Y_EfficientRecHits_inSegment
MuonServiceProxy * theService
std::vector< TH1F * > Phi_EfficientRecHits_inSegment
#define FirstCh
TH1F * TriggersFired
TH1F * segmentsPerEvent
TH1F * recHitsPerEvent
tuple cout
Definition: gather_cfg.py:41
TH1F * ALCTPerEvent
#define NumCh

Member Function Documentation

bool CSCEfficiency::applyTrigger ( edm::Handle< edm::TriggerResults > &  hltR,
const edm::TriggerNames triggerNames 
)
private

Definition at line 1557 of file CSCEfficiency.cc.

References ExpressReco_HICollisions_FallBack::andOr, gather_cfg::cout, edm::HandleBase::isValid(), and edm::TriggerNames::triggerNames().

1558  {
1559  bool triggerPassed = true;
1560  std::vector<std::string> hlNames=triggerNames.triggerNames();
1561  pointToTriggers.clear();
1562  for(size_t imyT = 0;imyT<myTriggers.size();++imyT){
1563  for (size_t iT=0; iT<hlNames.size(); ++iT) {
1564  //std::cout<<" iT = "<<iT<<" hlNames[iT] = "<<hlNames[iT]<<
1565  //" : wasrun = "<<hltR->wasrun(iT)<<" accept = "<<
1566  // hltR->accept(iT)<<" !error = "<<
1567  // !hltR->error(iT)<<std::endl;
1568  if(!imyT){
1569  if(hltR->wasrun(iT) &&
1570  hltR->accept(iT) &&
1571  !hltR->error(iT) ){
1572  TriggersFired->Fill(iT);
1573  }
1574  }
1575  if(hlNames[iT]==myTriggers[imyT]){
1576  pointToTriggers.push_back(iT);
1577  if(imyT){
1578  break;
1579  }
1580  }
1581  }
1582  }
1583  if(pointToTriggers.size()!=myTriggers.size()){
1584  pointToTriggers.clear();
1585  if(printalot){
1586  std::cout<<" Not all trigger names found - all trigger specifications will be ignored. Check your cfg file!"<<std::endl;
1587  }
1588  }
1589  else{
1590  if(pointToTriggers.size()){
1591  if(printalot){
1592  std::cout<<"The following triggers will be required in the event: "<<std::endl;
1593  for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
1594  std::cout<<" "<<hlNames[pointToTriggers[imyT]];
1595  }
1596  std::cout<<std::endl;
1597  std::cout<<" in condition (AND/OR) : "<<!andOr<<"/"<<andOr<<std::endl;
1598  }
1599  }
1600  }
1601 
1602  if (hltR.isValid()) {
1603  if(!pointToTriggers.size()){
1604  if(printalot){
1605  std::cout<<" No triggers specified in the configuration or all ignored - no trigger information will be considered"<<std::endl;
1606  }
1607  }
1608  for(size_t imyT =0; imyT <pointToTriggers.size();++imyT){
1609  if(hltR->wasrun(pointToTriggers[imyT]) &&
1610  hltR->accept(pointToTriggers[imyT]) &&
1611  !hltR->error(pointToTriggers[imyT]) ){
1612  triggerPassed = true;
1613  if(andOr){
1614  break;
1615  }
1616  }
1617  else{
1618  triggerPassed = false;
1619  if(!andOr){
1620  triggerPassed = false;
1621  break;
1622  }
1623  }
1624  }
1625  }
1626  else{
1627  if(printalot){
1628  std::cout<<" TriggerResults handle returns invalid state?! No trigger information will be considered"<<std::endl;
1629  }
1630  }
1631  if(printalot){
1632  std::cout<<" Trigger passed: "<<triggerPassed<<std::endl;
1633  }
1634  return triggerPassed;
1635 }
Strings const & triggerNames() const
Definition: TriggerNames.cc:24
std::vector< std::string > myTriggers
bool isValid() const
Definition: HandleBase.h:76
TH1F * TriggersFired
std::vector< int > pointToTriggers
tuple cout
Definition: gather_cfg.py:41
void CSCEfficiency::beginJob ( void  )
privatevirtual

Reimplemented from edm::EDFilter.

Definition at line 2048 of file CSCEfficiency.cc.

2049 {
2050 }
void CSCEfficiency::chamberCandidates ( int  station,
int  ring,
float  phi,
std::vector< int > &  coupleOfChambers 
)
private

Definition at line 1017 of file CSCEfficiency.cc.

References gather_cfg::cout, and M_PI.

1017  {
1018  coupleOfChambers.clear();
1019  // -pi< phi<+pi
1020  float phi_zero = 0.;// check! the phi at the "edge" of Ch 1
1021  float phi_const = 2.*M_PI/36.;
1022  int last_chamber = 36;
1023  int first_chamber = 1;
1024  if(1 != station && 1==ring){ // 18 chambers in the ring
1025  phi_const*=2;
1026  last_chamber /= 2;
1027  }
1028  if(phi<0.){
1029  if (printalot) std::cout<<" info: negative phi = "<<phi<<std::endl;
1030  phi += 2*M_PI;
1031  }
1032  float chamber_float = (phi - phi_zero)/phi_const;
1033  int chamber_int = int(chamber_float);
1034  if (chamber_float - float(chamber_int) -0.5 <0.){
1035  if(0!=chamber_int ){
1036  coupleOfChambers.push_back(chamber_int);
1037  }
1038  else{
1039  coupleOfChambers.push_back(last_chamber);
1040  }
1041  coupleOfChambers.push_back(chamber_int+1);
1042 
1043  }
1044  else{
1045  coupleOfChambers.push_back(chamber_int+1);
1046  if(last_chamber!=chamber_int+1){
1047  coupleOfChambers.push_back(chamber_int+2);
1048  }
1049  else{
1050  coupleOfChambers.push_back(first_chamber);
1051  }
1052  }
1053  if (printalot) std::cout<<" phi = "<<phi<<" phi_zero = "<<phi_zero<<" phi_const = "<<phi_const<<
1054  " candidate chambers: first ch = "<<coupleOfChambers[0]<<" second ch = "<<coupleOfChambers[1]<<std::endl;
1055 }
#define M_PI
Definition: BFit3D.cc:3
tuple cout
Definition: gather_cfg.py:41
Definition: DDAxes.h:10
bool CSCEfficiency::checkLocal ( double  yLocal,
double  yBoundary,
int  station,
int  ring 
)
private

Definition at line 605 of file CSCEfficiency.cc.

605  {
606 //---- check if it is in a good local region (sensitive area - geometrical and HV boundaries excluded)
607  bool pass = false;
608  std::vector <float> deadZoneCenter(6);
609  const float deadZoneHalf = 0.32*7/2;// wire spacing * (wires missing + 1)/2
610  float cutZone = deadZoneHalf + distanceFromDeadZone;//cm
611  //---- hardcoded... not good
612  if(station>1 && station<5){
613  if(2==ring){
614  deadZoneCenter[0]= -162.48 ;
615  deadZoneCenter[1] = -81.8744;
616  deadZoneCenter[2] = -21.18165;
617  deadZoneCenter[3] = 39.51105;
618  deadZoneCenter[4] = 100.2939;
619  deadZoneCenter[5] = 160.58;
620 
621  if(yLocal >yBoundary &&
622  ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
623  (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
624  (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone) ||
625  (yLocal> deadZoneCenter[3] + cutZone && yLocal< deadZoneCenter[4] - cutZone) ||
626  (yLocal> deadZoneCenter[4] + cutZone && yLocal< deadZoneCenter[5] - cutZone))){
627  pass = true;
628  }
629  }
630  else if(1==ring){
631  if(2==station){
632  deadZoneCenter[0]= -95.94 ;
633  deadZoneCenter[1] = -27.47;
634  deadZoneCenter[2] = 33.67;
635  deadZoneCenter[3] = 93.72;
636  }
637  else if(3==station){
638  deadZoneCenter[0]= -85.97 ;
639  deadZoneCenter[1] = -36.21;
640  deadZoneCenter[2] = 23.68;
641  deadZoneCenter[3] = 84.04;
642  }
643  else if(4==station){
644  deadZoneCenter[0]= -75.82;
645  deadZoneCenter[1] = -26.14;
646  deadZoneCenter[2] = 23.85;
647  deadZoneCenter[3] = 73.91;
648  }
649  if(yLocal >yBoundary &&
650  ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
651  (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
652  (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
653  pass = true;
654  }
655  }
656  }
657  else if(1==station){
658  if(3==ring){
659  deadZoneCenter[0]= -83.155 ;
660  deadZoneCenter[1] = -22.7401;
661  deadZoneCenter[2] = 27.86665;
662  deadZoneCenter[3] = 81.005;
663  if(yLocal > yBoundary &&
664  ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
665  (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
666  (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
667  pass = true;
668  }
669  }
670  else if(2==ring){
671  deadZoneCenter[0]= -86.285 ;
672  deadZoneCenter[1] = -32.88305;
673  deadZoneCenter[2] = 32.867423;
674  deadZoneCenter[3] = 88.205;
675  if(yLocal > (yBoundary) &&
676  ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) ||
677  (yLocal> deadZoneCenter[1] + cutZone && yLocal< deadZoneCenter[2] - cutZone) ||
678  (yLocal> deadZoneCenter[2] + cutZone && yLocal< deadZoneCenter[3] - cutZone))){
679  pass = true;
680  }
681  }
682  else{
683  deadZoneCenter[0]= -81.0;
684  deadZoneCenter[1] = 81.0;
685  if(yLocal > (yBoundary) &&
686  ((yLocal> deadZoneCenter[0] + cutZone && yLocal< deadZoneCenter[1] - cutZone) )){
687  pass = true;
688  }
689  }
690  }
691  return pass;
692 }
double distanceFromDeadZone
void CSCEfficiency::chooseDirection ( CLHEP::Hep3Vector &  innerPosition,
CLHEP::Hep3Vector &  outerPosition 
)
private

Definition at line 1499 of file CSCEfficiency.cc.

1499  {
1500 
1501  //---- Be careful with trigger conditions too
1502  if(!isIPdata){
1503  float dy = outerPosition.y() - innerPosition.y();
1504  float dz = outerPosition.z() - innerPosition.z();
1505  if(isBeamdata){
1506  if(dz>0){
1507  alongZ = true;
1508  }
1509  else{
1510  alongZ = false;
1511  }
1512  }
1513  else{//cosmics
1514  if(dy/dz>0){
1515  alongZ = false;
1516  }
1517  else{
1518  alongZ = true;
1519  }
1520  }
1521  }
1522 }
bool CSCEfficiency::efficienciesPerChamber ( CSCDetId id,
const CSCChamber cscChamber,
FreeTrajectoryState ftsChamber 
)
private

Definition at line 1058 of file CSCEfficiency.cc.

References gather_cfg::cout, FreeTrajectoryState::momentum(), dbtoconf::out, PV3DBase< T, PVType, FrameType >::phi(), PV3DBase< T, PVType, FrameType >::theta(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

1058  {
1059  int ec, st, rg, ch, secondRing;
1060  returnTypes(id, ec, st, rg, ch, secondRing);
1061 
1062  LocalVector localDir = cscChamber->toLocal(ftsChamber.momentum());
1063  if(printalot){
1064  std::cout<<" global dir = "<<ftsChamber.momentum()<<std::endl;
1065  std::cout<<" local dir = "<<localDir<<std::endl;
1066  std::cout<<" local theta = "<<localDir.theta()<<std::endl;
1067  }
1068  float dxdz = localDir.x()/localDir.z();
1069  float dydz = localDir.y()/localDir.z();
1070  if(2==st || 3==st){
1071  if(printalot){
1072  std::cout<<"st 3 or 4 ... flip dy/dz"<<std::endl;
1073  }
1074  dydz = - dydz;
1075  }
1076  if(printalot){
1077  std::cout<<"dy/dz = "<<dydz<<std::endl;
1078  }
1079  // Apply angle cut
1080  bool out = true;
1081  if(applyIPangleCuts){
1082  if(dydz>local_DY_DZ_Max || dydz<local_DY_DZ_Min || fabs(dxdz)>local_DX_DZ_Max){
1083  out = false;
1084  }
1085  }
1086 
1087  // Segments
1088  bool firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
1089  bool secondCondition = false;
1090  //---- ME1 is special as usual - ME1a and ME1b are actually one chamber
1091  if(secondRing>-1){
1092  secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
1093  }
1094  if(firstCondition || secondCondition){
1095  if(out){
1096  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(1);
1097  }
1098  }
1099  else{
1100  if(out){
1101  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(0);
1102  }
1103  }
1104 
1105  bool missingALCT = false;
1106  bool missingCLCT = false;
1107  if(useDigis){
1108  // ALCTs
1109  firstCondition = allALCT[ec][st][rg][ch];
1110  secondCondition = false;
1111  if(secondRing>-1){
1112  secondCondition = allALCT[ec][st][secondRing][ch];
1113  }
1114  if(firstCondition || secondCondition){
1115  if(out){
1116  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(3);
1117  }
1118  // always apply partial angle cuts for this kind of histos
1119  if(fabs(dxdz)<local_DX_DZ_Max){
1120  StHist[ec][st].EfficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
1121  ChHist[ec][st][rg][ch].EfficientALCT_dydz->Fill(dydz);
1122  }
1123  }
1124  else{
1125  missingALCT = true;
1126  if(out){
1127  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(2);
1128  }
1129  if(fabs(dxdz)<local_DX_DZ_Max){
1130  StHist[ec][st].InefficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
1131  ChHist[ec][st][rg][ch].InefficientALCT_dydz->Fill(dydz);
1132  }
1133  if(printalot){
1134  std::cout<<" missing ALCT (dy/dz = "<<dydz<<")";
1135  printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
1136  }
1137  }
1138 
1139  // CLCTs
1140  firstCondition = allCLCT[ec][st][rg][ch];
1141  secondCondition = false;
1142  if(secondRing>-1){
1143  secondCondition = allCLCT[ec][st][secondRing][ch];
1144  }
1145  if(firstCondition || secondCondition){
1146  if(out){
1147  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(5);
1148  }
1149  if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){
1150  StHist[ec][st].EfficientCLCT_momPhi->Fill(ftsChamber.momentum().phi() );// - phi chamber...
1151  ChHist[ec][st][rg][ch].EfficientCLCT_dxdz->Fill(dxdz);
1152  }
1153  }
1154  else{
1155  missingCLCT = true;
1156  if(out){
1157  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(4);
1158  }
1159  if(dydz<local_DY_DZ_Max && dydz>local_DY_DZ_Min){
1160  StHist[ec][st].InefficientCLCT_momPhi->Fill(ftsChamber.momentum().phi());// - phi chamber...
1161  ChHist[ec][st][rg][ch].InefficientCLCT_dxdz->Fill(dxdz);
1162  }
1163  if(printalot){
1164  std::cout<<" missing CLCT (dx/dz = "<<dxdz<<")";
1165  printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",ec+1,st+1,rg+1,ch+1);
1166  }
1167  }
1168  if(out){
1169  // CorrLCTs
1170  firstCondition = allCorrLCT[ec][st][rg][ch];
1171  secondCondition = false;
1172  if(secondRing>-1){
1173  secondCondition = allCorrLCT[ec][st][secondRing][ch];
1174  }
1175  if(firstCondition || secondCondition){
1176  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(7);
1177  }
1178  else{
1179  ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(6);
1180  }
1181  }
1182  }
1183  return out;
1184 }
struct CSCEfficiency::StationHistos StHist[2][4]
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
T y() const
Definition: PV3DBase.h:57
bool allCLCT[2][4][4][NumCh]
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:64
double local_DX_DZ_Max
void returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing)
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
double local_DY_DZ_Max
T z() const
Definition: PV3DBase.h:58
bool allALCT[2][4][4][NumCh]
GlobalVector momentum() const
double local_DY_DZ_Min
tuple out
Definition: dbtoconf.py:99
std::vector< std::pair< LocalPoint, LocalVector > > allSegments[2][4][4][NumCh]
tuple cout
Definition: gather_cfg.py:41
T x() const
Definition: PV3DBase.h:56
bool allCorrLCT[2][4][4][NumCh]
void CSCEfficiency::endJob ( void  )
privatevirtual

Reimplemented from edm::EDFilter.

Definition at line 2054 of file CSCEfficiency.cc.

2054  {
2055 }
double CSCEfficiency::extrapolate1D ( double  initPosition,
double  initDirection,
double  parameterOfTheLine 
)
private

Definition at line 1489 of file CSCEfficiency.cc.

1489  {
1490  double extrapolatedPosition = initPosition + initDirection*parameterOfTheLine;
1491  return extrapolatedPosition;
1492 }
void CSCEfficiency::fillDigiInfo ( edm::Handle< CSCALCTDigiCollection > &  alcts,
edm::Handle< CSCCLCTDigiCollection > &  clcts,
edm::Handle< CSCCorrelatedLCTDigiCollection > &  correlatedlcts,
edm::Handle< CSCWireDigiCollection > &  wires,
edm::Handle< CSCStripDigiCollection > &  strips,
edm::Handle< edm::PSimHitContainer > &  simhits,
edm::Handle< CSCRecHit2DCollection > &  rechits,
edm::Handle< CSCSegmentCollection > &  segments,
edm::ESHandle< CSCGeometry > &  cscGeom 
)
private

Definition at line 694 of file CSCEfficiency.cc.

References NumCh.

702  {
703  for(int iE=0;iE<2;iE++){
704  for(int iS=0;iS<4;iS++){
705  for(int iR=0;iR<4;iR++){
706  for(int iC=0;iC<NumCh;iC++){
707  allSegments[iE][iS][iR][iC].clear();
708  allCLCT[iE][iS][iR][iC] = allALCT[iE][iS][iR][iC] = allCorrLCT[iE][iS][iR][iC] = false;
709  for(int iL=0;iL<6;iL++){
710  allStrips[iE][iS][iR][iC][iL].clear();
711  allWG[iE][iS][iR][iC][iL].clear();
712  allRechits[iE][iS][iR][iC][iL].clear();
713  allSimhits[iE][iS][iR][iC][iL].clear();
714  }
715  }
716  }
717  }
718  }
719  //
720  if(useDigis){
721  fillLCT_info(alcts, clcts, correlatedlcts);
722  fillWG_info(wires, cscGeom);
723  fillStrips_info(strips);
724  }
725  fillRechitsSegments_info(rechits, segments, cscGeom);
726  if(!isData){
727  fillSimhit_info(simhits);
728  }
729 }
void fillWG_info(edm::Handle< CSCWireDigiCollection > &wires, edm::ESHandle< CSCGeometry > &cscGeom)
std::vector< std::pair< LocalPoint, int > > allSimhits[2][4][4][NumCh][6]
bool allCLCT[2][4][4][NumCh]
bool allALCT[2][4][4][NumCh]
std::vector< std::pair< std::pair< int, float >, int > > allWG[2][4][4][NumCh][6]
void fillLCT_info(edm::Handle< CSCALCTDigiCollection > &alcts, edm::Handle< CSCCLCTDigiCollection > &clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > &correlatedlcts)
void fillRechitsSegments_info(edm::Handle< CSCRecHit2DCollection > &rechits, edm::Handle< CSCSegmentCollection > &segments, edm::ESHandle< CSCGeometry > &cscGeom)
std::vector< std::pair< int, float > > allStrips[2][4][4][NumCh][6]
std::vector< std::pair< LocalPoint, LocalVector > > allSegments[2][4][4][NumCh]
void fillSimhit_info(edm::Handle< edm::PSimHitContainer > &simHits)
#define NumCh
std::vector< std::pair< LocalPoint, bool > > allRechits[2][4][4][NumCh][6]
bool allCorrLCT[2][4][4][NumCh]
void fillStrips_info(edm::Handle< CSCStripDigiCollection > &strips)
void CSCEfficiency::fillLCT_info ( edm::Handle< CSCALCTDigiCollection > &  alcts,
edm::Handle< CSCCLCTDigiCollection > &  clcts,
edm::Handle< CSCCorrelatedLCTDigiCollection > &  correlatedlcts 
)
private

Definition at line 732 of file CSCEfficiency.cc.

References FirstCh, j, and prof2calltree::last.

734  {
735  //---- ALCTDigis
736  int nSize = 0;
737  for (CSCALCTDigiCollection::DigiRangeIterator j=alcts->begin(); j!=alcts->end(); j++) {
738  ++nSize;
739  const CSCDetId& id = (*j).first;
740  const CSCALCTDigiCollection::Range& range =(*j).second;
742  range.first; digiIt!=range.second;
743  ++digiIt){
744  // Valid digi in the chamber (or in neighbouring chamber)
745  if((*digiIt).isValid()){
746  allALCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
747  }
748  }// for digis in layer
749  }// end of for (j=...
750  ALCTPerEvent->Fill(nSize);
751  //---- CLCTDigis
752  nSize = 0;
753  for (CSCCLCTDigiCollection::DigiRangeIterator j=clcts->begin(); j!=clcts->end(); j++) {
754  ++nSize;
755  const CSCDetId& id = (*j).first;
756  std::vector<CSCCLCTDigi>::const_iterator digiIt = (*j).second.first;
757  std::vector<CSCCLCTDigi>::const_iterator last = (*j).second.second;
758  for( ; digiIt != last; ++digiIt) {
759  // Valid digi in the chamber (or in neighbouring chamber)
760  if((*digiIt).isValid()){
761  allCLCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
762  }
763  }
764  }
765  CLCTPerEvent->Fill(nSize);
766  //---- CorrLCTDigis
767  for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j=correlatedlcts->begin(); j!=correlatedlcts->end(); j++) {
768  const CSCDetId& id = (*j).first;
769  std::vector<CSCCorrelatedLCTDigi>::const_iterator digiIt = (*j).second.first;
770  std::vector<CSCCorrelatedLCTDigi>::const_iterator last = (*j).second.second;
771  for( ; digiIt != last; ++digiIt) {
772  // Valid digi in the chamber (or in neighbouring chamber)
773  if((*digiIt).isValid()){
774  allCorrLCT[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] = true;
775  }
776  }
777  }
778 }
TH1F * CLCTPerEvent
bool allCLCT[2][4][4][NumCh]
bool allALCT[2][4][4][NumCh]
int j
Definition: DBlmapReader.cc:9
#define FirstCh
std::vector< DigiType >::const_iterator const_iterator
std::pair< const_iterator, const_iterator > Range
TH1F * ALCTPerEvent
bool allCorrLCT[2][4][4][NumCh]
void CSCEfficiency::fillRechitsSegments_info ( edm::Handle< CSCRecHit2DCollection > &  rechits,
edm::Handle< CSCSegmentCollection > &  segments,
edm::ESHandle< CSCGeometry > &  cscGeom 
)
private

Definition at line 863 of file CSCEfficiency.cc.

References CSCDetId::chamber(), gather_cfg::cout, CSCDetId, Reference_intrackfit_cff::endcap, CSCDetId::endcap(), FirstCh, CSCDetId::layer(), NumCh, relativeConstraints::ring, CSCDetId::ring(), findQualityFiles::size, relativeConstraints::station, CSCDetId::station(), GeomDet::toGlobal(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

866  {
867  //---- RECHITS AND SEGMENTS
868  //---- Loop over rechits
869  if (printalot){
870  //printf("\tGet the recHits collection.\t ");
871  printf(" The size of the rechit collection is %i\n",int(rechits->size()));
872  //printf("\t...start loop over rechits...\n");
873  }
874  recHitsPerEvent->Fill(rechits->size());
875  //---- Build iterator for rechits and loop :
877  for (recIt = rechits->begin(); recIt != rechits->end(); recIt++) {
878  //---- Find chamber with rechits in CSC
879  CSCDetId id = (CSCDetId)(*recIt).cscDetId();
880  if (printalot){
881  const CSCLayer* csclayer = cscGeom->layer( id);
882  LocalPoint rhitlocal = (*recIt).localPosition();
883  LocalError rerrlocal = (*recIt).localPositionError();
884  GlobalPoint rhitglobal= csclayer->toGlobal(rhitlocal);
885  printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
886  printf("\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n",
887  rhitlocal.x(), rhitlocal.y(), rhitlocal.z(), rerrlocal.xx(), rerrlocal.yy(), rerrlocal.xy(),
888  rhitglobal.x(), rhitglobal.y(), rhitglobal.z());
889  }
890  std::pair <LocalPoint, bool> recHitPos((*recIt).localPosition(), false);
891  allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][id.layer()-1].push_back(recHitPos);
892  }
893  //---- "Empty" chambers
894  for(int iE=0;iE<2;iE++){
895  for(int iS=0;iS<4;iS++){
896  for(int iR=0;iR<4;iR++){
897  for(int iC=0;iC<NumCh;iC++){
898  int numLayers = 0;
899  for(int iL=0;iL<6;iL++){
900  if(allRechits[iE][iS][iR][iC][iL].size()){
901  ++numLayers;
902  }
903  }
904  if(numLayers>1){
905  emptyChambers[iE][iS][iR][iC] = false;
906  }
907  else{
908  emptyChambers[iE][iS][iR][iC] = true;
909  }
910  }
911  }
912  }
913  }
914 
915  //
916  if (printalot){
917  printf(" The size of the segment collection is %i\n", int(segments->size()));
918  //printf("\t...start loop over segments...\n");
919  }
920  segmentsPerEvent->Fill(segments->size());
921  for(CSCSegmentCollection::const_iterator it = segments->begin(); it != segments->end(); it++) {
922  CSCDetId id = (CSCDetId)(*it).cscDetId();
923  StHist[id.endcap()-1][id.station()-1].segmentChi2_ndf->Fill((*it).chi2()/(*it).degreesOfFreedom());
924  StHist[id.endcap()-1][id.station()-1].hitsInSegment->Fill((*it).nRecHits());
925  if (printalot){
926  printf("\tendcap/station/ring/chamber: %i %i %i %i\n",
927  id.endcap(),id.station(),id.ring(),id.chamber());
928  std::cout<<"\tposition(loc) = "<<(*it).localPosition()<<" error(loc) = "<<(*it).localPositionError()<<std::endl;
929  std::cout<<"\t chi2/ndf = "<<(*it).chi2()/(*it).degreesOfFreedom()<<" nhits = "<<(*it).nRecHits() <<std::endl;
930 
931  }
932  allSegments[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh].push_back
933  (make_pair((*it).localPosition(), (*it).localDirection()));
934 
935 
936  //---- try to get the CSC recHits that contribute to this segment.
937  //if (printalot) printf("\tGet the recHits for this segment.\t");
938  std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
939  int nRH = (*it).nRecHits();
940  if (printalot){
941  printf("\tGet the recHits for this segment.\t");
942  printf(" nRH = %i\n",nRH);
943  }
944  //---- Find which of the rechits in the chamber is in the segment
945  int layerRH = 0;
946  for ( vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
947  ++layerRH;
948  CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
949  if(printalot){
950  printf("\t%i RH\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
951  layerRH,idRH.endcap(),idRH.station(),idRH.ring(),idRH.chamber(),idRH.layer());
952  }
953  for(size_t jRH = 0;
954  jRH<allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1].size();
955  ++jRH){
956  LocalPoint lp = allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first;
957  float xDiff = iRH->localPosition().x() -
958  allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.x();
959  float yDiff = iRH->localPosition().y() -
960  allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first.y();
961  if(fabs(xDiff)<0.0001 && fabs(yDiff)<0.0001){
962  std::pair <LocalPoint, bool>
963  recHitPos(allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH].first, true);
964  allRechits[idRH.endcap()-1][idRH.station()-1][idRH.ring()-1][idRH.chamber()-FirstCh][idRH.layer()-1][jRH] = recHitPos;
965  if(printalot){
966  std::cout<<" number of the rechit (from zero) in the segment = "<< jRH<<std::endl;
967  }
968  }
969  }
970  }
971  }
972 }
int chamber() const
Definition: CSCDetId.h:70
struct CSCEfficiency::StationHistos StHist[2][4]
float xx() const
Definition: LocalError.h:19
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
T y() const
Definition: PV3DBase.h:57
int layer() const
Definition: CSCDetId.h:63
int endcap() const
Definition: CSCDetId.h:95
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:45
float xy() const
Definition: LocalError.h:20
float yy() const
Definition: LocalError.h:21
T z() const
Definition: PV3DBase.h:58
int ring() const
Definition: CSCDetId.h:77
#define FirstCh
std::vector< std::pair< LocalPoint, LocalVector > > allSegments[2][4][4][NumCh]
TH1F * segmentsPerEvent
TH1F * recHitsPerEvent
bool emptyChambers[2][4][4][NumCh]
int station() const
Definition: CSCDetId.h:88
tuple cout
Definition: gather_cfg.py:41
#define NumCh
T x() const
Definition: PV3DBase.h:56
std::vector< std::pair< LocalPoint, bool > > allRechits[2][4][4][NumCh][6]
tuple size
Write out results.
void CSCEfficiency::fillSimhit_info ( edm::Handle< edm::PSimHitContainer > &  simHits)
private

Definition at line 852 of file CSCEfficiency.cc.

References CSCDetId::chamber(), CSCDetId, CSCDetId::endcap(), FirstCh, CSCDetId::layer(), CSCDetId::ring(), and CSCDetId::station().

852  {
853  //---- SIMHITS
854  edm::PSimHitContainer::const_iterator dSHsimIter;
855  for (dSHsimIter = simhits->begin(); dSHsimIter != simhits->end(); dSHsimIter++){
856  // Get DetID for this simHit:
857  CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
858  std::pair <LocalPoint, int> simHitPos((*dSHsimIter).localPosition(), (*dSHsimIter).particleType());
859  allSimhits[sId.endcap()-1][sId.station()-1][sId.ring()-1][sId.chamber()-FirstCh][sId.layer()-1].push_back(simHitPos);
860  }
861 }
int chamber() const
Definition: CSCDetId.h:70
std::vector< std::pair< LocalPoint, int > > allSimhits[2][4][4][NumCh][6]
int layer() const
Definition: CSCDetId.h:63
int endcap() const
Definition: CSCDetId.h:95
int ring() const
Definition: CSCDetId.h:77
#define FirstCh
int station() const
Definition: CSCDetId.h:88
void CSCEfficiency::fillStrips_info ( edm::Handle< CSCStripDigiCollection > &  strips)
private

Definition at line 808 of file CSCEfficiency.cc.

References CSCDetId, diffTreeTool::diff, j, prof2calltree::last, and crabWrap::threshold.

808  {
809  //---- STRIPS
810  for (CSCStripDigiCollection::DigiRangeIterator j=strips->begin(); j!=strips->end(); j++) {
811  CSCDetId id = (CSCDetId)(*j).first;
812  int largestADCValue = -1;
813  int largestStrip = -1;
814  std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
815  std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
816  for( ; digiItr != last; ++digiItr) {
817  int maxADC=largestADCValue;
818  int myStrip = digiItr->getStrip();
819  std::vector<int> myADCVals = digiItr->getADCCounts();
820  bool thisStripFired = false;
821  float thisPedestal = 0.5*(float)(myADCVals[0]+myADCVals[1]);
822  float threshold = 13.3 ;
823  float diff = 0.;
824  float peakADC = -1000.;
825  int peakTime = -1;
826  for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
827  diff = (float)myADCVals[iCount]-thisPedestal;
828  if (diff > threshold) {
829  thisStripFired = true;
830  if (myADCVals[iCount] > largestADCValue) {
831  largestADCValue = myADCVals[iCount];
832  largestStrip = myStrip;
833  }
834  }
835  if (diff > threshold && diff > peakADC) {
836  peakADC = diff;
837  peakTime = iCount;
838  }
839  }
840  if(largestADCValue>maxADC){// FIX IT!!!
841  maxADC = largestADCValue;
842  std::pair <int, float> LayerSignal (myStrip, peakADC);
843 
844  //---- AllStrips contains basic information about strips
845  //---- (strip number and peak signal for most significant strip in the layer)
846  allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].clear();
847  allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-1][id.layer()-1].push_back(LayerSignal);
848  }
849  }
850  }
851 }
float threshold
Definition: crabWrap.py:319
int j
Definition: DBlmapReader.cc:9
std::vector< std::pair< int, float > > allStrips[2][4][4][NumCh][6]
void CSCEfficiency::fillWG_info ( edm::Handle< CSCWireDigiCollection > &  wires,
edm::ESHandle< CSCGeometry > &  cscGeom 
)
private

Definition at line 780 of file CSCEfficiency.cc.

References CSCDetId, FirstCh, j, prof2calltree::last, TrapezoidalPlaneBounds::parameters(), and CSCLayerGeometry::yOfWireGroup().

780  {
781  //---- WIRE GROUPS
782  for (CSCWireDigiCollection::DigiRangeIterator j=wires->begin(); j!=wires->end(); j++) {
783  CSCDetId id = (CSCDetId)(*j).first;
784  const CSCLayer *layer_p = cscGeom->layer (id);
785  const CSCLayerGeometry *layerGeom = layer_p->geometry ();
786  //
787 
788  const std::vector<float> LayerBounds = layerGeom->parameters ();
789  std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
790  std::vector<CSCWireDigi>::const_iterator last = (*j).second.second;
791  //
792  for( ; digiItr != last; ++digiItr) {
793  std::pair < int, float > WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup()));
794  std::pair <std::pair < int, float >, int > LayerSignal(WG_pos, digiItr->getTimeBin());
795 
796  //---- AllWG contains basic information about WG (WG number and Y-position, time bin)
797  allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
798  [id.layer()-1].push_back(LayerSignal);
799  if(printalot){
800  //std::cout<<" WG check : "<<std::endl;
801  //printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
802  //std::cout<<" WG size = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
803  //[id.layer()-1].size()<<std::endl;
804  }
805  }
806  }
807 }
float yOfWireGroup(int wireGroup, float x=0.) const
virtual const std::vector< float > parameters() const
int j
Definition: DBlmapReader.cc:9
std::vector< std::pair< std::pair< int, float >, int > > allWG[2][4][4][NumCh][6]
#define FirstCh
bool CSCEfficiency::filter ( edm::Event event,
const edm::EventSetup eventSetup 
)
privatevirtual

Implements edm::EDFilter.

Definition at line 21 of file CSCEfficiency.cc.

References muon::caloCompatibility(), DeDxDiscriminatorTools::charge(), reco::TrackBase::charge(), reco::TrackBase::confirmed, gather_cfg::cout, CSCDetId, debug, deltaR(), MuonPatternRecoDumper::dumpFTS(), MuonPatternRecoDumper::dumpMuonId(), ExpressReco_HICollisions_FallBack::e, Reference_intrackfit_cff::endcap, reco::TrackBase::eta(), FirstCh, reco::Track::found(), TrajectoryStateOnSurface::freeState(), GeomDet::geographicalId(), edm::EventSetup::get(), reco::TrackBase::goodIterative, reco::TrackBase::highPurity, i, iEvent, reco::Track::innerDetId(), reco::Track::innerPosition(), TrajectoryStateOnSurface::isValid(), reco::TrackBase::loose, alignBH_cfg::minP, metsig::muon, ExpressReco_HICollisions_FallBack::muons, reco::TrackBase::normalizedChi2(), reco::Track::outerDetId(), reco::Track::outerMomentum(), reco::Track::outerPosition(), reco::TrackBase::p(), phi, reco::TrackBase::phi(), GloballyPositioned< T >::position(), reco::BeamSpot::position(), FreeTrajectoryState::position(), funct::pow(), edm::Handle< T >::product(), reco::TrackBase::pt(), reco::TrackBase::ptError(), reco::TrackBase::qoverp(), reco::TrackBase::qoverpError(), reco::TrackBase::quality(), reco::TrackBase::qualitySize, relativeConstraints::ring, reco::Muon::SegmentAndTrackArbitration, edm::View< T >::size(), mathSSE::sqrt(), relativeConstraints::station, RecoTauPiZeroBuilderPlugins_cfi::strips, reco::TrackBase::tight, ExpressReco_HICollisions_FallBack::track, ExpressReco_HICollisions_FallBack::trackCollection, edm::TriggerNames::triggerNames(), reco::TrackBase::undefQuality, and PV3DBase< T, PVType, FrameType >::z().

21  {
22  passTheEvent = false;
23  DataFlow->Fill(0.);
25 
26  //---- increment counter
28  // printalot debug output
30  int iRun = event.id().run();
31  int iEvent = event.id().event();
32  if(0==fmod(double (nEventsAnalyzed) ,double(1000) )){
33  if(printalot){
34  printf("\n==enter==CSCEfficiency===== run %i\tevent %i\tn Analyzed %i\n",iRun,iEvent,nEventsAnalyzed);
35  }
36  }
37  theService->update(eventSetup);
38  //---- These declarations create handles to the types of records that you want
39  //---- to retrieve from event "e".
40  if (printalot) printf("\tget handles for digi collections\n");
41 
42  //---- Pass the handle to the method "getByType", which is used to retrieve
43  //---- one and only one instance of the type in question out of event "e". If
44  //---- zero or more than one instance exists in the event an exception is thrown.
45  if (printalot) printf("\tpass handles\n");
53  //edm::Handle<reco::TrackCollection> saMuons;
54  edm::Handle<edm::View<reco::Track> > trackCollectionH;
56 
57  if(useDigis){
58  event.getByLabel(alctDigiTag_, alcts);
59  event.getByLabel(clctDigiTag_, clcts);
60  event.getByLabel(corrlctDigiTag_, correlatedlcts);
61 
62  event.getByLabel( stripDigiTag_, strips);
63  event.getByLabel( wireDigiTag_, wires);
64  }
65  if(!isData){
66  event.getByLabel(simHitTag, simhits);
67  }
68  event.getByLabel(rechitDigiTag_,rechits);
69  event.getByLabel(segmentDigiTag_, segments);
70  //event.getByLabel(saMuonTag,saMuons);
71  event.getByLabel(tracksTag,trackCollectionH);
72  const edm::View<reco::Track> trackCollection = *(trackCollectionH.product());
73 
74  //---- Get the CSC Geometry :
75  if (printalot) printf("\tget the CSC geometry.\n");
77  eventSetup.get<MuonGeometryRecord>().get(cscGeom);
78 
79  // use theTrackingGeometry instead of cscGeom?
80  ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
81  eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
82 
83  bool triggerPassed = true;
84  if(useTrigger){
85  // access the trigger information
86  // trigger names can be find in HLTrigger/Configuration/python/HLT_2E30_cff.py (or?)
87  // get hold of TriggerResults
89  event.getByLabel(hlTriggerResults_,hltR);
90  const edm::TriggerNames & triggerNames = event.triggerNames(*hltR);
91  triggerPassed = applyTrigger(hltR, triggerNames);
92  }
93  if(!triggerPassed){
94  return triggerPassed;
95  }
96  DataFlow->Fill(1.);
97  GlobalPoint gpZero(0.,0.,0.);
98  if(theService->magneticField()->inTesla(gpZero).mag2()<0.1){
99  magField = false;
100  }
101  else{
102  magField = true;
103  }
104 
105  //---- store info from digis
106  fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
107  //
109  edm::InputTag muonTag_("muons");
110  event.getByLabel(muonTag_,muons);
111 
112  edm::Handle<reco::BeamSpot> beamSpotHandle;
113  event.getByLabel("offlineBeamSpot", beamSpotHandle);
114  reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
115  //
116  std::vector <reco::MuonCollection::const_iterator> goodMuons_it;
117  unsigned int nPositiveZ = 0;
118  unsigned int nNegativeZ = 0;
119  float muonOuterZPosition = -99999.;
120  if(isIPdata){
121  if (printalot)std::cout<<" muons.size() = "<<muons->size() <<std::endl;
122  for ( reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon ) {
123  DataFlow->Fill(31.);
124  if (printalot) {
125  std::cout<<" iMuon = "<<muon-muons->begin()<<" charge = "<<muon->charge()<<" p = "<<muon->p()<<" pt = "<<muon->pt()<<
126  " eta = "<<muon->eta()<<" phi = "<<muon->phi()<<
127  " matches = "<<
128  muon->matches().size()<<" matched Seg = "<<muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration)<<" GLB/TR/STA = "<<
129  muon->isGlobalMuon()<<"/"<<muon->isTrackerMuon()<<"/"<<muon->isStandAloneMuon()<<std::endl;
130  }
131  if(!(muon->isTrackerMuon() && muon->isGlobalMuon())){
132  continue;
133  }
134  DataFlow->Fill(32.);
135  double relISO = ( muon->isolationR03().sumPt +
136  muon->isolationR03().emEt +
137  muon->isolationR03().hadEt)/muon->track()->pt();
138  if (printalot) {
139  std::cout<<" relISO = "<<relISO<<" emVetoEt = "<<muon->isolationR03().emVetoEt<<" caloComp = "<<
140  muon::caloCompatibility(*(muon))<<" dxy = "<<fabs(muon->track()->dxy(vertexBeamSpot.position()))<<std::endl;
141  }
142  if(
143  //relISO>0.1 || muon::caloCompatibility(*(muon))<.90 ||
144  fabs(muon->track()->dxy(vertexBeamSpot.position()))>0.2 || muon->pt()<6.){
145  continue;
146  }
147  DataFlow->Fill(33.);
148  if(muon->track()->hitPattern().numberOfValidPixelHits()<1 ||
149  muon->track()->hitPattern().numberOfValidTrackerHits()<11 ||
150  muon->combinedMuon()->hitPattern().numberOfValidMuonHits()<1 ||
151  muon->combinedMuon()->normalizedChi2()>10. ||
152  muon->numberOfMatches()<2){
153  continue;
154  }
155  DataFlow->Fill(34.);
156  float zOuter = muon->combinedMuon()->outerPosition().z();
157  float rhoOuter = muon->combinedMuon()->outerPosition().rho();
158  bool passDepth = true;
159  // barrel region
160  //if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 480.){
161  if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.){
162  passDepth = false;
163  }
164  // endcap region
165  //else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
166  else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
167  passDepth = false;
168  }
169  // overlap region
170  //else if ( fabs(zOuter) > 680. && fabs(zOuter) < 730. && rhoOuter < 480.){
171  else if ( fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.){
172  passDepth = false;
173  }
174  if(!passDepth){
175  continue;
176  }
177  DataFlow->Fill(35.);
178  goodMuons_it.push_back(muon);
179  if(muon->track()->momentum().z()>0.){
180  ++nPositiveZ;
181  }
182  if(muon->track()->momentum().z()<0.){
183  ++nNegativeZ;
184  }
185  }
186  }
187 
188  //
189 
190 
191  if (printalot) std::cout<<"Start track loop over "<<trackCollection.size()<<" tracks"<<std::endl;
192  for(edm::View<reco::Track>::size_type i=0; i<trackCollection.size(); ++i) {
193  DataFlow->Fill(2.);
194  edm::RefToBase<reco::Track> track(trackCollectionH, i);
195  //std::cout<<" iTR = "<<i<<" eta = "<<track->eta()<<" phi = "<<track->phi()<<std::cout<<" pt = "<<track->pt()<<std::endl;
196  if(isIPdata){
197  if (printalot){
198  std::cout<<" nNegativeZ = "<<nNegativeZ<<" nPositiveZ = "<<nPositiveZ<<std::endl;
199  }
200  if(nNegativeZ>1 || nPositiveZ>1){
201  break;
202  }
203  bool trackOK = false;
204  if (printalot){
205  std::cout<<" goodMuons_it.size() = "<<goodMuons_it.size()<<std::endl;
206  }
207  for(size_t iM=0;iM<goodMuons_it.size();++iM){
208  //std::cout<<" iM = "<<iM<<" eta = "<<goodMuons_it[iM]->track()->eta()<<
209  //" phi = "<<goodMuons_it[iM]->track()->phi()<<
210  //" pt = "<<goodMuons_it[iM]->track()->pt()<<std::endl;
211  float deltaR = pow(track->phi()-goodMuons_it[iM]->track()->phi(),2) +
212  pow(track->eta()-goodMuons_it[iM]->track()->eta(),2);
213  deltaR = sqrt(deltaR);
214  if (printalot){
215  std::cout<<" TR mu match to a tr: deltaR = "<<deltaR<<" dPt = "<<
216  track->pt()-goodMuons_it[iM]->track()->pt()<<std::endl;
217  }
218  if(deltaR>0.01 || fabs(track->pt()-goodMuons_it[iM]->track()->pt())>0.1 ){
219  continue;
220  }
221  else{
222  trackOK = true;
223  if (printalot){
224  std::cout<<" trackOK "<<std::endl;
225  }
226  muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
227  break;
228  //++nChosenTracks;
229  }
230  }
231  if(!trackOK){
232  if (printalot){
233  std::cout<<" failed: trackOK "<<std::endl;
234  }
235  continue;
236  }
237  }
238  else{
239  //---- Do we need a better "clean track" definition?
240  if(trackCollection.size()>2){
241  break;
242  }
243  DataFlow->Fill(3.);
244  if(!i && 2==trackCollection.size()){
246  edm::RefToBase<reco::Track> trackTwo(trackCollectionH, tType);
247  if(track->outerPosition().z()*trackTwo->outerPosition().z()>0){// in one and the same "endcap"
248  break;
249  }
250  }
251  }
252  DataFlow->Fill(4.);
253  if (printalot){
254  std::cout<<"i track = "<<i<<" P = "<<track->p()<<" chi2/ndf = "<<track->normalizedChi2()<<" nSeg = "<<segments->size()<<std::endl;
255  std::cout<<"quality undef/loose/tight/high/confirmed/goodIt/size "<<
256  track->quality(reco::Track::undefQuality)<<"/"<<
257  track->quality(reco::Track::loose)<<"/"<<
258  track->quality(reco::Track::tight)<<"/"<<
259  track->quality(reco::Track::highPurity)<<"/"<<
260  track->quality(reco::Track::confirmed)<<"/"<<
261  track->quality(reco::Track::goodIterative)<<"/"<<
262  track->quality(reco::Track::qualitySize)<<
263  std::endl;
264  std::cout<<" pt = "<< track->pt()<<" +-"<<track->ptError()<<" q/pt = "<<track->qoverp()<<" +- "<<track->qoverpError()<<std::endl;
265  //std::cout<<" const Pmin = "<<minTrackMomentum<<" pMax = "<<maxTrackMomentum<<" maxNormChi2 = "<<maxNormChi2<<std::endl;
266  std::cout<<" track inner position = "<<track->innerPosition()<<" outer position = "<<track->outerPosition()<<std::endl;
267  std::cout<<"track eta (outer) = "<<track->outerPosition().eta()<<" phi (outer) = "<<
268  track->outerPosition().phi()<<std::endl;
269  if(fabs(track->innerPosition().z())>500.){
270  DetId innerDetId(track->innerDetId());
271  std::cout<<" dump inner state MUON detid = "<<debug.dumpMuonId(innerDetId)<<std::endl;
272  }
273  if(fabs(track->outerPosition().z())>500.){
274  DetId outerDetId(track->outerDetId());
275  std::cout<<" dump outer state MUON detid = "<<debug.dumpMuonId(outerDetId)<<std::endl;
276  }
277 
278  std::cout<<" nHits = "<<track->found()<<std::endl;
279  /*
280  trackingRecHit_iterator rhbegin = track->recHitsBegin();
281  trackingRecHit_iterator rhend = track->recHitsEnd();
282  int iRH = 0;
283  for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){
284  const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
285  std::cout<<"hit "<<iRH<<" loc pos = " <<(*recHit)->localPosition()<<
286  " glob pos = " <<geomDet->toGlobal((*recHit)->localPosition())<<std::endl;
287  ++iRH;
288  }
289  */
290  }
291  float dpT_ov_pT = 0.;
292  if(fabs(track->pt())>0.001){
293  dpT_ov_pT = track->ptError()/ track->pt();
294  }
295  //---- These define a "good" track
296  if(track->normalizedChi2()>maxNormChi2){// quality
297  break;
298  }
299  DataFlow->Fill(5.);
300  if(track->found()<minTrackHits){// enough data points
301  break;
302  }
303  DataFlow->Fill(6.);
304  if(!segments->size()){// better have something in the CSC
305  break;
306  }
307  DataFlow->Fill(7.);
308  if(magField && (track->p()<minP || track->p()>maxP)){// proper energy range
309  break;
310  }
311  DataFlow->Fill(8.);
312  if(magField && (dpT_ov_pT >0.5) ){// not too crazy uncertainty
313  break;
314  }
315  DataFlow->Fill(9.);
316 
317  passTheEvent = true;
318  if (printalot) std::cout<<"good Track"<<std::endl;
319  CLHEP::Hep3Vector r3T_inner(track->innerPosition().x(),track->innerPosition().y(),track->innerPosition().z());
320  CLHEP::Hep3Vector r3T(track->outerPosition().x(),track->outerPosition().y(),track->outerPosition().z());
321  chooseDirection(r3T_inner, r3T);// for non-IP
322 
323  CLHEP::Hep3Vector p3T(track->outerMomentum().x(),track->outerMomentum().y(),track->outerMomentum().z());
324  CLHEP::Hep3Vector p3_propagated, r3_propagated;
325  AlgebraicSymMatrix66 cov_propagated, covT;
326  covT *= 1e-20;
327  cov_propagated *= 1e-20;
328  int charge = track->charge();
329  FreeTrajectoryState ftsStart = getFromCLHEP(p3T, r3T, charge, covT, &*(theService->magneticField()));
330  if (printalot){
331  std::cout<<" p = "<<track->p()<<" norm chi2 = "<<track->normalizedChi2()<<std::endl;
332  std::cout<<" dump the very first FTS = "<<debug.dumpFTS(ftsStart)<<std::endl;
333  }
334  TrajectoryStateOnSurface tSOSDest;
335  int endcap = 0;
336  //---- which endcap to look at
337  if(track->outerPosition().z()>0){
338  endcap = 1;
339  }
340  else{
341  endcap = 2;
342  }
343  int chamber = 1;
344  //---- a "refference" CSCDetId for each ring
345  std::vector< CSCDetId > refME;
346  for(int iS=1;iS<5;++iS){
347  for(int iR=1;iR<4;++iR){
348  if(1!=iS && iR>2){
349  continue;
350  }
351  else if(4==iS && iR>1){
352  continue;
353  }
354  refME.push_back( CSCDetId(endcap, iS, iR, chamber));
355  }
356  }
357  //---- loop over the "refference" CSCDetIds
358  for(size_t iSt = 0; iSt<refME.size();++iSt){
359  if (printalot){
360  std::cout<<"loop iStatation = "<<iSt<<std::endl;
361  std::cout<<"refME[iSt]: st = "<<refME[iSt].station()<<" rg = "<<refME[iSt].ring()<<std::endl;
362  }
363  std::map <std::string, bool> chamberTypes;
364  chamberTypes["ME11"] = false;
365  chamberTypes["ME12"] = false;
366  chamberTypes["ME13"] = false;
367  chamberTypes["ME21"] = false;
368  chamberTypes["ME22"] = false;
369  chamberTypes["ME31"] = false;
370  chamberTypes["ME32"] = false;
371  chamberTypes["ME41"] = false;
372  const CSCChamber* cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
373  DetId detId = cscChamber_base->geographicalId();
374  if (printalot){
375  std::cout<<" base iStation : eta = "<<cscGeom->idToDet(detId)->surface().position().eta()<<" phi = "<<
376  cscGeom->idToDet(detId)->surface().position().phi() << " y = " <<cscGeom->idToDet(detId)->surface().position().y()<<std::endl;
377  std::cout<<" dump base iStation detid = "<<debug.dumpMuonId(detId)<<std::endl;
378  std::cout<<" dump FTS start = "<<debug.dumpFTS(ftsStart)<<std::endl;
379  }
380  //---- propagate to this ME
381  tSOSDest = propagate(ftsStart, cscGeom->idToDet(detId)->surface());
382  if(tSOSDest.isValid()){
383  ftsStart = *tSOSDest.freeState();
384  if (printalot) std::cout<<" dump FTS end = "<<debug.dumpFTS(ftsStart)<<std::endl;
385  getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
386  float feta = fabs(r3_propagated.eta());
387  float phi = r3_propagated.phi();
388  //---- which rings are (possibly) penetrated
389  ringCandidates(refME[iSt].station(), feta, chamberTypes);
390 
391  map<std::string,bool>::iterator iter;
392  int iterations = 0;
393  //---- loop over ring candidates
394  for( iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++ ) {
395  ++iterations;
396  //---- is this ME a machinig candidate station
397  if(iter->second && (iterations-1)==int(iSt)){
398  if (printalot){
399  std::cout<<" Chamber type "<< iter->first<<" is a candidate..."<<std::endl;
400  std::cout<<" station() = "<< refME[iSt].station()<<" ring() = "<<refME[iSt].ring()<<" iSt = "<<iSt<<std::endl;
401  }
402  std::vector <int> coupleOfChambers;
403  //---- which chamber (and its closes neighbor) is penetrated by the track - candidates
404  chamberCandidates(refME[iSt].station(), refME[iSt].ring(), phi, coupleOfChambers);
405  //---- loop over the two chamber candidates
406  for(size_t iCh =0;iCh<coupleOfChambers.size();++iCh){
407  DataFlow->Fill(11.);
408  if (printalot) std::cout<<" Check chamber N = "<<coupleOfChambers.at(iCh)<<std::endl;;
409  if((!getAbsoluteEfficiency) && (true == emptyChambers
410  [refME[iSt].endcap()-1]
411  [refME[iSt].station()-1]
412  [refME[iSt].ring()-1]
413  [coupleOfChambers.at(iCh)-FirstCh])){
414  continue;
415  }
416  CSCDetId theCSCId(refME[iSt].endcap(), refME[iSt].station(), refME[iSt].ring(), coupleOfChambers.at(iCh));
417  const CSCChamber* cscChamber = cscGeom->chamber(theCSCId.chamberId());
418  const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
419  float zFTS = ftsStart.position().z();
420  float dz = fabs(bpCh.position().z() - zFTS);
421  float zDistInner = track->innerPosition().z() - bpCh.position().z();
422  float zDistOuter = track->outerPosition().z() - bpCh.position().z();
423  //---- only detectors between the inner and outer points of the track are considered for non IP-data
424  if(printalot){
425  std::cout<<" zIn = "<<track->innerPosition().z()<<" zOut = "<<track->outerPosition().z()<<" zSurf = "<<bpCh.position().z()<<std::endl;
426  }
427  if(!isIPdata && (zDistInner*zDistOuter>0. || fabs(zDistInner)<15. || fabs(zDistOuter)<15.)){ // for non IP-data
428  if(printalot){
429  std::cout<<" Not an intermediate (as defined) point... Skip."<<std::endl;
430  }
431  continue;
432  }
433  if(isIPdata && fabs(track->eta())<1.8){
434  if(fabs(muonOuterZPosition) - fabs(bpCh.position().z())<0 ||
435  fabs(muonOuterZPosition-bpCh.position().z())<15.){
436  continue;
437  }
438  }
439  DataFlow->Fill(13.);
440  //---- propagate to the chamber (from this ME) if it is a different surface (odd/even chambers)
441  if(dz>0.1){// i.e. non-zero (float 0 check is bad)
442  //if(fabs(zChanmber - zFTS ) > 0.1){
443  tSOSDest = propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface());
444  if(tSOSDest.isValid()){
445  ftsStart = *tSOSDest.freeState();
446  }
447  else{
448  if(printalot) std::cout<<"TSOS not valid! Break."<<std::endl;
449  break;
450  }
451  }
452  else{
453  if(printalot) std::cout<<" info: dz<0.1"<<std::endl;
454  }
455  DataFlow->Fill(15.);
456  FreeTrajectoryState ftsInit = ftsStart;
457  bool inDeadZone = false;
458  //---- loop over the 6 layers
459  for(int iLayer = 0;iLayer<6;++iLayer){
460  bool extrapolationPassed = true;
461  if (printalot){
462  std::cout<<" iLayer = "<<iLayer<<" dump FTS init = "<<debug.dumpFTS(ftsInit)<<std::endl;
463  std::cout<<" dump detid = "<<debug.dumpMuonId(cscChamber->geographicalId())<<std::endl;
464  std::cout<<"Surface to propagate to: pos = "<<cscChamber->layer(iLayer+1)->surface().position()<<" eta = "
465  <<cscChamber->layer(iLayer+1)->surface().position().eta()<<" phi = "
466  <<cscChamber->layer(iLayer+1)->surface().position().phi()<<std::endl;
467  }
468  //---- propagate to this layer
469  tSOSDest = propagate(ftsInit, cscChamber->layer(iLayer+1)->surface());
470  if(tSOSDest.isValid()){
471  ftsInit = *tSOSDest.freeState();
472  if (printalot) std::cout<<" Propagation between layers successful: dump FTS end = "<<debug.dumpFTS(ftsInit)<<std::endl;
473  getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
474  }
475  else{
476  if (printalot) std::cout<<"Propagation between layers not successful - notValid TSOS"<<std::endl;
477  extrapolationPassed = false;
478  inDeadZone = true;
479  }
480  //}
481  //---- Extrapolation passed? For each layer?
482  if(extrapolationPassed){
483  GlobalPoint theExtrapolationPoint(r3_propagated.x(),r3_propagated.y(),r3_propagated.z());
484  LocalPoint theLocalPoint = cscChamber->layer(iLayer+1)->toLocal(theExtrapolationPoint);
485  //std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<std::endl;
486  inDeadZone = ( inDeadZone ||
487  !inSensitiveLocalRegion(theLocalPoint.x(), theLocalPoint.y(),
488  refME[iSt].station(), refME[iSt].ring()));
489  if (printalot){
490  std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<"inDeadZone = "<<inDeadZone<<std::endl;
491  }
492  //---- break if in dead zone for any layer ("clean" tracks)
493  if(inDeadZone){
494  break;
495  }
496  }
497  else{
498  break;
499  }
500  }
501  DataFlow->Fill(17.);
502  //---- Is a track in a sensitive area for each layer?
503  if(!inDeadZone){//---- for any layer
504  DataFlow->Fill(19.);
505  if (printalot) std::cout<<"Do efficiencies..."<<std::endl;
506  //---- Do efficiencies
507  // angle cuts applied (if configured)
508  bool angle_flag = true; angle_flag = efficienciesPerChamber(theCSCId, cscChamber, ftsStart);
509  if(useDigis && angle_flag){
510  bool stripANDwire_flag; stripANDwire_flag = stripWire_Efficiencies(theCSCId, ftsStart);
511  }
512  if(angle_flag){
513  bool recHitANDsegment_flag; recHitANDsegment_flag = recHitSegment_Efficiencies(theCSCId, cscChamber, ftsStart);
514  if(!isData){
515  recSimHitEfficiency(theCSCId, ftsStart);
516  }
517  }
518  }
519  else{
520  if(printalot) std::cout<<" Not in active area for all layers"<<std::endl;
521  }
522  }
523  if(tSOSDest.isValid()){
524  ftsStart = *tSOSDest.freeState();
525  }
526  }
527  }
528  }
529  else{
530  if (printalot) std::cout<<" TSOS not valid..."<<std::endl;
531  }
532  }
533  }
534  //---- End
535  if (printalot) printf("==exit===CSCEfficiency===== run %i\tevent %i\n\n",iRun,iEvent);
536  return passTheEvent;
537 }
edm::InputTag corrlctDigiTag_
int i
Definition: DBlmapReader.cc:9
edm::InputTag segmentDigiTag_
bool applyTrigger(edm::Handle< edm::TriggerResults > &hltR, const edm::TriggerNames &triggerNames)
void chooseDirection(CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition)
void ringCandidates(int station, float absEta, std::map< std::string, bool > &chamberTypes)
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
double charge(const std::vector< uint8_t > &Ampls)
Strings const & triggerNames() const
Definition: TriggerNames.cc:24
bool getAbsoluteEfficiency
float caloCompatibility(const reco::Muon &muon)
FreeTrajectoryState getFromCLHEP(const CLHEP::Hep3Vector &p3, const CLHEP::Hep3Vector &r3, int charge, const AlgebraicSymMatrix66 &cov, const MagneticField *field)
std::string dumpMuonId(const DetId &id) const
std::string dumpFTS(const FreeTrajectoryState &fts) const
int iEvent
Definition: GenABIO.cc:243
MuonServiceProxy * theService
FreeTrajectoryState * freeState(bool withErrors=true) const
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
unsigned int size_type
Definition: View.h:90
edm::InputTag rechitDigiTag_
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:74
edm::InputTag hlTriggerResults_
void getFromFTS(const FreeTrajectoryState &fts, CLHEP::Hep3Vector &p3, CLHEP::Hep3Vector &r3, int &charge, AlgebraicSymMatrix66 &cov)
bool efficienciesPerChamber(CSCDetId &id, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
void chamberCandidates(int station, int ring, float phi, std::vector< int > &coupleOfChambers)
double maxNormChi2
bool stripWire_Efficiencies(CSCDetId &cscDetId, FreeTrajectoryState &ftsChamber)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
Definition: DetId.h:20
GlobalPoint position() const
void fillDigiInfo(edm::Handle< CSCALCTDigiCollection > &alcts, edm::Handle< CSCCLCTDigiCollection > &clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > &correlatedlcts, edm::Handle< CSCWireDigiCollection > &wires, edm::Handle< CSCStripDigiCollection > &strips, edm::Handle< edm::PSimHitContainer > &simhits, edm::Handle< CSCRecHit2DCollection > &rechits, edm::Handle< CSCSegmentCollection > &segments, edm::ESHandle< CSCGeometry > &cscGeom)
#define FirstCh
const T & get() const
Definition: EventSetup.h:55
edm::InputTag tracksTag
size_type size() const
bool recHitSegment_Efficiencies(CSCDetId &cscDetId, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
unsigned int printout_NEvents
T const * product() const
Definition: Handle.h:74
edm::InputTag clctDigiTag_
bool inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring)
bool emptyChambers[2][4][4][NumCh]
edm::InputTag simHitTag
tuple cout
Definition: gather_cfg.py:41
const Point & position() const
position
Definition: BeamSpot.h:63
edm::InputTag stripDigiTag_
#define debug
Definition: MEtoEDMFormat.h:34
edm::InputTag alctDigiTag_
edm::InputTag wireDigiTag_
const PositionType & position() const
unsigned int minTrackHits
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
TrajectoryStateOnSurface propagate(FreeTrajectoryState &ftsStart, const BoundPlane &bp)
Definition: DDAxes.h:10
FreeTrajectoryState CSCEfficiency::getFromCLHEP ( const CLHEP::Hep3Vector &  p3,
const CLHEP::Hep3Vector &  r3,
int  charge,
const AlgebraicSymMatrix66 cov,
const MagneticField field 
)
private

Definition at line 1466 of file CSCEfficiency.cc.

1468  {
1469 
1470  GlobalVector p3GV(p3.x(), p3.y(), p3.z());
1471  GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
1472  GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
1473 
1474  CartesianTrajectoryError tCov(cov);
1475 
1476  return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars) ;
1477 }
double charge(const std::vector< uint8_t > &Ampls)
double p3[4]
Definition: TauolaWrapper.h:91
void CSCEfficiency::getFromFTS ( const FreeTrajectoryState fts,
CLHEP::Hep3Vector &  p3,
CLHEP::Hep3Vector &  r3,
int &  charge,
AlgebraicSymMatrix66 cov 
)
private

Definition at line 1451 of file CSCEfficiency.cc.

References FreeTrajectoryState::cartesianError(), FreeTrajectoryState::charge(), FreeTrajectoryState::hasError(), CartesianTrajectoryError::matrix(), FreeTrajectoryState::momentum(), FreeTrajectoryState::position(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

1453  {
1454 
1455  GlobalVector p3GV = fts.momentum();
1456  GlobalPoint r3GP = fts.position();
1457 
1458  p3.set(p3GV.x(), p3GV.y(), p3GV.z());
1459  r3.set(r3GP.x(), r3GP.y(), r3GP.z());
1460 
1461  charge = fts.charge();
1462  cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
1463 
1464 }
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
T y() const
Definition: PV3DBase.h:57
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
T z() const
Definition: PV3DBase.h:58
GlobalVector momentum() const
const AlgebraicSymMatrix66 & matrix() const
GlobalPoint position() const
const CartesianTrajectoryError & cartesianError() const
T x() const
Definition: PV3DBase.h:56
double p3[4]
Definition: TauolaWrapper.h:91
bool CSCEfficiency::inSensitiveLocalRegion ( double  xLocal,
double  yLocal,
int  station,
int  ring 
)
private

Definition at line 540 of file CSCEfficiency.cc.

References abs.

540  {
541  //---- Good region means sensitive area of a chamber. "Local" stands for the local system
542  bool pass = false;
543  std::vector <double> chamberBounds(3);// the sensitive area
544  float y_center = 99999.;
545  //---- hardcoded... not good
546  if(station>1 && station<5){
547  if(2==ring){
548  chamberBounds[0] = 66.46/2; // (+-)x1 shorter
549  chamberBounds[1] = 127.15/2; // (+-)x2 longer
550  chamberBounds[2] = 323.06/2;
551  y_center = -0.95;
552  }
553  else{
554  if(2==station){
555  chamberBounds[0] = 54.00/2; // (+-)x1 shorter
556  chamberBounds[1] = 125.71/2; // (+-)x2 longer
557  chamberBounds[2] = 189.66/2;
558  y_center = -0.955;
559  }
560  else if(3==station){
561  chamberBounds[0] = 61.40/2; // (+-)x1 shorter
562  chamberBounds[1] = 125.71/2; // (+-)x2 longer
563  chamberBounds[2] = 169.70/2;
564  y_center = -0.97;
565  }
566  else if(4==station){
567  chamberBounds[0] = 69.01/2; // (+-)x1 shorter
568  chamberBounds[1] = 125.65/2; // (+-)x2 longer
569  chamberBounds[2] = 149.42/2;
570  y_center = -0.94;
571  }
572  }
573  }
574  else if(1==station){
575  if(3==ring){
576  chamberBounds[0] = 63.40/2; // (+-)x1 shorter
577  chamberBounds[1] = 92.10/2; // (+-)x2 longer
578  chamberBounds[2] = 164.16/2;
579  y_center = -1.075;
580  }
581  else if(2==ring){
582  chamberBounds[0] = 51.00/2; // (+-)x1 shorter
583  chamberBounds[1] = 83.74/2; // (+-)x2 longer
584  chamberBounds[2] = 174.49/2;
585  y_center = -0.96;
586  }
587  else{// to be investigated
588  chamberBounds[0] = 30./2;//40./2; // (+-)x1 shorter
589  chamberBounds[1] = 60./2;//100./2; // (+-)x2 longer
590  chamberBounds[2] = 160./2;//142./2;
591  y_center = 0.;
592  }
593  }
594  double yUp = chamberBounds[2] + y_center;
595  double yDown = - chamberBounds[2] + y_center;
596  double xBound1Shifted = chamberBounds[0]-distanceFromDeadZone;//
597  double xBound2Shifted = chamberBounds[1]-distanceFromDeadZone;//
598  double lineSlope = (yUp - yDown)/(xBound2Shifted-xBound1Shifted);
599  double lineConst = yUp - lineSlope*xBound2Shifted;
600  double yBoundary = lineSlope*abs(xLocal) + lineConst;
601  pass = checkLocal(yLocal, yBoundary, station, ring);
602  return pass;
603 }
#define abs(x)
Definition: mlp_lapack.h:159
double distanceFromDeadZone
bool checkLocal(double yLocal, double yBoundary, int station, int ring)
void CSCEfficiency::linearExtrapolation ( GlobalPoint  initialPosition,
GlobalVector  initialDirection,
float  zSurface,
std::vector< float > &  posZY 
)
private

Definition at line 1479 of file CSCEfficiency.cc.

References PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

1480  {
1481  double paramLine = lineParameter(initialPosition.z(), zSurface, initialDirection.z());
1482  double xPosition = extrapolate1D(initialPosition.x(), initialDirection.x(),paramLine);
1483  double yPosition = extrapolate1D(initialPosition.y(), initialDirection.y(),paramLine);
1484  posZY.clear();
1485  posZY.push_back(xPosition);
1486  posZY.push_back(yPosition);
1487 }
double lineParameter(double initZPosition, double destZPosition, double initZDirection)
T y() const
Definition: PV3DBase.h:57
T z() const
Definition: PV3DBase.h:58
double extrapolate1D(double initPosition, double initDirection, double parameterOfTheLine)
T x() const
Definition: PV3DBase.h:56
double CSCEfficiency::lineParameter ( double  initZPosition,
double  destZPosition,
double  initZDirection 
)
private

Definition at line 1494 of file CSCEfficiency.cc.

1494  {
1495  double paramLine = (destZPosition-initZPosition)/initZDirection;
1496  return paramLine;
1497 }
TrajectoryStateOnSurface CSCEfficiency::propagate ( FreeTrajectoryState ftsStart,
const BoundPlane bp 
)
private

Definition at line 1529 of file CSCEfficiency.cc.

References LargeD0_PixelPairStep_cff::propagator.

1529  {
1530  TrajectoryStateOnSurface tSOSDest;
1531  std::string propagatorName;
1532 /*
1533 // it would work if cosmic muons had properly assigned direction...
1534  bool dzPositive = bpDest.position().z() - ftsStart.position().z() > 0 ? true : false;
1535  //---- Be careful with trigger conditions too
1536  if(!isIPdata){
1537  bool rightDirection = !(alongZ^dzPositive);
1538  if(rightDirection){
1539  if(printalot) std::cout<<" propagate along momentum"<<std::endl;
1540  propagatorName = "SteppingHelixPropagatorAlong";
1541  }
1542  else{
1543  if(printalot) std::cout<<" propagate opposite momentum"<<std::endl;
1544  propagatorName = "SteppingHelixPropagatorOpposite";
1545  }
1546  }
1547  else{
1548  if(printalot) std::cout<<" propagate any (momentum)"<<std::endl;
1549  propagatorName = "SteppingHelixPropagatorAny";
1550  }
1551 */
1552  propagatorName = "SteppingHelixPropagatorAny";
1553  tSOSDest = propagator(propagatorName)->propagate(ftsStart, bpDest);
1554  return tSOSDest;
1555 }
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:9
const Propagator * propagator(std::string propagatorName) const
const Propagator * CSCEfficiency::propagator ( std::string  propagatorName) const
private

Definition at line 1524 of file CSCEfficiency.cc.

1524  {
1525  return &*theService->propagator(propagatorName);
1526 }
MuonServiceProxy * theService
bool CSCEfficiency::recHitSegment_Efficiencies ( CSCDetId cscDetId,
const CSCChamber cscChamber,
FreeTrajectoryState ftsChamber 
)
private

Definition at line 1310 of file CSCEfficiency.cc.

References gather_cfg::cout, Reference_intrackfit_cff::endcap, PV3DBase< T, PVType, FrameType >::eta(), first, FirstCh, CSCLayer::geometry(), CSCChamber::layer(), M_PI, CSCLayerGeometry::nearestStrip(), GloballyPositioned< T >::position(), FreeTrajectoryState::position(), funct::pow(), relativeConstraints::ring, edm::second(), findQualityFiles::size, mathSSE::sqrt(), relativeConstraints::station, CSCLayerGeometry::stripAngle(), GeomDet::surface(), GeomDet::toGlobal(), GeomDet::toLocal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

1310  {
1311  int ec, st, rg, ch, secondRing;
1312  returnTypes(id, ec, st, rg, ch, secondRing);
1313  bool firstCondition, secondCondition;
1314 
1315  std::vector <bool> missingLayers_rh(6);
1316  std::vector <int> usedInSegment(6);
1317  // Rechits
1318  if(printalot) std::cout<<"RecHits eff"<<std::endl;
1319  for(int iLayer=0;iLayer<6;++iLayer){
1320  firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
1321  secondCondition = false;
1322  int thisRing = rg;
1323  if(secondRing>-1){
1324  secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
1325  if(secondCondition){
1326  thisRing = secondRing;
1327  }
1328  }
1329  if(firstCondition || secondCondition){
1330  ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(iLayer+1);
1331  for(size_t iR=0;
1332  iR<allRechits[ec][st][thisRing][ch][iLayer].size();
1333  ++iR){
1334  if(allRechits[ec][st][thisRing][ch][iLayer][iR].second){
1335  usedInSegment[iLayer] = 1;
1336  break;
1337  }
1338  else{
1339  usedInSegment[iLayer] = -1;
1340  }
1341  }
1342  }
1343  else{
1344  missingLayers_rh[iLayer] = true;
1345  if(printalot){
1346  std::cout<<"missing rechits ";
1347  printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
1348  }
1349  }
1350  }
1351  GlobalVector globalDir;
1352  GlobalPoint globalPos;
1353  // Segments
1354  firstCondition = allSegments[ec][st][rg][ch].size() ? true : false;
1355  secondCondition = false;
1356  int secondSize = 0;
1357  int thisRing = rg;
1358  if(secondRing>-1){
1359  secondCondition = allSegments[ec][st][secondRing][ch].size() ? true : false;
1360  secondSize = allSegments[ec][st][secondRing][ch].size();
1361  if(secondCondition){
1362  thisRing = secondRing;
1363  }
1364  }
1365  if(firstCondition || secondCondition){
1366  if (printalot) std::cout<<"segments - start ec = "<<ec<<" st = "<<st<<" rg = "<<rg<<" ch = "<<ch<<std::endl;
1367  StHist[ec][st].EfficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
1368  if(1==allSegments[ec][st][rg][ch].size() + secondSize){
1369  globalDir = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].second);
1370  globalPos = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].first);
1371  StHist[ec][st].EfficientSegments_eta->Fill(fabs(ftsChamber.position().eta()));
1372  double residual = sqrt(pow(ftsChamber.position().x() - globalPos.x(),2)+
1373  pow(ftsChamber.position().y() - globalPos.y(),2)+
1374  pow(ftsChamber.position().z() - globalPos.z(),2));
1375  if (printalot) std::cout<<" fts.position() = "<<ftsChamber.position()<<" segPos = "<<globalPos<<" res = "<<residual<< std::endl;
1376  StHist[ec][st].ResidualSegments->Fill(residual);
1377  }
1378  for(int iLayer=0;iLayer<6;++iLayer){
1379  if(printalot) std::cout<<" iLayer = "<<iLayer<<" usedInSegment = "<<usedInSegment[iLayer]<<std::endl;
1380  if(0!=usedInSegment[iLayer]){
1381  if(-1==usedInSegment[iLayer]){
1382  ChHist[ec][st][rg][ch].InefficientSingleHits->Fill(iLayer+1);
1383  }
1384  ChHist[ec][st][rg][ch].AllSingleHits->Fill(iLayer+1);
1385  }
1386  firstCondition = allRechits[ec][st][rg][ch][iLayer].size() ? true : false;
1387  secondCondition = false;
1388  if(secondRing>-1){
1389  secondCondition = allRechits[ec][st][secondRing][ch][iLayer].size() ? true : false;
1390  }
1391  float stripAngle = 99999.;
1392  std::vector<float> posXY(2);
1393  bool oneSegment = false;
1394  if(1==allSegments[ec][st][rg][ch].size() + secondSize){
1395  oneSegment = true;
1396  const BoundPlane bp = cscChamber->layer(iLayer+1)->surface();
1397  linearExtrapolation(globalPos,globalDir, bp.position().z(), posXY);
1398  GlobalPoint gp_extrapol( posXY.at(0), posXY.at(1),bp.position().z());
1399  const LocalPoint lp_extrapol = cscChamber->layer(iLayer+1)->toLocal(gp_extrapol);
1400  posXY.at(0) = lp_extrapol.x();
1401  posXY.at(1) = lp_extrapol.y();
1402  int nearestStrip = cscChamber->layer(iLayer+1)->geometry()->nearestStrip(lp_extrapol);
1403  stripAngle = cscChamber->layer(iLayer+1)->geometry()->stripAngle(nearestStrip) - M_PI/2. ;
1404  }
1405  if(firstCondition || secondCondition){
1406  ChHist[ec][st][rg][ch].EfficientRechits_inSegment->Fill(iLayer+1);
1407  if(oneSegment){
1408  ChHist[ec][st][rg][ch].Y_EfficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
1409  ChHist[ec][st][rg][ch].Phi_EfficientRecHits_inSegment[iLayer]->Fill(stripAngle);
1410  }
1411  }
1412  else{
1413  if(oneSegment){
1414  ChHist[ec][st][rg][ch].Y_InefficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
1415  ChHist[ec][st][rg][ch].Phi_InefficientRecHits_inSegment[iLayer]->Fill(stripAngle);
1416  }
1417  }
1418  }
1419  }
1420  else{
1421  StHist[ec][st].InefficientSegments_XY->Fill(ftsChamber.position().x(),ftsChamber.position().y());
1422  if(printalot){
1423  std::cout<<"missing segment "<<std::endl;
1424  printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber());
1425  std::cout<<" fts.position() = "<<ftsChamber.position()<<std::endl;
1426  }
1427  }
1428  // Normalization
1429  ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(8);
1430  if(allSegments[ec][st][rg][ch].size()+secondSize<2){
1431  StHist[ec][st].AllSegments_eta->Fill(fabs(ftsChamber.position().eta()));
1432  }
1433  ChHist[ec][st][rg][id.chamber()-FirstCh].EfficientRechits_inSegment->Fill(9);
1434 
1435  return true;
1436 }
struct CSCEfficiency::StationHistos StHist[2][4]
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
std::vector< TH1F * > Y_InefficientRecHits_inSegment
T y() const
Definition: PV3DBase.h:57
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:64
void returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing)
U second(std::pair< T, U > const &p)
std::vector< TH1F * > Phi_InefficientRecHits_inSegment
std::vector< TH1F * > Y_EfficientRecHits_inSegment
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
Definition: CSCChamber.cc:41
void linearExtrapolation(GlobalPoint initialPosition, GlobalVector initialDirection, float zSurface, std::vector< float > &posZY)
bool first
Definition: L1TdeRCT.cc:79
std::vector< TH1F * > Phi_EfficientRecHits_inSegment
int nearestStrip(const LocalPoint &lp) const
GlobalPoint position() const
#define M_PI
Definition: BFit3D.cc:3
#define FirstCh
std::vector< std::pair< LocalPoint, LocalVector > > allSegments[2][4][4][NumCh]
T eta() const
Definition: PV3DBase.h:70
tuple cout
Definition: gather_cfg.py:41
T x() const
Definition: PV3DBase.h:56
const PositionType & position() const
std::vector< std::pair< LocalPoint, bool > > allRechits[2][4][4][NumCh][6]
float stripAngle(int strip) const
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
const CSCLayerGeometry * geometry() const
Definition: CSCLayer.h:47
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
bool CSCEfficiency::recSimHitEfficiency ( CSCDetId id,
FreeTrajectoryState ftsChamber 
)
private

Definition at line 1263 of file CSCEfficiency.cc.

References edm::second(), and findQualityFiles::size.

1263  {
1264  int ec, st, rg, ch, secondRing;
1265  returnTypes(id, ec, st, rg, ch, secondRing);
1266  bool firstCondition, secondCondition;
1267  for(int iLayer=0; iLayer<6;iLayer++){
1268  firstCondition = allSimhits[ec][st][rg][ch][iLayer].size() ? true : false;
1269  secondCondition = false;
1270  int thisRing = rg;
1271  if(secondRing>-1){
1272  secondCondition = allSimhits[ec][st][secondRing][ch][iLayer].size() ? true : false;
1273  if(secondCondition){
1274  thisRing = secondRing;
1275  }
1276  }
1277  if(firstCondition || secondCondition){
1278  for(size_t iSH=0;
1279  iSH<allSimhits[ec][st][thisRing][ch][iLayer].size();
1280  ++iSH){
1281  if(13 ==
1282  fabs(allSimhits[ec][st][thisRing][ch][iLayer][iSH].second)){
1283  ChHist[ec][st][rg][ch].SimSimhits->Fill(iLayer+1);
1284  if(allRechits[ec][st][thisRing][ch][iLayer].size()){
1285  ChHist[ec][st][rg][ch].SimRechits->Fill(iLayer+1);
1286  }
1287  break;
1288  }
1289  }
1290  //---- Next is not too usefull...
1291  /*
1292  for(unsigned int iSimHits=0;
1293  iSimHits<allSimhits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
1294  iSimHits++){
1295  ChHist[ec][st][rg][id.chamber()-FirstCh].SimSimhits_each->Fill(iLayer+1);
1296  }
1297  for(unsigned int iRecHits=0;
1298  iRecHits<allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
1299  iRecHits++){
1300  ChHist[ec][st][rg][id.chamber()-FirstCh].SimRechits_each->Fill(iLayer+1);
1301  }
1302  */
1303  //
1304  }
1305  }
1306  return true;
1307 }
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
std::vector< std::pair< LocalPoint, int > > allSimhits[2][4][4][NumCh][6]
void returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing)
U second(std::pair< T, U > const &p)
std::vector< std::pair< LocalPoint, bool > > allRechits[2][4][4][NumCh][6]
tuple size
Write out results.
void CSCEfficiency::returnTypes ( CSCDetId id,
int &  ec,
int &  st,
int &  rg,
int &  ch,
int &  secondRing 
)
private

Definition at line 1438 of file CSCEfficiency.cc.

References FirstCh, relativeConstraints::ring, and relativeConstraints::station.

1438  {
1439  ec = id.endcap()-1;
1440  st = id.station()-1;
1441  rg = id.ring()-1;
1442  secondRing = -1;
1443  if(1==id.station() && (4==id.ring() || 1==id.ring()) ){
1444  rg = 0;
1445  secondRing = 3;
1446  }
1447  ch = id.chamber()-FirstCh;
1448 }
#define FirstCh
void CSCEfficiency::ringCandidates ( int  station,
float  absEta,
std::map< std::string, bool > &  chamberTypes 
)
private

Definition at line 975 of file CSCEfficiency.cc.

975  {
976  // yeah, hardcoded again...
977  switch (station){
978  case 1:
979  if(feta>0.85 && feta<1.18){//ME13
980  chamberTypes["ME13"] = true;
981  }
982  if(feta>1.18 && feta<1.7){//ME12
983  chamberTypes["ME12"] = true;
984  }
985  if(feta>1.5 && feta<2.45){//ME11
986  chamberTypes["ME11"] = true;
987  }
988  break;
989  case 2:
990  if(feta>0.95 && feta<1.6){//ME22
991  chamberTypes["ME22"] = true;
992 
993  }
994  if(feta>1.55 && feta<2.45){//ME21
995  chamberTypes["ME21"] = true;
996  }
997  break;
998  case 3:
999  if(feta>1.08 && feta<1.72){//ME32
1000  chamberTypes["ME32"] = true;
1001 
1002  }
1003  if(feta>1.69 && feta<2.45){//ME31
1004  chamberTypes["ME31"] = true;
1005  }
1006  break;
1007  case 4:
1008  if(feta>1.78 && feta<2.45){//ME41
1009  chamberTypes["ME41"] = true;
1010  }
1011  break;
1012  default:
1013  break;
1014  }
1015 }
bool CSCEfficiency::stripWire_Efficiencies ( CSCDetId cscDetId,
FreeTrajectoryState ftsChamber 
)
private

Definition at line 1187 of file CSCEfficiency.cc.

References gather_cfg::cout, Reference_intrackfit_cff::endcap, FirstCh, FreeTrajectoryState::momentum(), relativeConstraints::ring, relativeConstraints::station, and PV3DBase< T, PVType, FrameType >::theta().

1187  {
1188  int ec, st, rg, ch, secondRing;
1189  returnTypes(id, ec, st, rg, ch, secondRing);
1190 
1191  bool firstCondition, secondCondition;
1192  int missingLayers_s = 0;
1193  int missingLayers_wg = 0;
1194  for(int iLayer=0;iLayer<6;iLayer++){
1195  //----Strips
1196  if(printalot){
1197  printf("\t%i swEff: \tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
1198  iLayer + 1,id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
1199  std::cout<<" size S = "<<allStrips[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<
1200  "size W = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size()<<std::endl;
1201 
1202  }
1203  firstCondition = allStrips[ec][st][rg][ch][iLayer].size() ? true : false;
1204  //allSegments[ec][st][rg][ch].size() ? true : false;
1205  secondCondition = false;
1206  if(secondRing>-1){
1207  secondCondition = allStrips[ec][st][secondRing][ch][iLayer].size() ? true : false;
1208  }
1209  if(firstCondition || secondCondition){
1210  ChHist[ec][st][rg][ch].EfficientStrips->Fill(iLayer+1);
1211  }
1212  else{
1213  if(printalot){
1214  std::cout<<"missing strips ";
1215  printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
1216  }
1217  }
1218  // Wires
1219  firstCondition = allWG[ec][st][rg][ch][iLayer].size() ? true : false;
1220  secondCondition = false;
1221  if(secondRing>-1){
1222  secondCondition = allWG[ec][st][secondRing][ch][iLayer].size() ? true : false;
1223  }
1224  if(firstCondition || secondCondition){
1225  ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(iLayer+1);
1226  }
1227  else{
1228  if(printalot){
1229  std::cout<<"missing wires ";
1230  printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),iLayer+1);
1231  }
1232  }
1233  }
1234  // Normalization
1235  if(6!=missingLayers_s){
1236  ChHist[ec][st][rg][ch].EfficientStrips->Fill(8);
1237  }
1238  if(6!=missingLayers_wg){
1239  ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(8);
1240  }
1241  ChHist[ec][st][rg][ch].EfficientStrips->Fill(9);
1242  ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(9);
1243 //
1244  ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(1);
1245  if(missingLayers_s!=missingLayers_wg){
1246  ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(2);
1247  if(6==missingLayers_wg){
1248  ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(3);
1249  ChHist[ec][st][rg][ch].NoWires_momTheta->Fill(ftsChamber.momentum().theta());
1250  }
1251  if(6==missingLayers_s){
1252  ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(4);
1253  ChHist[ec][st][rg][ch].NoStrips_momPhi->Fill(ftsChamber.momentum().theta());
1254  }
1255  }
1256  else if(6==missingLayers_s){
1257  ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(5);
1258  }
1259 
1260  return true;
1261 }
struct CSCEfficiency::ChamberHistos ChHist[2][4][3][LastCh-FirstCh+1]
void returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing)
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
std::vector< std::pair< std::pair< int, float >, int > > allWG[2][4][4][NumCh][6]
GlobalVector momentum() const
std::vector< std::pair< int, float > > allStrips[2][4][4][NumCh][6]
#define FirstCh
tuple cout
Definition: gather_cfg.py:41

Member Data Documentation

edm::InputTag CSCEfficiency::alctDigiTag_
private

Definition at line 135 of file CSCEfficiency.h.

TH1F* CSCEfficiency::ALCTPerEvent
private

Definition at line 267 of file CSCEfficiency.h.

bool CSCEfficiency::allALCT[2][4][4][NumCh]
private

Definition at line 189 of file CSCEfficiency.h.

bool CSCEfficiency::allCLCT[2][4][4][NumCh]
private

Definition at line 188 of file CSCEfficiency.h.

bool CSCEfficiency::allCorrLCT[2][4][4][NumCh]
private

Definition at line 190 of file CSCEfficiency.h.

std::vector<std::pair <LocalPoint, bool> > CSCEfficiency::allRechits[2][4][4][NumCh][6]
private

Definition at line 204 of file CSCEfficiency.h.

std::vector<std::pair <LocalPoint, LocalVector> > CSCEfficiency::allSegments[2][4][4][NumCh]
private

Definition at line 207 of file CSCEfficiency.h.

std::vector<std::pair <LocalPoint, int> > CSCEfficiency::allSimhits[2][4][4][NumCh][6]
private

Definition at line 200 of file CSCEfficiency.h.

std::vector<std::pair <int, float> > CSCEfficiency::allStrips[2][4][4][NumCh][6]
private

Definition at line 193 of file CSCEfficiency.h.

std::vector<std::pair <std::pair <int, float>, int> > CSCEfficiency::allWG[2][4][4][NumCh][6]
private

Definition at line 196 of file CSCEfficiency.h.

bool CSCEfficiency::alongZ
private

Definition at line 182 of file CSCEfficiency.h.

bool CSCEfficiency::andOr
private

Definition at line 168 of file CSCEfficiency.h.

bool CSCEfficiency::applyIPangleCuts
private

Definition at line 157 of file CSCEfficiency.h.

struct CSCEfficiency::ChamberHistos CSCEfficiency::ChHist[2][4][3][LastCh-FirstCh+1]
private
edm::InputTag CSCEfficiency::clctDigiTag_
private

Definition at line 136 of file CSCEfficiency.h.

TH1F* CSCEfficiency::CLCTPerEvent
private

Definition at line 268 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::corrlctDigiTag_
private

Definition at line 137 of file CSCEfficiency.h.

TH1F* CSCEfficiency::DataFlow
private

Definition at line 264 of file CSCEfficiency.h.

double CSCEfficiency::distanceFromDeadZone
private

Definition at line 151 of file CSCEfficiency.h.

bool CSCEfficiency::emptyChambers[2][4][4][NumCh]
private

Definition at line 210 of file CSCEfficiency.h.

bool CSCEfficiency::getAbsoluteEfficiency
private

Definition at line 149 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::hlTriggerResults_
private

Definition at line 165 of file CSCEfficiency.h.

bool CSCEfficiency::isBeamdata
private

Definition at line 148 of file CSCEfficiency.h.

bool CSCEfficiency::isData
private

Definition at line 146 of file CSCEfficiency.h.

bool CSCEfficiency::isIPdata
private

Definition at line 147 of file CSCEfficiency.h.

double CSCEfficiency::local_DX_DZ_Max
private

Definition at line 160 of file CSCEfficiency.h.

double CSCEfficiency::local_DY_DZ_Max
private

Definition at line 158 of file CSCEfficiency.h.

double CSCEfficiency::local_DY_DZ_Min
private

Definition at line 159 of file CSCEfficiency.h.

bool CSCEfficiency::magField
private

Definition at line 180 of file CSCEfficiency.h.

double CSCEfficiency::maxNormChi2
private

Definition at line 154 of file CSCEfficiency.h.

double CSCEfficiency::maxP
private

Definition at line 153 of file CSCEfficiency.h.

double CSCEfficiency::minP
private

Definition at line 152 of file CSCEfficiency.h.

unsigned int CSCEfficiency::minTrackHits
private

Definition at line 155 of file CSCEfficiency.h.

std::vector<std::string> CSCEfficiency::myTriggers
private

Definition at line 166 of file CSCEfficiency.h.

int CSCEfficiency::nEventsAnalyzed
private

Definition at line 178 of file CSCEfficiency.h.

bool CSCEfficiency::passTheEvent
private

Definition at line 184 of file CSCEfficiency.h.

std::vector<int> CSCEfficiency::pointToTriggers
private

Definition at line 167 of file CSCEfficiency.h.

bool CSCEfficiency::printalot
private

Definition at line 176 of file CSCEfficiency.h.

unsigned int CSCEfficiency::printout_NEvents
private

Definition at line 145 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::rechitDigiTag_
private

Definition at line 140 of file CSCEfficiency.h.

TH1F* CSCEfficiency::recHitsPerEvent
private

Definition at line 269 of file CSCEfficiency.h.

std::string CSCEfficiency::rootFileName
private

Definition at line 133 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::segmentDigiTag_
private

Definition at line 142 of file CSCEfficiency.h.

TH1F* CSCEfficiency::segmentsPerEvent
private

Definition at line 270 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::simHitTag
private

Definition at line 141 of file CSCEfficiency.h.

struct CSCEfficiency::StationHistos CSCEfficiency::StHist[2][4]
private
edm::InputTag CSCEfficiency::stripDigiTag_
private

Definition at line 138 of file CSCEfficiency.h.

TFile* CSCEfficiency::theFile
private

Definition at line 174 of file CSCEfficiency.h.

MuonServiceProxy* CSCEfficiency::theService
private

Definition at line 172 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::tracksTag
private

Definition at line 143 of file CSCEfficiency.h.

TH1F* CSCEfficiency::TriggersFired
private

Definition at line 265 of file CSCEfficiency.h.

bool CSCEfficiency::useDigis
private

Definition at line 150 of file CSCEfficiency.h.

bool CSCEfficiency::useTrigger
private

Definition at line 163 of file CSCEfficiency.h.

edm::InputTag CSCEfficiency::wireDigiTag_
private

Definition at line 139 of file CSCEfficiency.h.