TkMSParametrization/test/TestMS.cc

0001 // system includes
0002 #include <iostream>
0003
0004 // user includes
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
0016 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0017 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0018 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0019 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0020 #include "TrackingTools/GeomPropagators/interface/StraightLinePropagator.h"
0021 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0022 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0025
0026 // ROOT includes
0027 #include "TFile.h"
0028 #include "TH1D.h"
0029 #include "TProfile.h"
0030
0031 using namespace GeomDetEnumerators;
0032 using namespace std;
0033 using namespace edm;
0034
0035 template <class T>
0036 T sqr(T t) {
0037   return t * t;
0038 }
0039
0040 class TestMS : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0041 public:
0042   explicit TestMS(const edm::ParameterSet& conf);
0043   ~TestMS();
0044   virtual void beginRun(edm::Run const& run, const edm::EventSetup& es) override;
0045   virtual void analyze(const edm::Event& ev, const edm::EventSetup& es) override;
0046   virtual void endRun(edm::Run const& run, const edm::EventSetup& es) override;
0047
0048 private:
0049   edm::ESGetToken<GeometricSearchTracker, TrackerRecoGeometryRecord> trackerToken_;
0050   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bfieldToken_;
0051   TFile* rootFile;
0052   TH1F *hB1, *hB2, *hB3, *hF1, *hF2;
0053 };
0054
0055 TestMS::TestMS(const edm::ParameterSet& conf) : trackerToken_(esConsumes()), bfieldToken_(esConsumes()) {
0056   rootFile = new TFile("histos.root", "RECREATE");
0057   const int nch = 250;
0058   const float rmin = 0;
0059   const float zmin = 0;
0060   const float zmax = 25;
0061   const float rmax = 15;
0062
0063   edm::LogInfo("TestMS") << " CRATEING HISTOS";
0064   hF1 = new TH1F("hSigmaF1", "hSigmaF1", nch, rmin, rmax);
0065   hF2 = new TH1F("hSigmaF2", "hSigmaF2", nch, rmin, rmax);
0066   hB1 = new TH1F("hSigmaB1", "hSigmaB1", nch, zmin, zmax);
0067   hB2 = new TH1F("hSigmaB2", "hSigmaB2", nch, zmin, zmax);
0068   hB3 = new TH1F("hSigmaB3", "hSigmaB3", nch, zmin, zmax);
0069 }
0070
0071 TestMS::~TestMS() {
0072   edm::LogInfo("TestMS") << " DTOR";
0073   std::cout << "WRITING ROOT FILE" << std::endl;
0074   rootFile->Write();
0075   std::cout << "rootFile WRITTEN" << std::endl;
0076 }
0077
0078 void TestMS::analyze(const edm::Event& ev, const edm::EventSetup& es) {}
0079 void TestMS::endRun(edm::Run const& run, const edm::EventSetup& es) {}
0080
0081 void TestMS::beginRun(edm::Run const& run, const edm::EventSetup& es) {
0082   auto const& tracker = es.getData(trackerToken_);
0083   auto const& bfield = es.getData(bfieldToken_);
0084
0085   vector<BarrelDetLayer const*> barrel = tracker.barrelLayers();
0086   //  vector<ForwardDetLayer*> endcap=tracker->posForwardLayers();
0087   vector<ForwardDetLayer const*> endcap = tracker.negForwardLayers();
0088
0089   MultipleScatteringParametrisationMaker maker(tracker, bfield);
0090
0091   MultipleScatteringParametrisation sb1 = maker.parametrisation(barrel[0]);
0092   MultipleScatteringParametrisation sb2 = maker.parametrisation(barrel[1]);
0093   MultipleScatteringParametrisation sb3 = maker.parametrisation(barrel[2]);
0094
0095   MultipleScatteringParametrisation sf1 = maker.parametrisation(endcap[0]);
0096   MultipleScatteringParametrisation sf2 = maker.parametrisation(endcap[1]);
0097
0098   const int nch = 250;
0099   const float rmin = 0;
0100   const float zmin = 0;
0101   const float zmax = 25;
0102   const float rmax = 15;
0103   float r1 = 4.43;
0104   float r2 = 7.34;
0105   float r3 = 10.2;
0106   float z1 = 34.5;
0107   float z2 = 46.5;  // float z3 = 58.5;
0108   z1 = -z1;
0109   z2 = -z2;  // z3 = -z3;
0110
0111   cout << "HERE !!!! " << endl;
0112   for (int i = 0; i < nch; i++) {
0113     float r = rmin + (i + 0.5) * (rmax - rmin) / nch;
0114     float z = zmin + (i + 0.5) * (zmax - zmin) / nch;
0115     float cotThetaB1 = z / r1;
0116     float cotThetaB2 = z / r2;
0117     float cotThetaB3 = z / r3;
0118     float cotThetaF1 = z1 / r;
0119     float cotThetaF2 = z2 / r;
0120     float msB1 = sb1(1., cotThetaB1);
0121     float msB2 = sb2(1., cotThetaB2);
0122     float msB3 = sb3(1., cotThetaB3);
0123     float msF1 = sf1(1., cotThetaF1);
0124     float msF2 = sf2(1., cotThetaF2);
0125     //    std::cout <<"--------------------->"
0126     //            <<i<<" r="<<r<<" cot="<<cotThetaF2<<" ms="<<msF1<<" "<<msF2<<std::endl;
0127     //
0128     hB1->Fill(fabs(z), msB1);
0129     hB2->Fill(fabs(z), msB2);
0130     hB3->Fill(fabs(z), msB3);
0131     hF1->Fill(r, msF1);
0132     hF2->Fill(r, msF2);
0133   }
0134 }
0135
0136 DEFINE_FWK_MODULE(TestMS);