Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-08-09 23:47:33 

0001 /*! \brief   Checklist
0002  *  \details TTClusters and TTStubs
0003  *
0004  *  \author Nicola Pozzobon
0005  *  \author Sebastien Viret
0006  *  \date   2013, Jul 22
0007  *
0008  */
0009 
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/PluginManager/interface/ModuleDef.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0017 #include "MagneticField/Engine/interface/MagneticField.h"
0018 #include "DataFormats/L1TrackTrigger/interface/TTTypes.h"
0019 #include "DataFormats/L1TrackTrigger/interface/TTCluster.h"
0020 #include "DataFormats/L1TrackTrigger/interface/TTStub.h"
0021 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0022 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0023 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0024 #include "Geometry/CommonTopologies/interface/Topology.h"
0025 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0026 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0027 #include <TH1D.h>
0028 #include <TH2D.h>
0029 
0030 class AnalyzerClusterStub : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0031   /// Public methods
0032 public:
0033   /// Constructor/destructor
0034   explicit AnalyzerClusterStub(const edm::ParameterSet& iConfig);
0035   ~AnalyzerClusterStub() override;
0036   // Typical methods used on Loops over events
0037   void beginJob() override;
0038   void endJob() override;
0039   void beginRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0040   void endRun(const edm::Run& iEvent, const edm::EventSetup& iSetup) override {}
0041   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0042 
0043   /// Private methods and variables
0044 private:
0045   /// TrackingParticle and TrackingVertex
0046   TH2D* hSimVtx_XY;
0047   TH2D* hSimVtx_RZ;
0048 
0049   TH1D* hTPart_Pt;
0050   TH1D* hTPart_Eta_Pt10;
0051   TH1D* hTPart_Phi_Pt10;
0052 
0053   /// Global positions of TTClusters
0054   TH2D* hCluster_Barrel_XY;
0055   TH2D* hCluster_Barrel_XY_Zoom;
0056   TH2D* hCluster_Endcap_Fw_XY;
0057   TH2D* hCluster_Endcap_Bw_XY;
0058   TH2D* hCluster_RZ;
0059   TH2D* hCluster_Endcap_Fw_RZ_Zoom;
0060   TH2D* hCluster_Endcap_Bw_RZ_Zoom;
0061 
0062   TH1D* hCluster_IMem_Barrel;
0063   TH1D* hCluster_IMem_Endcap;
0064   TH1D* hCluster_OMem_Barrel;
0065   TH1D* hCluster_OMem_Endcap;
0066 
0067   TH1D* hCluster_Gen_Barrel;
0068   TH1D* hCluster_Unkn_Barrel;
0069   TH1D* hCluster_Comb_Barrel;
0070   TH1D* hCluster_Gen_Endcap;
0071   TH1D* hCluster_Unkn_Endcap;
0072   TH1D* hCluster_Comb_Endcap;
0073 
0074   TH1D* hCluster_Gen_Eta;
0075   TH1D* hCluster_Unkn_Eta;
0076   TH1D* hCluster_Comb_Eta;
0077 
0078   TH2D* hCluster_PID;
0079   TH2D* hCluster_W;
0080 
0081   TH1D* hTPart_Eta_INormalization;
0082   TH1D* hTPart_Eta_ICW_1;
0083   TH1D* hTPart_Eta_ICW_2;
0084   TH1D* hTPart_Eta_ICW_3;
0085   TH1D* hTPart_Eta_ONormalization;
0086   TH1D* hTPart_Eta_OCW_1;
0087   TH1D* hTPart_Eta_OCW_2;
0088   TH1D* hTPart_Eta_OCW_3;
0089 
0090   /// Global positions of TTStubs
0091   TH2D* hStub_Barrel_XY;
0092   TH2D* hStub_Barrel_XY_Zoom;
0093   TH2D* hStub_Endcap_Fw_XY;
0094   TH2D* hStub_Endcap_Bw_XY;
0095   TH2D* hStub_RZ;
0096   TH2D* hStub_Endcap_Fw_RZ_Zoom;
0097   TH2D* hStub_Endcap_Bw_RZ_Zoom;
0098 
0099   TH1D* hStub_Barrel;
0100   TH1D* hStub_Endcap;
0101 
0102   TH1D* hStub_Gen_Barrel;
0103   TH1D* hStub_Unkn_Barrel;
0104   TH1D* hStub_Comb_Barrel;
0105   TH1D* hStub_Gen_Endcap;
0106   TH1D* hStub_Unkn_Endcap;
0107   TH1D* hStub_Comb_Endcap;
0108 
0109   TH1D* hStub_Gen_Eta;
0110   TH1D* hStub_Unkn_Eta;
0111   TH1D* hStub_Comb_Eta;
0112 
0113   TH1D* hStub_PID;
0114   TH2D* hStub_Barrel_W;
0115   TH2D* hStub_Barrel_O;
0116   TH2D* hStub_Endcap_W;
0117   TH2D* hStub_Endcap_O;
0118 
0119   /// Stub finding coverage
0120   TH1D* hTPart_Eta_Pt10_Normalization;
0121   TH1D* hTPart_Eta_Pt10_NumPS;
0122   TH1D* hTPart_Eta_Pt10_Num2S;
0123 
0124   /// Denominator for Stub Prod Eff
0125   std::map<unsigned int, TH1D*> mapCluLayer_hTPart_Pt;
0126   std::map<unsigned int, TH1D*> mapCluLayer_hTPart_Eta_Pt10;
0127   std::map<unsigned int, TH1D*> mapCluLayer_hTPart_Phi_Pt10;
0128   std::map<unsigned int, TH1D*> mapCluDisk_hTPart_Pt;
0129   std::map<unsigned int, TH1D*> mapCluDisk_hTPart_Eta_Pt10;
0130   std::map<unsigned int, TH1D*> mapCluDisk_hTPart_Phi_Pt10;
0131   /// Numerator for Stub Prod Eff
0132   std::map<unsigned int, TH1D*> mapStubLayer_hTPart_Pt;
0133   std::map<unsigned int, TH1D*> mapStubLayer_hTPart_Eta_Pt10;
0134   std::map<unsigned int, TH1D*> mapStubLayer_hTPart_Phi_Pt10;
0135   std::map<unsigned int, TH1D*> mapStubDisk_hTPart_Pt;
0136   std::map<unsigned int, TH1D*> mapStubDisk_hTPart_Eta_Pt10;
0137   std::map<unsigned int, TH1D*> mapStubDisk_hTPart_Phi_Pt10;
0138 
0139   /// Comparison of Stubs to TrackingParticles
0140   std::map<unsigned int, TH2D*> mapStubLayer_hStub_InvPt_TPart_InvPt;
0141   std::map<unsigned int, TH2D*> mapStubLayer_hStub_Pt_TPart_Pt;
0142   std::map<unsigned int, TH2D*> mapStubLayer_hStub_Eta_TPart_Eta;
0143   std::map<unsigned int, TH2D*> mapStubLayer_hStub_Phi_TPart_Phi;
0144   std::map<unsigned int, TH2D*> mapStubDisk_hStub_InvPt_TPart_InvPt;
0145   std::map<unsigned int, TH2D*> mapStubDisk_hStub_Pt_TPart_Pt;
0146   std::map<unsigned int, TH2D*> mapStubDisk_hStub_Eta_TPart_Eta;
0147   std::map<unsigned int, TH2D*> mapStubDisk_hStub_Phi_TPart_Phi;
0148 
0149   /// Residuals
0150   std::map<unsigned int, TH2D*> mapStubLayer_hStub_InvPtRes_TPart_Eta;
0151   std::map<unsigned int, TH2D*> mapStubLayer_hStub_PtRes_TPart_Eta;
0152   std::map<unsigned int, TH2D*> mapStubLayer_hStub_EtaRes_TPart_Eta;
0153   std::map<unsigned int, TH2D*> mapStubLayer_hStub_PhiRes_TPart_Eta;
0154   std::map<unsigned int, TH2D*> mapStubDisk_hStub_InvPtRes_TPart_Eta;
0155   std::map<unsigned int, TH2D*> mapStubDisk_hStub_PtRes_TPart_Eta;
0156   std::map<unsigned int, TH2D*> mapStubDisk_hStub_EtaRes_TPart_Eta;
0157   std::map<unsigned int, TH2D*> mapStubDisk_hStub_PhiRes_TPart_Eta;
0158 
0159   /// Stub Width vs Pt
0160   std::map<unsigned int, TH2D*> mapStubLayer_hStub_W_TPart_Pt;
0161   std::map<unsigned int, TH2D*> mapStubLayer_hStub_W_TPart_InvPt;
0162   std::map<unsigned int, TH2D*> mapStubDisk_hStub_W_TPart_Pt;
0163   std::map<unsigned int, TH2D*> mapStubDisk_hStub_W_TPart_InvPt;
0164 
0165   /// Containers of parameters passed by python
0166   /// configuration file
0167   edm::ParameterSet config;
0168 
0169   bool testedGeometry;
0170   bool DebugMode;
0171 };
0172 
0173 //////////////////////////////////
0174 //                              //
0175 //     CLASS IMPLEMENTATION     //
0176 //                              //
0177 //////////////////////////////////
0178 
0179 //////////////
0180 // CONSTRUCTOR
0181 AnalyzerClusterStub::AnalyzerClusterStub(edm::ParameterSet const& iConfig) : config(iConfig) {
0182   /// Insert here what you need to initialize
0183   usesResource("TFileService");
0184   DebugMode = iConfig.getParameter<bool>("DebugMode");
0185 }
0186 
0187 /////////////
0188 // DESTRUCTOR
0189 AnalyzerClusterStub::~AnalyzerClusterStub() {
0190   /// Insert here what you need to delete
0191   /// when you close the class instance
0192 }
0193 
0194 //////////
0195 // END JOB
0196 void AnalyzerClusterStub::endJob()  //edm::Run& run, const edm::EventSetup& iSetup
0197 {
0198   /// Things to be done at the exit of the event Loop
0199   std::cerr << " AnalyzerClusterStub::endJob" << std::endl;
0200   /// End of things to be done at the exit from the event Loop
0201 }
0202 
0203 ////////////
0204 // BEGIN JOB
0205 void AnalyzerClusterStub::beginJob() {
0206   /// Initialize all slave variables
0207   /// mainly histogram ranges and resolution
0208   testedGeometry = false;
0209 
0210   std::ostringstream histoName;
0211   std::ostringstream histoTitle;
0212 
0213   /// Things to be done before entering the event Loop
0214   std::cerr << " AnalyzerClusterStub::beginJob" << std::endl;
0215 
0216   /// Book histograms etc
0217   edm::Service<TFileService> fs;
0218 
0219   /// Prepare for LogXY Plots
0220   int NumBins = 200;
0221   double MinPt = 0.0;
0222   double MaxPt = 100.0;
0223 
0224   double* BinVec = new double[NumBins + 1];
0225   for (int iBin = 0; iBin < NumBins + 1; iBin++) {
0226     double temp = pow(10, (-NumBins + iBin) / (MaxPt - MinPt));
0227     BinVec[iBin] = temp;
0228   }
0229 
0230   /// TrackingParticle and TrackingVertex
0231   hSimVtx_XY = fs->make<TH2D>("hSimVtx_XY", "SimVtx y vs. x", 200, -0.4, 0.4, 200, -0.4, 0.4);
0232   hSimVtx_RZ = fs->make<TH2D>("hSimVtx_RZ", "SimVtx #rho vs. z", 200, -50, 50, 200, 0, 0.4);
0233   hSimVtx_XY->Sumw2();
0234   hSimVtx_RZ->Sumw2();
0235 
0236   hTPart_Pt = fs->make<TH1D>("hTPart_Pt", "TPart p_{T}", 200, 0, 50);
0237   hTPart_Eta_Pt10 = fs->make<TH1D>("hTPart_Eta_Pt10", "TPart #eta (p_{T} > 10 GeV/c)", 180, -M_PI, M_PI);
0238   hTPart_Phi_Pt10 = fs->make<TH1D>("hTPart_Phi_Pt10", "TPart #phi (p_{T} > 10 GeV/c)", 180, -M_PI, M_PI);
0239   hTPart_Pt->Sumw2();
0240   hTPart_Eta_Pt10->Sumw2();
0241   hTPart_Phi_Pt10->Sumw2();
0242 
0243   /// Global position of TTCluster
0244   hCluster_Barrel_XY = fs->make<TH2D>("hCluster_Barrel_XY", "TTCluster Barrel y vs. x", 960, -120, 120, 960, -120, 120);
0245   hCluster_Barrel_XY_Zoom =
0246       fs->make<TH2D>("hCluster_Barrel_XY_Zoom", "TTCluster Barrel y vs. x", 960, 30, 60, 960, -15, 15);
0247   hCluster_Endcap_Fw_XY =
0248       fs->make<TH2D>("hCluster_Endcap_Fw_XY", "TTCluster Forward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0249   hCluster_Endcap_Bw_XY =
0250       fs->make<TH2D>("hCluster_Endcap_Bw_XY", "TTCluster Backward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0251   hCluster_RZ = fs->make<TH2D>("hCluster_RZ", "TTCluster #rho vs. z", 900, -300, 300, 480, 0, 120);
0252   hCluster_Endcap_Fw_RZ_Zoom =
0253       fs->make<TH2D>("hCluster_Endcap_Fw_RZ_Zoom", "TTCluster Forward Endcap #rho vs. z", 960, 140, 170, 960, 30, 60);
0254   hCluster_Endcap_Bw_RZ_Zoom = fs->make<TH2D>(
0255       "hCluster_Endcap_Bw_RZ_Zoom", "TTCluster Backward Endcap #rho vs. z", 960, -170, -140, 960, 70, 100);
0256   hCluster_Barrel_XY->Sumw2();
0257   hCluster_Barrel_XY_Zoom->Sumw2();
0258   hCluster_Endcap_Fw_XY->Sumw2();
0259   hCluster_Endcap_Bw_XY->Sumw2();
0260   hCluster_RZ->Sumw2();
0261   hCluster_Endcap_Fw_RZ_Zoom->Sumw2();
0262   hCluster_Endcap_Bw_RZ_Zoom->Sumw2();
0263 
0264   hCluster_IMem_Barrel = fs->make<TH1D>("hCluster_IMem_Barrel", "Inner TTCluster Stack", 12, -0.5, 11.5);
0265   hCluster_IMem_Endcap = fs->make<TH1D>("hCluster_IMem_Endcap", "Inner TTCluster Stack", 12, -0.5, 11.5);
0266   hCluster_OMem_Barrel = fs->make<TH1D>("hCluster_OMem_Barrel", "Outer TTCluster Stack", 12, -0.5, 11.5);
0267   hCluster_OMem_Endcap = fs->make<TH1D>("hCluster_OMem_Endcap", "Outer TTCluster Stack", 12, -0.5, 11.5);
0268   hCluster_IMem_Barrel->Sumw2();
0269   hCluster_IMem_Endcap->Sumw2();
0270   hCluster_OMem_Barrel->Sumw2();
0271   hCluster_OMem_Endcap->Sumw2();
0272 
0273   hCluster_Gen_Barrel = fs->make<TH1D>("hCluster_Gen_Barrel", "Genuine TTCluster Stack", 12, -0.5, 11.5);
0274   hCluster_Unkn_Barrel = fs->make<TH1D>("hCluster_Unkn_Barrel", "Unknown TTCluster Stack", 12, -0.5, 11.5);
0275   hCluster_Comb_Barrel = fs->make<TH1D>("hCluster_Comb_Barrel", "Combinatorial TTCluster Stack", 12, -0.5, 11.5);
0276   hCluster_Gen_Endcap = fs->make<TH1D>("hCluster_Gen_Endcap", "Genuine TTCluster Stack", 12, -0.5, 11.5);
0277   hCluster_Unkn_Endcap = fs->make<TH1D>("hCluster_Unkn_Endcap", "Unknown TTCluster Stack", 12, -0.5, 11.5);
0278   hCluster_Comb_Endcap = fs->make<TH1D>("hCluster_Comb_Endcap", "Combinatorial TTCluster Stack", 12, -0.5, 11.5);
0279   hCluster_Gen_Barrel->Sumw2();
0280   hCluster_Unkn_Barrel->Sumw2();
0281   hCluster_Comb_Barrel->Sumw2();
0282   hCluster_Gen_Endcap->Sumw2();
0283   hCluster_Unkn_Endcap->Sumw2();
0284   hCluster_Comb_Endcap->Sumw2();
0285 
0286   hCluster_Gen_Eta = fs->make<TH1D>("hCluster_Gen_Eta", "Genuine TTCluster #eta", 90, 0, M_PI);
0287   hCluster_Unkn_Eta = fs->make<TH1D>("hCluster_Unkn_Eta", "Unknown TTCluster #eta", 90, 0, M_PI);
0288   hCluster_Comb_Eta = fs->make<TH1D>("hCluster_Comb_Eta", "Combinatorial TTCluster #eta", 90, 0, M_PI);
0289   hCluster_Gen_Eta->Sumw2();
0290   hCluster_Unkn_Eta->Sumw2();
0291   hCluster_Comb_Eta->Sumw2();
0292 
0293   hCluster_PID = fs->make<TH2D>("hCluster_PID", "TTCluster PID (Member)", 501, -250.5, 250.5, 2, -0.5, 1.5);
0294   hCluster_W = fs->make<TH2D>("hCluster_W", "TTCluster Width (Member)", 10, -0.5, 9.5, 2, -0.5, 1.5);
0295   hCluster_PID->Sumw2();
0296   hCluster_W->Sumw2();
0297 
0298   hTPart_Eta_INormalization = fs->make<TH1D>("hTPart_Eta_INormalization", "TParticles vs. TPart #eta", 90, 0, M_PI);
0299   hTPart_Eta_ICW_1 = fs->make<TH1D>("hTPart_Eta_ICW_1", "CW 1 vs. TPart #eta", 90, 0, M_PI);
0300   hTPart_Eta_ICW_2 = fs->make<TH1D>("hTPart_Eta_ICW_2", "CW 2 vs. TPart #eta", 90, 0, M_PI);
0301   hTPart_Eta_ICW_3 = fs->make<TH1D>("hTPart_Eta_ICW_3", "CW 3 or more vs. TPart #eta", 90, 0, M_PI);
0302   hTPart_Eta_INormalization->Sumw2();
0303   hTPart_Eta_ICW_1->Sumw2();
0304   hTPart_Eta_ICW_2->Sumw2();
0305   hTPart_Eta_ICW_3->Sumw2();
0306 
0307   hTPart_Eta_ONormalization = fs->make<TH1D>("hTPart_Eta_ONormalization", "TParticles vs. TPart #eta", 90, 0, M_PI);
0308   hTPart_Eta_OCW_1 = fs->make<TH1D>("hTPart_Eta_OCW_1", "CW 1 vs. TPart #eta", 90, 0, M_PI);
0309   hTPart_Eta_OCW_2 = fs->make<TH1D>("hTPart_Eta_OCW_2", "CW 2 vs. TPart #eta", 90, 0, M_PI);
0310   hTPart_Eta_OCW_3 = fs->make<TH1D>("hTPart_Eta_OCW_3", "CW 3 or more vs. TPart #eta", 90, 0, M_PI);
0311   hTPart_Eta_ONormalization->Sumw2();
0312   hTPart_Eta_OCW_1->Sumw2();
0313   hTPart_Eta_OCW_2->Sumw2();
0314   hTPart_Eta_OCW_3->Sumw2();
0315 
0316   /// Global position of TTStub
0317   hStub_Barrel_XY = fs->make<TH2D>("hStub_Barrel_XY", "TTStub Barrel y vs. x", 960, -120, 120, 960, -120, 120);
0318   hStub_Barrel_XY_Zoom = fs->make<TH2D>("hStub_Barrel_XY_Zoom", "TTStub Barrel y vs. x", 960, 30, 60, 960, -15, 15);
0319   hStub_Endcap_Fw_XY =
0320       fs->make<TH2D>("hStub_Endcap_Fw_XY", "TTStub Forward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0321   hStub_Endcap_Bw_XY =
0322       fs->make<TH2D>("hStub_Endcap_Bw_XY", "TTStub Backward Endcap y vs. x", 960, -120, 120, 960, -120, 120);
0323   hStub_RZ = fs->make<TH2D>("hStub_RZ", "TTStub #rho vs. z", 900, -300, 300, 480, 0, 120);
0324   hStub_Endcap_Fw_RZ_Zoom =
0325       fs->make<TH2D>("hStub_Endcap_Fw_RZ_Zoom", "TTStub Forward Endcap #rho vs. z", 960, 140, 170, 960, 30, 60);
0326   hStub_Endcap_Bw_RZ_Zoom =
0327       fs->make<TH2D>("hStub_Endcap_Bw_RZ_Zoom", "TTStub Backward Endcap #rho vs. z", 960, -170, -140, 960, 70, 100);
0328   hStub_Barrel_XY->Sumw2();
0329   hStub_Barrel_XY_Zoom->Sumw2();
0330   hStub_Endcap_Fw_XY->Sumw2();
0331   hStub_Endcap_Bw_XY->Sumw2();
0332   hStub_RZ->Sumw2();
0333   hStub_Endcap_Fw_RZ_Zoom->Sumw2();
0334   hStub_Endcap_Bw_RZ_Zoom->Sumw2();
0335 
0336   hStub_Barrel = fs->make<TH1D>("hStub_Barrel", "TTStub Stack", 12, -0.5, 11.5);
0337   hStub_Endcap = fs->make<TH1D>("hStub_Endcap", "TTStub Stack", 12, -0.5, 11.5);
0338   hStub_Barrel->Sumw2();
0339   hStub_Endcap->Sumw2();
0340 
0341   hStub_Gen_Barrel = fs->make<TH1D>("hStub_Gen_Barrel", "Genuine TTStub Stack", 12, -0.5, 11.5);
0342   hStub_Unkn_Barrel = fs->make<TH1D>("hStub_Unkn_Barrel", "Unknown  TTStub Stack", 12, -0.5, 11.5);
0343   hStub_Comb_Barrel = fs->make<TH1D>("hStub_Comb_Barrel", "Combinatorial TTStub Stack", 12, -0.5, 11.5);
0344   hStub_Gen_Endcap = fs->make<TH1D>("hStub_Gen_Endcap", "Genuine TTStub Stack", 12, -0.5, 11.5);
0345   hStub_Unkn_Endcap = fs->make<TH1D>("hStub_Unkn_Endcap", "Unknown  TTStub Stack", 12, -0.5, 11.5);
0346   hStub_Comb_Endcap = fs->make<TH1D>("hStub_Comb_Endcap", "Combinatorial TTStub Stack", 12, -0.5, 11.5);
0347   hStub_Gen_Barrel->Sumw2();
0348   hStub_Unkn_Barrel->Sumw2();
0349   hStub_Comb_Barrel->Sumw2();
0350   hStub_Gen_Endcap->Sumw2();
0351   hStub_Unkn_Endcap->Sumw2();
0352   hStub_Comb_Endcap->Sumw2();
0353 
0354   hStub_Gen_Eta = fs->make<TH1D>("hStub_Gen_Eta", "Genuine TTStub #eta", 90, 0, M_PI);
0355   hStub_Unkn_Eta = fs->make<TH1D>("hStub_Unkn_Eta", "Unknown TTStub #eta", 90, 0, M_PI);
0356   hStub_Comb_Eta = fs->make<TH1D>("hStub_Comb_Eta", "Combinatorial TTStub #eta", 90, 0, M_PI);
0357   hStub_Gen_Eta->Sumw2();
0358   hStub_Unkn_Eta->Sumw2();
0359   hStub_Comb_Eta->Sumw2();
0360 
0361   hStub_PID = fs->make<TH1D>("hStub_PID", "TTStub PID", 501, -250.5, 250.5);
0362   hStub_Barrel_W =
0363       fs->make<TH2D>("hStub_Barrel_W", "TTStub Post-Corr Displacement (Layer)", 12, -0.5, 11.5, 43, -10.75, 10.75);
0364   hStub_Barrel_O = fs->make<TH2D>("hStub_Barrel_O", "TTStub Offset (Layer)", 12, -0.5, 11.5, 43, -10.75, 10.75);
0365   hStub_Endcap_W =
0366       fs->make<TH2D>("hStub_Endcap_W", "TTStub Post-Corr Displacement (Layer)", 12, -0.5, 11.5, 43, -10.75, 10.75);
0367   hStub_Endcap_O = fs->make<TH2D>("hStub_Endcap_O", "TTStub Offset (Layer)", 12, -0.5, 11.5, 43, -10.75, 10.75);
0368 
0369   hStub_PID->Sumw2();
0370   hStub_Barrel_W->Sumw2();
0371   hStub_Barrel_O->Sumw2();
0372   hStub_Endcap_W->Sumw2();
0373   hStub_Endcap_O->Sumw2();
0374 
0375   hTPart_Eta_Pt10_Normalization =
0376       fs->make<TH1D>("hTPart_Eta_Pt10_Normalization", "TParticles vs. TPart #eta", 90, 0, M_PI);
0377   hTPart_Eta_Pt10_NumPS = fs->make<TH1D>("hTPart_Eta_Pt10_NumPS", "PS Stubs vs. TPart #eta", 90, 0, M_PI);
0378   hTPart_Eta_Pt10_Num2S = fs->make<TH1D>("hTPart_Eta_Pt10_Num2S", "2S Stubs vs. TPart #eta", 90, 0, M_PI);
0379   hTPart_Eta_Pt10_Normalization->Sumw2();
0380   hTPart_Eta_Pt10_NumPS->Sumw2();
0381   hTPart_Eta_Pt10_Num2S->Sumw2();
0382 
0383   /// Stub Production Efficiency and comparison to TrackingParticle
0384   for (unsigned int stackIdx = 0; stackIdx < 12; stackIdx++) {
0385     /// BARREL
0386 
0387     /// Denominators
0388     histoName.str("");
0389     histoName << "hTPart_Pt_Clu_L" << stackIdx;
0390     histoTitle.str("");
0391     histoTitle << "TPart p_{T}, Cluster, Barrel Stack " << stackIdx;
0392     mapCluLayer_hTPart_Pt[stackIdx] = fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50);
0393     histoName.str("");
0394     histoName << "hTPart_Eta_Pt10_Clu_L" << stackIdx;
0395     histoTitle.str("");
0396     histoTitle << "TPart #eta (p_{T} > 10 GeV/c), Cluster, Barrel Stack " << stackIdx;
0397     mapCluLayer_hTPart_Eta_Pt10[stackIdx] =
0398         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0399     histoName.str("");
0400     histoName << "hTPart_Phi_Pt10_Clu_L" << stackIdx;
0401     histoTitle.str("");
0402     histoTitle << "TPart #phi (p_{T} > 10 GeV/c), Cluster, Barrel Stack " << stackIdx;
0403     mapCluLayer_hTPart_Phi_Pt10[stackIdx] =
0404         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0405     mapCluLayer_hTPart_Pt[stackIdx]->Sumw2();
0406     mapCluLayer_hTPart_Eta_Pt10[stackIdx]->Sumw2();
0407     mapCluLayer_hTPart_Phi_Pt10[stackIdx]->Sumw2();
0408 
0409     /// Numerators GeV/c
0410     histoName.str("");
0411     histoName << "hTPart_Pt_Stub_L" << stackIdx;
0412     histoTitle.str("");
0413     histoTitle << "TPart p_{T}, Stub, Barrel Stack " << stackIdx;
0414     mapStubLayer_hTPart_Pt[stackIdx] = fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50);
0415     histoName.str("");
0416     histoName << "hTPart_Eta_Pt10_Stub_L" << stackIdx;
0417     histoTitle.str("");
0418     histoTitle << "TPart #eta (p_{T} > 10 GeV/c), Stub, Barrel Stack " << stackIdx;
0419     mapStubLayer_hTPart_Eta_Pt10[stackIdx] =
0420         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0421     histoName.str("");
0422     histoName << "hTPart_Phi_Pt10_Stub_L" << stackIdx;
0423     histoTitle.str("");
0424     histoTitle << "TPart #phi (p_{T} > 10 GeV/c), Stub, Barrel Stack " << stackIdx;
0425     mapStubLayer_hTPart_Phi_Pt10[stackIdx] =
0426         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0427     mapStubLayer_hTPart_Pt[stackIdx]->Sumw2();
0428     mapStubLayer_hTPart_Eta_Pt10[stackIdx]->Sumw2();
0429     mapStubLayer_hTPart_Phi_Pt10[stackIdx]->Sumw2();
0430 
0431     /// Comparison to TrackingParticle
0432     histoName.str("");
0433     histoName << "hStub_InvPt_TPart_InvPt_L" << stackIdx;
0434     histoTitle.str("");
0435     histoTitle << "Stub p_{T}^{-1} vs. TPart p_{T}^{-1}, Barrel Stack " << stackIdx;
0436     mapStubLayer_hStub_InvPt_TPart_InvPt[stackIdx] =
0437         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0.0, 0.8, 200, 0.0, 0.8);
0438     mapStubLayer_hStub_InvPt_TPart_InvPt[stackIdx]->GetXaxis()->Set(NumBins, BinVec);
0439     mapStubLayer_hStub_InvPt_TPart_InvPt[stackIdx]->GetYaxis()->Set(NumBins, BinVec);
0440     mapStubLayer_hStub_InvPt_TPart_InvPt[stackIdx]->Sumw2();
0441 
0442     histoName.str("");
0443     histoName << "hStub_Pt_TPart_Pt_L" << stackIdx;
0444     histoTitle.str("");
0445     histoTitle << "Stub p_{T} vs. TPart p_{T}, Barrel Stack " << stackIdx;
0446     mapStubLayer_hStub_Pt_TPart_Pt[stackIdx] =
0447         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 100, 0, 50, 100, 0, 50);
0448     mapStubLayer_hStub_Pt_TPart_Pt[stackIdx]->Sumw2();
0449 
0450     histoName.str("");
0451     histoName << "hStub_Eta_TPart_Eta_L" << stackIdx;
0452     histoTitle.str("");
0453     histoTitle << "Stub #eta vs. TPart #eta, Barrel Stack " << stackIdx;
0454     mapStubLayer_hStub_Eta_TPart_Eta[stackIdx] =
0455         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 180, -M_PI, M_PI);
0456     mapStubLayer_hStub_Eta_TPart_Eta[stackIdx]->Sumw2();
0457 
0458     histoName.str("");
0459     histoName << "hStub_Phi_TPart_Phi_L" << stackIdx;
0460     histoTitle.str("");
0461     histoTitle << "Stub #phi vs. TPart #phi, Barrel Stack " << stackIdx;
0462     mapStubLayer_hStub_Phi_TPart_Phi[stackIdx] =
0463         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 180, -M_PI, M_PI);
0464     mapStubLayer_hStub_Phi_TPart_Phi[stackIdx]->Sumw2();
0465 
0466     /// Residuals
0467     histoName.str("");
0468     histoName << "hStub_InvPtRes_TPart_Eta_L" << stackIdx;
0469     histoTitle.str("");
0470     histoTitle << "Stub p_{T}^{-1} - TPart p_{T}^{-1} vs. TPart #eta, Barrel Stack " << stackIdx;
0471     mapStubLayer_hStub_InvPtRes_TPart_Eta[stackIdx] =
0472         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -2.0, 2.0);
0473     mapStubLayer_hStub_InvPtRes_TPart_Eta[stackIdx]->Sumw2();
0474 
0475     histoName.str("");
0476     histoName << "hStub_PtRes_TPart_Eta_L" << stackIdx;
0477     histoTitle.str("");
0478     histoTitle << "Stub p_{T} - TPart p_{T} vs. TPart #eta, Barrel Stack " << stackIdx;
0479     mapStubLayer_hStub_PtRes_TPart_Eta[stackIdx] =
0480         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -40, 40);
0481     mapStubLayer_hStub_PtRes_TPart_Eta[stackIdx]->Sumw2();
0482 
0483     histoName.str("");
0484     histoName << "hStub_EtaRes_TPart_Eta_L" << stackIdx;
0485     histoTitle.str("");
0486     histoTitle << "Stub #eta - TPart #eta vs. TPart #eta, Barrel Stack " << stackIdx;
0487     mapStubLayer_hStub_EtaRes_TPart_Eta[stackIdx] =
0488         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -2, 2);
0489     mapStubLayer_hStub_EtaRes_TPart_Eta[stackIdx]->Sumw2();
0490 
0491     histoName.str("");
0492     histoName << "hStub_PhiRes_TPart_Eta_L" << stackIdx;
0493     histoTitle.str("");
0494     histoTitle << "Stub #phi - TPart #phi vs. TPart #eta, Barrel Stack " << stackIdx;
0495     mapStubLayer_hStub_PhiRes_TPart_Eta[stackIdx] =
0496         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -0.5, 0.5);
0497     mapStubLayer_hStub_PhiRes_TPart_Eta[stackIdx]->Sumw2();
0498 
0499     /// Stub Width vs. Pt
0500     histoName.str("");
0501     histoName << "hStub_W_TPart_Pt_L" << stackIdx;
0502     histoTitle.str("");
0503     histoTitle << "Stub Width vs. TPart p_{T}, Barrel Stack " << stackIdx;
0504     mapStubLayer_hStub_W_TPart_Pt[stackIdx] =
0505         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50, 41, -10.25, 10.25);
0506     mapStubLayer_hStub_W_TPart_Pt[stackIdx]->Sumw2();
0507 
0508     histoName.str("");
0509     histoName << "hStub_W_TPart_InvPt_L" << stackIdx;
0510     histoTitle.str("");
0511     histoTitle << "Stub Width vs. TPart p_{T}^{-1}, Barrel Stack " << stackIdx;
0512     mapStubLayer_hStub_W_TPart_InvPt[stackIdx] =
0513         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 0.8, 41, -10.25, 10.25);
0514     mapStubLayer_hStub_W_TPart_InvPt[stackIdx]->GetXaxis()->Set(NumBins, BinVec);
0515     mapStubLayer_hStub_W_TPart_InvPt[stackIdx]->Sumw2();
0516 
0517     /// ENDCAP
0518 
0519     /// Denominators
0520     histoName.str("");
0521     histoName << "hTPart_Pt_Clu_D" << stackIdx;
0522     histoTitle.str("");
0523     histoTitle << "TPart p_{T}, Cluster, Endcap Stack " << stackIdx;
0524     mapCluDisk_hTPart_Pt[stackIdx] = fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50);
0525     histoName.str("");
0526     histoName << "hTPart_Eta_Pt10_Clu_D" << stackIdx;
0527     histoTitle.str("");
0528     histoTitle << "TPart #eta (p_{T} > 10 GeV/c), Cluster, Endcap Stack " << stackIdx;
0529     mapCluDisk_hTPart_Eta_Pt10[stackIdx] =
0530         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0531     histoName.str("");
0532     histoName << "hTPart_Phi_Pt10_Clu_D" << stackIdx;
0533     histoTitle.str("");
0534     histoTitle << "TPart #phi (p_{T} > 10 GeV/c), Cluster, Endcap Stack " << stackIdx;
0535     mapCluDisk_hTPart_Phi_Pt10[stackIdx] =
0536         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0537     mapCluDisk_hTPart_Pt[stackIdx]->Sumw2();
0538     mapCluDisk_hTPart_Eta_Pt10[stackIdx]->Sumw2();
0539     mapCluDisk_hTPart_Phi_Pt10[stackIdx]->Sumw2();
0540 
0541     /// Numerators GeV/c
0542     histoName.str("");
0543     histoName << "hTPart_Pt_Stub_D" << stackIdx;
0544     histoTitle.str("");
0545     histoTitle << "TPart p_{T}, Stub, Endcap Stack " << stackIdx;
0546     mapStubDisk_hTPart_Pt[stackIdx] = fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50);
0547     histoName.str("");
0548     histoName << "hTPart_Eta_Pt10_Stub_D" << stackIdx;
0549     histoTitle.str("");
0550     histoTitle << "TPart #eta (p_{T} > 10 GeV/c), Stub, Endcap Stack " << stackIdx;
0551     mapStubDisk_hTPart_Eta_Pt10[stackIdx] =
0552         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0553     histoName.str("");
0554     histoName << "hTPart_Phi_Pt10_Stub_D" << stackIdx;
0555     histoTitle.str("");
0556     histoTitle << "TPart #phi (p_{T} > 10 GeV/c), Stub, Endcap Stack " << stackIdx;
0557     mapStubDisk_hTPart_Phi_Pt10[stackIdx] =
0558         fs->make<TH1D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI);
0559     mapStubDisk_hTPart_Pt[stackIdx]->Sumw2();
0560     mapStubDisk_hTPart_Eta_Pt10[stackIdx]->Sumw2();
0561     mapStubDisk_hTPart_Phi_Pt10[stackIdx]->Sumw2();
0562 
0563     /// Comparison to TrackingParticle
0564     histoName.str("");
0565     histoName << "hStub_InvPt_TPart_InvPt_D" << stackIdx;
0566     histoTitle.str("");
0567     histoTitle << "Stub p_{T}^{-1} vs. TPart p_{T}^{-1}, Endcap Stack " << stackIdx;
0568     mapStubDisk_hStub_InvPt_TPart_InvPt[stackIdx] =
0569         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0.0, 0.8, 200, 0.0, 0.8);
0570     mapStubDisk_hStub_InvPt_TPart_InvPt[stackIdx]->GetXaxis()->Set(NumBins, BinVec);
0571     mapStubDisk_hStub_InvPt_TPart_InvPt[stackIdx]->GetYaxis()->Set(NumBins, BinVec);
0572     mapStubDisk_hStub_InvPt_TPart_InvPt[stackIdx]->Sumw2();
0573 
0574     histoName.str("");
0575     histoName << "hStub_Pt_TPart_Pt_D" << stackIdx;
0576     histoTitle.str("");
0577     histoTitle << "Stub p_{T} vs. TPart p_{T}, Endcap Stack " << stackIdx;
0578     mapStubDisk_hStub_Pt_TPart_Pt[stackIdx] =
0579         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 100, 0, 50, 100, 0, 50);
0580     mapStubDisk_hStub_Pt_TPart_Pt[stackIdx]->Sumw2();
0581 
0582     histoName.str("");
0583     histoName << "hStub_Eta_TPart_Eta_D" << stackIdx;
0584     histoTitle.str("");
0585     histoTitle << "Stub #eta vs. TPart #eta, Endcap Stack " << stackIdx;
0586     mapStubDisk_hStub_Eta_TPart_Eta[stackIdx] =
0587         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 180, -M_PI, M_PI);
0588     mapStubDisk_hStub_Eta_TPart_Eta[stackIdx]->Sumw2();
0589 
0590     histoName.str("");
0591     histoName << "hStub_Phi_TPart_Phi_D" << stackIdx;
0592     histoTitle.str("");
0593     histoTitle << "Stub #phi vs. TPart #phi, Endcap Stack " << stackIdx;
0594     mapStubDisk_hStub_Phi_TPart_Phi[stackIdx] =
0595         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 180, -M_PI, M_PI);
0596     mapStubDisk_hStub_Phi_TPart_Phi[stackIdx]->Sumw2();
0597 
0598     /// Residuals
0599     histoName.str("");
0600     histoName << "hStub_InvPtRes_TPart_Eta_D" << stackIdx;
0601     histoTitle.str("");
0602     histoTitle << "Stub p_{T}^{-1} - TPart p_{T}^{-1} vs. TPart #eta, Endcap Stack " << stackIdx;
0603     mapStubDisk_hStub_InvPtRes_TPart_Eta[stackIdx] =
0604         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -2.0, 2.0);
0605     mapStubDisk_hStub_InvPtRes_TPart_Eta[stackIdx]->Sumw2();
0606 
0607     histoName.str("");
0608     histoName << "hStub_PtRes_TPart_Eta_D" << stackIdx;
0609     histoTitle.str("");
0610     histoTitle << "Stub p_{T} - TPart p_{T} vs. TPart #eta, Endcap Stack " << stackIdx;
0611     mapStubDisk_hStub_PtRes_TPart_Eta[stackIdx] =
0612         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -40, 40);
0613     mapStubDisk_hStub_PtRes_TPart_Eta[stackIdx]->Sumw2();
0614 
0615     histoName.str("");
0616     histoName << "hStub_EtaRes_TPart_Eta_D" << stackIdx;
0617     histoTitle.str("");
0618     histoTitle << "Stub #eta - TPart #eta vs. TPart #eta, Endcap Stack " << stackIdx;
0619     mapStubDisk_hStub_EtaRes_TPart_Eta[stackIdx] =
0620         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -2, 2);
0621     mapStubDisk_hStub_EtaRes_TPart_Eta[stackIdx]->Sumw2();
0622 
0623     histoName.str("");
0624     histoName << "hStub_PhiRes_TPart_Eta_D" << stackIdx;
0625     histoTitle.str("");
0626     histoTitle << "Stub #phi - TPart #phi vs. TPart #eta, Endcap Stack " << stackIdx;
0627     mapStubDisk_hStub_PhiRes_TPart_Eta[stackIdx] =
0628         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 180, -M_PI, M_PI, 100, -0.5, 0.5);
0629     mapStubDisk_hStub_PhiRes_TPart_Eta[stackIdx]->Sumw2();
0630 
0631     /// Stub Width vs. Pt
0632     histoName.str("");
0633     histoName << "hStub_W_TPart_Pt_D" << stackIdx;
0634     histoTitle.str("");
0635     histoTitle << "Stub Width vs. TPart p_{T}, Endcap Stack " << stackIdx;
0636     mapStubDisk_hStub_W_TPart_Pt[stackIdx] =
0637         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 50, 41, -10.25, 10.25);
0638     mapStubDisk_hStub_W_TPart_Pt[stackIdx]->Sumw2();
0639 
0640     histoName.str("");
0641     histoName << "hStub_W_TPart_InvPt_D" << stackIdx;
0642     histoTitle.str("");
0643     histoTitle << "Stub Width vs. TPart p_{T}^{-1}, Endcap Stack " << stackIdx;
0644     mapStubDisk_hStub_W_TPart_InvPt[stackIdx] =
0645         fs->make<TH2D>(histoName.str().c_str(), histoTitle.str().c_str(), 200, 0, 0.8, 41, -10.25, 10.25);
0646     mapStubDisk_hStub_W_TPart_InvPt[stackIdx]->GetXaxis()->Set(NumBins, BinVec);
0647     mapStubDisk_hStub_W_TPart_InvPt[stackIdx]->Sumw2();
0648   }
0649 
0650   /// End of things to be done before entering the event Loop
0651 }
0652 
0653 //////////
0654 // ANALYZE
0655 void AnalyzerClusterStub::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0656   /*
0657   /// Geometry handles etc
0658   edm::ESHandle< TrackerGeometry >                GeometryHandle;
0659   edm::ESHandle< StackedTrackerGeometry >         StackedGeometryHandle;
0660   const StackedTrackerGeometry*                   theStackedGeometry;
0661   StackedTrackerGeometry::StackContainerIterator  StackedTrackerIterator;
0662 
0663   /// Geometry setup
0664   /// Set pointers to Geometry
0665   iSetup.get< TrackerDigiGeometryRecord >().get(GeometryHandle);
0666   /// Set pointers to Stacked Modules
0667   iSetup.get< StackedTrackerGeometryRecord >().get(StackedGeometryHandle);
0668   theStackedGeometry = StackedGeometryHandle.product(); /// Note this is different 
0669                                                         /// from the "global" geometry
0670 
0671   /// Magnetic Field
0672   edm::ESHandle< MagneticField > magneticFieldHandle;
0673   iSetup.get< IdealMagneticFieldRecord >().get(magneticFieldHandle);
0674   const MagneticField* theMagneticField = magneticFieldHandle.product();
0675   double mMagneticFieldStrength = theMagneticField->inTesla(GlobalPoint(0,0,0)).z();
0676 
0677   /// TrackingParticles
0678   edm::Handle< std::vector< TrackingParticle > > TrackingParticleHandle;
0679   iEvent.getByLabel( "mix", "MergedTrackTruth", TrackingParticleHandle );
0680   edm::Handle< std::vector< TrackingVertex > > TrackingVertexHandle;
0681   iEvent.getByLabel( "mix", "MergedTrackTruth", TrackingVertexHandle );
0682 
0683   /// Track Trigger
0684   edm::Handle< edmNew::DetSetVector< TTCluster< Ref_Phase2TrackerDigi_ > > > Phase2TrackerDigiTTClusterHandle;
0685   edm::Handle< edmNew::DetSetVector< TTStub< Ref_Phase2TrackerDigi_ > > >    Phase2TrackerDigiTTStubHandle;
0686   /// NOTE: the InputTag for the "Accepted" clusters is different from the "Inclusive" one
0687   /// "TTClustersFromPhase2TrackerDigis", "ClusterInclusive" BUT ...
0688   /// ... "TTStubsFromPhase2TrackerDigis", "ClusterAccepted"
0689   iEvent.getByLabel( "TTClustersFromPhase2TrackerDigis", "ClusterInclusive", Phase2TrackerDigiTTClusterHandle );
0690   iEvent.getByLabel( "TTStubsFromPhase2TrackerDigis", "StubAccepted",        Phase2TrackerDigiTTStubHandle );
0691 
0692   /// Track Trigger MC Truth
0693   edm::Handle< TTClusterAssociationMap< Ref_Phase2TrackerDigi_ > > MCTruthTTClusterHandle;
0694   edm::Handle< TTStubAssociationMap< Ref_Phase2TrackerDigi_ > >    MCTruthTTStubHandle;
0695   iEvent.getByLabel( "TTClusterAssociatorFromPhase2TrackerDigis", "ClusterInclusive", MCTruthTTClusterHandle );
0696   iEvent.getByLabel( "TTStubAssociatorFromPhase2TrackerDigis", "StubAccepted",        MCTruthTTStubHandle );
0697 
0698   ////////////////////////////////
0699   /// COLLECT STUB INFORMATION ///
0700   ////////////////////////////////
0701 
0702   /// Eta coverage
0703   /// Go on only if there are TrackingParticles
0704   if ( TrackingParticleHandle->size() > 0 )
0705   {
0706     /// Loop over TrackingParticles
0707     unsigned int tpCnt = 0;
0708     std::vector< TrackingParticle >::const_iterator iterTP;
0709     for ( iterTP = TrackingParticleHandle->begin();
0710           iterTP != TrackingParticleHandle->end();
0711           ++iterTP )
0712     {
0713       /// Make the pointer
0714       edm::Ptr< TrackingParticle > tempTPPtr( TrackingParticleHandle, tpCnt++ );
0715 
0716       /// Search the cluster MC map
0717       std::vector< edm::Ref< edmNew::DetSetVector< TTCluster< Ref_Phase2TrackerDigi_ > >, TTCluster< Ref_Phase2TrackerDigi_ > > > theseClusters = MCTruthTTClusterHandle->findTTClusterRefs( tempTPPtr );
0718 
0719       if ( theseClusters.size() > 0 )
0720       {
0721         bool normIClu = false;
0722         bool normOClu = false;
0723 
0724         /// Loop over the Clusters
0725         for ( unsigned int jc = 0; jc < theseClusters.size(); jc++ )
0726         {
0727           /// Check if it is good
0728           bool genuineClu = MCTruthTTClusterHandle->isGenuine( theseClusters.at(jc) );
0729           if ( !genuineClu )
0730             continue;
0731 
0732           unsigned int stackMember = theseClusters.at(jc)->getStackMember();
0733           unsigned int clusterWidth = theseClusters.at(jc)->findWidth();
0734 
0735           if ( stackMember == 0 )
0736           {
0737             if ( normIClu == false )
0738             {
0739               hTPart_Eta_INormalization->Fill( fabs( tempTPPtr->momentum().eta() ) );
0740               normIClu = true;
0741             }
0742 
0743             if ( clusterWidth == 1 )
0744             {
0745               hTPart_Eta_ICW_1->Fill( fabs( tempTPPtr->momentum().eta() ) );
0746             }
0747             else if ( clusterWidth == 2 )
0748             {
0749               hTPart_Eta_ICW_2->Fill( fabs( tempTPPtr->momentum().eta() ) );
0750             }
0751             else
0752             {
0753               hTPart_Eta_ICW_3->Fill( fabs( tempTPPtr->momentum().eta() ) );
0754             }
0755           }
0756           else if ( stackMember == 1 )
0757           {
0758             if ( normOClu == false )
0759             {
0760               hTPart_Eta_ONormalization->Fill( fabs( tempTPPtr->momentum().eta() ) );
0761               normOClu = true;
0762             }
0763 
0764             if ( clusterWidth == 1 )
0765             {
0766               hTPart_Eta_OCW_1->Fill( fabs( tempTPPtr->momentum().eta() ) );
0767             }
0768             else if ( clusterWidth == 2 )
0769             {
0770               hTPart_Eta_OCW_2->Fill( fabs( tempTPPtr->momentum().eta() ) );
0771             }
0772             else
0773             {
0774               hTPart_Eta_OCW_3->Fill( fabs( tempTPPtr->momentum().eta() ) );
0775             }
0776           }
0777         } /// End of loop over clusters
0778       }
0779 
0780       /// Search the stub MC truth map
0781       std::vector< edm::Ref< edmNew::DetSetVector< TTStub< Ref_Phase2TrackerDigi_ > >, TTStub< Ref_Phase2TrackerDigi_ > > > theseStubs = MCTruthTTStubHandle->findTTStubRefs( tempTPPtr );
0782 
0783       if ( tempTPPtr->p4().pt() <= 10 )
0784         continue; 
0785 
0786       if ( theseStubs.size() > 0 )
0787       {
0788         bool normStub = false;
0789 
0790         /// Loop over the Stubs
0791         for ( unsigned int js = 0; js < theseStubs.size(); js++ )
0792         {
0793           /// Check if it is good
0794           bool genuineStub = MCTruthTTStubHandle->isGenuine( theseStubs.at(js) );
0795           if ( !genuineStub )
0796             continue;
0797 
0798           if ( normStub == false )
0799           {
0800             hTPart_Eta_Pt10_Normalization->Fill( fabs( tempTPPtr->momentum().eta() ) );
0801             normStub = true;
0802           }
0803 
0804           /// Classify the stub
0805           StackedTrackerDetId stDetId( theseStubs.at(js)->getDetId() );
0806           /// Check if there are PS modules in seed or candidate
0807           const GeomDetUnit* det0 = theStackedGeometry->idToDetUnit( stDetId, 0 );
0808           const GeomDetUnit* det1 = theStackedGeometry->idToDetUnit( stDetId, 1 );
0809           /// Find pixel pitch and topology related information
0810           const PixelGeomDetUnit* pix0 = dynamic_cast< const PixelGeomDetUnit* >( det0 );
0811           const PixelGeomDetUnit* pix1 = dynamic_cast< const PixelGeomDetUnit* >( det1 );
0812           const PixelTopology* top0 = dynamic_cast< const PixelTopology* >( &(pix0->specificTopology()) );
0813           const PixelTopology* top1 = dynamic_cast< const PixelTopology* >( &(pix1->specificTopology()) );
0814           int cols0 = top0->ncolumns();
0815           int cols1 = top1->ncolumns();
0816           int ratio = cols0/cols1; /// This assumes the ratio is integer!
0817 
0818           if ( ratio == 1 ) /// 2S Modules
0819           {
0820             hTPart_Eta_Pt10_Num2S->Fill( fabs( tempTPPtr->momentum().eta() ) );
0821           }
0822           else /// PS
0823           {
0824             hTPart_Eta_Pt10_NumPS->Fill( fabs( tempTPPtr->momentum().eta() ) );
0825           }
0826         } /// End of loop over the Stubs generated by this TrackingParticle
0827       }
0828     } /// End of loop over TrackingParticles
0829   }
0830 
0831   /// Maps to store TrackingParticle information
0832   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > > tpPerLayer;
0833   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > > tpPerDisk;
0834 
0835   /// Loop over the input Clusters
0836   typename edmNew::DetSetVector< TTCluster< Ref_Phase2TrackerDigi_ > >::const_iterator inputIter;
0837   typename edmNew::DetSet< TTCluster< Ref_Phase2TrackerDigi_ > >::const_iterator contentIter;
0838   for ( inputIter = Phase2TrackerDigiTTClusterHandle->begin();
0839         inputIter != Phase2TrackerDigiTTClusterHandle->end();
0840         ++inputIter )
0841   {
0842     for ( contentIter = inputIter->begin();
0843           contentIter != inputIter->end();
0844           ++contentIter )
0845     {
0846       /// Make the reference to be put in the map
0847       edm::Ref< edmNew::DetSetVector< TTCluster< Ref_Phase2TrackerDigi_ > >, TTCluster< Ref_Phase2TrackerDigi_ > > tempCluRef = edmNew::makeRefTo( Phase2TrackerDigiTTClusterHandle, contentIter );
0848 
0849       StackedTrackerDetId detIdClu( tempCluRef->getDetId() );
0850       unsigned int memberClu = tempCluRef->getStackMember();
0851       bool genuineClu     = MCTruthTTClusterHandle->isGenuine( tempCluRef );
0852       bool combinClu      = MCTruthTTClusterHandle->isCombinatoric( tempCluRef );
0853       //bool unknownClu     = MCTruthTTClusterHandle->isUnknown( tempCluRef );
0854       int partClu         = 999999999;
0855       if ( genuineClu )
0856       {
0857         edm::Ptr< TrackingParticle > thisTP = MCTruthTTClusterHandle->findTrackingParticlePtr( tempCluRef );
0858         partClu = thisTP->pdgId();
0859       }
0860       unsigned int widClu = tempCluRef->findWidth();
0861       GlobalPoint posClu  = theStackedGeometry->findAverageGlobalPosition( &(*tempCluRef) );
0862       
0863       hCluster_RZ->Fill( posClu.z(), posClu.perp() );
0864 
0865       if ( detIdClu.isBarrel() )
0866       {
0867         if ( memberClu == 0 )
0868         {
0869           hCluster_IMem_Barrel->Fill( detIdClu.iLayer() );
0870         }
0871         else
0872         {
0873           hCluster_OMem_Barrel->Fill( detIdClu.iLayer() );
0874         }
0875 
0876         if ( genuineClu )
0877         {
0878           hCluster_Gen_Barrel->Fill( detIdClu.iLayer() );
0879         }
0880         else if ( combinClu )
0881         {
0882           hCluster_Comb_Barrel->Fill( detIdClu.iLayer() );
0883         }
0884         else
0885         {
0886           hCluster_Unkn_Barrel->Fill( detIdClu.iLayer() );
0887         }
0888 
0889         hCluster_Barrel_XY->Fill( posClu.x(), posClu.y() );
0890         hCluster_Barrel_XY_Zoom->Fill( posClu.x(), posClu.y() );
0891       }
0892       else if ( detIdClu.isEndcap() )
0893       {
0894         if ( memberClu == 0 )
0895         {
0896           hCluster_IMem_Endcap->Fill( detIdClu.iDisk() );
0897         }
0898         else
0899         {
0900           hCluster_OMem_Endcap->Fill( detIdClu.iDisk() );
0901         }
0902 
0903         if ( genuineClu )
0904         {
0905           hCluster_Gen_Endcap->Fill( detIdClu.iDisk() );
0906         }
0907         else if ( combinClu )
0908         {
0909           hCluster_Comb_Endcap->Fill( detIdClu.iDisk() );
0910         }
0911         else
0912         {
0913           hCluster_Unkn_Endcap->Fill( detIdClu.iDisk() );
0914         }
0915 
0916         if ( posClu.z() > 0 )
0917         {
0918           hCluster_Endcap_Fw_XY->Fill( posClu.x(), posClu.y() );
0919           hCluster_Endcap_Fw_RZ_Zoom->Fill( posClu.z(), posClu.perp() );
0920         }
0921         else
0922         {
0923           hCluster_Endcap_Bw_XY->Fill( posClu.x(), posClu.y() );
0924           hCluster_Endcap_Bw_RZ_Zoom->Fill( posClu.z(), posClu.perp() );
0925         }
0926       }
0927 
0928       /// Another way of looking at MC truth
0929       if ( genuineClu )
0930       {
0931         hCluster_Gen_Eta->Fill( fabs(posClu.eta()) );
0932       }
0933       else if ( combinClu )
0934       {
0935         hCluster_Comb_Eta->Fill( fabs(posClu.eta()) );
0936       }
0937       else
0938       {
0939         hCluster_Unkn_Eta->Fill( fabs(posClu.eta()) );
0940       }
0941 
0942       hCluster_PID->Fill( partClu, memberClu );
0943       hCluster_W->Fill( widClu, memberClu );
0944 
0945       /// Store Track information in maps, skip if the Cluster is not good
0946       if ( !genuineClu && !combinClu ) continue;
0947 
0948       std::vector< edm::Ptr< TrackingParticle > > theseTPs = MCTruthTTClusterHandle->findTrackingParticlePtrs( tempCluRef );
0949 
0950       for ( unsigned int i = 0; i < theseTPs.size(); i++ )
0951       {
0952         edm::Ptr< TrackingParticle > tpPtr = theseTPs.at(i);
0953 
0954         if ( tpPtr.isNull() )
0955           continue;
0956 
0957         /// Get the corresponding vertex and reject the track
0958         /// if its vertex is outside the beampipe
0959         if ( tpPtr->vertex().rho() >= 2.0 )
0960           continue;
0961 
0962         if ( detIdClu.isBarrel() )
0963         {
0964           if ( tpPerLayer.find( detIdClu.iLayer() ) == tpPerLayer.end() )
0965           {
0966             std::vector< edm::Ptr< TrackingParticle > > tempVec;
0967             tpPerLayer.insert( make_pair( detIdClu.iLayer(), tempVec ) );
0968           }
0969           tpPerLayer[detIdClu.iLayer()].push_back( tpPtr );
0970         }
0971         else if ( detIdClu.isEndcap() )
0972         {
0973           if ( tpPerDisk.find( detIdClu.iDisk() ) == tpPerDisk.end() )
0974           {
0975             std::vector< edm::Ptr< TrackingParticle > > tempVec;
0976             tpPerDisk.insert( make_pair( detIdClu.iDisk(), tempVec ) );
0977           }
0978           tpPerDisk[detIdClu.iDisk()].push_back( tpPtr );
0979         }
0980       }
0981     }
0982   } /// End of Loop over TTClusters
0983 
0984   /// Clean the maps for TrackingParticles and fill histograms
0985   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > >::iterator iterTPPerLayer;
0986   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > >::iterator iterTPPerDisk;
0987 
0988   for ( iterTPPerLayer = tpPerLayer.begin();
0989         iterTPPerLayer != tpPerLayer.end();
0990         ++iterTPPerLayer )
0991   {
0992     /// Remove duplicates, if any
0993     std::vector< edm::Ptr< TrackingParticle > > tempVec = iterTPPerLayer->second;
0994     std::sort( tempVec.begin(), tempVec.end() );
0995     tempVec.erase( std::unique( tempVec.begin(), tempVec.end() ), tempVec.end() );
0996 
0997     /// Loop over the TrackingParticles in this piece of the map
0998     for ( unsigned int i = 0; i < tempVec.size(); i++ )
0999     {
1000       if ( tempVec.at(i).isNull() ) continue;
1001       TrackingParticle thisTP = *(tempVec.at(i));
1002       mapCluLayer_hTPart_Pt[ iterTPPerLayer->first ]->Fill( thisTP.p4().pt() );
1003       if ( thisTP.p4().pt() > 10.0 )
1004       {
1005         mapCluLayer_hTPart_Eta_Pt10[ iterTPPerLayer->first ]->Fill( thisTP.momentum().eta() );
1006         mapCluLayer_hTPart_Phi_Pt10[ iterTPPerLayer->first ]->Fill( thisTP.momentum().phi() > M_PI ?
1007                                                                      thisTP.momentum().phi() - 2*M_PI :
1008                                                                      thisTP.momentum().phi() );    
1009       }
1010     }
1011   }
1012 
1013   for ( iterTPPerDisk = tpPerDisk.begin();
1014         iterTPPerDisk != tpPerDisk.end();
1015         ++iterTPPerDisk )
1016   {
1017     /// Remove duplicates, if any
1018     std::vector< edm::Ptr< TrackingParticle > > tempVec = iterTPPerDisk->second;
1019     std::sort( tempVec.begin(), tempVec.end() );
1020     tempVec.erase( std::unique( tempVec.begin(), tempVec.end() ), tempVec.end() );
1021 
1022     /// Loop over the TrackingParticles in this piece of the map
1023     for ( unsigned int i = 0; i < tempVec.size(); i++ )
1024     {
1025       if ( tempVec.at(i).isNull() ) continue;
1026       TrackingParticle thisTP = *(tempVec.at(i));
1027       mapCluDisk_hTPart_Pt[ iterTPPerDisk->first ]->Fill( thisTP.p4().pt() );
1028       if ( thisTP.p4().pt() > 10.0 )
1029       {
1030         mapCluDisk_hTPart_Eta_Pt10[ iterTPPerDisk->first ]->Fill( thisTP.momentum().eta() );
1031         mapCluDisk_hTPart_Phi_Pt10[ iterTPPerDisk->first ]->Fill( thisTP.momentum().phi() > M_PI ?
1032                                                                   thisTP.momentum().phi() - 2*M_PI :
1033                                                                   thisTP.momentum().phi() );    
1034       }
1035     }
1036   }
1037 
1038   ////////////////////////////////
1039   /// COLLECT STUB INFORMATION ///
1040   ////////////////////////////////
1041 
1042   /// Maps to store TrackingParticle information
1043   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > > tpPerStubLayer;
1044   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > > tpPerStubDisk;
1045 
1046   /// Loop over the input Stubs
1047   typename edmNew::DetSetVector< TTStub< Ref_Phase2TrackerDigi_ > >::const_iterator otherInputIter;
1048   typename edmNew::DetSet< TTStub< Ref_Phase2TrackerDigi_ > >::const_iterator otherContentIter;
1049   for ( otherInputIter = Phase2TrackerDigiTTStubHandle->begin();
1050         otherInputIter != Phase2TrackerDigiTTStubHandle->end();
1051         ++otherInputIter )
1052   {
1053     for ( otherContentIter = otherInputIter->begin();
1054           otherContentIter != otherInputIter->end();
1055           ++otherContentIter )
1056     {
1057       /// Make the reference to be put in the map
1058       edm::Ref< edmNew::DetSetVector< TTStub< Ref_Phase2TrackerDigi_ > >, TTStub< Ref_Phase2TrackerDigi_ > > tempStubRef = edmNew::makeRefTo( Phase2TrackerDigiTTStubHandle, otherContentIter );
1059 
1060       StackedTrackerDetId detIdStub( tempStubRef->getDetId() );
1061 
1062       bool genuineStub    = MCTruthTTStubHandle->isGenuine( tempStubRef );
1063       bool combinStub     = MCTruthTTStubHandle->isCombinatoric( tempStubRef );
1064       //bool unknownStub    = MCTruthTTStubHandle->isUnknown( tempStubRef );
1065       int partStub         = 999999999;
1066       if ( genuineStub )
1067       {
1068         edm::Ptr< TrackingParticle > thisTP = MCTruthTTStubHandle->findTrackingParticlePtr( tempStubRef );
1069         partStub = thisTP->pdgId();
1070       }
1071       double displStub    = tempStubRef->getRawBend();
1072       double offsetStub   = tempStubRef->getBendOffset();
1073       GlobalPoint posStub = theStackedGeometry->findGlobalPosition( &(*tempStubRef) );
1074 
1075       hStub_RZ->Fill( posStub.z(), posStub.perp() );
1076 
1077       if ( detIdStub.isBarrel() )
1078       {
1079         hStub_Barrel->Fill( detIdStub.iLayer() );
1080 
1081         if ( genuineStub )
1082         {
1083           hStub_Gen_Barrel->Fill( detIdStub.iLayer() );
1084         }
1085         else if ( combinStub )
1086         {
1087           hStub_Comb_Barrel->Fill( detIdStub.iLayer() );
1088         }
1089         else
1090         {
1091           hStub_Unkn_Barrel->Fill( detIdStub.iLayer() );
1092         }
1093 
1094         hStub_Barrel_XY->Fill( posStub.x(), posStub.y() );
1095         hStub_Barrel_XY_Zoom->Fill( posStub.x(), posStub.y() );
1096       }
1097       else if ( detIdStub.isEndcap() )
1098       {
1099         hStub_Endcap->Fill( detIdStub.iDisk() );
1100 
1101         if ( genuineStub )
1102         {
1103           hStub_Gen_Endcap->Fill( detIdStub.iDisk() );
1104         }
1105         else if ( combinStub )
1106         {
1107           hStub_Comb_Endcap->Fill( detIdStub.iDisk() );
1108         }
1109         else
1110         {
1111           hStub_Unkn_Endcap->Fill( detIdStub.iDisk() );
1112         }
1113 
1114         if ( posStub.z() > 0 ) 
1115         {
1116           hStub_Endcap_Fw_XY->Fill( posStub.x(), posStub.y() );
1117           hStub_Endcap_Fw_RZ_Zoom->Fill( posStub.z(), posStub.perp() );
1118         }
1119         else
1120         {
1121           hStub_Endcap_Bw_XY->Fill( posStub.x(), posStub.y() );
1122           hStub_Endcap_Bw_RZ_Zoom->Fill( posStub.z(), posStub.perp() );
1123         }
1124       }
1125 
1126       /// Another way of looking at MC truth
1127       if ( genuineStub )
1128       {
1129         hStub_Gen_Eta->Fill( fabs(posStub.eta()) );
1130       }
1131       else if ( combinStub )
1132       {
1133         hStub_Comb_Eta->Fill( fabs(posStub.eta()) );
1134       }
1135       else
1136       {
1137         hStub_Unkn_Eta->Fill( fabs(posStub.eta()) );
1138       }
1139 
1140       hStub_PID->Fill( partStub );
1141 
1142       /// Store Track information in maps, skip if the Cluster is not good
1143       if ( !genuineStub ) continue;
1144 
1145       edm::Ptr< TrackingParticle > tpPtr = MCTruthTTStubHandle->findTrackingParticlePtr( tempStubRef );
1146 
1147       /// Get the corresponding vertex and reject the track
1148       /// if its vertex is outside the beampipe
1149       if ( tpPtr->vertex().rho() >= 2.0 )
1150         continue;
1151 
1152       if ( detIdStub.isBarrel() )
1153       {
1154         if ( tpPerStubLayer.find( detIdStub.iLayer() ) == tpPerStubLayer.end() )
1155         {
1156           std::vector< edm::Ptr< TrackingParticle > > tempVec;
1157           tpPerStubLayer.insert( make_pair( detIdStub.iLayer(), tempVec ) );
1158         }
1159         tpPerStubLayer[detIdStub.iLayer()].push_back( tpPtr );
1160 
1161         hStub_Barrel_W->Fill( detIdStub.iLayer(), displStub - offsetStub );
1162         hStub_Barrel_O->Fill( detIdStub.iLayer(), offsetStub );
1163       }
1164       else if ( detIdStub.isEndcap() )
1165       {
1166         if ( tpPerStubDisk.find( detIdStub.iDisk() ) == tpPerStubDisk.end() )
1167         {
1168           std::vector< edm::Ptr< TrackingParticle > > tempVec;
1169           tpPerStubDisk.insert( make_pair( detIdStub.iDisk(), tempVec ) );
1170         }
1171         tpPerStubDisk[detIdStub.iDisk()].push_back( tpPtr );
1172 
1173         hStub_Endcap_W->Fill( detIdStub.iDisk(), displStub - offsetStub );
1174         hStub_Endcap_O->Fill( detIdStub.iDisk(), offsetStub );
1175       }
1176       
1177       /// Compare to TrackingParticle
1178 
1179       if ( tpPtr.isNull() ) continue; /// This prevents to fill the vector if the TrackingParticle is not found
1180       TrackingParticle thisTP = *tpPtr;
1181 
1182       double simPt = thisTP.p4().pt();
1183       double simEta = thisTP.momentum().eta();
1184       double simPhi = thisTP.momentum().phi();
1185       double recPt = theStackedGeometry->findRoughPt( mMagneticFieldStrength, &(*tempStubRef) );
1186       double recEta = theStackedGeometry->findGlobalDirection( &(*tempStubRef) ).eta();
1187       double recPhi = theStackedGeometry->findGlobalDirection( &(*tempStubRef) ).phi();
1188 
1189       if ( simPhi > M_PI )
1190       {
1191         simPhi -= 2*M_PI;
1192       }
1193       if ( recPhi > M_PI )
1194       {
1195         recPhi -= 2*M_PI;
1196       }
1197 
1198       if ( detIdStub.isBarrel() )
1199       {
1200         mapStubLayer_hStub_InvPt_TPart_InvPt[ detIdStub.iLayer() ]->Fill( 1./simPt, 1./recPt );
1201         mapStubLayer_hStub_Pt_TPart_Pt[ detIdStub.iLayer() ]->Fill( simPt, recPt );
1202         mapStubLayer_hStub_Eta_TPart_Eta[ detIdStub.iLayer() ]->Fill( simEta, recEta );
1203         mapStubLayer_hStub_Phi_TPart_Phi[ detIdStub.iLayer() ]->Fill( simPhi, recPhi );
1204 
1205         mapStubLayer_hStub_InvPtRes_TPart_Eta[ detIdStub.iLayer() ]->Fill( simEta, 1./recPt - 1./simPt );
1206         mapStubLayer_hStub_PtRes_TPart_Eta[ detIdStub.iLayer() ]->Fill( simEta, recPt - simPt );
1207         mapStubLayer_hStub_EtaRes_TPart_Eta[ detIdStub.iLayer() ]->Fill( simEta, recEta - simEta );
1208         mapStubLayer_hStub_PhiRes_TPart_Eta[ detIdStub.iLayer() ]->Fill( simEta, recPhi - simPhi );
1209 
1210         mapStubLayer_hStub_W_TPart_Pt[ detIdStub.iLayer() ]->Fill( simPt, displStub - offsetStub );
1211         mapStubLayer_hStub_W_TPart_InvPt[ detIdStub.iLayer() ]->Fill( 1./simPt, displStub - offsetStub );
1212       }
1213       else if ( detIdStub.isEndcap() )
1214       {
1215         mapStubDisk_hStub_InvPt_TPart_InvPt[ detIdStub.iDisk() ]->Fill( 1./simPt, 1./recPt );
1216         mapStubDisk_hStub_Pt_TPart_Pt[ detIdStub.iDisk() ]->Fill( simPt, recPt );
1217         mapStubDisk_hStub_Eta_TPart_Eta[ detIdStub.iDisk() ]->Fill( simEta, recEta );
1218         mapStubDisk_hStub_Phi_TPart_Phi[ detIdStub.iDisk() ]->Fill( simPhi, recPhi );
1219 
1220         mapStubDisk_hStub_InvPtRes_TPart_Eta[ detIdStub.iDisk() ]->Fill( simEta, 1./recPt - 1./simPt );
1221         mapStubDisk_hStub_PtRes_TPart_Eta[ detIdStub.iDisk() ]->Fill( simEta, recPt - simPt );
1222         mapStubDisk_hStub_EtaRes_TPart_Eta[ detIdStub.iDisk() ]->Fill( simEta, recEta - simEta );
1223         mapStubDisk_hStub_PhiRes_TPart_Eta[ detIdStub.iDisk() ]->Fill( simEta, recPhi - simPhi );
1224 
1225         mapStubDisk_hStub_W_TPart_Pt[ detIdStub.iDisk() ]->Fill( simPt, displStub - offsetStub );
1226         mapStubDisk_hStub_W_TPart_InvPt[ detIdStub.iDisk() ]->Fill( 1./simPt, displStub - offsetStub );
1227       }
1228     }
1229   } /// End of loop over TTStubs
1230 
1231   /// Clean the maps for TrackingParticles and fill histograms
1232   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > >::iterator iterTPPerStubLayer;
1233   std::map< unsigned int, std::vector< edm::Ptr< TrackingParticle > > >::iterator iterTPPerStubDisk;
1234 
1235   for ( iterTPPerStubLayer = tpPerStubLayer.begin();
1236         iterTPPerStubLayer != tpPerStubLayer.end();
1237         ++iterTPPerStubLayer ) 
1238   {
1239     /// Remove duplicates, if any
1240     std::vector< edm::Ptr< TrackingParticle > > tempVec = iterTPPerStubLayer->second;
1241     std::sort( tempVec.begin(), tempVec.end() );
1242     tempVec.erase( std::unique( tempVec.begin(), tempVec.end() ), tempVec.end() );
1243 
1244     /// Loop over the TrackingParticles in this piece of the map
1245     for ( unsigned int i = 0; i < tempVec.size(); i++ )
1246     {
1247       if ( tempVec.at(i).isNull() ) continue;
1248       TrackingParticle thisTP = *(tempVec.at(i));
1249       mapStubLayer_hTPart_Pt[ iterTPPerStubLayer->first ]->Fill( thisTP.p4().pt() );
1250       if ( thisTP.p4().pt() > 10.0 )
1251       {
1252         mapStubLayer_hTPart_Eta_Pt10[ iterTPPerStubLayer->first ]->Fill( thisTP.momentum().eta() );
1253         mapStubLayer_hTPart_Phi_Pt10[ iterTPPerStubLayer->first ]->Fill( thisTP.momentum().phi() > M_PI ?
1254                                                                           thisTP.momentum().phi() - 2*M_PI :
1255                                                                           thisTP.momentum().phi() );    
1256       }
1257     }
1258   }
1259 
1260   for ( iterTPPerStubDisk = tpPerStubDisk.begin();
1261         iterTPPerStubDisk != tpPerStubDisk.end();
1262         ++iterTPPerStubDisk ) 
1263   {
1264     /// Remove duplicates, if any
1265     std::vector< edm::Ptr< TrackingParticle > > tempVec = iterTPPerStubDisk->second;
1266     std::sort( tempVec.begin(), tempVec.end() );
1267     tempVec.erase( std::unique( tempVec.begin(), tempVec.end() ), tempVec.end() );
1268 
1269     /// Loop over the TrackingParticles in this piece of the map
1270     for ( unsigned int i = 0; i < tempVec.size(); i++ )
1271     {
1272       if ( tempVec.at(i).isNull() ) continue;
1273       TrackingParticle thisTP = *(tempVec.at(i));
1274       mapStubDisk_hTPart_Pt[ iterTPPerStubDisk->first ]->Fill( thisTP.p4().pt() );
1275       if ( thisTP.p4().pt() > 10.0 )
1276       {
1277         mapStubDisk_hTPart_Eta_Pt10[ iterTPPerStubDisk->first ]->Fill( thisTP.momentum().eta() );
1278         mapStubDisk_hTPart_Phi_Pt10[ iterTPPerStubDisk->first ]->Fill( thisTP.momentum().phi() > M_PI ?
1279                                                                         thisTP.momentum().phi() - 2*M_PI :
1280                                                                         thisTP.momentum().phi() );    
1281       }
1282     }
1283   }
1284 
1285   /// //////////////////////////
1286   /// SPECTRUM OF SIM TRACKS ///
1287   /// WITHIN PRIMARY VERTEX  ///
1288   /// CONSTRAINTS            ///
1289   /// //////////////////////////
1290 
1291   /// Go on only if there are TrackingParticles
1292   if ( TrackingParticleHandle->size() != 0 )
1293   {
1294     /// Loop over TrackingParticles
1295     std::vector< TrackingParticle >::const_iterator iterTrackingParticles;
1296     for ( iterTrackingParticles = TrackingParticleHandle->begin();
1297           iterTrackingParticles != TrackingParticleHandle->end();
1298           ++iterTrackingParticles )
1299     {
1300       /// Get the corresponding vertex
1301       /// Assume perfectly round beamspot
1302       /// Correct and get the correct TrackingParticle Vertex position wrt beam center
1303       if ( iterTrackingParticles->vertex().rho() >= 2 )
1304         continue;
1305 
1306       /// First of all, check beamspot and correction
1307       hSimVtx_XY->Fill( iterTrackingParticles->vertex().x(), iterTrackingParticles->vertex().y() );
1308       hSimVtx_RZ->Fill( iterTrackingParticles->vertex().z(), iterTrackingParticles->vertex().rho() );
1309 
1310       /// Here we have only tracks form primary vertices
1311       /// Check Pt spectrum and pseudorapidity for over-threshold tracks
1312       hTPart_Pt->Fill( iterTrackingParticles->p4().pt() );
1313       if ( iterTrackingParticles->p4().pt() > 10.0 )
1314       {
1315         hTPart_Eta_Pt10->Fill( iterTrackingParticles->momentum().eta() );
1316         hTPart_Phi_Pt10->Fill( iterTrackingParticles->momentum().phi() > M_PI ?
1317                                 iterTrackingParticles->momentum().phi() - 2*M_PI :
1318                                 iterTrackingParticles->momentum().phi() );
1319       }
1320     } /// End of Loop over TrackingParticles
1321   } /// End of if ( TrackingParticleHandle->size() != 0 )
1322 */
1323 }  /// End of analyze()
1324 
1325 ///////////////////////////
1326 // DEFINE THIS AS A PLUG-IN
1327 DEFINE_FWK_MODULE(AnalyzerClusterStub);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!