CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoJets/JetAnalyzers/src/myFastSimVal.cc

Go to the documentation of this file.
00001 // myFastSimVal.cc
00002 // Description:  Comparison between Jet Algorithms
00003 // Author: Frank Chlebana
00004 // Date:  08 - August - 2007
00005 // 
00006 #include "RecoJets/JetAnalyzers/interface/myFastSimVal.h"
00007 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
00008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00009 #include "DataFormats/JetReco/interface/CaloJet.h"
00010 #include "DataFormats/JetReco/interface/GenJet.h"
00011 #include "DataFormats/Math/interface/deltaR.h"
00012 #include "DataFormats/Math/interface/deltaPhi.h"
00013 // #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
00014 #include "DataFormats/Candidate/interface/Candidate.h"
00015 // #include "FWCore/Framework/interface/Handle.h"
00016 #include "FWCore/Framework/interface/Event.h"
00017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00018 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
00019 #include <TROOT.h>
00020 #include <TSystem.h>
00021 #include <TFile.h>
00022 #include <TCanvas.h>
00023 #include <cmath>
00024 using namespace edm;
00025 using namespace reco;
00026 using namespace std;
00027 
00028 #define DEBUG 1
00029 // #define MAXJETS 50
00030 #define MAXJETS 100
00031 
00032 // Get the algorithm of the jet collections we will read from the .cfg file 
00033 // which defines the value of the strings CaloJetAlgorithm and GenJetAlgorithm.
00034 myFastSimVal::myFastSimVal( const ParameterSet & cfg ) :
00035   CaloJetAlgorithm1( cfg.getParameter<string>( "CaloJetAlgorithm1" ) ), 
00036   CaloJetAlgorithm2( cfg.getParameter<string>( "CaloJetAlgorithm2" ) ), 
00037   CaloJetAlgorithm3( cfg.getParameter<string>( "CaloJetAlgorithm3" ) ), 
00038   CaloJetAlgorithm4( cfg.getParameter<string>( "CaloJetAlgorithm4" ) ), 
00039   GenJetAlgorithm1( cfg.getParameter<string>( "GenJetAlgorithm1" ) ),
00040   GenJetAlgorithm2( cfg.getParameter<string>( "GenJetAlgorithm2" ) ),
00041   GenJetAlgorithm3( cfg.getParameter<string>( "GenJetAlgorithm3" ) ),
00042   GenJetAlgorithm4( cfg.getParameter<string>( "GenJetAlgorithm4" ) ),
00043   JetCorrectionService( cfg.getParameter<string>( "JetCorrectionService" ) )
00044 {
00045 }
00046 
00047 
00048 
00049 int nEvent = 0;
00050 
00051 void myFastSimVal::beginJob( ) {
00052 
00053   // Open the histogram file and book some associated histograms
00054   m_file=new TFile("histo.root","RECREATE"); 
00055 
00056   tMassGen      =  TH1F("tMassGen","T Mass Gen",100,0,200);
00057   tbarMassGen   =  TH1F("tbarMassGen","Tbar Mass Gen",100,0,200);
00058 
00059   tMass         =  TH1F("tMass","T Mass",100,0,200);
00060   tbarMass      =  TH1F("tbarMass","Tbar Mass",100,0,200);
00061 
00062   topMassParton =  TH1F("topMassParton","Top Mass Parton",100,0,200);
00063   topMass1      =  TH1F("topMass1","Top Mass 1",100,0,200);
00064   topMass2      =  TH1F("topMass2","Top Mass 2",100,0,200);
00065   topMass3      =  TH1F("topMass3","Top Mass 3",100,0,200);
00066 
00067   ZpMass         =  TH1F("ZpMass","Generated Zp Mass",160,0,8000);
00068   ZpMassGen      =  TH1F("ZpMassGen","Gen Zp Mass",160,0,8000);
00069   ZpMassMatched1 =  TH1F("ZpMassMatched1","Calor Zp Mass 1",160,0,8000);
00070   ZpMassMatched2 =  TH1F("ZpMassMatched2","Calor Zp Mass 2",160,0,8000);
00071   ZpMassMatched3 =  TH1F("ZpMassMatched3","Calor Zp Mass 3",160,0,8000);
00072 
00073   ZpMassGen10      =  TH1F("ZpMassGen10","Gen Zp Mass",160,0,8000);
00074   ZpMassGen13      =  TH1F("ZpMassGen13","Gen Zp Mass",160,0,8000);
00075   ZpMassGen40      =  TH1F("ZpMassGen40","Gen Zp Mass",160,0,8000);
00076 
00077   ZpMass_700_10      =  TH1F("ZpMass_700_10","Parton Zp Mass",100,0,1000);
00078   ZpMass_700_13      =  TH1F("ZpMass_700_13","Parton Zp Mass",100,0,1000);
00079   ZpMass_700_40      =  TH1F("ZpMass_700_40","Parton Zp Mass",100,0,1000);
00080 
00081   ZpMassGen_700_10      =  TH1F("ZpMassGen_700_10","Gen Zp Mass",100,0,1000);
00082   ZpMassGen_700_13      =  TH1F("ZpMassGen_700_13","Gen Zp Mass",100,0,1000);
00083   ZpMassGen_700_40      =  TH1F("ZpMassGen_700_40","Gen Zp Mass",100,0,1000);
00084 
00085   ZpMassGen_2000_10      =  TH1F("ZpMassGen_2000_10","Gen Zp Mass",100,1500,2500);
00086   ZpMassGen_2000_13      =  TH1F("ZpMassGen_2000_13","Gen Zp Mass",100,1500,2500);
00087   ZpMassGen_2000_40      =  TH1F("ZpMassGen_2000_40","Gen Zp Mass",100,1500,2500);
00088 
00089   ZpMass_2000_10      =  TH1F("ZpMass_2000_10","Parton Zp Mass",100,1500,2500);
00090   ZpMass_2000_13      =  TH1F("ZpMass_2000_13","Parton Zp Mass",100,1500,2500);
00091   ZpMass_2000_40      =  TH1F("ZpMass_2000_40","Parton Zp Mass",100,1500,2500);
00092 
00093   ZpMassGen_5000_10      =  TH1F("ZpMassGen_5000_10","Gen Zp Mass",150,4000,5500);
00094   ZpMassGen_5000_13      =  TH1F("ZpMassGen_5000_13","Gen Zp Mass",150,4000,5500);
00095   ZpMassGen_5000_40      =  TH1F("ZpMassGen_5000_40","Gen Zp Mass",150,4000,5500);
00096 
00097   ZpMass_5000_10      =  TH1F("ZpMass_5000_10","Parton Zp Mass",150,4000,5500);
00098   ZpMass_5000_13      =  TH1F("ZpMass_5000_13","Parton Zp Mass",150,4000,5500);
00099   ZpMass_5000_40      =  TH1F("ZpMass_5000_40","Parton Zp Mass",150,4000,5500);
00100 
00101   ZpMassRes101     =  TH1F("ZpMassRes101","Zp Mass Resolution 1",100,-2,2);
00102   ZpMassRes102     =  TH1F("ZpMassRes102","Zp Mass Resolution 2",100,-2,2);
00103   ZpMassRes103     =  TH1F("ZpMassRes103","Zp Mass Resolution 3",100,-2,2);
00104 
00105   ZpMassRes131     =  TH1F("ZpMassRes131","Zp Mass Resolution 1",100,-2,2);
00106   ZpMassRes132     =  TH1F("ZpMassRes132","Zp Mass Resolution 2",100,-2,2);
00107   ZpMassRes133     =  TH1F("ZpMassRes133","Zp Mass Resolution 3",100,-2,2);
00108 
00109   ZpMassRes401     =  TH1F("ZpMassRes401","Zp Mass Resolution 1",100,-2,2);
00110   ZpMassRes402     =  TH1F("ZpMassRes402","Zp Mass Resolution 2",100,-2,2);
00111   ZpMassRes403     =  TH1F("ZpMassRes403","Zp Mass Resolution 3",100,-2,2);
00112 
00113   ZpMassResL101     =  TH1F("ZpMassResL101","Zp Mass Resolution Leading Jets 1",100,0,2);
00114   ZpMassResL102     =  TH1F("ZpMassResL102","Zp Mass Resolution Leading Jets 2",100,0,2);
00115   ZpMassResL103     =  TH1F("ZpMassResL103","Zp Mass Resolution Leading Jets 3",100,0,2);
00116 
00117   ZpMassResRL101     =  TH1F("ZpMassResRL101","Zp Mass Res. Ratio Leading Jets 1",100,0,2);
00118   ZpMassResRL102     =  TH1F("ZpMassResRL102","Zp Mass Res. Ratio Leading Jets 2",100,0,2);
00119   ZpMassResRL103     =  TH1F("ZpMassResRL103","Zp Mass Res. Ratio Leading Jets 3",100,0,2);
00120 
00121   ZpMassResRLoP101     =  TH1F("ZpMassResRLoP101","Zp Mass RLoP Ratio Leading Jets 1",100,0,2);
00122   ZpMassResRLoP102     =  TH1F("ZpMassResRLoP102","Zp Mass RLoP Ratio Leading Jets 2",100,0,2);
00123   ZpMassResRLoP103     =  TH1F("ZpMassResRLoP103","Zp Mass RLoP Ratio Leading Jets 3",100,0,2);
00124 
00125   ZpMassResPRL101     =  TH1F("ZpMassResPRL101","Zp Mass Res. P Ratio Leading Jets 1",100,0,2);
00126   ZpMassResPRL102     =  TH1F("ZpMassResPRL102","Zp Mass Res. P Ratio Leading Jets 2",100,0,2);
00127   ZpMassResPRL103     =  TH1F("ZpMassResPRL103","Zp Mass Res. P Ratio Leading Jets 3",100,0,2);
00128 
00129 
00130   ZpMassResL131     =  TH1F("ZpMassResL131","Zp Mass Resolution Leading Jets 1",100,0,2);
00131   ZpMassResL132     =  TH1F("ZpMassResL132","Zp Mass Resolution Leading Jets 2",100,0,2);
00132   ZpMassResL133     =  TH1F("ZpMassResL133","Zp Mass Resolution Leading Jets 3",100,0,2);
00133 
00134   ZpMassResRL131     =  TH1F("ZpMassResRL131","Zp Mass Res. Ratio Leading Jets 1",100,0,2);
00135   ZpMassResRL132     =  TH1F("ZpMassResRL132","Zp Mass Res. Ratio Leading Jets 2",100,0,2);
00136   ZpMassResRL133     =  TH1F("ZpMassResRL133","Zp Mass Res. Ratio Leading Jets 3",100,0,2);
00137 
00138   ZpMassResRLoP131     =  TH1F("ZpMassResRLoP131","Zp Mass RLoP Ratio Leading Jets 1",100,0,2);
00139   ZpMassResRLoP132     =  TH1F("ZpMassResRLoP132","Zp Mass RLoP Ratio Leading Jets 2",100,0,2);
00140   ZpMassResRLoP133     =  TH1F("ZpMassResRLoP133","Zp Mass RLoP Ratio Leading Jets 3",100,0,2);
00141 
00142   ZpMassResPRL131     =  TH1F("ZpMassResPRL131","Zp Mass Res. P Ratio Leading Jets 1",100,0,2);
00143   ZpMassResPRL132     =  TH1F("ZpMassResPRL132","Zp Mass Res. P Ratio Leading Jets 2",100,0,2);
00144   ZpMassResPRL133     =  TH1F("ZpMassResPRL133","Zp Mass Res. P Ratio Leading Jets 3",100,0,2);
00145 
00146 
00147   ZpMassResL401     =  TH1F("ZpMassResL401","Zp Mass Resolution Leading Jets 1",100,0,2);
00148   ZpMassResL402     =  TH1F("ZpMassResL402","Zp Mass Resolution Leading Jets 2",100,0,2);
00149   ZpMassResL403     =  TH1F("ZpMassResL403","Zp Mass Resolution Leading Jets 3",100,0,2);
00150 
00151   ZpMassResRL401     =  TH1F("ZpMassResRL401","Zp Mass Res. Ratio Leading Jets 1",100,0,2);
00152   ZpMassResRL402     =  TH1F("ZpMassResRL402","Zp Mass Res. Ratio Leading Jets 2",100,0,2);
00153   ZpMassResRL403     =  TH1F("ZpMassResRL403","Zp Mass Res. Ratio Leading Jets 3",100,0,2);
00154 
00155   ZpMassResRLoP401     =  TH1F("ZpMassResRLoP401","Zp Mass RLoP Ratio Leading Jets 1",100,0,2);
00156   ZpMassResRLoP402     =  TH1F("ZpMassResRLoP402","Zp Mass RLoP Ratio Leading Jets 2",100,0,2);
00157   ZpMassResRLoP403     =  TH1F("ZpMassResRLoP403","Zp Mass RLoP Ratio Leading Jets 3",100,0,2);
00158 
00159   ZpMassResPRL401     =  TH1F("ZpMassResPRL401","Zp Mass Res. P Ratio Leading Jets 1",100,0,2);
00160   ZpMassResPRL402     =  TH1F("ZpMassResPRL402","Zp Mass Res. P Ratio Leading Jets 2",100,0,2);
00161   ZpMassResPRL403     =  TH1F("ZpMassResPRL403","Zp Mass Res. P Ratio Leading Jets 3",100,0,2);
00162 
00163   dijetMass1 =  TH1F("dijetMass1","DiJet Mass 1",100,0,4000);
00164   dijetMass12 =  TH1F("dijetMass12","DiJet Mass 1 2",100,0,6000);
00165   dijetMass13 =  TH1F("dijetMass13","DiJet Mass 1 3",100,0,12000);
00166   dijetMass2 =  TH1F("dijetMass2","DiJet Mass 2",100,0,4000);
00167   dijetMass22 =  TH1F("dijetMass22","DiJet Mass 2 2",100,0,6000);
00168   dijetMass23 =  TH1F("dijetMass23","DiJet Mass 2 3",100,0,12000);
00169   dijetMass3 =  TH1F("dijetMass3","DiJet Mass 3",100,0,4000);
00170   dijetMass32 =  TH1F("dijetMass32","DiJet Mass 3 2",100,0,6000);
00171   dijetMass33 =  TH1F("dijetMass33","DiJet Mass 3 3",100,0,12000);
00172   dijetMass4 =  TH1F("dijetMass4","DiJet Mass 4",100,0,4000);
00173   dijetMass42 =  TH1F("dijetMass42","DiJet Mass 4 2",100,0,6000);
00174   dijetMass43 =  TH1F("dijetMass43","DiJet Mass 4 3",100,0,12000);
00175 
00176   dijetMass101 =  TH1F("dijetMass101","DiJet Mass 1",100,0,6000);
00177   dijetMass131 =  TH1F("dijetMass131","DiJet Mass 1",100,0,6000);
00178   dijetMass401 =  TH1F("dijetMass401","DiJet Mass 1",100,0,6000);
00179 
00180   dijetMass102 =  TH1F("dijetMass102","DiJet Mass 2",100,0,6000);
00181   dijetMass132 =  TH1F("dijetMass132","DiJet Mass 2",100,0,6000);
00182   dijetMass402 =  TH1F("dijetMass402","DiJet Mass 2",100,0,6000);
00183 
00184   dijetMass103 =  TH1F("dijetMass103","DiJet Mass 3",100,0,10000);
00185   dijetMass133 =  TH1F("dijetMass133","DiJet Mass 3",100,0,10000);
00186   dijetMass403 =  TH1F("dijetMass403","DiJet Mass 3",100,0,10000);
00187 
00188   dijetMass_700_101 =  TH1F("dijetMass_700_101","DiJet Mass 1",100,0,1000);
00189   dijetMass_700_131 =  TH1F("dijetMass_700_131","DiJet Mass 1",100,0,1000);
00190   dijetMass_700_401 =  TH1F("dijetMass_700_401","DiJet Mass 1",100,0,1000);
00191 
00192   dijetMass_2000_101 =  TH1F("dijetMass_2000_101","DiJet Mass 1",100,1500,2500);
00193   dijetMass_2000_131 =  TH1F("dijetMass_2000_131","DiJet Mass 1",100,1500,2500);
00194   dijetMass_2000_401 =  TH1F("dijetMass_2000_401","DiJet Mass 1",100,1500,2500);
00195 
00196   dijetMass_5000_101 =  TH1F("dijetMass_5000_101","DiJet Mass 1",150,4000,5500);
00197   dijetMass_5000_131 =  TH1F("dijetMass_5000_131","DiJet Mass 1",150,4000,5500);
00198   dijetMass_5000_401 =  TH1F("dijetMass_5000_401","DiJet Mass 1",150,4000,5500);
00199 
00200 
00201   dijetMassCor1   =  TH1F("dijetMassCor1","DiJet Mass 1",160,0,8000);
00202   dijetMassCor101 =  TH1F("dijetMassCor101","DiJet Mass Cor 101",160,0,8000);
00203   dijetMassCor131 =  TH1F("dijetMassCor131","DiJet Mass Cor 131",160,0,8000);
00204   dijetMassCor401 =  TH1F("dijetMassCor401","DiJet Mass Cor 401",160,0,8000);
00205 
00206   dijetMassCor_700_1   =  TH1F("dijetMassCor_700_1","DiJet Mass 1",100,0,1000);
00207   dijetMassCor_700_101 =  TH1F("dijetMassCor_700_101","DiJet Mass Cor 101",100,0,1000);
00208   dijetMassCor_700_131 =  TH1F("dijetMassCor_700_131","DiJet Mass Cor 131",100,0,1000);
00209   dijetMassCor_700_401 =  TH1F("dijetMassCor_700_401","DiJet Mass Cor 401",100,0,1000);
00210 
00211   dijetMassCor_2000_1   =  TH1F("dijetMassCor_2000_1","DiJet Mass 1",100,1500,2500);
00212   dijetMassCor_2000_101 =  TH1F("dijetMassCor_2000_101","DiJet Mass Cor 101",100,1500,2500);
00213   dijetMassCor_2000_131 =  TH1F("dijetMassCor_2000_131","DiJet Mass Cor 131",100,1500,2500);
00214   dijetMassCor_2000_401 =  TH1F("dijetMassCor_2000_401","DiJet Mass Cor 401",100,1500,2500);
00215 
00216   dijetMassCor_5000_1   =  TH1F("dijetMassCor_5000_1","DiJet Mass 1",150,4000,5500);
00217   dijetMassCor_5000_101 =  TH1F("dijetMassCor_5000_101","DiJet Mass Cor 101",150,4000,5500);
00218   dijetMassCor_5000_131 =  TH1F("dijetMassCor_5000_131","DiJet Mass Cor 131",150,4000,5500);
00219   dijetMassCor_5000_401 =  TH1F("dijetMassCor_5000_401","DiJet Mass Cor 401",150,4000,5500);
00220 
00221   dijetMassP1 =  TH1F("dijetMassP1","DiJet Mass P 1",160,0,8000);
00222   dijetMassP2 =  TH1F("dijetMassP2","DiJet Mass P 2",160,0,8000);
00223   dijetMassP3 =  TH1F("dijetMassP3","DiJet Mass P 3",160,0,8000);
00224 
00225 
00226   dijetMassP101 =  TH1F("dijetMassP101","DiJet Mass P 1",160,0,8000);
00227   dijetMassP131 =  TH1F("dijetMassP131","DiJet Mass P 1",160,0,8000);
00228   dijetMassP401 =  TH1F("dijetMassP401","DiJet Mass P 1",160,0,8000);
00229 
00230   dijetMassP_700_101 =  TH1F("dijetMassP_700_101","DiJet Mass P 1",100,0,1000);
00231   dijetMassP_700_131 =  TH1F("dijetMassP_700_131","DiJet Mass P 1",100,0,1000);
00232   dijetMassP_700_401 =  TH1F("dijetMassP_700_401","DiJet Mass P 1",100,0,1000);
00233 
00234   dijetMassP_2000_101 =  TH1F("dijetMassP_2000_101","DiJet Mass P 1",100,1500,2500);
00235   dijetMassP_2000_131 =  TH1F("dijetMassP_2000_131","DiJet Mass P 1",100,1500,2500);
00236   dijetMassP_2000_401 =  TH1F("dijetMassP_2000_401","DiJet Mass P 1",100,1500,2500);
00237 
00238   dijetMassP_5000_101 =  TH1F("dijetMassP_5000_101","DiJet Mass P 1",150,4000,5500);
00239   dijetMassP_5000_131 =  TH1F("dijetMassP_5000_131","DiJet Mass P 1",150,4000,5500);
00240   dijetMassP_5000_401 =  TH1F("dijetMassP_5000_401","DiJet Mass P 1",150,4000,5500);
00241 
00242   totEneLeadJetEta1_1 = TH1F("totEneLeadJetEta1_1","Total Energy Lead Jet Eta1 1",100,0,1500);
00243   totEneLeadJetEta2_1 = TH1F("totEneLeadJetEta2_1","Total Energy Lead Jet Eta2 1",100,0,1500);
00244   totEneLeadJetEta3_1 = TH1F("totEneLeadJetEta3_1","Total Energy Lead Jet Eta3 1",100,0,1500);
00245   hadEneLeadJetEta1_1 = TH1F("hadEneLeadJetEta1_1","Hadronic Energy Lead Jet Eta1 1",100,0,1500);
00246   hadEneLeadJetEta2_1 = TH1F("hadEneLeadJetEta2_1","Hadronic Energy Lead Jet Eta2 1",100,0,1500);
00247   hadEneLeadJetEta3_1 = TH1F("hadEneLeadJetEta3_1","Hadronic Energy Lead Jet Eta3 1",100,0,1500);
00248   emEneLeadJetEta1_1 = TH1F("emEneLeadJetEta1_1","EM Energy Lead Jet Eta1 1",100,0,1500);
00249   emEneLeadJetEta2_1 = TH1F("emEneLeadJetEta2_1","EM Energy Lead Jet Eta2 1",100,0,1500);
00250   emEneLeadJetEta3_1 = TH1F("emEneLeadJetEta3_1","EM Energy Lead Jet Eta3 1",100,0,1500);
00251 
00252   totEneLeadJetEta1_2 = TH1F("totEneLeadJetEta1_2","Total Energy Lead Jet Eta1 2",100,0,6000);
00253   totEneLeadJetEta2_2 = TH1F("totEneLeadJetEta2_2","Total Energy Lead Jet Eta2 2",100,0,6000);
00254   totEneLeadJetEta3_2 = TH1F("totEneLeadJetEta3_2","Total Energy Lead Jet Eta3 2",100,0,6000);
00255   hadEneLeadJetEta1_2 = TH1F("hadEneLeadJetEta1_2","Hadronic Energy Lead Jet Eta1 2",100,0,6000);
00256   hadEneLeadJetEta2_2 = TH1F("hadEneLeadJetEta2_2","Hadronic Energy Lead Jet Eta2 2",100,0,6000);
00257   hadEneLeadJetEta3_2 = TH1F("hadEneLeadJetEta3_2","Hadronic Energy Lead Jet Eta3 2",100,0,6000);
00258   emEneLeadJetEta1_2 = TH1F("emEneLeadJetEta1_2","EM Energy Lead Jet Eta1 2",100,0,5000);
00259   emEneLeadJetEta2_2 = TH1F("emEneLeadJetEta2_2","EM Energy Lead Jet Eta2 2",100,0,5000);
00260   emEneLeadJetEta3_2 = TH1F("emEneLeadJetEta3_2","EM Energy Lead Jet Eta3 2",100,0,5000);
00261 
00262   hadEneLeadJet1 = TH1F("hadEneLeadJet1","Hadronic Energy Lead Jet 1",100,0,3000);
00263   hadEneLeadJet12 = TH1F("hadEneLeadJet12","Hadronic Energy Lead Jet 1 2",100,0,4000);
00264   hadEneLeadJet13 = TH1F("hadEneLeadJet13","Hadronic Energy Lead Jet 1 3",100,0,6000);
00265   hadEneLeadJet2 = TH1F("hadEneLeadJet2","Hadronic Energy Lead Jet 2",100,0,3000);
00266   hadEneLeadJet22 = TH1F("hadEneLeadJet22","Hadronic Energy Lead Jet 2 2",100,0,4000);
00267   hadEneLeadJet23 = TH1F("hadEneLeadJet23","Hadronic Energy Lead Jet 2 3",100,0,6000);
00268   hadEneLeadJet3 = TH1F("hadEneLeadJet3","Hadronic Energy Lead Jet 3",100,0,3000);
00269   hadEneLeadJet32 = TH1F("hadEneLeadJet32","Hadronic Energy Lead Jet 3 2",100,0,4000);
00270   hadEneLeadJet33 = TH1F("hadEneLeadJet33","Hadronic Energy Lead Jet 3 3",100,0,6000);
00271 
00272   emEneLeadJet1 = TH1F("emEneLeadJet1","EM Energy Lead Jet 1",100,0,1500);
00273   emEneLeadJet12 = TH1F("emEneLeadJet12","EM Energy Lead Jet 1 2",100,0,3000);
00274   emEneLeadJet13 = TH1F("emEneLeadJet13","EM Energy Lead Jet 1 3",100,0,5000);
00275   emEneLeadJet2 = TH1F("emEneLeadJet2","EM Energy Lead Jet 2",100,0,1500);
00276   emEneLeadJet22 = TH1F("emEneLeadJet22","EM Energy Lead Jet 2 2",100,0,3000);
00277   emEneLeadJet23 = TH1F("emEneLeadJet23","EM Energy Lead Jet 2 3",100,0,5000);
00278   emEneLeadJet3 = TH1F("emEneLeadJet3","EM Energy Lead Jet 3",100,0,1500);
00279   emEneLeadJet32 = TH1F("emEneLeadJet32","EM Energy Lead Jet 3 2",100,0,3000);
00280   emEneLeadJet33 = TH1F("emEneLeadJet33","EM Energy Lead Jet 3 3",100,0,5000);
00281 
00282   hadFracEta11 = TH1F("hadFracEta11","Hadronic Fraction Eta1 Jet 1",100,0,1);
00283   hadFracEta21 = TH1F("hadFracEta21","Hadronic Fraction Eta2 Jet 1",100,0,1);
00284   hadFracEta31 = TH1F("hadFracEta31","Hadronic Fraction Eta3 Jet 1",100,0,1);
00285 
00286   hadFracEta12 = TH1F("hadFracEta12","Hadronic Fraction Eta1 Jet 2",100,0,1);
00287   hadFracEta22 = TH1F("hadFracEta22","Hadronic Fraction Eta2 Jet 2",100,0,1);
00288   hadFracEta32 = TH1F("hadFracEta32","Hadronic Fraction Eta3 Jet 2",100,0,1);
00289 
00290   hadFracEta13 = TH1F("hadFracEta13","Hadronic Fraction Eta1 Jet 3",100,0,1);
00291   hadFracEta23 = TH1F("hadFracEta23","Hadronic Fraction Eta2 Jet 3",100,0,1);
00292   hadFracEta33 = TH1F("hadFracEta33","Hadronic Fraction Eta3 Jet 3",100,0,1);
00293 
00294   hadFracLeadJet1 = TH1F("hadFracLeadJet1","Hadronic Fraction Lead Jet 1",100,0,1);
00295   hadFracLeadJet2 = TH1F("hadFracLeadJet2","Hadronic Fraction Lead Jet 2",100,0,1);
00296   hadFracLeadJet3 = TH1F("hadFracLeadJet3","Hadronic Fraction Lead Jet 3",100,0,1);
00297 
00298   SumEt1 = TH1F("SumEt1","SumEt 1",100,0,1000);
00299   SumEt12 = TH1F("SumEt12","SumEt 1 2",100,0,4000);
00300   SumEt13 = TH1F("SumEt13","SumEt 1 3",100,0,15000);
00301 
00302   MET1   = TH1F("MET1",  "MET 1",100,0,200);
00303   MET12   = TH1F("MET12",  "MET 1 2",100,0,1000);
00304   MET13   = TH1F("MET13",  "MET 1 3",100,0,3000);
00305 
00306   nTowersLeadJet1  = TH1F("nTowersLeadJet1","Number of Towers Lead Jet 1",100,0,100);
00307   nTowersLeadJet2  = TH1F("nTowersLeadJet2","Number of Towers Lead Jet 2",100,0,100);
00308   nTowersLeadJet3  = TH1F("nTowersLeadJet3","Number of Towers Lead Jet 3",100,0,100);
00309 
00310   nTowersSecondJet1  = TH1F("nTowersSecondJet1","Number of Towers Second Jet 1",100,0,100);
00311   nTowersSecondJet2  = TH1F("nTowersSecondJet2","Number of Towers Second Jet 2",100,0,100);
00312   nTowersSecondJet3  = TH1F("nTowersSecondJet3","Number of Towers Second Jet 3",100,0,100);
00313 
00314   hf_PtResponse1   = TProfile("PtResponse1","Pt Response 1", 100, -5, 5, 0, 10);
00315   hf_PtResponse2   = TProfile("PtResponse2","Pt Response 2", 100, -5, 5, 0, 10);
00316   hf_PtResponse3   = TProfile("PtResponse3","Pt Response 3", 100, -5, 5, 0, 10);
00317   hf_PtResponse4   = TProfile("PtResponse4","Pt Response 4", 100, -5, 5, 0, 10);
00318 
00319   hf_TowerDelR1   = TProfile("hf_TowerDelR1","Tower Del R 1", 100, 0, 2, 0, 10);
00320   hf_TowerDelR12   = TProfile("hf_TowerDelR12","Tower Del R 1", 80, 0, 0.8, 0, 10);
00321   hf_TowerDelR2   = TProfile("hf_TowerDelR2","Tower Del R 2", 100, 0, 2, 0, 10);
00322   hf_TowerDelR22   = TProfile("hf_TowerDelR22","Tower Del R 2", 80, 0, 0.8, 0, 10);
00323   hf_TowerDelR3   = TProfile("hf_TowerDelR3","Tower Del R 3", 100, 0, 2, 0, 10);
00324   hf_TowerDelR32   = TProfile("hf_TowerDelR32","Tower Del R 3", 80, 0, 0.8, 0, 10);
00325 
00326   hf_sumTowerAllEx = TH1F("sumTowerAllEx","Tower Ex",100,-1000,1000);
00327   hf_sumTowerAllEy = TH1F("sumTowerAllEy","Tower Ey",100,-1000,1000);
00328 
00329   hf_TowerJetEt1  = TH1F("TowerJetEt1","Tower/Jet Et 1",50,0,1);
00330 
00331   nTowers1  = TH1F("nTowers1","Number of Towers pt 0.5",100,0,500);
00332   nTowers2  = TH1F("nTowers2","Number of Towers pt 1.0",100,0,500);
00333   nTowers3  = TH1F("nTowers3","Number of Towers pt 1.5",100,0,500);
00334   nTowers4  = TH1F("nTowers4","Number of Towers pt 2.0",100,0,500);
00335 
00336   nTowersLeadJetPt1  = TH1F("nTowersLeadJetPt1","Number of Towers in Lead Jet pt 0.5",100,0,200);
00337   nTowersLeadJetPt2  = TH1F("nTowersLeadJetPt2","Number of Towers in Lead Jet pt 1.0",100,0,200);
00338   nTowersLeadJetPt3  = TH1F("nTowersLeadJetPt3","Number of Towers in Lead Jet pt 1.5",100,0,200);
00339   nTowersLeadJetPt4  = TH1F("nTowersLeadJetPt4","Number of Towers in Lead Jet pt 2.0",100,0,200);
00340 
00341   TowerEtLeadJet1 = TH1F("TowerEtLeadJet1","Towers Et Lead Jet 1",100,0,2000);
00342   TowerEtLeadJet12 = TH1F("TowerEtLeadJet12","Towers Et Lead Jet 1 2",100,0,6000);
00343   TowerEtLeadJet13 = TH1F("TowerEtLeadJet13","Towers Et Lead Jet 1 3",100,0,300);
00344   TowerEtLeadJet2 = TH1F("TowerEtLeadJet2","Towers Et Lead Jet 2",100,0,2000);
00345   TowerEtLeadJet22 = TH1F("TowerEtLeadJet22","Towers Et Lead Jet 2 2",100,0,6000);
00346   TowerEtLeadJet23 = TH1F("TowerEtLeadJet23","Towers Et Lead Jet 2 3",100,0,300);
00347   TowerEtLeadJet3 = TH1F("TowerEtLeadJet3","Towers Et Lead Jet 3",100,0,2000);
00348   TowerEtLeadJet32 = TH1F("TowerEtLeadJet32","Towers Et Lead Jet 3 2",100,0,6000);
00349   TowerEtLeadJet33 = TH1F("TowerEtLeadJet33","Towers Et Lead Jet 3 3",100,0,300);
00350 
00351   hf_nJet1 = TProfile("hf_nJet1", "Num Jets 1", 100, 0, 5000, 0, 50);
00352   hf_nJet2 = TProfile("hf_nJet2", "Num Jets 2", 100, 0, 5000, 0, 50);
00353   hf_nJet3 = TProfile("hf_nJet3", "Num Jets 3", 100, 0, 5000, 0, 50);
00354   hf_nJet4 = TProfile("hf_nJet4", "Num Jets 4", 100, 0, 5000, 0, 50);
00355 
00356   hf_nJet1s = TProfile("hf_nJet1s", "Num Jets 1", 100, 0, 200, 0, 50);
00357   hf_nJet2s = TProfile("hf_nJet2s", "Num Jets 2", 100, 0, 200, 0, 50);
00358   hf_nJet3s = TProfile("hf_nJet3s", "Num Jets 3", 100, 0, 200, 0, 50);
00359   hf_nJet4s = TProfile("hf_nJet4s", "Num Jets 4", 100, 0, 200, 0, 50);
00360 
00361   hf_nJet11 = TProfile("hf_nJet11", "Num Jets 1 1", 100, 0, 3000, 0, 50);
00362   hf_nJet21 = TProfile("hf_nJet21", "Num Jets 2 1", 100, 0, 3000, 0, 50);
00363   hf_nJet31 = TProfile("hf_nJet31", "Num Jets 3 1", 100, 0, 3000, 0, 50);
00364   hf_nJet41 = TProfile("hf_nJet41", "Num Jets 4 1", 100, 0, 3000, 0, 50);
00365 
00366   dRPar1   = TH1F("dRPar1","Parton dR with matched CaloJet1",100,0,0.5);
00367   dPhiPar1 = TH1F("dPhiPar1","Parton dPhi with matched CaloJet1",200,-0.5,0.5);
00368   dEtaPar1 = TH1F("dEtaPar1","Parton dEta with matched CaloJet1",200,-0.5,0.5);
00369   dPtPar1  = TH1F("dPtPar1","Parton dPt with matched CaloJet1",200,-50,50);
00370 
00371   dRPar2   = TH1F("dRPar2","Parton dR with matched CaloJet2",100,0,0.5);
00372   dPhiPar2 = TH1F("dPhiPar2","Parton dPhi with matched CaloJet2",200,-0.5,0.5);
00373   dEtaPar2 = TH1F("dEtaPar2","Parton dEta with matched CaloJet2",200,-0.5,0.5);
00374   dPtPar2  = TH1F("dPtPar2","Parton dPt with matched CaloJet2",200,-50,50);
00375 
00376   dRPar3   = TH1F("dRPar3","Parton dR with matched CaloJet3",100,0,0.5);
00377   dPhiPar3 = TH1F("dPhiPar3","Parton dPhi with matched CaloJet3",200,-0.5,0.5);
00378   dEtaPar3 = TH1F("dEtaPar3","Parton dEta with matched CaloJet3",200,-0.5,0.5);
00379   dPtPar3  = TH1F("dPtPar3","Parton dPt with matched CaloJet3",200,-50,50);
00380 
00381   dRPar4   = TH1F("dRPar4","Parton dR with matched CaloJet4",100,0,0.5);
00382   dPhiPar4 = TH1F("dPhiPar4","Parton dPhi with matched CaloJet4",200,-0.5,0.5);
00383   dEtaPar4 = TH1F("dEtaPar4","Parton dEta with matched CaloJet4",200,-0.5,0.5);
00384   dPtPar4  = TH1F("dPtPar4","Parton dPt with matched CaloJet4",200,-50,50);
00385 
00386   dRParton    = TH1F("dRParton","dR Parton",100,0,10.0);
00387   dRPartonMin = TH1F("dRPartonMin","Min dR Parton",100,0,2.0);
00388 
00389   dR1   = TH1F("dR1","GenJets dR with matched CaloJet",100,0,0.5);
00390   dPhi1 = TH1F("dPhi1","GenJets dPhi with matched CaloJet",200,-0.5,0.5);
00391   dEta1 = TH1F("dEta1","GenJets dEta with matched CaloJet",200,-0.5,0.5);
00392   dPt1  = TH1F("dPt1","GenJets dPt with matched CaloJet",200,-100,100);
00393   dPtFrac1  = TH1F("dPtFrac1","GenJets dPt frac with matched CaloJet",100,-1,1);
00394   dPt20Frac1   = TH1F("dPt20Frac1","GenJets dPt frac with matched CaloJet",100,-1,1);
00395   dPt40Frac1   = TH1F("dPt40Frac1","GenJets dPt frac with matched CaloJet",100,-1,1);
00396   dPt80Frac1   = TH1F("dPt80Frac1","GenJets dPt frac with matched CaloJet",100,-1,1);
00397   dPt100Frac1  = TH1F("dPt100Frac1","GenJets dPt frac with matched CaloJet",100,-1,1);
00398 
00399   dR2   = TH1F("dR2","GenJets dR with matched CaloJet",100,0,0.5);
00400   dPhi2 = TH1F("dPhi2","GenJets dPhi with matched CaloJet",200,-0.5,0.5);
00401   dEta2 = TH1F("dEta2","GenJets dEta with matched CaloJet",200,-0.5,0.5);
00402   dPt2  = TH1F("dPt2","GenJets dPt with matched CaloJet",200,-100,100);
00403   dPtFrac2  = TH1F("dPtFrac2","GenJets dPt frac with matched CaloJet",100,-1,1);
00404   dPt20Frac2   = TH1F("dPt20Frac2","GenJets dPt frac with matched CaloJet",100,-1,1);
00405   dPt40Frac2   = TH1F("dPt40Frac2","GenJets dPt frac with matched CaloJet",100,-1,1);
00406   dPt80Frac2   = TH1F("dPt80Frac2","GenJets dPt frac with matched CaloJet",100,-1,1);
00407   dPt100Frac2  = TH1F("dPt100Frac2","GenJets dPt frac with matched CaloJet",100,-1,1);
00408 
00409   dR3   = TH1F("dR3","GenJets dR with matched CaloJet",100,0,0.5);
00410   dPhi3 = TH1F("dPhi3","GenJets dPhi with matched CaloJet",200,-0.5,0.5);
00411   dEta3 = TH1F("dEta3","GenJets dEta with matched CaloJet",200,-0.5,0.5);
00412   dPt3  = TH1F("dPt3","GenJets dPt with matched CaloJet",200,-100,100);
00413   dPtFrac3  = TH1F("dPtFrac3","GenJets dPt frac with matched CaloJet",100,-1,1);
00414   dPt20Frac3   = TH1F("dPt20Frac3","GenJets dPt frac with matched CaloJet",100,-1,1);
00415   dPt40Frac3   = TH1F("dPt40Frac3","GenJets dPt frac with matched CaloJet",100,-1,1);
00416   dPt80Frac3   = TH1F("dPt80Frac3","GenJets dPt frac with matched CaloJet",100,-1,1);
00417   dPt100Frac3  = TH1F("dPt100Frac3","GenJets dPt frac with matched CaloJet",100,-1,1);
00418 
00419   dR4   = TH1F("dR4","GenJets dR with matched CaloJet",100,0,0.5);
00420   dPhi4 = TH1F("dPhi4","GenJets dPhi with matched CaloJet",200,-0.5,0.5);
00421   dEta4 = TH1F("dEta4","GenJets dEta with matched CaloJet",200,-0.5,0.5);
00422   dPt4  = TH1F("dPt4","GenJets dPt with matched CaloJet",200,-100,100);
00423   dPtFrac4  = TH1F("dPtFrac4","GenJets dPt frac with matched CaloJet",100,-1,1);
00424   dPt20Frac4   = TH1F("dPt20Frac4","GenJets dPt frac with matched CaloJet",100,-1,1);
00425   dPt40Frac4   = TH1F("dPt40Frac4","GenJets dPt frac with matched CaloJet",100,-1,1);
00426   dPt80Frac4   = TH1F("dPt80Frac4","GenJets dPt frac with matched CaloJet",100,-1,1);
00427   dPt100Frac4  = TH1F("dPt100Frac4","GenJets dPt frac with matched CaloJet",100,-1,1);
00428 
00429   dR12   = TH1F("dR12","dR MidPoint - SISCone",100,0,0.5);
00430   dPhi12 = TH1F("dPhi12","dPhi MidPoint - SISCone",200,-0.5,0.5);
00431   dEta12 = TH1F("dEta12","dEta MidPoint - SISCone",200,-0.5,0.5);
00432   dPt12  = TH1F("dPt12","dPt MidPoint - SISCone",200,-100,100);
00433 
00434 
00435 
00436   h_nCalJets1 =  TH1F( "nCalJets1",  "Number of CalJets1", 20, 0, 20 );
00437   h_nCalJets2 =  TH1F( "nCalJets2",  "Number of CalJets2", 20, 0, 20 );
00438   h_nCalJets3 =  TH1F( "nCalJets3",  "Number of CalJets3", 20, 0, 20 );
00439   h_nCalJets4 =  TH1F( "nCalJets4",  "Number of CalJets4", 20, 0, 20 );
00440 
00441   h_lowPtCal11 =  TH1F( "lowPtCal11",  "Low p_{T} of CalJet1 Eta1", 100, 0, 100 );
00442   h_lowPtCal21 =  TH1F( "lowPtCal21",  "Low p_{T} of CalJet2 Eta1", 100, 0, 100 );
00443   h_lowPtCal31 =  TH1F( "lowPtCal31",  "Low p_{T} of CalJet3 Eta1", 100, 0, 100 );
00444   h_lowPtCal41 =  TH1F( "lowPtCal41",  "Low p_{T} of CalJet4 Eta1", 100, 0, 100 );
00445 
00446   h_lowPtCal12 =  TH1F( "lowPtCal12",  "Low p_{T} of CalJet1 Eta2", 100, 0, 100 );
00447   h_lowPtCal22 =  TH1F( "lowPtCal22",  "Low p_{T} of CalJet2 Eta2", 100, 0, 100 );
00448   h_lowPtCal32 =  TH1F( "lowPtCal32",  "Low p_{T} of CalJet3 Eta2", 100, 0, 100 );
00449   h_lowPtCal42 =  TH1F( "lowPtCal42",  "Low p_{T} of CalJet4 Eta2", 100, 0, 100 );
00450 
00451   h_lowPtCal13 =  TH1F( "lowPtCal13",  "Low p_{T} of CalJet1 Eta3", 100, 0, 100 );
00452   h_lowPtCal23 =  TH1F( "lowPtCal23",  "Low p_{T} of CalJet2 Eta3", 100, 0, 100 );
00453   h_lowPtCal33 =  TH1F( "lowPtCal33",  "Low p_{T} of CalJet3 Eta3", 100, 0, 100 );
00454   h_lowPtCal43 =  TH1F( "lowPtCal43",  "Low p_{T} of CalJet4 Eta3", 100, 0, 100 );
00455 
00456   h_lowPtCal1c11 =  TH1F( "lowPtCal1c11",  "Low p_{T} of CalJet1 c1 Eta1", 100, 0, 100 );
00457   h_lowPtCal2c11 =  TH1F( "lowPtCal2c11",  "Low p_{T} of CalJet2 c1 Eta1", 100, 0, 100 );
00458   h_lowPtCal3c11 =  TH1F( "lowPtCal3c11",  "Low p_{T} of CalJet3 c1 Eta1", 100, 0, 100 );
00459   h_lowPtCal4c11 =  TH1F( "lowPtCal4c11",  "Low p_{T} of CalJet4 c1 Eta1", 100, 0, 100 );
00460 
00461   h_lowPtCal1c12 =  TH1F( "lowPtCal1c12",  "Low p_{T} of CalJet1 c1 Eta2", 100, 0, 100 );
00462   h_lowPtCal2c12 =  TH1F( "lowPtCal2c12",  "Low p_{T} of CalJet2 c1 Eta2", 100, 0, 100 );
00463   h_lowPtCal3c12 =  TH1F( "lowPtCal3c12",  "Low p_{T} of CalJet3 c1 Eta2", 100, 0, 100 );
00464   h_lowPtCal4c12 =  TH1F( "lowPtCal4c12",  "Low p_{T} of CalJet4 c1 Eta2", 100, 0, 100 );
00465 
00466   h_lowPtCal1c13 =  TH1F( "lowPtCal1c13",  "Low p_{T} of CalJet1 c1 Eta3", 100, 0, 100 );
00467   h_lowPtCal2c13 =  TH1F( "lowPtCal2c13",  "Low p_{T} of CalJet2 c1 Eta3", 100, 0, 100 );
00468   h_lowPtCal3c13 =  TH1F( "lowPtCal3c13",  "Low p_{T} of CalJet3 c1 Eta3", 100, 0, 100 );
00469   h_lowPtCal4c13 =  TH1F( "lowPtCal4c13",  "Low p_{T} of CalJet4 c1 Eta3", 100, 0, 100 );
00470 
00471 
00472   matchedAllPt11 =  TH1F( "matchedAllPt11",  "p_{T} of CalJet1 Eta1", 50, 0, 250 );
00473   matchedAllPt12 =  TH1F( "matchedAllPt12",  "p_{T} of CalJet1 Eta2", 50, 0, 250 );
00474   matchedAllPt13 =  TH1F( "matchedAllPt13",  "p_{T} of CalJet1 Eta3", 50, 0, 250 );
00475   matchedPt11    =  TH1F( "matchedPt11",  "p_{T} of CalJet1 Eta1", 50, 0, 250 );
00476   matchedPt12    =  TH1F( "matchedPt12",  "p_{T} of CalJet1 Eta2", 50, 0, 250 );
00477   matchedPt13    =  TH1F( "matchedPt13",  "p_{T} of CalJet1 Eta3", 50, 0, 250 );
00478 
00479   matchedAllPt21 =  TH1F( "matchedAllPt21",  "p_{T} of CalJet2 Eta1", 50, 0, 250 );
00480   matchedAllPt22 =  TH1F( "matchedAllPt22",  "p_{T} of CalJet2 Eta2", 50, 0, 250 );
00481   matchedAllPt23 =  TH1F( "matchedAllPt23",  "p_{T} of CalJet2 Eta3", 50, 0, 250 );
00482   matchedPt21    =  TH1F( "matchedPt21",  "p_{T} of CalJet2 Eta1", 50, 0, 250 );
00483   matchedPt22    =  TH1F( "matchedPt22",  "p_{T} of CalJet2 Eta2", 50, 0, 250 );
00484   matchedPt23    =  TH1F( "matchedPt23",  "p_{T} of CalJet2 Eta3", 50, 0, 250 );
00485 
00486   matchedAllPt31 =  TH1F( "matchedAllPt31",  "p_{T} of CalJet3 Eta1", 50, 0, 250 );
00487   matchedAllPt32 =  TH1F( "matchedAllPt32",  "p_{T} of CalJet3 Eta2", 50, 0, 250 );
00488   matchedAllPt33 =  TH1F( "matchedAllPt33",  "p_{T} of CalJet3 Eta3", 50, 0, 250 );
00489   matchedPt31    =  TH1F( "matchedPt31",  "p_{T} of CalJet3 Eta1", 50, 0, 250 );
00490   matchedPt32    =  TH1F( "matchedPt32",  "p_{T} of CalJet3 Eta2", 50, 0, 250 );
00491   matchedPt33    =  TH1F( "matchedPt33",  "p_{T} of CalJet3 Eta3", 50, 0, 250 );
00492 
00493   matchedAllPt41 =  TH1F( "matchedAllPt41",  "p_{T} of CalJet4 Eta1", 50, 0, 250 );
00494   matchedAllPt42 =  TH1F( "matchedAllPt42",  "p_{T} of CalJet4 Eta2", 50, 0, 250 );
00495   matchedAllPt43 =  TH1F( "matchedAllPt43",  "p_{T} of CalJet4 Eta3", 50, 0, 250 );
00496   matchedPt41    =  TH1F( "matchedPt41",  "p_{T} of CalJet4 Eta1", 50, 0, 250 );
00497   matchedPt42    =  TH1F( "matchedPt42",  "p_{T} of CalJet4 Eta2", 50, 0, 250 );
00498   matchedPt43    =  TH1F( "matchedPt43",  "p_{T} of CalJet4 Eta3", 50, 0, 250 );
00499 
00500 
00501 
00502   h_ptCal1 =  TH1F( "ptCal1",  "p_{T} of CalJet1", 50, 0, 1000 );
00503   h_ptCal12 =  TH1F( "ptCal12",  "p_{T} of CalJet1 2", 50, 0, 6000 );
00504   h_ptCal13 =  TH1F( "ptCal13",  "p_{T} of CalJet1 2", 50, 0, 300 );
00505 
00506   h_ptCal2 =  TH1F( "ptCal2",  "p_{T} of CalJet2", 50, 0, 1000 );
00507   h_ptCal22 =  TH1F( "ptCal22",  "p_{T} of CalJet2 2", 50, 0, 6000 );
00508   h_ptCal23 =  TH1F( "ptCal23",  "p_{T} of CalJet2 2", 50, 0, 300 );
00509 
00510   h_ptCal3 =  TH1F( "ptCal3",  "p_{T} of CalJet3", 50, 0, 1000 );
00511   h_ptCal32 =  TH1F( "ptCal32",  "p_{T} of CalJet3 2", 50, 0, 6000 );
00512   h_ptCal33 =  TH1F( "ptCal33",  "p_{T} of CalJet3 2", 50, 0, 300 );
00513 
00514   h_ptCal4 =  TH1F( "ptCal4",  "p_{T} of CalJet4", 50, 0, 1000 );
00515   h_ptCal42 =  TH1F( "ptCal42",  "p_{T} of CalJet4 2", 50, 0, 6000 );
00516   h_ptCal43 =  TH1F( "ptCal43",  "p_{T} of CalJet4 2", 50, 0, 300 );
00517 
00518   h_etaCal1 = TH1F( "etaCal1", "#eta of  CalJet1", 100, -4, 4 );
00519   h_etaCal2 = TH1F( "etaCal2", "#eta of  CalJet2", 100, -4, 4 );
00520   h_etaCal3 = TH1F( "etaCal3", "#eta of  CalJet3", 100, -4, 4 );
00521   h_etaCal4 = TH1F( "etaCal4", "#eta of  CalJet4", 100, -4, 4 );
00522 
00523   h_phiCal1 = TH1F( "phiCal1", "#phi of  CalJet1", 50, -M_PI, M_PI );
00524   h_phiCal2 = TH1F( "phiCal2", "#phi of  CalJet2", 50, -M_PI, M_PI );
00525   h_phiCal3 = TH1F( "phiCal3", "#phi of  CalJet3", 50, -M_PI, M_PI );
00526   h_phiCal4 = TH1F( "phiCal4", "#phi of  CalJet4", 50, -M_PI, M_PI );
00527 
00528   h_ptCalL1 =  TH1F( "ptCalL1",  "p_{T} of CalJetL1", 50, 0, 300 );
00529   h_ptCalL12 =  TH1F( "ptCalL12",  "p_{T} of CalJetL1 2", 50, 0, 1200 );
00530   h_ptCalL13 =  TH1F( "ptCalL13",  "p_{T} of CalJetL1 3", 50, 0, 6000 );
00531   h_ptCalL2 =  TH1F( "ptCalL2",  "p_{T} of CalJetL2", 50, 0, 300 );
00532   h_ptCalL22 =  TH1F( "ptCalL22",  "p_{T} of CalJetL2 2", 50, 0, 1200 );
00533   h_ptCalL23 =  TH1F( "ptCalL23",  "p_{T} of CalJetL2 3", 50, 0, 6000 );
00534   h_ptCalL3 =  TH1F( "ptCalL3",  "p_{T} of CalJetL3", 50, 0, 300 );
00535   h_ptCalL32 =  TH1F( "ptCalL32",  "p_{T} of CalJetL3 2", 50, 0, 1200 );
00536   h_ptCalL33 =  TH1F( "ptCalL33",  "p_{T} of CalJetL3 3", 50, 0, 6000 );
00537   h_ptCalL4 =  TH1F( "ptCalL4",  "p_{T} of CalJetL4", 50, 0, 300 );
00538   h_ptCalL42 =  TH1F( "ptCalL42",  "p_{T} of CalJetL4 2", 50, 0, 1200 );
00539   h_ptCalL43 =  TH1F( "ptCalL43",  "p_{T} of CalJetL4 3", 50, 0, 6000 );
00540 
00541 
00542   h_etaCalL1 = TH1F( "etaCalL1", "#eta of  CalJetL1", 100, -4, 4 );
00543   h_etaCalL2 = TH1F( "etaCalL2", "#eta of  CalJetL2", 100, -4, 4 );
00544   h_etaCalL3 = TH1F( "etaCalL3", "#eta of  CalJetL3", 100, -4, 4 );
00545   h_etaCalL4 = TH1F( "etaCalL4", "#eta of  CalJetL4", 100, -4, 4 );
00546   h_phiCalL1 = TH1F( "phiCalL1", "#phi of  CalJetL1", 50, -M_PI, M_PI );
00547   h_phiCalL2 = TH1F( "phiCalL2", "#phi of  CalJetL2", 50, -M_PI, M_PI );
00548   h_phiCalL3 = TH1F( "phiCalL3", "#phi of  CalJetL3", 50, -M_PI, M_PI );
00549   h_phiCalL4 = TH1F( "phiCalL4", "#phi of  CalJetL4", 50, -M_PI, M_PI );
00550 
00551   h_nGenJets1 =  TH1F( "nGenJets1",  "Number of GenJets1", 20, 0, 20 );
00552   h_nGenJets2 =  TH1F( "nGenJets2",  "Number of GenJets2", 20, 0, 20 );
00553   h_nGenJets3 =  TH1F( "nGenJets3",  "Number of GenJets3", 20, 0, 20 );
00554   h_nGenJets4 =  TH1F( "nGenJets4",  "Number of GenJets4", 20, 0, 20 );
00555 
00556   h_ptGen1 =  TH1F( "ptGen1",  "p_{T} of GenJet1", 50, 0, 1000 );
00557   h_ptGen12 =  TH1F( "ptGen12",  "p_{T} of GenJet1 2", 50, 0, 6000 );
00558   h_ptGen13 =  TH1F( "ptGen13",  "p_{T} of GenJet1 3", 50, 0, 300 );
00559 
00560   h_ptGen2 =  TH1F( "ptGen2",  "p_{T} of GenJet2", 50, 0, 1000 );
00561   h_ptGen22 =  TH1F( "ptGen22",  "p_{T} of GenJet2 2", 50, 0, 6000 );
00562   h_ptGen23 =  TH1F( "ptGen23",  "p_{T} of GenJet2 3", 50, 0, 300 );
00563 
00564   h_ptGen3 =  TH1F( "ptGen3",  "p_{T} of GenJet3", 50, 0, 1000 );
00565   h_ptGen32 =  TH1F( "ptGen32",  "p_{T} of GenJet3 2", 50, 0, 6000 );
00566   h_ptGen33 =  TH1F( "ptGen33",  "p_{T} of GenJet3 3", 50, 0, 300 );
00567 
00568   h_ptGen4 =  TH1F( "ptGen4",  "p_{T} of GenJet4", 50, 0, 1000 );
00569   h_ptGen42 =  TH1F( "ptGen42",  "p_{T} of GenJet4 2", 50, 0, 6000 );
00570   h_ptGen43 =  TH1F( "ptGen43",  "p_{T} of GenJet4 3", 50, 0, 300 );
00571 
00572 
00573   h_etaGen1 = TH1F( "etaGen1", "#eta of GenJet1", 100, -4, 4 );
00574   h_etaGen2 = TH1F( "etaGen2", "#eta of GenJet2", 100, -4, 4 );
00575   h_etaGen3 = TH1F( "etaGen3", "#eta of GenJet3", 100, -4, 4 );
00576   h_phiGen1 = TH1F( "phiGen1", "#phi of GenJet1", 50, -M_PI, M_PI );
00577   h_phiGen2 = TH1F( "phiGen2", "#phi of GenJet2", 50, -M_PI, M_PI );
00578   h_phiGen3 = TH1F( "phiGen3", "#phi of GenJet3", 50, -M_PI, M_PI );
00579 
00580   h_ptGenL1 =  TH1F( "ptGenL1",  "p_{T} of GenJetL1", 50, 0, 300 );
00581   h_ptGenL12 =  TH1F( "ptGenL12",  "p_{T} of GenJetL1 2", 50, 0, 1200 );
00582   h_ptGenL13 =  TH1F( "ptGenL13",  "p_{T} of GenJetL1 3", 50, 0, 6000 );
00583   h_ptGenL2 =  TH1F( "ptGenL2",  "p_{T} of GenJetL2", 50, 0, 300 );
00584   h_ptGenL22 =  TH1F( "ptGenL22",  "p_{T} of GenJetL2 2", 50, 0, 1200 );
00585   h_ptGenL23 =  TH1F( "ptGenL23",  "p_{T} of GenJetL2 3", 50, 0, 6000 );
00586   h_ptGenL3 =  TH1F( "ptGenL3",  "p_{T} of GenJetL3", 50, 0, 300 );
00587   h_ptGenL32 =  TH1F( "ptGenL32",  "p_{T} of GenJetL32", 50, 0, 1200 );
00588   h_ptGenL33 =  TH1F( "ptGenL33",  "p_{T} of GenJetL33", 50, 0, 6000 );
00589 
00590 
00591   h_etaGenL1 = TH1F( "etaGenL1", "#eta of GenJetL1", 100, -4, 4 );
00592   h_etaGenL2 = TH1F( "etaGenL2", "#eta of GenJetL2", 100, -4, 4 );
00593   h_etaGenL3 = TH1F( "etaGenL3", "#eta of GenJetL3", 100, -4, 4 );
00594   h_phiGenL1 = TH1F( "phiGenL1", "#phi of GenJetL1", 50, -M_PI, M_PI );
00595   h_phiGenL2 = TH1F( "phiGenL2", "#phi of GenJetL2", 50, -M_PI, M_PI );
00596   h_phiGenL3 = TH1F( "phiGenL3", "#phi of GenJetL3", 50, -M_PI, M_PI );
00597 
00598   h_jetEt1 = TH1F( "jetEt1", "Total Jet Et", 100, 0, 3000 );
00599   h_jetEt2 = TH1F( "jetEt2", "Total Jet Et", 100, 0, 3000 );
00600   h_jetEt3 = TH1F( "jetEt3", "Total Jet Et", 100, 0, 3000 );
00601 
00602   h_jet1Pt1 = TH1F( "jet1Pt1", "Jet Pt", 100, 0, 3000 );
00603   h_jet2Pt1 = TH1F( "jet2Pt1", "Jet Pt", 100, 0, 3000 );
00604   h_jet3Pt1 = TH1F( "jet3Pt1", "Jet Pt", 100, 0, 3000 );
00605   h_jet4Pt1 = TH1F( "jet4Pt1", "Jet Pt", 100, 0, 3000 );
00606   h_jet5Pt1 = TH1F( "jet5Pt1", "Jet Pt", 100, 0, 3000 );
00607   h_jet6Pt1 = TH1F( "jet6Pt1", "Jet Pt", 100, 0, 3000 );
00608   h_jet7Pt1 = TH1F( "jet7Pt1", "Jet Pt", 100, 0, 3000 );
00609 
00610   h_jet1Pt2 = TH1F( "jet1Pt2", "Jet Pt", 100, 0, 3000 );
00611   h_jet2Pt2 = TH1F( "jet2Pt2", "Jet Pt", 100, 0, 3000 );
00612   h_jet3Pt2 = TH1F( "jet3Pt2", "Jet Pt", 100, 0, 3000 );
00613   h_jet4Pt2 = TH1F( "jet4Pt2", "Jet Pt", 100, 0, 3000 );
00614   h_jet5Pt2 = TH1F( "jet5Pt2", "Jet Pt", 100, 0, 3000 );
00615   h_jet6Pt2 = TH1F( "jet6Pt2", "Jet Pt", 100, 0, 3000 );
00616   h_jet7Pt2 = TH1F( "jet7Pt2", "Jet Pt", 100, 0, 3000 );
00617 
00618   h_jet1Pt3 = TH1F( "jet1Pt3", "Jet Pt", 100, 0, 3000 );
00619   h_jet2Pt3 = TH1F( "jet2Pt3", "Jet Pt", 100, 0, 3000 );
00620   h_jet3Pt3 = TH1F( "jet3Pt3", "Jet Pt", 100, 0, 3000 );
00621   h_jet4Pt3 = TH1F( "jet4Pt3", "Jet Pt", 100, 0, 3000 );
00622   h_jet5Pt3 = TH1F( "jet5Pt3", "Jet Pt", 100, 0, 3000 );
00623   h_jet6Pt3 = TH1F( "jet6Pt3", "Jet Pt", 100, 0, 3000 );
00624   h_jet7Pt3 = TH1F( "jet7Pt3", "Jet Pt", 100, 0, 3000 );
00625 
00626   h_jet1Pt4 = TH1F( "jet1Pt4", "Jet Pt", 100, 0, 3000 );
00627   h_jet2Pt4 = TH1F( "jet2Pt4", "Jet Pt", 100, 0, 3000 );
00628   h_jet3Pt4 = TH1F( "jet3Pt4", "Jet Pt", 100, 0, 3000 );
00629   h_jet4Pt4 = TH1F( "jet4Pt4", "Jet Pt", 100, 0, 3000 );
00630   h_jet5Pt4 = TH1F( "jet5Pt4", "Jet Pt", 100, 0, 3000 );
00631   h_jet6Pt4 = TH1F( "jet6Pt4", "Jet Pt", 100, 0, 3000 );
00632   h_jet7Pt4 = TH1F( "jet7Pt4", "Jet Pt", 100, 0, 3000 );
00633 
00634 
00635   h_totMissEt1 = TH1F( "totMissEt1", "Total Unclustered Et", 100, 0, 500 );
00636   h_totMissEt2 = TH1F( "totMissEt2", "Total Unclustered Et", 100, 0, 500 );
00637   h_totMissEt3 = TH1F( "totMissEt3", "Total Unclustered Et", 100, 0, 500 );
00638 
00639   h_missEt1 = TH1F( "missEt1", "Unclustered Et", 100, 0, 50 );
00640   h_missEt2 = TH1F( "missEt2", "Unclustered Et", 100, 0, 50 );
00641   h_missEt3 = TH1F( "missEt3", "Unclustered Et", 100, 0, 50 );
00642 
00643   h_missEt1s = TH1F( "missEt1s", "Unclustered Et", 100, 0, 2 );
00644   h_missEt2s = TH1F( "missEt2s", "Unclustered Et", 100, 0, 2 );
00645   h_missEt3s = TH1F( "missEt3s", "Unclustered Et", 100, 0, 2 );
00646 
00647   ParMatch1 = TH1F( "ParMatch1", "Number of Matched Jets 1", 10, 0, 10 );
00648   ParMatch2 = TH1F( "ParMatch2", "Number of Matched Jets 2", 10, 0, 10 );
00649   ParMatch3 = TH1F( "ParMatch3", "Number of Matched Jets 3", 10, 0, 10 );
00650 
00651 }
00652 
00653 
00654 void myFastSimVal::analyze( const Event& evt, const EventSetup& es ) {
00655 
00656   int EtaOk10, EtaOk13, EtaOk40;
00657 
00658   double ZpM, ZpMG, ZpMM;
00659   double LeadMass1, LeadMass2, LeadMass3, LeadMass4;
00660   double LeadMassP1, LeadMassP2, LeadMassP3;
00661 
00662   float pt1, pt2, pt3;
00663 
00664   float minJetPt = 30.;
00665   float minJetPt10 = 10.;
00666   int jetInd, allJetInd;
00667   int usedInd = -1;  
00668   //  double matchedDelR = 0.1;
00669     double matchedDelR = 0.3;
00670 
00671   ZpMG = 0;
00672   LeadMass1 = -1;
00673   LeadMass2 = -1;
00674   LeadMass3 = -1;
00675 
00676   math::XYZTLorentzVector p4tmp[2], p4cortmp[2];
00677   nEvent++;
00678 
00679   // ********************************
00680   // **** Get the CaloJet1 collection
00681   // ********************************
00682 
00683   Handle<CaloJetCollection> caloJets1;
00684   evt.getByLabel( CaloJetAlgorithm1, caloJets1 );
00685 
00686   // Count Jets above Pt cut
00687   for (int istep = 0; istep < 100; ++istep) {
00688     int     njet = 0;
00689     float ptStep = (istep * (5000./100.));
00690     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00691       if ( cal->pt() > ptStep ) njet++;      
00692     }
00693 
00694     hf_nJet1.Fill( ptStep, njet );
00695   }
00696 
00697   // Count Jets above Pt cut
00698   for (int istep = 0; istep < 100; ++istep) {
00699     int     njet = 0;
00700     float ptStep = (istep * (200./100.));
00701     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00702       if ( cal->pt() > ptStep ) njet++;      
00703     }
00704 
00705     hf_nJet1s.Fill( ptStep, njet );
00706   }
00707 
00708   // Count Jets above Pt cut
00709   for (int istep = 0; istep < 100; ++istep) {
00710     int     njet = 0;
00711     float ptStep = (istep * (3000./100.));
00712     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00713       if ( cal->pt() > ptStep ) njet++;      
00714     }
00715 
00716     hf_nJet11.Fill( ptStep, njet );
00717   }
00718 
00719 
00720   //Loop over the two leading CaloJets and fill some histograms
00721   jetInd    = 0;
00722   allJetInd = 0;
00723   EtaOk10 = 0;
00724   EtaOk13 = 0;
00725   EtaOk40 = 0;
00726 
00727   //  const JetCorrector* corrector = 
00728   //    JetCorrector::getJetCorrector (JetCorrectionService, es);
00729 
00730   double highestPt;
00731   double nextPt;
00732 
00733   highestPt = 0.0;
00734   nextPt    = 0.0;
00735   
00736 
00737   for( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {
00738     
00739     //    double scale = corrector->correction (*cal);
00740     double scale = 1.0;
00741     double corPt = scale*cal->pt();
00742     //    double corPt = cal->pt();
00743 
00744     
00745     if (corPt>highestPt) {
00746       nextPt      = highestPt;
00747       p4cortmp[1] = p4cortmp[0]; 
00748       highestPt   = corPt;
00749       p4cortmp[0] = scale*cal->p4();
00750     } else if (corPt>nextPt) {
00751       nextPt      = corPt;
00752       p4cortmp[1] = scale*cal->p4();
00753     }
00754 
00755     /***
00756     std::cout << ">>> Corr Jet: corPt = " 
00757               << corPt << ", scale = " << scale 
00758               << " pt = " << cal->pt()
00759               << " highestPt = " << highestPt 
00760               << " nextPt = "    << nextPt 
00761               << std::endl;  
00762     ****/
00763 
00764     allJetInd++;
00765     if (allJetInd == 1) {
00766       h_jet1Pt1.Fill( cal->pt() );
00767       pt1 = cal->pt();
00768       p4tmp[0] = cal->p4();
00769       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
00770       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
00771       if ( fabs(cal->eta()) < 4.0) EtaOk40++;
00772     }
00773     if (allJetInd == 2) {
00774       h_jet2Pt1.Fill( cal->pt() );
00775       p4tmp[1] = cal->p4();
00776       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
00777       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
00778       if ( fabs(cal->eta()) < 4.0) EtaOk40++;
00779     }
00780     if ( (allJetInd == 1) || (allJetInd == 2) ) {
00781 
00782       h_ptCalL1.Fill( cal->pt() );   
00783       h_ptCalL12.Fill( cal->pt() );   
00784       h_ptCalL13.Fill( cal->pt() );   
00785 
00786       h_etaCalL1.Fill( cal->eta() );
00787       h_phiCalL1.Fill( cal->phi() );
00788     }
00789 
00790     if (allJetInd == 3) h_jet3Pt1.Fill( cal->pt() );
00791     if (allJetInd == 4) h_jet4Pt1.Fill( cal->pt() );
00792     if (allJetInd == 5) h_jet5Pt1.Fill( cal->pt() );
00793     if (allJetInd == 6) h_jet6Pt1.Fill( cal->pt() );
00794     if (allJetInd == 7) h_jet7Pt1.Fill( cal->pt() );
00795 
00796     if ( fabs(cal->eta()) < 1.3) {
00797       h_lowPtCal11.Fill( cal->pt() );   
00798       if ( cal->pt() > 10.) {
00799         h_lowPtCal1c11.Fill( cal->pt() );   
00800       }
00801     }
00802     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
00803       h_lowPtCal12.Fill( cal->pt() );   
00804       if ( cal->pt() > 10.) {
00805         h_lowPtCal1c12.Fill( cal->pt() );   
00806       }
00807     }
00808     if ( fabs(cal->eta()) > 3.0) {
00809       h_lowPtCal13.Fill( cal->pt() );   
00810       if ( cal->pt() > 10.) {
00811         h_lowPtCal1c13.Fill( cal->pt() );   
00812       }
00813     }
00814 
00815 
00816 
00817     if ( cal->pt() > minJetPt) {
00818       //    std::cout << "CALO JET1 #" << jetInd << std::endl << cal->print() << std::endl;
00819       h_ptCal1.Fill( cal->pt() );   
00820       h_ptCal12.Fill( cal->pt() );   
00821       h_ptCal13.Fill( cal->pt() );   
00822 
00823       h_etaCal1.Fill( cal->eta() );
00824       h_phiCal1.Fill( cal->phi() );
00825       jetInd++;
00826     }
00827   }
00828 
00829   //  h_nCalJets1.Fill( caloJets1->size() );   
00830   h_nCalJets1.Fill( jetInd ); 
00831   if (jetInd > 1) {
00832     LeadMass1 = (p4tmp[0]+p4tmp[1]).mass();
00833     dijetMass1.Fill( LeadMass1 );    
00834     dijetMass12.Fill( LeadMass1 );    
00835     dijetMass13.Fill( LeadMass1 );    
00836     if (EtaOk10 == 2) {
00837       dijetMass101.Fill( LeadMass1 );        
00838       dijetMass102.Fill( LeadMass1 );  
00839       dijetMass103.Fill( LeadMass1 );  
00840       dijetMass_700_101.Fill( LeadMass1 );  
00841       dijetMass_2000_101.Fill( LeadMass1 );  
00842       dijetMass_5000_101.Fill( LeadMass1 );  
00843     }
00844     if (EtaOk13 == 2) {
00845       dijetMass131.Fill( LeadMass1 );  
00846       dijetMass132.Fill( LeadMass1 );  
00847       dijetMass133.Fill( LeadMass1 );  
00848       dijetMass_700_131.Fill( LeadMass1 );  
00849       dijetMass_2000_131.Fill( LeadMass1 );  
00850       dijetMass_5000_131.Fill( LeadMass1 );  
00851     }
00852     if (EtaOk40 == 2) {
00853       dijetMass401.Fill( LeadMass1 );        
00854       dijetMass402.Fill( LeadMass1 );  
00855       dijetMass403.Fill( LeadMass1 );  
00856       dijetMass_700_401.Fill( LeadMass1 );  
00857       dijetMass_2000_401.Fill( LeadMass1 );  
00858       dijetMass_5000_401.Fill( LeadMass1 );  
00859     }
00860 
00861     LeadMass1 = (p4cortmp[0]+p4cortmp[1]).mass();
00862 
00863     /****
00864     if (LeadMass1 < 30.) {
00865       std::cout << " XXX Low Mass " 
00866                 << (p4tmp[0]+p4tmp[1]).mass() 
00867                 << " / " 
00868                 << (p4cortmp[0]+p4cortmp[1]).mass() 
00869                 << std::endl;
00870 
00871       std::cout << " p4 1 = " << p4tmp[0]
00872                 << " p4 2 = " << p4tmp[1]
00873                 << " p4 cor 1 = " << p4cortmp[0]
00874                 << " p4 cor 2 = " << p4cortmp[0]
00875                 << endl;
00876 
00877     }
00878     ****/
00879 
00880     /****
00881     dijetMassCor1.Fill( LeadMass1 );    
00882     dijetMassCor_700_1.Fill( LeadMass1 );    
00883     dijetMassCor_2000_1.Fill( LeadMass1 );    
00884     dijetMassCor_5000_1.Fill( LeadMass1 );    
00885 
00886     if (EtaOk10 == 2) {
00887       dijetMassCor101.Fill( LeadMass1 );  
00888       dijetMassCor_700_101.Fill( LeadMass1 );  
00889       dijetMassCor_2000_101.Fill( LeadMass1 );  
00890       dijetMassCor_5000_101.Fill( LeadMass1 );  
00891     }
00892     if (EtaOk13 == 2) {
00893       dijetMassCor131.Fill( LeadMass1 );  
00894       dijetMassCor_700_131.Fill( LeadMass1 );  
00895       dijetMassCor_2000_131.Fill( LeadMass1 );  
00896       dijetMassCor_5000_131.Fill( LeadMass1 );  
00897     }
00898     if (EtaOk40 == 2) {
00899       dijetMassCor401.Fill( LeadMass1 ); 
00900       dijetMassCor_700_401.Fill( LeadMass1 ); 
00901       dijetMassCor_2000_401.Fill( LeadMass1 ); 
00902       dijetMassCor_5000_401.Fill( LeadMass1 ); 
00903     }
00904     ****/
00905 
00906   }
00907 
00908 
00909   // ********************************
00910   // **** Get the CaloJet2 collection
00911   // ********************************
00912   Handle<CaloJetCollection> caloJets2;
00913   evt.getByLabel( CaloJetAlgorithm2, caloJets2 );
00914 
00915   // Count Jets above Pt cut
00916   for (int istep = 0; istep < 100; ++istep) {
00917     int     njet = 0;
00918     float ptStep = (istep * (5000./100.));
00919 
00920     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00921       if ( cal->pt() > ptStep ) njet++;      
00922     
00923     hf_nJet2.Fill( ptStep, njet );
00924   }
00925 
00926   for (int istep = 0; istep < 100; ++istep) {
00927     int     njet = 0;
00928     float ptStep = (istep * (200./100.));
00929 
00930     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00931       if ( cal->pt() > ptStep ) njet++;      
00932     
00933     hf_nJet2s.Fill( ptStep, njet );
00934   }
00935 
00936 
00937   for (int istep = 0; istep < 100; ++istep) {
00938     int     njet = 0;
00939     float ptStep = (istep * (3000./100.));
00940 
00941     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00942       if ( cal->pt() > ptStep ) njet++;      
00943     
00944     hf_nJet21.Fill( ptStep, njet );
00945   }
00946 
00947 
00948 
00949 
00950 
00951   //Loop over the two leading CaloJets and fill some histograms
00952   jetInd = 0;
00953   allJetInd = 0;
00954   for( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal ) {
00955 
00956     allJetInd++;
00957     if (allJetInd == 1) {
00958       h_jet1Pt2.Fill( cal->pt() );
00959       pt2 = cal->pt();
00960       p4tmp[0] = cal->p4();
00961     }
00962     if (allJetInd == 2) {
00963       h_jet2Pt2.Fill( cal->pt() );
00964       p4tmp[1] = cal->p4();
00965     }
00966     if ( (allJetInd == 1) || (allJetInd == 2) ) {
00967       h_ptCalL2.Fill( cal->pt() );   
00968       h_ptCalL22.Fill( cal->pt() );   
00969       h_ptCalL23.Fill( cal->pt() );   
00970 
00971       h_etaCalL2.Fill( cal->eta() );
00972       h_phiCalL2.Fill( cal->phi() );
00973     }
00974     if (allJetInd == 3) h_jet3Pt2.Fill( cal->pt() );
00975     if (allJetInd == 4) h_jet4Pt2.Fill( cal->pt() );
00976     if (allJetInd == 5) h_jet5Pt2.Fill( cal->pt() );
00977     if (allJetInd == 6) h_jet6Pt2.Fill( cal->pt() );
00978     if (allJetInd == 7) h_jet7Pt2.Fill( cal->pt() );
00979 
00980     if ( fabs(cal->eta()) < 1.3) {
00981       h_lowPtCal21.Fill( cal->pt() );   
00982       if ( cal->pt() > 10.) {
00983         h_lowPtCal2c11.Fill( cal->pt() );   
00984       }
00985     }
00986     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
00987       h_lowPtCal22.Fill( cal->pt() );   
00988       if ( cal->pt() > 10.) {
00989         h_lowPtCal2c12.Fill( cal->pt() );   
00990       }
00991     }
00992     if ( fabs(cal->eta()) > 3.0) {
00993       h_lowPtCal23.Fill( cal->pt() );   
00994       if ( cal->pt() > 10.) {
00995         h_lowPtCal2c13.Fill( cal->pt() );   
00996       }
00997     }
00998 
00999 
01000     if ( cal->pt() > minJetPt) {
01001       h_ptCal2.Fill( cal->pt() );   
01002       h_ptCal22.Fill( cal->pt() );   
01003       h_ptCal23.Fill( cal->pt() );   
01004 
01005       h_etaCal2.Fill( cal->eta() );
01006       h_phiCal2.Fill( cal->phi() );
01007       jetInd++;
01008     }
01009   }
01010   //  h_nCalJets2.Fill( caloJets2->size() ); 
01011   h_nCalJets2.Fill( jetInd ); 
01012   if (jetInd > 1) {
01013     LeadMass2 = (p4tmp[0]+p4tmp[1]).mass();
01014     dijetMass2.Fill( LeadMass2 );
01015     dijetMass22.Fill( LeadMass2 );
01016     dijetMass23.Fill( LeadMass2 );
01017 
01018 
01019     dijetMassCor1.Fill( LeadMass2 );    
01020     dijetMassCor_700_1.Fill( LeadMass2 );    
01021     dijetMassCor_2000_1.Fill( LeadMass2 );    
01022     dijetMassCor_5000_1.Fill( LeadMass2 );    
01023 
01024     if (EtaOk10 == 2) {
01025       dijetMassCor101.Fill( LeadMass2 );  
01026       dijetMassCor_700_101.Fill( LeadMass2 );  
01027       dijetMassCor_2000_101.Fill( LeadMass2 );  
01028       dijetMassCor_5000_101.Fill( LeadMass2 );  
01029     }
01030     if (EtaOk13 == 2) {
01031       dijetMassCor131.Fill( LeadMass2 );  
01032       dijetMassCor_700_131.Fill( LeadMass2 );  
01033       dijetMassCor_2000_131.Fill( LeadMass2 );  
01034       dijetMassCor_5000_131.Fill( LeadMass2 );  
01035     }
01036     if (EtaOk40 == 2) {
01037       dijetMassCor401.Fill( LeadMass2 ); 
01038       dijetMassCor_700_401.Fill( LeadMass2 ); 
01039       dijetMassCor_2000_401.Fill( LeadMass2 ); 
01040       dijetMassCor_5000_401.Fill( LeadMass2 ); 
01041     }
01042 
01043 
01044   }
01045 
01046 
01047 
01048   // ********************************
01049   // **** Get the CaloJet3 collection
01050   // ********************************
01051   Handle<CaloJetCollection> caloJets3;
01052   evt.getByLabel( CaloJetAlgorithm3, caloJets3 );
01053 
01054   //Loop over the two leading CaloJets and fill some histograms
01055   jetInd = 0;
01056   allJetInd = 0;
01057 
01058   // Count Jets above Pt cut
01059   for (int istep = 0; istep < 100; ++istep) {
01060     int     njet = 0;
01061     float ptStep = (istep * (5000./100.));
01062 
01063     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01064       if ( cal->pt() > ptStep ) njet++;      
01065     
01066 
01067     hf_nJet3.Fill( ptStep, njet );
01068   }
01069 
01070   for (int istep = 0; istep < 100; ++istep) {
01071     int     njet = 0;
01072     float ptStep = (istep * (200./100.));
01073 
01074     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01075       if ( cal->pt() > ptStep ) njet++;          
01076 
01077     hf_nJet3s.Fill( ptStep, njet );
01078   }
01079 
01080   for (int istep = 0; istep < 100; ++istep) {
01081     int     njet = 0;
01082     float ptStep = (istep * (3000./100.));
01083 
01084     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01085       if ( cal->pt() > ptStep ) njet++;          
01086 
01087     hf_nJet31.Fill( ptStep, njet );
01088   }
01089 
01090 
01091   for( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal ) {
01092 
01093     allJetInd++;
01094     if (allJetInd == 1) {
01095       h_jet1Pt3.Fill( cal->pt() );
01096       pt3 = cal->pt();
01097       p4tmp[0] = cal->p4();
01098     }
01099     if (allJetInd == 2) {
01100       h_jet2Pt3.Fill( cal->pt() );
01101       p4tmp[1] = cal->p4();
01102     }
01103     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01104       h_ptCalL3.Fill( cal->pt() );   
01105       h_ptCalL32.Fill( cal->pt() );   
01106       h_ptCalL33.Fill( cal->pt() );   
01107 
01108       h_etaCalL3.Fill( cal->eta() );
01109       h_phiCalL3.Fill( cal->phi() );
01110     }
01111     if (allJetInd == 3) h_jet3Pt3.Fill( cal->pt() );
01112     if (allJetInd == 4) h_jet4Pt3.Fill( cal->pt() );
01113     if (allJetInd == 5) h_jet5Pt3.Fill( cal->pt() );
01114     if (allJetInd == 6) h_jet6Pt3.Fill( cal->pt() );
01115     if (allJetInd == 7) h_jet7Pt3.Fill( cal->pt() );
01116 
01117 
01118     if ( fabs(cal->eta()) < 1.3) {
01119       h_lowPtCal31.Fill( cal->pt() );   
01120       if ( cal->pt() > 10.) {
01121         h_lowPtCal3c11.Fill( cal->pt() );   
01122       }
01123     }
01124     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
01125       h_lowPtCal32.Fill( cal->pt() );   
01126       if ( cal->pt() > 10.) {
01127         h_lowPtCal3c12.Fill( cal->pt() );   
01128       }
01129     }
01130     if ( fabs(cal->eta()) > 3.0) {
01131       h_lowPtCal33.Fill( cal->pt() );   
01132       if ( cal->pt() > 10.) {
01133         h_lowPtCal3c13.Fill( cal->pt() );   
01134       }
01135     }
01136 
01137 
01138 
01139     if ( cal->pt() > minJetPt) {
01140       //    std::cout << "CALO JET3 #" << jetInd << std::endl << cal->print() << std::endl;
01141       h_ptCal3.Fill( cal->pt() );   
01142       h_ptCal32.Fill( cal->pt() );   
01143       h_ptCal33.Fill( cal->pt() );   
01144 
01145       h_etaCal3.Fill( cal->eta() );
01146       h_phiCal3.Fill( cal->phi() );
01147       jetInd++;
01148     }
01149   }
01150   //  h_nCalJets3.Fill( caloJets3->size() ); 
01151   h_nCalJets3.Fill( jetInd ); 
01152   if (jetInd > 1) {
01153     LeadMass3 = (p4tmp[0]+p4tmp[1]).mass();
01154     dijetMass3.Fill( LeadMass3 );
01155     dijetMass32.Fill( LeadMass3 );
01156     dijetMass33.Fill( LeadMass3 );
01157   }
01158 
01159 
01160 
01161   // ********************************
01162   // **** Get the CaloJet4 collection
01163   // ********************************
01164 
01165   Handle<CaloJetCollection> caloJets4;
01166   evt.getByLabel( CaloJetAlgorithm4, caloJets4 );
01167 
01168   //Loop over the two leading CaloJets and fill some histograms
01169   jetInd = 0;
01170   allJetInd = 0;
01171 
01172   // Count Jets above Pt cut
01173   for (int istep = 0; istep < 100; ++istep) {
01174     int     njet = 0;
01175     float ptStep = (istep * (5000./100.));
01176 
01177     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01178       if ( cal->pt() > ptStep ) njet++;      
01179     
01180 
01181     hf_nJet4.Fill( ptStep, njet );
01182   }
01183 
01184   for (int istep = 0; istep < 100; ++istep) {
01185     int     njet = 0;
01186     float ptStep = (istep * (200./100.));
01187 
01188     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01189       if ( cal->pt() > ptStep ) njet++;          
01190 
01191     hf_nJet4s.Fill( ptStep, njet );
01192   }
01193 
01194   for (int istep = 0; istep < 100; ++istep) {
01195     int     njet = 0;
01196     float ptStep = (istep * (3000./100.));
01197 
01198     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01199       if ( cal->pt() > ptStep ) njet++;          
01200 
01201     hf_nJet41.Fill( ptStep, njet );
01202   }
01203 
01204 
01205   for( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal ) {
01206 
01207     allJetInd++;
01208     if (allJetInd == 1) {
01209       h_jet1Pt4.Fill( cal->pt() );
01210       //      pt3 = cal->pt();
01211       p4tmp[0] = cal->p4();
01212     }
01213     if (allJetInd == 2) {
01214       h_jet2Pt4.Fill( cal->pt() );
01215       p4tmp[1] = cal->p4();
01216     }
01217     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01218       h_ptCalL4.Fill(  cal->pt() );   
01219       h_ptCalL42.Fill( cal->pt() );   
01220       h_ptCalL43.Fill( cal->pt() );   
01221 
01222       h_etaCalL4.Fill( cal->eta() );
01223       h_phiCalL4.Fill( cal->phi() );
01224     }
01225     if (allJetInd == 3) h_jet3Pt4.Fill( cal->pt() );
01226     if (allJetInd == 4) h_jet4Pt4.Fill( cal->pt() );
01227     if (allJetInd == 5) h_jet5Pt4.Fill( cal->pt() );
01228     if (allJetInd == 6) h_jet6Pt4.Fill( cal->pt() );
01229     if (allJetInd == 7) h_jet7Pt4.Fill( cal->pt() );
01230 
01231 
01232     if ( fabs(cal->eta()) < 1.3) {
01233       h_lowPtCal41.Fill( cal->pt() );   
01234       if ( cal->pt() > 10.) {
01235         h_lowPtCal4c11.Fill( cal->pt() );   
01236       }
01237     }
01238     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
01239       h_lowPtCal42.Fill( cal->pt() );   
01240       if ( cal->pt() > 10.) {
01241         h_lowPtCal4c12.Fill( cal->pt() );   
01242       }
01243     }
01244     if ( fabs(cal->eta()) > 3.0) {
01245       h_lowPtCal43.Fill( cal->pt() );   
01246       if ( cal->pt() > 10.) {
01247         h_lowPtCal4c13.Fill( cal->pt() );   
01248       }
01249     }
01250 
01251 
01252     if ( cal->pt() > minJetPt) {
01253       //    std::cout << "CALO JET4 #" << jetInd << std::endl << cal->print() << std::endl;
01254       h_ptCal4.Fill( cal->pt() );   
01255       h_ptCal42.Fill( cal->pt() );   
01256       h_ptCal43.Fill( cal->pt() );   
01257 
01258       h_etaCal4.Fill( cal->eta() );
01259       h_phiCal4.Fill( cal->phi() );
01260       jetInd++;
01261     }
01262   }
01263   //  h_nCalJets4.Fill( caloJets4->size() ); 
01264   h_nCalJets4.Fill( jetInd ); 
01265   if (jetInd > 1) {
01266     LeadMass4 = (p4tmp[0]+p4tmp[1]).mass();
01267     dijetMass4.Fill( LeadMass4 );
01268     dijetMass42.Fill( LeadMass4 );
01269     dijetMass43.Fill( LeadMass4 );
01270   }
01271 
01272 
01273 
01274   // *********************************************
01275   // *********************************************
01276 
01277 
01278   //**** Get the GenJet1 collection
01279   Handle<GenJetCollection> genJets1;
01280   evt.getByLabel( GenJetAlgorithm1, genJets1 );
01281 
01282 
01283   //Loop over the two leading GenJets and fill some histograms
01284   jetInd    = 0;
01285   allJetInd = 0;
01286   for( GenJetCollection::const_iterator gen = genJets1->begin(); gen != genJets1->end(); ++ gen ) {
01287     allJetInd++;
01288     if (allJetInd == 1) {
01289       p4tmp[0] = gen->p4();
01290     }
01291     if (allJetInd == 2) {
01292       p4tmp[1] = gen->p4();
01293     }
01294 
01295     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01296       h_ptGenL1.Fill( gen->pt() );   
01297       h_ptGenL12.Fill( gen->pt() );   
01298       h_ptGenL13.Fill( gen->pt() );   
01299 
01300       h_etaGenL1.Fill( gen->eta() );
01301       h_phiGenL1.Fill( gen->phi() );
01302     }
01303 
01304     if ( gen->pt() > minJetPt) {
01305       // std::cout << "GEN JET1 #" << jetInd << std::endl << gen->print() << std::endl;
01306       h_ptGen1.Fill( gen->pt() );   
01307       h_ptGen12.Fill( gen->pt() );   
01308       h_ptGen13.Fill( gen->pt() );   
01309 
01310       h_etaGen1.Fill( gen->eta() );
01311       h_phiGen1.Fill( gen->phi() );
01312       jetInd++;
01313     }
01314   }
01315 
01316   LeadMassP1 = (p4tmp[0]+p4tmp[1]).mass();
01317   dijetMassP1.Fill( LeadMassP1 ); 
01318   if (EtaOk10 == 2) {
01319     dijetMassP101.Fill( LeadMassP1 ); 
01320     dijetMassP_700_101.Fill( LeadMassP1 ); 
01321     dijetMassP_2000_101.Fill( LeadMassP1 ); 
01322     dijetMassP_5000_101.Fill( LeadMassP1 ); 
01323    }
01324   if (EtaOk13 == 2) {
01325     dijetMassP131.Fill( LeadMassP1 ); 
01326     dijetMassP_700_131.Fill( LeadMassP1 ); 
01327     dijetMassP_2000_131.Fill( LeadMassP1 ); 
01328     dijetMassP_5000_131.Fill( LeadMassP1 ); 
01329 
01330    }
01331   if (EtaOk40 == 2) {
01332     dijetMassP401.Fill( LeadMassP1 );
01333     dijetMassP_5000_401.Fill( LeadMassP1 );
01334     dijetMassP_5000_401.Fill( LeadMassP1 );
01335     dijetMassP_5000_401.Fill( LeadMassP1 );
01336   } 
01337 
01338   //  h_nGenJets1.Fill( genJets1->size() ); 
01339   h_nGenJets1.Fill( jetInd ); 
01340 
01341   //**** Get the GenJet2 collection
01342   Handle<GenJetCollection> genJets2;
01343   evt.getByLabel( GenJetAlgorithm2, genJets2 );
01344 
01345   //Loop over the two leading GenJets and fill some histograms
01346   jetInd    = 0;
01347   allJetInd = 0;
01348   for( GenJetCollection::const_iterator gen = genJets2->begin(); gen != genJets2->end(); ++ gen ) {
01349     allJetInd++;
01350     if (allJetInd == 1) {
01351       p4tmp[0] = gen->p4();
01352     }
01353     if (allJetInd == 2) {
01354       p4tmp[1] = gen->p4();
01355     }
01356     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01357       h_ptGenL2.Fill( gen->pt() );   
01358       h_ptGenL22.Fill( gen->pt() );   
01359       h_ptGenL23.Fill( gen->pt() );   
01360 
01361       h_etaGenL2.Fill( gen->eta() );
01362       h_phiGenL2.Fill( gen->phi() );
01363     }
01364 
01365     if ( gen->pt() > minJetPt) {
01366       // std::cout << "GEN JET2 #" << jetInd << std::endl << gen->print() << std::endl;
01367       h_ptGen2.Fill( gen->pt() );   
01368       h_ptGen22.Fill( gen->pt() );   
01369       h_ptGen23.Fill( gen->pt() );   
01370 
01371       h_etaGen2.Fill( gen->eta() );
01372       h_phiGen2.Fill( gen->phi() );
01373       jetInd++;
01374     }
01375   }
01376   
01377   LeadMassP2 = (p4tmp[0]+p4tmp[1]).mass();
01378   dijetMassP2.Fill( LeadMassP2 ); 
01379 
01380 
01381   //  h_nGenJets2.Fill( genJets2->size() ); 
01382   h_nGenJets2.Fill( jetInd ); 
01383 
01384   //**** Get the GenJet3 collection
01385   Handle<GenJetCollection> genJets3;
01386   evt.getByLabel( GenJetAlgorithm3, genJets3 );
01387 
01388   //Loop over the two leading GenJets and fill some histograms
01389   jetInd    = 0;
01390   allJetInd = 0;
01391   for( GenJetCollection::const_iterator gen = genJets3->begin(); gen != genJets3->end(); ++ gen ) {
01392     allJetInd++;
01393     if (allJetInd == 1) {
01394       p4tmp[0] = gen->p4();
01395     }
01396     if (allJetInd == 2) {
01397       p4tmp[1] = gen->p4();
01398     }    
01399     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01400       h_ptGenL3.Fill( gen->pt() );   
01401       h_ptGenL32.Fill( gen->pt() );   
01402       h_ptGenL33.Fill( gen->pt() );   
01403 
01404       h_etaGenL3.Fill( gen->eta() );
01405       h_phiGenL3.Fill( gen->phi() );
01406     }
01407 
01408     if ( gen->pt() > minJetPt) {
01409       // std::cout << "GEN JET3 #" << jetInd << std::endl << gen->print() << std::endl;
01410       h_ptGen3.Fill( gen->pt() );   
01411       h_ptGen32.Fill( gen->pt() );   
01412       h_ptGen33.Fill( gen->pt() );   
01413 
01414       h_etaGen3.Fill( gen->eta() );
01415       h_phiGen3.Fill( gen->phi() );
01416       jetInd++;
01417     }
01418   }
01419 
01420   LeadMassP3 = (p4tmp[0]+p4tmp[1]).mass();
01421   dijetMassP3.Fill( LeadMassP3 ); 
01422 
01423 
01424   //  h_nGenJets3.Fill( genJets3->size() ); 
01425   h_nGenJets3.Fill( jetInd ); 
01426 
01427 
01428   // *********************
01429   // MidPoint Jet Matching
01430   
01431   Handle<GenJetCollection>  genJets;
01432   Handle<CaloJetCollection> caloJets;
01433 
01434   //  evt.getByLabel( "midPointCone5GenJets",  genJets );
01435   //  evt.getByLabel( "midPointCone5CaloJets", caloJets );
01436   evt.getByLabel( GenJetAlgorithm1, genJets );
01437   evt.getByLabel( CaloJetAlgorithm1, caloJets );
01438 
01439 
01440   int maxJets = MAXJETS;
01441 
01442   jetInd = 0;
01443   double dRmin[MAXJETS];
01444   math::XYZTLorentzVector  p4gen[MAXJETS], p4cal[MAXJETS],  
01445     p4par[MAXJETS], p4Zp[MAXJETS], p4part[MAXJETS];
01446 
01447   int used[MAXJETS];
01448   int nj;
01449 
01450   for( int i=0; i<maxJets; ++i ) used[i] = 0;  
01451 
01452   //  cout << ">>>>>>>>> " << endl;
01453 
01454 
01455 
01456   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01457        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01458 
01459     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01460     dRmin[jetInd] = 1000.0;
01461 
01462     nj      = 0;
01463     usedInd = -1;
01464     
01465     for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) { 
01466  
01467       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01468  
01469       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01470         dRmin[jetInd] = delR;        // delta R of match
01471         p4cal[jetInd] = cal->p4();   // Matched Cal 4-vector
01472         usedInd       = nj;     
01473       }
01474 
01475       nj++;
01476     }
01477 
01478     if (fabs(p4gen[jetInd].eta()) < 1.3) {
01479       matchedAllPt11.Fill(p4gen[jetInd].pt());
01480     }
01481     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01482       matchedAllPt12.Fill(p4gen[jetInd].pt());
01483     }
01484     if (fabs(p4gen[jetInd].eta()) > 3.) {
01485       matchedAllPt13.Fill(p4gen[jetInd].pt());     
01486     }
01487 
01488     if (usedInd != -1) {
01489 
01490       used[usedInd] = 1;
01491 
01492       if (p4cal[jetInd].pt() > minJetPt10)
01493         hf_PtResponse1.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01494 
01495       if (fabs(p4gen[jetInd].eta()) < 1.3) {
01496         matchedPt11.Fill(p4gen[jetInd].pt());
01497       } 
01498       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01499         matchedPt12.Fill(p4gen[jetInd].pt());
01500       }
01501       if (fabs(p4gen[jetInd].eta()) > 3.) {
01502         matchedPt13.Fill(p4gen[jetInd].pt());
01503       }
01504       /***
01505       cout << " >>>SISCone "     
01506            << jetInd                 << " " 
01507            << p4gen[jetInd].pt()     << " " 
01508            << p4gen[jetInd].eta()    << " " 
01509            << p4gen[jetInd].phi()    << " " 
01510            << p4cal[jetInd].pt()     << " " 
01511            << p4cal[jetInd].eta()    << " " 
01512            << p4cal[jetInd].phi()    
01513            << endl;
01514       ***/
01515 
01516       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01517       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01518         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01519           dPt20Frac1.Fill(dpt/p4gen[jetInd].pt());
01520         }
01521         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01522           dPt40Frac1.Fill(dpt/p4gen[jetInd].pt());
01523         }
01524         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01525           dPt80Frac1.Fill(dpt/p4gen[jetInd].pt());
01526         }
01527         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01528           dPt100Frac1.Fill(dpt/p4gen[jetInd].pt());
01529         }
01530       }
01531       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) ) {
01532 
01533         dR1.Fill(dRmin[jetInd]);
01534         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01535         dPhi1.Fill(dphi);
01536         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01537         dEta1.Fill(deta);
01538         dPt1.Fill(dpt);         
01539         dPtFrac1.Fill(dpt/p4gen[jetInd].pt());
01540 
01541         if ( ( (dpt/p4gen[jetInd].pt()) < -0.5 ) &&  ( fabs(dpt) > 90. ) ) {
01542 
01543           cout << " deltaR min = "     << dRmin[jetInd]
01544                << " Ind = " << jetInd  << " / " << usedInd << " / " << used[nj]
01545                << " Del pt / frac  = " << dpt
01546                << " / "                << dpt/p4gen[jetInd].pt()
01547                << " cal/gen pt   = "   << p4cal[jetInd].pt()
01548                << " / "                << p4gen[jetInd].pt()
01549                << " cal/gen eta  = "   << p4cal[jetInd].eta()
01550                << " / "                << p4gen[jetInd].eta()
01551                << " cal/gen phi  = "   << p4cal[jetInd].phi()
01552                << " / "                << p4gen[jetInd].phi()
01553                << endl;
01554         }
01555 
01556       }
01557       
01558       jetInd++;    
01559 
01560     }
01561 
01562   }
01563 
01564 
01565 
01566   // *********************
01567   // Seedless Jet Matching
01568   
01569   //  Handle<GenJetCollection>  genJets;
01570   //  Handle<CaloJetCollection> caloJets;
01571 
01572   //  evt.getByLabel( "sisCone5GenJets",  genJets );
01573   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
01574   evt.getByLabel( GenJetAlgorithm2, genJets );
01575   evt.getByLabel( CaloJetAlgorithm2, caloJets );
01576 
01577   // int maxJets = 20;
01578   jetInd = 0;
01579   //  double dRmin[20];
01580   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01581 
01582   for( int i=0; i<maxJets; ++i ) used[i] = 0;    
01583   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01584        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01585     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01586     dRmin[jetInd] = 1000.0;
01587 
01588     nj      = 0;
01589     usedInd = -1;
01590 
01591     for( CaloJetCollection::const_iterator cal = caloJets->begin(); 
01592          cal != caloJets->end(); ++ cal ) { 
01593       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01594 
01595       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01596         dRmin[jetInd] = delR;       // delta R of match
01597         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
01598         usedInd       = nj;
01599       }
01600       nj++;
01601     }
01602 
01603     if (fabs(p4gen[jetInd].eta()) < 1.3)   {      
01604       matchedAllPt21.Fill(p4gen[jetInd].pt());
01605     }
01606     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01607       matchedAllPt22.Fill(p4gen[jetInd].pt());
01608     }
01609     if (fabs(p4gen[jetInd].eta()) > 3.) {
01610       matchedAllPt23.Fill(p4gen[jetInd].pt());
01611     }
01612 
01613 
01614     if (usedInd != -1) {
01615 
01616       used[usedInd] = 1;
01617 
01618       if (p4cal[jetInd].pt() > minJetPt10)
01619         hf_PtResponse2.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01620 
01621       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01622         matchedPt21.Fill(p4gen[jetInd].pt());
01623       }
01624       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01625         matchedPt22.Fill(p4gen[jetInd].pt());
01626       }
01627       if (fabs(p4gen[jetInd].eta()) > 3.) {
01628         matchedPt23.Fill(p4gen[jetInd].pt());
01629       }
01630       /***
01631       cout << " >>>IterCone "     
01632            << jetInd                 << " " 
01633            << p4gen[jetInd].pt()     << " " 
01634            << p4gen[jetInd].eta()    << " " 
01635            << p4gen[jetInd].phi()    << " " 
01636            << p4cal[jetInd].pt()     << " " 
01637            << p4cal[jetInd].eta()    << " " 
01638            << p4cal[jetInd].phi()    
01639            << endl;
01640       ***/
01641 
01642 
01643       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01644       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01645         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01646           dPt20Frac2.Fill(dpt/p4gen[jetInd].pt());
01647         }
01648         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01649           dPt40Frac2.Fill(dpt/p4gen[jetInd].pt());
01650         }
01651         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01652           dPt80Frac2.Fill(dpt/p4gen[jetInd].pt());
01653         }
01654         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01655           dPt100Frac2.Fill(dpt/p4gen[jetInd].pt());
01656         }
01657       }
01658       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) )  {
01659 
01660         dR2.Fill(dRmin[jetInd]);
01661         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01662         dPhi2.Fill(dphi);
01663         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01664         dEta2.Fill(deta);
01665         dPt2.Fill(dpt);
01666         dPtFrac2.Fill(dpt/p4gen[jetInd].pt());       
01667 
01668       }      
01669 
01670       jetInd++;    
01671     }
01672 
01673   }
01674 
01675   // *********************
01676   // Kt Jet Matching
01677   
01678   //  Handle<GenJetCollection>  genJets;
01679   //  Handle<CaloJetCollection> caloJets;
01680 
01681   //  evt.getByLabel( "sisCone5GenJets",  genJets );
01682   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
01683   evt.getByLabel( GenJetAlgorithm3, genJets );
01684   evt.getByLabel( CaloJetAlgorithm3, caloJets );
01685 
01686   // int maxJets = 20;
01687   jetInd = 0;
01688   //  double dRmin[20];
01689   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01690 
01691   for( int i=0; i<maxJets; ++i ) used[i] = 0;
01692   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01693        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01694     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01695     dRmin[jetInd] = 1000.0;
01696 
01697     nj = 0;
01698     usedInd = -1;
01699 
01700     for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) { 
01701       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01702 
01703       if ( (delR<dRmin[jetInd]) &&  (delR < matchedDelR) && (used[nj] == 0) ) {
01704         dRmin[jetInd] = delR;        // delta R of match
01705         p4cal[jetInd] = cal->p4();   // Matched Cal 4-vector
01706         usedInd       = nj;
01707       }
01708       nj++;
01709     }
01710 
01711     if (fabs(p4gen[jetInd].eta()) < 1.3) {
01712       matchedAllPt31.Fill(p4gen[jetInd].pt());
01713     }
01714     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01715       matchedAllPt32.Fill(p4gen[jetInd].pt());
01716     }
01717    if (fabs(p4gen[jetInd].eta()) > 3.) {
01718       matchedAllPt33.Fill(p4gen[jetInd].pt());
01719     }
01720 
01721 
01722     if (usedInd != -1) {
01723       used[usedInd] = 1;
01724 
01725       if (p4cal[jetInd].pt() > minJetPt10)
01726         hf_PtResponse3.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01727 
01728       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01729         matchedPt31.Fill(p4gen[jetInd].pt());
01730       }
01731       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01732         matchedPt32.Fill(p4gen[jetInd].pt());
01733       }
01734       if (fabs(p4gen[jetInd].eta()) > 3.) {
01735         matchedPt33.Fill(p4gen[jetInd].pt());
01736       }
01737       /***
01738       cout << " >>>MidPoint "     
01739            << jetInd                 << " " 
01740            << p4gen[jetInd].pt()     << " " 
01741            << p4gen[jetInd].eta()    << " " 
01742            << p4gen[jetInd].phi()    << " " 
01743            << p4cal[jetInd].pt()     << " " 
01744            << p4cal[jetInd].eta()    << " " 
01745            << p4cal[jetInd].phi()    
01746            << endl;
01747       ***/
01748 
01749       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01750       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01751         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01752           dPt20Frac3.Fill(dpt/p4gen[jetInd].pt());
01753         }
01754         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01755           dPt40Frac3.Fill(dpt/p4gen[jetInd].pt());
01756         }
01757         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01758           dPt80Frac3.Fill(dpt/p4gen[jetInd].pt());
01759         }
01760         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01761           dPt100Frac3.Fill(dpt/p4gen[jetInd].pt());
01762         }
01763       }
01764 
01765       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) ) {
01766 
01767         dR3.Fill(dRmin[jetInd]);
01768         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01769         dPhi3.Fill(dphi);
01770         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01771         dEta3.Fill(deta);
01772         dPt3.Fill(dpt);
01773         dPtFrac3.Fill(dpt/p4gen[jetInd].pt());
01774 
01775 
01776       }
01777       
01778       jetInd++;    
01779     }
01780   }
01781 
01782   // *********************
01783   // Jet Algorithm 4
01784   
01785   evt.getByLabel( GenJetAlgorithm4,  genJets );
01786   evt.getByLabel( CaloJetAlgorithm4, caloJets );
01787 
01788   // int maxJets = 20;
01789   jetInd = 0;
01790   //  double dRmin[20];
01791   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01792 
01793   for( int i=0; i<maxJets; ++i ) used[i] = 0;    
01794   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01795        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01796     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01797     dRmin[jetInd] = 1000.0;
01798 
01799     nj      = 0;
01800     usedInd = -1;
01801 
01802     for( CaloJetCollection::const_iterator cal = caloJets->begin(); 
01803          cal != caloJets->end(); ++ cal ) { 
01804       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01805 
01806       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01807         dRmin[jetInd] = delR;       // delta R of match
01808         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
01809         usedInd       = nj;
01810       }
01811       nj++;
01812     }
01813 
01814     if (fabs(p4gen[jetInd].eta()) < 1.3)   {      
01815       matchedAllPt41.Fill(p4gen[jetInd].pt());
01816     }
01817     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01818       matchedAllPt42.Fill(p4gen[jetInd].pt());
01819     }
01820     if (fabs(p4gen[jetInd].eta()) > 3.) {
01821       matchedAllPt43.Fill(p4gen[jetInd].pt());
01822     }
01823 
01824 
01825     if (usedInd != -1) {
01826 
01827       used[usedInd] = 1;
01828 
01829       if (p4cal[jetInd].pt() > minJetPt10)
01830         hf_PtResponse4.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01831 
01832       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01833         matchedPt41.Fill(p4gen[jetInd].pt());
01834       }
01835       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01836         matchedPt42.Fill(p4gen[jetInd].pt());
01837       }
01838       if (fabs(p4gen[jetInd].eta()) > 3.) {
01839         matchedPt43.Fill(p4gen[jetInd].pt());
01840       }
01841 
01842       /***
01843       cout << " >>>KtClus "     
01844            << jetInd                 << " " 
01845            << p4gen[jetInd].pt()     << " " 
01846            << p4gen[jetInd].eta()    << " " 
01847            << p4gen[jetInd].phi()    << " " 
01848            << p4cal[jetInd].pt()     << " " 
01849            << p4cal[jetInd].eta()    << " " 
01850            << p4cal[jetInd].phi()    
01851            << endl;
01852       ***/
01853 
01854 
01855       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01856       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01857         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01858           dPt20Frac4.Fill(dpt/p4gen[jetInd].pt());
01859         }
01860         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01861           dPt40Frac4.Fill(dpt/p4gen[jetInd].pt());
01862         }
01863         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01864           dPt80Frac4.Fill(dpt/p4gen[jetInd].pt());
01865         }
01866         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01867           dPt100Frac4.Fill(dpt/p4gen[jetInd].pt());
01868         }
01869       }
01870 
01871       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) )  {
01872 
01873         dR4.Fill(dRmin[jetInd]);
01874         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01875         dPhi4.Fill(dphi);
01876         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01877         dEta4.Fill(deta);
01878         dPt4.Fill(dpt);
01879         dPtFrac4.Fill(dpt/p4gen[jetInd].pt());       
01880 
01881       }      
01882 
01883       jetInd++;    
01884     }
01885 
01886   }
01887 
01888 
01889 
01890   // *********************
01891   // MidPoint - Seedless Jet Matching
01892   
01893   Handle<CaloJetCollection> calo1Jets;
01894   Handle<CaloJetCollection> calo2Jets;
01895   Handle<CaloJetCollection> calo3Jets;
01896 
01897   evt.getByLabel( CaloJetAlgorithm1, calo1Jets );
01898   evt.getByLabel( CaloJetAlgorithm2, calo2Jets );
01899   evt.getByLabel( CaloJetAlgorithm3, calo3Jets );
01900 
01901   jetInd = 0;
01902 
01903   for( int i=0; i<maxJets; ++i ) used[i] = 0;
01904   for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
01905        cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
01906 
01907     p4gen[jetInd] = cal1->p4();    //Gen 4-vector
01908     dRmin[jetInd] = 1000.0;
01909 
01910     nj = 0;
01911     for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); cal2 != calo2Jets->end(); ++cal2 ) { 
01912 
01913       double delR = deltaR( cal1->eta(), cal1->phi(), cal2->eta(), cal2->phi() ); 
01914       if ( (delR<dRmin[jetInd]) && (used[nj] == 0) ) {
01915         dRmin[jetInd] = delR;        // delta R of match
01916         p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
01917         usedInd       = nj;
01918       }
01919       nj++;
01920     }
01921     used[usedInd] = 1;
01922 
01923     if (p4gen[jetInd].pt() > minJetPt) {
01924       dR12.Fill(dRmin[jetInd]);
01925       double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01926       dPhi12.Fill(dphi);
01927       double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01928       dEta12.Fill(deta);
01929       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01930       dPt12.Fill(dpt);
01931     }
01932 
01933     jetInd++;    
01934   }
01935 
01936   // ******************************************
01937   // ******************************************
01938 
01939 
01940   Handle<CandidateCollection> genParticles;
01941   evt.getByLabel("genParticleCandidates",genParticles);
01942 
01943 
01944   // *********************
01945   // Partons (Z')
01946 
01947   int nPart = 0;
01948   for (size_t i =0;i< genParticles->size(); i++) {
01949 
01950     const Candidate &p = (*genParticles)[i];
01951     //    int Status =  p.status();
01952     //    bool ParticleIsStable = Status==1;
01953     int id = p.pdgId();
01954 
01955     if (id == 32) {
01956 
01957       if (p.numberOfDaughters() != 0) {
01958         /***
01959         cout << "Z': part = "   << i << " id = " << id 
01960              << " daughters = " << p.numberOfDaughters() 
01961              << " mass = "      << p.mass()
01962              << endl;
01963         ***/
01964         ZpMG =  p.mass();
01965         ZpMassGen.Fill( ZpMG );
01966         if (EtaOk10 == 2) {
01967           ZpMassGen10.Fill( ZpMG );
01968           ZpMassGen_700_10.Fill( ZpMG );
01969           ZpMassGen_2000_10.Fill( ZpMG );
01970           ZpMassGen_5000_10.Fill( ZpMG );
01971         }
01972         if (EtaOk13 == 2) {
01973           ZpMassGen13.Fill( ZpMG );
01974           ZpMassGen_700_13.Fill( ZpMG );
01975           ZpMassGen_2000_13.Fill( ZpMG );
01976           ZpMassGen_5000_13.Fill( ZpMG );
01977         }
01978         if (EtaOk40 == 2) {
01979           ZpMassGen40.Fill( ZpMG );
01980           ZpMassGen_700_40.Fill( ZpMG );
01981           ZpMassGen_2000_40.Fill( ZpMG );
01982           ZpMassGen_5000_40.Fill( ZpMG );
01983         }
01984       }
01985 
01986       for( int id1=0, nd1=p.numberOfDaughters(); id1 < nd1; ++id1 ) {
01987 
01988         const Candidate * d1 = p.daughter(id1);
01989 
01990         if ( abs(d1->pdgId()) != 32) {  
01991           math::XYZTLorentzVector momentum=d1->p4();
01992           p4Zp[nPart] = momentum=d1->p4();
01993           nPart++;
01994         }
01995 
01996       }
01997     }
01998 
01999   }
02000 
02001   // *********************
02002   // Match jets to Zp
02003   int genInd = 0;
02004 
02005   if (nPart == 2) {
02006     
02007     ZpM = (p4Zp[0]+p4Zp[1]).mass();
02008     ZpMass.Fill( ZpM );
02009 
02010     if (EtaOk10 == 2) {
02011       ZpMass_700_10.Fill( ZpM );
02012       ZpMass_2000_10.Fill( ZpM );
02013       ZpMass_5000_10.Fill( ZpM );
02014     }
02015     if (EtaOk13 == 2) {
02016       ZpMass_700_13.Fill( ZpM );
02017       ZpMass_2000_13.Fill( ZpM );
02018       ZpMass_5000_13.Fill( ZpM );
02019     }
02020     if (EtaOk40 == 2) {
02021       ZpMass_700_40.Fill( ZpM );
02022       ZpMass_2000_40.Fill( ZpM );
02023       ZpMass_5000_40.Fill( ZpM );
02024     }
02025 
02026     int usedInd;   
02027 
02028     // ***********
02029     // **** Calor1
02030     usedInd = -1;
02031     jetInd  = 0;
02032 
02033     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02034     for( int i=0; i<2; ++i ) { 
02035       
02036       dRmin[jetInd]  = 1000.0;
02037 
02038       int nj = 0;
02039       for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
02040            cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
02041         
02042         double delR = deltaR( cal1->eta(), cal1->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02043 
02044         //      if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02045         if ( (delR < dRmin[jetInd]) && (used[nj] == 0) ) {
02046           dRmin[jetInd] = delR;        // delta R of match
02047           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
02048           usedInd       = nj;
02049           genInd        = i;
02050         }
02051 
02052         /****
02053         cout << "Delta R = " << delR
02054              << " deltaR min = " << dRmin[jetInd]
02055              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02056              << " cal1 eta  = " << cal1->eta()
02057              << " p4par eta = " << p4Zp[i].eta()
02058              << " cal1 phi  = " << cal1->phi()
02059              << " p4par phi = " << p4Zp[i].phi()
02060              << endl;
02061         cout << "    " 
02062              << " p4par = " << p4Zp[i]
02063              << " p4cal = " << cal1->p4()
02064              << endl;
02065         ***/
02066                 
02067         nj++;
02068       }
02069 
02070       // Found matched jet
02071       if (usedInd != -1) {
02072         used[usedInd] = 1;
02073         jetInd++;    
02074       }
02075 
02076     }
02077     
02078     ZpMM = (p4cal[0]+p4cal[1]).mass();
02079     ZpMassMatched1.Fill( ZpMM );
02080 
02081     if ((ZpMG != 0) && (EtaOk40 == 2)) {
02082       ZpMassRes401.Fill( (ZpMM - ZpMG) / ZpMG );
02083 
02084       ZpMassResL401.Fill( (LeadMass1 - ZpMG) / ZpMG );
02085       ZpMassResL402.Fill( (LeadMass2 - ZpMG) / ZpMG );
02086       ZpMassResL403.Fill( (LeadMass3 - ZpMG) / ZpMG );
02087       
02088       ZpMassResRL401.Fill( LeadMass1 / ZpMG );
02089       ZpMassResRL402.Fill( LeadMass2 / ZpMG );
02090       ZpMassResRL403.Fill( LeadMass3 / ZpMG );
02091 
02092       ZpMassResRLoP401.Fill( LeadMass1 / LeadMassP1 );
02093       ZpMassResRLoP402.Fill( LeadMass2 / LeadMassP2 );
02094       ZpMassResRLoP403.Fill( LeadMass3 / LeadMassP2 );
02095 
02096       ZpMassResPRL401.Fill( LeadMassP1 / ZpMG );
02097       ZpMassResPRL402.Fill( LeadMassP2 / ZpMG );
02098       ZpMassResPRL403.Fill( LeadMassP3 / ZpMG );
02099       
02100     }
02101 
02102     if ((ZpMG != 0) && (EtaOk10 == 2)) {
02103       ZpMassRes101.Fill( (ZpMM - ZpMG) / ZpMG );
02104 
02105       ZpMassResL101.Fill( (LeadMass1 - ZpMG) / ZpMG );
02106       ZpMassResL102.Fill( (LeadMass2 - ZpMG) / ZpMG );
02107       ZpMassResL103.Fill( (LeadMass3 - ZpMG) / ZpMG );
02108       
02109       ZpMassResRL101.Fill( LeadMass1 / ZpMG );
02110       ZpMassResRL102.Fill( LeadMass2 / ZpMG );
02111       ZpMassResRL103.Fill( LeadMass3 / ZpMG );
02112 
02113       ZpMassResRLoP101.Fill( LeadMass1 / LeadMassP1 );
02114       ZpMassResRLoP102.Fill( LeadMass2 / LeadMassP2 );
02115       ZpMassResRLoP103.Fill( LeadMass3 / LeadMassP2 );
02116 
02117       ZpMassResPRL101.Fill( LeadMassP1 / ZpMG );
02118       ZpMassResPRL102.Fill( LeadMassP2 / ZpMG );
02119       ZpMassResPRL103.Fill( LeadMassP3 / ZpMG );
02120       
02121     }
02122 
02123     if ((ZpMG != 0) && (EtaOk13 == 2)) {
02124       ZpMassRes131.Fill( (ZpMM - ZpMG) / ZpMG );
02125 
02126       ZpMassResL131.Fill( (LeadMass1 - ZpMG) / ZpMG );
02127       ZpMassResL132.Fill( (LeadMass2 - ZpMG) / ZpMG );
02128       ZpMassResL133.Fill( (LeadMass3 - ZpMG) / ZpMG );
02129       
02130       ZpMassResRL131.Fill( LeadMass1 / ZpMG );
02131       ZpMassResRL132.Fill( LeadMass2 / ZpMG );
02132       ZpMassResRL133.Fill( LeadMass3 / ZpMG );
02133 
02134       ZpMassResRLoP131.Fill( LeadMass1 / LeadMassP1 );
02135       ZpMassResRLoP132.Fill( LeadMass2 / LeadMassP2 );
02136       ZpMassResRLoP133.Fill( LeadMass3 / LeadMassP2 );
02137 
02138       ZpMassResPRL131.Fill( LeadMassP1 / ZpMG );
02139       ZpMassResPRL132.Fill( LeadMassP2 / ZpMG );
02140       ZpMassResPRL133.Fill( LeadMassP3 / ZpMG );
02141       
02142     }
02143 
02144 
02145 
02146     // ***********
02147     // **** Calor2
02148     usedInd = -1;
02149     jetInd  = 0;
02150 
02151     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02152     for( int i=0; i<2; ++i ) { 
02153       
02154       dRmin[jetInd]  = 1000.0;
02155 
02156       int nj = 0;
02157       for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); 
02158            cal2 != calo2Jets->end() && jetInd<maxJets; ++cal2 ) { 
02159         
02160         double delR = deltaR( cal2->eta(), cal2->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02161 
02162         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02163           dRmin[jetInd] = delR;        // delta R of match
02164           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
02165           usedInd       = nj;
02166         }
02167 
02168         /****   
02169         cout << "Delta R = " << delR
02170              << " deltaR min = " << dRmin[jetInd]
02171              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02172              << " p4par = " << p4par[i]
02173              << " p4cal = " << cal1->p4()
02174              << " cal1 eta  = " << cal1->eta()
02175              << " p4par eta = " << p4par[i].eta()
02176              << endl;
02177         ****/
02178                 
02179         nj++;
02180       }
02181 
02182       // Found matched jet
02183       if (usedInd != -1) {
02184         used[usedInd] = 1;
02185         jetInd++;    
02186       }
02187 
02188     }
02189     
02190     ZpMM = (p4cal[0]+p4cal[1]).mass();
02191     ZpMassMatched2.Fill( ZpMM );
02192     ZpMassRes402.Fill( (ZpMM - ZpM) / ZpM );
02193 
02194 
02195     // ***********
02196     // **** Calor3
02197     usedInd = -1;
02198     jetInd  = 0;
02199 
02200     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02201     for( int i=0; i<2; ++i ) { 
02202       
02203       dRmin[jetInd]  = 1000.0;
02204 
02205       int nj = 0;
02206       for( CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); 
02207            cal3 != calo3Jets->end() && jetInd<maxJets; ++cal3 ) { 
02208         
02209         double delR = deltaR( cal3->eta(), cal3->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02210 
02211         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02212           dRmin[jetInd] = delR;        // delta R of match
02213           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
02214           usedInd       = nj;
02215         }
02216 
02217         /****   
02218         cout << "Delta R = " << delR
02219              << " deltaR min = " << dRmin[jetInd]
02220              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02221              << " p4par = " << p4par[i]
02222              << " p4cal = " << cal1->p4()
02223              << " cal1 eta  = " << cal1->eta()
02224              << " p4par eta = " << p4par[i].eta()
02225              << endl;
02226         ****/
02227                 
02228         nj++;
02229       }
02230 
02231 
02232       // Found matched jet
02233       if (usedInd != -1) {
02234         used[usedInd] = 1;
02235         jetInd++;    
02236       }
02237 
02238     }
02239     
02240     ZpMM = (p4cal[0]+p4cal[1]).mass();
02241     ZpMassMatched3.Fill( ZpMM );
02242     ZpMassRes403.Fill( (ZpMM - ZpM) / ZpM );
02243     
02244   } else {
02245     cout << "Z' (3): nPart = " << nPart << endl;
02246   }
02247 
02248 
02249 
02250   // *********************
02251   // Partons (ttbar) Jet Matching
02252 
02253   //  cout << ">>> Begin MC list" << endl;
02254   int nJet = 0;
02255 
02256 
02257   int ii = 1;
02258   int jj = 4;
02259 
02260   for (size_t i =0;i< genParticles->size(); i++) {
02261 
02262     const Candidate &p = (*genParticles)[i];
02263     //    int Status =  p.status();
02264     //    bool ParticleIsStable = Status==1;
02265     int id = p.pdgId();
02266 
02267     // Top Quarks
02268     if (abs(id) == 6) {
02269       cout << "TOP: id = " << id << " mass = " << p.mass() << endl;
02270 
02271       topMassParton.Fill(p.mass());
02272 
02273       if (id == 6)  tMassGen.Fill(p.mass());
02274       if (id == -6) tbarMassGen.Fill(p.mass());
02275 
02276       for( size_t id1=0, nd1=p.numberOfDaughters(); id1 < nd1; ++id1 ) {
02277 
02278         const Candidate * d1 = p.daughter(id1);
02279 
02280         // b - quark
02281         if ( abs(d1->pdgId()) == 5) {
02282 
02283           math::XYZTLorentzVector momentum=d1->p4();
02284           p4par[nJet++] = momentum=d1->p4();
02285 
02286           cout << "Daughter1: id = " << d1->pdgId() 
02287                << " daughters = " << d1->numberOfDaughters() 
02288                << " mother 1   = " << (d1->mother())->pdgId() 
02289                << " Momentum " << momentum << " GeV/c"
02290                << endl;   
02291 
02292           
02293           if ( (d1->mother())->pdgId() == 6 ) {
02294             p4part[0] = momentum=d1->p4();
02295             cout << ">>> part0 = " << p4part[0] << endl;
02296           }
02297           if ( (d1->mother())->pdgId() == -6) {
02298             p4part[3] = momentum=d1->p4();
02299             cout << ">>> part3 = " << p4part[3] << endl;
02300           }
02301 
02302         }
02303 
02304         
02305         // W
02306         // Check for fully hadronic decay
02307 
02308         if ( abs(d1->pdgId()) == 24) {
02309 
02310           for( size_t id2=0, nd2=d1->numberOfDaughters(); id2 < nd2; ++id2 ) {
02311 
02312             const Candidate * d2 = d1->daughter(id2);
02313 
02314             if (abs(d2->pdgId()) < 9) {
02315 
02316               math::XYZVector vertex(d2->vx(),d2->vy(),d2->vz());
02317               math::XYZTLorentzVector momentum=d2->p4();
02318               p4par[nJet++] = momentum=d2->p4();
02319 
02320               if ( (d1->mother())->pdgId() == 6 ) {
02321                 p4part[ii] = momentum=d2->p4();
02322                 cout << ">>> part" << ii << " = " << p4part[ii] << endl;
02323                 ii++;
02324               }
02325               if ( (d1->mother())->pdgId() == -6 ) {
02326                 p4part[jj] = momentum=d2->p4();
02327                 cout << ">>> part" << jj << " = " << p4part[jj] << endl;
02328                 jj++;
02329               }
02330 
02331               cout << "Daughter2: id = " << d2->pdgId() 
02332                    << " daughters = " << d2->numberOfDaughters() 
02333                    << " mother 2   = " << (d2->mother())->pdgId() 
02334                    << " Momentum " << momentum << " GeV/c"
02335                    << endl;
02336             }
02337 
02338           }
02339         }
02340 
02341         //      if ( pdgId == d->pdgId() && d->status() == 1 ) {        
02342         //      }
02343       }
02344     }
02345 
02346   }
02347   //  cout << ">>> N Jets = " << nJet << endl;
02348     
02349   if (nJet == 6) {
02350 
02351     double tmass    = (p4part[0]+p4part[1]+p4part[2]).mass();
02352     double tbarmass = (p4part[3]+p4part[4]+p4part[5]).mass();
02353 
02354     tMass.Fill(tmass);
02355     tbarMass.Fill(tbarmass);
02356     
02357     cout << ">>> T Mass = " << tmass << " / " << tbarmass << endl;
02358 
02359     double mindR = 1000.;
02360     for( size_t i=0; i<6; ++i ) { 
02361       for( size_t j=0; j<6; ++j ) { 
02362         if (j > i) {
02363           double delR = deltaR( p4par[i].eta(), p4par[i].phi(), p4par[j].eta(), p4par[j].phi() ); 
02364           if (delR < mindR) mindR = delR;
02365           dRParton.Fill(delR);
02366         }
02367       }
02368     }
02369     dRPartonMin.Fill(mindR);
02370 
02371     int usedInd;
02372     usedInd = -1;
02373     jetInd  = 0;
02374     
02375     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02376     for( int i=0; i<6; ++i ) { 
02377       
02378       dRmin[jetInd]  = 1000.0;
02379 
02380       int nj = 0;
02381       for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
02382            cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
02383         
02384         double delR = deltaR( cal1->eta(), cal1->phi(), p4par[i].eta(), p4par[i].phi() ); 
02385 
02386         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02387           dRmin[jetInd] = delR;        // delta R of match
02388           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
02389           usedInd       = nj;
02390           genInd        = i;
02391         }
02392 
02393         /****   
02394         cout << "Delta R = " << delR
02395              << " deltaR min = " << dRmin[jetInd]
02396              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02397              << " p4par = " << p4par[i]
02398              << " p4cal = " << cal1->p4()
02399              << " cal1 eta  = " << cal1->eta()
02400              << " p4par eta = " << p4par[i].eta()
02401              << endl;
02402         ****/
02403 
02404         
02405         nj++;
02406       }
02407 
02408 
02409       // Found matched jet
02410       if (usedInd != -1) {
02411         used[usedInd] = 1;
02412             
02413         dRPar1.Fill(dRmin[jetInd]);
02414         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02415         dPhiPar1.Fill(dphi);
02416         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02417         dEtaPar1.Fill(deta);
02418         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02419         dPtPar1.Fill(dpt);
02420         jetInd++;    
02421       }
02422 
02423     }
02424     ParMatch1.Fill(jetInd);
02425     if (jetInd == 6) {
02426       topMass1.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02427       topMass1.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02428     }
02429 
02430     /***
02431     cout << "Collection Size = " <<  calo1Jets->size() 
02432          << " / " << jetInd
02433          << endl;
02434     ***/
02435 
02436     // ***********************
02437     jetInd = 0;
02438     usedInd = -1;
02439 
02440     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02441     for( int i=0; i<6; ++i ) { 
02442       
02443       dRmin[jetInd]  = 1000.0;
02444 
02445       int nj = 0;      
02446       for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); 
02447            cal2 != calo2Jets->end() && jetInd<maxJets; ++cal2 ) { 
02448         
02449         double delR = deltaR( cal2->eta(), cal2->phi(), p4par[i].eta(), p4par[i].phi() ); 
02450 
02451         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02452           dRmin[jetInd] = delR;        // delta R of match
02453           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
02454           usedInd       = nj;
02455           genInd        = i;
02456         }
02457 
02458         /****
02459         cout << "Delta R = " << delR
02460              << " deltaR min = " << dRmin[jetInd]
02461              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02462              << " cal2 eta  = " << cal2->eta()
02463              << " p4par eta = " << p4par[i].eta()
02464              << endl;
02465         ****/
02466         
02467         nj++;   
02468       }
02469       if (usedInd != -1) {
02470         used[usedInd] = 1;
02471 
02472         dRPar2.Fill(dRmin[jetInd]);
02473         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02474         dPhiPar2.Fill(dphi);
02475         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02476         dEtaPar2.Fill(deta);
02477         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02478         dPtPar2.Fill(dpt);
02479         
02480         jetInd++;    
02481       }
02482 
02483     }
02484     ParMatch2.Fill(jetInd);
02485     if (jetInd == 6) {
02486       topMass2.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02487       topMass2.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02488     }
02489 
02490 
02491     /***
02492     cout << "Collection Size = " <<  calo2Jets->size() 
02493          << " / " << jetInd
02494          << endl;
02495     ***/
02496 
02497 
02498     // ***********************
02499     jetInd = 0;
02500     usedInd = -1;
02501 
02502     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02503     for( int i=0; i<6; ++i ) { 
02504       
02505       dRmin[jetInd]  = 1000.0;
02506 
02507       int nj = 0;
02508       for( CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); 
02509            cal3 != calo3Jets->end() && jetInd<maxJets; ++cal3 ) { 
02510         
02511         double delR = deltaR( cal3->eta(), cal3->phi(), p4par[i].eta(), p4par[i].phi() ); 
02512 
02513         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) )  {
02514           dRmin[jetInd] = delR;        // delta R of match
02515           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
02516           usedInd       = nj;
02517           genInd        = i;
02518         }
02519         /****
02520         cout << "Delta R = " << delR
02521              << " deltaR min = " << dRmin[jetInd]
02522              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02523              << " cal3 eta  = " << cal3->eta()
02524              << " p4par eta = " << p4par[i].eta()
02525              << endl;
02526         ****/
02527         
02528         nj++;
02529       }
02530       if (usedInd != -1) {
02531         used[usedInd] = 1;
02532 
02533         dRPar3.Fill(dRmin[jetInd]);
02534         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02535         dPhiPar3.Fill(dphi);
02536         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02537         dEtaPar3.Fill(deta);
02538         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02539         dPtPar3.Fill(dpt);
02540         
02541         jetInd++;    
02542       }
02543     }
02544     ParMatch3.Fill(jetInd);
02545     if (jetInd == 6) {
02546       topMass3.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02547       topMass3.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02548     }
02549 
02550     /***
02551     cout << "Collection Size = " <<  calo3Jets->size() 
02552          << " / " << jetInd
02553          << endl;
02554     ***/
02555 
02556   }
02557 
02558   int nTow1, nTow2, nTow3, nTow4;
02559 
02560   Handle<CaloJetCollection> jets;
02561 
02562   // *********************
02563   // Jet Properties
02564   // *********************
02565 
02566 
02567 
02568   // --- Loop over jets and make a list of all the used towers
02569   evt.getByLabel( CaloJetAlgorithm1, jets );
02570   int jjet = 0;
02571   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02572     jjet++;
02573 
02574     float hadEne  = ijet->hadEnergyInHB() + ijet->hadEnergyInHO() + 
02575                     ijet->hadEnergyInHE() + ijet->hadEnergyInHF();                   
02576     float emEne   = ijet->emEnergyInEB() + ijet->emEnergyInEE() + ijet->emEnergyInHF();
02577     float had     = ijet->energyFractionHadronic();    
02578 
02579     float j_et = ijet->et();
02580 
02581     if (fabs(ijet->eta()) < 1.3) {
02582       totEneLeadJetEta1_1.Fill(hadEne+emEne); 
02583       hadEneLeadJetEta1_1.Fill(hadEne); 
02584       emEneLeadJetEta1_1.Fill(emEne);       
02585 
02586       totEneLeadJetEta1_2.Fill(hadEne+emEne); 
02587       hadEneLeadJetEta1_2.Fill(hadEne); 
02588       emEneLeadJetEta1_2.Fill(emEne);       
02589 
02590       if (ijet->pt() > minJetPt10) 
02591         hadFracEta11.Fill(had);
02592     }
02593     if ((fabs(ijet->eta()) > 1.3) && (fabs(ijet->eta()) < 3.) ) {
02594 
02595       totEneLeadJetEta2_1.Fill(hadEne+emEne); 
02596       hadEneLeadJetEta2_1.Fill(hadEne); 
02597       emEneLeadJetEta2_1.Fill(emEne);   
02598     
02599       totEneLeadJetEta2_2.Fill(hadEne+emEne); 
02600       hadEneLeadJetEta2_2.Fill(hadEne); 
02601       emEneLeadJetEta2_2.Fill(emEne);       
02602 
02603       if (ijet->pt() > minJetPt10) 
02604         hadFracEta21.Fill(had);
02605     }
02606     if (fabs(ijet->eta()) > 3.) {
02607 
02608       totEneLeadJetEta3_1.Fill(hadEne+emEne); 
02609       hadEneLeadJetEta3_1.Fill(hadEne); 
02610       emEneLeadJetEta3_1.Fill(emEne); 
02611 
02612       totEneLeadJetEta3_2.Fill(hadEne+emEne); 
02613       hadEneLeadJetEta3_2.Fill(hadEne); 
02614       emEneLeadJetEta3_2.Fill(emEne); 
02615 
02616       if (ijet->pt() > minJetPt10) 
02617         hadFracEta31.Fill(had);
02618     }
02619 
02620     
02621     if (jjet == 1) {
02622       hadFracLeadJet1.Fill(had);
02623       hadEneLeadJet1.Fill(hadEne);
02624       hadEneLeadJet12.Fill(hadEne);
02625       hadEneLeadJet13.Fill(hadEne);
02626       emEneLeadJet1.Fill(emEne);
02627       emEneLeadJet12.Fill(emEne);
02628       emEneLeadJet13.Fill(emEne);
02629     }
02630 
02631     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02632     int nConstituents = jetCaloRefs.size();
02633 
02634     if (jjet == 1) {
02635 
02636       nTow1 = nTow2 = nTow3 = nTow4 = 0;
02637       for (int i = 0; i <nConstituents ; i++){
02638 
02639         float et  = jetCaloRefs[i]->et();
02640 
02641         if (et > 0.5) nTow1++;
02642         if (et > 1.0) nTow2++;
02643         if (et > 1.5) nTow3++;
02644         if (et > 2.0) nTow4++;
02645         
02646         hf_TowerJetEt1.Fill(et/j_et);
02647 
02648       }
02649 
02650       nTowersLeadJetPt1.Fill(nTow1);
02651       nTowersLeadJetPt2.Fill(nTow2);
02652       nTowersLeadJetPt3.Fill(nTow3);
02653       nTowersLeadJetPt4.Fill(nTow4);
02654 
02655     }
02656 
02657     if ( (jjet == 1) && (fabs(ijet->eta()) < 1.3) ) {
02658       nTowersLeadJet1.Fill( nConstituents );    
02659       
02660       for (int i = 0; i <nConstituents ; i++){
02661         float t_et  = jetCaloRefs[i]->et();
02662         double delR = deltaR( ijet->eta(), ijet->phi(), jetCaloRefs[i]->eta(), jetCaloRefs[i]->phi() );
02663         hf_TowerDelR1.Fill( delR,  t_et/j_et);
02664         hf_TowerDelR12.Fill( delR,  t_et/j_et);
02665         TowerEtLeadJet1.Fill( t_et );
02666         TowerEtLeadJet12.Fill( t_et );
02667         TowerEtLeadJet13.Fill( t_et );
02668       }
02669     }
02670 
02671 
02672     if ( (jjet == 2) && (fabs(ijet->eta()) < 1.3) ) {
02673       nTowersSecondJet1.Fill( nConstituents );    
02674     }
02675 
02676 
02677 
02678   }
02679 
02680   // *********************
02681   // Unclustered Energy
02682   // *********************
02683 
02684   double SumPtJet(0);
02685 
02686   double SumEtNotJets(0);
02687   double SumEtJets(0);
02688   double SumEtTowers(0);
02689 
02690   double sumJetPx(0);
02691   double sumJetPy(0);  
02692 
02693   double sumTowerAllPx(0);
02694   double sumTowerAllPy(0);
02695 
02696   double sumTowerAllEx(0);
02697   double sumTowerAllEy(0);
02698 
02699   std::vector<CaloTowerPtr>   UsedTowerList;
02700   std::vector<CaloTower>      TowerUsedInJets;
02701   std::vector<CaloTower>      TowerNotUsedInJets;
02702 
02703 
02704   // *********************
02705   // Towers
02706 
02707   Handle<CaloTowerCollection> caloTowers;
02708   evt.getByLabel( "towerMaker", caloTowers );
02709 
02710 
02711 
02712   nTow1 = nTow2 = nTow3 = nTow4 = 0;
02713 
02714   double sum_et = 0.0;
02715   double sum_ex = 0.0;
02716   double sum_ey = 0.0;
02717   //  double sum_ez = 0.0;
02718 
02719   // --- Loop over towers and make a lists of used and unused towers
02720   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02721        tower != caloTowers->end(); tower++) {
02722     
02723     Double_t  et = tower->et();
02724     
02725     if (et > 0.5) nTow1++;
02726     if (et > 1.0) nTow2++;
02727     if (et > 1.5) nTow3++;
02728     if (et > 2.0) nTow4++;
02729 
02730     if(et>0.5) {
02731 
02732       // ********
02733       double phix   = tower->phi();
02734       //      double theta = tower->theta();
02735       //      double e     = tower->energy();
02736       //      double et    = e*sin(theta);
02737       //      double et    = tower->emEt() + tower->hadEt();
02738       double et    = tower->et();
02739 
02740       //      sum_ez += e*cos(theta);
02741       sum_et += et;
02742       sum_ex += et*cos(phix);
02743       sum_ey += et*sin(phix);
02744       // ********
02745 
02746       Double_t phi = tower->phi();
02747       SumEtTowers += tower->et();
02748       
02749       sumTowerAllEx += et*cos(phi);
02750       sumTowerAllEy += et*sin(phi);
02751 
02752     }
02753         
02754   }
02755 
02756   SumEt1.Fill(sum_et);
02757   SumEt12.Fill(sum_et);
02758   SumEt13.Fill(sum_et);
02759 
02760   MET1.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02761   MET12.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02762   MET13.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02763 
02764   //  met->mex   = -sum_ex;
02765   //  met->mey   = -sum_ey;
02766   //  met->mez   = -sum_ez;
02767   //  met->met   = sqrt( sum_ex*sum_ex + sum_ey*sum_ey );
02768   // cout << "MET = " << met->met << endl;
02769   //  met->sumet = sum_et;
02770   //  met->phi   = atan2( -sum_ey, -sum_ex ); 
02771 
02772   
02773 
02774   hf_sumTowerAllEx.Fill(sumTowerAllEx);
02775   hf_sumTowerAllEy.Fill(sumTowerAllEy);
02776 
02777   nTowers1.Fill(nTow1);
02778   nTowers2.Fill(nTow2);
02779   nTowers3.Fill(nTow3);
02780   nTowers4.Fill(nTow4);
02781 
02782   // *********************
02783   // MidPoint 
02784   //
02785   UsedTowerList.clear();
02786   TowerUsedInJets.clear();
02787   TowerNotUsedInJets.clear();
02788 
02789   // --- Loop over jets and make a list of all the used towers
02790   evt.getByLabel( CaloJetAlgorithm1, jets );
02791   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02792 
02793     Double_t jetPt  = ijet->pt();
02794     Double_t jetPhi = ijet->phi();
02795     
02796     //    if (jetPt>5.0) {
02797 
02798       Double_t jetPx = jetPt*cos(jetPhi);
02799       Double_t jetPy = jetPt*sin(jetPhi);
02800       
02801       sumJetPx +=jetPx;
02802       sumJetPy +=jetPy;
02803 
02804       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02805       int nConstituents = jetCaloRefs.size();
02806       for (int i = 0; i <nConstituents ; i++){
02807         UsedTowerList.push_back(jetCaloRefs[i]);
02808       }
02809 
02810       SumPtJet +=jetPt;
02811     //    }
02812 
02813   }
02814 
02815 
02816   int NTowersUsed = UsedTowerList.size();
02817       
02818   // --- Loop over towers and make a lists of used and unused towers
02819   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02820        tower != caloTowers->end(); tower++) {
02821     
02822     CaloTower  t = *tower;
02823     Double_t  et = tower->et();
02824 
02825     if(et>0) {
02826 
02827       Double_t phi = tower->phi();
02828       SumEtTowers += tower->et();
02829       
02830       sumTowerAllPx += et*cos(phi);
02831       sumTowerAllPy += et*sin(phi);
02832 
02833       bool used = false;
02834 
02835       for(int i=0; i<NTowersUsed; i++){
02836         if(tower->id() == UsedTowerList[i]->id()){
02837           used=true;
02838           break;
02839         }
02840       }
02841       
02842       if (used) {
02843         TowerUsedInJets.push_back(t);
02844       } else {
02845         TowerNotUsedInJets.push_back(t);
02846       }
02847 
02848     }
02849 
02850   }
02851 
02852   int nUsed    = TowerUsedInJets.size();
02853   int nNotUsed = TowerNotUsedInJets.size();
02854   
02855   SumEtJets    = 0;
02856   SumEtNotJets = 0;
02857 
02858   for(int i=0;i<nUsed;i++){
02859     SumEtJets += TowerUsedInJets[i].et();
02860   }
02861   h_jetEt1.Fill(SumEtJets);
02862 
02863   for(int i=0;i<nNotUsed;i++){
02864     if (TowerNotUsedInJets[i].et() > 0.5) 
02865       SumEtNotJets += TowerNotUsedInJets[i].et();
02866     h_missEt1.Fill(TowerNotUsedInJets[i].et());
02867     h_missEt1s.Fill(TowerNotUsedInJets[i].et());
02868   }
02869   h_totMissEt1.Fill(SumEtNotJets);
02870 
02871 
02872 
02873   // *********************
02874   // SISCone
02875   //
02876   UsedTowerList.clear();
02877   TowerUsedInJets.clear();
02878   TowerNotUsedInJets.clear();
02879 
02880   // --- Loop over jets and make a list of all the used towers
02881   evt.getByLabel( CaloJetAlgorithm2, jets );
02882   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02883 
02884     Double_t jetPt  = ijet->pt();
02885     Double_t jetPhi = ijet->phi();
02886     
02887     //    if (jetPt>5.0) {
02888 
02889       Double_t jetPx = jetPt*cos(jetPhi);
02890       Double_t jetPy = jetPt*sin(jetPhi);
02891       
02892       sumJetPx +=jetPx;
02893       sumJetPy +=jetPy;
02894 
02895       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02896       int nConstituents = jetCaloRefs.size();
02897       for (int i = 0; i <nConstituents ; i++){
02898         UsedTowerList.push_back(jetCaloRefs[i]);
02899       }
02900 
02901       SumPtJet +=jetPt;
02902     //    }
02903 
02904   }
02905 
02906 
02907   //  Handle<CaloTowerCollection> caloTowers;
02908   //  evt.getByLabel( "towerMaker", caloTowers );
02909 
02910   NTowersUsed = UsedTowerList.size();
02911       
02912   // --- Loop over towers and make a lists of used and unused towers
02913   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02914        tower != caloTowers->end(); tower++) {
02915     
02916     CaloTower  t = *tower;
02917     Double_t  et = tower->et();
02918 
02919     if(et>0) {
02920 
02921       Double_t phi = tower->phi();
02922 
02923       SumEtTowers += tower->et();
02924       
02925       sumTowerAllPx += et*cos(phi);
02926       sumTowerAllPy += et*sin(phi);
02927 
02928       bool used = false;
02929 
02930       for(int i=0; i<NTowersUsed; i++){
02931         if(tower->id() == UsedTowerList[i]->id()){
02932           used=true;
02933           break;
02934         }
02935       }
02936       
02937       if (used) {
02938         TowerUsedInJets.push_back(t);
02939       } else {
02940         TowerNotUsedInJets.push_back(t);
02941       }
02942 
02943     }
02944 
02945   }
02946 
02947   nUsed    = TowerUsedInJets.size();
02948   nNotUsed = TowerNotUsedInJets.size();
02949   
02950   SumEtJets    = 0;
02951   SumEtNotJets = 0;
02952 
02953   for(int i=0;i<nUsed;i++){
02954     SumEtJets += TowerUsedInJets[i].et();
02955   }
02956   h_jetEt2.Fill(SumEtJets);
02957 
02958   for(int i=0;i<nNotUsed;i++){
02959     if (TowerNotUsedInJets[i].et() > 0.5) 
02960       SumEtNotJets += TowerNotUsedInJets[i].et();
02961     h_missEt2.Fill(TowerNotUsedInJets[i].et());
02962     h_missEt2s.Fill(TowerNotUsedInJets[i].et());
02963   }
02964   h_totMissEt2.Fill(SumEtNotJets);
02965 
02966 
02967   // *********************
02968   // KtClus
02969   //
02970   UsedTowerList.clear();
02971   TowerUsedInJets.clear();
02972   TowerNotUsedInJets.clear();
02973 
02974   // --- Loop over jets and make a list of all the used towers
02975   evt.getByLabel( CaloJetAlgorithm3, jets );
02976   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02977 
02978     Double_t jetPt  = ijet->pt();
02979     Double_t jetPhi = ijet->phi();
02980     
02981     //    if (jetPt>5.0) {
02982 
02983       Double_t jetPx = jetPt*cos(jetPhi);
02984       Double_t jetPy = jetPt*sin(jetPhi);
02985       
02986       sumJetPx +=jetPx;
02987       sumJetPy +=jetPy;
02988 
02989       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02990       int nConstituents = jetCaloRefs.size();
02991       for (int i = 0; i <nConstituents ; i++){
02992         UsedTowerList.push_back(jetCaloRefs[i]);
02993       }
02994 
02995       SumPtJet +=jetPt;
02996     //    }
02997 
02998   }
02999 
03000 
03001   //  Handle<CaloTowerCollection> caloTowers;
03002   //  evt.getByLabel( "towerMaker", caloTowers );
03003 
03004   NTowersUsed = UsedTowerList.size();
03005       
03006   // --- Loop over towers and make a lists of used and unused towers
03007   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
03008        tower != caloTowers->end(); tower++) {
03009     
03010     CaloTower  t = *tower;
03011     Double_t  et = tower->et();
03012 
03013     if(et>0) {
03014 
03015       //      Double_t phi = tower->phi();
03016 
03017       //      SumEtTowers   += tower->et();      
03018       //      sumTowerAllPx += et*cos(phi);
03019       //      sumTowerAllPy += et*sin(phi);
03020 
03021       bool used = false;
03022 
03023       for(int i=0; i<NTowersUsed; i++){
03024         if(tower->id() == UsedTowerList[i]->id()){
03025           used=true;
03026           break;
03027         }
03028       }
03029       
03030       if (used) {
03031         TowerUsedInJets.push_back(t);
03032       } else {
03033         TowerNotUsedInJets.push_back(t);
03034       }
03035 
03036     }
03037 
03038   }
03039 
03040   nUsed    = TowerUsedInJets.size();
03041   nNotUsed = TowerNotUsedInJets.size();
03042   
03043   SumEtJets    = 0;
03044   SumEtNotJets = 0;
03045 
03046   for(int i=0;i<nUsed;i++){
03047     SumEtJets += TowerUsedInJets[i].et();
03048   }
03049   h_jetEt3.Fill(SumEtJets);
03050 
03051   for(int i=0;i<nNotUsed;i++){
03052     if (TowerNotUsedInJets[i].et() > 0.5) 
03053       SumEtNotJets += TowerNotUsedInJets[i].et();
03054     h_missEt3.Fill(TowerNotUsedInJets[i].et());
03055     h_missEt3s.Fill(TowerNotUsedInJets[i].et());
03056   }
03057   h_totMissEt3.Fill(SumEtNotJets);
03058 
03059 }
03060 
03061 
03062 
03063 void myFastSimVal::endJob() {
03064 
03065   //Write out the histogram file.
03066   m_file->Write(); 
03067 
03068 }
03069 #include "FWCore/Framework/interface/MakerMacros.h"
03070 DEFINE_FWK_MODULE(myFastSimVal);