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.13 by naodell, Fri Jul 23 12:45:32 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 155 | Line 161 | void MPIntuple::analyze(const edm::Event
161    eventNumber  = iEvent.id().event();
162    runNumber    = iEvent.id().run();
163    lumiSection  = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
164 <  //  isMC = iEvent.???();
165 <  
166 <  int PFcount  = 0;
167 <  int Gencount = 0;
168 <  int Vtxcount = 0;
164 >  bunchCross   = (unsigned int)iEvent.bunchCrossing();
165 >  isRealData   = iEvent.isRealData();
166 >
167 >  int PFcount   = 0;
168 >  int Gencount  = 0;
169 >  int Vtxcount  = 0;
170  
171    if(doPFJets_){
172  
173 +    edm::Handle<reco::JetTagCollection> bTagHandle1;
174 +    iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
175 +    const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
176 +    reco::JetTagCollection::const_iterator jet_it_1;
177 +
178 +    edm::Handle<reco::JetTagCollection> bTagHandle2;
179 +    iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
180 +    const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
181 +    reco::JetTagCollection::const_iterator jet_it_2;
182 +
183      const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
184      const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
185  
# Line 175 | Line 192 | void MPIntuple::analyze(const edm::Event
192        
193        reco::PFJet myJet = reco::PFJet(*jet_iter);
194  
195 <      float scale2 = correctorL2->correction(myJet);
196 <      float scale3 = correctorL3->correction(myJet);
195 > //      float scale2 = correctorL2->correction(myJet);
196 > //      float scale3 = correctorL3->correction(myJet);
197  
198        if(myJet.et() > 5){
199 +
200 +        for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
201 +           if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
202 +        }
203 +
204 +        for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
205 +           if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
206 +        }
207                  
208          TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
209        
# Line 195 | Line 220 | void MPIntuple::analyze(const edm::Event
220  
221          jetCon->SetNumChPart(myJet.chargedMultiplicity());
222  
223 +        if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
224 +        if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
225 +
226 +        float scale2 = correctorL2->correction(myJet);
227 +        myJet.scaleEnergy(scale2);
228 +        float scale3 = correctorL3->correction(myJet);
229 +        myJet.scaleEnergy(scale3);
230 +
231 +        //more corrections?
232 +
233          jetCon->SetJetCorr(2, scale2);
234          jetCon->SetJetCorr(3, scale3);
200                
235  
236          ++PFcount;
237        }      
238      }  
239    }
240    
241 <  if(doGenJets_){
241 >  if(!isRealData){
242      
243 +
244      Handle<reco::GenJetCollection> GenJets;
245      iEvent.getByLabel(GenJetHandle_, GenJets);
246      
# Line 221 | Line 256 | void MPIntuple::analyze(const edm::Event
256          ++Gencount;
257          
258        }
259 <    }
259 >    }    
260    }
261 <
262 <
261 >  
262 >  
263    Handle<reco::VertexCollection> primaryVtcs;
264    iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
265 <
265 >  
266    for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
267      
268      reco::Vertex myVtx = reco::Vertex(*vtx_iter);
269 <    
269 >    if(!myVtx.isValid() || myVtx.isFake()) continue;
270      TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
271        
272      vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
273      vtxCon->SetNDof(myVtx.ndof());
274      vtxCon->SetChi2(myVtx.chi2());
275 <            
275 >    vtxCon->SetNTrks(myVtx.tracksSize());
276 >    vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
277      ++Vtxcount;
278      
279    }
280  
281  
282 +  if(triggerHLT_){
283  
284 <  //----------  Filling HLT trigger bits!  ------------
248 <
249 <  edm::Handle<TriggerResults> hltR;
250 <  triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
251 <  iEvent.getByLabel(triggerResultsTag_,hltR);
284 >    //----------  Filling HLT trigger bits!  ------------
285  
286 <  const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
287 <  hlNames=triggerNames.triggerNames();
288 <
289 <  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"};
290 <
291 <  bool triggerPassed = false;
292 <  triggerStatus = 0x0;
293 <  
261 <  for (uint iT=0; iT<hlNames.size(); ++iT) {
262 <
263 <    triggerPassed = triggerDecision(hltR, iT);
286 >    edm::Handle<TriggerResults> hltR;
287 >    triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
288 >    iEvent.getByLabel(triggerResultsTag_,hltR);
289 >    
290 >    const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
291 >    hlNames=triggerNames.triggerNames();
292 >    
293 >    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"};
294      
295 <    if(triggerPassed){
295 >    bool triggerPassed = false;
296 >    triggerStatus = 0x0;
297 >    
298 >    for (uint i=0; i<hlNames.size(); ++i) {
299 >      
300 >      triggerPassed = triggerDecision(hltR, i);
301        
302 <      for (int j = 0; j != 32; ++j){
302 >      if(triggerPassed){
303          
304 <        if (hlNames[iT] == MPI_TriggerNames[j])
305 <          {
306 <            cout<<"trigger name: "<<hlNames[iT]<<"  status: "<<triggerPassed<<endl;
307 <            triggerStatus |= 0x01 << j;
308 <          }
304 >        for (uint j = 0; j != 32; ++j){
305 >          
306 >          if (hlNames[i] == MPI_TriggerNames[j])
307 >            {
308 >              //            cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
309 >              triggerStatus |= 0x01 << j;
310 >              
311 >            }
312 >        }
313        }
314 <    }
315 <  }
314 >    }
315 >  }
316 >  //  cout<< "total status: "<<hex<<triggerStatus<<endl;
317    
318 <
279 <  if(PFcount > 0 || Gencount > 0 || Vtxcount > 0);  sTree -> Fill();
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