CSCLayerGeometryInside

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
/// CSCLayerGeometryInside.cc
/// Test CSCLayerGeometry::inside fiducial function
/// Tim Cox - 06.05.2009

#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Geometry/CSCGeometry/interface/CSCGeometry.h"
#include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
#include "Geometry/CSCGeometry/interface/CSCLayer.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"

#include "TH2F.h"
#include "TFile.h"
#include "TRandom3.h"

#include <string>

class CSCLayerGeometryInside : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit CSCLayerGeometryInside(const edm::ParameterSet&);
  ~CSCLayerGeometryInside() override;

  void beginJob() override {}
  void analyze(edm::Event const&, edm::EventSetup const&) override;
  void endJob() override {}

  const std::string& myName() { return myName_; }

private:
  const int dashedLineWidth_;
  const std::string dashedLine_;
  const std::string myName_;

  const int ntries;  // no. of points to generate (set in config)

  // No,. of bins and range of x and y generation (set in config)
  const int nbx;
  const int nby;
  const double xlo;
  const double xhi;
  const double ylo;
  const double yhi;

  const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> tokGeom_;

  //      TH1F* h1as;
  //      TH1F* h1bs;

  TH2F* h1ai;
  TH2F* h1ao;
  TH2F* h1bi;
  TH2F* h1bo;
  TH2F* hall;

  TRandom3* tro;
};

CSCLayerGeometryInside::CSCLayerGeometryInside(const edm::ParameterSet& ps)
    : dashedLineWidth_(194),
      dashedLine_(std::string(dashedLineWidth_, '-')),
      myName_("CSCLayerGeometryInside"),
      ntries(ps.getUntrackedParameter<int>("ntries")),
      nbx(ps.getUntrackedParameter<int>("nbx")),
      nby(ps.getUntrackedParameter<int>("nby")),
      xlo(ps.getUntrackedParameter<double>("xlo")),
      xhi(ps.getUntrackedParameter<double>("xhi")),
      ylo(ps.getUntrackedParameter<double>("ylo")),
      yhi(ps.getUntrackedParameter<double>("yhi")),
      tokGeom_(esConsumes()) {
  usesResource("TFileService");

  std::cout << myName_ << " constructor:" << std::endl;
  std::cout << "x range is " << nbx << " bins from " << xlo << " to " << xhi << std::endl;
  std::cout << "y range is " << nby << " bins from " << ylo << " to " << yhi << std::endl;

  edm::Service<TFileService> fs;

  tro = new TRandom3();

  //  h1as = fs->make<TH1F>( "h1as", "ME1a strip", 800, 0., 80. );
  //  h1bs = fs->make<TH1F>( "h1bs", "ME1b strip", 800, 0., 80. );

  hall = fs->make<TH2F>("hall", "All      y vs x", nbx, xlo, xhi, nby, ylo, yhi);
  h1ai = fs->make<TH2F>("h1ai", "ME1a  in y vs x", nbx, xlo, xhi, nby, ylo, yhi);
  h1bi = fs->make<TH2F>("h1bi", "ME1b  in y vs x", nbx, xlo, xhi, nby, ylo, yhi);
  h1ao = fs->make<TH2F>("h1ao", "ME1a out y vs x", nbx, xlo, xhi, nby, ylo, yhi);
  h1bo = fs->make<TH2F>("h1bo", "ME1b out y vs x", nbx, xlo, xhi, nby, ylo, yhi);
}

CSCLayerGeometryInside::~CSCLayerGeometryInside() { delete tro; }

void CSCLayerGeometryInside::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  std::cout << myName() << ": Analyzer..." << std::endl;
  std::cout << "start " << dashedLine_ << std::endl;

  // Build the geometry
  const edm::ESHandle<CSCGeometry>& pgeom = iSetup.getHandle(tokGeom_);

  std::cout << "CSCGeometry contains " << pgeom->chambers().size() << " chambers" << std::endl;
  std::cout << "CSCGeometry contains " << pgeom->layers().size() << " layers" << std::endl;

  std::cout << dashedLine_ << std::endl;

  // Construct an ME1a and a ME1b layer
  CSCDetId id1a(1, 1, 4, 1, 1);
  CSCDetId id1b(1, 1, 1, 1, 1);
  // Get a pointer to each layer
  const CSCLayer* p1a = pgeom->layer(id1a);
  const CSCLayer* p1b = pgeom->layer(id1b);
  // Get a pointer to each layer's layergeometry
  const CSCLayerGeometry* lg1a = p1a->geometry();
  const CSCLayerGeometry* lg1b = p1b->geometry();

  // Now generate a whole load of random LocalPoints and check whether they'e inside or outside the strip region

  //   float fstrips1a = lg1a->numberOfStrips();
  //   float fstrips1b = lg1b->numberOfStrips();
  //   float epsilon = 1.0e-06;

  for (int i = 1; i <= ntries; ++i) {
    float x = tro->Uniform(xlo, xhi);
    float y = tro->Uniform(ylo, yhi);
    LocalPoint lp(x, y);

    // Examine the float strip value returned by the StripTopology
    //     float strip1a = lg1a->topology()->strip( lp );
    //     float strip1b = lg1b->topology()->strip( lp );

    //     if ( strip1a > epsilon && strip1a < fstrips1a-epsilon ) h1as->Fill( strip1a );
    //     if ( strip1b > epsilon && strip1b < fstrips1b-epsilon ) h1bs->Fill( strip1b );

    // histogram each generated point
    hall->Fill(x, y);

    // histogram according to inside/outside 1a and 1b
    if (lg1a->inside(lp))
      h1ai->Fill(x, y);
    else
      h1ao->Fill(x, y);

    if (lg1b->inside(lp))
      h1bi->Fill(x, y);
    else
      h1bo->Fill(x, y);
  }

  std::cout << dashedLine_ << " end" << std::endl;
}

// This is a plug-in
DEFINE_FWK_MODULE(CSCLayerGeometryInside);