Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:25:27

0001 // myFastSimVal.cc
0002 // Description:  Comparison between Jet Algorithms
0003 // Author: Frank Chlebana
0004 // Date:  08 - August - 2007
0005 //
0006 #include "RecoJets/JetAnalyzers/interface/myFastSimVal.h"
0007 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
0008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0009 #include "DataFormats/JetReco/interface/CaloJet.h"
0010 #include "DataFormats/JetReco/interface/GenJet.h"
0011 #include "DataFormats/Math/interface/deltaR.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013 // #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
0014 #include "DataFormats/Candidate/interface/Candidate.h"
0015 // #include "FWCore/Framework/interface/Handle.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include <TROOT.h>
0019 #include <TSystem.h>
0020 #include <TFile.h>
0021 #include <TCanvas.h>
0022 #include <cmath>
0023 using namespace edm;
0024 using namespace reco;
0025 using namespace std;
0026 
0027 #define DEBUG 1
0028 // #define MAXJETS 50
0029 #define MAXJETS 100
0030 
0031 // Get the algorithm of the jet collections we will read from the .cfg file
0032 // which defines the value of the strings CaloJetAlgorithm and GenJetAlgorithm.
0033 myFastSimVal::myFastSimVal(const ParameterSet& cfg)
0034     : CaloJetAlgorithm1(cfg.getParameter<string>("CaloJetAlgorithm1")),
0035       CaloJetAlgorithm2(cfg.getParameter<string>("CaloJetAlgorithm2")),
0036       CaloJetAlgorithm3(cfg.getParameter<string>("CaloJetAlgorithm3")),
0037       CaloJetAlgorithm4(cfg.getParameter<string>("CaloJetAlgorithm4")),
0038       GenJetAlgorithm1(cfg.getParameter<string>("GenJetAlgorithm1")),
0039       GenJetAlgorithm2(cfg.getParameter<string>("GenJetAlgorithm2")),
0040       GenJetAlgorithm3(cfg.getParameter<string>("GenJetAlgorithm3")),
0041       GenJetAlgorithm4(cfg.getParameter<string>("GenJetAlgorithm4")),
0042       JetCorrectionService(cfg.getParameter<string>("JetCorrectionService")) {}
0043 
0044 int nEvent = 0;
0045 
0046 void myFastSimVal::beginJob() {
0047   // Open the histogram file and book some associated histograms
0048   m_file = new TFile("histo.root", "RECREATE");
0049 
0050   tMassGen = TH1F("tMassGen", "T Mass Gen", 100, 0, 200);
0051   tbarMassGen = TH1F("tbarMassGen", "Tbar Mass Gen", 100, 0, 200);
0052 
0053   tMass = TH1F("tMass", "T Mass", 100, 0, 200);
0054   tbarMass = TH1F("tbarMass", "Tbar Mass", 100, 0, 200);
0055 
0056   topMassParton = TH1F("topMassParton", "Top Mass Parton", 100, 0, 200);
0057   topMass1 = TH1F("topMass1", "Top Mass 1", 100, 0, 200);
0058   topMass2 = TH1F("topMass2", "Top Mass 2", 100, 0, 200);
0059   topMass3 = TH1F("topMass3", "Top Mass 3", 100, 0, 200);
0060 
0061   ZpMass = TH1F("ZpMass", "Generated Zp Mass", 160, 0, 8000);
0062   ZpMassGen = TH1F("ZpMassGen", "Gen Zp Mass", 160, 0, 8000);
0063   ZpMassMatched1 = TH1F("ZpMassMatched1", "Calor Zp Mass 1", 160, 0, 8000);
0064   ZpMassMatched2 = TH1F("ZpMassMatched2", "Calor Zp Mass 2", 160, 0, 8000);
0065   ZpMassMatched3 = TH1F("ZpMassMatched3", "Calor Zp Mass 3", 160, 0, 8000);
0066 
0067   ZpMassGen10 = TH1F("ZpMassGen10", "Gen Zp Mass", 160, 0, 8000);
0068   ZpMassGen13 = TH1F("ZpMassGen13", "Gen Zp Mass", 160, 0, 8000);
0069   ZpMassGen40 = TH1F("ZpMassGen40", "Gen Zp Mass", 160, 0, 8000);
0070 
0071   ZpMass_700_10 = TH1F("ZpMass_700_10", "Parton Zp Mass", 100, 0, 1000);
0072   ZpMass_700_13 = TH1F("ZpMass_700_13", "Parton Zp Mass", 100, 0, 1000);
0073   ZpMass_700_40 = TH1F("ZpMass_700_40", "Parton Zp Mass", 100, 0, 1000);
0074 
0075   ZpMassGen_700_10 = TH1F("ZpMassGen_700_10", "Gen Zp Mass", 100, 0, 1000);
0076   ZpMassGen_700_13 = TH1F("ZpMassGen_700_13", "Gen Zp Mass", 100, 0, 1000);
0077   ZpMassGen_700_40 = TH1F("ZpMassGen_700_40", "Gen Zp Mass", 100, 0, 1000);
0078 
0079   ZpMassGen_2000_10 = TH1F("ZpMassGen_2000_10", "Gen Zp Mass", 100, 1500, 2500);
0080   ZpMassGen_2000_13 = TH1F("ZpMassGen_2000_13", "Gen Zp Mass", 100, 1500, 2500);
0081   ZpMassGen_2000_40 = TH1F("ZpMassGen_2000_40", "Gen Zp Mass", 100, 1500, 2500);
0082 
0083   ZpMass_2000_10 = TH1F("ZpMass_2000_10", "Parton Zp Mass", 100, 1500, 2500);
0084   ZpMass_2000_13 = TH1F("ZpMass_2000_13", "Parton Zp Mass", 100, 1500, 2500);
0085   ZpMass_2000_40 = TH1F("ZpMass_2000_40", "Parton Zp Mass", 100, 1500, 2500);
0086 
0087   ZpMassGen_5000_10 = TH1F("ZpMassGen_5000_10", "Gen Zp Mass", 150, 4000, 5500);
0088   ZpMassGen_5000_13 = TH1F("ZpMassGen_5000_13", "Gen Zp Mass", 150, 4000, 5500);
0089   ZpMassGen_5000_40 = TH1F("ZpMassGen_5000_40", "Gen Zp Mass", 150, 4000, 5500);
0090 
0091   ZpMass_5000_10 = TH1F("ZpMass_5000_10", "Parton Zp Mass", 150, 4000, 5500);
0092   ZpMass_5000_13 = TH1F("ZpMass_5000_13", "Parton Zp Mass", 150, 4000, 5500);
0093   ZpMass_5000_40 = TH1F("ZpMass_5000_40", "Parton Zp Mass", 150, 4000, 5500);
0094 
0095   ZpMassRes101 = TH1F("ZpMassRes101", "Zp Mass Resolution 1", 100, -2, 2);
0096   ZpMassRes102 = TH1F("ZpMassRes102", "Zp Mass Resolution 2", 100, -2, 2);
0097   ZpMassRes103 = TH1F("ZpMassRes103", "Zp Mass Resolution 3", 100, -2, 2);
0098 
0099   ZpMassRes131 = TH1F("ZpMassRes131", "Zp Mass Resolution 1", 100, -2, 2);
0100   ZpMassRes132 = TH1F("ZpMassRes132", "Zp Mass Resolution 2", 100, -2, 2);
0101   ZpMassRes133 = TH1F("ZpMassRes133", "Zp Mass Resolution 3", 100, -2, 2);
0102 
0103   ZpMassRes401 = TH1F("ZpMassRes401", "Zp Mass Resolution 1", 100, -2, 2);
0104   ZpMassRes402 = TH1F("ZpMassRes402", "Zp Mass Resolution 2", 100, -2, 2);
0105   ZpMassRes403 = TH1F("ZpMassRes403", "Zp Mass Resolution 3", 100, -2, 2);
0106 
0107   ZpMassResL101 = TH1F("ZpMassResL101", "Zp Mass Resolution Leading Jets 1", 100, 0, 2);
0108   ZpMassResL102 = TH1F("ZpMassResL102", "Zp Mass Resolution Leading Jets 2", 100, 0, 2);
0109   ZpMassResL103 = TH1F("ZpMassResL103", "Zp Mass Resolution Leading Jets 3", 100, 0, 2);
0110 
0111   ZpMassResRL101 = TH1F("ZpMassResRL101", "Zp Mass Res. Ratio Leading Jets 1", 100, 0, 2);
0112   ZpMassResRL102 = TH1F("ZpMassResRL102", "Zp Mass Res. Ratio Leading Jets 2", 100, 0, 2);
0113   ZpMassResRL103 = TH1F("ZpMassResRL103", "Zp Mass Res. Ratio Leading Jets 3", 100, 0, 2);
0114 
0115   ZpMassResRLoP101 = TH1F("ZpMassResRLoP101", "Zp Mass RLoP Ratio Leading Jets 1", 100, 0, 2);
0116   ZpMassResRLoP102 = TH1F("ZpMassResRLoP102", "Zp Mass RLoP Ratio Leading Jets 2", 100, 0, 2);
0117   ZpMassResRLoP103 = TH1F("ZpMassResRLoP103", "Zp Mass RLoP Ratio Leading Jets 3", 100, 0, 2);
0118 
0119   ZpMassResPRL101 = TH1F("ZpMassResPRL101", "Zp Mass Res. P Ratio Leading Jets 1", 100, 0, 2);
0120   ZpMassResPRL102 = TH1F("ZpMassResPRL102", "Zp Mass Res. P Ratio Leading Jets 2", 100, 0, 2);
0121   ZpMassResPRL103 = TH1F("ZpMassResPRL103", "Zp Mass Res. P Ratio Leading Jets 3", 100, 0, 2);
0122 
0123   ZpMassResL131 = TH1F("ZpMassResL131", "Zp Mass Resolution Leading Jets 1", 100, 0, 2);
0124   ZpMassResL132 = TH1F("ZpMassResL132", "Zp Mass Resolution Leading Jets 2", 100, 0, 2);
0125   ZpMassResL133 = TH1F("ZpMassResL133", "Zp Mass Resolution Leading Jets 3", 100, 0, 2);
0126 
0127   ZpMassResRL131 = TH1F("ZpMassResRL131", "Zp Mass Res. Ratio Leading Jets 1", 100, 0, 2);
0128   ZpMassResRL132 = TH1F("ZpMassResRL132", "Zp Mass Res. Ratio Leading Jets 2", 100, 0, 2);
0129   ZpMassResRL133 = TH1F("ZpMassResRL133", "Zp Mass Res. Ratio Leading Jets 3", 100, 0, 2);
0130 
0131   ZpMassResRLoP131 = TH1F("ZpMassResRLoP131", "Zp Mass RLoP Ratio Leading Jets 1", 100, 0, 2);
0132   ZpMassResRLoP132 = TH1F("ZpMassResRLoP132", "Zp Mass RLoP Ratio Leading Jets 2", 100, 0, 2);
0133   ZpMassResRLoP133 = TH1F("ZpMassResRLoP133", "Zp Mass RLoP Ratio Leading Jets 3", 100, 0, 2);
0134 
0135   ZpMassResPRL131 = TH1F("ZpMassResPRL131", "Zp Mass Res. P Ratio Leading Jets 1", 100, 0, 2);
0136   ZpMassResPRL132 = TH1F("ZpMassResPRL132", "Zp Mass Res. P Ratio Leading Jets 2", 100, 0, 2);
0137   ZpMassResPRL133 = TH1F("ZpMassResPRL133", "Zp Mass Res. P Ratio Leading Jets 3", 100, 0, 2);
0138 
0139   ZpMassResL401 = TH1F("ZpMassResL401", "Zp Mass Resolution Leading Jets 1", 100, 0, 2);
0140   ZpMassResL402 = TH1F("ZpMassResL402", "Zp Mass Resolution Leading Jets 2", 100, 0, 2);
0141   ZpMassResL403 = TH1F("ZpMassResL403", "Zp Mass Resolution Leading Jets 3", 100, 0, 2);
0142 
0143   ZpMassResRL401 = TH1F("ZpMassResRL401", "Zp Mass Res. Ratio Leading Jets 1", 100, 0, 2);
0144   ZpMassResRL402 = TH1F("ZpMassResRL402", "Zp Mass Res. Ratio Leading Jets 2", 100, 0, 2);
0145   ZpMassResRL403 = TH1F("ZpMassResRL403", "Zp Mass Res. Ratio Leading Jets 3", 100, 0, 2);
0146 
0147   ZpMassResRLoP401 = TH1F("ZpMassResRLoP401", "Zp Mass RLoP Ratio Leading Jets 1", 100, 0, 2);
0148   ZpMassResRLoP402 = TH1F("ZpMassResRLoP402", "Zp Mass RLoP Ratio Leading Jets 2", 100, 0, 2);
0149   ZpMassResRLoP403 = TH1F("ZpMassResRLoP403", "Zp Mass RLoP Ratio Leading Jets 3", 100, 0, 2);
0150 
0151   ZpMassResPRL401 = TH1F("ZpMassResPRL401", "Zp Mass Res. P Ratio Leading Jets 1", 100, 0, 2);
0152   ZpMassResPRL402 = TH1F("ZpMassResPRL402", "Zp Mass Res. P Ratio Leading Jets 2", 100, 0, 2);
0153   ZpMassResPRL403 = TH1F("ZpMassResPRL403", "Zp Mass Res. P Ratio Leading Jets 3", 100, 0, 2);
0154 
0155   dijetMass1 = TH1F("dijetMass1", "DiJet Mass 1", 100, 0, 4000);
0156   dijetMass12 = TH1F("dijetMass12", "DiJet Mass 1 2", 100, 0, 6000);
0157   dijetMass13 = TH1F("dijetMass13", "DiJet Mass 1 3", 100, 0, 12000);
0158   dijetMass2 = TH1F("dijetMass2", "DiJet Mass 2", 100, 0, 4000);
0159   dijetMass22 = TH1F("dijetMass22", "DiJet Mass 2 2", 100, 0, 6000);
0160   dijetMass23 = TH1F("dijetMass23", "DiJet Mass 2 3", 100, 0, 12000);
0161   dijetMass3 = TH1F("dijetMass3", "DiJet Mass 3", 100, 0, 4000);
0162   dijetMass32 = TH1F("dijetMass32", "DiJet Mass 3 2", 100, 0, 6000);
0163   dijetMass33 = TH1F("dijetMass33", "DiJet Mass 3 3", 100, 0, 12000);
0164   dijetMass4 = TH1F("dijetMass4", "DiJet Mass 4", 100, 0, 4000);
0165   dijetMass42 = TH1F("dijetMass42", "DiJet Mass 4 2", 100, 0, 6000);
0166   dijetMass43 = TH1F("dijetMass43", "DiJet Mass 4 3", 100, 0, 12000);
0167 
0168   dijetMass101 = TH1F("dijetMass101", "DiJet Mass 1", 100, 0, 6000);
0169   dijetMass131 = TH1F("dijetMass131", "DiJet Mass 1", 100, 0, 6000);
0170   dijetMass401 = TH1F("dijetMass401", "DiJet Mass 1", 100, 0, 6000);
0171 
0172   dijetMass102 = TH1F("dijetMass102", "DiJet Mass 2", 100, 0, 6000);
0173   dijetMass132 = TH1F("dijetMass132", "DiJet Mass 2", 100, 0, 6000);
0174   dijetMass402 = TH1F("dijetMass402", "DiJet Mass 2", 100, 0, 6000);
0175 
0176   dijetMass103 = TH1F("dijetMass103", "DiJet Mass 3", 100, 0, 10000);
0177   dijetMass133 = TH1F("dijetMass133", "DiJet Mass 3", 100, 0, 10000);
0178   dijetMass403 = TH1F("dijetMass403", "DiJet Mass 3", 100, 0, 10000);
0179 
0180   dijetMass_700_101 = TH1F("dijetMass_700_101", "DiJet Mass 1", 100, 0, 1000);
0181   dijetMass_700_131 = TH1F("dijetMass_700_131", "DiJet Mass 1", 100, 0, 1000);
0182   dijetMass_700_401 = TH1F("dijetMass_700_401", "DiJet Mass 1", 100, 0, 1000);
0183 
0184   dijetMass_2000_101 = TH1F("dijetMass_2000_101", "DiJet Mass 1", 100, 1500, 2500);
0185   dijetMass_2000_131 = TH1F("dijetMass_2000_131", "DiJet Mass 1", 100, 1500, 2500);
0186   dijetMass_2000_401 = TH1F("dijetMass_2000_401", "DiJet Mass 1", 100, 1500, 2500);
0187 
0188   dijetMass_5000_101 = TH1F("dijetMass_5000_101", "DiJet Mass 1", 150, 4000, 5500);
0189   dijetMass_5000_131 = TH1F("dijetMass_5000_131", "DiJet Mass 1", 150, 4000, 5500);
0190   dijetMass_5000_401 = TH1F("dijetMass_5000_401", "DiJet Mass 1", 150, 4000, 5500);
0191 
0192   dijetMassCor1 = TH1F("dijetMassCor1", "DiJet Mass 1", 160, 0, 8000);
0193   dijetMassCor101 = TH1F("dijetMassCor101", "DiJet Mass Cor 101", 160, 0, 8000);
0194   dijetMassCor131 = TH1F("dijetMassCor131", "DiJet Mass Cor 131", 160, 0, 8000);
0195   dijetMassCor401 = TH1F("dijetMassCor401", "DiJet Mass Cor 401", 160, 0, 8000);
0196 
0197   dijetMassCor_700_1 = TH1F("dijetMassCor_700_1", "DiJet Mass 1", 100, 0, 1000);
0198   dijetMassCor_700_101 = TH1F("dijetMassCor_700_101", "DiJet Mass Cor 101", 100, 0, 1000);
0199   dijetMassCor_700_131 = TH1F("dijetMassCor_700_131", "DiJet Mass Cor 131", 100, 0, 1000);
0200   dijetMassCor_700_401 = TH1F("dijetMassCor_700_401", "DiJet Mass Cor 401", 100, 0, 1000);
0201 
0202   dijetMassCor_2000_1 = TH1F("dijetMassCor_2000_1", "DiJet Mass 1", 100, 1500, 2500);
0203   dijetMassCor_2000_101 = TH1F("dijetMassCor_2000_101", "DiJet Mass Cor 101", 100, 1500, 2500);
0204   dijetMassCor_2000_131 = TH1F("dijetMassCor_2000_131", "DiJet Mass Cor 131", 100, 1500, 2500);
0205   dijetMassCor_2000_401 = TH1F("dijetMassCor_2000_401", "DiJet Mass Cor 401", 100, 1500, 2500);
0206 
0207   dijetMassCor_5000_1 = TH1F("dijetMassCor_5000_1", "DiJet Mass 1", 150, 4000, 5500);
0208   dijetMassCor_5000_101 = TH1F("dijetMassCor_5000_101", "DiJet Mass Cor 101", 150, 4000, 5500);
0209   dijetMassCor_5000_131 = TH1F("dijetMassCor_5000_131", "DiJet Mass Cor 131", 150, 4000, 5500);
0210   dijetMassCor_5000_401 = TH1F("dijetMassCor_5000_401", "DiJet Mass Cor 401", 150, 4000, 5500);
0211 
0212   dijetMassP1 = TH1F("dijetMassP1", "DiJet Mass P 1", 160, 0, 8000);
0213   dijetMassP2 = TH1F("dijetMassP2", "DiJet Mass P 2", 160, 0, 8000);
0214   dijetMassP3 = TH1F("dijetMassP3", "DiJet Mass P 3", 160, 0, 8000);
0215 
0216   dijetMassP101 = TH1F("dijetMassP101", "DiJet Mass P 1", 160, 0, 8000);
0217   dijetMassP131 = TH1F("dijetMassP131", "DiJet Mass P 1", 160, 0, 8000);
0218   dijetMassP401 = TH1F("dijetMassP401", "DiJet Mass P 1", 160, 0, 8000);
0219 
0220   dijetMassP_700_101 = TH1F("dijetMassP_700_101", "DiJet Mass P 1", 100, 0, 1000);
0221   dijetMassP_700_131 = TH1F("dijetMassP_700_131", "DiJet Mass P 1", 100, 0, 1000);
0222   dijetMassP_700_401 = TH1F("dijetMassP_700_401", "DiJet Mass P 1", 100, 0, 1000);
0223 
0224   dijetMassP_2000_101 = TH1F("dijetMassP_2000_101", "DiJet Mass P 1", 100, 1500, 2500);
0225   dijetMassP_2000_131 = TH1F("dijetMassP_2000_131", "DiJet Mass P 1", 100, 1500, 2500);
0226   dijetMassP_2000_401 = TH1F("dijetMassP_2000_401", "DiJet Mass P 1", 100, 1500, 2500);
0227 
0228   dijetMassP_5000_101 = TH1F("dijetMassP_5000_101", "DiJet Mass P 1", 150, 4000, 5500);
0229   dijetMassP_5000_131 = TH1F("dijetMassP_5000_131", "DiJet Mass P 1", 150, 4000, 5500);
0230   dijetMassP_5000_401 = TH1F("dijetMassP_5000_401", "DiJet Mass P 1", 150, 4000, 5500);
0231 
0232   totEneLeadJetEta1_1 = TH1F("totEneLeadJetEta1_1", "Total Energy Lead Jet Eta1 1", 100, 0, 1500);
0233   totEneLeadJetEta2_1 = TH1F("totEneLeadJetEta2_1", "Total Energy Lead Jet Eta2 1", 100, 0, 1500);
0234   totEneLeadJetEta3_1 = TH1F("totEneLeadJetEta3_1", "Total Energy Lead Jet Eta3 1", 100, 0, 1500);
0235   hadEneLeadJetEta1_1 = TH1F("hadEneLeadJetEta1_1", "Hadronic Energy Lead Jet Eta1 1", 100, 0, 1500);
0236   hadEneLeadJetEta2_1 = TH1F("hadEneLeadJetEta2_1", "Hadronic Energy Lead Jet Eta2 1", 100, 0, 1500);
0237   hadEneLeadJetEta3_1 = TH1F("hadEneLeadJetEta3_1", "Hadronic Energy Lead Jet Eta3 1", 100, 0, 1500);
0238   emEneLeadJetEta1_1 = TH1F("emEneLeadJetEta1_1", "EM Energy Lead Jet Eta1 1", 100, 0, 1500);
0239   emEneLeadJetEta2_1 = TH1F("emEneLeadJetEta2_1", "EM Energy Lead Jet Eta2 1", 100, 0, 1500);
0240   emEneLeadJetEta3_1 = TH1F("emEneLeadJetEta3_1", "EM Energy Lead Jet Eta3 1", 100, 0, 1500);
0241 
0242   totEneLeadJetEta1_2 = TH1F("totEneLeadJetEta1_2", "Total Energy Lead Jet Eta1 2", 100, 0, 6000);
0243   totEneLeadJetEta2_2 = TH1F("totEneLeadJetEta2_2", "Total Energy Lead Jet Eta2 2", 100, 0, 6000);
0244   totEneLeadJetEta3_2 = TH1F("totEneLeadJetEta3_2", "Total Energy Lead Jet Eta3 2", 100, 0, 6000);
0245   hadEneLeadJetEta1_2 = TH1F("hadEneLeadJetEta1_2", "Hadronic Energy Lead Jet Eta1 2", 100, 0, 6000);
0246   hadEneLeadJetEta2_2 = TH1F("hadEneLeadJetEta2_2", "Hadronic Energy Lead Jet Eta2 2", 100, 0, 6000);
0247   hadEneLeadJetEta3_2 = TH1F("hadEneLeadJetEta3_2", "Hadronic Energy Lead Jet Eta3 2", 100, 0, 6000);
0248   emEneLeadJetEta1_2 = TH1F("emEneLeadJetEta1_2", "EM Energy Lead Jet Eta1 2", 100, 0, 5000);
0249   emEneLeadJetEta2_2 = TH1F("emEneLeadJetEta2_2", "EM Energy Lead Jet Eta2 2", 100, 0, 5000);
0250   emEneLeadJetEta3_2 = TH1F("emEneLeadJetEta3_2", "EM Energy Lead Jet Eta3 2", 100, 0, 5000);
0251 
0252   hadEneLeadJet1 = TH1F("hadEneLeadJet1", "Hadronic Energy Lead Jet 1", 100, 0, 3000);
0253   hadEneLeadJet12 = TH1F("hadEneLeadJet12", "Hadronic Energy Lead Jet 1 2", 100, 0, 4000);
0254   hadEneLeadJet13 = TH1F("hadEneLeadJet13", "Hadronic Energy Lead Jet 1 3", 100, 0, 6000);
0255   hadEneLeadJet2 = TH1F("hadEneLeadJet2", "Hadronic Energy Lead Jet 2", 100, 0, 3000);
0256   hadEneLeadJet22 = TH1F("hadEneLeadJet22", "Hadronic Energy Lead Jet 2 2", 100, 0, 4000);
0257   hadEneLeadJet23 = TH1F("hadEneLeadJet23", "Hadronic Energy Lead Jet 2 3", 100, 0, 6000);
0258   hadEneLeadJet3 = TH1F("hadEneLeadJet3", "Hadronic Energy Lead Jet 3", 100, 0, 3000);
0259   hadEneLeadJet32 = TH1F("hadEneLeadJet32", "Hadronic Energy Lead Jet 3 2", 100, 0, 4000);
0260   hadEneLeadJet33 = TH1F("hadEneLeadJet33", "Hadronic Energy Lead Jet 3 3", 100, 0, 6000);
0261 
0262   emEneLeadJet1 = TH1F("emEneLeadJet1", "EM Energy Lead Jet 1", 100, 0, 1500);
0263   emEneLeadJet12 = TH1F("emEneLeadJet12", "EM Energy Lead Jet 1 2", 100, 0, 3000);
0264   emEneLeadJet13 = TH1F("emEneLeadJet13", "EM Energy Lead Jet 1 3", 100, 0, 5000);
0265   emEneLeadJet2 = TH1F("emEneLeadJet2", "EM Energy Lead Jet 2", 100, 0, 1500);
0266   emEneLeadJet22 = TH1F("emEneLeadJet22", "EM Energy Lead Jet 2 2", 100, 0, 3000);
0267   emEneLeadJet23 = TH1F("emEneLeadJet23", "EM Energy Lead Jet 2 3", 100, 0, 5000);
0268   emEneLeadJet3 = TH1F("emEneLeadJet3", "EM Energy Lead Jet 3", 100, 0, 1500);
0269   emEneLeadJet32 = TH1F("emEneLeadJet32", "EM Energy Lead Jet 3 2", 100, 0, 3000);
0270   emEneLeadJet33 = TH1F("emEneLeadJet33", "EM Energy Lead Jet 3 3", 100, 0, 5000);
0271 
0272   hadFracEta11 = TH1F("hadFracEta11", "Hadronic Fraction Eta1 Jet 1", 100, 0, 1);
0273   hadFracEta21 = TH1F("hadFracEta21", "Hadronic Fraction Eta2 Jet 1", 100, 0, 1);
0274   hadFracEta31 = TH1F("hadFracEta31", "Hadronic Fraction Eta3 Jet 1", 100, 0, 1);
0275 
0276   hadFracEta12 = TH1F("hadFracEta12", "Hadronic Fraction Eta1 Jet 2", 100, 0, 1);
0277   hadFracEta22 = TH1F("hadFracEta22", "Hadronic Fraction Eta2 Jet 2", 100, 0, 1);
0278   hadFracEta32 = TH1F("hadFracEta32", "Hadronic Fraction Eta3 Jet 2", 100, 0, 1);
0279 
0280   hadFracEta13 = TH1F("hadFracEta13", "Hadronic Fraction Eta1 Jet 3", 100, 0, 1);
0281   hadFracEta23 = TH1F("hadFracEta23", "Hadronic Fraction Eta2 Jet 3", 100, 0, 1);
0282   hadFracEta33 = TH1F("hadFracEta33", "Hadronic Fraction Eta3 Jet 3", 100, 0, 1);
0283 
0284   hadFracLeadJet1 = TH1F("hadFracLeadJet1", "Hadronic Fraction Lead Jet 1", 100, 0, 1);
0285   hadFracLeadJet2 = TH1F("hadFracLeadJet2", "Hadronic Fraction Lead Jet 2", 100, 0, 1);
0286   hadFracLeadJet3 = TH1F("hadFracLeadJet3", "Hadronic Fraction Lead Jet 3", 100, 0, 1);
0287 
0288   SumEt1 = TH1F("SumEt1", "SumEt 1", 100, 0, 1000);
0289   SumEt12 = TH1F("SumEt12", "SumEt 1 2", 100, 0, 4000);
0290   SumEt13 = TH1F("SumEt13", "SumEt 1 3", 100, 0, 15000);
0291 
0292   MET1 = TH1F("MET1", "MET 1", 100, 0, 200);
0293   MET12 = TH1F("MET12", "MET 1 2", 100, 0, 1000);
0294   MET13 = TH1F("MET13", "MET 1 3", 100, 0, 3000);
0295 
0296   nTowersLeadJet1 = TH1F("nTowersLeadJet1", "Number of Towers Lead Jet 1", 100, 0, 100);
0297   nTowersLeadJet2 = TH1F("nTowersLeadJet2", "Number of Towers Lead Jet 2", 100, 0, 100);
0298   nTowersLeadJet3 = TH1F("nTowersLeadJet3", "Number of Towers Lead Jet 3", 100, 0, 100);
0299 
0300   nTowersSecondJet1 = TH1F("nTowersSecondJet1", "Number of Towers Second Jet 1", 100, 0, 100);
0301   nTowersSecondJet2 = TH1F("nTowersSecondJet2", "Number of Towers Second Jet 2", 100, 0, 100);
0302   nTowersSecondJet3 = TH1F("nTowersSecondJet3", "Number of Towers Second Jet 3", 100, 0, 100);
0303 
0304   hf_PtResponse1 = TProfile("PtResponse1", "Pt Response 1", 100, -5, 5, 0, 10);
0305   hf_PtResponse2 = TProfile("PtResponse2", "Pt Response 2", 100, -5, 5, 0, 10);
0306   hf_PtResponse3 = TProfile("PtResponse3", "Pt Response 3", 100, -5, 5, 0, 10);
0307   hf_PtResponse4 = TProfile("PtResponse4", "Pt Response 4", 100, -5, 5, 0, 10);
0308 
0309   hf_TowerDelR1 = TProfile("hf_TowerDelR1", "Tower Del R 1", 100, 0, 2, 0, 10);
0310   hf_TowerDelR12 = TProfile("hf_TowerDelR12", "Tower Del R 1", 80, 0, 0.8, 0, 10);
0311   hf_TowerDelR2 = TProfile("hf_TowerDelR2", "Tower Del R 2", 100, 0, 2, 0, 10);
0312   hf_TowerDelR22 = TProfile("hf_TowerDelR22", "Tower Del R 2", 80, 0, 0.8, 0, 10);
0313   hf_TowerDelR3 = TProfile("hf_TowerDelR3", "Tower Del R 3", 100, 0, 2, 0, 10);
0314   hf_TowerDelR32 = TProfile("hf_TowerDelR32", "Tower Del R 3", 80, 0, 0.8, 0, 10);
0315 
0316   hf_sumTowerAllEx = TH1F("sumTowerAllEx", "Tower Ex", 100, -1000, 1000);
0317   hf_sumTowerAllEy = TH1F("sumTowerAllEy", "Tower Ey", 100, -1000, 1000);
0318 
0319   hf_TowerJetEt1 = TH1F("TowerJetEt1", "Tower/Jet Et 1", 50, 0, 1);
0320 
0321   nTowers1 = TH1F("nTowers1", "Number of Towers pt 0.5", 100, 0, 500);
0322   nTowers2 = TH1F("nTowers2", "Number of Towers pt 1.0", 100, 0, 500);
0323   nTowers3 = TH1F("nTowers3", "Number of Towers pt 1.5", 100, 0, 500);
0324   nTowers4 = TH1F("nTowers4", "Number of Towers pt 2.0", 100, 0, 500);
0325 
0326   nTowersLeadJetPt1 = TH1F("nTowersLeadJetPt1", "Number of Towers in Lead Jet pt 0.5", 100, 0, 200);
0327   nTowersLeadJetPt2 = TH1F("nTowersLeadJetPt2", "Number of Towers in Lead Jet pt 1.0", 100, 0, 200);
0328   nTowersLeadJetPt3 = TH1F("nTowersLeadJetPt3", "Number of Towers in Lead Jet pt 1.5", 100, 0, 200);
0329   nTowersLeadJetPt4 = TH1F("nTowersLeadJetPt4", "Number of Towers in Lead Jet pt 2.0", 100, 0, 200);
0330 
0331   TowerEtLeadJet1 = TH1F("TowerEtLeadJet1", "Towers Et Lead Jet 1", 100, 0, 2000);
0332   TowerEtLeadJet12 = TH1F("TowerEtLeadJet12", "Towers Et Lead Jet 1 2", 100, 0, 6000);
0333   TowerEtLeadJet13 = TH1F("TowerEtLeadJet13", "Towers Et Lead Jet 1 3", 100, 0, 300);
0334   TowerEtLeadJet2 = TH1F("TowerEtLeadJet2", "Towers Et Lead Jet 2", 100, 0, 2000);
0335   TowerEtLeadJet22 = TH1F("TowerEtLeadJet22", "Towers Et Lead Jet 2 2", 100, 0, 6000);
0336   TowerEtLeadJet23 = TH1F("TowerEtLeadJet23", "Towers Et Lead Jet 2 3", 100, 0, 300);
0337   TowerEtLeadJet3 = TH1F("TowerEtLeadJet3", "Towers Et Lead Jet 3", 100, 0, 2000);
0338   TowerEtLeadJet32 = TH1F("TowerEtLeadJet32", "Towers Et Lead Jet 3 2", 100, 0, 6000);
0339   TowerEtLeadJet33 = TH1F("TowerEtLeadJet33", "Towers Et Lead Jet 3 3", 100, 0, 300);
0340 
0341   hf_nJet1 = TProfile("hf_nJet1", "Num Jets 1", 100, 0, 5000, 0, 50);
0342   hf_nJet2 = TProfile("hf_nJet2", "Num Jets 2", 100, 0, 5000, 0, 50);
0343   hf_nJet3 = TProfile("hf_nJet3", "Num Jets 3", 100, 0, 5000, 0, 50);
0344   hf_nJet4 = TProfile("hf_nJet4", "Num Jets 4", 100, 0, 5000, 0, 50);
0345 
0346   hf_nJet1s = TProfile("hf_nJet1s", "Num Jets 1", 100, 0, 200, 0, 50);
0347   hf_nJet2s = TProfile("hf_nJet2s", "Num Jets 2", 100, 0, 200, 0, 50);
0348   hf_nJet3s = TProfile("hf_nJet3s", "Num Jets 3", 100, 0, 200, 0, 50);
0349   hf_nJet4s = TProfile("hf_nJet4s", "Num Jets 4", 100, 0, 200, 0, 50);
0350 
0351   hf_nJet11 = TProfile("hf_nJet11", "Num Jets 1 1", 100, 0, 3000, 0, 50);
0352   hf_nJet21 = TProfile("hf_nJet21", "Num Jets 2 1", 100, 0, 3000, 0, 50);
0353   hf_nJet31 = TProfile("hf_nJet31", "Num Jets 3 1", 100, 0, 3000, 0, 50);
0354   hf_nJet41 = TProfile("hf_nJet41", "Num Jets 4 1", 100, 0, 3000, 0, 50);
0355 
0356   dRPar1 = TH1F("dRPar1", "Parton dR with matched CaloJet1", 100, 0, 0.5);
0357   dPhiPar1 = TH1F("dPhiPar1", "Parton dPhi with matched CaloJet1", 200, -0.5, 0.5);
0358   dEtaPar1 = TH1F("dEtaPar1", "Parton dEta with matched CaloJet1", 200, -0.5, 0.5);
0359   dPtPar1 = TH1F("dPtPar1", "Parton dPt with matched CaloJet1", 200, -50, 50);
0360 
0361   dRPar2 = TH1F("dRPar2", "Parton dR with matched CaloJet2", 100, 0, 0.5);
0362   dPhiPar2 = TH1F("dPhiPar2", "Parton dPhi with matched CaloJet2", 200, -0.5, 0.5);
0363   dEtaPar2 = TH1F("dEtaPar2", "Parton dEta with matched CaloJet2", 200, -0.5, 0.5);
0364   dPtPar2 = TH1F("dPtPar2", "Parton dPt with matched CaloJet2", 200, -50, 50);
0365 
0366   dRPar3 = TH1F("dRPar3", "Parton dR with matched CaloJet3", 100, 0, 0.5);
0367   dPhiPar3 = TH1F("dPhiPar3", "Parton dPhi with matched CaloJet3", 200, -0.5, 0.5);
0368   dEtaPar3 = TH1F("dEtaPar3", "Parton dEta with matched CaloJet3", 200, -0.5, 0.5);
0369   dPtPar3 = TH1F("dPtPar3", "Parton dPt with matched CaloJet3", 200, -50, 50);
0370 
0371   dRPar4 = TH1F("dRPar4", "Parton dR with matched CaloJet4", 100, 0, 0.5);
0372   dPhiPar4 = TH1F("dPhiPar4", "Parton dPhi with matched CaloJet4", 200, -0.5, 0.5);
0373   dEtaPar4 = TH1F("dEtaPar4", "Parton dEta with matched CaloJet4", 200, -0.5, 0.5);
0374   dPtPar4 = TH1F("dPtPar4", "Parton dPt with matched CaloJet4", 200, -50, 50);
0375 
0376   dRParton = TH1F("dRParton", "dR Parton", 100, 0, 10.0);
0377   dRPartonMin = TH1F("dRPartonMin", "Min dR Parton", 100, 0, 2.0);
0378 
0379   dR1 = TH1F("dR1", "GenJets dR with matched CaloJet", 100, 0, 0.5);
0380   dPhi1 = TH1F("dPhi1", "GenJets dPhi with matched CaloJet", 200, -0.5, 0.5);
0381   dEta1 = TH1F("dEta1", "GenJets dEta with matched CaloJet", 200, -0.5, 0.5);
0382   dPt1 = TH1F("dPt1", "GenJets dPt with matched CaloJet", 200, -100, 100);
0383   dPtFrac1 = TH1F("dPtFrac1", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0384   dPt20Frac1 = TH1F("dPt20Frac1", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0385   dPt40Frac1 = TH1F("dPt40Frac1", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0386   dPt80Frac1 = TH1F("dPt80Frac1", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0387   dPt100Frac1 = TH1F("dPt100Frac1", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0388 
0389   dR2 = TH1F("dR2", "GenJets dR with matched CaloJet", 100, 0, 0.5);
0390   dPhi2 = TH1F("dPhi2", "GenJets dPhi with matched CaloJet", 200, -0.5, 0.5);
0391   dEta2 = TH1F("dEta2", "GenJets dEta with matched CaloJet", 200, -0.5, 0.5);
0392   dPt2 = TH1F("dPt2", "GenJets dPt with matched CaloJet", 200, -100, 100);
0393   dPtFrac2 = TH1F("dPtFrac2", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0394   dPt20Frac2 = TH1F("dPt20Frac2", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0395   dPt40Frac2 = TH1F("dPt40Frac2", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0396   dPt80Frac2 = TH1F("dPt80Frac2", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0397   dPt100Frac2 = TH1F("dPt100Frac2", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0398 
0399   dR3 = TH1F("dR3", "GenJets dR with matched CaloJet", 100, 0, 0.5);
0400   dPhi3 = TH1F("dPhi3", "GenJets dPhi with matched CaloJet", 200, -0.5, 0.5);
0401   dEta3 = TH1F("dEta3", "GenJets dEta with matched CaloJet", 200, -0.5, 0.5);
0402   dPt3 = TH1F("dPt3", "GenJets dPt with matched CaloJet", 200, -100, 100);
0403   dPtFrac3 = TH1F("dPtFrac3", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0404   dPt20Frac3 = TH1F("dPt20Frac3", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0405   dPt40Frac3 = TH1F("dPt40Frac3", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0406   dPt80Frac3 = TH1F("dPt80Frac3", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0407   dPt100Frac3 = TH1F("dPt100Frac3", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0408 
0409   dR4 = TH1F("dR4", "GenJets dR with matched CaloJet", 100, 0, 0.5);
0410   dPhi4 = TH1F("dPhi4", "GenJets dPhi with matched CaloJet", 200, -0.5, 0.5);
0411   dEta4 = TH1F("dEta4", "GenJets dEta with matched CaloJet", 200, -0.5, 0.5);
0412   dPt4 = TH1F("dPt4", "GenJets dPt with matched CaloJet", 200, -100, 100);
0413   dPtFrac4 = TH1F("dPtFrac4", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0414   dPt20Frac4 = TH1F("dPt20Frac4", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0415   dPt40Frac4 = TH1F("dPt40Frac4", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0416   dPt80Frac4 = TH1F("dPt80Frac4", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0417   dPt100Frac4 = TH1F("dPt100Frac4", "GenJets dPt frac with matched CaloJet", 100, -1, 1);
0418 
0419   dR12 = TH1F("dR12", "dR MidPoint - SISCone", 100, 0, 0.5);
0420   dPhi12 = TH1F("dPhi12", "dPhi MidPoint - SISCone", 200, -0.5, 0.5);
0421   dEta12 = TH1F("dEta12", "dEta MidPoint - SISCone", 200, -0.5, 0.5);
0422   dPt12 = TH1F("dPt12", "dPt MidPoint - SISCone", 200, -100, 100);
0423 
0424   h_nCalJets1 = TH1F("nCalJets1", "Number of CalJets1", 20, 0, 20);
0425   h_nCalJets2 = TH1F("nCalJets2", "Number of CalJets2", 20, 0, 20);
0426   h_nCalJets3 = TH1F("nCalJets3", "Number of CalJets3", 20, 0, 20);
0427   h_nCalJets4 = TH1F("nCalJets4", "Number of CalJets4", 20, 0, 20);
0428 
0429   h_lowPtCal11 = TH1F("lowPtCal11", "Low p_{T} of CalJet1 Eta1", 100, 0, 100);
0430   h_lowPtCal21 = TH1F("lowPtCal21", "Low p_{T} of CalJet2 Eta1", 100, 0, 100);
0431   h_lowPtCal31 = TH1F("lowPtCal31", "Low p_{T} of CalJet3 Eta1", 100, 0, 100);
0432   h_lowPtCal41 = TH1F("lowPtCal41", "Low p_{T} of CalJet4 Eta1", 100, 0, 100);
0433 
0434   h_lowPtCal12 = TH1F("lowPtCal12", "Low p_{T} of CalJet1 Eta2", 100, 0, 100);
0435   h_lowPtCal22 = TH1F("lowPtCal22", "Low p_{T} of CalJet2 Eta2", 100, 0, 100);
0436   h_lowPtCal32 = TH1F("lowPtCal32", "Low p_{T} of CalJet3 Eta2", 100, 0, 100);
0437   h_lowPtCal42 = TH1F("lowPtCal42", "Low p_{T} of CalJet4 Eta2", 100, 0, 100);
0438 
0439   h_lowPtCal13 = TH1F("lowPtCal13", "Low p_{T} of CalJet1 Eta3", 100, 0, 100);
0440   h_lowPtCal23 = TH1F("lowPtCal23", "Low p_{T} of CalJet2 Eta3", 100, 0, 100);
0441   h_lowPtCal33 = TH1F("lowPtCal33", "Low p_{T} of CalJet3 Eta3", 100, 0, 100);
0442   h_lowPtCal43 = TH1F("lowPtCal43", "Low p_{T} of CalJet4 Eta3", 100, 0, 100);
0443 
0444   h_lowPtCal1c11 = TH1F("lowPtCal1c11", "Low p_{T} of CalJet1 c1 Eta1", 100, 0, 100);
0445   h_lowPtCal2c11 = TH1F("lowPtCal2c11", "Low p_{T} of CalJet2 c1 Eta1", 100, 0, 100);
0446   h_lowPtCal3c11 = TH1F("lowPtCal3c11", "Low p_{T} of CalJet3 c1 Eta1", 100, 0, 100);
0447   h_lowPtCal4c11 = TH1F("lowPtCal4c11", "Low p_{T} of CalJet4 c1 Eta1", 100, 0, 100);
0448 
0449   h_lowPtCal1c12 = TH1F("lowPtCal1c12", "Low p_{T} of CalJet1 c1 Eta2", 100, 0, 100);
0450   h_lowPtCal2c12 = TH1F("lowPtCal2c12", "Low p_{T} of CalJet2 c1 Eta2", 100, 0, 100);
0451   h_lowPtCal3c12 = TH1F("lowPtCal3c12", "Low p_{T} of CalJet3 c1 Eta2", 100, 0, 100);
0452   h_lowPtCal4c12 = TH1F("lowPtCal4c12", "Low p_{T} of CalJet4 c1 Eta2", 100, 0, 100);
0453 
0454   h_lowPtCal1c13 = TH1F("lowPtCal1c13", "Low p_{T} of CalJet1 c1 Eta3", 100, 0, 100);
0455   h_lowPtCal2c13 = TH1F("lowPtCal2c13", "Low p_{T} of CalJet2 c1 Eta3", 100, 0, 100);
0456   h_lowPtCal3c13 = TH1F("lowPtCal3c13", "Low p_{T} of CalJet3 c1 Eta3", 100, 0, 100);
0457   h_lowPtCal4c13 = TH1F("lowPtCal4c13", "Low p_{T} of CalJet4 c1 Eta3", 100, 0, 100);
0458 
0459   matchedAllPt11 = TH1F("matchedAllPt11", "p_{T} of CalJet1 Eta1", 50, 0, 250);
0460   matchedAllPt12 = TH1F("matchedAllPt12", "p_{T} of CalJet1 Eta2", 50, 0, 250);
0461   matchedAllPt13 = TH1F("matchedAllPt13", "p_{T} of CalJet1 Eta3", 50, 0, 250);
0462   matchedPt11 = TH1F("matchedPt11", "p_{T} of CalJet1 Eta1", 50, 0, 250);
0463   matchedPt12 = TH1F("matchedPt12", "p_{T} of CalJet1 Eta2", 50, 0, 250);
0464   matchedPt13 = TH1F("matchedPt13", "p_{T} of CalJet1 Eta3", 50, 0, 250);
0465 
0466   matchedAllPt21 = TH1F("matchedAllPt21", "p_{T} of CalJet2 Eta1", 50, 0, 250);
0467   matchedAllPt22 = TH1F("matchedAllPt22", "p_{T} of CalJet2 Eta2", 50, 0, 250);
0468   matchedAllPt23 = TH1F("matchedAllPt23", "p_{T} of CalJet2 Eta3", 50, 0, 250);
0469   matchedPt21 = TH1F("matchedPt21", "p_{T} of CalJet2 Eta1", 50, 0, 250);
0470   matchedPt22 = TH1F("matchedPt22", "p_{T} of CalJet2 Eta2", 50, 0, 250);
0471   matchedPt23 = TH1F("matchedPt23", "p_{T} of CalJet2 Eta3", 50, 0, 250);
0472 
0473   matchedAllPt31 = TH1F("matchedAllPt31", "p_{T} of CalJet3 Eta1", 50, 0, 250);
0474   matchedAllPt32 = TH1F("matchedAllPt32", "p_{T} of CalJet3 Eta2", 50, 0, 250);
0475   matchedAllPt33 = TH1F("matchedAllPt33", "p_{T} of CalJet3 Eta3", 50, 0, 250);
0476   matchedPt31 = TH1F("matchedPt31", "p_{T} of CalJet3 Eta1", 50, 0, 250);
0477   matchedPt32 = TH1F("matchedPt32", "p_{T} of CalJet3 Eta2", 50, 0, 250);
0478   matchedPt33 = TH1F("matchedPt33", "p_{T} of CalJet3 Eta3", 50, 0, 250);
0479 
0480   matchedAllPt41 = TH1F("matchedAllPt41", "p_{T} of CalJet4 Eta1", 50, 0, 250);
0481   matchedAllPt42 = TH1F("matchedAllPt42", "p_{T} of CalJet4 Eta2", 50, 0, 250);
0482   matchedAllPt43 = TH1F("matchedAllPt43", "p_{T} of CalJet4 Eta3", 50, 0, 250);
0483   matchedPt41 = TH1F("matchedPt41", "p_{T} of CalJet4 Eta1", 50, 0, 250);
0484   matchedPt42 = TH1F("matchedPt42", "p_{T} of CalJet4 Eta2", 50, 0, 250);
0485   matchedPt43 = TH1F("matchedPt43", "p_{T} of CalJet4 Eta3", 50, 0, 250);
0486 
0487   h_ptCal1 = TH1F("ptCal1", "p_{T} of CalJet1", 50, 0, 1000);
0488   h_ptCal12 = TH1F("ptCal12", "p_{T} of CalJet1 2", 50, 0, 6000);
0489   h_ptCal13 = TH1F("ptCal13", "p_{T} of CalJet1 2", 50, 0, 300);
0490 
0491   h_ptCal2 = TH1F("ptCal2", "p_{T} of CalJet2", 50, 0, 1000);
0492   h_ptCal22 = TH1F("ptCal22", "p_{T} of CalJet2 2", 50, 0, 6000);
0493   h_ptCal23 = TH1F("ptCal23", "p_{T} of CalJet2 2", 50, 0, 300);
0494 
0495   h_ptCal3 = TH1F("ptCal3", "p_{T} of CalJet3", 50, 0, 1000);
0496   h_ptCal32 = TH1F("ptCal32", "p_{T} of CalJet3 2", 50, 0, 6000);
0497   h_ptCal33 = TH1F("ptCal33", "p_{T} of CalJet3 2", 50, 0, 300);
0498 
0499   h_ptCal4 = TH1F("ptCal4", "p_{T} of CalJet4", 50, 0, 1000);
0500   h_ptCal42 = TH1F("ptCal42", "p_{T} of CalJet4 2", 50, 0, 6000);
0501   h_ptCal43 = TH1F("ptCal43", "p_{T} of CalJet4 2", 50, 0, 300);
0502 
0503   h_etaCal1 = TH1F("etaCal1", "#eta of  CalJet1", 100, -4, 4);
0504   h_etaCal2 = TH1F("etaCal2", "#eta of  CalJet2", 100, -4, 4);
0505   h_etaCal3 = TH1F("etaCal3", "#eta of  CalJet3", 100, -4, 4);
0506   h_etaCal4 = TH1F("etaCal4", "#eta of  CalJet4", 100, -4, 4);
0507 
0508   h_phiCal1 = TH1F("phiCal1", "#phi of  CalJet1", 50, -M_PI, M_PI);
0509   h_phiCal2 = TH1F("phiCal2", "#phi of  CalJet2", 50, -M_PI, M_PI);
0510   h_phiCal3 = TH1F("phiCal3", "#phi of  CalJet3", 50, -M_PI, M_PI);
0511   h_phiCal4 = TH1F("phiCal4", "#phi of  CalJet4", 50, -M_PI, M_PI);
0512 
0513   h_ptCalL1 = TH1F("ptCalL1", "p_{T} of CalJetL1", 50, 0, 300);
0514   h_ptCalL12 = TH1F("ptCalL12", "p_{T} of CalJetL1 2", 50, 0, 1200);
0515   h_ptCalL13 = TH1F("ptCalL13", "p_{T} of CalJetL1 3", 50, 0, 6000);
0516   h_ptCalL2 = TH1F("ptCalL2", "p_{T} of CalJetL2", 50, 0, 300);
0517   h_ptCalL22 = TH1F("ptCalL22", "p_{T} of CalJetL2 2", 50, 0, 1200);
0518   h_ptCalL23 = TH1F("ptCalL23", "p_{T} of CalJetL2 3", 50, 0, 6000);
0519   h_ptCalL3 = TH1F("ptCalL3", "p_{T} of CalJetL3", 50, 0, 300);
0520   h_ptCalL32 = TH1F("ptCalL32", "p_{T} of CalJetL3 2", 50, 0, 1200);
0521   h_ptCalL33 = TH1F("ptCalL33", "p_{T} of CalJetL3 3", 50, 0, 6000);
0522   h_ptCalL4 = TH1F("ptCalL4", "p_{T} of CalJetL4", 50, 0, 300);
0523   h_ptCalL42 = TH1F("ptCalL42", "p_{T} of CalJetL4 2", 50, 0, 1200);
0524   h_ptCalL43 = TH1F("ptCalL43", "p_{T} of CalJetL4 3", 50, 0, 6000);
0525 
0526   h_etaCalL1 = TH1F("etaCalL1", "#eta of  CalJetL1", 100, -4, 4);
0527   h_etaCalL2 = TH1F("etaCalL2", "#eta of  CalJetL2", 100, -4, 4);
0528   h_etaCalL3 = TH1F("etaCalL3", "#eta of  CalJetL3", 100, -4, 4);
0529   h_etaCalL4 = TH1F("etaCalL4", "#eta of  CalJetL4", 100, -4, 4);
0530   h_phiCalL1 = TH1F("phiCalL1", "#phi of  CalJetL1", 50, -M_PI, M_PI);
0531   h_phiCalL2 = TH1F("phiCalL2", "#phi of  CalJetL2", 50, -M_PI, M_PI);
0532   h_phiCalL3 = TH1F("phiCalL3", "#phi of  CalJetL3", 50, -M_PI, M_PI);
0533   h_phiCalL4 = TH1F("phiCalL4", "#phi of  CalJetL4", 50, -M_PI, M_PI);
0534 
0535   h_nGenJets1 = TH1F("nGenJets1", "Number of GenJets1", 20, 0, 20);
0536   h_nGenJets2 = TH1F("nGenJets2", "Number of GenJets2", 20, 0, 20);
0537   h_nGenJets3 = TH1F("nGenJets3", "Number of GenJets3", 20, 0, 20);
0538   h_nGenJets4 = TH1F("nGenJets4", "Number of GenJets4", 20, 0, 20);
0539 
0540   h_ptGen1 = TH1F("ptGen1", "p_{T} of GenJet1", 50, 0, 1000);
0541   h_ptGen12 = TH1F("ptGen12", "p_{T} of GenJet1 2", 50, 0, 6000);
0542   h_ptGen13 = TH1F("ptGen13", "p_{T} of GenJet1 3", 50, 0, 300);
0543 
0544   h_ptGen2 = TH1F("ptGen2", "p_{T} of GenJet2", 50, 0, 1000);
0545   h_ptGen22 = TH1F("ptGen22", "p_{T} of GenJet2 2", 50, 0, 6000);
0546   h_ptGen23 = TH1F("ptGen23", "p_{T} of GenJet2 3", 50, 0, 300);
0547 
0548   h_ptGen3 = TH1F("ptGen3", "p_{T} of GenJet3", 50, 0, 1000);
0549   h_ptGen32 = TH1F("ptGen32", "p_{T} of GenJet3 2", 50, 0, 6000);
0550   h_ptGen33 = TH1F("ptGen33", "p_{T} of GenJet3 3", 50, 0, 300);
0551 
0552   h_ptGen4 = TH1F("ptGen4", "p_{T} of GenJet4", 50, 0, 1000);
0553   h_ptGen42 = TH1F("ptGen42", "p_{T} of GenJet4 2", 50, 0, 6000);
0554   h_ptGen43 = TH1F("ptGen43", "p_{T} of GenJet4 3", 50, 0, 300);
0555 
0556   h_etaGen1 = TH1F("etaGen1", "#eta of GenJet1", 100, -4, 4);
0557   h_etaGen2 = TH1F("etaGen2", "#eta of GenJet2", 100, -4, 4);
0558   h_etaGen3 = TH1F("etaGen3", "#eta of GenJet3", 100, -4, 4);
0559   h_phiGen1 = TH1F("phiGen1", "#phi of GenJet1", 50, -M_PI, M_PI);
0560   h_phiGen2 = TH1F("phiGen2", "#phi of GenJet2", 50, -M_PI, M_PI);
0561   h_phiGen3 = TH1F("phiGen3", "#phi of GenJet3", 50, -M_PI, M_PI);
0562 
0563   h_ptGenL1 = TH1F("ptGenL1", "p_{T} of GenJetL1", 50, 0, 300);
0564   h_ptGenL12 = TH1F("ptGenL12", "p_{T} of GenJetL1 2", 50, 0, 1200);
0565   h_ptGenL13 = TH1F("ptGenL13", "p_{T} of GenJetL1 3", 50, 0, 6000);
0566   h_ptGenL2 = TH1F("ptGenL2", "p_{T} of GenJetL2", 50, 0, 300);
0567   h_ptGenL22 = TH1F("ptGenL22", "p_{T} of GenJetL2 2", 50, 0, 1200);
0568   h_ptGenL23 = TH1F("ptGenL23", "p_{T} of GenJetL2 3", 50, 0, 6000);
0569   h_ptGenL3 = TH1F("ptGenL3", "p_{T} of GenJetL3", 50, 0, 300);
0570   h_ptGenL32 = TH1F("ptGenL32", "p_{T} of GenJetL32", 50, 0, 1200);
0571   h_ptGenL33 = TH1F("ptGenL33", "p_{T} of GenJetL33", 50, 0, 6000);
0572 
0573   h_etaGenL1 = TH1F("etaGenL1", "#eta of GenJetL1", 100, -4, 4);
0574   h_etaGenL2 = TH1F("etaGenL2", "#eta of GenJetL2", 100, -4, 4);
0575   h_etaGenL3 = TH1F("etaGenL3", "#eta of GenJetL3", 100, -4, 4);
0576   h_phiGenL1 = TH1F("phiGenL1", "#phi of GenJetL1", 50, -M_PI, M_PI);
0577   h_phiGenL2 = TH1F("phiGenL2", "#phi of GenJetL2", 50, -M_PI, M_PI);
0578   h_phiGenL3 = TH1F("phiGenL3", "#phi of GenJetL3", 50, -M_PI, M_PI);
0579 
0580   h_jetEt1 = TH1F("jetEt1", "Total Jet Et", 100, 0, 3000);
0581   h_jetEt2 = TH1F("jetEt2", "Total Jet Et", 100, 0, 3000);
0582   h_jetEt3 = TH1F("jetEt3", "Total Jet Et", 100, 0, 3000);
0583 
0584   h_jet1Pt1 = TH1F("jet1Pt1", "Jet Pt", 100, 0, 3000);
0585   h_jet2Pt1 = TH1F("jet2Pt1", "Jet Pt", 100, 0, 3000);
0586   h_jet3Pt1 = TH1F("jet3Pt1", "Jet Pt", 100, 0, 3000);
0587   h_jet4Pt1 = TH1F("jet4Pt1", "Jet Pt", 100, 0, 3000);
0588   h_jet5Pt1 = TH1F("jet5Pt1", "Jet Pt", 100, 0, 3000);
0589   h_jet6Pt1 = TH1F("jet6Pt1", "Jet Pt", 100, 0, 3000);
0590   h_jet7Pt1 = TH1F("jet7Pt1", "Jet Pt", 100, 0, 3000);
0591 
0592   h_jet1Pt2 = TH1F("jet1Pt2", "Jet Pt", 100, 0, 3000);
0593   h_jet2Pt2 = TH1F("jet2Pt2", "Jet Pt", 100, 0, 3000);
0594   h_jet3Pt2 = TH1F("jet3Pt2", "Jet Pt", 100, 0, 3000);
0595   h_jet4Pt2 = TH1F("jet4Pt2", "Jet Pt", 100, 0, 3000);
0596   h_jet5Pt2 = TH1F("jet5Pt2", "Jet Pt", 100, 0, 3000);
0597   h_jet6Pt2 = TH1F("jet6Pt2", "Jet Pt", 100, 0, 3000);
0598   h_jet7Pt2 = TH1F("jet7Pt2", "Jet Pt", 100, 0, 3000);
0599 
0600   h_jet1Pt3 = TH1F("jet1Pt3", "Jet Pt", 100, 0, 3000);
0601   h_jet2Pt3 = TH1F("jet2Pt3", "Jet Pt", 100, 0, 3000);
0602   h_jet3Pt3 = TH1F("jet3Pt3", "Jet Pt", 100, 0, 3000);
0603   h_jet4Pt3 = TH1F("jet4Pt3", "Jet Pt", 100, 0, 3000);
0604   h_jet5Pt3 = TH1F("jet5Pt3", "Jet Pt", 100, 0, 3000);
0605   h_jet6Pt3 = TH1F("jet6Pt3", "Jet Pt", 100, 0, 3000);
0606   h_jet7Pt3 = TH1F("jet7Pt3", "Jet Pt", 100, 0, 3000);
0607 
0608   h_jet1Pt4 = TH1F("jet1Pt4", "Jet Pt", 100, 0, 3000);
0609   h_jet2Pt4 = TH1F("jet2Pt4", "Jet Pt", 100, 0, 3000);
0610   h_jet3Pt4 = TH1F("jet3Pt4", "Jet Pt", 100, 0, 3000);
0611   h_jet4Pt4 = TH1F("jet4Pt4", "Jet Pt", 100, 0, 3000);
0612   h_jet5Pt4 = TH1F("jet5Pt4", "Jet Pt", 100, 0, 3000);
0613   h_jet6Pt4 = TH1F("jet6Pt4", "Jet Pt", 100, 0, 3000);
0614   h_jet7Pt4 = TH1F("jet7Pt4", "Jet Pt", 100, 0, 3000);
0615 
0616   h_totMissEt1 = TH1F("totMissEt1", "Total Unclustered Et", 100, 0, 500);
0617   h_totMissEt2 = TH1F("totMissEt2", "Total Unclustered Et", 100, 0, 500);
0618   h_totMissEt3 = TH1F("totMissEt3", "Total Unclustered Et", 100, 0, 500);
0619 
0620   h_missEt1 = TH1F("missEt1", "Unclustered Et", 100, 0, 50);
0621   h_missEt2 = TH1F("missEt2", "Unclustered Et", 100, 0, 50);
0622   h_missEt3 = TH1F("missEt3", "Unclustered Et", 100, 0, 50);
0623 
0624   h_missEt1s = TH1F("missEt1s", "Unclustered Et", 100, 0, 2);
0625   h_missEt2s = TH1F("missEt2s", "Unclustered Et", 100, 0, 2);
0626   h_missEt3s = TH1F("missEt3s", "Unclustered Et", 100, 0, 2);
0627 
0628   ParMatch1 = TH1F("ParMatch1", "Number of Matched Jets 1", 10, 0, 10);
0629   ParMatch2 = TH1F("ParMatch2", "Number of Matched Jets 2", 10, 0, 10);
0630   ParMatch3 = TH1F("ParMatch3", "Number of Matched Jets 3", 10, 0, 10);
0631 }
0632 
0633 void myFastSimVal::analyze(const Event& evt, const EventSetup& es) {
0634   int EtaOk10, EtaOk13, EtaOk40;
0635 
0636   double ZpM, ZpMG, ZpMM;
0637   double LeadMass1, LeadMass2, LeadMass3, LeadMass4;
0638   double LeadMassP1, LeadMassP2, LeadMassP3;
0639 
0640   float minJetPt = 30.;
0641   float minJetPt10 = 10.;
0642   int jetInd, allJetInd;
0643   int usedInd = -1;
0644   //  double matchedDelR = 0.1;
0645   double matchedDelR = 0.3;
0646 
0647   ZpMG = 0;
0648   LeadMass1 = -1;
0649   LeadMass2 = -1;
0650   LeadMass3 = -1;
0651 
0652   math::XYZTLorentzVector p4tmp[2], p4cortmp[2];
0653   nEvent++;
0654 
0655   // ********************************
0656   // **** Get the CaloJet1 collection
0657   // ********************************
0658 
0659   Handle<CaloJetCollection> caloJets1;
0660   evt.getByLabel(CaloJetAlgorithm1, caloJets1);
0661 
0662   // Count Jets above Pt cut
0663   for (int istep = 0; istep < 100; ++istep) {
0664     int njet = 0;
0665     float ptStep = (istep * (5000. / 100.));
0666     for (CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++cal) {
0667       if (cal->pt() > ptStep)
0668         njet++;
0669     }
0670 
0671     hf_nJet1.Fill(ptStep, njet);
0672   }
0673 
0674   // Count Jets above Pt cut
0675   for (int istep = 0; istep < 100; ++istep) {
0676     int njet = 0;
0677     float ptStep = (istep * (200. / 100.));
0678     for (CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++cal) {
0679       if (cal->pt() > ptStep)
0680         njet++;
0681     }
0682 
0683     hf_nJet1s.Fill(ptStep, njet);
0684   }
0685 
0686   // Count Jets above Pt cut
0687   for (int istep = 0; istep < 100; ++istep) {
0688     int njet = 0;
0689     float ptStep = (istep * (3000. / 100.));
0690     for (CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++cal) {
0691       if (cal->pt() > ptStep)
0692         njet++;
0693     }
0694 
0695     hf_nJet11.Fill(ptStep, njet);
0696   }
0697 
0698   //Loop over the two leading CaloJets and fill some histograms
0699   jetInd = 0;
0700   allJetInd = 0;
0701   EtaOk10 = 0;
0702   EtaOk13 = 0;
0703   EtaOk40 = 0;
0704 
0705   double highestPt;
0706   double nextPt;
0707 
0708   highestPt = 0.0;
0709   nextPt = 0.0;
0710 
0711   for (CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++cal) {
0712     //    double scale = corrector->correction (*cal);
0713     double scale = 1.0;
0714     double corPt = scale * cal->pt();
0715     //    double corPt = cal->pt();
0716 
0717     if (corPt > highestPt) {
0718       nextPt = highestPt;
0719       p4cortmp[1] = p4cortmp[0];
0720       highestPt = corPt;
0721       p4cortmp[0] = scale * cal->p4();
0722     } else if (corPt > nextPt) {
0723       nextPt = corPt;
0724       p4cortmp[1] = scale * cal->p4();
0725     }
0726 
0727     /***
0728     std::cout << ">>> Corr Jet: corPt = " 
0729           << corPt << ", scale = " << scale 
0730           << " pt = " << cal->pt()
0731           << " highestPt = " << highestPt 
0732           << " nextPt = "    << nextPt 
0733           << std::endl;  
0734     ****/
0735 
0736     allJetInd++;
0737     if (allJetInd == 1) {
0738       h_jet1Pt1.Fill(cal->pt());
0739       p4tmp[0] = cal->p4();
0740       if (fabs(cal->eta()) < 1.0)
0741         EtaOk10++;
0742       if (fabs(cal->eta()) < 1.3)
0743         EtaOk13++;
0744       if (fabs(cal->eta()) < 4.0)
0745         EtaOk40++;
0746     }
0747     if (allJetInd == 2) {
0748       h_jet2Pt1.Fill(cal->pt());
0749       p4tmp[1] = cal->p4();
0750       if (fabs(cal->eta()) < 1.0)
0751         EtaOk10++;
0752       if (fabs(cal->eta()) < 1.3)
0753         EtaOk13++;
0754       if (fabs(cal->eta()) < 4.0)
0755         EtaOk40++;
0756     }
0757     if ((allJetInd == 1) || (allJetInd == 2)) {
0758       h_ptCalL1.Fill(cal->pt());
0759       h_ptCalL12.Fill(cal->pt());
0760       h_ptCalL13.Fill(cal->pt());
0761 
0762       h_etaCalL1.Fill(cal->eta());
0763       h_phiCalL1.Fill(cal->phi());
0764     }
0765 
0766     if (allJetInd == 3)
0767       h_jet3Pt1.Fill(cal->pt());
0768     if (allJetInd == 4)
0769       h_jet4Pt1.Fill(cal->pt());
0770     if (allJetInd == 5)
0771       h_jet5Pt1.Fill(cal->pt());
0772     if (allJetInd == 6)
0773       h_jet6Pt1.Fill(cal->pt());
0774     if (allJetInd == 7)
0775       h_jet7Pt1.Fill(cal->pt());
0776 
0777     if (fabs(cal->eta()) < 1.3) {
0778       h_lowPtCal11.Fill(cal->pt());
0779       if (cal->pt() > 10.) {
0780         h_lowPtCal1c11.Fill(cal->pt());
0781       }
0782     }
0783     if ((fabs(cal->eta()) > 1.3) && (fabs(cal->eta()) < 3.)) {
0784       h_lowPtCal12.Fill(cal->pt());
0785       if (cal->pt() > 10.) {
0786         h_lowPtCal1c12.Fill(cal->pt());
0787       }
0788     }
0789     if (fabs(cal->eta()) > 3.0) {
0790       h_lowPtCal13.Fill(cal->pt());
0791       if (cal->pt() > 10.) {
0792         h_lowPtCal1c13.Fill(cal->pt());
0793       }
0794     }
0795 
0796     if (cal->pt() > minJetPt) {
0797       //    std::cout << "CALO JET1 #" << jetInd << std::endl << cal->print() << std::endl;
0798       h_ptCal1.Fill(cal->pt());
0799       h_ptCal12.Fill(cal->pt());
0800       h_ptCal13.Fill(cal->pt());
0801 
0802       h_etaCal1.Fill(cal->eta());
0803       h_phiCal1.Fill(cal->phi());
0804       jetInd++;
0805     }
0806   }
0807 
0808   //  h_nCalJets1.Fill( caloJets1->size() );
0809   h_nCalJets1.Fill(jetInd);
0810   if (jetInd > 1) {
0811     LeadMass1 = (p4tmp[0] + p4tmp[1]).mass();
0812     dijetMass1.Fill(LeadMass1);
0813     dijetMass12.Fill(LeadMass1);
0814     dijetMass13.Fill(LeadMass1);
0815     if (EtaOk10 == 2) {
0816       dijetMass101.Fill(LeadMass1);
0817       dijetMass102.Fill(LeadMass1);
0818       dijetMass103.Fill(LeadMass1);
0819       dijetMass_700_101.Fill(LeadMass1);
0820       dijetMass_2000_101.Fill(LeadMass1);
0821       dijetMass_5000_101.Fill(LeadMass1);
0822     }
0823     if (EtaOk13 == 2) {
0824       dijetMass131.Fill(LeadMass1);
0825       dijetMass132.Fill(LeadMass1);
0826       dijetMass133.Fill(LeadMass1);
0827       dijetMass_700_131.Fill(LeadMass1);
0828       dijetMass_2000_131.Fill(LeadMass1);
0829       dijetMass_5000_131.Fill(LeadMass1);
0830     }
0831     if (EtaOk40 == 2) {
0832       dijetMass401.Fill(LeadMass1);
0833       dijetMass402.Fill(LeadMass1);
0834       dijetMass403.Fill(LeadMass1);
0835       dijetMass_700_401.Fill(LeadMass1);
0836       dijetMass_2000_401.Fill(LeadMass1);
0837       dijetMass_5000_401.Fill(LeadMass1);
0838     }
0839 
0840     LeadMass1 = (p4cortmp[0] + p4cortmp[1]).mass();
0841 
0842     /****
0843     if (LeadMass1 < 30.) {
0844       std::cout << " XXX Low Mass " 
0845         << (p4tmp[0]+p4tmp[1]).mass() 
0846         << " / " 
0847         << (p4cortmp[0]+p4cortmp[1]).mass() 
0848         << std::endl;
0849 
0850       std::cout << " p4 1 = " << p4tmp[0]
0851         << " p4 2 = " << p4tmp[1]
0852         << " p4 cor 1 = " << p4cortmp[0]
0853         << " p4 cor 2 = " << p4cortmp[0]
0854         << endl;
0855 
0856     }
0857     ****/
0858 
0859     /****
0860     dijetMassCor1.Fill( LeadMass1 );    
0861     dijetMassCor_700_1.Fill( LeadMass1 );    
0862     dijetMassCor_2000_1.Fill( LeadMass1 );    
0863     dijetMassCor_5000_1.Fill( LeadMass1 );    
0864 
0865     if (EtaOk10 == 2) {
0866       dijetMassCor101.Fill( LeadMass1 );  
0867       dijetMassCor_700_101.Fill( LeadMass1 );  
0868       dijetMassCor_2000_101.Fill( LeadMass1 );  
0869       dijetMassCor_5000_101.Fill( LeadMass1 );  
0870     }
0871     if (EtaOk13 == 2) {
0872       dijetMassCor131.Fill( LeadMass1 );  
0873       dijetMassCor_700_131.Fill( LeadMass1 );  
0874       dijetMassCor_2000_131.Fill( LeadMass1 );  
0875       dijetMassCor_5000_131.Fill( LeadMass1 );  
0876     }
0877     if (EtaOk40 == 2) {
0878       dijetMassCor401.Fill( LeadMass1 ); 
0879       dijetMassCor_700_401.Fill( LeadMass1 ); 
0880       dijetMassCor_2000_401.Fill( LeadMass1 ); 
0881       dijetMassCor_5000_401.Fill( LeadMass1 ); 
0882     }
0883     ****/
0884   }
0885 
0886   // ********************************
0887   // **** Get the CaloJet2 collection
0888   // ********************************
0889   Handle<CaloJetCollection> caloJets2;
0890   evt.getByLabel(CaloJetAlgorithm2, caloJets2);
0891 
0892   // Count Jets above Pt cut
0893   for (int istep = 0; istep < 100; ++istep) {
0894     int njet = 0;
0895     float ptStep = (istep * (5000. / 100.));
0896 
0897     for (CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal)
0898       if (cal->pt() > ptStep)
0899         njet++;
0900 
0901     hf_nJet2.Fill(ptStep, njet);
0902   }
0903 
0904   for (int istep = 0; istep < 100; ++istep) {
0905     int njet = 0;
0906     float ptStep = (istep * (200. / 100.));
0907 
0908     for (CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal)
0909       if (cal->pt() > ptStep)
0910         njet++;
0911 
0912     hf_nJet2s.Fill(ptStep, njet);
0913   }
0914 
0915   for (int istep = 0; istep < 100; ++istep) {
0916     int njet = 0;
0917     float ptStep = (istep * (3000. / 100.));
0918 
0919     for (CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal)
0920       if (cal->pt() > ptStep)
0921         njet++;
0922 
0923     hf_nJet21.Fill(ptStep, njet);
0924   }
0925 
0926   //Loop over the two leading CaloJets and fill some histograms
0927   jetInd = 0;
0928   allJetInd = 0;
0929   for (CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal) {
0930     allJetInd++;
0931     if (allJetInd == 1) {
0932       h_jet1Pt2.Fill(cal->pt());
0933       p4tmp[0] = cal->p4();
0934     }
0935     if (allJetInd == 2) {
0936       h_jet2Pt2.Fill(cal->pt());
0937       p4tmp[1] = cal->p4();
0938     }
0939     if ((allJetInd == 1) || (allJetInd == 2)) {
0940       h_ptCalL2.Fill(cal->pt());
0941       h_ptCalL22.Fill(cal->pt());
0942       h_ptCalL23.Fill(cal->pt());
0943 
0944       h_etaCalL2.Fill(cal->eta());
0945       h_phiCalL2.Fill(cal->phi());
0946     }
0947     if (allJetInd == 3)
0948       h_jet3Pt2.Fill(cal->pt());
0949     if (allJetInd == 4)
0950       h_jet4Pt2.Fill(cal->pt());
0951     if (allJetInd == 5)
0952       h_jet5Pt2.Fill(cal->pt());
0953     if (allJetInd == 6)
0954       h_jet6Pt2.Fill(cal->pt());
0955     if (allJetInd == 7)
0956       h_jet7Pt2.Fill(cal->pt());
0957 
0958     if (fabs(cal->eta()) < 1.3) {
0959       h_lowPtCal21.Fill(cal->pt());
0960       if (cal->pt() > 10.) {
0961         h_lowPtCal2c11.Fill(cal->pt());
0962       }
0963     }
0964     if ((fabs(cal->eta()) > 1.3) && (fabs(cal->eta()) < 3.)) {
0965       h_lowPtCal22.Fill(cal->pt());
0966       if (cal->pt() > 10.) {
0967         h_lowPtCal2c12.Fill(cal->pt());
0968       }
0969     }
0970     if (fabs(cal->eta()) > 3.0) {
0971       h_lowPtCal23.Fill(cal->pt());
0972       if (cal->pt() > 10.) {
0973         h_lowPtCal2c13.Fill(cal->pt());
0974       }
0975     }
0976 
0977     if (cal->pt() > minJetPt) {
0978       h_ptCal2.Fill(cal->pt());
0979       h_ptCal22.Fill(cal->pt());
0980       h_ptCal23.Fill(cal->pt());
0981 
0982       h_etaCal2.Fill(cal->eta());
0983       h_phiCal2.Fill(cal->phi());
0984       jetInd++;
0985     }
0986   }
0987   //  h_nCalJets2.Fill( caloJets2->size() );
0988   h_nCalJets2.Fill(jetInd);
0989   if (jetInd > 1) {
0990     LeadMass2 = (p4tmp[0] + p4tmp[1]).mass();
0991     dijetMass2.Fill(LeadMass2);
0992     dijetMass22.Fill(LeadMass2);
0993     dijetMass23.Fill(LeadMass2);
0994 
0995     dijetMassCor1.Fill(LeadMass2);
0996     dijetMassCor_700_1.Fill(LeadMass2);
0997     dijetMassCor_2000_1.Fill(LeadMass2);
0998     dijetMassCor_5000_1.Fill(LeadMass2);
0999 
1000     if (EtaOk10 == 2) {
1001       dijetMassCor101.Fill(LeadMass2);
1002       dijetMassCor_700_101.Fill(LeadMass2);
1003       dijetMassCor_2000_101.Fill(LeadMass2);
1004       dijetMassCor_5000_101.Fill(LeadMass2);
1005     }
1006     if (EtaOk13 == 2) {
1007       dijetMassCor131.Fill(LeadMass2);
1008       dijetMassCor_700_131.Fill(LeadMass2);
1009       dijetMassCor_2000_131.Fill(LeadMass2);
1010       dijetMassCor_5000_131.Fill(LeadMass2);
1011     }
1012     if (EtaOk40 == 2) {
1013       dijetMassCor401.Fill(LeadMass2);
1014       dijetMassCor_700_401.Fill(LeadMass2);
1015       dijetMassCor_2000_401.Fill(LeadMass2);
1016       dijetMassCor_5000_401.Fill(LeadMass2);
1017     }
1018   }
1019 
1020   // ********************************
1021   // **** Get the CaloJet3 collection
1022   // ********************************
1023   Handle<CaloJetCollection> caloJets3;
1024   evt.getByLabel(CaloJetAlgorithm3, caloJets3);
1025 
1026   //Loop over the two leading CaloJets and fill some histograms
1027   jetInd = 0;
1028   allJetInd = 0;
1029 
1030   // Count Jets above Pt cut
1031   for (int istep = 0; istep < 100; ++istep) {
1032     int njet = 0;
1033     float ptStep = (istep * (5000. / 100.));
1034 
1035     for (CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++cal)
1036       if (cal->pt() > ptStep)
1037         njet++;
1038 
1039     hf_nJet3.Fill(ptStep, njet);
1040   }
1041 
1042   for (int istep = 0; istep < 100; ++istep) {
1043     int njet = 0;
1044     float ptStep = (istep * (200. / 100.));
1045 
1046     for (CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++cal)
1047       if (cal->pt() > ptStep)
1048         njet++;
1049 
1050     hf_nJet3s.Fill(ptStep, njet);
1051   }
1052 
1053   for (int istep = 0; istep < 100; ++istep) {
1054     int njet = 0;
1055     float ptStep = (istep * (3000. / 100.));
1056 
1057     for (CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++cal)
1058       if (cal->pt() > ptStep)
1059         njet++;
1060 
1061     hf_nJet31.Fill(ptStep, njet);
1062   }
1063 
1064   for (CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++cal) {
1065     allJetInd++;
1066     if (allJetInd == 1) {
1067       h_jet1Pt3.Fill(cal->pt());
1068       p4tmp[0] = cal->p4();
1069     }
1070     if (allJetInd == 2) {
1071       h_jet2Pt3.Fill(cal->pt());
1072       p4tmp[1] = cal->p4();
1073     }
1074     if ((allJetInd == 1) || (allJetInd == 2)) {
1075       h_ptCalL3.Fill(cal->pt());
1076       h_ptCalL32.Fill(cal->pt());
1077       h_ptCalL33.Fill(cal->pt());
1078 
1079       h_etaCalL3.Fill(cal->eta());
1080       h_phiCalL3.Fill(cal->phi());
1081     }
1082     if (allJetInd == 3)
1083       h_jet3Pt3.Fill(cal->pt());
1084     if (allJetInd == 4)
1085       h_jet4Pt3.Fill(cal->pt());
1086     if (allJetInd == 5)
1087       h_jet5Pt3.Fill(cal->pt());
1088     if (allJetInd == 6)
1089       h_jet6Pt3.Fill(cal->pt());
1090     if (allJetInd == 7)
1091       h_jet7Pt3.Fill(cal->pt());
1092 
1093     if (fabs(cal->eta()) < 1.3) {
1094       h_lowPtCal31.Fill(cal->pt());
1095       if (cal->pt() > 10.) {
1096         h_lowPtCal3c11.Fill(cal->pt());
1097       }
1098     }
1099     if ((fabs(cal->eta()) > 1.3) && (fabs(cal->eta()) < 3.)) {
1100       h_lowPtCal32.Fill(cal->pt());
1101       if (cal->pt() > 10.) {
1102         h_lowPtCal3c12.Fill(cal->pt());
1103       }
1104     }
1105     if (fabs(cal->eta()) > 3.0) {
1106       h_lowPtCal33.Fill(cal->pt());
1107       if (cal->pt() > 10.) {
1108         h_lowPtCal3c13.Fill(cal->pt());
1109       }
1110     }
1111 
1112     if (cal->pt() > minJetPt) {
1113       //    std::cout << "CALO JET3 #" << jetInd << std::endl << cal->print() << std::endl;
1114       h_ptCal3.Fill(cal->pt());
1115       h_ptCal32.Fill(cal->pt());
1116       h_ptCal33.Fill(cal->pt());
1117 
1118       h_etaCal3.Fill(cal->eta());
1119       h_phiCal3.Fill(cal->phi());
1120       jetInd++;
1121     }
1122   }
1123   //  h_nCalJets3.Fill( caloJets3->size() );
1124   h_nCalJets3.Fill(jetInd);
1125   if (jetInd > 1) {
1126     LeadMass3 = (p4tmp[0] + p4tmp[1]).mass();
1127     dijetMass3.Fill(LeadMass3);
1128     dijetMass32.Fill(LeadMass3);
1129     dijetMass33.Fill(LeadMass3);
1130   }
1131 
1132   // ********************************
1133   // **** Get the CaloJet4 collection
1134   // ********************************
1135 
1136   Handle<CaloJetCollection> caloJets4;
1137   evt.getByLabel(CaloJetAlgorithm4, caloJets4);
1138 
1139   //Loop over the two leading CaloJets and fill some histograms
1140   jetInd = 0;
1141   allJetInd = 0;
1142 
1143   // Count Jets above Pt cut
1144   for (int istep = 0; istep < 100; ++istep) {
1145     int njet = 0;
1146     float ptStep = (istep * (5000. / 100.));
1147 
1148     for (CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++cal)
1149       if (cal->pt() > ptStep)
1150         njet++;
1151 
1152     hf_nJet4.Fill(ptStep, njet);
1153   }
1154 
1155   for (int istep = 0; istep < 100; ++istep) {
1156     int njet = 0;
1157     float ptStep = (istep * (200. / 100.));
1158 
1159     for (CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++cal)
1160       if (cal->pt() > ptStep)
1161         njet++;
1162 
1163     hf_nJet4s.Fill(ptStep, njet);
1164   }
1165 
1166   for (int istep = 0; istep < 100; ++istep) {
1167     int njet = 0;
1168     float ptStep = (istep * (3000. / 100.));
1169 
1170     for (CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++cal)
1171       if (cal->pt() > ptStep)
1172         njet++;
1173 
1174     hf_nJet41.Fill(ptStep, njet);
1175   }
1176 
1177   for (CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++cal) {
1178     allJetInd++;
1179     if (allJetInd == 1) {
1180       h_jet1Pt4.Fill(cal->pt());
1181       p4tmp[0] = cal->p4();
1182     }
1183     if (allJetInd == 2) {
1184       h_jet2Pt4.Fill(cal->pt());
1185       p4tmp[1] = cal->p4();
1186     }
1187     if ((allJetInd == 1) || (allJetInd == 2)) {
1188       h_ptCalL4.Fill(cal->pt());
1189       h_ptCalL42.Fill(cal->pt());
1190       h_ptCalL43.Fill(cal->pt());
1191 
1192       h_etaCalL4.Fill(cal->eta());
1193       h_phiCalL4.Fill(cal->phi());
1194     }
1195     if (allJetInd == 3)
1196       h_jet3Pt4.Fill(cal->pt());
1197     if (allJetInd == 4)
1198       h_jet4Pt4.Fill(cal->pt());
1199     if (allJetInd == 5)
1200       h_jet5Pt4.Fill(cal->pt());
1201     if (allJetInd == 6)
1202       h_jet6Pt4.Fill(cal->pt());
1203     if (allJetInd == 7)
1204       h_jet7Pt4.Fill(cal->pt());
1205 
1206     if (fabs(cal->eta()) < 1.3) {
1207       h_lowPtCal41.Fill(cal->pt());
1208       if (cal->pt() > 10.) {
1209         h_lowPtCal4c11.Fill(cal->pt());
1210       }
1211     }
1212     if ((fabs(cal->eta()) > 1.3) && (fabs(cal->eta()) < 3.)) {
1213       h_lowPtCal42.Fill(cal->pt());
1214       if (cal->pt() > 10.) {
1215         h_lowPtCal4c12.Fill(cal->pt());
1216       }
1217     }
1218     if (fabs(cal->eta()) > 3.0) {
1219       h_lowPtCal43.Fill(cal->pt());
1220       if (cal->pt() > 10.) {
1221         h_lowPtCal4c13.Fill(cal->pt());
1222       }
1223     }
1224 
1225     if (cal->pt() > minJetPt) {
1226       //    std::cout << "CALO JET4 #" << jetInd << std::endl << cal->print() << std::endl;
1227       h_ptCal4.Fill(cal->pt());
1228       h_ptCal42.Fill(cal->pt());
1229       h_ptCal43.Fill(cal->pt());
1230 
1231       h_etaCal4.Fill(cal->eta());
1232       h_phiCal4.Fill(cal->phi());
1233       jetInd++;
1234     }
1235   }
1236   //  h_nCalJets4.Fill( caloJets4->size() );
1237   h_nCalJets4.Fill(jetInd);
1238   if (jetInd > 1) {
1239     LeadMass4 = (p4tmp[0] + p4tmp[1]).mass();
1240     dijetMass4.Fill(LeadMass4);
1241     dijetMass42.Fill(LeadMass4);
1242     dijetMass43.Fill(LeadMass4);
1243   }
1244 
1245   // *********************************************
1246   // *********************************************
1247 
1248   //**** Get the GenJet1 collection
1249   Handle<GenJetCollection> genJets1;
1250   evt.getByLabel(GenJetAlgorithm1, genJets1);
1251 
1252   //Loop over the two leading GenJets and fill some histograms
1253   jetInd = 0;
1254   allJetInd = 0;
1255   for (GenJetCollection::const_iterator gen = genJets1->begin(); gen != genJets1->end(); ++gen) {
1256     allJetInd++;
1257     if (allJetInd == 1) {
1258       p4tmp[0] = gen->p4();
1259     }
1260     if (allJetInd == 2) {
1261       p4tmp[1] = gen->p4();
1262     }
1263 
1264     if ((allJetInd == 1) || (allJetInd == 2)) {
1265       h_ptGenL1.Fill(gen->pt());
1266       h_ptGenL12.Fill(gen->pt());
1267       h_ptGenL13.Fill(gen->pt());
1268 
1269       h_etaGenL1.Fill(gen->eta());
1270       h_phiGenL1.Fill(gen->phi());
1271     }
1272 
1273     if (gen->pt() > minJetPt) {
1274       // std::cout << "GEN JET1 #" << jetInd << std::endl << gen->print() << std::endl;
1275       h_ptGen1.Fill(gen->pt());
1276       h_ptGen12.Fill(gen->pt());
1277       h_ptGen13.Fill(gen->pt());
1278 
1279       h_etaGen1.Fill(gen->eta());
1280       h_phiGen1.Fill(gen->phi());
1281       jetInd++;
1282     }
1283   }
1284 
1285   LeadMassP1 = (p4tmp[0] + p4tmp[1]).mass();
1286   dijetMassP1.Fill(LeadMassP1);
1287   if (EtaOk10 == 2) {
1288     dijetMassP101.Fill(LeadMassP1);
1289     dijetMassP_700_101.Fill(LeadMassP1);
1290     dijetMassP_2000_101.Fill(LeadMassP1);
1291     dijetMassP_5000_101.Fill(LeadMassP1);
1292   }
1293   if (EtaOk13 == 2) {
1294     dijetMassP131.Fill(LeadMassP1);
1295     dijetMassP_700_131.Fill(LeadMassP1);
1296     dijetMassP_2000_131.Fill(LeadMassP1);
1297     dijetMassP_5000_131.Fill(LeadMassP1);
1298   }
1299   if (EtaOk40 == 2) {
1300     dijetMassP401.Fill(LeadMassP1);
1301     dijetMassP_5000_401.Fill(LeadMassP1);
1302     dijetMassP_5000_401.Fill(LeadMassP1);
1303     dijetMassP_5000_401.Fill(LeadMassP1);
1304   }
1305 
1306   //  h_nGenJets1.Fill( genJets1->size() );
1307   h_nGenJets1.Fill(jetInd);
1308 
1309   //**** Get the GenJet2 collection
1310   Handle<GenJetCollection> genJets2;
1311   evt.getByLabel(GenJetAlgorithm2, genJets2);
1312 
1313   //Loop over the two leading GenJets and fill some histograms
1314   jetInd = 0;
1315   allJetInd = 0;
1316   for (GenJetCollection::const_iterator gen = genJets2->begin(); gen != genJets2->end(); ++gen) {
1317     allJetInd++;
1318     if (allJetInd == 1) {
1319       p4tmp[0] = gen->p4();
1320     }
1321     if (allJetInd == 2) {
1322       p4tmp[1] = gen->p4();
1323     }
1324     if ((allJetInd == 1) || (allJetInd == 2)) {
1325       h_ptGenL2.Fill(gen->pt());
1326       h_ptGenL22.Fill(gen->pt());
1327       h_ptGenL23.Fill(gen->pt());
1328 
1329       h_etaGenL2.Fill(gen->eta());
1330       h_phiGenL2.Fill(gen->phi());
1331     }
1332 
1333     if (gen->pt() > minJetPt) {
1334       // std::cout << "GEN JET2 #" << jetInd << std::endl << gen->print() << std::endl;
1335       h_ptGen2.Fill(gen->pt());
1336       h_ptGen22.Fill(gen->pt());
1337       h_ptGen23.Fill(gen->pt());
1338 
1339       h_etaGen2.Fill(gen->eta());
1340       h_phiGen2.Fill(gen->phi());
1341       jetInd++;
1342     }
1343   }
1344 
1345   LeadMassP2 = (p4tmp[0] + p4tmp[1]).mass();
1346   dijetMassP2.Fill(LeadMassP2);
1347 
1348   //  h_nGenJets2.Fill( genJets2->size() );
1349   h_nGenJets2.Fill(jetInd);
1350 
1351   //**** Get the GenJet3 collection
1352   Handle<GenJetCollection> genJets3;
1353   evt.getByLabel(GenJetAlgorithm3, genJets3);
1354 
1355   //Loop over the two leading GenJets and fill some histograms
1356   jetInd = 0;
1357   allJetInd = 0;
1358   for (GenJetCollection::const_iterator gen = genJets3->begin(); gen != genJets3->end(); ++gen) {
1359     allJetInd++;
1360     if (allJetInd == 1) {
1361       p4tmp[0] = gen->p4();
1362     }
1363     if (allJetInd == 2) {
1364       p4tmp[1] = gen->p4();
1365     }
1366     if ((allJetInd == 1) || (allJetInd == 2)) {
1367       h_ptGenL3.Fill(gen->pt());
1368       h_ptGenL32.Fill(gen->pt());
1369       h_ptGenL33.Fill(gen->pt());
1370 
1371       h_etaGenL3.Fill(gen->eta());
1372       h_phiGenL3.Fill(gen->phi());
1373     }
1374 
1375     if (gen->pt() > minJetPt) {
1376       // std::cout << "GEN JET3 #" << jetInd << std::endl << gen->print() << std::endl;
1377       h_ptGen3.Fill(gen->pt());
1378       h_ptGen32.Fill(gen->pt());
1379       h_ptGen33.Fill(gen->pt());
1380 
1381       h_etaGen3.Fill(gen->eta());
1382       h_phiGen3.Fill(gen->phi());
1383       jetInd++;
1384     }
1385   }
1386 
1387   LeadMassP3 = (p4tmp[0] + p4tmp[1]).mass();
1388   dijetMassP3.Fill(LeadMassP3);
1389 
1390   //  h_nGenJets3.Fill( genJets3->size() );
1391   h_nGenJets3.Fill(jetInd);
1392 
1393   // *********************
1394   // MidPoint Jet Matching
1395 
1396   Handle<GenJetCollection> genJets;
1397   Handle<CaloJetCollection> caloJets;
1398 
1399   //  evt.getByLabel( "midPointCone5GenJets",  genJets );
1400   //  evt.getByLabel( "midPointCone5CaloJets", caloJets );
1401   evt.getByLabel(GenJetAlgorithm1, genJets);
1402   evt.getByLabel(CaloJetAlgorithm1, caloJets);
1403 
1404   int maxJets = MAXJETS;
1405 
1406   jetInd = 0;
1407   double dRmin[MAXJETS];
1408   math::XYZTLorentzVector p4gen[MAXJETS], p4cal[MAXJETS], p4par[MAXJETS], p4Zp[MAXJETS], p4part[MAXJETS];
1409 
1410   int used[MAXJETS];
1411   int nj;
1412 
1413   for (int i = 0; i < maxJets; ++i)
1414     used[i] = 0;
1415 
1416   //  cout << ">>>>>>>>> " << endl;
1417 
1418   for (GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd < maxJets; ++gen) {
1419     p4gen[jetInd] = gen->p4();  //Gen 4-vector
1420     dRmin[jetInd] = 1000.0;
1421 
1422     nj = 0;
1423     usedInd = -1;
1424 
1425     for (CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal) {
1426       double delR = deltaR(cal->eta(), cal->phi(), gen->eta(), gen->phi());
1427 
1428       if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
1429         dRmin[jetInd] = delR;       // delta R of match
1430         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
1431         usedInd = nj;
1432       }
1433 
1434       nj++;
1435     }
1436 
1437     if (fabs(p4gen[jetInd].eta()) < 1.3) {
1438       matchedAllPt11.Fill(p4gen[jetInd].pt());
1439     }
1440     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1441       matchedAllPt12.Fill(p4gen[jetInd].pt());
1442     }
1443     if (fabs(p4gen[jetInd].eta()) > 3.) {
1444       matchedAllPt13.Fill(p4gen[jetInd].pt());
1445     }
1446 
1447     if (usedInd != -1) {
1448       used[usedInd] = 1;
1449 
1450       if (p4cal[jetInd].pt() > minJetPt10)
1451         hf_PtResponse1.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt() / p4gen[jetInd].pt());
1452 
1453       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1454         matchedPt11.Fill(p4gen[jetInd].pt());
1455       }
1456       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1457         matchedPt12.Fill(p4gen[jetInd].pt());
1458       }
1459       if (fabs(p4gen[jetInd].eta()) > 3.) {
1460         matchedPt13.Fill(p4gen[jetInd].pt());
1461       }
1462       /***
1463       cout << " >>>SISCone "     
1464        << jetInd                 << " " 
1465        << p4gen[jetInd].pt()     << " " 
1466        << p4gen[jetInd].eta()    << " " 
1467        << p4gen[jetInd].phi()    << " " 
1468        << p4cal[jetInd].pt()     << " " 
1469        << p4cal[jetInd].eta()    << " " 
1470        << p4cal[jetInd].phi()    
1471        << endl;
1472       ***/
1473 
1474       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
1475       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1476         if ((p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40)) {
1477           dPt20Frac1.Fill(dpt / p4gen[jetInd].pt());
1478         }
1479         if ((p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60)) {
1480           dPt40Frac1.Fill(dpt / p4gen[jetInd].pt());
1481         }
1482         if ((p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100)) {
1483           dPt80Frac1.Fill(dpt / p4gen[jetInd].pt());
1484         }
1485         if ((p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120)) {
1486           dPt100Frac1.Fill(dpt / p4gen[jetInd].pt());
1487         }
1488       }
1489       if ((p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3)) {
1490         dR1.Fill(dRmin[jetInd]);
1491         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
1492         dPhi1.Fill(dphi);
1493         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
1494         dEta1.Fill(deta);
1495         dPt1.Fill(dpt);
1496         dPtFrac1.Fill(dpt / p4gen[jetInd].pt());
1497 
1498         if (((dpt / p4gen[jetInd].pt()) < -0.5) && (fabs(dpt) > 90.)) {
1499           cout << " deltaR min = " << dRmin[jetInd] << " Ind = " << jetInd << " / " << usedInd << " / " << used[nj]
1500                << " Del pt / frac  = " << dpt << " / " << dpt / p4gen[jetInd].pt()
1501                << " cal/gen pt   = " << p4cal[jetInd].pt() << " / " << p4gen[jetInd].pt()
1502                << " cal/gen eta  = " << p4cal[jetInd].eta() << " / " << p4gen[jetInd].eta()
1503                << " cal/gen phi  = " << p4cal[jetInd].phi() << " / " << p4gen[jetInd].phi() << endl;
1504         }
1505       }
1506 
1507       jetInd++;
1508     }
1509   }
1510 
1511   // *********************
1512   // Seedless Jet Matching
1513 
1514   //  Handle<GenJetCollection>  genJets;
1515   //  Handle<CaloJetCollection> caloJets;
1516 
1517   //  evt.getByLabel( "sisCone5GenJets",  genJets );
1518   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
1519   evt.getByLabel(GenJetAlgorithm2, genJets);
1520   evt.getByLabel(CaloJetAlgorithm2, caloJets);
1521 
1522   // int maxJets = 20;
1523   jetInd = 0;
1524   //  double dRmin[20];
1525   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
1526 
1527   for (int i = 0; i < maxJets; ++i)
1528     used[i] = 0;
1529   for (GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd < maxJets; ++gen) {
1530     p4gen[jetInd] = gen->p4();  //Gen 4-vector
1531     dRmin[jetInd] = 1000.0;
1532 
1533     nj = 0;
1534     usedInd = -1;
1535 
1536     for (CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal) {
1537       double delR = deltaR(cal->eta(), cal->phi(), gen->eta(), gen->phi());
1538 
1539       if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
1540         dRmin[jetInd] = delR;       // delta R of match
1541         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
1542         usedInd = nj;
1543       }
1544       nj++;
1545     }
1546 
1547     if (fabs(p4gen[jetInd].eta()) < 1.3) {
1548       matchedAllPt21.Fill(p4gen[jetInd].pt());
1549     }
1550     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1551       matchedAllPt22.Fill(p4gen[jetInd].pt());
1552     }
1553     if (fabs(p4gen[jetInd].eta()) > 3.) {
1554       matchedAllPt23.Fill(p4gen[jetInd].pt());
1555     }
1556 
1557     if (usedInd != -1) {
1558       used[usedInd] = 1;
1559 
1560       if (p4cal[jetInd].pt() > minJetPt10)
1561         hf_PtResponse2.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt() / p4gen[jetInd].pt());
1562 
1563       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1564         matchedPt21.Fill(p4gen[jetInd].pt());
1565       }
1566       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1567         matchedPt22.Fill(p4gen[jetInd].pt());
1568       }
1569       if (fabs(p4gen[jetInd].eta()) > 3.) {
1570         matchedPt23.Fill(p4gen[jetInd].pt());
1571       }
1572       /***
1573       cout << " >>>IterCone "     
1574        << jetInd                 << " " 
1575        << p4gen[jetInd].pt()     << " " 
1576        << p4gen[jetInd].eta()    << " " 
1577        << p4gen[jetInd].phi()    << " " 
1578        << p4cal[jetInd].pt()     << " " 
1579        << p4cal[jetInd].eta()    << " " 
1580        << p4cal[jetInd].phi()    
1581        << endl;
1582       ***/
1583 
1584       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
1585       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1586         if ((p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40)) {
1587           dPt20Frac2.Fill(dpt / p4gen[jetInd].pt());
1588         }
1589         if ((p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60)) {
1590           dPt40Frac2.Fill(dpt / p4gen[jetInd].pt());
1591         }
1592         if ((p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100)) {
1593           dPt80Frac2.Fill(dpt / p4gen[jetInd].pt());
1594         }
1595         if ((p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120)) {
1596           dPt100Frac2.Fill(dpt / p4gen[jetInd].pt());
1597         }
1598       }
1599       if ((p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3)) {
1600         dR2.Fill(dRmin[jetInd]);
1601         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
1602         dPhi2.Fill(dphi);
1603         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
1604         dEta2.Fill(deta);
1605         dPt2.Fill(dpt);
1606         dPtFrac2.Fill(dpt / p4gen[jetInd].pt());
1607       }
1608 
1609       jetInd++;
1610     }
1611   }
1612 
1613   // *********************
1614   // Kt Jet Matching
1615 
1616   //  Handle<GenJetCollection>  genJets;
1617   //  Handle<CaloJetCollection> caloJets;
1618 
1619   //  evt.getByLabel( "sisCone5GenJets",  genJets );
1620   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
1621   evt.getByLabel(GenJetAlgorithm3, genJets);
1622   evt.getByLabel(CaloJetAlgorithm3, caloJets);
1623 
1624   // int maxJets = 20;
1625   jetInd = 0;
1626   //  double dRmin[20];
1627   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
1628 
1629   for (int i = 0; i < maxJets; ++i)
1630     used[i] = 0;
1631   for (GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd < maxJets; ++gen) {
1632     p4gen[jetInd] = gen->p4();  //Gen 4-vector
1633     dRmin[jetInd] = 1000.0;
1634 
1635     nj = 0;
1636     usedInd = -1;
1637 
1638     for (CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal) {
1639       double delR = deltaR(cal->eta(), cal->phi(), gen->eta(), gen->phi());
1640 
1641       if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
1642         dRmin[jetInd] = delR;       // delta R of match
1643         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
1644         usedInd = nj;
1645       }
1646       nj++;
1647     }
1648 
1649     if (fabs(p4gen[jetInd].eta()) < 1.3) {
1650       matchedAllPt31.Fill(p4gen[jetInd].pt());
1651     }
1652     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1653       matchedAllPt32.Fill(p4gen[jetInd].pt());
1654     }
1655     if (fabs(p4gen[jetInd].eta()) > 3.) {
1656       matchedAllPt33.Fill(p4gen[jetInd].pt());
1657     }
1658 
1659     if (usedInd != -1) {
1660       used[usedInd] = 1;
1661 
1662       if (p4cal[jetInd].pt() > minJetPt10)
1663         hf_PtResponse3.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt() / p4gen[jetInd].pt());
1664 
1665       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1666         matchedPt31.Fill(p4gen[jetInd].pt());
1667       }
1668       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1669         matchedPt32.Fill(p4gen[jetInd].pt());
1670       }
1671       if (fabs(p4gen[jetInd].eta()) > 3.) {
1672         matchedPt33.Fill(p4gen[jetInd].pt());
1673       }
1674       /***
1675       cout << " >>>MidPoint "     
1676        << jetInd                 << " " 
1677        << p4gen[jetInd].pt()     << " " 
1678        << p4gen[jetInd].eta()    << " " 
1679        << p4gen[jetInd].phi()    << " " 
1680        << p4cal[jetInd].pt()     << " " 
1681        << p4cal[jetInd].eta()    << " " 
1682        << p4cal[jetInd].phi()    
1683        << endl;
1684       ***/
1685 
1686       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
1687       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1688         if ((p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40)) {
1689           dPt20Frac3.Fill(dpt / p4gen[jetInd].pt());
1690         }
1691         if ((p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60)) {
1692           dPt40Frac3.Fill(dpt / p4gen[jetInd].pt());
1693         }
1694         if ((p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100)) {
1695           dPt80Frac3.Fill(dpt / p4gen[jetInd].pt());
1696         }
1697         if ((p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120)) {
1698           dPt100Frac3.Fill(dpt / p4gen[jetInd].pt());
1699         }
1700       }
1701 
1702       if ((p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3)) {
1703         dR3.Fill(dRmin[jetInd]);
1704         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
1705         dPhi3.Fill(dphi);
1706         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
1707         dEta3.Fill(deta);
1708         dPt3.Fill(dpt);
1709         dPtFrac3.Fill(dpt / p4gen[jetInd].pt());
1710       }
1711 
1712       jetInd++;
1713     }
1714   }
1715 
1716   // *********************
1717   // Jet Algorithm 4
1718 
1719   evt.getByLabel(GenJetAlgorithm4, genJets);
1720   evt.getByLabel(CaloJetAlgorithm4, caloJets);
1721 
1722   // int maxJets = 20;
1723   jetInd = 0;
1724   //  double dRmin[20];
1725   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
1726 
1727   for (int i = 0; i < maxJets; ++i)
1728     used[i] = 0;
1729   for (GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && jetInd < maxJets; ++gen) {
1730     p4gen[jetInd] = gen->p4();  //Gen 4-vector
1731     dRmin[jetInd] = 1000.0;
1732 
1733     nj = 0;
1734     usedInd = -1;
1735 
1736     for (CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal) {
1737       double delR = deltaR(cal->eta(), cal->phi(), gen->eta(), gen->phi());
1738 
1739       if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
1740         dRmin[jetInd] = delR;       // delta R of match
1741         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
1742         usedInd = nj;
1743       }
1744       nj++;
1745     }
1746 
1747     if (fabs(p4gen[jetInd].eta()) < 1.3) {
1748       matchedAllPt41.Fill(p4gen[jetInd].pt());
1749     }
1750     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1751       matchedAllPt42.Fill(p4gen[jetInd].pt());
1752     }
1753     if (fabs(p4gen[jetInd].eta()) > 3.) {
1754       matchedAllPt43.Fill(p4gen[jetInd].pt());
1755     }
1756 
1757     if (usedInd != -1) {
1758       used[usedInd] = 1;
1759 
1760       if (p4cal[jetInd].pt() > minJetPt10)
1761         hf_PtResponse4.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt() / p4gen[jetInd].pt());
1762 
1763       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1764         matchedPt41.Fill(p4gen[jetInd].pt());
1765       }
1766       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.)) {
1767         matchedPt42.Fill(p4gen[jetInd].pt());
1768       }
1769       if (fabs(p4gen[jetInd].eta()) > 3.) {
1770         matchedPt43.Fill(p4gen[jetInd].pt());
1771       }
1772 
1773       /***
1774       cout << " >>>KtClus "     
1775        << jetInd                 << " " 
1776        << p4gen[jetInd].pt()     << " " 
1777        << p4gen[jetInd].eta()    << " " 
1778        << p4gen[jetInd].phi()    << " " 
1779        << p4cal[jetInd].pt()     << " " 
1780        << p4cal[jetInd].eta()    << " " 
1781        << p4cal[jetInd].phi()    
1782        << endl;
1783       ***/
1784 
1785       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
1786       if (fabs(p4gen[jetInd].eta()) < 1.3) {
1787         if ((p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40)) {
1788           dPt20Frac4.Fill(dpt / p4gen[jetInd].pt());
1789         }
1790         if ((p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60)) {
1791           dPt40Frac4.Fill(dpt / p4gen[jetInd].pt());
1792         }
1793         if ((p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100)) {
1794           dPt80Frac4.Fill(dpt / p4gen[jetInd].pt());
1795         }
1796         if ((p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120)) {
1797           dPt100Frac4.Fill(dpt / p4gen[jetInd].pt());
1798         }
1799       }
1800 
1801       if ((p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3)) {
1802         dR4.Fill(dRmin[jetInd]);
1803         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
1804         dPhi4.Fill(dphi);
1805         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
1806         dEta4.Fill(deta);
1807         dPt4.Fill(dpt);
1808         dPtFrac4.Fill(dpt / p4gen[jetInd].pt());
1809       }
1810 
1811       jetInd++;
1812     }
1813   }
1814 
1815   // *********************
1816   // MidPoint - Seedless Jet Matching
1817 
1818   Handle<CaloJetCollection> calo1Jets;
1819   Handle<CaloJetCollection> calo2Jets;
1820   Handle<CaloJetCollection> calo3Jets;
1821 
1822   evt.getByLabel(CaloJetAlgorithm1, calo1Jets);
1823   evt.getByLabel(CaloJetAlgorithm2, calo2Jets);
1824   evt.getByLabel(CaloJetAlgorithm3, calo3Jets);
1825 
1826   jetInd = 0;
1827 
1828   for (int i = 0; i < maxJets; ++i)
1829     used[i] = 0;
1830   for (CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); cal1 != calo1Jets->end() && jetInd < maxJets;
1831        ++cal1) {
1832     p4gen[jetInd] = cal1->p4();  //Gen 4-vector
1833     dRmin[jetInd] = 1000.0;
1834 
1835     nj = 0;
1836     for (CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); cal2 != calo2Jets->end(); ++cal2) {
1837       double delR = deltaR(cal1->eta(), cal1->phi(), cal2->eta(), cal2->phi());
1838       if ((delR < dRmin[jetInd]) && (used[nj] == 0)) {
1839         dRmin[jetInd] = delR;        // delta R of match
1840         p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
1841         usedInd = nj;
1842       }
1843       nj++;
1844     }
1845     used[usedInd] = 1;
1846 
1847     if (p4gen[jetInd].pt() > minJetPt) {
1848       dR12.Fill(dRmin[jetInd]);
1849       double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
1850       dPhi12.Fill(dphi);
1851       double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
1852       dEta12.Fill(deta);
1853       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
1854       dPt12.Fill(dpt);
1855     }
1856 
1857     jetInd++;
1858   }
1859 
1860   // ******************************************
1861   // ******************************************
1862 
1863   Handle<CandidateCollection> genParticles;
1864   evt.getByLabel("genParticleCandidates", genParticles);
1865 
1866   // *********************
1867   // Partons (Z')
1868 
1869   int nPart = 0;
1870   for (size_t i = 0; i < genParticles->size(); i++) {
1871     const Candidate& p = (*genParticles)[i];
1872     //    int Status =  p.status();
1873     //    bool ParticleIsStable = Status==1;
1874     int id = p.pdgId();
1875 
1876     if (id == 32) {
1877       if (p.numberOfDaughters() != 0) {
1878         /***
1879     cout << "Z': part = "   << i << " id = " << id 
1880          << " daughters = " << p.numberOfDaughters() 
1881          << " mass = "      << p.mass()
1882          << endl;
1883     ***/
1884         ZpMG = p.mass();
1885         ZpMassGen.Fill(ZpMG);
1886         if (EtaOk10 == 2) {
1887           ZpMassGen10.Fill(ZpMG);
1888           ZpMassGen_700_10.Fill(ZpMG);
1889           ZpMassGen_2000_10.Fill(ZpMG);
1890           ZpMassGen_5000_10.Fill(ZpMG);
1891         }
1892         if (EtaOk13 == 2) {
1893           ZpMassGen13.Fill(ZpMG);
1894           ZpMassGen_700_13.Fill(ZpMG);
1895           ZpMassGen_2000_13.Fill(ZpMG);
1896           ZpMassGen_5000_13.Fill(ZpMG);
1897         }
1898         if (EtaOk40 == 2) {
1899           ZpMassGen40.Fill(ZpMG);
1900           ZpMassGen_700_40.Fill(ZpMG);
1901           ZpMassGen_2000_40.Fill(ZpMG);
1902           ZpMassGen_5000_40.Fill(ZpMG);
1903         }
1904       }
1905 
1906       for (int id1 = 0, nd1 = p.numberOfDaughters(); id1 < nd1; ++id1) {
1907         const Candidate* d1 = p.daughter(id1);
1908 
1909         if (abs(d1->pdgId()) != 32) {
1910           math::XYZTLorentzVector momentum = d1->p4();
1911           p4Zp[nPart] = momentum = d1->p4();
1912           nPart++;
1913         }
1914       }
1915     }
1916   }
1917 
1918   // *********************
1919   // Match jets to Zp
1920   int genInd = 0;
1921 
1922   if (nPart == 2) {
1923     ZpM = (p4Zp[0] + p4Zp[1]).mass();
1924     ZpMass.Fill(ZpM);
1925 
1926     if (EtaOk10 == 2) {
1927       ZpMass_700_10.Fill(ZpM);
1928       ZpMass_2000_10.Fill(ZpM);
1929       ZpMass_5000_10.Fill(ZpM);
1930     }
1931     if (EtaOk13 == 2) {
1932       ZpMass_700_13.Fill(ZpM);
1933       ZpMass_2000_13.Fill(ZpM);
1934       ZpMass_5000_13.Fill(ZpM);
1935     }
1936     if (EtaOk40 == 2) {
1937       ZpMass_700_40.Fill(ZpM);
1938       ZpMass_2000_40.Fill(ZpM);
1939       ZpMass_5000_40.Fill(ZpM);
1940     }
1941 
1942     int usedInd;
1943 
1944     // ***********
1945     // **** Calor1
1946     usedInd = -1;
1947     jetInd = 0;
1948 
1949     for (int i = 0; i < maxJets; ++i)
1950       used[i] = 0;
1951     for (int i = 0; i < 2; ++i) {
1952       dRmin[jetInd] = 1000.0;
1953 
1954       int nj = 0;
1955       for (CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); cal1 != calo1Jets->end() && jetInd < maxJets;
1956            ++cal1) {
1957         double delR = deltaR(cal1->eta(), cal1->phi(), p4Zp[i].eta(), p4Zp[i].phi());
1958 
1959         //  if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
1960         if ((delR < dRmin[jetInd]) && (used[nj] == 0)) {
1961           dRmin[jetInd] = delR;        // delta R of match
1962           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
1963           usedInd = nj;
1964           genInd = i;
1965         }
1966 
1967         /****
1968     cout << "Delta R = " << delR
1969          << " deltaR min = " << dRmin[jetInd]
1970          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
1971          << " cal1 eta  = " << cal1->eta()
1972          << " p4par eta = " << p4Zp[i].eta()
1973          << " cal1 phi  = " << cal1->phi()
1974          << " p4par phi = " << p4Zp[i].phi()
1975          << endl;
1976     cout << "    " 
1977          << " p4par = " << p4Zp[i]
1978          << " p4cal = " << cal1->p4()
1979          << endl;
1980     ***/
1981 
1982         nj++;
1983       }
1984 
1985       // Found matched jet
1986       if (usedInd != -1) {
1987         used[usedInd] = 1;
1988         jetInd++;
1989       }
1990     }
1991 
1992     ZpMM = (p4cal[0] + p4cal[1]).mass();
1993     ZpMassMatched1.Fill(ZpMM);
1994 
1995     if ((ZpMG != 0) && (EtaOk40 == 2)) {
1996       ZpMassRes401.Fill((ZpMM - ZpMG) / ZpMG);
1997 
1998       ZpMassResL401.Fill((LeadMass1 - ZpMG) / ZpMG);
1999       ZpMassResL402.Fill((LeadMass2 - ZpMG) / ZpMG);
2000       ZpMassResL403.Fill((LeadMass3 - ZpMG) / ZpMG);
2001 
2002       ZpMassResRL401.Fill(LeadMass1 / ZpMG);
2003       ZpMassResRL402.Fill(LeadMass2 / ZpMG);
2004       ZpMassResRL403.Fill(LeadMass3 / ZpMG);
2005 
2006       ZpMassResRLoP401.Fill(LeadMass1 / LeadMassP1);
2007       ZpMassResRLoP402.Fill(LeadMass2 / LeadMassP2);
2008       ZpMassResRLoP403.Fill(LeadMass3 / LeadMassP2);
2009 
2010       ZpMassResPRL401.Fill(LeadMassP1 / ZpMG);
2011       ZpMassResPRL402.Fill(LeadMassP2 / ZpMG);
2012       ZpMassResPRL403.Fill(LeadMassP3 / ZpMG);
2013     }
2014 
2015     if ((ZpMG != 0) && (EtaOk10 == 2)) {
2016       ZpMassRes101.Fill((ZpMM - ZpMG) / ZpMG);
2017 
2018       ZpMassResL101.Fill((LeadMass1 - ZpMG) / ZpMG);
2019       ZpMassResL102.Fill((LeadMass2 - ZpMG) / ZpMG);
2020       ZpMassResL103.Fill((LeadMass3 - ZpMG) / ZpMG);
2021 
2022       ZpMassResRL101.Fill(LeadMass1 / ZpMG);
2023       ZpMassResRL102.Fill(LeadMass2 / ZpMG);
2024       ZpMassResRL103.Fill(LeadMass3 / ZpMG);
2025 
2026       ZpMassResRLoP101.Fill(LeadMass1 / LeadMassP1);
2027       ZpMassResRLoP102.Fill(LeadMass2 / LeadMassP2);
2028       ZpMassResRLoP103.Fill(LeadMass3 / LeadMassP2);
2029 
2030       ZpMassResPRL101.Fill(LeadMassP1 / ZpMG);
2031       ZpMassResPRL102.Fill(LeadMassP2 / ZpMG);
2032       ZpMassResPRL103.Fill(LeadMassP3 / ZpMG);
2033     }
2034 
2035     if ((ZpMG != 0) && (EtaOk13 == 2)) {
2036       ZpMassRes131.Fill((ZpMM - ZpMG) / ZpMG);
2037 
2038       ZpMassResL131.Fill((LeadMass1 - ZpMG) / ZpMG);
2039       ZpMassResL132.Fill((LeadMass2 - ZpMG) / ZpMG);
2040       ZpMassResL133.Fill((LeadMass3 - ZpMG) / ZpMG);
2041 
2042       ZpMassResRL131.Fill(LeadMass1 / ZpMG);
2043       ZpMassResRL132.Fill(LeadMass2 / ZpMG);
2044       ZpMassResRL133.Fill(LeadMass3 / ZpMG);
2045 
2046       ZpMassResRLoP131.Fill(LeadMass1 / LeadMassP1);
2047       ZpMassResRLoP132.Fill(LeadMass2 / LeadMassP2);
2048       ZpMassResRLoP133.Fill(LeadMass3 / LeadMassP2);
2049 
2050       ZpMassResPRL131.Fill(LeadMassP1 / ZpMG);
2051       ZpMassResPRL132.Fill(LeadMassP2 / ZpMG);
2052       ZpMassResPRL133.Fill(LeadMassP3 / ZpMG);
2053     }
2054 
2055     // ***********
2056     // **** Calor2
2057     usedInd = -1;
2058     jetInd = 0;
2059 
2060     for (int i = 0; i < maxJets; ++i)
2061       used[i] = 0;
2062     for (int i = 0; i < 2; ++i) {
2063       dRmin[jetInd] = 1000.0;
2064 
2065       int nj = 0;
2066       for (CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); cal2 != calo2Jets->end() && jetInd < maxJets;
2067            ++cal2) {
2068         double delR = deltaR(cal2->eta(), cal2->phi(), p4Zp[i].eta(), p4Zp[i].phi());
2069 
2070         if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
2071           dRmin[jetInd] = delR;        // delta R of match
2072           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
2073           usedInd = nj;
2074         }
2075 
2076         /****   
2077     cout << "Delta R = " << delR
2078          << " deltaR min = " << dRmin[jetInd]
2079          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
2080          << " p4par = " << p4par[i]
2081          << " p4cal = " << cal1->p4()
2082          << " cal1 eta  = " << cal1->eta()
2083          << " p4par eta = " << p4par[i].eta()
2084          << endl;
2085     ****/
2086 
2087         nj++;
2088       }
2089 
2090       // Found matched jet
2091       if (usedInd != -1) {
2092         used[usedInd] = 1;
2093         jetInd++;
2094       }
2095     }
2096 
2097     ZpMM = (p4cal[0] + p4cal[1]).mass();
2098     ZpMassMatched2.Fill(ZpMM);
2099     ZpMassRes402.Fill((ZpMM - ZpM) / ZpM);
2100 
2101     // ***********
2102     // **** Calor3
2103     usedInd = -1;
2104     jetInd = 0;
2105 
2106     for (int i = 0; i < maxJets; ++i)
2107       used[i] = 0;
2108     for (int i = 0; i < 2; ++i) {
2109       dRmin[jetInd] = 1000.0;
2110 
2111       int nj = 0;
2112       for (CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); cal3 != calo3Jets->end() && jetInd < maxJets;
2113            ++cal3) {
2114         double delR = deltaR(cal3->eta(), cal3->phi(), p4Zp[i].eta(), p4Zp[i].phi());
2115 
2116         if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
2117           dRmin[jetInd] = delR;        // delta R of match
2118           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
2119           usedInd = nj;
2120         }
2121 
2122         /****   
2123     cout << "Delta R = " << delR
2124          << " deltaR min = " << dRmin[jetInd]
2125          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
2126          << " p4par = " << p4par[i]
2127          << " p4cal = " << cal1->p4()
2128          << " cal1 eta  = " << cal1->eta()
2129          << " p4par eta = " << p4par[i].eta()
2130          << endl;
2131     ****/
2132 
2133         nj++;
2134       }
2135 
2136       // Found matched jet
2137       if (usedInd != -1) {
2138         used[usedInd] = 1;
2139         jetInd++;
2140       }
2141     }
2142 
2143     ZpMM = (p4cal[0] + p4cal[1]).mass();
2144     ZpMassMatched3.Fill(ZpMM);
2145     ZpMassRes403.Fill((ZpMM - ZpM) / ZpM);
2146 
2147   } else {
2148     cout << "Z' (3): nPart = " << nPart << endl;
2149   }
2150 
2151   // *********************
2152   // Partons (ttbar) Jet Matching
2153 
2154   //  cout << ">>> Begin MC list" << endl;
2155   int nJet = 0;
2156 
2157   int ii = 1;
2158   int jj = 4;
2159 
2160   for (size_t i = 0; i < genParticles->size(); i++) {
2161     const Candidate& p = (*genParticles)[i];
2162     //    int Status =  p.status();
2163     //    bool ParticleIsStable = Status==1;
2164     int id = p.pdgId();
2165 
2166     // Top Quarks
2167     if (abs(id) == 6) {
2168       cout << "TOP: id = " << id << " mass = " << p.mass() << endl;
2169 
2170       topMassParton.Fill(p.mass());
2171 
2172       if (id == 6)
2173         tMassGen.Fill(p.mass());
2174       if (id == -6)
2175         tbarMassGen.Fill(p.mass());
2176 
2177       for (size_t id1 = 0, nd1 = p.numberOfDaughters(); id1 < nd1; ++id1) {
2178         const Candidate* d1 = p.daughter(id1);
2179 
2180         // b - quark
2181         if (abs(d1->pdgId()) == 5) {
2182           math::XYZTLorentzVector momentum = d1->p4();
2183           p4par[nJet++] = momentum = d1->p4();
2184 
2185           cout << "Daughter1: id = " << d1->pdgId() << " daughters = " << d1->numberOfDaughters()
2186                << " mother 1   = " << (d1->mother())->pdgId() << " Momentum " << momentum << " GeV/c" << endl;
2187 
2188           if ((d1->mother())->pdgId() == 6) {
2189             p4part[0] = momentum = d1->p4();
2190             cout << ">>> part0 = " << p4part[0] << endl;
2191           }
2192           if ((d1->mother())->pdgId() == -6) {
2193             p4part[3] = momentum = d1->p4();
2194             cout << ">>> part3 = " << p4part[3] << endl;
2195           }
2196         }
2197 
2198         // W
2199         // Check for fully hadronic decay
2200 
2201         if (abs(d1->pdgId()) == 24) {
2202           for (size_t id2 = 0, nd2 = d1->numberOfDaughters(); id2 < nd2; ++id2) {
2203             const Candidate* d2 = d1->daughter(id2);
2204 
2205             if (abs(d2->pdgId()) < 9) {
2206               math::XYZVector vertex(d2->vx(), d2->vy(), d2->vz());
2207               math::XYZTLorentzVector momentum = d2->p4();
2208               p4par[nJet++] = momentum = d2->p4();
2209 
2210               if ((d1->mother())->pdgId() == 6) {
2211                 p4part[ii] = momentum = d2->p4();
2212                 cout << ">>> part" << ii << " = " << p4part[ii] << endl;
2213                 ii++;
2214               }
2215               if ((d1->mother())->pdgId() == -6) {
2216                 p4part[jj] = momentum = d2->p4();
2217                 cout << ">>> part" << jj << " = " << p4part[jj] << endl;
2218                 jj++;
2219               }
2220 
2221               cout << "Daughter2: id = " << d2->pdgId() << " daughters = " << d2->numberOfDaughters()
2222                    << " mother 2   = " << (d2->mother())->pdgId() << " Momentum " << momentum << " GeV/c" << endl;
2223             }
2224           }
2225         }
2226 
2227         //  if ( pdgId == d->pdgId() && d->status() == 1 ) {
2228         //  }
2229       }
2230     }
2231   }
2232   //  cout << ">>> N Jets = " << nJet << endl;
2233 
2234   if (nJet == 6) {
2235     double tmass = (p4part[0] + p4part[1] + p4part[2]).mass();
2236     double tbarmass = (p4part[3] + p4part[4] + p4part[5]).mass();
2237 
2238     tMass.Fill(tmass);
2239     tbarMass.Fill(tbarmass);
2240 
2241     cout << ">>> T Mass = " << tmass << " / " << tbarmass << endl;
2242 
2243     double mindR = 1000.;
2244     for (size_t i = 0; i < 6; ++i) {
2245       for (size_t j = 0; j < 6; ++j) {
2246         if (j > i) {
2247           double delR = deltaR(p4par[i].eta(), p4par[i].phi(), p4par[j].eta(), p4par[j].phi());
2248           if (delR < mindR)
2249             mindR = delR;
2250           dRParton.Fill(delR);
2251         }
2252       }
2253     }
2254     dRPartonMin.Fill(mindR);
2255 
2256     int usedInd;
2257     usedInd = -1;
2258     jetInd = 0;
2259 
2260     for (int i = 0; i < maxJets; ++i)
2261       used[i] = 0;
2262     for (int i = 0; i < 6; ++i) {
2263       dRmin[jetInd] = 1000.0;
2264 
2265       int nj = 0;
2266       for (CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); cal1 != calo1Jets->end() && jetInd < maxJets;
2267            ++cal1) {
2268         double delR = deltaR(cal1->eta(), cal1->phi(), p4par[i].eta(), p4par[i].phi());
2269 
2270         if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
2271           dRmin[jetInd] = delR;        // delta R of match
2272           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
2273           usedInd = nj;
2274           genInd = i;
2275         }
2276 
2277         /****   
2278     cout << "Delta R = " << delR
2279          << " deltaR min = " << dRmin[jetInd]
2280          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
2281          << " p4par = " << p4par[i]
2282          << " p4cal = " << cal1->p4()
2283          << " cal1 eta  = " << cal1->eta()
2284          << " p4par eta = " << p4par[i].eta()
2285          << endl;
2286     ****/
2287 
2288         nj++;
2289       }
2290 
2291       // Found matched jet
2292       if (usedInd != -1) {
2293         used[usedInd] = 1;
2294 
2295         dRPar1.Fill(dRmin[jetInd]);
2296         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
2297         dPhiPar1.Fill(dphi);
2298         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
2299         dEtaPar1.Fill(deta);
2300         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
2301         dPtPar1.Fill(dpt);
2302         jetInd++;
2303       }
2304     }
2305     ParMatch1.Fill(jetInd);
2306     if (jetInd == 6) {
2307       topMass1.Fill((p4cal[0] + p4cal[1] + p4cal[2]).mass());
2308       topMass1.Fill((p4cal[3] + p4cal[4] + p4cal[5]).mass());
2309     }
2310 
2311     /***
2312     cout << "Collection Size = " <<  calo1Jets->size() 
2313      << " / " << jetInd
2314      << endl;
2315     ***/
2316 
2317     // ***********************
2318     jetInd = 0;
2319     usedInd = -1;
2320 
2321     for (int i = 0; i < maxJets; ++i)
2322       used[i] = 0;
2323     for (int i = 0; i < 6; ++i) {
2324       dRmin[jetInd] = 1000.0;
2325 
2326       int nj = 0;
2327       for (CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); cal2 != calo2Jets->end() && jetInd < maxJets;
2328            ++cal2) {
2329         double delR = deltaR(cal2->eta(), cal2->phi(), p4par[i].eta(), p4par[i].phi());
2330 
2331         if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
2332           dRmin[jetInd] = delR;        // delta R of match
2333           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
2334           usedInd = nj;
2335           genInd = i;
2336         }
2337 
2338         /****
2339     cout << "Delta R = " << delR
2340          << " deltaR min = " << dRmin[jetInd]
2341          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
2342          << " cal2 eta  = " << cal2->eta()
2343          << " p4par eta = " << p4par[i].eta()
2344          << endl;
2345     ****/
2346 
2347         nj++;
2348       }
2349       if (usedInd != -1) {
2350         used[usedInd] = 1;
2351 
2352         dRPar2.Fill(dRmin[jetInd]);
2353         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
2354         dPhiPar2.Fill(dphi);
2355         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
2356         dEtaPar2.Fill(deta);
2357         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
2358         dPtPar2.Fill(dpt);
2359 
2360         jetInd++;
2361       }
2362     }
2363     ParMatch2.Fill(jetInd);
2364     if (jetInd == 6) {
2365       topMass2.Fill((p4cal[0] + p4cal[1] + p4cal[2]).mass());
2366       topMass2.Fill((p4cal[3] + p4cal[4] + p4cal[5]).mass());
2367     }
2368 
2369     /***
2370     cout << "Collection Size = " <<  calo2Jets->size() 
2371      << " / " << jetInd
2372      << endl;
2373     ***/
2374 
2375     // ***********************
2376     jetInd = 0;
2377     usedInd = -1;
2378 
2379     for (int i = 0; i < maxJets; ++i)
2380       used[i] = 0;
2381     for (int i = 0; i < 6; ++i) {
2382       dRmin[jetInd] = 1000.0;
2383 
2384       int nj = 0;
2385       for (CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); cal3 != calo3Jets->end() && jetInd < maxJets;
2386            ++cal3) {
2387         double delR = deltaR(cal3->eta(), cal3->phi(), p4par[i].eta(), p4par[i].phi());
2388 
2389         if ((delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0)) {
2390           dRmin[jetInd] = delR;        // delta R of match
2391           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
2392           usedInd = nj;
2393           genInd = i;
2394         }
2395         /****
2396     cout << "Delta R = " << delR
2397          << " deltaR min = " << dRmin[jetInd]
2398          << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
2399          << " cal3 eta  = " << cal3->eta()
2400          << " p4par eta = " << p4par[i].eta()
2401          << endl;
2402     ****/
2403 
2404         nj++;
2405       }
2406       if (usedInd != -1) {
2407         used[usedInd] = 1;
2408 
2409         dRPar3.Fill(dRmin[jetInd]);
2410         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
2411         dPhiPar3.Fill(dphi);
2412         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
2413         dEtaPar3.Fill(deta);
2414         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
2415         dPtPar3.Fill(dpt);
2416 
2417         jetInd++;
2418       }
2419     }
2420     ParMatch3.Fill(jetInd);
2421     if (jetInd == 6) {
2422       topMass3.Fill((p4cal[0] + p4cal[1] + p4cal[2]).mass());
2423       topMass3.Fill((p4cal[3] + p4cal[4] + p4cal[5]).mass());
2424     }
2425 
2426     /***
2427     cout << "Collection Size = " <<  calo3Jets->size() 
2428      << " / " << jetInd
2429      << endl;
2430     ***/
2431   }
2432 
2433   int nTow1, nTow2, nTow3, nTow4;
2434 
2435   Handle<CaloJetCollection> jets;
2436 
2437   // *********************
2438   // Jet Properties
2439   // *********************
2440 
2441   // --- Loop over jets and make a list of all the used towers
2442   evt.getByLabel(CaloJetAlgorithm1, jets);
2443   int jjet = 0;
2444   for (CaloJetCollection::const_iterator ijet = jets->begin(); ijet != jets->end(); ijet++) {
2445     jjet++;
2446 
2447     float hadEne = ijet->hadEnergyInHB() + ijet->hadEnergyInHO() + ijet->hadEnergyInHE() + ijet->hadEnergyInHF();
2448     float emEne = ijet->emEnergyInEB() + ijet->emEnergyInEE() + ijet->emEnergyInHF();
2449     float had = ijet->energyFractionHadronic();
2450 
2451     float j_et = ijet->et();
2452 
2453     if (fabs(ijet->eta()) < 1.3) {
2454       totEneLeadJetEta1_1.Fill(hadEne + emEne);
2455       hadEneLeadJetEta1_1.Fill(hadEne);
2456       emEneLeadJetEta1_1.Fill(emEne);
2457 
2458       totEneLeadJetEta1_2.Fill(hadEne + emEne);
2459       hadEneLeadJetEta1_2.Fill(hadEne);
2460       emEneLeadJetEta1_2.Fill(emEne);
2461 
2462       if (ijet->pt() > minJetPt10)
2463         hadFracEta11.Fill(had);
2464     }
2465     if ((fabs(ijet->eta()) > 1.3) && (fabs(ijet->eta()) < 3.)) {
2466       totEneLeadJetEta2_1.Fill(hadEne + emEne);
2467       hadEneLeadJetEta2_1.Fill(hadEne);
2468       emEneLeadJetEta2_1.Fill(emEne);
2469 
2470       totEneLeadJetEta2_2.Fill(hadEne + emEne);
2471       hadEneLeadJetEta2_2.Fill(hadEne);
2472       emEneLeadJetEta2_2.Fill(emEne);
2473 
2474       if (ijet->pt() > minJetPt10)
2475         hadFracEta21.Fill(had);
2476     }
2477     if (fabs(ijet->eta()) > 3.) {
2478       totEneLeadJetEta3_1.Fill(hadEne + emEne);
2479       hadEneLeadJetEta3_1.Fill(hadEne);
2480       emEneLeadJetEta3_1.Fill(emEne);
2481 
2482       totEneLeadJetEta3_2.Fill(hadEne + emEne);
2483       hadEneLeadJetEta3_2.Fill(hadEne);
2484       emEneLeadJetEta3_2.Fill(emEne);
2485 
2486       if (ijet->pt() > minJetPt10)
2487         hadFracEta31.Fill(had);
2488     }
2489 
2490     if (jjet == 1) {
2491       hadFracLeadJet1.Fill(had);
2492       hadEneLeadJet1.Fill(hadEne);
2493       hadEneLeadJet12.Fill(hadEne);
2494       hadEneLeadJet13.Fill(hadEne);
2495       emEneLeadJet1.Fill(emEne);
2496       emEneLeadJet12.Fill(emEne);
2497       emEneLeadJet13.Fill(emEne);
2498     }
2499 
2500     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
2501     int nConstituents = jetCaloRefs.size();
2502 
2503     if (jjet == 1) {
2504       nTow1 = nTow2 = nTow3 = nTow4 = 0;
2505       for (int i = 0; i < nConstituents; i++) {
2506         float et = jetCaloRefs[i]->et();
2507 
2508         if (et > 0.5)
2509           nTow1++;
2510         if (et > 1.0)
2511           nTow2++;
2512         if (et > 1.5)
2513           nTow3++;
2514         if (et > 2.0)
2515           nTow4++;
2516 
2517         hf_TowerJetEt1.Fill(et / j_et);
2518       }
2519 
2520       nTowersLeadJetPt1.Fill(nTow1);
2521       nTowersLeadJetPt2.Fill(nTow2);
2522       nTowersLeadJetPt3.Fill(nTow3);
2523       nTowersLeadJetPt4.Fill(nTow4);
2524     }
2525 
2526     if ((jjet == 1) && (fabs(ijet->eta()) < 1.3)) {
2527       nTowersLeadJet1.Fill(nConstituents);
2528 
2529       for (int i = 0; i < nConstituents; i++) {
2530         float t_et = jetCaloRefs[i]->et();
2531         double delR = deltaR(ijet->eta(), ijet->phi(), jetCaloRefs[i]->eta(), jetCaloRefs[i]->phi());
2532         hf_TowerDelR1.Fill(delR, t_et / j_et);
2533         hf_TowerDelR12.Fill(delR, t_et / j_et);
2534         TowerEtLeadJet1.Fill(t_et);
2535         TowerEtLeadJet12.Fill(t_et);
2536         TowerEtLeadJet13.Fill(t_et);
2537       }
2538     }
2539 
2540     if ((jjet == 2) && (fabs(ijet->eta()) < 1.3)) {
2541       nTowersSecondJet1.Fill(nConstituents);
2542     }
2543   }
2544 
2545   // *********************
2546   // Unclustered Energy
2547   // *********************
2548 
2549   double SumEtNotJets(0);
2550   double SumEtJets(0);
2551 
2552   double sumTowerAllEx(0);
2553   double sumTowerAllEy(0);
2554 
2555   std::vector<CaloTowerPtr> UsedTowerList;
2556   std::vector<CaloTower> TowerUsedInJets;
2557   std::vector<CaloTower> TowerNotUsedInJets;
2558 
2559   // *********************
2560   // Towers
2561 
2562   Handle<CaloTowerCollection> caloTowers;
2563   evt.getByLabel("towerMaker", caloTowers);
2564 
2565   nTow1 = nTow2 = nTow3 = nTow4 = 0;
2566 
2567   double sum_et = 0.0;
2568   double sum_ex = 0.0;
2569   double sum_ey = 0.0;
2570   //  double sum_ez = 0.0;
2571 
2572   // --- Loop over towers and make a lists of used and unused towers
2573   for (CaloTowerCollection::const_iterator tower = caloTowers->begin(); tower != caloTowers->end(); tower++) {
2574     Double_t et = tower->et();
2575 
2576     if (et > 0.5)
2577       nTow1++;
2578     if (et > 1.0)
2579       nTow2++;
2580     if (et > 1.5)
2581       nTow3++;
2582     if (et > 2.0)
2583       nTow4++;
2584 
2585     if (et > 0.5) {
2586       // ********
2587       double phix = tower->phi();
2588       //      double theta = tower->theta();
2589       //      double e     = tower->energy();
2590       //      double et    = e*sin(theta);
2591       //      double et    = tower->emEt() + tower->hadEt();
2592       double et = tower->et();
2593 
2594       //      sum_ez += e*cos(theta);
2595       sum_et += et;
2596       sum_ex += et * cos(phix);
2597       sum_ey += et * sin(phix);
2598       // ********
2599 
2600       Double_t phi = tower->phi();
2601 
2602       sumTowerAllEx += et * cos(phi);
2603       sumTowerAllEy += et * sin(phi);
2604     }
2605   }
2606 
2607   SumEt1.Fill(sum_et);
2608   SumEt12.Fill(sum_et);
2609   SumEt13.Fill(sum_et);
2610 
2611   MET1.Fill(sqrt(sum_ex * sum_ex + sum_ey * sum_ey));
2612   MET12.Fill(sqrt(sum_ex * sum_ex + sum_ey * sum_ey));
2613   MET13.Fill(sqrt(sum_ex * sum_ex + sum_ey * sum_ey));
2614 
2615   //  met->mex   = -sum_ex;
2616   //  met->mey   = -sum_ey;
2617   //  met->mez   = -sum_ez;
2618   //  met->met   = sqrt( sum_ex*sum_ex + sum_ey*sum_ey );
2619   // cout << "MET = " << met->met << endl;
2620   //  met->sumet = sum_et;
2621   //  met->phi   = atan2( -sum_ey, -sum_ex );
2622 
2623   hf_sumTowerAllEx.Fill(sumTowerAllEx);
2624   hf_sumTowerAllEy.Fill(sumTowerAllEy);
2625 
2626   nTowers1.Fill(nTow1);
2627   nTowers2.Fill(nTow2);
2628   nTowers3.Fill(nTow3);
2629   nTowers4.Fill(nTow4);
2630 
2631   // *********************
2632   // MidPoint
2633   //
2634   UsedTowerList.clear();
2635   TowerUsedInJets.clear();
2636   TowerNotUsedInJets.clear();
2637 
2638   // --- Loop over jets and make a list of all the used towers
2639   evt.getByLabel(CaloJetAlgorithm1, jets);
2640   for (CaloJetCollection::const_iterator ijet = jets->begin(); ijet != jets->end(); ijet++) {
2641     //    if (jetPt>5.0) {
2642 
2643     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
2644     int nConstituents = jetCaloRefs.size();
2645     for (int i = 0; i < nConstituents; i++) {
2646       UsedTowerList.push_back(jetCaloRefs[i]);
2647     }
2648 
2649     //    }
2650   }
2651 
2652   int NTowersUsed = UsedTowerList.size();
2653 
2654   // --- Loop over towers and make a lists of used and unused towers
2655   for (CaloTowerCollection::const_iterator tower = caloTowers->begin(); tower != caloTowers->end(); tower++) {
2656     const CaloTower& t = *tower;
2657     Double_t et = tower->et();
2658 
2659     if (et > 0) {
2660       bool used = false;
2661 
2662       for (int i = 0; i < NTowersUsed; i++) {
2663         if (tower->id() == UsedTowerList[i]->id()) {
2664           used = true;
2665           break;
2666         }
2667       }
2668 
2669       if (used) {
2670         TowerUsedInJets.push_back(t);
2671       } else {
2672         TowerNotUsedInJets.push_back(t);
2673       }
2674     }
2675   }
2676 
2677   int nUsed = TowerUsedInJets.size();
2678   int nNotUsed = TowerNotUsedInJets.size();
2679 
2680   SumEtJets = 0;
2681   SumEtNotJets = 0;
2682 
2683   for (int i = 0; i < nUsed; i++) {
2684     SumEtJets += TowerUsedInJets[i].et();
2685   }
2686   h_jetEt1.Fill(SumEtJets);
2687 
2688   for (int i = 0; i < nNotUsed; i++) {
2689     if (TowerNotUsedInJets[i].et() > 0.5)
2690       SumEtNotJets += TowerNotUsedInJets[i].et();
2691     h_missEt1.Fill(TowerNotUsedInJets[i].et());
2692     h_missEt1s.Fill(TowerNotUsedInJets[i].et());
2693   }
2694   h_totMissEt1.Fill(SumEtNotJets);
2695 
2696   // *********************
2697   // SISCone
2698   //
2699   UsedTowerList.clear();
2700   TowerUsedInJets.clear();
2701   TowerNotUsedInJets.clear();
2702 
2703   // --- Loop over jets and make a list of all the used towers
2704   evt.getByLabel(CaloJetAlgorithm2, jets);
2705   for (CaloJetCollection::const_iterator ijet = jets->begin(); ijet != jets->end(); ijet++) {
2706     //    if (jetPt>5.0) {
2707 
2708     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
2709     int nConstituents = jetCaloRefs.size();
2710     for (int i = 0; i < nConstituents; i++) {
2711       UsedTowerList.push_back(jetCaloRefs[i]);
2712     }
2713 
2714     //    }
2715   }
2716 
2717   //  Handle<CaloTowerCollection> caloTowers;
2718   //  evt.getByLabel( "towerMaker", caloTowers );
2719 
2720   NTowersUsed = UsedTowerList.size();
2721 
2722   // --- Loop over towers and make a lists of used and unused towers
2723   for (CaloTowerCollection::const_iterator tower = caloTowers->begin(); tower != caloTowers->end(); tower++) {
2724     const CaloTower& t = *tower;
2725     Double_t et = tower->et();
2726 
2727     if (et > 0) {
2728       bool used = false;
2729 
2730       for (int i = 0; i < NTowersUsed; i++) {
2731         if (tower->id() == UsedTowerList[i]->id()) {
2732           used = true;
2733           break;
2734         }
2735       }
2736 
2737       if (used) {
2738         TowerUsedInJets.push_back(t);
2739       } else {
2740         TowerNotUsedInJets.push_back(t);
2741       }
2742     }
2743   }
2744 
2745   nUsed = TowerUsedInJets.size();
2746   nNotUsed = TowerNotUsedInJets.size();
2747 
2748   SumEtJets = 0;
2749   SumEtNotJets = 0;
2750 
2751   for (int i = 0; i < nUsed; i++) {
2752     SumEtJets += TowerUsedInJets[i].et();
2753   }
2754   h_jetEt2.Fill(SumEtJets);
2755 
2756   for (int i = 0; i < nNotUsed; i++) {
2757     if (TowerNotUsedInJets[i].et() > 0.5)
2758       SumEtNotJets += TowerNotUsedInJets[i].et();
2759     h_missEt2.Fill(TowerNotUsedInJets[i].et());
2760     h_missEt2s.Fill(TowerNotUsedInJets[i].et());
2761   }
2762   h_totMissEt2.Fill(SumEtNotJets);
2763 
2764   // *********************
2765   // KtClus
2766   //
2767   UsedTowerList.clear();
2768   TowerUsedInJets.clear();
2769   TowerNotUsedInJets.clear();
2770 
2771   // --- Loop over jets and make a list of all the used towers
2772   evt.getByLabel(CaloJetAlgorithm3, jets);
2773   for (CaloJetCollection::const_iterator ijet = jets->begin(); ijet != jets->end(); ijet++) {
2774     //    if (jetPt>5.0) {
2775 
2776     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
2777     int nConstituents = jetCaloRefs.size();
2778     for (int i = 0; i < nConstituents; i++) {
2779       UsedTowerList.push_back(jetCaloRefs[i]);
2780     }
2781 
2782     //    }
2783   }
2784 
2785   //  Handle<CaloTowerCollection> caloTowers;
2786   //  evt.getByLabel( "towerMaker", caloTowers );
2787 
2788   NTowersUsed = UsedTowerList.size();
2789 
2790   // --- Loop over towers and make a lists of used and unused towers
2791   for (CaloTowerCollection::const_iterator tower = caloTowers->begin(); tower != caloTowers->end(); tower++) {
2792     const CaloTower& t = *tower;
2793     Double_t et = tower->et();
2794 
2795     if (et > 0) {
2796       bool used = false;
2797 
2798       for (int i = 0; i < NTowersUsed; i++) {
2799         if (tower->id() == UsedTowerList[i]->id()) {
2800           used = true;
2801           break;
2802         }
2803       }
2804 
2805       if (used) {
2806         TowerUsedInJets.push_back(t);
2807       } else {
2808         TowerNotUsedInJets.push_back(t);
2809       }
2810     }
2811   }
2812 
2813   nUsed = TowerUsedInJets.size();
2814   nNotUsed = TowerNotUsedInJets.size();
2815 
2816   SumEtJets = 0;
2817   SumEtNotJets = 0;
2818 
2819   for (int i = 0; i < nUsed; i++) {
2820     SumEtJets += TowerUsedInJets[i].et();
2821   }
2822   h_jetEt3.Fill(SumEtJets);
2823 
2824   for (int i = 0; i < nNotUsed; i++) {
2825     if (TowerNotUsedInJets[i].et() > 0.5)
2826       SumEtNotJets += TowerNotUsedInJets[i].et();
2827     h_missEt3.Fill(TowerNotUsedInJets[i].et());
2828     h_missEt3s.Fill(TowerNotUsedInJets[i].et());
2829   }
2830   h_totMissEt3.Fill(SumEtNotJets);
2831 }
2832 
2833 void myFastSimVal::endJob() {
2834   //Write out the histogram file.
2835   m_file->Write();
2836 }
2837 #include "FWCore/Framework/interface/MakerMacros.h"
2838 DEFINE_FWK_MODULE(myFastSimVal);