Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:31:02

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