CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoJets/JetAnalyzers/src/myJetAna.cc

Go to the documentation of this file.
00001 // myJetAna.cc
00002 // Description:  Access Cruzet Data
00003 // Author: Frank Chlebana
00004 // Date:  24 - July - 2008
00005 // 
00006 #include "RecoJets/JetAnalyzers/interface/myJetAna.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/METReco/interface/CaloMETCollection.h"
00012 #include "DataFormats/METReco/interface/CaloMET.h"
00013 
00014 #include "DataFormats/Common/interface/Handle.h"
00015 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"
00016 
00017 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00018 #include "DataFormats/HcalRecHit/interface/HcalSourcePositionData.h"
00019 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00020 
00021 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00022 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00023 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00024 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00025 
00026 #include "FWCore/Framework/interface/ESHandle.h"
00027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00028 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00029 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00030 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00031 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00032 
00033 #include "DataFormats/TrackReco/interface/Track.h"
00034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00035 #include "DataFormats/TrackReco/interface/Track.h"
00036 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00037 
00038 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
00039 
00040 #include "DataFormats/MuonReco/interface/Muon.h"
00041 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00042 
00043 #include "DataFormats/VertexReco/interface/Vertex.h"
00044 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00045 
00046 // #include "DataFormats/PhotonReco/interface/PhotonFwd.h"
00047 // #include "DataFormats/PhotonReco/interface/Photon.h"
00048 
00049 
00050 #include "DataFormats/Math/interface/deltaR.h"
00051 #include "DataFormats/Math/interface/deltaPhi.h"
00052 // #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
00053 #include "DataFormats/Candidate/interface/Candidate.h"
00054 // #include "FWCore/Framework/interface/Handle.h"
00055 #include "FWCore/Framework/interface/Event.h"
00056 // #include "FWCore/Framework/interface/EventSetup.h"
00057 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00058 
00059 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00060 
00061 #include "DataFormats/HLTReco/interface/TriggerObject.h"
00062 // #include "FWCore/Common/interface/TriggerNames.h"
00063 #include "DataFormats/Common/interface/TriggerResults.h"
00064 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
00065 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
00066 
00067 // #include "DataFormats/Scalers/interface/DcsStatus.h"
00068 
00069 // include files
00070 #include "CommonTools/RecoAlgos/interface/HBHENoiseFilter.h"
00071 #include "DataFormats/METReco/interface/HcalNoiseSummary.h"
00072 
00073 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCaloFlagLabels.h"  
00074 
00075 
00076 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00077 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
00078 #include <TROOT.h>
00079 #include <TSystem.h>
00080 #include <TFile.h>
00081 #include <TCanvas.h>
00082 #include <cmath>
00083 
00084 using namespace edm;
00085 using namespace reco;
00086 using namespace std;
00087 
00088 
00089 #define INVALID 9999.
00090 #define DEBUG false
00091 #define MAXJETS 100
00092 
00093 typedef struct RBX_struct {
00094   double et;
00095   double hadEnergy;
00096   double emEnergy;
00097   float  hcalTime;
00098   float  ecalTime;
00099   int    nTowers;
00100 } RBX ;
00101 
00102 typedef struct HPD_struct {
00103   double et;
00104   double hadEnergy;
00105   double emEnergy;
00106   double time;
00107   float  hcalTime;
00108   float  ecalTime;
00109   int    nTowers;
00110 } HPD ;
00111 
00112 
00113 
00114 float totBNC, nBNC[4000];
00115 
00116 RBX RBXColl[36];
00117 HPD HPDColl[144];
00118 
00119 // ************************
00120 // ************************
00121 
00122 // Get the algorithm of the jet collections we will read from the .cfg file 
00123 // which defines the value of the strings CaloJetAlgorithm and GenJetAlgorithm.
00124 
00125 myJetAna::myJetAna( const ParameterSet & cfg ) :
00126   CaloJetAlgorithm( cfg.getParameter<string>( "CaloJetAlgorithm" ) ), 
00127   GenJetAlgorithm( cfg.getParameter<string>( "GenJetAlgorithm" ) )  
00128 {
00129   theTriggerResultsLabel = cfg.getParameter<edm::InputTag>("TriggerResultsLabel");
00130 }
00131 
00132 
00133 // ************************
00134 // ************************
00135 
00136 void myJetAna::beginJob( ) {
00137 
00138 
00139 
00140   edm::Service<TFileService> fs;
00141 
00142   // --- passed selection cuts 
00143   h_pt     = fs->make<TH1F>( "pt",  "Jet p_{T}", 100, 0, 50 );
00144   h_ptRBX  = fs->make<TH1F>( "ptRBX",  "RBX: Jet p_{T}", 100, 0, 50 );
00145   h_ptHPD  = fs->make<TH1F>( "ptHPD",  "HPD: Jet p_{T}", 100, 0, 50 );
00146   h_ptTower  = fs->make<TH1F>( "ptTower",  "Jet p_{T}", 100, 0, 50 );
00147   h_et     = fs->make<TH1F>( "et",  "Jet E_{T}", 100, 0, 50 );
00148   h_eta    = fs->make<TH1F>( "eta", "Jet #eta", 100, -5, 5 );
00149   h_phi    = fs->make<TH1F>( "phi", "Jet #phi", 50, -M_PI, M_PI );
00150   // ---
00151 
00152   hitEtaEt  = fs->make<TH1F>( "hitEtaEt", "RecHit #eta", 90, -45, 45 );
00153   hitEta    = fs->make<TH1F>( "hitEta", "RecHit #eta", 90, -45, 45 );
00154   hitPhi    = fs->make<TH1F>( "hitPhi", "RecHit #phi", 73, 0, 73 );
00155 
00156   caloEtaEt  = fs->make<TH1F>( "caloEtaEt", "CaloTower #eta", 100, -4, 4 );
00157   caloEta    = fs->make<TH1F>( "caloEta", "CaloTower #eta", 100, -4, 4 );
00158   caloPhi    = fs->make<TH1F>( "caloPhi", "CaloTower #phi", 50, -M_PI, M_PI );
00159 
00160   dijetMass  =  fs->make<TH1F>("dijetMass","DiJet Mass",100,0,100);
00161 
00162   totEneLeadJetEta1 = fs->make<TH1F>("totEneLeadJetEta1","Total Energy Lead Jet Eta1 1",100,0,100);
00163   totEneLeadJetEta2 = fs->make<TH1F>("totEneLeadJetEta2","Total Energy Lead Jet Eta2 1",150,0,150);
00164   totEneLeadJetEta3 = fs->make<TH1F>("totEneLeadJetEta3","Total Energy Lead Jet Eta3 1",150,0,150);
00165 
00166   hadEneLeadJetEta1 = fs->make<TH1F>("hadEneLeadJetEta1","Hadronic Energy Lead Jet Eta1 1",50,0,50);
00167   hadEneLeadJetEta2 = fs->make<TH1F>("hadEneLeadJetEta2","Hadronic Energy Lead Jet Eta2 1",100,0,100);
00168   hadEneLeadJetEta3 = fs->make<TH1F>("hadEneLeadJetEta3","Hadronic Energy Lead Jet Eta3 1",100,0,100);
00169   emEneLeadJetEta1  = fs->make<TH1F>("emEneLeadJetEta1","EM Energy Lead Jet Eta1 1",50,0,50);
00170   emEneLeadJetEta2  = fs->make<TH1F>("emEneLeadJetEta2","EM Energy Lead Jet Eta2 1",100,0,100);
00171   emEneLeadJetEta3  = fs->make<TH1F>("emEneLeadJetEta3","EM Energy Lead Jet Eta3 1",100,0,100);
00172 
00173 
00174   hadFracEta1 = fs->make<TH1F>("hadFracEta11","Hadronic Fraction Eta1 Jet 1",100,0,1);
00175   hadFracEta2 = fs->make<TH1F>("hadFracEta21","Hadronic Fraction Eta2 Jet 1",100,0,1);
00176   hadFracEta3 = fs->make<TH1F>("hadFracEta31","Hadronic Fraction Eta3 Jet 1",100,0,1);
00177 
00178   HFSumEt  = fs->make<TH1F>("HFSumEt","HFSumEt",100,0,100);
00179   HFMET    = fs->make<TH1F>("HFMET",  "HFMET",120,0,120);
00180 
00181   SumEt  = fs->make<TH1F>("SumEt","SumEt",100,0,100);
00182   MET    = fs->make<TH1F>("MET",  "MET",120,0,120);
00183   OERMET    = fs->make<TH1F>("OERMET",  "OERMET",120,0,120);
00184   METSig = fs->make<TH1F>("METSig",  "METSig",100,0,50);
00185   MEx    = fs->make<TH1F>("MEx",  "MEx",100,-20,20);
00186   MEy    = fs->make<TH1F>("MEy",  "MEy",100,-20,20);
00187   METPhi = fs->make<TH1F>("METPhi",  "METPhi",315,0,3.15);
00188   MET_RBX    = fs->make<TH1F>("MET_RBX",  "MET",100,0,1000);
00189   MET_HPD    = fs->make<TH1F>("MET_HPD",  "MET",100,0,1000);
00190   MET_Tower  = fs->make<TH1F>("MET_Tower",  "MET",100,0,1000);
00191 
00192   SiClusters = fs->make<TH1F>("SiClusters",  "SiClusters",150,0,1500);
00193 
00194   h_Vx     = fs->make<TH1F>("Vx",  "Vx",100,-0.5,0.5);
00195   h_Vy     = fs->make<TH1F>("Vy",  "Vy",100,-0.5,0.5);
00196   h_Vz     = fs->make<TH1F>("Vz",  "Vz",100,-20,20);
00197   h_VNTrks = fs->make<TH1F>("VNTrks",  "VNTrks",10,1,100);
00198 
00199   h_Trk_pt   = fs->make<TH1F>("Trk_pt",  "Trk_pt",100,0,20);
00200   h_Trk_NTrk = fs->make<TH1F>("Trk_NTrk",  "Trk_NTrk",150,0,150);
00201 
00202   hf_sumTowerAllEx = fs->make<TH1F>("sumTowerAllEx","Tower Ex",100,-1000,1000);
00203   hf_sumTowerAllEy = fs->make<TH1F>("sumTowerAllEy","Tower Ey",100,-1000,1000);
00204 
00205   hf_TowerJetEt   = fs->make<TH1F>("TowerJetEt","Tower/Jet Et 1",50,0,1);
00206 
00207   ETime = fs->make<TH1F>("ETime","Ecal Time",200,-200,200);
00208   HTime = fs->make<TH1F>("HTime","Hcal Time",200,-200,200);
00209 
00210   towerHadEnHB    = fs->make<TH1F>("towerHadEnHB" ,"HB: Calo Tower HAD Energy",210,-1,20);
00211   towerHadEnHE    = fs->make<TH1F>("towerHadEnHE" ,"HE: Calo Tower HAD Energy",510,-1,50);
00212   towerHadEnHF    = fs->make<TH1F>("towerHadEnHF" ,"HF: Calo Tower HAD Energy",510,-1,50);
00213 
00214   towerEmEnHB    = fs->make<TH1F>("towerEmEnHB" ,"HB: Calo Tower EM Energy",210,-1,20);
00215   towerEmEnHE    = fs->make<TH1F>("towerEmEnHE" ,"HE: Calo Tower EM Energy",510,-1,50);
00216   towerEmEnHF    = fs->make<TH1F>("towerEmEnHF" ,"HF: Calo Tower EM Energy",510,-1,50);
00217 
00218   towerHadEn    = fs->make<TH1F>("towerHadEn" ,"Hadronic Energy in Calo Tower",2000,-100,100);
00219   towerEmEn     = fs->make<TH1F>("towerEmEn"  ,"EM Energy in Calo Tower",2000,-100,100);
00220   towerOuterEn  = fs->make<TH1F>("towerOuterEn"  ,"HO Energy in Calo Tower",2000,-100,100);
00221 
00222   towerEmFrac   = fs->make<TH1F>("towerEmFrac","EM Fraction of Energy in Calo Tower",100,-1.,1.);
00223 
00224   RBX_et        = fs->make<TH1F>("RBX_et","ET in RBX",1000,-20,100);
00225   RBX_hadEnergy = fs->make<TH1F>("RBX_hadEnergy","Hcal Energy in RBX",1000,-20,100);
00226   RBX_hcalTime  = fs->make<TH1F>("RBX_hcalTime","Hcal Time in RBX",200,-200,200);
00227   RBX_nTowers   = fs->make<TH1F>("RBX_nTowers","Number of Towers in RBX",75,0,75);
00228   RBX_N         = fs->make<TH1F>("RBX_N","Number of RBX",10,0,10);
00229 
00230   HPD_et        = fs->make<TH1F>("HPD_et","ET in HPD",1000,-20,100);
00231   HPD_hadEnergy = fs->make<TH1F>("HPD_hadEnergy","Hcal Energy in HPD",1000,-20,100);
00232   HPD_hcalTime  = fs->make<TH1F>("HPD_hcalTime","Hcal Time in HPD",200,-200,200);
00233   HPD_nTowers   = fs->make<TH1F>("HPD_nTowers","Number of Towers in HPD",20,0,20);
00234   HPD_N         = fs->make<TH1F>("HPD_N","Number of HPD",10,0,10);
00235   
00236   nTowers1  = fs->make<TH1F>("nTowers1","Number of Towers pt 0.5",100,0,200);
00237   nTowers2  = fs->make<TH1F>("nTowers2","Number of Towers pt 1.0",100,0,200);
00238   nTowers3  = fs->make<TH1F>("nTowers3","Number of Towers pt 1.5",100,0,200);
00239   nTowers4  = fs->make<TH1F>("nTowers4","Number of Towers pt 2.0",100,0,200);
00240 
00241   nTowersLeadJetPt1  = fs->make<TH1F>("nTowersLeadJetPt1","Number of Towers in Lead Jet pt 0.5",100,0,100);
00242   nTowersLeadJetPt2  = fs->make<TH1F>("nTowersLeadJetPt2","Number of Towers in Lead Jet pt 1.0",100,0,100);
00243   nTowersLeadJetPt3  = fs->make<TH1F>("nTowersLeadJetPt3","Number of Towers in Lead Jet pt 1.5",100,0,100);
00244   nTowersLeadJetPt4  = fs->make<TH1F>("nTowersLeadJetPt4","Number of Towers in Lead Jet pt 2.0",100,0,100);
00245 
00246   h_nCalJets  =  fs->make<TH1F>( "nCalJets",  "Number of CalJets", 20, 0, 20 );
00247 
00248   HBEneOOT  = fs->make<TH1F>( "HBEneOOT",  "HBEneOOT", 200, -5, 10 );
00249   HEEneOOT  = fs->make<TH1F>( "HEEneOOT",  "HEEneOOT", 200, -5, 10 );
00250   HFEneOOT  = fs->make<TH1F>( "HFEneOOT",  "HFEneOOT", 200, -5, 10 );
00251   HOEneOOT  = fs->make<TH1F>( "HOEneOOT",  "HOEneOOT", 200, -5, 10 );
00252 
00253   HBEneOOTTh  = fs->make<TH1F>( "HBEneOOTTh",  "HBEneOOTTh", 200, -5, 10 );
00254   HEEneOOTTh  = fs->make<TH1F>( "HEEneOOTTh",  "HEEneOOTTh", 200, -5, 10 );
00255   HFEneOOTTh  = fs->make<TH1F>( "HFEneOOTTh",  "HFEneOOTTh", 200, -5, 10 );
00256   HOEneOOTTh  = fs->make<TH1F>( "HOEneOOTTh",  "HOEneOOTTh", 200, -5, 10 );
00257 
00258   HBEneOOTTh1  = fs->make<TH1F>( "HBEneOOTTh1",  "HBEneOOT", 200, -5, 10 );
00259   HEEneOOTTh1  = fs->make<TH1F>( "HEEneOOTTh1",  "HEEneOOT", 200, -5, 10 );
00260   HFEneOOTTh1  = fs->make<TH1F>( "HFEneOOTTh1",  "HFEneOOT", 200, -5, 10 );
00261   HOEneOOTTh1  = fs->make<TH1F>( "HOEneOOTTh1",  "HOEneOOT", 200, -5, 10 );
00262 
00263   HBEneTThr     = fs->make<TH1F>( "HBEneTThr",  "HBEneTThr", 105, -5, 100 );
00264   HEEneTThr     = fs->make<TH1F>( "HEEneTThr",  "HEEneTThr", 105, -5, 100 );
00265   HFEneTThr     = fs->make<TH1F>( "HFEneTThr",  "HFEneTThr", 105, -5, 100 );
00266 
00267 
00268   HBEne     = fs->make<TH1F>( "HBEne",  "HBEne", 205, -5, 200 );
00269   HBEneTh   = fs->make<TH1F>( "HBEneTh",  "HBEneTh", 205, -5, 200 );
00270   HBEneTh1  = fs->make<TH1F>( "HBEneTh1",  "HBEneTh1", 205, -5, 200 );
00271   HBEneX    = fs->make<TH1F>( "HBEneX",  "HBEneX", 200, -5, 10 );
00272   HBEneY    = fs->make<TH1F>( "HBEneY",  "HBEnedY", 200, -5, 10 );
00273   HBTime    = fs->make<TH1F>( "HBTime", "HBTime", 200, -100, 100 );
00274   HBTimeTh  = fs->make<TH1F>( "HBTimeTh", "HBTimeTh", 200, -100, 100 );
00275   HBTimeTh1  = fs->make<TH1F>( "HBTimeTh1", "HBTimeTh1", 200, -100, 100 );
00276   HBTimeTh2  = fs->make<TH1F>( "HBTimeTh2", "HBTimeTh2", 200, -100, 100 );
00277   HBTimeTh3  = fs->make<TH1F>( "HBTimeTh3", "HBTimeTh3", 200, -100, 100 );
00278   HBTimeThR  = fs->make<TH1F>( "HBTimeThR", "HBTimeThR", 200, -100, 100 );
00279   HBTimeTh1R  = fs->make<TH1F>( "HBTimeTh1R", "HBTimeTh1R", 200, -100, 100 );
00280   HBTimeTh2R  = fs->make<TH1F>( "HBTimeTh2R", "HBTimeTh2R", 200, -100, 100 );
00281   HBTimeTh3R  = fs->make<TH1F>( "HBTimeTh3R", "HBTimeTh3R", 200, -100, 100 );
00282 
00283   HBTimeFlagged    = fs->make<TH1F>( "HBTimeFlagged", "HBTimeFlagged", 200, -100, 100 );
00284   HBTimeThFlagged    = fs->make<TH1F>( "HBTimeThFlagged", "HBTimeThFlagged", 200, -100, 100 );
00285   HBTimeTh1Flagged    = fs->make<TH1F>( "HBTimeTh1Flagged", "HBTimeTh1Flagged", 200, -100, 100 );
00286   HBTimeTh2Flagged    = fs->make<TH1F>( "HBTimeTh2Flagged", "HBTimeTh2Flagged", 200, -100, 100 );
00287 
00288   HBTimeFlagged2    = fs->make<TH1F>( "HBTimeFlagged2", "HBTimeFlagged2", 200, -100, 100 );
00289   HBTimeThFlagged2    = fs->make<TH1F>( "HBTimeThFlagged2", "HBTimeThFlagged2", 200, -100, 100 );
00290   HBTimeTh1Flagged2    = fs->make<TH1F>( "HBTimeTh1Flagged2", "HBTimeTh1Flagged2", 200, -100, 100 );
00291   HBTimeTh2Flagged2    = fs->make<TH1F>( "HBTimeTh2Flagged2", "HBTimeTh2Flagged2", 200, -100, 100 );
00292 
00293   HBTimeX   = fs->make<TH1F>( "HBTimeX", "HBTimeX", 200, -100, 100 );
00294   HBTimeY   = fs->make<TH1F>( "HBTimeY", "HBTimeY", 200, -100, 100 );
00295   HEEne     = fs->make<TH1F>( "HEEne",  "HEEne", 205, -5, 200 );
00296   HEEneTh   = fs->make<TH1F>( "HEEneTh",  "HEEneTh", 205, -5, 200 );
00297   HEEneTh1  = fs->make<TH1F>( "HEEneTh1",  "HEEneTh1", 205, -5, 200 );
00298   HEEneX    = fs->make<TH1F>( "HEEneX",  "HEEneX", 200, -5, 10 );
00299   HEEneY    = fs->make<TH1F>( "HEEneY",  "HEEneY", 200, -5, 10 );
00300   HEposEne  = fs->make<TH1F>( "HEposEne",  "HEposEne", 200, -5, 10 );
00301   HEnegEne  = fs->make<TH1F>( "HEnegEne",  "HEnegEne", 200, -5, 10 );
00302   HETime    = fs->make<TH1F>( "HETime", "HETime", 200, -100, 100 );
00303   HETimeTh  = fs->make<TH1F>( "HETimeTh", "HETimeTh", 200, -100, 100 );
00304   HETimeTh1  = fs->make<TH1F>( "HETimeTh1", "HETimeTh1", 200, -100, 100 );
00305   HETimeTh2  = fs->make<TH1F>( "HETimeTh2", "HETimeTh2", 200, -100, 100 );
00306   HETimeTh3  = fs->make<TH1F>( "HETimeTh3", "HETimeTh3", 200, -100, 100 );
00307   HETimeThR  = fs->make<TH1F>( "HETimeThR", "HETimeThR", 200, -100, 100 );
00308   HETimeTh1R  = fs->make<TH1F>( "HETimeTh1R", "HETimeTh1R", 200, -100, 100 );
00309   HETimeTh2R  = fs->make<TH1F>( "HETimeTh2R", "HETimeTh2R", 200, -100, 100 );
00310   HETimeTh3R  = fs->make<TH1F>( "HETimeTh3R", "HETimeTh3R", 200, -100, 100 );
00311 
00312   HETimeFlagged    = fs->make<TH1F>( "HETimeFlagged", "HETimeFlagged", 200, -100, 100 );
00313   HETimeThFlagged    = fs->make<TH1F>( "HETimeThFlagged", "HETimeThFlagged", 200, -100, 100 );
00314   HETimeTh1Flagged    = fs->make<TH1F>( "HETimeTh1Flagged", "HETimeTh1Flagged", 200, -100, 100 );
00315   HETimeTh2Flagged    = fs->make<TH1F>( "HETimeTh2Flagged", "HETimeTh2Flagged", 200, -100, 100 );
00316 
00317   HETimeFlagged2    = fs->make<TH1F>( "HETimeFlagged2", "HETimeFlagged2", 200, -100, 100 );
00318   HETimeThFlagged2    = fs->make<TH1F>( "HETimeThFlagged2", "HETimeThFlagged2", 200, -100, 100 );
00319   HETimeTh1Flagged2    = fs->make<TH1F>( "HETimeTh1Flagged2", "HETimeTh1Flagged2", 200, -100, 100 );
00320   HETimeTh2Flagged2    = fs->make<TH1F>( "HETimeTh2Flagged2", "HETimeTh2Flagged2", 200, -100, 100 );
00321 
00322   HETimeX   = fs->make<TH1F>( "HETimeX", "HETimeX", 200, -100, 100 );
00323   HETimeY   = fs->make<TH1F>( "HETimeY", "HETimeY", 200, -100, 100 );
00324   HEposTime = fs->make<TH1F>( "HEposTime",  "HEposTime", 200, -100, 100 );
00325   HEnegTime = fs->make<TH1F>( "HEnegTime",  "HEnegTime", 200, -100, 100 );
00326   HOEne     = fs->make<TH1F>( "HOEne",  "HOEne", 200, -5, 10 );
00327   HOEneTh   = fs->make<TH1F>( "HOEneTh",  "HOEneTh", 200, -5, 10 );
00328   HOEneTh1  = fs->make<TH1F>( "HOEneTh1",  "HOEneTh1", 200, -5, 10 );
00329   HOTime    = fs->make<TH1F>( "HOTime", "HOTime", 200, -100, 100 );
00330   HOTimeTh  = fs->make<TH1F>( "HOTimeTh", "HOTimeTh", 200, -100, 100 );
00331 
00332   // Histos for separating SiPMs and HPDs in HO:
00333   HOSEne     = fs->make<TH1F>( "HOSEne",  "HOSEne", 12000, -20, 100 );
00334   HOSTime    = fs->make<TH1F>( "HOSTime", "HOSTime", 200, -100, 100 );
00335   HOHEne     = fs->make<TH1F>( "HOHEne",  "HOHEne", 12000, -20, 100 );
00336   HOHTime    = fs->make<TH1F>( "HOHTime", "HOHTime", 200, -100, 100 );
00337 
00338   HOHr0Ene      = fs->make<TH1F>( "HOHr0Ene"  ,   "HOHr0Ene", 12000, -20 , 100 );
00339   HOHr0Time     = fs->make<TH1F>( "HOHr0Time" ,  "HOHr0Time",   200, -200, 200 );
00340   HOHrm1Ene     = fs->make<TH1F>( "HOHrm1Ene" ,  "HOHrm1Ene", 12000, -20 , 100 );
00341   HOHrm1Time    = fs->make<TH1F>( "HOHrm1Time", "HOHrm1Time",   200, -200, 200 );
00342   HOHrm2Ene     = fs->make<TH1F>( "HOHrm2Ene" ,  "HOHrm2Ene", 12000, -20 , 100 );
00343   HOHrm2Time    = fs->make<TH1F>( "HOHrm2Time", "HOHrm2Time",   200, -200, 200 );
00344   HOHrp1Ene     = fs->make<TH1F>( "HOHrp1Ene" ,  "HOHrp1Ene", 12000, -20 , 100 );
00345   HOHrp1Time    = fs->make<TH1F>( "HOHrp1Time", "HOHrp1Time",   200, -200, 200 );
00346   HOHrp2Ene     = fs->make<TH1F>( "HOHrp2Ene" ,  "HOHrp2Ene", 12000, -20 , 100 );
00347   HOHrp2Time    = fs->make<TH1F>( "HOHrp2Time", "HOHrp2Time",   200, -200, 200 );
00348 
00349   HBTvsE    = fs->make<TH2F>( "HBTvsE", "HBTvsE",305, -5, 300, 100, -100, 100);
00350   HETvsE    = fs->make<TH2F>( "HETvsE", "HETvsE",305, -5, 300, 100, -100, 100);
00351 
00352   HFTvsE            = fs->make<TH2F>( "HFTvsE", "HFTvsE",305, -5, 300, 100, -100, 100);
00353   HFTvsEFlagged     = fs->make<TH2F>( "HFTvsEFlagged", "HFTvsEFlagged",305, -5, 300, 100, -100, 100);  
00354   HFTvsEFlagged2    = fs->make<TH2F>( "HFTvsEFlagged2", "HFTvsEFlagged2",305, -5, 300, 100, -100, 100);
00355 
00356   HFTvsEThr    = fs->make<TH2F>( "HFTvsEThr", "HFTvsEThr",305, -5, 300, 100, -100, 100);
00357   HFTvsEFlaggedThr    = fs->make<TH2F>( "HFTvsEFlaggedThr", "HFTvsEFlaggedThr",305, -5, 300, 100, -100, 100);  
00358   HFTvsEFlagged2Thr    = fs->make<TH2F>( "HFTvsEFlagged2Thr", "HFTvsEFlagged2Thr",305, -5, 300, 100, -100, 100);
00359 
00360   HOTvsE    = fs->make<TH2F>( "HOTvsE", "HOTvsE",305, -5, 300, 100, -100, 100);
00361 
00362   HFvsZ    = fs->make<TH2F>( "HFvsZ", "HFvsZ",100,-50,50,100,-50,50);
00363 
00364 
00365 
00366   HOocc    = fs->make<TH2F>( "HOocc", "HOocc",85,-42.5,42.5,70,0.5,70.5);
00367   HBocc    = fs->make<TH2F>( "HBocc", "HBocc",85,-42.5,42.5,70,0.5,70.5);
00368   HEocc    = fs->make<TH2F>( "HEocc", "HEocc",85,-42.5,42.5,70,0.5,70.5);
00369   HFocc    = fs->make<TH2F>( "HFocc", "HFocc",85,-42.5,42.5,70,0.5,70.5);
00370   HFoccTime    = fs->make<TH2F>( "HFoccTime", "HFoccTime",85,-42.5,42.5,70,0.5,70.5);
00371   HFoccFlagged    = fs->make<TH2F>( "HFoccFlagged", "HFoccFlagged",85,-42.5,42.5,70,0.5,70.5);
00372   HFoccFlagged2    = fs->make<TH2F>( "HFoccFlagged2", "HFoccFlagged2",85,-42.5,42.5,70,0.5,70.5);
00373 
00374   HFEtaPhiNFlagged    = fs->make<TH2F>( "HFEtaPhiNFlagged", "HFEtaPhiNFlagged",85,-42.5,42.5,70,0.5,70.5);
00375 
00376   //  HFEtaFlagged    = fs->make<TProfile>( "HFEtaFlagged", "HFEtaFlagged",85,-42.5,42.5,0, 10000);
00377   HFEtaFlagged     = fs->make<TH1F>( "HFEtaFlagged", "HFEtaFlagged",85,-42.5,42.5);
00378   HFEtaFlaggedL    = fs->make<TH1F>( "HFEtaFlaggedL", "HFEtaFlaggedL",85,-42.5,42.5);
00379   HFEtaFlaggedLN    = fs->make<TH1F>( "HFEtaFlaggedLN", "HFEtaFlaggedLN",85,-42.5,42.5);
00380   HFEtaFlaggedS    = fs->make<TH1F>( "HFEtaFlaggedS", "HFEtaFlaggedS",85,-42.5,42.5);
00381   HFEtaFlaggedSN    = fs->make<TH1F>( "HFEtaFlaggedSN", "HFEtaFlaggedSN",85,-42.5,42.5);
00382 
00383   HFEtaNFlagged    = fs->make<TProfile>( "HFEtaNFlagged", "HFEtaNFlagged",85,-42.5,42.5,0, 10000);
00384 
00385   HOoccOOT    = fs->make<TH2F>( "HOoccOOT", "HOoccOOT",85,-42.5,42.5,70,0.5,70.5);
00386   HBoccOOT    = fs->make<TH2F>( "HBoccOOT", "HBoccOOT",85,-42.5,42.5,70,0.5,70.5);
00387   HEoccOOT    = fs->make<TH2F>( "HEoccOOT", "HEoccOOT",85,-42.5,42.5,70,0.5,70.5);
00388   HFoccOOT    = fs->make<TH2F>( "HFoccOOT", "HFoccOOT",85,-42.5,42.5,70,0.5,70.5);
00389 
00390   HFEnePMT0     = fs->make<TH1F>( "HFEnePMT0",  "HFEnePMT0", 210, -10, 200 );
00391   HFEnePMT1     = fs->make<TH1F>( "HFEnePMT1",  "HFEnePMT1", 210, -10, 200 );
00392   HFEnePMT2     = fs->make<TH1F>( "HFEnePMT2",  "HFEnePMT2", 210, -10, 200 );
00393   HFTimePMT0    = fs->make<TH1F>( "HFTimePMT0", "HFTimePMT0", 200, -100, 100 );
00394   HFTimePMT1    = fs->make<TH1F>( "HFTimePMT1", "HFTimePMT1", 200, -100, 100 );
00395   HFTimePMT2    = fs->make<TH1F>( "HFTimePMT2", "HFTimePMT2", 200, -100, 100 );
00396 
00397   HFEne     = fs->make<TH1F>( "HFEne",  "HFEne", 210, -10, 200 );
00398   HFEneFlagged     = fs->make<TH1F>( "HFEneFlagged",  "HFEneFlagged", 210, -10, 200 );
00399   HFEneFlagged2     = fs->make<TH1F>( "HFEneFlagged2",  "HFEneFlagged2", 210, -10, 200 );
00400   HFEneTh   = fs->make<TH1F>( "HFEneTh",  "HFEneTh", 210, -10, 200 );
00401   HFEneTh1  = fs->make<TH1F>( "HFEneTh1",  "HFEneTh1", 210, -10, 200 );
00402   HFEneP    = fs->make<TH1F>( "HFEneP",  "HFEneP", 200, -5, 10 );
00403   HFEneM    = fs->make<TH1F>( "HFEneM",  "HFEneM", 200, -5, 10 );
00404   HFTime    = fs->make<TH1F>( "HFTime", "HFTime", 200, -100, 100 );
00405   PMTHits   = fs->make<TH1F>( "PMTHits", "PMTHits", 10, 0, 10 );
00406   HFTimeFlagged  = fs->make<TH1F>( "HFTimeFlagged", "HFTimeFlagged", 200, -100, 100 );
00407 
00408   HFTimeFlagged2  = fs->make<TH1F>( "HFTimeFlagged2", "HFTimeFlagged2", 200, -100, 100 );
00409   HFTimeThFlagged2  = fs->make<TH1F>( "HFTimeThFlagged2", "HFTimeThFlagged2", 200, -100, 100 );
00410   HFTimeTh1Flagged2  = fs->make<TH1F>( "HFTimeTh1Flagged2", "HFTimeTh1Flagged2", 200, -100, 100 );
00411   HFTimeTh2Flagged2  = fs->make<TH1F>( "HFTimeTh2Flagged2", "HFTimeTh2Flagged2", 200, -100, 100 );
00412   HFTimeTh3Flagged2  = fs->make<TH1F>( "HFTimeTh3Flagged2", "HFTimeTh3Flagged2", 200, -100, 100 );
00413 
00414   HFTimeFlagged3  = fs->make<TH1F>( "HFTimeFlagged3", "HFTimeFlagged3", 200, -100, 100 );
00415   HFTimeThFlagged3  = fs->make<TH1F>( "HFTimeThFlagged3", "HFTimeThFlagged3", 200, -100, 100 );
00416   HFTimeTh1Flagged3  = fs->make<TH1F>( "HFTimeTh1Flagged3", "HFTimeTh1Flagged3", 200, -100, 100 );
00417   HFTimeTh2Flagged3  = fs->make<TH1F>( "HFTimeTh2Flagged3", "HFTimeTh2Flagged3", 200, -100, 100 );
00418   HFTimeTh3Flagged3  = fs->make<TH1F>( "HFTimeTh3Flagged3", "HFTimeTh3Flagged3", 200, -100, 100 );
00419 
00420   HFTimeThFlagged  = fs->make<TH1F>( "HFTimeThFlagged", "HFTimeThFlagged", 200, -100, 100 );
00421   HFTimeTh2Flagged  = fs->make<TH1F>( "HFTimeTh2Flagged", "HFTimeTh2Flagged", 200, -100, 100 );
00422   HFTimeTh3Flagged  = fs->make<TH1F>( "HFTimeTh3Flagged", "HFTimeTh3Flagged", 200, -100, 100 );
00423 
00424   HFTimeThFlaggedR  = fs->make<TH1F>( "HFTimeThFlaggedR", "HFTimeThFlaggedR", 200, -100, 100 );
00425   HFTimeThFlaggedR1  = fs->make<TH1F>( "HFTimeThFlaggedR1", "HFTimeThFlaggedR1", 200, -100, 100 );
00426   HFTimeThFlaggedR2  = fs->make<TH1F>( "HFTimeThFlaggedR2", "HFTimeThFlaggedR2", 200, -100, 100 );
00427   HFTimeThFlaggedR3  = fs->make<TH1F>( "HFTimeThFlaggedR3", "HFTimeThFlaggedR3", 200, -100, 100 );
00428   HFTimeThFlaggedR4  = fs->make<TH1F>( "HFTimeThFlaggedR4", "HFTimeThFlaggedR4", 200, -100, 100 );
00429   HFTimeThFlaggedRM  = fs->make<TH1F>( "HFTimeThFlaggedRM", "HFTimeThFlaggedRM", 200, -100, 100 );
00430   TrkMultFlagged0  = fs->make<TH1F>( "TrkMultFlagged0", "TrkMultFlagged0", 100, 0, 100 );
00431   TrkMultFlagged1  = fs->make<TH1F>( "TrkMultFlagged1", "TrkMultFlagged1", 100, 0, 100 );
00432   TrkMultFlagged2  = fs->make<TH1F>( "TrkMultFlagged2", "TrkMultFlagged2", 100, 0, 100 );
00433   TrkMultFlagged3  = fs->make<TH1F>( "TrkMultFlagged3", "TrkMultFlagged3", 100, 0, 100 );
00434   TrkMultFlagged4  = fs->make<TH1F>( "TrkMultFlagged4", "TrkMultFlagged4", 100, 0, 100 );
00435   TrkMultFlaggedM  = fs->make<TH1F>( "TrkMultFlaggedM", "TrkMultFlaggedM", 100, 0, 100 );
00436   HFTimeTh  = fs->make<TH1F>( "HFTimeTh", "HFTimeTh", 200, -100, 100 );
00437   HFTimeTh1  = fs->make<TH1F>( "HFTimeTh1", "HFTimeTh1", 200, -100, 100 );
00438   HFTimeTh2  = fs->make<TH1F>( "HFTimeTh2", "HFTimeTh2", 200, -100, 100 );
00439   HFTimeTh3  = fs->make<TH1F>( "HFTimeTh3", "HFTimeTh3", 200, -100, 100 );
00440   HFTimeThR  = fs->make<TH1F>( "HFTimeThR", "HFTimeThR", 200, -100, 100 );
00441   HFTimeTh1R  = fs->make<TH1F>( "HFTimeTh1R", "HFTimeTh1R", 200, -100, 100 );
00442   HFTimeTh2R  = fs->make<TH1F>( "HFTimeTh2R", "HFTimeTh2R", 200, -100, 100 );
00443   HFTimeTh3R  = fs->make<TH1F>( "HFTimeTh3R", "HFTimeTh3R", 200, -100, 100 );
00444   HFTimeP   = fs->make<TH1F>( "HFTimeP", "HFTimeP", 100, -100, 50 );
00445   HFTimeM   = fs->make<TH1F>( "HFTimeM", "HFTimeM", 100, -100, 50 );
00446   HFTimePMa = fs->make<TH1F>( "HFTimePMa", "HFTimePMa", 100, -100, 100 );
00447   HFTimePM  = fs->make<TH1F>( "HFTimePM", "HFTimePM", 100, -100, 100 );
00448 
00449   // Histos for separating HF long/short fibers:
00450   HFLEneAll  = fs->make<TH1F>( "HFLEneAll",  "HFLEneAll", 210, -10, 200 );
00451   HFLEneAllF = fs->make<TH1F>( "HFLEneAllF",  "HFLEneAllF", 210, -10, 200 );
00452   HFSEneAll  = fs->make<TH1F>( "HFSEneAll",  "HFSEneAll", 210, -10, 200 );
00453   HFSEneAllF = fs->make<TH1F>( "HFSEneAllF",  "HFSEneAllF", 210, -10, 200 );
00454   HFLEne     = fs->make<TH1F>( "HFLEne",  "HFLEne", 200, -5, 10 );
00455   HFLTime    = fs->make<TH1F>( "HFLTime", "HFLTime", 200, -100, 100 );
00456   HFSEne     = fs->make<TH1F>( "HFSEne",  "HFSEne", 200, -5, 10 );
00457   HFSTime    = fs->make<TH1F>( "HFSTime", "HFSTime", 200, -100, 100 );
00458   HFLSRatio  = fs->make<TH1F>( "HFLSRatio",  "HFLSRatio", 220, -1.1, 1.1 );
00459 
00460   HFOERatio  = fs->make<TH1F>( "HFOERatio",  "HFOERatio", 2200, -1.1, 1.1 );
00461 
00462   HFLvsS     = fs->make<TH2F>( "HFLvsS", "HFLvsS",220,-20,200,220,-20,200);
00463   HFLEneNoS  = fs->make<TH1F>( "HFLEneNoS",  "HFLEneNoS", 205, -5, 200 );
00464   HFSEneNoL  = fs->make<TH1F>( "HFSEneNoL",  "HFSEneNoL", 205, -5, 200 );
00465   HFLEneNoSFlagged  = fs->make<TH1F>( "HFLEneNoSFlagged",  "HFLEneNoSFlagged", 205, -5, 200 );
00466   HFSEneNoLFlagged  = fs->make<TH1F>( "HFSEneNoLFlagged",  "HFSEneNoLFlagged", 205, -5, 200 );
00467   HFLEneNoSFlaggedN  = fs->make<TH1F>( "HFLEneNoSFlaggedN",  "HFLEneNoSFlaggedN", 205, -5, 200 );
00468   HFSEneNoLFlaggedN  = fs->make<TH1F>( "HFSEneNoLFlaggedN",  "HFSEneNoLFlaggedN", 205, -5, 200 );
00469 
00470 
00471   EBEne     = fs->make<TH1F>( "EBEne",  "EBEne", 200, -5, 10 );
00472   EBEneTh   = fs->make<TH1F>( "EBEneTh",  "EBEneTh", 200, -5, 10 );
00473   EBEneX    = fs->make<TH1F>( "EBEneX",  "EBEneX", 200, -5, 10 );
00474   EBEneY    = fs->make<TH1F>( "EBEneY",  "EBEneY", 200, -5, 10 );
00475   EBTime    = fs->make<TH1F>( "EBTime", "EBTime", 200, -100, 100 );
00476   EBTimeTh  = fs->make<TH1F>( "EBTimeTh", "EBTimeTh", 200, -100, 100 );
00477   EBTimeX   = fs->make<TH1F>( "EBTimeX", "EBTimeX", 200, -100, 100 );
00478   EBTimeY   = fs->make<TH1F>( "EBTimeY", "EBTimeY", 200, -100, 100 );
00479   EEEne     = fs->make<TH1F>( "EEEne",  "EEEne", 200, -5, 10 );
00480   EEEneTh   = fs->make<TH1F>( "EEEneTh",  "EEEneTh", 200, -5, 10 );
00481   EEEneX    = fs->make<TH1F>( "EEEneX",  "EEEneX", 200, -5, 10 );
00482   EEEneY    = fs->make<TH1F>( "EEEneY",  "EEEneY", 200, -5, 10 );
00483   EEnegEne  = fs->make<TH1F>( "EEnegEne",  "EEnegEne", 200, -5, 10 );
00484   EEposEne  = fs->make<TH1F>( "EEposEne",  "EEposEne", 200, -5, 10 );
00485   EETime    = fs->make<TH1F>( "EETime", "EETime", 200, -100, 100 );
00486   EETimeTh  = fs->make<TH1F>( "EETimeTh", "EETimeTh", 200, -100, 100 );
00487   EETimeX   = fs->make<TH1F>( "EETimeX", "EETimeX", 200, -100, 100 );
00488   EETimeY   = fs->make<TH1F>( "EETimeY", "EETimeY", 200, -100, 100 );
00489   EEnegTime = fs->make<TH1F>( "EEnegTime", "EEnegTime", 200, -100, 100 );
00490   EEposTime = fs->make<TH1F>( "EEposTime", "EEposTime", 200, -100, 100 );
00491 
00492   h_nTowersCal = fs->make<TH1F>( "nTowersCal",  "N Towers in Jet", 100, 0, 50 );
00493   h_EMFracCal  = fs->make<TH1F>( "EMFracCal",  "EM Fraction in Jet", 100, -1.1, 1.1 );
00494   h_ptCal      = fs->make<TH1F>( "ptCal",  "p_{T} of CalJet", 100, 0, 50 );
00495   h_etaCal     = fs->make<TH1F>( "etaCal", "#eta of  CalJet", 100, -4, 4 );
00496   h_phiCal     = fs->make<TH1F>( "phiCal", "#phi of  CalJet", 50, -M_PI, M_PI );
00497 
00498   h_nGenJets  =  fs->make<TH1F>( "nGenJets",  "Number of GenJets", 20, 0, 20 );
00499 
00500   h_ptGen     =  fs->make<TH1F>( "ptGen",  "p_{T} of GenJet", 100, 0, 50 );
00501   h_etaGen    =  fs->make<TH1F>( "etaGen", "#eta of GenJet", 100, -4, 4 );
00502   h_phiGen    =  fs->make<TH1F>( "phiGen", "#phi of GenJet", 50, -M_PI, M_PI );
00503 
00504   h_ptGenL    =  fs->make<TH1F>( "ptGenL",  "p_{T} of GenJetL", 100, 0, 50 );
00505   h_etaGenL   =  fs->make<TH1F>( "etaGenL", "#eta of GenJetL", 100, -4, 4 );
00506   h_phiGenL   =  fs->make<TH1F>( "phiGenL", "#phi of GenJetL", 50, -M_PI, M_PI );
00507 
00508   h_jetEt      = fs->make<TH1F>( "jetEt", "Total Jet Et", 100, 0, 3000 );
00509 
00510   h_jet1Pt       = fs->make<TH1F>( "jet1Pt", "Jet1 Pt", 100, 0, 1000 );
00511   h_jet2Pt       = fs->make<TH1F>( "jet2Pt", "Jet2 Pt", 100, 0, 1000 );
00512   h_jet1Eta       = fs->make<TH1F>( "jet1Eta", "Jet1 Eta", 50, -5, 5 );
00513   h_jet2Eta       = fs->make<TH1F>( "jet2Eta", "Jet2 Eta", 50, -5, 5 );
00514   h_jet1PtHLT    = fs->make<TH1F>( "jet1PtHLT", "Jet1 Pt HLT", 100, 0, 1000 );
00515 
00516   h_TotalUnclusteredEt = fs->make<TH1F>( "TotalUnclusteredEt", "Total Unclustered Et", 100, 0, 500 );
00517   h_UnclusteredEt      = fs->make<TH1F>( "UnclusteredEt", "Unclustered Et", 100, 0, 50 );
00518   h_UnclusteredEts     = fs->make<TH1F>( "UnclusteredEts", "Unclustered Et", 100, 0, 2 );
00519 
00520   h_ClusteredE         = fs->make<TH1F>( "ClusteredE", "Clustered E", 200, 0, 20 );
00521   h_TotalClusteredE    = fs->make<TH1F>( "TotalClusteredE", "Total Clustered E", 200, 0, 100 );
00522   h_UnclusteredE       = fs->make<TH1F>( "UnclusteredE", "Unclustered E", 200, 0, 20 );
00523   h_TotalUnclusteredE  = fs->make<TH1F>( "TotalUnclusteredE", "Total Unclustered E", 200, 0, 100 );
00524 
00525   jetHOEne              = fs->make<TH1F>("jetHOEne" ,"HO Energy in Jet",100, 0,100);
00526   jetEMFraction         = fs->make<TH1F>( "jetEMFraction", "Jet EM Fraction", 100, -1.1, 1.1 );
00527   NTowers              = fs->make<TH1F>( "NTowers", "Number of Towers", 100, 0, 100 );
00528 
00529 
00530   h_EmEnergy   = fs->make<TH2F>( "EmEnergy",  "Em Energy",  90, -45, 45, 73, 0, 73 );
00531   h_HadEnergy  = fs->make<TH2F>( "HadEnergy", "Had Energy", 90, -45, 45, 73, 0, 73 );
00532 
00533   st_Pt            = fs->make<TH1F>( "st_Pt", "Pt", 200, 0, 200 );
00534   st_Constituents  = fs->make<TH1F>( "st_Constituents", "Constituents", 200, 0, 200 );
00535   st_Energy        = fs->make<TH1F>( "st_Energy", "Tower Energy", 200, 0, 200 );
00536   st_EmEnergy      = fs->make<TH1F>( "st_EmEnergy", "Tower EmEnergy", 200, 0, 200 );
00537   st_HadEnergy     = fs->make<TH1F>( "st_HadEnergy", "Tower HadEnergy", 200, 0, 200 );
00538   st_OuterEnergy   = fs->make<TH1F>( "st_OuterEnergy", "Tower OuterEnergy", 200, 0, 200 );
00539   st_Eta           = fs->make<TH1F>( "st_Eta", "Eta", 100, -4, 4 );
00540   st_Phi           = fs->make<TH1F>( "st_Phi", "Phi", 50, -M_PI, M_PI );
00541   st_iEta          = fs->make<TH1F>( "st_iEta", "iEta", 60, -30, 30 );
00542   st_iPhi          = fs->make<TH1F>( "st_iPhi", "iPhi", 80, 0, 80 );
00543   st_Frac          = fs->make<TH1F>( "st_Frac", "Frac", 100, 0, 1 );
00544 
00545 
00546   EBvHB           = fs->make<TH2F>( "EBvHB", "EB vs HB",1000,0,4500000.,1000,0,1000000.);
00547   EEvHE           = fs->make<TH2F>( "EEvHE", "EE vs HE",1000,0,4500000.,1000,0,200000.);
00548 
00549   ECALvHCAL       = fs->make<TH2F>( "ECALvHCAL", "ECAL vs HCAL",100,0,20000000.,100,-500000,500000.);
00550   ECALvHCALEta1   = fs->make<TH2F>( "ECALvHCALEta1", "ECAL vs HCALEta1",100,0,20000000.,100,-500000,500000.);
00551   ECALvHCALEta2   = fs->make<TH2F>( "ECALvHCALEta2", "ECAL vs HCALEta2",100,0,20000000.,100,-500000,500000.);
00552   ECALvHCALEta3   = fs->make<TH2F>( "ECALvHCALEta3", "ECAL vs HCALEta3",100,0,20000000.,100,-500000,500000.);
00553 
00554   EMF_Eta   = fs->make<TProfile>("EMF_Eta","EMF Eta", 100, -50, 50, 0, 10);
00555   EMF_Phi   = fs->make<TProfile>("EMF_Phi","EMF Phi", 100, 0, 100, 0, 10);
00556   EMF_EtaX  = fs->make<TProfile>("EMF_EtaX","EMF EtaX", 100, -50, 50, 0, 10);
00557   EMF_PhiX  = fs->make<TProfile>("EMF_PhiX","EMF PhiX", 100, 0, 100, 0, 10);
00558 
00559   HFTimeVsiEtaP  = fs->make<TProfile>("HFTimeVsiEtaP","HFTimeVsiEtaP", 13, 28.5, 41.5, -100, 100);
00560   HFTimeVsiEtaM  = fs->make<TProfile>("HFTimeVsiEtaM","HFTimeVsiEtaM", 13, -41.5, -28.5, -100, 100);
00561 
00562   HFTimeVsiEtaP5  = fs->make<TProfile>("HFTimeVsiEtaP5","HFTimeVsiEtaP5", 13, 28.5, 41.5, -100, 100);
00563   HFTimeVsiEtaM5  = fs->make<TProfile>("HFTimeVsiEtaM5","HFTimeVsiEtaM5", 13, -41.5, -28.5, -100, 100);
00564 
00565   HFTimeVsiEtaP20  = fs->make<TProfile>("HFTimeVsiEtaP20","HFTimeVsiEtaP20", 13, 28.5, 41.5, -100, 100);
00566   HFTimeVsiEtaM20  = fs->make<TProfile>("HFTimeVsiEtaM20","HFTimeVsiEtaM20", 13, -41.5, -28.5, -100, 100);
00567 
00568   NPass          = fs->make<TH1F>( "NPass", "NPass", 3, -1, 1 );
00569   NTotal         = fs->make<TH1F>( "NTotal", "NTotal", 3, -1, 1 );
00570   NTime         = fs->make<TH1F>( "NTime", "NTime", 10, 0, 10 );
00571 
00572 
00573   HFRecHitEne   = fs->make<TH1F>( "HFRecHitEne", "HFRecHitEne", 300, 0, 3000 );
00574   HFRecHitEneClean   = fs->make<TH1F>( "HFRecHitEneClean", "HFRecHitEneClean", 300, 0, 3000 );
00575   HFRecHitTime  = fs->make<TH1F>( "HFRecHitTime", "HFRecHitTime", 120, -60, 60 );
00576 
00577 
00578   HFLongShortPhi  = fs->make<TH1F>( "HFLongShortPhi", "HFLongShortPhi", 73, 0, 73 );
00579   HFLongShortEta  = fs->make<TH1F>( "HFLongShortEta", "HFLongShortEta", 90, -45, 45 );
00580   HFLongShortEne  = fs->make<TH1F>( "HFLongShortEne", "HFLongShortEne", 300, 0, 3000 );
00581   HFLongShortTime = fs->make<TH1F>( "HFLongShortTime", "HFLongShortTime", 120, -60, 60 );
00582   HFLongShortNHits = fs->make<TH1F>( "HFLongShortNHits", "HFLongShortNHits", 30, 0, 30 );
00583 
00584   HFDigiTimePhi  = fs->make<TH1F>( "HFDigiTimePhi", "HFDigiTimePhi", 73, 0, 73 );
00585   HFDigiTimeEta  = fs->make<TH1F>( "HFDigiTimeEta", "HFDigiTimeEta", 90, -45, 45 );
00586   HFDigiTimeEne  = fs->make<TH1F>( "HFDigiTimeEne", "HFDigiTimeEne", 300, 0, 3000 );
00587   HFDigiTimeTime = fs->make<TH1F>( "HFDigiTimeTime", "HFDigiTimeTime", 120, -60, 60 );
00588   HFDigiTimeNHits = fs->make<TH1F>( "HFDigiTimeNHits", "HFDigiTimeNHits", 30, 0, 30 );
00589 
00590 
00591   totBNC = 0;
00592   for (int i=0; i<4000; i++)  nBNC[i] = 0;
00593 
00594 }
00595 
00596 // ************************
00597 // ************************
00598 void myJetAna::analyze( const edm::Event& evt, const edm::EventSetup& es ) {
00599  
00600   using namespace edm;
00601 
00602   bool Pass, Pass_HFTime, Pass_DiJet, Pass_BunchCrossing, Pass_Vertex;
00603 
00604   int EtaOk10, EtaOk13, EtaOk40;
00605 
00606   double LeadMass;
00607 
00608   double HFRecHit[100][100][2];
00609   int HFRecHitFlag[100][100][2];
00610 
00611   double towerEtCut, towerECut, towerE;
00612 
00613   towerEtCut = 1.0;
00614   towerECut  = 1.0;
00615 
00616   unsigned int StableRun = 123732;
00617 
00618   double HBHEThreshold = 2.0;
00619   double HFThreshold   = 2.0;
00620   double HOThreshold   = 2.0;
00621   double EBEEThreshold = 2.0;
00622 
00623   double HBHEThreshold1 = 4.0;
00624   double HFThreshold1   = 4.0;
00625   double HOThreshold1   = 4.0;
00626   //double EBEEThreshold1 = 4.0;
00627 
00628   double HBHEThreshold2 = 10.0;
00629   double HFThreshold2   = 10.0;
00630   //double HOThreshold2   = 10.0;
00631   //double EBEEThreshold2 = 10.0;
00632 
00633   double HBHEThreshold3 = 40.0;
00634   double HFThreshold3   = 40.0;
00635   //double HOThreshold3   = 40.0;
00636   //double EBEEThreshold3 = 40.0;
00637 
00638   float minJetPt = 20.;
00639   float minJetPt10 = 10.;
00640   int jetInd, allJetInd;
00641   LeadMass = -1;
00642 
00643   //  Handle<DcsStatusCollection> dcsStatus;
00644   //  evt.getByLabel("scalersRawToDigi", dcsStatus);
00645   //  std::cout << dcsStatus << std::endl;
00646   //  if (dcsStatus.isValid()) {
00647   //  }
00648 
00649   //  DcsStatus dcsStatus;
00650   //  Handle<DcsStatus> dcsStatus;
00651   //  evt.getByLabel("dcsStatus", dcsStatus);
00652 
00653 
00654   math::XYZTLorentzVector p4tmp[2], p4cortmp[2];
00655 
00656   // --------------------------------------------------------------
00657   // --------------------------------------------------------------
00658 
00659   /***
00660   std::cout << ">>>> ANA: Run = "    << evt.id().run() 
00661             << " Event = " << evt.id().event()
00662             << " Bunch Crossing = " << evt.bunchCrossing() 
00663             << " Orbit Number = "   << evt.orbitNumber()
00664             << " Luminosity Block = "  << evt.luminosityBlock()
00665             <<  std::endl;
00666   ***/
00667 
00668   // *********************
00669   // *** Filter Event
00670   // *********************
00671   Pass = false;
00672 
00673   /***
00674   if (evt.bunchCrossing()== 100) {
00675     Pass = true;
00676   } else {
00677     Pass = false;
00678   }
00679   ***/
00680 
00681   // ***********************
00682   // ***  Pass Trigger
00683   // ***********************
00684 
00685 
00686   // **** Get the TriggerResults container
00687   Handle<TriggerResults> triggerResults;
00688   evt.getByLabel(theTriggerResultsLabel, triggerResults);
00689   //  evt.getByLabel("TriggerResults::HLT", triggerResults);
00690 
00691   if (triggerResults.isValid()) {
00692     if (DEBUG) std::cout << "trigger valid " << std::endl;
00693     //    edm::TriggerNames triggerNames;    // TriggerNames class
00694     //    triggerNames.init(*triggerResults);
00695     unsigned int n = triggerResults->size();
00696     for (unsigned int i=0; i!=n; i++) {
00697 
00698       /***
00699       std::cout << ">>> Trigger Name (" <<  i << ") = " << triggerNames.triggerName(i)
00700                 << " Accept = " << triggerResults->accept(i)
00701                 << std::endl;
00702       ***/
00703       /****
00704       if (triggerResults->accept(i) == 1) {
00705         std::cout << "+++ Trigger Name (" <<  i << ") = " << triggerNames.triggerName(i)
00706                   << " Accept = " << triggerResults->accept(i)
00707                   << std::endl;
00708       }
00709       ****/
00710 
00711       //      if (DEBUG) std::cout <<  triggerNames.triggerName(i) << std::endl;
00712 
00713       //      if ( (triggerNames.triggerName(i) == "HLT_ZeroBias")  || 
00714       //           (triggerNames.triggerName(i) == "HLT_MinBias")   || 
00715       //           (triggerNames.triggerName(i) == "HLT_MinBiasHcal") )  {
00716 
00717     }
00718       
00719   } else {
00720 
00721     edm::Handle<TriggerResults> *tr = new edm::Handle<TriggerResults>;
00722     triggerResults = (*tr);
00723 
00724     //     std::cout << "triggerResults is not valid" << std::endl;
00725     //     std::cout << triggerResults << std::endl;
00726     //     std::cout << triggerResults.isValid() << std::endl;
00727     
00728     if (DEBUG) std::cout << "trigger not valid " << std::endl;
00729     edm::LogInfo("myJetAna") << "TriggerResults::HLT not found, "
00730       "automatically select events";
00731 
00732     //return;
00733   }
00734 
00735 
00736   
00737   /***
00738   Handle<L1GlobalTriggerReadoutRecord> gtRecord;
00739   evt.getByLabel("gtDigis",gtRecord);
00740   const TechnicalTriggerWord tWord = gtRecord->technicalTriggerWord();
00741 
00742   ***/
00743 
00744 
00745   // *************************
00746   // ***  Pass Bunch Crossing
00747   // *************************
00748 
00749   // *** Check Luminosity Section
00750   if (evt.id().run() == 122294)
00751     if ( (evt.luminosityBlock() >= 37) && (evt.luminosityBlock() <= 43) ) 
00752       Pass = true;
00753   if (evt.id().run() == 122314)
00754     if ( (evt.luminosityBlock() >= 24) && (evt.luminosityBlock() <= 37) ) 
00755       Pass = true;
00756   if (evt.id().run() == 123575)
00757       Pass = true;
00758   if (evt.id().run() == 123596)
00759       Pass = true;
00760 
00761   // ***********
00762   if (evt.id().run() == 124009) {
00763     if ( (evt.bunchCrossing() == 51) ||
00764          (evt.bunchCrossing() == 151) ||
00765          (evt.bunchCrossing() == 2824) ) {
00766       Pass = true;
00767     }
00768   }
00769 
00770   if (evt.id().run() == 124020) {
00771     if ( (evt.bunchCrossing() == 51) ||
00772          (evt.bunchCrossing() == 151) ||
00773          (evt.bunchCrossing() == 2824) ) {
00774       Pass = true;
00775     }
00776   }
00777 
00778   if (evt.id().run() == 124024) {
00779     if ( (evt.bunchCrossing() == 51) ||
00780          (evt.bunchCrossing() == 151) ||
00781          (evt.bunchCrossing() == 2824) ) {
00782       Pass = true;
00783     }
00784   }
00785 
00786   if ( (evt.bunchCrossing() == 51)   ||
00787        (evt.bunchCrossing() == 151)  ||
00788        (evt.bunchCrossing() == 2596) || 
00789        (evt.bunchCrossing() == 2724) || 
00790        (evt.bunchCrossing() == 2824) ||
00791        (evt.bunchCrossing() == 3487) ) {
00792     Pass_BunchCrossing = true;
00793   } else {
00794     Pass_BunchCrossing = false;
00795   }
00796   
00797 
00798   // ***********************
00799   // ***  Pass HF Timing 
00800   // ***********************
00801 
00802   double HFM_ETime, HFP_ETime;
00803   double HFM_E, HFP_E;
00804   double HF_PMM;
00805 
00806   HFM_ETime = 0.;
00807   HFM_E = 0.;
00808   HFP_ETime = 0.;
00809   HFP_E = 0.;
00810 
00811   for (int i=0; i<100; i++) {
00812     for (int j=0; j<100; j++) {
00813       HFRecHit[i][j][0] = -10.;
00814       HFRecHit[i][j][1] = -10.;
00815       HFRecHitFlag[i][j][0]  = 0;
00816       HFRecHitFlag[i][j][1]  = 0;
00817     }
00818   }
00819 
00820 
00821   int nTime = 0;
00822   int NHFLongShortHits;
00823   int NHFDigiTimeHits;
00824   NHFLongShortHits = 0;
00825   NHFDigiTimeHits = 0;
00826 
00827   //  edm::Handle<reco::VertexCollection> vertexCollection;
00828 
00829   try {
00830     std::vector<edm::Handle<HFRecHitCollection> > colls;
00831     evt.getManyByType(colls);
00832 
00833     std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
00834     for (i=colls.begin(); i!=colls.end(); i++) {
00835       
00836       for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
00837         if (j->id().subdet() == HcalForward) {
00838 
00839           HFRecHitEne->Fill(j->energy());
00840           if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 0) && 
00841                (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 0) ) {
00842             HFRecHitEneClean->Fill(j->energy());
00843           }
00844  
00845           HFRecHitTime->Fill(j->time());
00846 
00847           int myFlag;
00848           myFlag= j->flagField(HcalCaloFlagLabels::HFLongShort);
00849           if (myFlag==1) {
00850             NHFLongShortHits++;
00851             HFLongShortPhi->Fill(j->id().iphi());
00852             HFLongShortEta->Fill(j->id().ieta());
00853             HFLongShortEne->Fill(j->energy());
00854             HFLongShortTime->Fill(j->time());
00855           }
00856 
00857           myFlag= j->flagField(HcalCaloFlagLabels::HFDigiTime);
00858           if (myFlag==1) {
00859             NHFDigiTimeHits++;
00860             HFDigiTimePhi->Fill(j->id().iphi());
00861             HFDigiTimeEta->Fill(j->id().ieta());
00862             HFDigiTimeEne->Fill(j->energy());
00863             HFDigiTimeTime->Fill(j->time());
00864           }
00865 
00866           
00867           float en = j->energy();
00868           float time = j->time();
00869           if ((en > 20.) && (time>20.)) {
00870             HFoccTime->Fill(j->id().ieta(),j->id().iphi());
00871             nTime++;
00872           }
00873           HcalDetId id(j->detid().rawId());
00874           int ieta = id.ieta();
00875           int iphi = id.iphi();
00876           int depth = id.depth();
00877 
00878 
00879           // Long:  depth = 1
00880           // Short: depth = 2
00881           HFRecHit[ieta+41][iphi][depth-1] = en;
00882           HFRecHitFlag[ieta+41][iphi][depth-1] = j->flagField(0);
00883 
00884           /****
00885           std::cout << "RecHit Flag = " 
00886                     << j->flagField(0)a
00887                     << std::endl;
00888           ***/
00889 
00890           if (j->id().ieta()<0) {
00891             if (j->energy() > HFThreshold) {
00892               HFM_ETime += j->energy()*j->time(); 
00893               HFM_E     += j->energy();
00894             }
00895           } else {
00896             if (j->energy() > HFThreshold) {
00897               HFP_ETime += j->energy()*j->time(); 
00898               HFP_E     += j->energy();
00899             }
00900           }
00901 
00902         }
00903       }
00904       break;
00905     }
00906   } catch (...) {
00907     cout << "No HF RecHits." << endl;
00908   }
00909 
00910   cout << "N HF Hits" << NHFLongShortHits << " " << NHFDigiTimeHits << endl;
00911   HFLongShortNHits->Fill(NHFLongShortHits);
00912   HFDigiTimeNHits->Fill(NHFDigiTimeHits);
00913 
00914   NTime->Fill(nTime);
00915 
00916   double OER = 0, OddEne, EvenEne;
00917   int nOdd, nEven;
00918 
00919   for (int iphi=0; iphi<100; iphi++) {
00920     OddEne = EvenEne = 0.;
00921     nOdd  = 0;
00922     nEven = 0;
00923     for (int ieta=0; ieta<100; ieta++) {
00924       if (HFRecHit[ieta][iphi][0] > 1.0) {
00925         if (ieta%2 == 0) {
00926           EvenEne += HFRecHit[ieta][iphi][0]; 
00927           nEven++;
00928         } else {
00929           OddEne  += HFRecHit[ieta][iphi][0];
00930           nOdd++;
00931         }
00932       }
00933       if (HFRecHit[ieta][iphi][1] > 1.0) {
00934         if (ieta%2 == 0) {
00935           EvenEne += HFRecHit[ieta][iphi][1]; 
00936           nEven++;
00937         } else {
00938           OddEne  += HFRecHit[ieta][iphi][1]; 
00939           nOdd++;
00940         }
00941       }
00942     }
00943     if (((OddEne + EvenEne) > 10.) && (nOdd > 1) && (nEven > 1)) {
00944       OER = (OddEne - EvenEne) / (OddEne + EvenEne);
00945       HFOERatio->Fill(OER);
00946     }
00947   }
00948 
00949   if ((HFP_E > 0.) && (HFM_E > 0.)) {
00950     HF_PMM = (HFP_ETime / HFP_E) - (HFM_ETime / HFM_E);
00951     HFTimePMa->Fill(HF_PMM); 
00952   } else {
00953     HF_PMM = INVALID;
00954   }
00955 
00956   
00957   if (fabs(HF_PMM) < 10.)  {
00958     Pass_HFTime = true;
00959   } else {
00960     Pass_HFTime = false;
00961   }
00962 
00963 
00964   // **************************
00965   // ***  Pass DiJet Criteria
00966   // **************************
00967   double highestPt;
00968   double nextPt;
00969   // double dphi;
00970   int    nDiJet, nJet;
00971 
00972   nJet      = 0;
00973   nDiJet    = 0;
00974   highestPt = 0.0;
00975   nextPt    = 0.0;
00976 
00977   allJetInd = 0;
00978   Handle<CaloJetCollection> caloJets;
00979   evt.getByLabel( CaloJetAlgorithm, caloJets );
00980   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) {
00981 
00982     // TODO: verify first two jets are the leading jets
00983     if (nJet == 0) p4tmp[0] = cal->p4();
00984     if (nJet == 1) p4tmp[1] = cal->p4();
00985 
00986     if ( (cal->pt() > 3.) && 
00987          (fabs(cal->eta()) < 3.0) ) { 
00988       nDiJet++;
00989     }
00990     nJet++;
00991 
00992   }  
00993   
00994 
00995   if (nDiJet > 1) {
00996     //dphi = deltaPhi(p4tmp[0].phi(), p4tmp[1].phi());
00997     Pass_DiJet = true;
00998   } else {
00999     // dphi = INVALID;
01000     Pass_DiJet = false;
01001   }
01002       
01003 
01004   // **************************
01005   // ***  Pass Vertex
01006   // **************************
01007   double VTX;
01008   int nVTX;
01009 
01010   edm::Handle<reco::VertexCollection> vertexCollection;
01011   evt.getByLabel("offlinePrimaryVertices", vertexCollection);
01012   const reco::VertexCollection vC = *(vertexCollection.product());
01013 
01014   //  std::cout << "Reconstructed "<< vC.size() << " vertices" << std::endl ;
01015 
01016   nVTX = vC.size();
01017   for (reco::VertexCollection::const_iterator vertex=vC.begin(); vertex!=vC.end(); vertex++){
01018     VTX  = vertex->z();
01019   }
01020 
01021   if ( (fabs(VTX) < 20.) && (nVTX > 0) ){
01022     Pass_Vertex = true;
01023   } else {
01024     Pass_Vertex = false;
01025   }
01026 
01027   // ***********************
01028   // ***********************
01029 
01030 
01031   nBNC[evt.bunchCrossing()]++;
01032   totBNC++;
01033     
01034   //  Pass = true;
01035 
01036   // *** Check for tracks
01037   //  edm::Handle<reco::TrackCollection> trackCollection;
01038   //  evt.getByLabel("generalTracks", trackCollection);
01039   //  const reco::TrackCollection tC = *(trackCollection.product());
01040   //  if ((Pass) && (tC.size()>1)) {
01041   //  } else {
01042   //    Pass = false;
01043   //  } 
01044 
01045 
01046   // ********************************
01047   // *** Pixel Clusters
01048   // ********************************
01049   edm::Handle< edmNew::DetSetVector<SiPixelCluster> > hClusterColl;
01050   evt.getByLabel("siPixelClusters", hClusterColl);
01051   const edmNew::DetSetVector<SiPixelCluster> clustColl = *(hClusterColl.product());
01052 
01053   edm::Handle<reco::TrackCollection> trackCollection;
01054   evt.getByLabel("generalTracks", trackCollection);
01055   const reco::TrackCollection tC = *(trackCollection.product());
01056 
01057 
01058   // **************************
01059   // *** Event Passed Selection
01060   // **************************
01061 
01062 
01063   if (evt.id().run() == 1) {
01064     if ( (Pass_DiJet)         &&
01065          (Pass_Vertex) ) {
01066       Pass = true;
01067     } else {
01068       Pass = false;
01069     }
01070     Pass = true;
01071 
01072   } else {
01073     if ( (Pass_BunchCrossing) && 
01074          (Pass_HFTime)        &&
01075          (Pass_Vertex) ) {
01076       Pass = true;
01077     } else {
01078       Pass = false;
01079     }
01080   }
01081 
01082   /***
01083   std::cout << "+++ Result " 
01084             << " Event = " 
01085             << evt.id().run()
01086             << " LS = "
01087             << evt.luminosityBlock()
01088             << " dphi = "
01089             << dphi
01090             << " Pass = " 
01091             << Pass
01092             << std::endl;
01093   ***/
01094 
01095   NTotal->Fill(0);
01096   
01097   Pass = false;
01098   if ((tC.size() > 100) && (clustColl.size() > 1000)) Pass = true;
01099   Pass = true;
01100 
01101   /****
01102   if (Pass_HFTime) {
01103     Pass = true;
01104   } else {
01105     Pass = false;
01106   }
01107   ****/
01108 
01109   // **************************
01110   // *** Noise Summary Object
01111   // **************************
01112 
01113   edm::Handle<HcalNoiseSummary> summary_h;
01114   evt.getByType(summary_h);
01115   if(!summary_h.isValid()) {
01116     throw edm::Exception(edm::errors::ProductNotFound) << " could not find HcalNoiseSummary.\n";
01117     //    return true;
01118   }
01119 
01120   const HcalNoiseSummary summary = *summary_h;
01121 
01122   bool Pass_NoiseSummary;
01123   Pass_NoiseSummary = true;
01124   if(summary.minE2Over10TS()<0.7) {
01125     Pass_NoiseSummary = false;
01126   }
01127   if(summary.maxE2Over10TS()>0.96) {
01128     Pass_NoiseSummary = false;
01129   }
01130   if(summary.maxHPDHits()>=17) {
01131     Pass_NoiseSummary = false;
01132   }
01133   if(summary.maxRBXHits()>=999) {
01134     Pass_NoiseSummary = false;
01135   }
01136   if(summary.maxHPDNoOtherHits()>=10) {
01137     Pass_NoiseSummary = false;
01138   }
01139   if(summary.maxZeros()>=10) {
01140     Pass_NoiseSummary = false;
01141   }
01142   if(summary.min25GeVHitTime()<-9999.0) {
01143     Pass_NoiseSummary = false;
01144   }
01145   if(summary.max25GeVHitTime()>9999.0) {
01146     Pass_NoiseSummary = false;
01147   }
01148   if(summary.minRBXEMF()<0.01) {
01149   }
01150 
01151   if (Pass_NoiseSummary) {
01152     Pass = false;
01153   } else {
01154     Pass = true;
01155   }
01156 
01157 
01158   Pass = true;
01159   if (Pass) {
01160 
01161     NPass->Fill(0);
01162 
01163   // *********************
01164   // *** Classify Event
01165   // *********************
01166   int evtType = 0;
01167 
01168   Handle<CaloTowerCollection> caloTowers;
01169   evt.getByLabel( "towerMaker", caloTowers );
01170 
01171   for (int i=0;i<36;i++) {
01172     RBXColl[i].et        = 0;
01173     RBXColl[i].hadEnergy = 0;
01174     RBXColl[i].emEnergy  = 0;
01175     RBXColl[i].hcalTime  = 0;
01176     RBXColl[i].ecalTime  = 0;
01177     RBXColl[i].nTowers   = 0;
01178   }  
01179   for (int i=0;i<144;i++) {
01180     HPDColl[i].et        = 0;
01181     HPDColl[i].hadEnergy = 0;
01182     HPDColl[i].emEnergy  = 0;
01183     HPDColl[i].hcalTime  = 0;
01184     HPDColl[i].ecalTime  = 0;
01185     HPDColl[i].nTowers   = 0;
01186   }  
01187 
01188   double ETotal, emFrac;
01189   double HCALTotalCaloTowerE, ECALTotalCaloTowerE;
01190   double HCALTotalCaloTowerE_Eta1, ECALTotalCaloTowerE_Eta1;
01191   double HCALTotalCaloTowerE_Eta2, ECALTotalCaloTowerE_Eta2;
01192   double HCALTotalCaloTowerE_Eta3, ECALTotalCaloTowerE_Eta3;
01193 
01194   ETotal = 0.;
01195   emFrac = 0.;
01196     
01197   HCALTotalCaloTowerE = 0;
01198   ECALTotalCaloTowerE = 0;
01199   HCALTotalCaloTowerE_Eta1 = 0.;
01200   ECALTotalCaloTowerE_Eta1 = 0.;
01201   HCALTotalCaloTowerE_Eta2 = 0.; 
01202   ECALTotalCaloTowerE_Eta2 = 0.;
01203   HCALTotalCaloTowerE_Eta3 = 0.;
01204   ECALTotalCaloTowerE_Eta3 = 0.;
01205 
01206   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
01207        tower != caloTowers->end(); tower++) {
01208     ETotal += tower->hadEnergy();
01209     ETotal += tower->emEnergy();
01210   }
01211 
01212   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
01213        tower != caloTowers->end(); tower++) {
01214 
01215     // Raw tower energy without grouping or thresholds
01216     if (abs(tower->ieta()) < 100) EMF_Eta->Fill(tower->ieta(), emFrac);
01217 
01218     if (abs(tower->ieta()) < 15) {
01219       towerHadEnHB->Fill(tower->hadEnergy());
01220       towerEmEnHB->Fill(tower->emEnergy());
01221     }
01222     if ( (abs(tower->ieta()) > 17) && ((abs(tower->ieta()) < 30)) ){
01223       towerHadEnHE->Fill(tower->hadEnergy());
01224       towerEmEnHE->Fill(tower->emEnergy());
01225     }
01226     if (abs(tower->ieta()) > 29) {
01227       towerHadEnHF->Fill(tower->hadEnergy());
01228       towerEmEnHF->Fill(tower->emEnergy());
01229     }
01230 
01231     towerHadEn->Fill(tower->hadEnergy());
01232     towerEmEn->Fill(tower->emEnergy());
01233     towerOuterEn->Fill(tower->outerEnergy());
01234 
01235     //    towerHadEt->Fill(tower->hadEt());
01236     //    towerEmEt->Fill(tower->emEt());
01237     //    towerOuterEt->Fill(tower->outerEt());
01238 
01239     if ((tower->emEnergy()+tower->hadEnergy()) != 0) {
01240       emFrac = tower->emEnergy()/(tower->emEnergy()+tower->hadEnergy());
01241       towerEmFrac->Fill(emFrac);
01242     } else {
01243       emFrac = 0.;
01244     }
01245 
01246     /***
01247       std::cout << "ETotal = " << ETotal 
01248                 << " EMF = " << emFrac
01249                 << " EM = " << tower->emEnergy()
01250                 << " Tot = " << tower->emEnergy()+tower->hadEnergy()
01251                 << " ieta/iphi = " <<  tower->ieta() << " / "  << tower->iphi() 
01252                 << std::endl;
01253     ***/
01254 
01255     if (abs(tower->iphi()) < 100) EMF_Phi->Fill(tower->iphi(), emFrac);
01256     if (abs(tower->ieta()) < 100) EMF_Eta->Fill(tower->ieta(), emFrac);
01257     if ( (evt.id().run() == 120020) && (evt.id().event() == 453) ) {
01258       std::cout << "Bunch Crossing = " << evt.bunchCrossing() 
01259                 << " Orbit Number = "  << evt.orbitNumber()
01260                 <<  std::endl;
01261 
01262       if (abs(tower->iphi()) < 100) EMF_PhiX->Fill(tower->iphi(), emFrac);
01263       if (abs(tower->ieta()) < 100) EMF_EtaX->Fill(tower->ieta(), emFrac);
01264     }
01265     
01266     HCALTotalCaloTowerE += tower->hadEnergy();
01267     ECALTotalCaloTowerE += tower->emEnergy();
01268 
01269     towerE = tower->hadEnergy() + tower->emEnergy();
01270     if (tower->et() > towerEtCut) caloEtaEt->Fill(tower->eta());
01271     if (towerE      > towerECut)  caloEta->Fill(tower->eta());
01272     caloPhi->Fill(tower->phi());
01273 
01274     if (fabs(tower->eta()) < 1.3) {
01275       HCALTotalCaloTowerE_Eta1 += tower->hadEnergy();
01276       ECALTotalCaloTowerE_Eta1 += tower->emEnergy();
01277     }
01278     if ((fabs(tower->eta()) >= 1.3) && (fabs(tower->eta()) < 2.5)) {
01279       HCALTotalCaloTowerE_Eta2 += tower->hadEnergy();
01280       ECALTotalCaloTowerE_Eta2 += tower->emEnergy();
01281     }
01282     if (fabs(tower->eta()) > 2.5) {
01283       HCALTotalCaloTowerE_Eta3 += tower->hadEnergy();
01284       ECALTotalCaloTowerE_Eta3 += tower->emEnergy();
01285     }
01286 
01287     /***
01288     std::cout << "had = "  << tower->hadEnergy()
01289               << " em = "  << tower->emEnergy()
01290               << " fabs(eta) = " << fabs(tower->eta())
01291               << " ieta/iphi = " <<  tower->ieta() << " / "  << tower->iphi() 
01292               << std::endl;
01293     ***/
01294 
01295     if ((tower->hadEnergy() + tower->emEnergy()) > 2.0) {
01296 
01297       int iRBX = tower->iphi();
01298       iRBX = iRBX-2;
01299       if (iRBX == 0)  iRBX = 17;
01300       if (iRBX == -1) iRBX = 18;
01301       iRBX = (iRBX-1)/4;
01302 
01303       if (tower->ieta() < 0) iRBX += 18;
01304       if (iRBX < 36) {
01305         RBXColl[iRBX].et        += tower->et(); 
01306         RBXColl[iRBX].hadEnergy += tower->hadEnergy(); 
01307         RBXColl[iRBX].emEnergy  += tower->emEnergy(); 
01308         RBXColl[iRBX].hcalTime  += tower->hcalTime(); 
01309         RBXColl[iRBX].ecalTime  += tower->ecalTime(); 
01310         RBXColl[iRBX].nTowers++;
01311       }
01312       /***
01313       std::cout << "iRBX = " << iRBX << " "     
01314                 << "ieta/iphi = " <<  tower->ieta() << " / "  << tower->iphi() 
01315                 << " et = " << tower->et()
01316                 << std::endl;
01317       ***/
01318       int iHPD = tower->iphi();
01319       if (tower->ieta() < 0) iHPD = iHPD + 72;
01320       if (iHPD < 144) {
01321         HPDColl[iHPD].et        += tower->et(); 
01322         HPDColl[iHPD].hadEnergy += tower->hadEnergy(); 
01323         HPDColl[iHPD].emEnergy  += tower->emEnergy(); 
01324         HPDColl[iHPD].hcalTime  += tower->hcalTime(); 
01325         HPDColl[iHPD].ecalTime  += tower->ecalTime(); 
01326         HPDColl[iHPD].nTowers++;
01327       }
01328       /***
01329       std::cout << "iHPD = " << iHPD << " "     
01330                 << "ieta/iphi = " <<  tower->ieta() << " / "  << tower->iphi() 
01331                 << " et = " << tower->et()
01332                 << std::endl;
01333       ***/
01334 
01335     }
01336 
01337   }
01338 
01339   ECALvHCAL->Fill(HCALTotalCaloTowerE, ECALTotalCaloTowerE);
01340   ECALvHCALEta1->Fill(HCALTotalCaloTowerE_Eta1, ECALTotalCaloTowerE_Eta1);
01341   ECALvHCALEta2->Fill(HCALTotalCaloTowerE_Eta2, ECALTotalCaloTowerE_Eta2);
01342   ECALvHCALEta3->Fill(HCALTotalCaloTowerE_Eta3, ECALTotalCaloTowerE_Eta3);
01343 
01344   /***
01345   std::cout << " Total CaloTower Energy :  "
01346             << " ETotal= " << ETotal 
01347             << " HCAL= " << HCALTotalCaloTowerE 
01348             << " ECAL= " << ECALTotalCaloTowerE
01349             << std::endl;
01350   ***/
01351 
01352   /***
01353             << " HCAL Eta1 = "  << HCALTotalCaloTowerE_Eta1
01354             << " ECAL= " << ECALTotalCaloTowerE_Eta1
01355             << " HCAL Eta2 = " << HCALTotalCaloTowerE_Eta2
01356             << " ECAL= " << ECALTotalCaloTowerE_Eta2
01357             << " HCAL Eta3 = " << HCALTotalCaloTowerE_Eta3
01358             << " ECAL= " << ECALTotalCaloTowerE_Eta3
01359             << std::endl;
01360   ***/
01361 
01362 
01363   // Loop over the RBX Collection
01364   int nRBX = 0;
01365   int nTowers = 0;
01366   for (int i=0;i<36;i++) {
01367     RBX_et->Fill(RBXColl[i].et);
01368     RBX_hadEnergy->Fill(RBXColl[i].hadEnergy);
01369     RBX_hcalTime->Fill(RBXColl[i].hcalTime / RBXColl[i].nTowers);
01370     RBX_nTowers->Fill(RBXColl[i].nTowers);
01371     if (RBXColl[i].hadEnergy > 3.0) {
01372       nRBX++;
01373       nTowers = RBXColl[i].nTowers;
01374     }
01375   }
01376   RBX_N->Fill(nRBX);
01377   if ( (nRBX == 1) && (nTowers > 24) ) {
01378     evtType = 1;
01379   }
01380 
01381   // Loop over the HPD Collection
01382   int nHPD = 0;
01383   for (int i=0;i<144;i++) {
01384     HPD_et->Fill(HPDColl[i].et);
01385     HPD_hadEnergy->Fill(HPDColl[i].hadEnergy);
01386     HPD_hcalTime->Fill(HPDColl[i].hcalTime / HPDColl[i].nTowers);
01387     HPD_nTowers->Fill(HPDColl[i].nTowers);
01388     if (HPDColl[i].hadEnergy > 3.0) {
01389       nHPD++;
01390       nTowers = HPDColl[i].nTowers;
01391     }
01392   }
01393   HPD_N->Fill(nHPD);
01394   if ( (nHPD == 1) && (nTowers > 6) ) {
01395     evtType = 2;
01396     cout << " nHPD = "   << nHPD 
01397          << " Towers = " << nTowers
01398          << " Type = "   << evtType 
01399          << endl; 
01400   }
01401  
01402   // **************************************************************
01403   // ** Access Trigger Information
01404   // **************************************************************
01405 
01406   // **** Get the TriggerResults container
01407   Handle<TriggerResults> triggerResults;
01408   evt.getByLabel(theTriggerResultsLabel, triggerResults);
01409 
01410   Int_t JetLoPass = 0;
01411   
01412   if (triggerResults.isValid()) {
01413     if (DEBUG) std::cout << "trigger valid " << std::endl;
01414     //    edm::TriggerNames triggerNames;    // TriggerNames class
01415     //    triggerNames.init(*triggerResults);
01416     unsigned int n = triggerResults->size();
01417     for (unsigned int i=0; i!=n; i++) {
01418 
01419       /***
01420       std::cout << "   Trigger Name = " << triggerNames.triggerName(i)
01421                 << " Accept = " << triggerResults->accept(i)
01422                 << std::endl;
01423       ***/
01424 
01425       //      if (DEBUG) std::cout <<  triggerNames.triggerName(i) << std::endl;
01426 
01427       /***
01428       if ( triggerNames.triggerName(i) == "HLT_Jet30" ) {
01429         JetLoPass =  triggerResults->accept(i);
01430         if (DEBUG) std::cout << "Found  HLT_Jet30 " 
01431                              << JetLoPass
01432                              << std::endl;
01433       }
01434       ***/
01435 
01436     }
01437       
01438   } else {
01439 
01440     edm::Handle<TriggerResults> *tr = new edm::Handle<TriggerResults>;
01441     triggerResults = (*tr);
01442 
01443     //     std::cout << "triggerResults is not valid" << std::endl;
01444     //     std::cout << triggerResults << std::endl;
01445     //     std::cout << triggerResults.isValid() << std::endl;
01446     
01447     if (DEBUG) std::cout << "trigger not valid " << std::endl;
01448     edm::LogInfo("myJetAna") << "TriggerResults::HLT not found, "
01449       "automatically select events";
01450     //return;
01451   }
01452 
01453   /****
01454   Handle <L1GlobalTriggerReadoutRecord> gtRecord_h;
01455   evt.getByType (gtRecord_h); // assume only one L1 trigger record here
01456   const L1GlobalTriggerReadoutRecord* gtRecord = gtRecord_h.failedToGet () ? 0 : &*gtRecord_h;
01457   
01458   if (gtRecord) { // object is available
01459     for (int l1bit = 0; l1bit < 128; ++l1bit) {
01460       if (gtRecord->decisionWord() [l1bit]) h_L1TrigBit->Fill (l1bit);
01461     }
01462   }
01463   ****/
01464 
01465 
01466 
01467 
01468   // **************************************************************
01469   // ** Loop over the two leading CaloJets and fill some histograms
01470   // **************************************************************
01471   Handle<CaloJetCollection> caloJets;
01472   evt.getByLabel( CaloJetAlgorithm, caloJets );
01473 
01474 
01475   jetInd    = 0;
01476   allJetInd = 0;
01477 
01478   EtaOk10 = 0;
01479   EtaOk13 = 0;
01480   EtaOk40 = 0;
01481 
01482   //  const JetCorrector* corrector = 
01483   //    JetCorrector::getJetCorrector (JetCorrectionService, es);
01484 
01485 
01486   highestPt = 0.0;
01487   nextPt    = 0.0;
01488   
01489   for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++ cal ) {
01490     
01491     //    double scale = corrector->correction (*cal);
01492     double scale = 1.0;
01493     double corPt = scale*cal->pt();
01494     //    double corPt = cal->pt();
01495     //    cout << "Pt = " << cal->pt() << endl;
01496     
01497     if (corPt>highestPt) {
01498       nextPt      = highestPt;
01499       p4cortmp[1] = p4cortmp[0]; 
01500       highestPt   = corPt;
01501       p4cortmp[0] = scale*cal->p4();
01502     } else if (corPt>nextPt) {
01503       nextPt      = corPt;
01504       p4cortmp[1] = scale*cal->p4();
01505     }
01506 
01507     allJetInd++;
01508     if (allJetInd == 1) {
01509       h_jet1Pt->Fill( cal->pt() );
01510       h_jet1Eta->Fill( cal->eta() );
01511       if (JetLoPass != 0) h_jet1PtHLT->Fill( cal->pt() );
01512       p4tmp[0] = cal->p4();
01513       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
01514       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
01515       if ( fabs(cal->eta()) < 4.0) EtaOk40++;            
01516     }
01517     if (allJetInd == 2) {
01518       h_jet2Pt->Fill( cal->pt() );
01519       h_jet2Eta->Fill( cal->eta() );
01520       p4tmp[1] = cal->p4();
01521       if ( fabs(cal->eta()) < 1.0) EtaOk10++;
01522       if ( fabs(cal->eta()) < 1.3) EtaOk13++;
01523       if ( fabs(cal->eta()) < 4.0) EtaOk40++;
01524     }
01525 
01526     if ( cal->pt() > minJetPt) {
01527       const std::vector<CaloTowerPtr> jetCaloRefs = cal->getCaloConstituents();
01528       int nConstituents = jetCaloRefs.size();
01529       h_nTowersCal->Fill(nConstituents);
01530       h_EMFracCal->Fill(cal->emEnergyFraction());    
01531       h_ptCal->Fill( cal->pt() );         
01532       h_etaCal->Fill( cal->eta() );
01533       h_phiCal->Fill( cal->phi() );
01534       jetInd++;
01535     }
01536   }
01537 
01538   h_nCalJets->Fill( jetInd ); 
01539 
01540   if (jetInd > 1) {
01541     LeadMass = (p4tmp[0]+p4tmp[1]).mass();
01542     dijetMass->Fill( LeadMass );    
01543   }
01544 
01545 
01546   // ******************
01547   // *** Jet Properties
01548   // ******************
01549 
01550   int nTow1, nTow2, nTow3, nTow4;
01551   //  Handle<CaloJetCollection> jets;
01552   //  evt.getByLabel( CaloJetAlgorithm, jets );
01553 
01554   // *********************************************************
01555   // --- Loop over jets and make a list of all the used towers
01556   int jjet = 0;
01557   for ( CaloJetCollection::const_iterator ijet=caloJets->begin(); ijet!=caloJets->end(); ijet++) {
01558     jjet++;
01559 
01560     float hadEne  = ijet->hadEnergyInHB() + ijet->hadEnergyInHO() + 
01561                     ijet->hadEnergyInHE() + ijet->hadEnergyInHF();                   
01562     float emEne   = ijet->emEnergyInEB() + ijet->emEnergyInEE() + ijet->emEnergyInHF();
01563     float had     = ijet->energyFractionHadronic();    
01564     float j_et = ijet->et();
01565 
01566     // *** Barrel
01567     if (fabs(ijet->eta()) < 1.3) {
01568       totEneLeadJetEta1->Fill(hadEne+emEne); 
01569       hadEneLeadJetEta1->Fill(ijet->hadEnergyInHB()); 
01570       emEneLeadJetEta1->Fill(ijet->emEnergyInEB());       
01571       if (ijet->pt() > minJetPt10) hadFracEta1->Fill(had);
01572     }
01573 
01574     // *** EndCap
01575     if ((fabs(ijet->eta()) > 1.3) && (fabs(ijet->eta()) < 3.) ) {
01576       totEneLeadJetEta2->Fill(hadEne+emEne); 
01577       hadEneLeadJetEta2->Fill(ijet->hadEnergyInHE()); 
01578       emEneLeadJetEta2->Fill(ijet->emEnergyInEE());       
01579       if (ijet->pt() > minJetPt10) hadFracEta2->Fill(had);
01580     }
01581 
01582     // *** Forward
01583     if (fabs(ijet->eta()) > 3.) {
01584       totEneLeadJetEta3->Fill(hadEne+emEne); 
01585       hadEneLeadJetEta3->Fill(hadEne); 
01586       emEneLeadJetEta3->Fill(emEne); 
01587       if (ijet->pt() > minJetPt10) hadFracEta3->Fill(had);
01588     }
01589 
01590     // *** CaloTowers in Jet
01591     const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
01592     int nConstituents = jetCaloRefs.size();
01593     NTowers->Fill(nConstituents);
01594 
01595     if (jjet == 1) {
01596 
01597       nTow1 = nTow2 = nTow3 = nTow4 = 0;
01598       for (int i = 0; i <nConstituents ; i++){
01599 
01600         float et  = jetCaloRefs[i]->et();
01601 
01602         if (et > 0.5) nTow1++;
01603         if (et > 1.0) nTow2++;
01604         if (et > 1.5) nTow3++;
01605         if (et > 2.0) nTow4++;
01606         
01607         hf_TowerJetEt->Fill(et/j_et);
01608 
01609       }
01610 
01611       nTowersLeadJetPt1->Fill(nTow1);
01612       nTowersLeadJetPt2->Fill(nTow2);
01613       nTowersLeadJetPt3->Fill(nTow3);
01614       nTowersLeadJetPt4->Fill(nTow4);
01615 
01616     }
01617 
01618   }
01619 
01620 
01621   // **********************
01622   // *** Unclustered Energy
01623   // **********************
01624 
01625   double SumPtJet(0);
01626 
01627   double SumEtNotJets(0);
01628   double SumEtJets(0);
01629   double SumEtTowers(0);
01630   double TotalClusteredE(0);
01631   double TotalUnclusteredE(0);
01632 
01633   double sumJetPx(0);
01634   double sumJetPy(0);
01635 
01636   double sumTowerAllPx(0);
01637   double sumTowerAllPy(0);
01638 
01639   double sumTowerAllEx(0);
01640   double sumTowerAllEy(0);
01641 
01642   // double HCALTotalE;
01643   double HBTotalE, HETotalE, HOTotalE, HFTotalE;
01644   // double ECALTotalE;
01645   double EBTotalE, EETotalE;
01646 
01647   std::vector<CaloTowerPtr>   UsedTowerList;
01648   std::vector<CaloTower>      TowerUsedInJets;
01649   std::vector<CaloTower>      TowerNotUsedInJets;
01650 
01651   // *********************
01652   // *** Hcal recHits
01653   // *********************
01654 
01655   edm::Handle<HcalSourcePositionData> spd;
01656 
01657   // HCALTotalE = 0.;
01658   HBTotalE = HETotalE = HOTotalE = HFTotalE = 0.;
01659   try {
01660     std::vector<edm::Handle<HBHERecHitCollection> > colls;
01661     evt.getManyByType(colls);
01662     std::vector<edm::Handle<HBHERecHitCollection> >::iterator i;
01663     for (i=colls.begin(); i!=colls.end(); i++) {
01664 
01665 
01666       for (HBHERecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
01667         //      std::cout << *j << std::endl;
01668         if (j->id().subdet() == HcalBarrel) {
01669           HBEne->Fill(j->energy()); 
01670           HBTime->Fill(j->time()); 
01671           if (!Pass_NoiseSummary) HBTimeFlagged2->Fill(j->time()); 
01672           if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HBTimeFlagged->Fill(j->time()); 
01673           HBTvsE->Fill(j->energy(), j->time());
01674 
01675           if (j->time() > 20.) HBEneTThr->Fill(j->energy()); 
01676           
01677           if ((j->time()<-25.) || (j->time()>75.)) {
01678             HBEneOOT->Fill(j->energy()); 
01679             if (j->energy() > HBHEThreshold)  HBEneOOTTh->Fill(j->energy()); 
01680             if (j->energy() > HBHEThreshold1) HBEneOOTTh1->Fill(j->energy()); 
01681           }
01682           if (j->energy() > HBHEThreshold) {
01683             HBEneTh->Fill(j->energy()); 
01684             HBTimeTh->Fill(j->time()); 
01685             if (!Pass_NoiseSummary) HBTimeThFlagged2->Fill(j->time()); 
01686             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HBTimeThFlagged->Fill(j->time()); 
01687 
01688             if (evt.id().run() >= StableRun) HBTimeThR->Fill(j->time()); 
01689             HBTotalE += j->energy();
01690             HBocc->Fill(j->id().ieta(),j->id().iphi());
01691             hitEta->Fill(j->id().ieta());
01692             hitPhi->Fill(j->id().iphi());
01693           }
01694           if (j->energy() > HBHEThreshold1) {
01695             HBEneTh1->Fill(j->energy()); 
01696             HBTimeTh1->Fill(j->time()); 
01697             if (!Pass_NoiseSummary) HBTimeTh1Flagged2->Fill(j->time()); 
01698             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HBTimeTh1Flagged->Fill(j->time()); 
01699 
01700             if (evt.id().run() >= StableRun) HBTimeTh1R->Fill(j->time()); 
01701             if ((j->time()<-25.) || (j->time()>75.)) {
01702               HBoccOOT->Fill(j->id().ieta(),j->id().iphi());
01703             }
01704           }
01705           if (j->energy() > HBHEThreshold2) {
01706             HBTimeTh2->Fill(j->time()); 
01707             if (!Pass_NoiseSummary) HBTimeTh2Flagged2->Fill(j->time()); 
01708             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HBTimeTh2Flagged->Fill(j->time()); 
01709 
01710             if (evt.id().run() >= StableRun) HBTimeTh2R->Fill(j->time()); 
01711           }
01712           if (j->energy() > HBHEThreshold3) {
01713             HBTimeTh3->Fill(j->time()); 
01714             if (evt.id().run() >= StableRun) HBTimeTh3R->Fill(j->time()); 
01715           }
01716           if ( (evt.id().run() == 120020) && (evt.id().event() == 453) ) {
01717             HBEneX->Fill(j->energy()); 
01718             if (j->energy() > HBHEThreshold) HBTimeX->Fill(j->time()); 
01719           }
01720           if ( (evt.id().run() == 120020) && (evt.id().event() == 457) ) {
01721             HBEneY->Fill(j->energy()); 
01722             if (j->energy() > HBHEThreshold) HBTimeY->Fill(j->time()); 
01723           }
01724         }
01725         if (j->id().subdet() == HcalEndcap) {
01726           HEEne->Fill(j->energy()); 
01727           HETime->Fill(j->time()); 
01728           if (!Pass_NoiseSummary) HETimeFlagged2->Fill(j->time()); 
01729           if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HETimeFlagged->Fill(j->time()); 
01730           HETvsE->Fill(j->energy(), j->time());
01731 
01732           if (j->time() > 20.) HEEneTThr->Fill(j->energy()); 
01733 
01734           if ((j->time()<-25.) || (j->time()>75.)) {
01735             HEEneOOT->Fill(j->energy()); 
01736             if (j->energy() > HBHEThreshold)  HEEneOOTTh->Fill(j->energy());  
01737             if (j->energy() > HBHEThreshold1) HEEneOOTTh1->Fill(j->energy());  
01738           }
01739 
01740           if (j->energy() > HBHEThreshold) {
01741             HEEneTh->Fill(j->energy()); 
01742             HETimeTh->Fill(j->time()); 
01743             if (!Pass_NoiseSummary) HETimeThFlagged2->Fill(j->time()); 
01744             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0)  HETimeThFlagged->Fill(j->time()); 
01745 
01746             if (evt.id().run() >= StableRun) HETimeThR->Fill(j->time()); 
01747             HETotalE += j->energy();
01748             HEocc->Fill(j->id().ieta(),j->id().iphi());
01749             hitEta->Fill(j->id().ieta());
01750             hitPhi->Fill(j->id().iphi());
01751           }
01752           if (j->energy() > HBHEThreshold1) {
01753             HEEneTh1->Fill(j->energy()); 
01754             HETimeTh1->Fill(j->time()); 
01755             if (!Pass_NoiseSummary) HETimeTh1Flagged2->Fill(j->time()); 
01756             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HETimeTh1Flagged->Fill(j->time()); 
01757             if (evt.id().run() >= StableRun) HETimeTh1R->Fill(j->time()); 
01758             if ((j->time()<-25.) || (j->time()>75.)) {
01759               HEoccOOT->Fill(j->id().ieta(),j->id().iphi());
01760             }
01761           }
01762           if (j->energy() > HBHEThreshold2) {
01763             HETimeTh2->Fill(j->time()); 
01764             if (!Pass_NoiseSummary) HETimeTh2Flagged2->Fill(j->time()); 
01765             if (j->flagField(HcalCaloFlagLabels::HBHETimingShapedCutsBits) != 0) HETimeTh2Flagged->Fill(j->time()); 
01766             if (evt.id().run() >= StableRun) HETimeTh2R->Fill(j->time()); 
01767           }
01768           if (j->energy() > HBHEThreshold3) {
01769             HETimeTh3->Fill(j->time()); 
01770             if (evt.id().run() >= StableRun) HETimeTh3R->Fill(j->time()); 
01771           }
01772 
01773           if ( (evt.id().run() == 120020) && (evt.id().event() == 453) ) {
01774             HEEneX->Fill(j->energy()); 
01775             if (j->energy() > HBHEThreshold) HETimeX->Fill(j->time()); 
01776           }
01777           if ( (evt.id().run() == 120020) && (evt.id().event() == 457) ) {
01778             HEEneY->Fill(j->energy()); 
01779             if (j->energy() > HBHEThreshold) HETimeY->Fill(j->time()); 
01780           }
01781 
01782           // Fill +-HE separately
01783           if (j->id().ieta()<0) {
01784             HEnegEne->Fill(j->energy()); 
01785             if (j->energy() > HBHEThreshold) {
01786               HEnegTime->Fill(j->time()); 
01787             }
01788           } else {
01789             HEposEne->Fill(j->energy()); 
01790             if (j->energy() > HBHEThreshold) {
01791               HEposTime->Fill(j->time()); 
01792             }
01793           }
01794           
01795         }
01796 
01797         /***
01798         std::cout << j->id()     << " "
01799                   << j->id().subdet() << " "
01800                   << j->id().ieta()   << " "
01801                   << j->id().iphi()   << " "
01802                   << j->id().depth()  << " "
01803                   << j->energy() << " "
01804                   << j->time()   << std::endl;
01805         ****/
01806       }
01807     }
01808   } catch (...) {
01809     cout << "No HB/HE RecHits." << endl;
01810   }
01811 
01812 
01813   HFM_ETime = 0.;
01814   HFM_E = 0.;
01815   HFP_ETime = 0.;
01816   HFP_E = 0.;
01817 
01818   int NPMTHits;
01819   NPMTHits = 0;
01820   try {
01821     std::vector<edm::Handle<HFRecHitCollection> > colls;
01822     evt.getManyByType(colls);
01823     std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
01824     for (i=colls.begin(); i!=colls.end(); i++) {
01825       for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
01826         if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) ||
01827              (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01828              NPMTHits++;
01829         }
01830       }
01831       break;
01832     }
01833   } catch (...) {
01834     cout << "No HF RecHits." << endl;
01835   }
01836 
01837 
01838   PMTHits->Fill(NPMTHits); 
01839 
01840   try {
01841     std::vector<edm::Handle<HFRecHitCollection> > colls;
01842     evt.getManyByType(colls);
01843     std::vector<edm::Handle<HFRecHitCollection> >::iterator i;
01844     for (i=colls.begin(); i!=colls.end(); i++) {
01845       for (HFRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
01846 
01847         /****
01848         float en = j->energy();
01849         HcalDetId id(j->detid().rawId());
01850         int ieta = id.ieta();
01851         int iphi = id.iphi();
01852         int depth = id.depth();
01853         *****/
01854 
01855         //  std::cout << *j << std::endl;
01856 
01857         if (j->id().subdet() == HcalForward) {
01858 
01859           if (NPMTHits == 1) {
01860             if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) ||
01861                  (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01862               HFEtaFlagged->Fill(j->id().ieta());
01863               if (j->id().depth() == 1) HFEtaFlaggedL->Fill(j->id().ieta());
01864               if (j->id().depth() == 2) HFEtaFlaggedS->Fill(j->id().ieta());
01865             } else {
01866               HFEtaNFlagged->Fill(j->id().ieta(), j->energy());
01867               HFEtaPhiNFlagged->Fill(j->id().ieta(),j->id().iphi(),j->energy());
01868             }
01869           }
01870           if (j->energy() > 20.) {
01871             if (NPMTHits == 0) {
01872               HFEnePMT0->Fill(j->energy()); 
01873               HFTimePMT0->Fill(j->time()); 
01874             }
01875             if (NPMTHits == 1) {
01876               if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) ||
01877                    (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01878                 HFEnePMT1->Fill(j->energy()); 
01879                 HFTimePMT1->Fill(j->time()); 
01880               }
01881             }
01882             if (NPMTHits > 1) {
01883               if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) ||
01884                    (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01885                 HFEnePMT2->Fill(j->energy()); 
01886                 HFTimePMT2->Fill(j->time()); 
01887               }
01888             }
01889           }
01890 
01891           HFTimeVsiEtaP->Fill(j->id().ieta(), j->time());
01892           HFTimeVsiEtaM->Fill(j->id().ieta(), j->time());
01893 
01894           if (j->energy() > 5.) { 
01895             HFTimeVsiEtaP5->Fill(j->id().ieta(), j->time());
01896             HFTimeVsiEtaM5->Fill(j->id().ieta(), j->time());
01897           }       
01898 
01899           if (j->energy() > 20.) { 
01900             HFTimeVsiEtaP20->Fill(j->id().ieta(), j->time());
01901             HFTimeVsiEtaM20->Fill(j->id().ieta(), j->time());
01902           }       
01903 
01904           HFEne->Fill(j->energy()); 
01905           HFTime->Fill(j->time()); 
01906           HFTvsE->Fill(j->energy(), j->time());
01907 
01908           if (j->time() > 20.) HFEneTThr->Fill(j->energy()); 
01909          
01910           if (j->energy() > 10.) HFTvsEThr->Fill(j->energy(), j->time());
01911 
01912           if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1)|| 
01913                (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01914             HFEneFlagged->Fill(j->energy());
01915             HFoccFlagged->Fill(j->id().ieta(),j->id().iphi());
01916             HFTimeFlagged->Fill(j->time()); 
01917             HFTvsEFlagged->Fill(j->energy(), j->time());
01918 
01919             //      std::cout << "Flagged:  " << j->energy() << " "
01920             //                << j->time()
01921             //                << std::endl;
01922           }
01923 
01924 
01925           if (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) {
01926             HFEneFlagged2->Fill(j->energy());
01927             HFoccFlagged2->Fill(j->id().ieta(),j->id().iphi());
01928             HFTimeFlagged2->Fill(j->time()); 
01929             HFTvsEFlagged2->Fill(j->energy(), j->time());
01930             if (j->energy() > 10.) HFTvsEFlagged2Thr->Fill(j->energy(), j->time());
01931           }
01932 
01933           if (j->flagField(HcalCaloFlagLabels::HFDigiTime) == 1) {
01934             HFTimeFlagged3->Fill(j->time()); 
01935           }
01936 
01937           if (j->energy() > HFThreshold) {
01938             HFEneTh->Fill(j->energy()); 
01939             HFTimeTh->Fill(j->time()); 
01940             if (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) HFTimeThFlagged2->Fill(j->time()); 
01941             if (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) HFTimeThFlagged3->Fill(j->time()); 
01942 
01943             if (evt.id().run() >= StableRun) HFTimeThR->Fill(j->time()); 
01944             if ( (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1)|| 
01945                  (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) ) {
01946 
01947               HFTimeThFlagged->Fill(j->time()); 
01948 
01949               if (j->energy() > HFThreshold2) HFTimeTh2Flagged->Fill(j->time()); 
01950               if (j->energy() > HFThreshold3) HFTimeTh3Flagged->Fill(j->time()); 
01951 
01952               if (evt.id().run() >= StableRun) {
01953                 HFTimeThFlaggedR->Fill(j->time()); 
01954                 if (NPMTHits == 1) HFTimeThFlaggedR1->Fill(j->time()); 
01955                 if (NPMTHits == 2) HFTimeThFlaggedR2->Fill(j->time()); 
01956                 if (NPMTHits == 3) HFTimeThFlaggedR3->Fill(j->time()); 
01957                 if (NPMTHits == 4) HFTimeThFlaggedR4->Fill(j->time()); 
01958                 if (NPMTHits > 1) HFTimeThFlaggedRM->Fill(j->time()); 
01959               }
01960             }
01961             HFTotalE += j->energy();
01962             HFocc->Fill(j->id().ieta(),j->id().iphi());
01963             hitEta->Fill(j->id().ieta());
01964             hitPhi->Fill(j->id().iphi());
01965           }
01966 
01967           if (j->energy() > HFThreshold1) {
01968             HFEneTh1->Fill(j->energy());
01969             HFTimeTh1->Fill(j->time()); 
01970             if (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) HFTimeTh1Flagged2->Fill(j->time()); 
01971             if (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) HFTimeTh1Flagged3->Fill(j->time()); 
01972             if (evt.id().run() >= StableRun) HFTimeTh1R->Fill(j->time()); 
01973             if ((j->time()<-20.) || (j->time()>20.)) {
01974               HFoccOOT->Fill(j->id().ieta(),j->id().iphi());
01975             }
01976           } 
01977           if (j->energy() > HFThreshold2) {
01978             HFTimeTh2->Fill(j->time()); 
01979             if (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) HFTimeTh2Flagged2->Fill(j->time()); 
01980             if (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) HFTimeTh2Flagged3->Fill(j->time()); 
01981             if (evt.id().run() >= StableRun) HFTimeTh2R->Fill(j->time()); 
01982           }
01983           if (j->energy() > HFThreshold3) {
01984             HFTimeTh3->Fill(j->time()); 
01985             if (j->flagField(HcalCaloFlagLabels::HFLongShort) == 1) HFTimeTh3Flagged2->Fill(j->time()); 
01986             if (j->flagField(HcalCaloFlagLabels::HFDigiTime)  == 1) HFTimeTh3Flagged3->Fill(j->time()); 
01987             if (evt.id().run() >= StableRun) HFTimeTh3R->Fill(j->time()); 
01988           }
01989 
01990           if (j->id().ieta()<0) {
01991             if (j->energy() > HFThreshold) {
01992               //              HFTimeM->Fill(j->time()); 
01993               HFEneM->Fill(j->energy()); 
01994               HFM_ETime += j->energy()*j->time(); 
01995               HFM_E     += j->energy();
01996             }
01997           } else {
01998             if (j->energy() > HFThreshold) {
01999               //              HFTimeP->Fill(j->time()); 
02000               HFEneP->Fill(j->energy()); 
02001               HFP_ETime += j->energy()*j->time(); 
02002               HFP_E     += j->energy();
02003             }
02004           }
02005 
02006           // Long and short fibers
02007           if (j->id().depth() == 1){
02008             HFLEne->Fill(j->energy()); 
02009             if (j->energy() > HFThreshold) HFLTime->Fill(j->time());
02010           } else {
02011             HFSEne->Fill(j->energy()); 
02012             if (j->energy() > HFThreshold) HFSTime->Fill(j->time());
02013           }
02014         }
02015       }
02016       break;
02017 
02018     }
02019 
02020   } catch (...) {
02021     cout << "No HF RecHits." << endl;
02022   }
02023 
02024 
02025 
02026   for (int ieta=0; ieta<100; ieta++) {
02027      for (int iphi=0; iphi<100; iphi++) {
02028        double longF, shortF;
02029        if (HFRecHit[ieta][iphi][0] == -10.) {
02030          longF = 0.;
02031        } else {
02032          longF = HFRecHit[ieta][iphi][0];
02033        }
02034        if (HFRecHit[ieta][iphi][1] == -10.) {
02035          shortF = 0.;
02036        } else {
02037          shortF = HFRecHit[ieta][iphi][1];
02038        }
02039        //       if ((longF > HFThreshold) || (shortF > HFThreshold)) HFLSRatio->Fill((longF-shortF)/(longF+shortF));
02040 
02041        if (longF > 0.) HFLEneAll->Fill(longF);
02042        if (shortF > 0.) HFSEneAll->Fill(shortF);
02043 
02044 
02045        if ((longF > 20.) || (shortF > 20.)) {
02046          double R = (longF-shortF)/(longF+shortF);
02047          HFLSRatio->Fill(R);
02048          if (fabs(R) > 0.995) {    
02049 
02050            //      if (longF > 110.)  {
02051            //      if (longF > 50.)  {
02052            if (longF > (162.4-10.19*abs(ieta-41)+.21*abs(ieta-41)*abs(ieta-41)) )  {
02053              HFEtaFlaggedLN->Fill(ieta-41);
02054 
02055              HFLEneAllF->Fill(longF);
02056 
02057              if (shortF == 0.) HFLEneNoSFlaggedN->Fill(longF);
02058            }
02059            //      if (shortF > 70.)  {
02060            //      if (shortF > 50.)  {
02061            if (shortF > (129.9-6.61*abs(ieta-41)+0.1153*abs(ieta-41)*abs(ieta-41)) ) {
02062              HFEtaFlaggedSN->Fill(ieta-41);
02063 
02064              HFSEneAllF->Fill(shortF);
02065 
02066              if (longF == 0.) HFSEneNoLFlaggedN->Fill(shortF);
02067            }
02068          }
02069        }
02070        /***
02071        cout << "HF LS Ratio long= " 
02072             << longF
02073             << " short= "
02074             << shortF
02075             << endl;
02076        ***/
02077 
02078        HFLvsS->Fill(HFRecHit[ieta][iphi][1], HFRecHit[ieta][iphi][0]);         
02079        if ( (HFRecHit[ieta][iphi][1] == -10.) && (HFRecHit[ieta][iphi][0] != -10.) ) {
02080          HFLEneNoS->Fill(HFRecHit[ieta][iphi][0]);
02081          if (HFRecHitFlag[ieta][iphi][0] !=0 ) HFLEneNoSFlagged->Fill(HFRecHit[ieta][iphi][0]);
02082        }
02083        if ( (HFRecHit[ieta][iphi][0] == -10.) && (HFRecHit[ieta][iphi][1] != -10.) ) {
02084          HFSEneNoL->Fill(HFRecHit[ieta][iphi][1]);
02085          if (HFRecHitFlag[ieta][iphi][1] !=0 ) HFSEneNoLFlagged->Fill(HFRecHit[ieta][iphi][1]);
02086        }
02087 
02088      }
02089   }
02090 
02091   if (HFP_E > 0.) HFTimeP->Fill(HFP_ETime / HFP_E);
02092   if (HFM_E > 0.) HFTimeM->Fill(HFM_ETime / HFM_E);
02093 
02094   if ((HFP_E > 0.) && (HFM_E > 0.)) {
02095     HF_PMM = (HFP_ETime / HFP_E) - (HFM_ETime / HFM_E);
02096     HFTimePM->Fill(HF_PMM); 
02097   } else {
02098     HF_PMM = INVALID;
02099   }
02100 
02101 
02102 
02103   try {
02104     std::vector<edm::Handle<HORecHitCollection> > colls;
02105     evt.getManyByType(colls);
02106     std::vector<edm::Handle<HORecHitCollection> >::iterator i;
02107     for (i=colls.begin(); i!=colls.end(); i++) {
02108       for (HORecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
02109         if (j->id().subdet() == HcalOuter) {
02110           HOEne->Fill(j->energy()); 
02111           HOTime->Fill(j->time());
02112           HOTvsE->Fill(j->energy(), j->time());
02113           if (j->energy() > HOThreshold1) {
02114             HOEneTh1->Fill(j->energy()); 
02115           }
02116           if (j->energy() > HOThreshold) {
02117             HOEneTh->Fill(j->energy()); 
02118             HOTimeTh->Fill(j->time());
02119             HOTotalE += j->energy();
02120             HOocc->Fill(j->id().ieta(),j->id().iphi());
02121           }
02122 
02123           // Separate SiPMs and HPDs:
02124           if (((j->id().iphi()>=59 && j->id().iphi()<=70 && 
02125                 j->id().ieta()>=11 && j->id().ieta()<=15) || 
02126                (j->id().iphi()>=47 && j->id().iphi()<=58 && 
02127                 j->id().ieta()>=5 && j->id().ieta()<=10)))
02128           {  
02129             HOSEne->Fill(j->energy());
02130             if (j->energy() > HOThreshold) HOSTime->Fill(j->time());
02131           } else if ((j->id().iphi()<59 || j->id().iphi()>70 || 
02132                       j->id().ieta()<11 || j->id().ieta()>15) && 
02133                      (j->id().iphi()<47 || j->id().iphi()>58 ||
02134                       j->id().ieta()<5  || j->id().ieta()>10))
02135           {
02136             HOHEne->Fill(j->energy());
02137             if (j->energy() > HOThreshold) HOHTime->Fill(j->time());
02138             // Separate rings -1,-2,0,1,2 in HPDs:
02139             if (j->id().ieta()<= -11){
02140               HOHrm2Ene->Fill(j->energy());
02141               if (j->energy() > HOThreshold) HOHrm2Time->Fill(j->time());
02142             } else if (j->id().ieta()>= -10 && j->id().ieta() <= -5) {
02143               HOHrm1Ene->Fill(j->energy());
02144               if (j->energy() > HOThreshold) HOHrm1Time->Fill(j->time());
02145             } else if (j->id().ieta()>= -4 && j->id().ieta() <= 4) {
02146               HOHr0Ene->Fill(j->energy());
02147               if (j->energy() > HOThreshold) HOHr0Time->Fill(j->time());
02148             } else if (j->id().ieta()>= 5 && j->id().ieta() <= 10) {
02149               HOHrp1Ene->Fill(j->energy());
02150               if (j->energy() > HOThreshold) HOHrp1Time->Fill(j->time());
02151             } else if (j->id().ieta()>= 11) {
02152               HOHrp2Ene->Fill(j->energy());
02153               if (j->energy() > HOThreshold) HOHrp2Time->Fill(j->time());
02154             } else {
02155               std::cout << "Finding events that are in no ring !?!" << std::endl;
02156               std::cout << "eta = " << j->id().ieta() << std::endl;
02157               
02158             }
02159           } else {
02160             std::cout << "Finding events that are neither SiPM nor HPD!?" << std::endl;     
02161           }
02162 
02163           
02164 
02165         }
02166         //      std::cout << *j << std::endl;
02167       }
02168     }
02169   } catch (...) {
02170     cout << "No HO RecHits." << endl;
02171   }
02172 
02173   // HCALTotalE = HBTotalE + HETotalE + HFTotalE + HOTotalE;
02174   // ECALTotalE = 0.;
02175   EBTotalE = EETotalE = 0.;
02176 
02177 
02178   try {
02179     std::vector<edm::Handle<EcalRecHitCollection> > colls;
02180     evt.getManyByType(colls);
02181     std::vector<edm::Handle<EcalRecHitCollection> >::iterator i;
02182     for (i=colls.begin(); i!=colls.end(); i++) {
02183       for (EcalRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
02184         if (j->id().subdetId() == EcalBarrel) {
02185           EBEne->Fill(j->energy()); 
02186           EBTime->Fill(j->time()); 
02187           if (j->energy() > EBEEThreshold) {
02188             EBEneTh->Fill(j->energy()); 
02189             EBTimeTh->Fill(j->time()); 
02190           }
02191           if ( (evt.id().run() == 120020) && (evt.id().event() == 453) ) {
02192             EBEneX->Fill(j->energy()); 
02193             EBTimeX->Fill(j->time()); 
02194           }
02195           if ( (evt.id().run() == 120020) && (evt.id().event() == 457) ) {
02196             EBEneY->Fill(j->energy()); 
02197             EBTimeY->Fill(j->time()); 
02198           }
02199           EBTotalE += j->energy();
02200         }
02201         if (j->id().subdetId() == EcalEndcap) {
02202           EEEne->Fill(j->energy()); 
02203           EETime->Fill(j->time());
02204           if (j->energy() > EBEEThreshold) {
02205             EEEneTh->Fill(j->energy()); 
02206             EETimeTh->Fill(j->time()); 
02207           }
02208           if ( (evt.id().run() == 120020) && (evt.id().event() == 453) ) {
02209             EEEneX->Fill(j->energy()); 
02210             EETimeX->Fill(j->time()); 
02211           }
02212           if ( (evt.id().run() == 120020) && (evt.id().event() == 457 ) ) {
02213             EEEneY->Fill(j->energy()); 
02214             EETimeY->Fill(j->time()); 
02215           }
02216           EETotalE += j->energy();
02217         }
02218         //      std::cout << *j << std::endl;
02219         //      std::cout << "EB ID = " << j->id().subdetId() << "/" << EcalBarrel << std::endl;
02220       }
02221     }
02222   } catch (...) {
02223     cout << "No ECAL RecHits." << endl;
02224   }
02225 
02226   EBvHB->Fill(HBTotalE, EBTotalE);
02227   EEvHE->Fill(HETotalE, EETotalE);
02228 
02229   /*****
02230   try {
02231     std::vector<edm::Handle<EBRecHitCollection> > colls;
02232     evt.getManyByType(colls);
02233     std::vector<edm::Handle<EBRecHitCollection> >::iterator i;
02234 
02235     for (i=colls.begin(); i!=colls.end(); i++) {
02236       for (EBRecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
02237         //      if (j->id().subdetId() == EcalBarrel) {
02238           EBEne->Fill(j->energy()); 
02239           EBTime->Fill(j->time()); 
02240           //      EBTotalE = j->energy();
02241           //    }
02242           //    std::cout << *j << std::endl;
02243           //    std::cout << "EB ID = " << j->id().subdetId() << "/" << EcalBarrel << std::endl;
02244       }
02245     }
02246   } catch (...) {
02247     cout << "No EB RecHits." << endl;
02248   }
02249 
02250   try {
02251     std::vector<edm::Handle<EERecHitCollection> > colls;
02252     evt.getManyByType(colls);
02253     std::vector<edm::Handle<EERecHitCollection> >::iterator i;
02254     for (i=colls.begin(); i!=colls.end(); i++) {
02255       for (EERecHitCollection::const_iterator j=(*i)->begin(); j!=(*i)->end(); j++) {
02256         //      if (j->id().subdetId() == EcalEndcap) {
02257           EEEne->Fill(j->energy()); 
02258           EETime->Fill(j->time());
02259           //      EETotalE = j->energy();
02260           // Separate +-EE;
02261           EEDetId EEid = EEDetId(j->id());
02262           if (!EEid.positiveZ()) 
02263           {
02264             EEnegEne->Fill(j->energy()); 
02265             EEnegTime->Fill(j->time()); 
02266           }else{
02267             EEposEne->Fill(j->energy()); 
02268             EEposTime->Fill(j->time()); 
02269           }
02270           //    }
02271         //      std::cout << *j << std::endl;
02272       }
02273     }
02274   } catch (...) {
02275     cout << "No EE RecHits." << endl;
02276   }
02277   ******/
02278 
02279   // ECALTotalE = EBTotalE + EETotalE;
02280 
02281   if ( (EBTotalE > 320000)  && (EBTotalE < 330000) && 
02282        (HBTotalE > 2700000) && (HBTotalE < 2800000) ) {
02283 
02284     std::cout << ">>> Off Axis! " 
02285               << std::endl;
02286     
02287   }
02288 
02289   /***
02290   std::cout << " Rechits: Total Energy :  " 
02291             << " HCAL= " << HCALTotalE 
02292             << " ECAL= " << ECALTotalE
02293             << " HB = " << HBTotalE
02294             << " EB = " << EBTotalE
02295             << std::endl;
02296   ***/
02297 
02298 
02299   // *********************
02300   // *** CaloTowers
02301   // *********************
02302   //  Handle<CaloTowerCollection> caloTowers;
02303   //  evt.getByLabel( "towerMaker", caloTowers );
02304 
02305   nTow1 = nTow2 = nTow3 = nTow4 = 0;
02306 
02307   double sum_et = 0.0;
02308   double sum_ex = 0.0;
02309   double sum_ey = 0.0;
02310 
02311   double HFsum_et = 0.0;
02312   double HFsum_ex = 0.0;
02313   double HFsum_ey = 0.0;
02314   //  double sum_ez = 0.0;
02315 
02316 
02317   //  std::cout<<">>>> Run " << evt.id().run() << " Event " << evt.id().event() << std::endl;
02318   // --- Loop over towers and make a lists of used and unused towers
02319   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02320        tower != caloTowers->end(); tower++) {
02321 
02322     Double_t  et   = tower->et();
02323     Double_t  phix = tower->phi();
02324     
02325     if (et > 0.5) nTow1++;
02326     if (et > 1.0) nTow2++;
02327     if (et > 1.5) nTow3++;
02328     if (et > 2.0) nTow4++;
02329 
02330     //    if ( (fabs(tower->ieta() > 42)) ||  (fabs(tower->iphi()) > 72) ) {
02331     //      std::cout << "ieta/iphi = " <<  tower->ieta() << " / "  << tower->iphi() << std::endl;
02332     //    }
02333 
02334     if (tower->emEnergy() > 2.0) {
02335       h_EmEnergy->Fill (tower->ieta(), tower->iphi(), tower->emEnergy());
02336     }
02337     if (tower->hadEnergy() > 2.0) {
02338       h_HadEnergy->Fill (tower->ieta(), tower->iphi(), tower->hadEnergy());
02339     }
02340 
02341     if (fabs(tower->ieta()) > 29) {
02342       HFsum_et += et;
02343       HFsum_ex += et*cos(phix);
02344       HFsum_ey += et*sin(phix);
02345     }
02346 
02347 
02348     if (et>0.5) {
02349 
02350       ETime->Fill(tower->ecalTime());
02351       HTime->Fill(tower->hcalTime());
02352 
02353       // ********
02354       //      double theta = tower->theta();
02355       //      double e     = tower->energy();
02356       //      double et    = e*sin(theta);
02357       //      double et    = tower->emEt() + tower->hadEt();
02358       //      sum_ez += e*cos(theta);
02359       sum_et += et;
02360       sum_ex += et*cos(phix);
02361       sum_ey += et*sin(phix);
02362       // ********
02363 
02364       Double_t phi = tower->phi();
02365       SumEtTowers += tower->et();
02366 
02367       sumTowerAllEx += et*cos(phi);
02368       sumTowerAllEy += et*sin(phi);
02369 
02370     }
02371 
02372   }
02373 
02374   //  SumEt->Fill(sum_et);
02375   //  MET->Fill(sqrt( sum_ex*sum_ex + sum_ey*sum_ey));
02376 
02377   HFSumEt->Fill(HFsum_et);
02378   HFMET->Fill(sqrt( HFsum_ex*HFsum_ex + HFsum_ey*HFsum_ey));
02379 
02380   hf_sumTowerAllEx->Fill(sumTowerAllEx);
02381   hf_sumTowerAllEy->Fill(sumTowerAllEy);
02382 
02383   nTowers1->Fill(nTow1);
02384   nTowers2->Fill(nTow2);
02385   nTowers3->Fill(nTow3);
02386   nTowers4->Fill(nTow4);
02387 
02388 
02389   // *********************
02390   // *********************
02391 
02392   UsedTowerList.clear();
02393   TowerUsedInJets.clear();
02394   TowerNotUsedInJets.clear();
02395 
02396   // --- Loop over jets and make a list of all the used towers
02397   //  evt.getByLabel( CaloJetAlgorithm, jets );
02398   for ( CaloJetCollection::const_iterator ijet=caloJets->begin(); ijet!=caloJets->end(); ijet++) {
02399 
02400     Double_t jetPt  = ijet->pt();
02401     Double_t jetEta = ijet->eta();
02402     Double_t jetPhi = ijet->phi();
02403 
02404     //    if (jetPt>5.0) {
02405 
02406       Double_t jetPx = jetPt*cos(jetPhi);
02407       Double_t jetPy = jetPt*sin(jetPhi);
02408 
02409       sumJetPx +=jetPx;
02410       sumJetPy +=jetPy;
02411 
02412       const std::vector<CaloTowerPtr> jetCaloRefs = ijet->getCaloConstituents();
02413       int nConstituents = jetCaloRefs.size();
02414       for (int i = 0; i <nConstituents ; i++){
02415         
02416         UsedTowerList.push_back(jetCaloRefs[i]);
02417       }
02418       
02419       SumPtJet +=jetPt;
02420 
02421     //    }
02422 
02423       if ( (jetPt>80.0) && (fabs(jetEta) < 1.3) ){
02424         st_Pt->Fill( jetPt );
02425         int nConstituents = ijet->getCaloConstituents().size();
02426         st_Constituents->Fill( nConstituents );
02427         
02428         float maxEne = 0.;
02429         float totEne = 0.;
02430           
02431         for(unsigned twr=0; twr<ijet->getCaloConstituents().size(); ++twr){
02432           CaloTowerPtr tower = (ijet->getCaloConstituents())[twr];
02433           //      CaloTowerDetId id = tower->id();     
02434           if( tower->et()>0. ){
02435 
02436             if (tower->energy() > maxEne) maxEne = tower->energy();
02437             totEne += tower->energy();
02438 
02439             st_Energy->Fill( tower->energy() );
02440             st_EmEnergy->Fill( tower->emEnergy() );
02441             st_HadEnergy->Fill( tower->hadEnergy() );
02442             st_OuterEnergy->Fill( tower->outerEnergy() );
02443 
02444             st_Eta->Fill( tower->eta() );
02445             st_Phi->Fill( tower->phi() );
02446 
02447             st_iEta->Fill( tower->ieta() );
02448             st_iPhi->Fill( tower->iphi() );
02449 
02450             /****
02451             std::cout << ">>> Towers :  " 
02452                       << " " << tower->energy() 
02453                       << " " << tower->emEnergy()
02454                       << " " << tower->hadEnergy()
02455                       << " " << tower->outerEnergy()
02456                       << " " << tower->et()
02457                       << " " << tower->emEt() 
02458                       << " " << tower->hadEt() 
02459                       << " " << tower->outerEt() 
02460                       << " " << tower->eta() 
02461                       << " " << tower->phi()        
02462                       << std::endl;
02463             ****/
02464           }
02465         }
02466         st_Frac->Fill( maxEne / totEne );
02467 
02468       }
02469 
02470   }
02471 
02472   int NTowersUsed = UsedTowerList.size();
02473 
02474   // --- Loop over towers and make a lists of used and unused towers
02475   for (CaloTowerCollection::const_iterator tower = caloTowers->begin();
02476        tower != caloTowers->end(); tower++) {
02477 
02478     CaloTower  t = *tower;
02479     Double_t  et = tower->et();
02480 
02481     if(et>0) {
02482 
02483       Double_t phi = tower->phi();
02484       SumEtTowers += tower->et();
02485 
02486       sumTowerAllPx += et*cos(phi);
02487       sumTowerAllPy += et*sin(phi);
02488 
02489       bool used = false;
02490 
02491       for(int i=0; i<NTowersUsed; i++){
02492         if(tower->id() == UsedTowerList[i]->id()){
02493           used=true;
02494           break;
02495         }
02496       }
02497 
02498       if (used) {
02499         TowerUsedInJets.push_back(t);
02500       } else {
02501         TowerNotUsedInJets.push_back(t);
02502       }
02503     }
02504   }
02505 
02506   int nUsed    = TowerUsedInJets.size();
02507   int nNotUsed = TowerNotUsedInJets.size();
02508 
02509   SumEtJets    = 0;
02510   SumEtNotJets = 0;
02511   TotalClusteredE   = 0;
02512   TotalUnclusteredE = 0;
02513 
02514   for(int i=0;i<nUsed;i++){
02515     SumEtJets += TowerUsedInJets[i].et();
02516     h_ClusteredE->Fill(TowerUsedInJets[i].energy());
02517     if (TowerUsedInJets[i].energy() > 1.0) 
02518       TotalClusteredE += TowerUsedInJets[i].energy();
02519   }
02520   h_jetEt->Fill(SumEtJets);
02521 
02522   for(int i=0;i<nNotUsed;i++){
02523     if (TowerNotUsedInJets[i].et() > 0.5)
02524       SumEtNotJets += TowerNotUsedInJets[i].et();
02525     h_UnclusteredEt->Fill(TowerNotUsedInJets[i].et());
02526     h_UnclusteredEts->Fill(TowerNotUsedInJets[i].et());
02527     h_UnclusteredE->Fill(TowerNotUsedInJets[i].energy());
02528     if (TowerNotUsedInJets[i].energy() > 1.0)  
02529       TotalUnclusteredE += TowerNotUsedInJets[i].energy();
02530   }
02531 
02532   h_TotalClusteredE->Fill(TotalClusteredE);
02533   h_TotalUnclusteredE->Fill(TotalUnclusteredE);
02534   h_TotalUnclusteredEt->Fill(SumEtNotJets);
02535 
02536   // ********************************
02537   // *** CaloMET
02538   // ********************************
02539 
02540   edm::Handle<reco::CaloMETCollection> calometcoll;
02541   evt.getByLabel("met", calometcoll);
02542   if (calometcoll.isValid()) {
02543     const CaloMETCollection *calometcol = calometcoll.product();
02544     const CaloMET *calomet;
02545     calomet = &(calometcol->front());
02546 
02547     double caloSumET  = calomet->sumEt();
02548     double caloMET    = calomet->pt();
02549     double caloMETSig = calomet->mEtSig();
02550     double caloMEx    = calomet->px();
02551     double caloMEy    = calomet->py();
02552     double caloMETPhi = calomet->phi();
02553 
02554     SumEt->Fill(caloSumET);
02555     MET->Fill(caloMET);
02556     if (std::abs(OER) > 0.8) OERMET->Fill(caloMET);
02557 
02558     if (evtType == 0) MET_Tower->Fill(caloMET);
02559     if (evtType == 1) MET_RBX->Fill(caloMET);
02560     if (evtType == 2) MET_HPD->Fill(caloMET);
02561     METSig->Fill(caloMETSig);
02562     MEx->Fill(caloMEx);
02563     MEy->Fill(caloMEy);
02564     METPhi->Fill(caloMETPhi);
02565 
02566     /***
02567     double caloEz     = calomet->e_longitudinal();
02568 
02569     double caloMaxEtInEMTowers    = calomet->maxEtInEmTowers();
02570     double caloMaxEtInHadTowers   = calomet->maxEtInHadTowers();
02571     double caloEtFractionHadronic = calomet->etFractionHadronic();
02572     double caloEmEtFraction       = calomet->emEtFraction();
02573 
02574     double caloHadEtInHB = calomet->hadEtInHB();
02575     double caloHadEtInHO = calomet->hadEtInHO();
02576     double caloHadEtInHE = calomet->hadEtInHE();
02577     double caloHadEtInHF = calomet->hadEtInHF();
02578     double caloEmEtInEB  = calomet->emEtInEB();
02579     double caloEmEtInEE  = calomet->emEtInEE();
02580     double caloEmEtInHF  = calomet->emEtInHF();
02581     ****/
02582   }
02583 
02584   // ********************************
02585   // *** Vertex
02586   // ********************************
02587   VTX  = INVALID;
02588   nVTX = 0;
02589 
02590   edm::Handle<reco::VertexCollection> vertexCollection;
02591   evt.getByLabel("offlinePrimaryVertices", vertexCollection);
02592   const reco::VertexCollection vC = *(vertexCollection.product());
02593 
02594   //  std::cout << "Reconstructed "<< vC.size() << " vertices" << std::endl ;
02595   nVTX = vC.size();
02596   //double vertex_numTrks;
02597   for (reco::VertexCollection::const_iterator vertex=vC.begin(); vertex!=vC.end(); vertex++){
02598 
02599     h_Vx->Fill(vertex->x());
02600     h_Vy->Fill(vertex->y());
02601     h_Vz->Fill(vertex->z());
02602     VTX  = vertex->z();
02603     //    vertex_numTrks = vertex->tracksSize();
02604     //    h_VNTrks->Fill(vertex_numTrks);
02605 
02606   }
02607 
02608   if ((HF_PMM != INVALID) || (nVTX > 0)) {
02609     HFvsZ->Fill(HF_PMM,VTX);
02610   }
02611 
02612   // ********************************
02613   // *** Pixel Clusters
02614   // ********************************
02615   //  edm::Handle< edmNew::DetSetVector<SiPixelCluster> > hClusterColl;
02616   //  evt.getByLabel("siPixelClusters", hClusterColl);
02617   //  const edmNew::DetSetVector<SiPixelCluster> clustColl = *(hClusterColl.product());
02618 
02619   SiClusters->Fill(clustColl.size());
02620 
02621   // ********************************
02622   // *** Tracks
02623   // ********************************
02624   //  edm::Handle<reco::TrackCollection> trackCollection;
02625   //  evt.getByLabel("ctfWithMaterialTracks", trackCollection);
02626   //  evt.getByLabel("generalTracks", trackCollection);
02627   //  const reco::TrackCollection tC = *(trackCollection.product());
02628 
02629   //  std::cout << "ANA: Reconstructed "<< tC.size() << " tracks" << std::endl ;
02630 
02631   // *************************************
02632   /*****
02633   //Get the Vertex Collection
02634   edm::Handle<std::vector<reco::Vertex> > verticies;  evt.getByLabel("offlinePrimaryVertices", verticies);
02635 
02636   //Fill the variables
02637   int _ntracksw5 = 0;
02638   for (std::vector<reco::Vertex>::const_iterator it = verticies->begin(); it != verticies->end(); ++it) {
02639 
02640     //    ntracks->push_back(int(it->tracksSize())); //all tracks considered for vertexing
02641     //    isvalid->push_back(int(it->isValid()));
02642     //    isfake->push_back(int(it->isFake()));
02643 
02644     if(it->tracksSize() > 0) {
02645       std::vector<TrackBaseRef>::const_iterator trackIt;
02646       for( trackIt = it->tracks_begin(); trackIt != it->tracks_end(); trackIt++) {
02647         if(fabs((**trackIt).charge()) <= 1.)  { 
02648           //tracks that contribute with more than 0.5 weight in vertex reconstruction
02649           if (it->trackWeight(*trackIt) >= 0.5 ) 
02650             _ntracksw5++;
02651         }
02652       }
02653     }
02654   }
02655   *****/
02656   // *************************************
02657   
02658 
02659   h_Trk_NTrk->Fill(tC.size());
02660   if (NPMTHits == 0) TrkMultFlagged0->Fill(tC.size());
02661   if (NPMTHits == 1) TrkMultFlagged1->Fill(tC.size());
02662   if (NPMTHits == 2) TrkMultFlagged2->Fill(tC.size());
02663   if (NPMTHits == 3) TrkMultFlagged3->Fill(tC.size());
02664   if (NPMTHits == 4) TrkMultFlagged4->Fill(tC.size());
02665   if (NPMTHits > 1) TrkMultFlaggedM->Fill(tC.size());
02666   for (reco::TrackCollection::const_iterator track=tC.begin(); track!=tC.end(); track++){
02667     h_Trk_pt->Fill(track->pt());
02668   }
02669 
02670 
02671     /****
02672     std::cout << "Track number "<< i << std::endl ;
02673     std::cout << "\tmomentum: " << track->momentum()<< std::endl;
02674     std::cout << "\tPT: " << track->pt()<< std::endl;
02675     std::cout << "\tvertex: " << track->vertex()<< std::endl;
02676     std::cout << "\timpact parameter: " << track->d0()<< std::endl;
02677     std::cout << "\tcharge: " << track->charge()<< std::endl;
02678     std::cout << "\tnormalizedChi2: " << track->normalizedChi2()<< std::endl;
02679 
02680     cout<<"\tFrom EXTRA : "<<endl;
02681     cout<<"\t\touter PT "<< track->outerPt()<<endl;
02682     std::cout << "\t direction: " << track->seedDirection() << std::endl;
02683     ****/
02684 
02685   // ********************************
02686   // *** Photons
02687   // ********************************
02688   /***
02689   edm::Handle<reco::PhotonCollection> photonCollection;
02690   evt.getByLabel("photons", photonCollection);
02691   const reco::PhotonCollection pC = *(photonCollection.product());
02692 
02693   std::cout << "Reconstructed "<< pC.size() << " photons" << std::endl ;
02694   for (reco::PhotonCollection::const_iterator photon=pC.begin(); photon!=pC.end(); photon++){
02695   }
02696   ***/
02697 
02698   // ********************************
02699   // *** Muons
02700   // ********************************
02701   /***
02702   edm::Handle<reco::MuonCollection> muonCollection;
02703   evt.getByLabel("muons", muonCollection);
02704 
02705   const reco::MuonCollection mC = *(muonCollection.product());
02706 
02707   std::cout << "Reconstructed "<< mC.size() << " muons" << std::endl ;
02708   for (reco::MuonCollection::const_iterator muon=mC.begin(); muon!=mC.end(); muon++){
02709   }
02710   ***/
02711 
02712 
02713 
02714 
02715   // ********************************
02716   // *** Events passing seletion cuts
02717   // ********************************
02718 
02719   // --- Cosmic Cleanup
02720   // --- Vertex
02721   // --- Eta 
02722 
02723   int iJet; 
02724   iJet = 0;
02725   for( CaloJetCollection::const_iterator ijet = caloJets->begin(); ijet != caloJets->end(); ++ ijet ) {
02726     
02727     //    if ( (fabs(ijet->eta()) < 1.3) && 
02728     //   (fabs(ijet->pt())  > 20.) ) {
02729 
02730          //      (ijet->emEnergyFraction() > 0.01) &&
02731          //      (ijet->emEnergyFraction() > 0.99) ) {
02732 
02733     iJet++; 
02734     //    if (iJet == 1) {
02735     //      cout << " CaloJet: Event Type = "   << evtType 
02736     //     << " pt = " << ijet->pt()
02737     //     << endl; 
02738     //    }
02739     h_pt->Fill(ijet->pt());
02740     if (evtType == 0) h_ptTower->Fill(ijet->pt());
02741     if (evtType == 1) h_ptRBX->Fill(ijet->pt());
02742     if (evtType == 2) h_ptHPD->Fill(ijet->pt());
02743     h_et->Fill(ijet->et());
02744     h_eta->Fill(ijet->eta());
02745     h_phi->Fill(ijet->phi());
02746 
02747     jetHOEne->Fill(ijet->hadEnergyInHO());    
02748     jetEMFraction->Fill(ijet->emEnergyFraction());    
02749       
02750     //    }    
02751   }
02752 
02753 
02754 
02755   //*****************************
02756   //*** Get the GenJet collection
02757   //*****************************
02758 
02759       /**************
02760   Handle<GenJetCollection> genJets;
02761   evt.getByLabel( GenJetAlgorithm, genJets );
02762 
02763   //Loop over the two leading GenJets and fill some histograms
02764   jetInd    = 0;
02765   allJetInd = 0;
02766   for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end(); ++ gen ) {
02767     allJetInd++;
02768     if (allJetInd == 1) {
02769       p4tmp[0] = gen->p4();
02770     }
02771     if (allJetInd == 2) {
02772       p4tmp[1] = gen->p4();
02773     }
02774 
02775     if ( (allJetInd == 1) || (allJetInd == 2) ) {
02776       h_ptGenL->Fill( gen->pt() );
02777       h_etaGenL->Fill( gen->eta() );
02778       h_phiGenL->Fill( gen->phi() );
02779     }
02780 
02781     if ( gen->pt() > minJetPt) {
02782       // std::cout << "GEN JET1 #" << jetInd << std::endl << gen->print() << std::endl;
02783       h_ptGen->Fill( gen->pt() );
02784       h_etaGen->Fill( gen->eta() );
02785       h_phiGen->Fill( gen->phi() );
02786       jetInd++;
02787     }
02788   }
02789 
02790   h_nGenJets->Fill( jetInd );
02791       *******/
02792   }
02793 
02794 }
02795 
02796 // ***********************************
02797 // ***********************************
02798 void myJetAna::endJob() {
02799 
02800   for (int i=0; i<4000; i++) {
02801     if ((nBNC[i]/totBNC) > 0.05) {
02802       std::cout << "+++ " << i << " " 
02803                 << (nBNC[i]/totBNC) << " "
02804                 << nBNC[i]          << " " 
02805                 << totBNC           << " " 
02806                 << std::endl;      
02807     }
02808   }
02809 
02810 
02811 }
02812 #include "FWCore/Framework/interface/MakerMacros.h"
02813 DEFINE_FWK_MODULE(myJetAna);