Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:11 

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