ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
(Generate patch)

Comparing UserCode/MPIAnalyzer/src/MPIntuple.cc (file contents):
Revision 1.6 by naodell, Sun May 23 21:56:11 2010 UTC vs.
Revision 1.14 by naodell, Tue Oct 5 21:27:59 2010 UTC

# Line 1 | Line 1
1 /*
2 Description: <one line class summary>
3
4 Implementation:
5     n-tuple creator for jet studies
6 */
7 //
1   // Original Author:  "Nathaniel Odell"
2   // Secondary Author: Steven Won
3   // With contributions from: Andrey Pozdnyakov
4   //         Created:  Thurs April 22  2010
5   // $Id$
13 //
14 //
15
6  
7   // system include files
8   #include <memory>
# Line 25 | Line 15
15   #include "FWCore/Framework/interface/EventSetup.h"
16   #include "FWCore/Framework/interface/Event.h"
17   #include "FWCore/Framework/interface/MakerMacros.h"
28
18   #include "FWCore/ParameterSet/interface/ParameterSet.h"
30
19   #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
20   #include "Geometry/Records/interface/CaloGeometryRecord.h"
21 + #include "DataFormats/Math/interface/deltaPhi.h"
22  
23   // Jet and vertex functions
24   #include "DataFormats/JetReco/interface/CaloJet.h"
# Line 41 | Line 30
30   #include "DataFormats/JetReco/interface/Jet.h"
31   #include "DataFormats/VertexReco/interface/Vertex.h"
32   #include "DataFormats/VertexReco/interface/VertexFwd.h"
33 + #include "DataFormats/BTauReco/interface/JetTag.h"
34 + //#include "RecoVertex/PrimaryVertexProducer/interface/VertexHigherPtSquared.h"
35 +
36 + //GenParticles
37 + //#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
38  
39   // Need to get correctors
40   #include "JetMETCorrections/Objects/interface/JetCorrector.h"
# Line 86 | Line 80 | class MPIntuple : public edm::EDAnalyzer
80        virtual void endJob() ;
81  
82    bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger);
83 +  double sumPtSquared(const Vertex & v);
84  
85        // ----------member data ---------------------------
86  
# Line 93 | Line 88 | class MPIntuple : public edm::EDAnalyzer
88    edm::InputTag GenJetHandle_;
89    edm::InputTag PrimaryVtxHandle_;
90    edm::InputTag triggerResultsTag_;
96
97
91    //  Counting variables   //
92  
93 <  int eventNumber, runNumber, lumiSection;
93 >  int eventNumber, runNumber, lumiSection, bunchCross;
94  
95    TTree* sTree;
96    TFile* ntupleFile;
# Line 108 | Line 101 | class MPIntuple : public edm::EDAnalyzer
101  
102    bool doGenJets_;
103    bool doPFJets_;
104 <  bool isMC;
105 <
104 >  bool triggerHLT_;
105 >  bool isRealData;
106    string rootfilename;
107  
108    //Triggers
# Line 129 | Line 122 | MPIntuple::MPIntuple(const edm::Paramete
122    PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
123    doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")),
124    doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")),
125 +  triggerHLT_(iConfig.getUntrackedParameter<bool>("triggerHLT")),
126    rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")),
127    hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
128    hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT"))
129   {
130 <  //edm::TriggerNames
137 <  // triggerNames(iConfig);
130 >
131  
132   }
133  
# Line 144 | Line 137 | MPIntuple::~MPIntuple()
137  
138   }
139  
140 + double MPIntuple::sumPtSquared(const Vertex & v)
141 + {
142 +  double sum = 0.;
143 +  double pT;
144 +  for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
145 +    pT = (**it).pt();
146 +    double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
147 +
148 +    sum += pT*pT;
149 +  }
150 +  return sum;
151 + }
152 +
153   //
154   // member functions
155   //
# Line 152 | Line 158 | MPIntuple::~MPIntuple()
158   void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
159   {
160  
161 +
162    eventNumber  = iEvent.id().event();
163    runNumber    = iEvent.id().run();
164    lumiSection  = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
165 <  //  isMC = iEvent.???();
166 <  
167 <  int PFcount  = 0;
168 <  int Gencount = 0;
169 <  int Vtxcount = 0;
165 >  bunchCross   = (unsigned int)iEvent.bunchCrossing();
166 >  isRealData   = iEvent.isRealData();
167 >
168 >  int PFcount   = 0;
169 >  int Gencount  = 0;
170 >  int Vtxcount  = 0;
171  
172    if(doPFJets_){
173  
174 +    edm::Handle<reco::JetTagCollection> bTagHandle1;
175 +    iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
176 +    const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
177 +    reco::JetTagCollection::const_iterator jet_it_1;
178 +
179 +    edm::Handle<reco::JetTagCollection> bTagHandle2;
180 +    iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
181 +    const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
182 +    reco::JetTagCollection::const_iterator jet_it_2;
183 +
184 +
185      const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
186      const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
187  
# Line 174 | Line 193 | void MPIntuple::analyze(const edm::Event
193      for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
194        
195        reco::PFJet myJet = reco::PFJet(*jet_iter);
177
178      float scale2 = correctorL2->correction(myJet);
179      float scale3 = correctorL3->correction(myJet);
180
196        if(myJet.et() > 5){
197 +
198 +        for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
199 +           if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
200 +        }
201 +
202 +        for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
203 +           if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
204 +        }
205                  
206          TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
207        
# Line 195 | Line 218 | void MPIntuple::analyze(const edm::Event
218  
219          jetCon->SetNumChPart(myJet.chargedMultiplicity());
220  
221 +        if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
222 +        if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
223 +
224 +        float scale2 = correctorL2->correction(myJet);
225 +        myJet.scaleEnergy(scale2);
226 +        float scale3 = correctorL3->correction(myJet);
227 +        myJet.scaleEnergy(scale3);
228 +
229 +        //more corrections?
230 +
231          jetCon->SetJetCorr(2, scale2);
232          jetCon->SetJetCorr(3, scale3);
200                
233  
234          ++PFcount;
235        }      
236      }  
237    }
238    
239 <  if(doGenJets_){
239 >  if(!isRealData){
240      
241 +
242      Handle<reco::GenJetCollection> GenJets;
243      iEvent.getByLabel(GenJetHandle_, GenJets);
244      
# Line 221 | Line 254 | void MPIntuple::analyze(const edm::Event
254          ++Gencount;
255          
256        }
257 <    }
257 >    }    
258    }
259 <
260 <
259 >  
260 >  
261    Handle<reco::VertexCollection> primaryVtcs;
262    iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
263 <
263 >  
264    for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
265      
266      reco::Vertex myVtx = reco::Vertex(*vtx_iter);
267 <    
267 >    if(!myVtx.isValid() || myVtx.isFake()) continue;
268      TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
269        
270      vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
271      vtxCon->SetNDof(myVtx.ndof());
272      vtxCon->SetChi2(myVtx.chi2());
273 <            
273 >    vtxCon->SetNTrks(myVtx.tracksSize());
274 >    vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
275      ++Vtxcount;
276      
277    }
278  
279  
280 +  if(triggerHLT_){
281  
282 <  //----------  Filling HLT trigger bits!  ------------
282 >    //----------  Filling HLT trigger bits!  ------------
283  
284 <  edm::Handle<TriggerResults> hltR;
285 <  triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
286 <  iEvent.getByLabel(triggerResultsTag_,hltR);
287 <
288 <  const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
289 <  hlNames=triggerNames.triggerNames();
255 <
256 <  string MPI_TriggerNames[32] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"};
257 <
258 <  bool triggerPassed = false;
259 <  triggerStatus = 0x0;
260 <  
261 <  for (uint iT=0; iT<hlNames.size(); ++iT) {
262 <
263 <    triggerPassed = triggerDecision(hltR, iT);
284 >    edm::Handle<TriggerResults> hltR;
285 >    triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
286 >    iEvent.getByLabel(triggerResultsTag_,hltR);
287 >    
288 >    const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
289 >    hlNames=triggerNames.triggerNames();
290      
291 <    if(triggerPassed){
291 >    string MPI_TriggerNames[] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"};
292 >    
293 >    bool triggerPassed = false;
294 >    triggerStatus = 0x0;
295 >    
296 >    for (uint i=0; i<hlNames.size(); ++i) {
297 >      
298 >      triggerPassed = triggerDecision(hltR, i);
299        
300 <      for (int j = 0; j != 32; ++j){
300 >      if(triggerPassed){
301          
302 <        if (hlNames[iT] == MPI_TriggerNames[j])
303 <          {
304 <            cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
305 <            triggerStatus |= 0x01 << j;
306 <          }
302 >        for (uint j = 0; j != 32; ++j){
303 >          
304 >          if (hlNames[i] == MPI_TriggerNames[j])
305 >            {
306 >              //            cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
307 >              triggerStatus |= 0x01 << j;
308 >              
309 >            }
310 >        }
311        }
312 <    }
313 <  }
314 <  
312 >    }
313 >  }
314 >  //  cout<< "total status: "<<hex<<triggerStatus<<endl;
315 >
316  
317 <  if(PFcount > 0 || Gencount > 0 || Vtxcount > 0);  sTree -> Fill();
317 >
318 >  if((PFcount > 3 || Gencount > 3) && Vtxcount > 0)  sTree -> Fill();
319    
320    jet_ak5PF->Clear();
321    primaryVtx->Clear();
# Line 304 | Line 343 | void  MPIntuple::beginJob()
343    sTree->Branch("runNumber",&runNumber, "runNumber/I");
344    sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
345    sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
346 <  sTree->Branch("isMC",&isMC,"isMC/O");
347 <
346 >  sTree->Branch("isRealData",&isRealData, "isRealData/i");
347 >  sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
348  
349   }
350  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines