CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/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 
00663   float minJetPt = 30.;
00664   float minJetPt10 = 10.;
00665   int jetInd, allJetInd;
00666   int usedInd = -1;  
00667   //  double matchedDelR = 0.1;
00668     double matchedDelR = 0.3;
00669 
00670   ZpMG = 0;
00671   LeadMass1 = -1;
00672   LeadMass2 = -1;
00673   LeadMass3 = -1;
00674 
00675   math::XYZTLorentzVector p4tmp[2], p4cortmp[2];
00676   nEvent++;
00677 
00678   // ********************************
00679   // **** Get the CaloJet1 collection
00680   // ********************************
00681 
00682   Handle<CaloJetCollection> caloJets1;
00683   evt.getByLabel( CaloJetAlgorithm1, caloJets1 );
00684 
00685   // Count Jets above Pt cut
00686   for (int istep = 0; istep < 100; ++istep) {
00687     int     njet = 0;
00688     float ptStep = (istep * (5000./100.));
00689     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00690       if ( cal->pt() > ptStep ) njet++;      
00691     }
00692 
00693     hf_nJet1.Fill( ptStep, njet );
00694   }
00695 
00696   // Count Jets above Pt cut
00697   for (int istep = 0; istep < 100; ++istep) {
00698     int     njet = 0;
00699     float ptStep = (istep * (200./100.));
00700     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00701       if ( cal->pt() > ptStep ) njet++;      
00702     }
00703 
00704     hf_nJet1s.Fill( ptStep, njet );
00705   }
00706 
00707   // Count Jets above Pt cut
00708   for (int istep = 0; istep < 100; ++istep) {
00709     int     njet = 0;
00710     float ptStep = (istep * (3000./100.));
00711     for ( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {          
00712       if ( cal->pt() > ptStep ) njet++;      
00713     }
00714 
00715     hf_nJet11.Fill( ptStep, njet );
00716   }
00717 
00718 
00719   //Loop over the two leading CaloJets and fill some histograms
00720   jetInd    = 0;
00721   allJetInd = 0;
00722   EtaOk10 = 0;
00723   EtaOk13 = 0;
00724   EtaOk40 = 0;
00725 
00726   //  const JetCorrector* corrector = 
00727   //    JetCorrector::getJetCorrector (JetCorrectionService, es);
00728 
00729   double highestPt;
00730   double nextPt;
00731 
00732   highestPt = 0.0;
00733   nextPt    = 0.0;
00734   
00735 
00736   for( CaloJetCollection::const_iterator cal = caloJets1->begin(); cal != caloJets1->end(); ++ cal ) {
00737     
00738     //    double scale = corrector->correction (*cal);
00739     double scale = 1.0;
00740     double corPt = scale*cal->pt();
00741     //    double corPt = cal->pt();
00742 
00743     
00744     if (corPt>highestPt) {
00745       nextPt      = highestPt;
00746       p4cortmp[1] = p4cortmp[0]; 
00747       highestPt   = corPt;
00748       p4cortmp[0] = scale*cal->p4();
00749     } else if (corPt>nextPt) {
00750       nextPt      = corPt;
00751       p4cortmp[1] = scale*cal->p4();
00752     }
00753 
00754     /***
00755     std::cout << ">>> Corr Jet: corPt = " 
00756               << corPt << ", scale = " << scale 
00757               << " pt = " << cal->pt()
00758               << " highestPt = " << highestPt 
00759               << " nextPt = "    << nextPt 
00760               << std::endl;  
00761     ****/
00762 
00763     allJetInd++;
00764     if (allJetInd == 1) {
00765       h_jet1Pt1.Fill( cal->pt() );
00766       p4tmp[0] = cal->p4();
00767       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
00768       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
00769       if ( fabs(cal->eta()) < 4.0) EtaOk40++;
00770     }
00771     if (allJetInd == 2) {
00772       h_jet2Pt1.Fill( cal->pt() );
00773       p4tmp[1] = cal->p4();
00774       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
00775       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
00776       if ( fabs(cal->eta()) < 4.0) EtaOk40++;
00777     }
00778     if ( (allJetInd == 1) || (allJetInd == 2) ) {
00779 
00780       h_ptCalL1.Fill( cal->pt() );   
00781       h_ptCalL12.Fill( cal->pt() );   
00782       h_ptCalL13.Fill( cal->pt() );   
00783 
00784       h_etaCalL1.Fill( cal->eta() );
00785       h_phiCalL1.Fill( cal->phi() );
00786     }
00787 
00788     if (allJetInd == 3) h_jet3Pt1.Fill( cal->pt() );
00789     if (allJetInd == 4) h_jet4Pt1.Fill( cal->pt() );
00790     if (allJetInd == 5) h_jet5Pt1.Fill( cal->pt() );
00791     if (allJetInd == 6) h_jet6Pt1.Fill( cal->pt() );
00792     if (allJetInd == 7) h_jet7Pt1.Fill( cal->pt() );
00793 
00794     if ( fabs(cal->eta()) < 1.3) {
00795       h_lowPtCal11.Fill( cal->pt() );   
00796       if ( cal->pt() > 10.) {
00797         h_lowPtCal1c11.Fill( cal->pt() );   
00798       }
00799     }
00800     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
00801       h_lowPtCal12.Fill( cal->pt() );   
00802       if ( cal->pt() > 10.) {
00803         h_lowPtCal1c12.Fill( cal->pt() );   
00804       }
00805     }
00806     if ( fabs(cal->eta()) > 3.0) {
00807       h_lowPtCal13.Fill( cal->pt() );   
00808       if ( cal->pt() > 10.) {
00809         h_lowPtCal1c13.Fill( cal->pt() );   
00810       }
00811     }
00812 
00813 
00814 
00815     if ( cal->pt() > minJetPt) {
00816       //    std::cout << "CALO JET1 #" << jetInd << std::endl << cal->print() << std::endl;
00817       h_ptCal1.Fill( cal->pt() );   
00818       h_ptCal12.Fill( cal->pt() );   
00819       h_ptCal13.Fill( cal->pt() );   
00820 
00821       h_etaCal1.Fill( cal->eta() );
00822       h_phiCal1.Fill( cal->phi() );
00823       jetInd++;
00824     }
00825   }
00826 
00827   //  h_nCalJets1.Fill( caloJets1->size() );   
00828   h_nCalJets1.Fill( jetInd ); 
00829   if (jetInd > 1) {
00830     LeadMass1 = (p4tmp[0]+p4tmp[1]).mass();
00831     dijetMass1.Fill( LeadMass1 );    
00832     dijetMass12.Fill( LeadMass1 );    
00833     dijetMass13.Fill( LeadMass1 );    
00834     if (EtaOk10 == 2) {
00835       dijetMass101.Fill( LeadMass1 );        
00836       dijetMass102.Fill( LeadMass1 );  
00837       dijetMass103.Fill( LeadMass1 );  
00838       dijetMass_700_101.Fill( LeadMass1 );  
00839       dijetMass_2000_101.Fill( LeadMass1 );  
00840       dijetMass_5000_101.Fill( LeadMass1 );  
00841     }
00842     if (EtaOk13 == 2) {
00843       dijetMass131.Fill( LeadMass1 );  
00844       dijetMass132.Fill( LeadMass1 );  
00845       dijetMass133.Fill( LeadMass1 );  
00846       dijetMass_700_131.Fill( LeadMass1 );  
00847       dijetMass_2000_131.Fill( LeadMass1 );  
00848       dijetMass_5000_131.Fill( LeadMass1 );  
00849     }
00850     if (EtaOk40 == 2) {
00851       dijetMass401.Fill( LeadMass1 );        
00852       dijetMass402.Fill( LeadMass1 );  
00853       dijetMass403.Fill( LeadMass1 );  
00854       dijetMass_700_401.Fill( LeadMass1 );  
00855       dijetMass_2000_401.Fill( LeadMass1 );  
00856       dijetMass_5000_401.Fill( LeadMass1 );  
00857     }
00858 
00859     LeadMass1 = (p4cortmp[0]+p4cortmp[1]).mass();
00860 
00861     /****
00862     if (LeadMass1 < 30.) {
00863       std::cout << " XXX Low Mass " 
00864                 << (p4tmp[0]+p4tmp[1]).mass() 
00865                 << " / " 
00866                 << (p4cortmp[0]+p4cortmp[1]).mass() 
00867                 << std::endl;
00868 
00869       std::cout << " p4 1 = " << p4tmp[0]
00870                 << " p4 2 = " << p4tmp[1]
00871                 << " p4 cor 1 = " << p4cortmp[0]
00872                 << " p4 cor 2 = " << p4cortmp[0]
00873                 << endl;
00874 
00875     }
00876     ****/
00877 
00878     /****
00879     dijetMassCor1.Fill( LeadMass1 );    
00880     dijetMassCor_700_1.Fill( LeadMass1 );    
00881     dijetMassCor_2000_1.Fill( LeadMass1 );    
00882     dijetMassCor_5000_1.Fill( LeadMass1 );    
00883 
00884     if (EtaOk10 == 2) {
00885       dijetMassCor101.Fill( LeadMass1 );  
00886       dijetMassCor_700_101.Fill( LeadMass1 );  
00887       dijetMassCor_2000_101.Fill( LeadMass1 );  
00888       dijetMassCor_5000_101.Fill( LeadMass1 );  
00889     }
00890     if (EtaOk13 == 2) {
00891       dijetMassCor131.Fill( LeadMass1 );  
00892       dijetMassCor_700_131.Fill( LeadMass1 );  
00893       dijetMassCor_2000_131.Fill( LeadMass1 );  
00894       dijetMassCor_5000_131.Fill( LeadMass1 );  
00895     }
00896     if (EtaOk40 == 2) {
00897       dijetMassCor401.Fill( LeadMass1 ); 
00898       dijetMassCor_700_401.Fill( LeadMass1 ); 
00899       dijetMassCor_2000_401.Fill( LeadMass1 ); 
00900       dijetMassCor_5000_401.Fill( LeadMass1 ); 
00901     }
00902     ****/
00903 
00904   }
00905 
00906 
00907   // ********************************
00908   // **** Get the CaloJet2 collection
00909   // ********************************
00910   Handle<CaloJetCollection> caloJets2;
00911   evt.getByLabel( CaloJetAlgorithm2, caloJets2 );
00912 
00913   // Count Jets above Pt cut
00914   for (int istep = 0; istep < 100; ++istep) {
00915     int     njet = 0;
00916     float ptStep = (istep * (5000./100.));
00917 
00918     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00919       if ( cal->pt() > ptStep ) njet++;      
00920     
00921     hf_nJet2.Fill( ptStep, njet );
00922   }
00923 
00924   for (int istep = 0; istep < 100; ++istep) {
00925     int     njet = 0;
00926     float ptStep = (istep * (200./100.));
00927 
00928     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00929       if ( cal->pt() > ptStep ) njet++;      
00930     
00931     hf_nJet2s.Fill( ptStep, njet );
00932   }
00933 
00934 
00935   for (int istep = 0; istep < 100; ++istep) {
00936     int     njet = 0;
00937     float ptStep = (istep * (3000./100.));
00938 
00939     for ( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++ cal )
00940       if ( cal->pt() > ptStep ) njet++;      
00941     
00942     hf_nJet21.Fill( ptStep, njet );
00943   }
00944 
00945 
00946 
00947 
00948 
00949   //Loop over the two leading CaloJets and fill some histograms
00950   jetInd = 0;
00951   allJetInd = 0;
00952   for( CaloJetCollection::const_iterator cal = caloJets2->begin(); cal != caloJets2->end(); ++cal ) {
00953 
00954     allJetInd++;
00955     if (allJetInd == 1) {
00956       h_jet1Pt2.Fill( cal->pt() );
00957       p4tmp[0] = cal->p4();
00958     }
00959     if (allJetInd == 2) {
00960       h_jet2Pt2.Fill( cal->pt() );
00961       p4tmp[1] = cal->p4();
00962     }
00963     if ( (allJetInd == 1) || (allJetInd == 2) ) {
00964       h_ptCalL2.Fill( cal->pt() );   
00965       h_ptCalL22.Fill( cal->pt() );   
00966       h_ptCalL23.Fill( cal->pt() );   
00967 
00968       h_etaCalL2.Fill( cal->eta() );
00969       h_phiCalL2.Fill( cal->phi() );
00970     }
00971     if (allJetInd == 3) h_jet3Pt2.Fill( cal->pt() );
00972     if (allJetInd == 4) h_jet4Pt2.Fill( cal->pt() );
00973     if (allJetInd == 5) h_jet5Pt2.Fill( cal->pt() );
00974     if (allJetInd == 6) h_jet6Pt2.Fill( cal->pt() );
00975     if (allJetInd == 7) h_jet7Pt2.Fill( cal->pt() );
00976 
00977     if ( fabs(cal->eta()) < 1.3) {
00978       h_lowPtCal21.Fill( cal->pt() );   
00979       if ( cal->pt() > 10.) {
00980         h_lowPtCal2c11.Fill( cal->pt() );   
00981       }
00982     }
00983     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
00984       h_lowPtCal22.Fill( cal->pt() );   
00985       if ( cal->pt() > 10.) {
00986         h_lowPtCal2c12.Fill( cal->pt() );   
00987       }
00988     }
00989     if ( fabs(cal->eta()) > 3.0) {
00990       h_lowPtCal23.Fill( cal->pt() );   
00991       if ( cal->pt() > 10.) {
00992         h_lowPtCal2c13.Fill( cal->pt() );   
00993       }
00994     }
00995 
00996 
00997     if ( cal->pt() > minJetPt) {
00998       h_ptCal2.Fill( cal->pt() );   
00999       h_ptCal22.Fill( cal->pt() );   
01000       h_ptCal23.Fill( cal->pt() );   
01001 
01002       h_etaCal2.Fill( cal->eta() );
01003       h_phiCal2.Fill( cal->phi() );
01004       jetInd++;
01005     }
01006   }
01007   //  h_nCalJets2.Fill( caloJets2->size() ); 
01008   h_nCalJets2.Fill( jetInd ); 
01009   if (jetInd > 1) {
01010     LeadMass2 = (p4tmp[0]+p4tmp[1]).mass();
01011     dijetMass2.Fill( LeadMass2 );
01012     dijetMass22.Fill( LeadMass2 );
01013     dijetMass23.Fill( LeadMass2 );
01014 
01015 
01016     dijetMassCor1.Fill( LeadMass2 );    
01017     dijetMassCor_700_1.Fill( LeadMass2 );    
01018     dijetMassCor_2000_1.Fill( LeadMass2 );    
01019     dijetMassCor_5000_1.Fill( LeadMass2 );    
01020 
01021     if (EtaOk10 == 2) {
01022       dijetMassCor101.Fill( LeadMass2 );  
01023       dijetMassCor_700_101.Fill( LeadMass2 );  
01024       dijetMassCor_2000_101.Fill( LeadMass2 );  
01025       dijetMassCor_5000_101.Fill( LeadMass2 );  
01026     }
01027     if (EtaOk13 == 2) {
01028       dijetMassCor131.Fill( LeadMass2 );  
01029       dijetMassCor_700_131.Fill( LeadMass2 );  
01030       dijetMassCor_2000_131.Fill( LeadMass2 );  
01031       dijetMassCor_5000_131.Fill( LeadMass2 );  
01032     }
01033     if (EtaOk40 == 2) {
01034       dijetMassCor401.Fill( LeadMass2 ); 
01035       dijetMassCor_700_401.Fill( LeadMass2 ); 
01036       dijetMassCor_2000_401.Fill( LeadMass2 ); 
01037       dijetMassCor_5000_401.Fill( LeadMass2 ); 
01038     }
01039 
01040 
01041   }
01042 
01043 
01044 
01045   // ********************************
01046   // **** Get the CaloJet3 collection
01047   // ********************************
01048   Handle<CaloJetCollection> caloJets3;
01049   evt.getByLabel( CaloJetAlgorithm3, caloJets3 );
01050 
01051   //Loop over the two leading CaloJets and fill some histograms
01052   jetInd = 0;
01053   allJetInd = 0;
01054 
01055   // Count Jets above Pt cut
01056   for (int istep = 0; istep < 100; ++istep) {
01057     int     njet = 0;
01058     float ptStep = (istep * (5000./100.));
01059 
01060     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01061       if ( cal->pt() > ptStep ) njet++;      
01062     
01063 
01064     hf_nJet3.Fill( ptStep, njet );
01065   }
01066 
01067   for (int istep = 0; istep < 100; ++istep) {
01068     int     njet = 0;
01069     float ptStep = (istep * (200./100.));
01070 
01071     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01072       if ( cal->pt() > ptStep ) njet++;          
01073 
01074     hf_nJet3s.Fill( ptStep, njet );
01075   }
01076 
01077   for (int istep = 0; istep < 100; ++istep) {
01078     int     njet = 0;
01079     float ptStep = (istep * (3000./100.));
01080 
01081     for ( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal )
01082       if ( cal->pt() > ptStep ) njet++;          
01083 
01084     hf_nJet31.Fill( ptStep, njet );
01085   }
01086 
01087 
01088   for( CaloJetCollection::const_iterator cal = caloJets3->begin(); cal != caloJets3->end(); ++ cal ) {
01089 
01090     allJetInd++;
01091     if (allJetInd == 1) {
01092       h_jet1Pt3.Fill( cal->pt() );
01093       p4tmp[0] = cal->p4();
01094     }
01095     if (allJetInd == 2) {
01096       h_jet2Pt3.Fill( cal->pt() );
01097       p4tmp[1] = cal->p4();
01098     }
01099     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01100       h_ptCalL3.Fill( cal->pt() );   
01101       h_ptCalL32.Fill( cal->pt() );   
01102       h_ptCalL33.Fill( cal->pt() );   
01103 
01104       h_etaCalL3.Fill( cal->eta() );
01105       h_phiCalL3.Fill( cal->phi() );
01106     }
01107     if (allJetInd == 3) h_jet3Pt3.Fill( cal->pt() );
01108     if (allJetInd == 4) h_jet4Pt3.Fill( cal->pt() );
01109     if (allJetInd == 5) h_jet5Pt3.Fill( cal->pt() );
01110     if (allJetInd == 6) h_jet6Pt3.Fill( cal->pt() );
01111     if (allJetInd == 7) h_jet7Pt3.Fill( cal->pt() );
01112 
01113 
01114     if ( fabs(cal->eta()) < 1.3) {
01115       h_lowPtCal31.Fill( cal->pt() );   
01116       if ( cal->pt() > 10.) {
01117         h_lowPtCal3c11.Fill( cal->pt() );   
01118       }
01119     }
01120     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
01121       h_lowPtCal32.Fill( cal->pt() );   
01122       if ( cal->pt() > 10.) {
01123         h_lowPtCal3c12.Fill( cal->pt() );   
01124       }
01125     }
01126     if ( fabs(cal->eta()) > 3.0) {
01127       h_lowPtCal33.Fill( cal->pt() );   
01128       if ( cal->pt() > 10.) {
01129         h_lowPtCal3c13.Fill( cal->pt() );   
01130       }
01131     }
01132 
01133 
01134 
01135     if ( cal->pt() > minJetPt) {
01136       //    std::cout << "CALO JET3 #" << jetInd << std::endl << cal->print() << std::endl;
01137       h_ptCal3.Fill( cal->pt() );   
01138       h_ptCal32.Fill( cal->pt() );   
01139       h_ptCal33.Fill( cal->pt() );   
01140 
01141       h_etaCal3.Fill( cal->eta() );
01142       h_phiCal3.Fill( cal->phi() );
01143       jetInd++;
01144     }
01145   }
01146   //  h_nCalJets3.Fill( caloJets3->size() ); 
01147   h_nCalJets3.Fill( jetInd ); 
01148   if (jetInd > 1) {
01149     LeadMass3 = (p4tmp[0]+p4tmp[1]).mass();
01150     dijetMass3.Fill( LeadMass3 );
01151     dijetMass32.Fill( LeadMass3 );
01152     dijetMass33.Fill( LeadMass3 );
01153   }
01154 
01155 
01156 
01157   // ********************************
01158   // **** Get the CaloJet4 collection
01159   // ********************************
01160 
01161   Handle<CaloJetCollection> caloJets4;
01162   evt.getByLabel( CaloJetAlgorithm4, caloJets4 );
01163 
01164   //Loop over the two leading CaloJets and fill some histograms
01165   jetInd = 0;
01166   allJetInd = 0;
01167 
01168   // Count Jets above Pt cut
01169   for (int istep = 0; istep < 100; ++istep) {
01170     int     njet = 0;
01171     float ptStep = (istep * (5000./100.));
01172 
01173     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01174       if ( cal->pt() > ptStep ) njet++;      
01175     
01176 
01177     hf_nJet4.Fill( ptStep, njet );
01178   }
01179 
01180   for (int istep = 0; istep < 100; ++istep) {
01181     int     njet = 0;
01182     float ptStep = (istep * (200./100.));
01183 
01184     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01185       if ( cal->pt() > ptStep ) njet++;          
01186 
01187     hf_nJet4s.Fill( ptStep, njet );
01188   }
01189 
01190   for (int istep = 0; istep < 100; ++istep) {
01191     int     njet = 0;
01192     float ptStep = (istep * (3000./100.));
01193 
01194     for ( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal )
01195       if ( cal->pt() > ptStep ) njet++;          
01196 
01197     hf_nJet41.Fill( ptStep, njet );
01198   }
01199 
01200 
01201   for( CaloJetCollection::const_iterator cal = caloJets4->begin(); cal != caloJets4->end(); ++ cal ) {
01202 
01203     allJetInd++;
01204     if (allJetInd == 1) {
01205       h_jet1Pt4.Fill( cal->pt() );
01206       p4tmp[0] = cal->p4();
01207     }
01208     if (allJetInd == 2) {
01209       h_jet2Pt4.Fill( cal->pt() );
01210       p4tmp[1] = cal->p4();
01211     }
01212     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01213       h_ptCalL4.Fill(  cal->pt() );   
01214       h_ptCalL42.Fill( cal->pt() );   
01215       h_ptCalL43.Fill( cal->pt() );   
01216 
01217       h_etaCalL4.Fill( cal->eta() );
01218       h_phiCalL4.Fill( cal->phi() );
01219     }
01220     if (allJetInd == 3) h_jet3Pt4.Fill( cal->pt() );
01221     if (allJetInd == 4) h_jet4Pt4.Fill( cal->pt() );
01222     if (allJetInd == 5) h_jet5Pt4.Fill( cal->pt() );
01223     if (allJetInd == 6) h_jet6Pt4.Fill( cal->pt() );
01224     if (allJetInd == 7) h_jet7Pt4.Fill( cal->pt() );
01225 
01226 
01227     if ( fabs(cal->eta()) < 1.3) {
01228       h_lowPtCal41.Fill( cal->pt() );   
01229       if ( cal->pt() > 10.) {
01230         h_lowPtCal4c11.Fill( cal->pt() );   
01231       }
01232     }
01233     if ( (fabs(cal->eta())> 1.3) && ( fabs(cal->eta()) < 3.) )  {
01234       h_lowPtCal42.Fill( cal->pt() );   
01235       if ( cal->pt() > 10.) {
01236         h_lowPtCal4c12.Fill( cal->pt() );   
01237       }
01238     }
01239     if ( fabs(cal->eta()) > 3.0) {
01240       h_lowPtCal43.Fill( cal->pt() );   
01241       if ( cal->pt() > 10.) {
01242         h_lowPtCal4c13.Fill( cal->pt() );   
01243       }
01244     }
01245 
01246 
01247     if ( cal->pt() > minJetPt) {
01248       //    std::cout << "CALO JET4 #" << jetInd << std::endl << cal->print() << std::endl;
01249       h_ptCal4.Fill( cal->pt() );   
01250       h_ptCal42.Fill( cal->pt() );   
01251       h_ptCal43.Fill( cal->pt() );   
01252 
01253       h_etaCal4.Fill( cal->eta() );
01254       h_phiCal4.Fill( cal->phi() );
01255       jetInd++;
01256     }
01257   }
01258   //  h_nCalJets4.Fill( caloJets4->size() ); 
01259   h_nCalJets4.Fill( jetInd ); 
01260   if (jetInd > 1) {
01261     LeadMass4 = (p4tmp[0]+p4tmp[1]).mass();
01262     dijetMass4.Fill( LeadMass4 );
01263     dijetMass42.Fill( LeadMass4 );
01264     dijetMass43.Fill( LeadMass4 );
01265   }
01266 
01267 
01268 
01269   // *********************************************
01270   // *********************************************
01271 
01272 
01273   //**** Get the GenJet1 collection
01274   Handle<GenJetCollection> genJets1;
01275   evt.getByLabel( GenJetAlgorithm1, genJets1 );
01276 
01277 
01278   //Loop over the two leading GenJets and fill some histograms
01279   jetInd    = 0;
01280   allJetInd = 0;
01281   for( GenJetCollection::const_iterator gen = genJets1->begin(); gen != genJets1->end(); ++ gen ) {
01282     allJetInd++;
01283     if (allJetInd == 1) {
01284       p4tmp[0] = gen->p4();
01285     }
01286     if (allJetInd == 2) {
01287       p4tmp[1] = gen->p4();
01288     }
01289 
01290     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01291       h_ptGenL1.Fill( gen->pt() );   
01292       h_ptGenL12.Fill( gen->pt() );   
01293       h_ptGenL13.Fill( gen->pt() );   
01294 
01295       h_etaGenL1.Fill( gen->eta() );
01296       h_phiGenL1.Fill( gen->phi() );
01297     }
01298 
01299     if ( gen->pt() > minJetPt) {
01300       // std::cout << "GEN JET1 #" << jetInd << std::endl << gen->print() << std::endl;
01301       h_ptGen1.Fill( gen->pt() );   
01302       h_ptGen12.Fill( gen->pt() );   
01303       h_ptGen13.Fill( gen->pt() );   
01304 
01305       h_etaGen1.Fill( gen->eta() );
01306       h_phiGen1.Fill( gen->phi() );
01307       jetInd++;
01308     }
01309   }
01310 
01311   LeadMassP1 = (p4tmp[0]+p4tmp[1]).mass();
01312   dijetMassP1.Fill( LeadMassP1 ); 
01313   if (EtaOk10 == 2) {
01314     dijetMassP101.Fill( LeadMassP1 ); 
01315     dijetMassP_700_101.Fill( LeadMassP1 ); 
01316     dijetMassP_2000_101.Fill( LeadMassP1 ); 
01317     dijetMassP_5000_101.Fill( LeadMassP1 ); 
01318    }
01319   if (EtaOk13 == 2) {
01320     dijetMassP131.Fill( LeadMassP1 ); 
01321     dijetMassP_700_131.Fill( LeadMassP1 ); 
01322     dijetMassP_2000_131.Fill( LeadMassP1 ); 
01323     dijetMassP_5000_131.Fill( LeadMassP1 ); 
01324 
01325    }
01326   if (EtaOk40 == 2) {
01327     dijetMassP401.Fill( LeadMassP1 );
01328     dijetMassP_5000_401.Fill( LeadMassP1 );
01329     dijetMassP_5000_401.Fill( LeadMassP1 );
01330     dijetMassP_5000_401.Fill( LeadMassP1 );
01331   } 
01332 
01333   //  h_nGenJets1.Fill( genJets1->size() ); 
01334   h_nGenJets1.Fill( jetInd ); 
01335 
01336   //**** Get the GenJet2 collection
01337   Handle<GenJetCollection> genJets2;
01338   evt.getByLabel( GenJetAlgorithm2, genJets2 );
01339 
01340   //Loop over the two leading GenJets and fill some histograms
01341   jetInd    = 0;
01342   allJetInd = 0;
01343   for( GenJetCollection::const_iterator gen = genJets2->begin(); gen != genJets2->end(); ++ gen ) {
01344     allJetInd++;
01345     if (allJetInd == 1) {
01346       p4tmp[0] = gen->p4();
01347     }
01348     if (allJetInd == 2) {
01349       p4tmp[1] = gen->p4();
01350     }
01351     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01352       h_ptGenL2.Fill( gen->pt() );   
01353       h_ptGenL22.Fill( gen->pt() );   
01354       h_ptGenL23.Fill( gen->pt() );   
01355 
01356       h_etaGenL2.Fill( gen->eta() );
01357       h_phiGenL2.Fill( gen->phi() );
01358     }
01359 
01360     if ( gen->pt() > minJetPt) {
01361       // std::cout << "GEN JET2 #" << jetInd << std::endl << gen->print() << std::endl;
01362       h_ptGen2.Fill( gen->pt() );   
01363       h_ptGen22.Fill( gen->pt() );   
01364       h_ptGen23.Fill( gen->pt() );   
01365 
01366       h_etaGen2.Fill( gen->eta() );
01367       h_phiGen2.Fill( gen->phi() );
01368       jetInd++;
01369     }
01370   }
01371   
01372   LeadMassP2 = (p4tmp[0]+p4tmp[1]).mass();
01373   dijetMassP2.Fill( LeadMassP2 ); 
01374 
01375 
01376   //  h_nGenJets2.Fill( genJets2->size() ); 
01377   h_nGenJets2.Fill( jetInd ); 
01378 
01379   //**** Get the GenJet3 collection
01380   Handle<GenJetCollection> genJets3;
01381   evt.getByLabel( GenJetAlgorithm3, genJets3 );
01382 
01383   //Loop over the two leading GenJets and fill some histograms
01384   jetInd    = 0;
01385   allJetInd = 0;
01386   for( GenJetCollection::const_iterator gen = genJets3->begin(); gen != genJets3->end(); ++ gen ) {
01387     allJetInd++;
01388     if (allJetInd == 1) {
01389       p4tmp[0] = gen->p4();
01390     }
01391     if (allJetInd == 2) {
01392       p4tmp[1] = gen->p4();
01393     }    
01394     if ( (allJetInd == 1) || (allJetInd == 2) ) {
01395       h_ptGenL3.Fill( gen->pt() );   
01396       h_ptGenL32.Fill( gen->pt() );   
01397       h_ptGenL33.Fill( gen->pt() );   
01398 
01399       h_etaGenL3.Fill( gen->eta() );
01400       h_phiGenL3.Fill( gen->phi() );
01401     }
01402 
01403     if ( gen->pt() > minJetPt) {
01404       // std::cout << "GEN JET3 #" << jetInd << std::endl << gen->print() << std::endl;
01405       h_ptGen3.Fill( gen->pt() );   
01406       h_ptGen32.Fill( gen->pt() );   
01407       h_ptGen33.Fill( gen->pt() );   
01408 
01409       h_etaGen3.Fill( gen->eta() );
01410       h_phiGen3.Fill( gen->phi() );
01411       jetInd++;
01412     }
01413   }
01414 
01415   LeadMassP3 = (p4tmp[0]+p4tmp[1]).mass();
01416   dijetMassP3.Fill( LeadMassP3 ); 
01417 
01418 
01419   //  h_nGenJets3.Fill( genJets3->size() ); 
01420   h_nGenJets3.Fill( jetInd ); 
01421 
01422 
01423   // *********************
01424   // MidPoint Jet Matching
01425   
01426   Handle<GenJetCollection>  genJets;
01427   Handle<CaloJetCollection> caloJets;
01428 
01429   //  evt.getByLabel( "midPointCone5GenJets",  genJets );
01430   //  evt.getByLabel( "midPointCone5CaloJets", caloJets );
01431   evt.getByLabel( GenJetAlgorithm1, genJets );
01432   evt.getByLabel( CaloJetAlgorithm1, caloJets );
01433 
01434 
01435   int maxJets = MAXJETS;
01436 
01437   jetInd = 0;
01438   double dRmin[MAXJETS];
01439   math::XYZTLorentzVector  p4gen[MAXJETS], p4cal[MAXJETS],  
01440     p4par[MAXJETS], p4Zp[MAXJETS], p4part[MAXJETS];
01441 
01442   int used[MAXJETS];
01443   int nj;
01444 
01445   for( int i=0; i<maxJets; ++i ) used[i] = 0;  
01446 
01447   //  cout << ">>>>>>>>> " << endl;
01448 
01449 
01450 
01451   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01452        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01453 
01454     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01455     dRmin[jetInd] = 1000.0;
01456 
01457     nj      = 0;
01458     usedInd = -1;
01459     
01460     for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) { 
01461  
01462       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01463  
01464       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01465         dRmin[jetInd] = delR;        // delta R of match
01466         p4cal[jetInd] = cal->p4();   // Matched Cal 4-vector
01467         usedInd       = nj;     
01468       }
01469 
01470       nj++;
01471     }
01472 
01473     if (fabs(p4gen[jetInd].eta()) < 1.3) {
01474       matchedAllPt11.Fill(p4gen[jetInd].pt());
01475     }
01476     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01477       matchedAllPt12.Fill(p4gen[jetInd].pt());
01478     }
01479     if (fabs(p4gen[jetInd].eta()) > 3.) {
01480       matchedAllPt13.Fill(p4gen[jetInd].pt());     
01481     }
01482 
01483     if (usedInd != -1) {
01484 
01485       used[usedInd] = 1;
01486 
01487       if (p4cal[jetInd].pt() > minJetPt10)
01488         hf_PtResponse1.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01489 
01490       if (fabs(p4gen[jetInd].eta()) < 1.3) {
01491         matchedPt11.Fill(p4gen[jetInd].pt());
01492       } 
01493       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01494         matchedPt12.Fill(p4gen[jetInd].pt());
01495       }
01496       if (fabs(p4gen[jetInd].eta()) > 3.) {
01497         matchedPt13.Fill(p4gen[jetInd].pt());
01498       }
01499       /***
01500       cout << " >>>SISCone "     
01501            << jetInd                 << " " 
01502            << p4gen[jetInd].pt()     << " " 
01503            << p4gen[jetInd].eta()    << " " 
01504            << p4gen[jetInd].phi()    << " " 
01505            << p4cal[jetInd].pt()     << " " 
01506            << p4cal[jetInd].eta()    << " " 
01507            << p4cal[jetInd].phi()    
01508            << endl;
01509       ***/
01510 
01511       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01512       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01513         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01514           dPt20Frac1.Fill(dpt/p4gen[jetInd].pt());
01515         }
01516         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01517           dPt40Frac1.Fill(dpt/p4gen[jetInd].pt());
01518         }
01519         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01520           dPt80Frac1.Fill(dpt/p4gen[jetInd].pt());
01521         }
01522         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01523           dPt100Frac1.Fill(dpt/p4gen[jetInd].pt());
01524         }
01525       }
01526       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) ) {
01527 
01528         dR1.Fill(dRmin[jetInd]);
01529         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01530         dPhi1.Fill(dphi);
01531         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01532         dEta1.Fill(deta);
01533         dPt1.Fill(dpt);         
01534         dPtFrac1.Fill(dpt/p4gen[jetInd].pt());
01535 
01536         if ( ( (dpt/p4gen[jetInd].pt()) < -0.5 ) &&  ( fabs(dpt) > 90. ) ) {
01537 
01538           cout << " deltaR min = "     << dRmin[jetInd]
01539                << " Ind = " << jetInd  << " / " << usedInd << " / " << used[nj]
01540                << " Del pt / frac  = " << dpt
01541                << " / "                << dpt/p4gen[jetInd].pt()
01542                << " cal/gen pt   = "   << p4cal[jetInd].pt()
01543                << " / "                << p4gen[jetInd].pt()
01544                << " cal/gen eta  = "   << p4cal[jetInd].eta()
01545                << " / "                << p4gen[jetInd].eta()
01546                << " cal/gen phi  = "   << p4cal[jetInd].phi()
01547                << " / "                << p4gen[jetInd].phi()
01548                << endl;
01549         }
01550 
01551       }
01552       
01553       jetInd++;    
01554 
01555     }
01556 
01557   }
01558 
01559 
01560 
01561   // *********************
01562   // Seedless Jet Matching
01563   
01564   //  Handle<GenJetCollection>  genJets;
01565   //  Handle<CaloJetCollection> caloJets;
01566 
01567   //  evt.getByLabel( "sisCone5GenJets",  genJets );
01568   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
01569   evt.getByLabel( GenJetAlgorithm2, genJets );
01570   evt.getByLabel( CaloJetAlgorithm2, caloJets );
01571 
01572   // int maxJets = 20;
01573   jetInd = 0;
01574   //  double dRmin[20];
01575   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01576 
01577   for( int i=0; i<maxJets; ++i ) used[i] = 0;    
01578   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01579        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01580     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01581     dRmin[jetInd] = 1000.0;
01582 
01583     nj      = 0;
01584     usedInd = -1;
01585 
01586     for( CaloJetCollection::const_iterator cal = caloJets->begin(); 
01587          cal != caloJets->end(); ++ cal ) { 
01588       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01589 
01590       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01591         dRmin[jetInd] = delR;       // delta R of match
01592         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
01593         usedInd       = nj;
01594       }
01595       nj++;
01596     }
01597 
01598     if (fabs(p4gen[jetInd].eta()) < 1.3)   {      
01599       matchedAllPt21.Fill(p4gen[jetInd].pt());
01600     }
01601     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01602       matchedAllPt22.Fill(p4gen[jetInd].pt());
01603     }
01604     if (fabs(p4gen[jetInd].eta()) > 3.) {
01605       matchedAllPt23.Fill(p4gen[jetInd].pt());
01606     }
01607 
01608 
01609     if (usedInd != -1) {
01610 
01611       used[usedInd] = 1;
01612 
01613       if (p4cal[jetInd].pt() > minJetPt10)
01614         hf_PtResponse2.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01615 
01616       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01617         matchedPt21.Fill(p4gen[jetInd].pt());
01618       }
01619       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01620         matchedPt22.Fill(p4gen[jetInd].pt());
01621       }
01622       if (fabs(p4gen[jetInd].eta()) > 3.) {
01623         matchedPt23.Fill(p4gen[jetInd].pt());
01624       }
01625       /***
01626       cout << " >>>IterCone "     
01627            << jetInd                 << " " 
01628            << p4gen[jetInd].pt()     << " " 
01629            << p4gen[jetInd].eta()    << " " 
01630            << p4gen[jetInd].phi()    << " " 
01631            << p4cal[jetInd].pt()     << " " 
01632            << p4cal[jetInd].eta()    << " " 
01633            << p4cal[jetInd].phi()    
01634            << endl;
01635       ***/
01636 
01637 
01638       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01639       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01640         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01641           dPt20Frac2.Fill(dpt/p4gen[jetInd].pt());
01642         }
01643         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01644           dPt40Frac2.Fill(dpt/p4gen[jetInd].pt());
01645         }
01646         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01647           dPt80Frac2.Fill(dpt/p4gen[jetInd].pt());
01648         }
01649         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01650           dPt100Frac2.Fill(dpt/p4gen[jetInd].pt());
01651         }
01652       }
01653       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) )  {
01654 
01655         dR2.Fill(dRmin[jetInd]);
01656         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01657         dPhi2.Fill(dphi);
01658         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01659         dEta2.Fill(deta);
01660         dPt2.Fill(dpt);
01661         dPtFrac2.Fill(dpt/p4gen[jetInd].pt());       
01662 
01663       }      
01664 
01665       jetInd++;    
01666     }
01667 
01668   }
01669 
01670   // *********************
01671   // Kt Jet Matching
01672   
01673   //  Handle<GenJetCollection>  genJets;
01674   //  Handle<CaloJetCollection> caloJets;
01675 
01676   //  evt.getByLabel( "sisCone5GenJets",  genJets );
01677   //  evt.getByLabel( "sisCone5CaloJets", caloJets );
01678   evt.getByLabel( GenJetAlgorithm3, genJets );
01679   evt.getByLabel( CaloJetAlgorithm3, caloJets );
01680 
01681   // int maxJets = 20;
01682   jetInd = 0;
01683   //  double dRmin[20];
01684   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01685 
01686   for( int i=0; i<maxJets; ++i ) used[i] = 0;
01687   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01688        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01689     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01690     dRmin[jetInd] = 1000.0;
01691 
01692     nj = 0;
01693     usedInd = -1;
01694 
01695     for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) { 
01696       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01697 
01698       if ( (delR<dRmin[jetInd]) &&  (delR < matchedDelR) && (used[nj] == 0) ) {
01699         dRmin[jetInd] = delR;        // delta R of match
01700         p4cal[jetInd] = cal->p4();   // Matched Cal 4-vector
01701         usedInd       = nj;
01702       }
01703       nj++;
01704     }
01705 
01706     if (fabs(p4gen[jetInd].eta()) < 1.3) {
01707       matchedAllPt31.Fill(p4gen[jetInd].pt());
01708     }
01709     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01710       matchedAllPt32.Fill(p4gen[jetInd].pt());
01711     }
01712    if (fabs(p4gen[jetInd].eta()) > 3.) {
01713       matchedAllPt33.Fill(p4gen[jetInd].pt());
01714     }
01715 
01716 
01717     if (usedInd != -1) {
01718       used[usedInd] = 1;
01719 
01720       if (p4cal[jetInd].pt() > minJetPt10)
01721         hf_PtResponse3.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01722 
01723       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01724         matchedPt31.Fill(p4gen[jetInd].pt());
01725       }
01726       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01727         matchedPt32.Fill(p4gen[jetInd].pt());
01728       }
01729       if (fabs(p4gen[jetInd].eta()) > 3.) {
01730         matchedPt33.Fill(p4gen[jetInd].pt());
01731       }
01732       /***
01733       cout << " >>>MidPoint "     
01734            << jetInd                 << " " 
01735            << p4gen[jetInd].pt()     << " " 
01736            << p4gen[jetInd].eta()    << " " 
01737            << p4gen[jetInd].phi()    << " " 
01738            << p4cal[jetInd].pt()     << " " 
01739            << p4cal[jetInd].eta()    << " " 
01740            << p4cal[jetInd].phi()    
01741            << endl;
01742       ***/
01743 
01744       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01745       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01746         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01747           dPt20Frac3.Fill(dpt/p4gen[jetInd].pt());
01748         }
01749         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01750           dPt40Frac3.Fill(dpt/p4gen[jetInd].pt());
01751         }
01752         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01753           dPt80Frac3.Fill(dpt/p4gen[jetInd].pt());
01754         }
01755         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01756           dPt100Frac3.Fill(dpt/p4gen[jetInd].pt());
01757         }
01758       }
01759 
01760       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) ) {
01761 
01762         dR3.Fill(dRmin[jetInd]);
01763         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01764         dPhi3.Fill(dphi);
01765         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01766         dEta3.Fill(deta);
01767         dPt3.Fill(dpt);
01768         dPtFrac3.Fill(dpt/p4gen[jetInd].pt());
01769 
01770 
01771       }
01772       
01773       jetInd++;    
01774     }
01775   }
01776 
01777   // *********************
01778   // Jet Algorithm 4
01779   
01780   evt.getByLabel( GenJetAlgorithm4,  genJets );
01781   evt.getByLabel( CaloJetAlgorithm4, caloJets );
01782 
01783   // int maxJets = 20;
01784   jetInd = 0;
01785   //  double dRmin[20];
01786   //  math::XYZTLorentzVector p4jet[20], p4gen[20], p4cal[20], p4cor[20];
01787 
01788   for( int i=0; i<maxJets; ++i ) used[i] = 0;    
01789   for( GenJetCollection::const_iterator gen = genJets->begin(); 
01790        gen != genJets->end() && jetInd<maxJets; ++ gen ) { 
01791     p4gen[jetInd] = gen->p4();    //Gen 4-vector
01792     dRmin[jetInd] = 1000.0;
01793 
01794     nj      = 0;
01795     usedInd = -1;
01796 
01797     for( CaloJetCollection::const_iterator cal = caloJets->begin(); 
01798          cal != caloJets->end(); ++ cal ) { 
01799       double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
01800 
01801       if ( (delR<dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
01802         dRmin[jetInd] = delR;       // delta R of match
01803         p4cal[jetInd] = cal->p4();  // Matched Cal 4-vector
01804         usedInd       = nj;
01805       }
01806       nj++;
01807     }
01808 
01809     if (fabs(p4gen[jetInd].eta()) < 1.3)   {      
01810       matchedAllPt41.Fill(p4gen[jetInd].pt());
01811     }
01812     if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01813       matchedAllPt42.Fill(p4gen[jetInd].pt());
01814     }
01815     if (fabs(p4gen[jetInd].eta()) > 3.) {
01816       matchedAllPt43.Fill(p4gen[jetInd].pt());
01817     }
01818 
01819 
01820     if (usedInd != -1) {
01821 
01822       used[usedInd] = 1;
01823 
01824       if (p4cal[jetInd].pt() > minJetPt10)
01825         hf_PtResponse4.Fill(p4cal[jetInd].eta(), p4cal[jetInd].pt()/p4gen[jetInd].pt());
01826 
01827       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01828         matchedPt41.Fill(p4gen[jetInd].pt());
01829       }
01830       if ((fabs(p4gen[jetInd].eta()) > 1.3) && (fabs(p4gen[jetInd].eta()) < 3.) ) {
01831         matchedPt42.Fill(p4gen[jetInd].pt());
01832       }
01833       if (fabs(p4gen[jetInd].eta()) > 3.) {
01834         matchedPt43.Fill(p4gen[jetInd].pt());
01835       }
01836 
01837       /***
01838       cout << " >>>KtClus "     
01839            << jetInd                 << " " 
01840            << p4gen[jetInd].pt()     << " " 
01841            << p4gen[jetInd].eta()    << " " 
01842            << p4gen[jetInd].phi()    << " " 
01843            << p4cal[jetInd].pt()     << " " 
01844            << p4cal[jetInd].eta()    << " " 
01845            << p4cal[jetInd].phi()    
01846            << endl;
01847       ***/
01848 
01849 
01850       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01851       if (fabs(p4gen[jetInd].eta()) < 1.3)  {
01852         if ( (p4gen[jetInd].pt() > 20) && (p4gen[jetInd].pt() < 40) ) {
01853           dPt20Frac4.Fill(dpt/p4gen[jetInd].pt());
01854         }
01855         if ( (p4gen[jetInd].pt() > 40) && (p4gen[jetInd].pt() < 60) ) {
01856           dPt40Frac4.Fill(dpt/p4gen[jetInd].pt());
01857         }
01858         if ( (p4gen[jetInd].pt() > 80) && (p4gen[jetInd].pt() < 100) ) {
01859           dPt80Frac4.Fill(dpt/p4gen[jetInd].pt());
01860         }
01861         if ( (p4gen[jetInd].pt() > 100) && (p4gen[jetInd].pt() < 120) ) {
01862           dPt100Frac4.Fill(dpt/p4gen[jetInd].pt());
01863         }
01864       }
01865 
01866       if ( (p4gen[jetInd].pt() > minJetPt10) && (fabs(p4gen[jetInd].eta()) < 1.3) )  {
01867 
01868         dR4.Fill(dRmin[jetInd]);
01869         double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01870         dPhi4.Fill(dphi);
01871         double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01872         dEta4.Fill(deta);
01873         dPt4.Fill(dpt);
01874         dPtFrac4.Fill(dpt/p4gen[jetInd].pt());       
01875 
01876       }      
01877 
01878       jetInd++;    
01879     }
01880 
01881   }
01882 
01883 
01884 
01885   // *********************
01886   // MidPoint - Seedless Jet Matching
01887   
01888   Handle<CaloJetCollection> calo1Jets;
01889   Handle<CaloJetCollection> calo2Jets;
01890   Handle<CaloJetCollection> calo3Jets;
01891 
01892   evt.getByLabel( CaloJetAlgorithm1, calo1Jets );
01893   evt.getByLabel( CaloJetAlgorithm2, calo2Jets );
01894   evt.getByLabel( CaloJetAlgorithm3, calo3Jets );
01895 
01896   jetInd = 0;
01897 
01898   for( int i=0; i<maxJets; ++i ) used[i] = 0;
01899   for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
01900        cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
01901 
01902     p4gen[jetInd] = cal1->p4();    //Gen 4-vector
01903     dRmin[jetInd] = 1000.0;
01904 
01905     nj = 0;
01906     for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); cal2 != calo2Jets->end(); ++cal2 ) { 
01907 
01908       double delR = deltaR( cal1->eta(), cal1->phi(), cal2->eta(), cal2->phi() ); 
01909       if ( (delR<dRmin[jetInd]) && (used[nj] == 0) ) {
01910         dRmin[jetInd] = delR;        // delta R of match
01911         p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
01912         usedInd       = nj;
01913       }
01914       nj++;
01915     }
01916     used[usedInd] = 1;
01917 
01918     if (p4gen[jetInd].pt() > minJetPt) {
01919       dR12.Fill(dRmin[jetInd]);
01920       double dphi = deltaPhi(p4cal[jetInd].phi(), p4gen[jetInd].phi());
01921       dPhi12.Fill(dphi);
01922       double deta = p4cal[jetInd].eta() - p4gen[jetInd].eta();
01923       dEta12.Fill(deta);
01924       double dpt = p4cal[jetInd].pt() - p4gen[jetInd].pt();
01925       dPt12.Fill(dpt);
01926     }
01927 
01928     jetInd++;    
01929   }
01930 
01931   // ******************************************
01932   // ******************************************
01933 
01934 
01935   Handle<CandidateCollection> genParticles;
01936   evt.getByLabel("genParticleCandidates",genParticles);
01937 
01938 
01939   // *********************
01940   // Partons (Z')
01941 
01942   int nPart = 0;
01943   for (size_t i =0;i< genParticles->size(); i++) {
01944 
01945     const Candidate &p = (*genParticles)[i];
01946     //    int Status =  p.status();
01947     //    bool ParticleIsStable = Status==1;
01948     int id = p.pdgId();
01949 
01950     if (id == 32) {
01951 
01952       if (p.numberOfDaughters() != 0) {
01953         /***
01954         cout << "Z': part = "   << i << " id = " << id 
01955              << " daughters = " << p.numberOfDaughters() 
01956              << " mass = "      << p.mass()
01957              << endl;
01958         ***/
01959         ZpMG =  p.mass();
01960         ZpMassGen.Fill( ZpMG );
01961         if (EtaOk10 == 2) {
01962           ZpMassGen10.Fill( ZpMG );
01963           ZpMassGen_700_10.Fill( ZpMG );
01964           ZpMassGen_2000_10.Fill( ZpMG );
01965           ZpMassGen_5000_10.Fill( ZpMG );
01966         }
01967         if (EtaOk13 == 2) {
01968           ZpMassGen13.Fill( ZpMG );
01969           ZpMassGen_700_13.Fill( ZpMG );
01970           ZpMassGen_2000_13.Fill( ZpMG );
01971           ZpMassGen_5000_13.Fill( ZpMG );
01972         }
01973         if (EtaOk40 == 2) {
01974           ZpMassGen40.Fill( ZpMG );
01975           ZpMassGen_700_40.Fill( ZpMG );
01976           ZpMassGen_2000_40.Fill( ZpMG );
01977           ZpMassGen_5000_40.Fill( ZpMG );
01978         }
01979       }
01980 
01981       for( int id1=0, nd1=p.numberOfDaughters(); id1 < nd1; ++id1 ) {
01982 
01983         const Candidate * d1 = p.daughter(id1);
01984 
01985         if ( abs(d1->pdgId()) != 32) {  
01986           math::XYZTLorentzVector momentum=d1->p4();
01987           p4Zp[nPart] = momentum=d1->p4();
01988           nPart++;
01989         }
01990 
01991       }
01992     }
01993 
01994   }
01995 
01996   // *********************
01997   // Match jets to Zp
01998   int genInd = 0;
01999 
02000   if (nPart == 2) {
02001     
02002     ZpM = (p4Zp[0]+p4Zp[1]).mass();
02003     ZpMass.Fill( ZpM );
02004 
02005     if (EtaOk10 == 2) {
02006       ZpMass_700_10.Fill( ZpM );
02007       ZpMass_2000_10.Fill( ZpM );
02008       ZpMass_5000_10.Fill( ZpM );
02009     }
02010     if (EtaOk13 == 2) {
02011       ZpMass_700_13.Fill( ZpM );
02012       ZpMass_2000_13.Fill( ZpM );
02013       ZpMass_5000_13.Fill( ZpM );
02014     }
02015     if (EtaOk40 == 2) {
02016       ZpMass_700_40.Fill( ZpM );
02017       ZpMass_2000_40.Fill( ZpM );
02018       ZpMass_5000_40.Fill( ZpM );
02019     }
02020 
02021     int usedInd;   
02022 
02023     // ***********
02024     // **** Calor1
02025     usedInd = -1;
02026     jetInd  = 0;
02027 
02028     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02029     for( int i=0; i<2; ++i ) { 
02030       
02031       dRmin[jetInd]  = 1000.0;
02032 
02033       int nj = 0;
02034       for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
02035            cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
02036         
02037         double delR = deltaR( cal1->eta(), cal1->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02038 
02039         //      if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02040         if ( (delR < dRmin[jetInd]) && (used[nj] == 0) ) {
02041           dRmin[jetInd] = delR;        // delta R of match
02042           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
02043           usedInd       = nj;
02044           genInd        = i;
02045         }
02046 
02047         /****
02048         cout << "Delta R = " << delR
02049              << " deltaR min = " << dRmin[jetInd]
02050              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02051              << " cal1 eta  = " << cal1->eta()
02052              << " p4par eta = " << p4Zp[i].eta()
02053              << " cal1 phi  = " << cal1->phi()
02054              << " p4par phi = " << p4Zp[i].phi()
02055              << endl;
02056         cout << "    " 
02057              << " p4par = " << p4Zp[i]
02058              << " p4cal = " << cal1->p4()
02059              << endl;
02060         ***/
02061                 
02062         nj++;
02063       }
02064 
02065       // Found matched jet
02066       if (usedInd != -1) {
02067         used[usedInd] = 1;
02068         jetInd++;    
02069       }
02070 
02071     }
02072     
02073     ZpMM = (p4cal[0]+p4cal[1]).mass();
02074     ZpMassMatched1.Fill( ZpMM );
02075 
02076     if ((ZpMG != 0) && (EtaOk40 == 2)) {
02077       ZpMassRes401.Fill( (ZpMM - ZpMG) / ZpMG );
02078 
02079       ZpMassResL401.Fill( (LeadMass1 - ZpMG) / ZpMG );
02080       ZpMassResL402.Fill( (LeadMass2 - ZpMG) / ZpMG );
02081       ZpMassResL403.Fill( (LeadMass3 - ZpMG) / ZpMG );
02082       
02083       ZpMassResRL401.Fill( LeadMass1 / ZpMG );
02084       ZpMassResRL402.Fill( LeadMass2 / ZpMG );
02085       ZpMassResRL403.Fill( LeadMass3 / ZpMG );
02086 
02087       ZpMassResRLoP401.Fill( LeadMass1 / LeadMassP1 );
02088       ZpMassResRLoP402.Fill( LeadMass2 / LeadMassP2 );
02089       ZpMassResRLoP403.Fill( LeadMass3 / LeadMassP2 );
02090 
02091       ZpMassResPRL401.Fill( LeadMassP1 / ZpMG );
02092       ZpMassResPRL402.Fill( LeadMassP2 / ZpMG );
02093       ZpMassResPRL403.Fill( LeadMassP3 / ZpMG );
02094       
02095     }
02096 
02097     if ((ZpMG != 0) && (EtaOk10 == 2)) {
02098       ZpMassRes101.Fill( (ZpMM - ZpMG) / ZpMG );
02099 
02100       ZpMassResL101.Fill( (LeadMass1 - ZpMG) / ZpMG );
02101       ZpMassResL102.Fill( (LeadMass2 - ZpMG) / ZpMG );
02102       ZpMassResL103.Fill( (LeadMass3 - ZpMG) / ZpMG );
02103       
02104       ZpMassResRL101.Fill( LeadMass1 / ZpMG );
02105       ZpMassResRL102.Fill( LeadMass2 / ZpMG );
02106       ZpMassResRL103.Fill( LeadMass3 / ZpMG );
02107 
02108       ZpMassResRLoP101.Fill( LeadMass1 / LeadMassP1 );
02109       ZpMassResRLoP102.Fill( LeadMass2 / LeadMassP2 );
02110       ZpMassResRLoP103.Fill( LeadMass3 / LeadMassP2 );
02111 
02112       ZpMassResPRL101.Fill( LeadMassP1 / ZpMG );
02113       ZpMassResPRL102.Fill( LeadMassP2 / ZpMG );
02114       ZpMassResPRL103.Fill( LeadMassP3 / ZpMG );
02115       
02116     }
02117 
02118     if ((ZpMG != 0) && (EtaOk13 == 2)) {
02119       ZpMassRes131.Fill( (ZpMM - ZpMG) / ZpMG );
02120 
02121       ZpMassResL131.Fill( (LeadMass1 - ZpMG) / ZpMG );
02122       ZpMassResL132.Fill( (LeadMass2 - ZpMG) / ZpMG );
02123       ZpMassResL133.Fill( (LeadMass3 - ZpMG) / ZpMG );
02124       
02125       ZpMassResRL131.Fill( LeadMass1 / ZpMG );
02126       ZpMassResRL132.Fill( LeadMass2 / ZpMG );
02127       ZpMassResRL133.Fill( LeadMass3 / ZpMG );
02128 
02129       ZpMassResRLoP131.Fill( LeadMass1 / LeadMassP1 );
02130       ZpMassResRLoP132.Fill( LeadMass2 / LeadMassP2 );
02131       ZpMassResRLoP133.Fill( LeadMass3 / LeadMassP2 );
02132 
02133       ZpMassResPRL131.Fill( LeadMassP1 / ZpMG );
02134       ZpMassResPRL132.Fill( LeadMassP2 / ZpMG );
02135       ZpMassResPRL133.Fill( LeadMassP3 / ZpMG );
02136       
02137     }
02138 
02139 
02140 
02141     // ***********
02142     // **** Calor2
02143     usedInd = -1;
02144     jetInd  = 0;
02145 
02146     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02147     for( int i=0; i<2; ++i ) { 
02148       
02149       dRmin[jetInd]  = 1000.0;
02150 
02151       int nj = 0;
02152       for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); 
02153            cal2 != calo2Jets->end() && jetInd<maxJets; ++cal2 ) { 
02154         
02155         double delR = deltaR( cal2->eta(), cal2->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02156 
02157         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02158           dRmin[jetInd] = delR;        // delta R of match
02159           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
02160           usedInd       = nj;
02161         }
02162 
02163         /****   
02164         cout << "Delta R = " << delR
02165              << " deltaR min = " << dRmin[jetInd]
02166              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02167              << " p4par = " << p4par[i]
02168              << " p4cal = " << cal1->p4()
02169              << " cal1 eta  = " << cal1->eta()
02170              << " p4par eta = " << p4par[i].eta()
02171              << endl;
02172         ****/
02173                 
02174         nj++;
02175       }
02176 
02177       // Found matched jet
02178       if (usedInd != -1) {
02179         used[usedInd] = 1;
02180         jetInd++;    
02181       }
02182 
02183     }
02184     
02185     ZpMM = (p4cal[0]+p4cal[1]).mass();
02186     ZpMassMatched2.Fill( ZpMM );
02187     ZpMassRes402.Fill( (ZpMM - ZpM) / ZpM );
02188 
02189 
02190     // ***********
02191     // **** Calor3
02192     usedInd = -1;
02193     jetInd  = 0;
02194 
02195     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02196     for( int i=0; i<2; ++i ) { 
02197       
02198       dRmin[jetInd]  = 1000.0;
02199 
02200       int nj = 0;
02201       for( CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); 
02202            cal3 != calo3Jets->end() && jetInd<maxJets; ++cal3 ) { 
02203         
02204         double delR = deltaR( cal3->eta(), cal3->phi(), p4Zp[i].eta(), p4Zp[i].phi() ); 
02205 
02206         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02207           dRmin[jetInd] = delR;        // delta R of match
02208           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
02209           usedInd       = nj;
02210         }
02211 
02212         /****   
02213         cout << "Delta R = " << delR
02214              << " deltaR min = " << dRmin[jetInd]
02215              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02216              << " p4par = " << p4par[i]
02217              << " p4cal = " << cal1->p4()
02218              << " cal1 eta  = " << cal1->eta()
02219              << " p4par eta = " << p4par[i].eta()
02220              << endl;
02221         ****/
02222                 
02223         nj++;
02224       }
02225 
02226 
02227       // Found matched jet
02228       if (usedInd != -1) {
02229         used[usedInd] = 1;
02230         jetInd++;    
02231       }
02232 
02233     }
02234     
02235     ZpMM = (p4cal[0]+p4cal[1]).mass();
02236     ZpMassMatched3.Fill( ZpMM );
02237     ZpMassRes403.Fill( (ZpMM - ZpM) / ZpM );
02238     
02239   } else {
02240     cout << "Z' (3): nPart = " << nPart << endl;
02241   }
02242 
02243 
02244 
02245   // *********************
02246   // Partons (ttbar) Jet Matching
02247 
02248   //  cout << ">>> Begin MC list" << endl;
02249   int nJet = 0;
02250 
02251 
02252   int ii = 1;
02253   int jj = 4;
02254 
02255   for (size_t i =0;i< genParticles->size(); i++) {
02256 
02257     const Candidate &p = (*genParticles)[i];
02258     //    int Status =  p.status();
02259     //    bool ParticleIsStable = Status==1;
02260     int id = p.pdgId();
02261 
02262     // Top Quarks
02263     if (abs(id) == 6) {
02264       cout << "TOP: id = " << id << " mass = " << p.mass() << endl;
02265 
02266       topMassParton.Fill(p.mass());
02267 
02268       if (id == 6)  tMassGen.Fill(p.mass());
02269       if (id == -6) tbarMassGen.Fill(p.mass());
02270 
02271       for( size_t id1=0, nd1=p.numberOfDaughters(); id1 < nd1; ++id1 ) {
02272 
02273         const Candidate * d1 = p.daughter(id1);
02274 
02275         // b - quark
02276         if ( abs(d1->pdgId()) == 5) {
02277 
02278           math::XYZTLorentzVector momentum=d1->p4();
02279           p4par[nJet++] = momentum=d1->p4();
02280 
02281           cout << "Daughter1: id = " << d1->pdgId() 
02282                << " daughters = " << d1->numberOfDaughters() 
02283                << " mother 1   = " << (d1->mother())->pdgId() 
02284                << " Momentum " << momentum << " GeV/c"
02285                << endl;   
02286 
02287           
02288           if ( (d1->mother())->pdgId() == 6 ) {
02289             p4part[0] = momentum=d1->p4();
02290             cout << ">>> part0 = " << p4part[0] << endl;
02291           }
02292           if ( (d1->mother())->pdgId() == -6) {
02293             p4part[3] = momentum=d1->p4();
02294             cout << ">>> part3 = " << p4part[3] << endl;
02295           }
02296 
02297         }
02298 
02299         
02300         // W
02301         // Check for fully hadronic decay
02302 
02303         if ( abs(d1->pdgId()) == 24) {
02304 
02305           for( size_t id2=0, nd2=d1->numberOfDaughters(); id2 < nd2; ++id2 ) {
02306 
02307             const Candidate * d2 = d1->daughter(id2);
02308 
02309             if (abs(d2->pdgId()) < 9) {
02310 
02311               math::XYZVector vertex(d2->vx(),d2->vy(),d2->vz());
02312               math::XYZTLorentzVector momentum=d2->p4();
02313               p4par[nJet++] = momentum=d2->p4();
02314 
02315               if ( (d1->mother())->pdgId() == 6 ) {
02316                 p4part[ii] = momentum=d2->p4();
02317                 cout << ">>> part" << ii << " = " << p4part[ii] << endl;
02318                 ii++;
02319               }
02320               if ( (d1->mother())->pdgId() == -6 ) {
02321                 p4part[jj] = momentum=d2->p4();
02322                 cout << ">>> part" << jj << " = " << p4part[jj] << endl;
02323                 jj++;
02324               }
02325 
02326               cout << "Daughter2: id = " << d2->pdgId() 
02327                    << " daughters = " << d2->numberOfDaughters() 
02328                    << " mother 2   = " << (d2->mother())->pdgId() 
02329                    << " Momentum " << momentum << " GeV/c"
02330                    << endl;
02331             }
02332 
02333           }
02334         }
02335 
02336         //      if ( pdgId == d->pdgId() && d->status() == 1 ) {        
02337         //      }
02338       }
02339     }
02340 
02341   }
02342   //  cout << ">>> N Jets = " << nJet << endl;
02343     
02344   if (nJet == 6) {
02345 
02346     double tmass    = (p4part[0]+p4part[1]+p4part[2]).mass();
02347     double tbarmass = (p4part[3]+p4part[4]+p4part[5]).mass();
02348 
02349     tMass.Fill(tmass);
02350     tbarMass.Fill(tbarmass);
02351     
02352     cout << ">>> T Mass = " << tmass << " / " << tbarmass << endl;
02353 
02354     double mindR = 1000.;
02355     for( size_t i=0; i<6; ++i ) { 
02356       for( size_t j=0; j<6; ++j ) { 
02357         if (j > i) {
02358           double delR = deltaR( p4par[i].eta(), p4par[i].phi(), p4par[j].eta(), p4par[j].phi() ); 
02359           if (delR < mindR) mindR = delR;
02360           dRParton.Fill(delR);
02361         }
02362       }
02363     }
02364     dRPartonMin.Fill(mindR);
02365 
02366     int usedInd;
02367     usedInd = -1;
02368     jetInd  = 0;
02369     
02370     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02371     for( int i=0; i<6; ++i ) { 
02372       
02373       dRmin[jetInd]  = 1000.0;
02374 
02375       int nj = 0;
02376       for( CaloJetCollection::const_iterator cal1 = calo1Jets->begin(); 
02377            cal1 != calo1Jets->end() && jetInd<maxJets; ++cal1 ) { 
02378         
02379         double delR = deltaR( cal1->eta(), cal1->phi(), p4par[i].eta(), p4par[i].phi() ); 
02380 
02381         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02382           dRmin[jetInd] = delR;        // delta R of match
02383           p4cal[jetInd] = cal1->p4();  // Matched Cal 4-vector
02384           usedInd       = nj;
02385           genInd        = i;
02386         }
02387 
02388         /****   
02389         cout << "Delta R = " << delR
02390              << " deltaR min = " << dRmin[jetInd]
02391              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02392              << " p4par = " << p4par[i]
02393              << " p4cal = " << cal1->p4()
02394              << " cal1 eta  = " << cal1->eta()
02395              << " p4par eta = " << p4par[i].eta()
02396              << endl;
02397         ****/
02398 
02399         
02400         nj++;
02401       }
02402 
02403 
02404       // Found matched jet
02405       if (usedInd != -1) {
02406         used[usedInd] = 1;
02407             
02408         dRPar1.Fill(dRmin[jetInd]);
02409         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02410         dPhiPar1.Fill(dphi);
02411         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02412         dEtaPar1.Fill(deta);
02413         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02414         dPtPar1.Fill(dpt);
02415         jetInd++;    
02416       }
02417 
02418     }
02419     ParMatch1.Fill(jetInd);
02420     if (jetInd == 6) {
02421       topMass1.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02422       topMass1.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02423     }
02424 
02425     /***
02426     cout << "Collection Size = " <<  calo1Jets->size() 
02427          << " / " << jetInd
02428          << endl;
02429     ***/
02430 
02431     // ***********************
02432     jetInd = 0;
02433     usedInd = -1;
02434 
02435     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02436     for( int i=0; i<6; ++i ) { 
02437       
02438       dRmin[jetInd]  = 1000.0;
02439 
02440       int nj = 0;      
02441       for( CaloJetCollection::const_iterator cal2 = calo2Jets->begin(); 
02442            cal2 != calo2Jets->end() && jetInd<maxJets; ++cal2 ) { 
02443         
02444         double delR = deltaR( cal2->eta(), cal2->phi(), p4par[i].eta(), p4par[i].phi() ); 
02445 
02446         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) ) {
02447           dRmin[jetInd] = delR;        // delta R of match
02448           p4cal[jetInd] = cal2->p4();  // Matched Cal 4-vector
02449           usedInd       = nj;
02450           genInd        = i;
02451         }
02452 
02453         /****
02454         cout << "Delta R = " << delR
02455              << " deltaR min = " << dRmin[jetInd]
02456              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02457              << " cal2 eta  = " << cal2->eta()
02458              << " p4par eta = " << p4par[i].eta()
02459              << endl;
02460         ****/
02461         
02462         nj++;   
02463       }
02464       if (usedInd != -1) {
02465         used[usedInd] = 1;
02466 
02467         dRPar2.Fill(dRmin[jetInd]);
02468         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02469         dPhiPar2.Fill(dphi);
02470         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02471         dEtaPar2.Fill(deta);
02472         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02473         dPtPar2.Fill(dpt);
02474         
02475         jetInd++;    
02476       }
02477 
02478     }
02479     ParMatch2.Fill(jetInd);
02480     if (jetInd == 6) {
02481       topMass2.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02482       topMass2.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02483     }
02484 
02485 
02486     /***
02487     cout << "Collection Size = " <<  calo2Jets->size() 
02488          << " / " << jetInd
02489          << endl;
02490     ***/
02491 
02492 
02493     // ***********************
02494     jetInd = 0;
02495     usedInd = -1;
02496 
02497     for( int i=0; i<maxJets; ++i ) used[i] = 0;
02498     for( int i=0; i<6; ++i ) { 
02499       
02500       dRmin[jetInd]  = 1000.0;
02501 
02502       int nj = 0;
02503       for( CaloJetCollection::const_iterator cal3 = calo3Jets->begin(); 
02504            cal3 != calo3Jets->end() && jetInd<maxJets; ++cal3 ) { 
02505         
02506         double delR = deltaR( cal3->eta(), cal3->phi(), p4par[i].eta(), p4par[i].phi() ); 
02507 
02508         if ( (delR < dRmin[jetInd]) && (delR < matchedDelR) && (used[nj] == 0) )  {
02509           dRmin[jetInd] = delR;        // delta R of match
02510           p4cal[jetInd] = cal3->p4();  // Matched Cal 4-vector
02511           usedInd       = nj;
02512           genInd        = i;
02513         }
02514         /****
02515         cout << "Delta R = " << delR
02516              << " deltaR min = " << dRmin[jetInd]
02517              << " Ind = " << jetInd << " / " << nj << " / " << used[nj]
02518              << " cal3 eta  = " << cal3->eta()
02519              << " p4par eta = " << p4par[i].eta()
02520              << endl;
02521         ****/
02522         
02523         nj++;
02524       }
02525       if (usedInd != -1) {
02526         used[usedInd] = 1;
02527 
02528         dRPar3.Fill(dRmin[jetInd]);
02529         double dphi = deltaPhi(p4cal[jetInd].phi(), p4par[genInd].phi());
02530         dPhiPar3.Fill(dphi);
02531         double deta = p4cal[jetInd].eta() - p4par[genInd].eta();
02532         dEtaPar3.Fill(deta);
02533         double dpt = p4cal[jetInd].pt() - p4par[genInd].pt();
02534         dPtPar3.Fill(dpt);
02535         
02536         jetInd++;    
02537       }
02538     }
02539     ParMatch3.Fill(jetInd);
02540     if (jetInd == 6) {
02541       topMass3.Fill( (p4cal[0]+p4cal[1]+p4cal[2]).mass() );
02542       topMass3.Fill( (p4cal[3]+p4cal[4]+p4cal[5]).mass() );
02543     }
02544 
02545     /***
02546     cout << "Collection Size = " <<  calo3Jets->size() 
02547          << " / " << jetInd
02548          << endl;
02549     ***/
02550 
02551   }
02552 
02553   int nTow1, nTow2, nTow3, nTow4;
02554 
02555   Handle<CaloJetCollection> jets;
02556 
02557   // *********************
02558   // Jet Properties
02559   // *********************
02560 
02561 
02562 
02563   // --- Loop over jets and make a list of all the used towers
02564   evt.getByLabel( CaloJetAlgorithm1, jets );
02565   int jjet = 0;
02566   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02567     jjet++;
02568 
02569     float hadEne  = ijet->hadEnergyInHB() + ijet->hadEnergyInHO() + 
02570                     ijet->hadEnergyInHE() + ijet->hadEnergyInHF();                   
02571     float emEne   = ijet->emEnergyInEB() + ijet->emEnergyInEE() + ijet->emEnergyInHF();
02572     float had     = ijet->energyFractionHadronic();    
02573 
02574     float j_et = ijet->et();
02575 
02576     if (fabs(ijet->eta()) < 1.3) {
02577       totEneLeadJetEta1_1.Fill(hadEne+emEne); 
02578       hadEneLeadJetEta1_1.Fill(hadEne); 
02579       emEneLeadJetEta1_1.Fill(emEne);       
02580 
02581       totEneLeadJetEta1_2.Fill(hadEne+emEne); 
02582       hadEneLeadJetEta1_2.Fill(hadEne); 
02583       emEneLeadJetEta1_2.Fill(emEne);       
02584 
02585       if (ijet->pt() > minJetPt10) 
02586         hadFracEta11.Fill(had);
02587     }
02588     if ((fabs(ijet->eta()) > 1.3) && (fabs(ijet->eta()) < 3.) ) {
02589 
02590       totEneLeadJetEta2_1.Fill(hadEne+emEne); 
02591       hadEneLeadJetEta2_1.Fill(hadEne); 
02592       emEneLeadJetEta2_1.Fill(emEne);   
02593     
02594       totEneLeadJetEta2_2.Fill(hadEne+emEne); 
02595       hadEneLeadJetEta2_2.Fill(hadEne); 
02596       emEneLeadJetEta2_2.Fill(emEne);       
02597 
02598       if (ijet->pt() > minJetPt10) 
02599         hadFracEta21.Fill(had);
02600     }
02601     if (fabs(ijet->eta()) > 3.) {
02602 
02603       totEneLeadJetEta3_1.Fill(hadEne+emEne); 
02604       hadEneLeadJetEta3_1.Fill(hadEne); 
02605       emEneLeadJetEta3_1.Fill(emEne); 
02606 
02607       totEneLeadJetEta3_2.Fill(hadEne+emEne); 
02608       hadEneLeadJetEta3_2.Fill(hadEne); 
02609       emEneLeadJetEta3_2.Fill(emEne); 
02610 
02611       if (ijet->pt() > minJetPt10) 
02612         hadFracEta31.Fill(had);
02613     }
02614 
02615     
02616     if (jjet == 1) {
02617       hadFracLeadJet1.Fill(had);
02618       hadEneLeadJet1.Fill(hadEne);
02619       hadEneLeadJet12.Fill(hadEne);
02620       hadEneLeadJet13.Fill(hadEne);
02621       emEneLeadJet1.Fill(emEne);
02622       emEneLeadJet12.Fill(emEne);
02623       emEneLeadJet13.Fill(emEne);
02624     }
02625 
02626     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02627     int nConstituents = jetCaloRefs.size();
02628 
02629     if (jjet == 1) {
02630 
02631       nTow1 = nTow2 = nTow3 = nTow4 = 0;
02632       for (int i = 0; i <nConstituents ; i++){
02633 
02634         float et  = jetCaloRefs[i]->et();
02635 
02636         if (et > 0.5) nTow1++;
02637         if (et > 1.0) nTow2++;
02638         if (et > 1.5) nTow3++;
02639         if (et > 2.0) nTow4++;
02640         
02641         hf_TowerJetEt1.Fill(et/j_et);
02642 
02643       }
02644 
02645       nTowersLeadJetPt1.Fill(nTow1);
02646       nTowersLeadJetPt2.Fill(nTow2);
02647       nTowersLeadJetPt3.Fill(nTow3);
02648       nTowersLeadJetPt4.Fill(nTow4);
02649 
02650     }
02651 
02652     if ( (jjet == 1) && (fabs(ijet->eta()) < 1.3) ) {
02653       nTowersLeadJet1.Fill( nConstituents );    
02654       
02655       for (int i = 0; i <nConstituents ; i++){
02656         float t_et  = jetCaloRefs[i]->et();
02657         double delR = deltaR( ijet->eta(), ijet->phi(), jetCaloRefs[i]->eta(), jetCaloRefs[i]->phi() );
02658         hf_TowerDelR1.Fill( delR,  t_et/j_et);
02659         hf_TowerDelR12.Fill( delR,  t_et/j_et);
02660         TowerEtLeadJet1.Fill( t_et );
02661         TowerEtLeadJet12.Fill( t_et );
02662         TowerEtLeadJet13.Fill( t_et );
02663       }
02664     }
02665 
02666 
02667     if ( (jjet == 2) && (fabs(ijet->eta()) < 1.3) ) {
02668       nTowersSecondJet1.Fill( nConstituents );    
02669     }
02670 
02671 
02672 
02673   }
02674 
02675   // *********************
02676   // Unclustered Energy
02677   // *********************
02678 
02679   double SumPtJet(0);
02680 
02681   double SumEtNotJets(0);
02682   double SumEtJets(0);
02683   double SumEtTowers(0);
02684 
02685   double sumJetPx(0);
02686   double sumJetPy(0);  
02687 
02688   double sumTowerAllPx(0);
02689   double sumTowerAllPy(0);
02690 
02691   double sumTowerAllEx(0);
02692   double sumTowerAllEy(0);
02693 
02694   std::vector<CaloTowerPtr>   UsedTowerList;
02695   std::vector<CaloTower>      TowerUsedInJets;
02696   std::vector<CaloTower>      TowerNotUsedInJets;
02697 
02698 
02699   // *********************
02700   // Towers
02701 
02702   Handle<CaloTowerCollection> caloTowers;
02703   evt.getByLabel( "towerMaker", caloTowers );
02704 
02705 
02706 
02707   nTow1 = nTow2 = nTow3 = nTow4 = 0;
02708 
02709   double sum_et = 0.0;
02710   double sum_ex = 0.0;
02711   double sum_ey = 0.0;
02712   //  double sum_ez = 0.0;
02713 
02714   // --- Loop over towers and make a lists of used and unused towers
02715   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02716        tower != caloTowers->end(); tower++) {
02717     
02718     Double_t  et = tower->et();
02719     
02720     if (et > 0.5) nTow1++;
02721     if (et > 1.0) nTow2++;
02722     if (et > 1.5) nTow3++;
02723     if (et > 2.0) nTow4++;
02724 
02725     if(et>0.5) {
02726 
02727       // ********
02728       double phix   = tower->phi();
02729       //      double theta = tower->theta();
02730       //      double e     = tower->energy();
02731       //      double et    = e*sin(theta);
02732       //      double et    = tower->emEt() + tower->hadEt();
02733       double et    = tower->et();
02734 
02735       //      sum_ez += e*cos(theta);
02736       sum_et += et;
02737       sum_ex += et*cos(phix);
02738       sum_ey += et*sin(phix);
02739       // ********
02740 
02741       Double_t phi = tower->phi();
02742       SumEtTowers += tower->et();
02743       
02744       sumTowerAllEx += et*cos(phi);
02745       sumTowerAllEy += et*sin(phi);
02746 
02747     }
02748         
02749   }
02750 
02751   SumEt1.Fill(sum_et);
02752   SumEt12.Fill(sum_et);
02753   SumEt13.Fill(sum_et);
02754 
02755   MET1.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02756   MET12.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02757   MET13.Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02758 
02759   //  met->mex   = -sum_ex;
02760   //  met->mey   = -sum_ey;
02761   //  met->mez   = -sum_ez;
02762   //  met->met   = sqrt( sum_ex*sum_ex + sum_ey*sum_ey );
02763   // cout << "MET = " << met->met << endl;
02764   //  met->sumet = sum_et;
02765   //  met->phi   = atan2( -sum_ey, -sum_ex ); 
02766 
02767   
02768 
02769   hf_sumTowerAllEx.Fill(sumTowerAllEx);
02770   hf_sumTowerAllEy.Fill(sumTowerAllEy);
02771 
02772   nTowers1.Fill(nTow1);
02773   nTowers2.Fill(nTow2);
02774   nTowers3.Fill(nTow3);
02775   nTowers4.Fill(nTow4);
02776 
02777   // *********************
02778   // MidPoint 
02779   //
02780   UsedTowerList.clear();
02781   TowerUsedInJets.clear();
02782   TowerNotUsedInJets.clear();
02783 
02784   // --- Loop over jets and make a list of all the used towers
02785   evt.getByLabel( CaloJetAlgorithm1, jets );
02786   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02787 
02788     Double_t jetPt  = ijet->pt();
02789     Double_t jetPhi = ijet->phi();
02790     
02791     //    if (jetPt>5.0) {
02792 
02793       Double_t jetPx = jetPt*cos(jetPhi);
02794       Double_t jetPy = jetPt*sin(jetPhi);
02795       
02796       sumJetPx +=jetPx;
02797       sumJetPy +=jetPy;
02798 
02799       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02800       int nConstituents = jetCaloRefs.size();
02801       for (int i = 0; i <nConstituents ; i++){
02802         UsedTowerList.push_back(jetCaloRefs[i]);
02803       }
02804 
02805       SumPtJet +=jetPt;
02806     //    }
02807 
02808   }
02809 
02810 
02811   int NTowersUsed = UsedTowerList.size();
02812       
02813   // --- Loop over towers and make a lists of used and unused towers
02814   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02815        tower != caloTowers->end(); tower++) {
02816     
02817     CaloTower  t = *tower;
02818     Double_t  et = tower->et();
02819 
02820     if(et>0) {
02821 
02822       Double_t phi = tower->phi();
02823       SumEtTowers += tower->et();
02824       
02825       sumTowerAllPx += et*cos(phi);
02826       sumTowerAllPy += et*sin(phi);
02827 
02828       bool used = false;
02829 
02830       for(int i=0; i<NTowersUsed; i++){
02831         if(tower->id() == UsedTowerList[i]->id()){
02832           used=true;
02833           break;
02834         }
02835       }
02836       
02837       if (used) {
02838         TowerUsedInJets.push_back(t);
02839       } else {
02840         TowerNotUsedInJets.push_back(t);
02841       }
02842 
02843     }
02844 
02845   }
02846 
02847   int nUsed    = TowerUsedInJets.size();
02848   int nNotUsed = TowerNotUsedInJets.size();
02849   
02850   SumEtJets    = 0;
02851   SumEtNotJets = 0;
02852 
02853   for(int i=0;i<nUsed;i++){
02854     SumEtJets += TowerUsedInJets[i].et();
02855   }
02856   h_jetEt1.Fill(SumEtJets);
02857 
02858   for(int i=0;i<nNotUsed;i++){
02859     if (TowerNotUsedInJets[i].et() > 0.5) 
02860       SumEtNotJets += TowerNotUsedInJets[i].et();
02861     h_missEt1.Fill(TowerNotUsedInJets[i].et());
02862     h_missEt1s.Fill(TowerNotUsedInJets[i].et());
02863   }
02864   h_totMissEt1.Fill(SumEtNotJets);
02865 
02866 
02867 
02868   // *********************
02869   // SISCone
02870   //
02871   UsedTowerList.clear();
02872   TowerUsedInJets.clear();
02873   TowerNotUsedInJets.clear();
02874 
02875   // --- Loop over jets and make a list of all the used towers
02876   evt.getByLabel( CaloJetAlgorithm2, jets );
02877   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02878 
02879     Double_t jetPt  = ijet->pt();
02880     Double_t jetPhi = ijet->phi();
02881     
02882     //    if (jetPt>5.0) {
02883 
02884       Double_t jetPx = jetPt*cos(jetPhi);
02885       Double_t jetPy = jetPt*sin(jetPhi);
02886       
02887       sumJetPx +=jetPx;
02888       sumJetPy +=jetPy;
02889 
02890       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02891       int nConstituents = jetCaloRefs.size();
02892       for (int i = 0; i <nConstituents ; i++){
02893         UsedTowerList.push_back(jetCaloRefs[i]);
02894       }
02895 
02896       SumPtJet +=jetPt;
02897     //    }
02898 
02899   }
02900 
02901 
02902   //  Handle<CaloTowerCollection> caloTowers;
02903   //  evt.getByLabel( "towerMaker", caloTowers );
02904 
02905   NTowersUsed = UsedTowerList.size();
02906       
02907   // --- Loop over towers and make a lists of used and unused towers
02908   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02909        tower != caloTowers->end(); tower++) {
02910     
02911     CaloTower  t = *tower;
02912     Double_t  et = tower->et();
02913 
02914     if(et>0) {
02915 
02916       Double_t phi = tower->phi();
02917 
02918       SumEtTowers += tower->et();
02919       
02920       sumTowerAllPx += et*cos(phi);
02921       sumTowerAllPy += et*sin(phi);
02922 
02923       bool used = false;
02924 
02925       for(int i=0; i<NTowersUsed; i++){
02926         if(tower->id() == UsedTowerList[i]->id()){
02927           used=true;
02928           break;
02929         }
02930       }
02931       
02932       if (used) {
02933         TowerUsedInJets.push_back(t);
02934       } else {
02935         TowerNotUsedInJets.push_back(t);
02936       }
02937 
02938     }
02939 
02940   }
02941 
02942   nUsed    = TowerUsedInJets.size();
02943   nNotUsed = TowerNotUsedInJets.size();
02944   
02945   SumEtJets    = 0;
02946   SumEtNotJets = 0;
02947 
02948   for(int i=0;i<nUsed;i++){
02949     SumEtJets += TowerUsedInJets[i].et();
02950   }
02951   h_jetEt2.Fill(SumEtJets);
02952 
02953   for(int i=0;i<nNotUsed;i++){
02954     if (TowerNotUsedInJets[i].et() > 0.5) 
02955       SumEtNotJets += TowerNotUsedInJets[i].et();
02956     h_missEt2.Fill(TowerNotUsedInJets[i].et());
02957     h_missEt2s.Fill(TowerNotUsedInJets[i].et());
02958   }
02959   h_totMissEt2.Fill(SumEtNotJets);
02960 
02961 
02962   // *********************
02963   // KtClus
02964   //
02965   UsedTowerList.clear();
02966   TowerUsedInJets.clear();
02967   TowerNotUsedInJets.clear();
02968 
02969   // --- Loop over jets and make a list of all the used towers
02970   evt.getByLabel( CaloJetAlgorithm3, jets );
02971   for ( CaloJetCollection::const_iterator ijet=jets->begin(); ijet!=jets->end(); ijet++) {
02972 
02973     Double_t jetPt  = ijet->pt();
02974     Double_t jetPhi = ijet->phi();
02975     
02976     //    if (jetPt>5.0) {
02977 
02978       Double_t jetPx = jetPt*cos(jetPhi);
02979       Double_t jetPy = jetPt*sin(jetPhi);
02980       
02981       sumJetPx +=jetPx;
02982       sumJetPy +=jetPy;
02983 
02984       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02985       int nConstituents = jetCaloRefs.size();
02986       for (int i = 0; i <nConstituents ; i++){
02987         UsedTowerList.push_back(jetCaloRefs[i]);
02988       }
02989 
02990       SumPtJet +=jetPt;
02991     //    }
02992 
02993   }
02994 
02995 
02996   //  Handle<CaloTowerCollection> caloTowers;
02997   //  evt.getByLabel( "towerMaker", caloTowers );
02998 
02999   NTowersUsed = UsedTowerList.size();
03000       
03001   // --- Loop over towers and make a lists of used and unused towers
03002   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
03003        tower != caloTowers->end(); tower++) {
03004     
03005     CaloTower  t = *tower;
03006     Double_t  et = tower->et();
03007 
03008     if(et>0) {
03009 
03010       //      Double_t phi = tower->phi();
03011 
03012       //      SumEtTowers   += tower->et();      
03013       //      sumTowerAllPx += et*cos(phi);
03014       //      sumTowerAllPy += et*sin(phi);
03015 
03016       bool used = false;
03017 
03018       for(int i=0; i<NTowersUsed; i++){
03019         if(tower->id() == UsedTowerList[i]->id()){
03020           used=true;
03021           break;
03022         }
03023       }
03024       
03025       if (used) {
03026         TowerUsedInJets.push_back(t);
03027       } else {
03028         TowerNotUsedInJets.push_back(t);
03029       }
03030 
03031     }
03032 
03033   }
03034 
03035   nUsed    = TowerUsedInJets.size();
03036   nNotUsed = TowerNotUsedInJets.size();
03037   
03038   SumEtJets    = 0;
03039   SumEtNotJets = 0;
03040 
03041   for(int i=0;i<nUsed;i++){
03042     SumEtJets += TowerUsedInJets[i].et();
03043   }
03044   h_jetEt3.Fill(SumEtJets);
03045 
03046   for(int i=0;i<nNotUsed;i++){
03047     if (TowerNotUsedInJets[i].et() > 0.5) 
03048       SumEtNotJets += TowerNotUsedInJets[i].et();
03049     h_missEt3.Fill(TowerNotUsedInJets[i].et());
03050     h_missEt3s.Fill(TowerNotUsedInJets[i].et());
03051   }
03052   h_totMissEt3.Fill(SumEtNotJets);
03053 
03054 }
03055 
03056 
03057 
03058 void myFastSimVal::endJob() {
03059 
03060   //Write out the histogram file.
03061   m_file->Write(); 
03062 
03063 }
03064 #include "FWCore/Framework/interface/MakerMacros.h"
03065 DEFINE_FWK_MODULE(myFastSimVal);