Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:57

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