ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Updater.cc
Revision: 1.1
Committed: Mon Oct 31 11:05:23 2011 UTC (13 years, 6 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Log Message:
first working version of the updater to recompute the trigger weights without redoing the whole step2

File Contents

# User Rev Content
1 arizzi 1.1 #include <TH1F.h>
2     #include <TH3F.h>
3     #include <TH2F.h>
4     #include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
5     #include <TROOT.h>
6     #include <TFile.h>
7     #include <TTree.h>
8     #include <TSystem.h>
9     #include "FWCore/ParameterSet/interface/ProcessDesc.h"
10     #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
11     #include "DataFormats/Math/interface/deltaR.h"
12     #include "DataFormats/Math/interface/deltaPhi.h"
13    
14     //btagging
15     //#include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
16     //#include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
17     #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
18    
19     #include <sstream>
20     #include <string>
21    
22     #define MAXJ 30
23     #define MAXL 10
24     #define MAXB 10
25     typedef struct
26     {
27     float et;
28     float sumet;
29     float sig;
30     float phi;
31     } METInfo;
32    
33    
34     int main(int argc, char* argv[])
35     {
36     gROOT->Reset();
37    
38     // parse arguments
39     if ( argc < 2 ) {
40     return 0;
41     }
42     // get the python configuration
43     PythonProcessDesc builder(argv[1]);
44     const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
45     const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
46     const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
47     // now get each parameter
48     // int maxEvents_( in.getParameter<int>("maxEvents") );
49     // int skipEvents_( in.getParameter<int>("skipEvents") );
50     // unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
51     std::string outputFile_( out.getParameter<std::string>("fileName" ) );
52     std::string inputFile_( in.getParameter<std::string>("fileName" ) );
53     bool replaceWeights( ana.getParameter<bool>("replaceWeights" ) );
54     TriggerWeight triggerWeight(ana);
55    
56    
57    
58     //Get old file, old tree and set top branch address
59     TFile *oldfile = new TFile(inputFile_.c_str());
60     TTree *oldtree = (TTree*)oldfile->Get("tree");
61     TH1F * count = (TH1F*)oldfile->Get("Count");
62     TH1F * countWithPU = (TH1F*)oldfile->Get("CountWithPU");
63     TH1F * countWithPU2011B = (TH1F*)oldfile->Get("CountWithPU2011B");
64     TH3F * input3DPU = (TH3F*)oldfile->Get("Input3DPU");
65    
66    
67     Int_t nentries = (Int_t)oldtree->GetEntries();
68     Int_t nvlep;
69     Float_t vLepton_pt[50];
70     Float_t vLepton_eta[50];
71     Int_t nhJets;
72     Float_t hJet_pt[50];
73     Float_t hJet_eta[50];
74     Int_t naJets;
75     Float_t aJet_pt[50];
76     Float_t aJet_eta[50];
77     Int_t Vtype;
78     METInfo MET;
79    
80     oldtree->SetBranchAddress("nvlep", &nvlep);
81     oldtree->SetBranchAddress("nhJets", &nhJets);
82     oldtree->SetBranchAddress("naJets", &naJets);
83    
84     oldtree->SetBranchAddress("vLepton_pt", vLepton_pt);
85     oldtree->SetBranchAddress("vLepton_eta", vLepton_eta);
86    
87     oldtree->SetBranchAddress("hJet_pt", hJet_pt);
88     oldtree->SetBranchAddress("hJet_eta", hJet_eta);
89    
90     oldtree->SetBranchAddress("aJet_pt", aJet_pt);
91     oldtree->SetBranchAddress("aJet_eta", aJet_eta);
92    
93     oldtree->SetBranchAddress("Vtype", &Vtype);
94     oldtree->SetBranchAddress("MET", &MET);
95    
96    
97     Float_t weightTrig;
98     Float_t weightTrigMay;
99     Float_t weightTrigV4;
100     Float_t weightTrigMET;
101     Float_t weightTrigOrMu30;
102     Float_t weightEleRecoAndId;
103     Float_t weightEleTrigJetMETPart;
104     Float_t weightEleTrigElePart;
105     if(replaceWeights)
106     {
107     std::cout << "Replacing the weights in the same branch names" << std::endl;
108     oldtree->SetBranchAddress("weightTrig", &weightTrig);
109     oldtree->SetBranchAddress("weightTrigMay", &weightTrigMay);
110     oldtree->SetBranchAddress("weightTrigV4", &weightTrigV4);
111     oldtree->SetBranchAddress("weightTrigMET", &weightTrigMET);
112     oldtree->SetBranchAddress("weightTrigOrMu30", &weightTrigOrMu30);
113     oldtree->SetBranchAddress("weightEleRecoAndId", &weightEleRecoAndId);
114     oldtree->SetBranchAddress("weightEleTrigJetMETPart", &weightEleTrigJetMETPart);
115     oldtree->SetBranchAddress("weightEleTrigElePart", &weightEleTrigElePart);
116     }
117     //Create a new file + a clone of old tree in new file + clone of the histos + additional/updated weights
118    
119     TFile *newfile = new TFile(outputFile_.c_str(),"RECREATE");
120     TTree *newtree = oldtree->CloneTree(0);
121     if(count) count->Clone()->Write();
122     if(countWithPU) countWithPU->Clone()->Write();
123     if(countWithPU2011B) countWithPU2011B->Clone()->Write();
124     if(input3DPU) input3DPU->Clone()->Write();
125    
126    
127     if(!replaceWeights)
128     {
129     std::cout << "Creating new branch names with _up postfix" << std::endl;
130     newtree->Branch("weightTrig_up", &weightTrig, "weightTrig_up/F");
131     newtree->Branch("weightTrigMay_up", &weightTrigMay,"weightTrigMay/F");
132     newtree->Branch("weightTrigV4_up", &weightTrigV4,"weightTrigV4/F");
133     newtree->Branch("weightTrigMET_up", &weightTrigMET,"weightTrigMET/F");
134     newtree->Branch("weightTrigOrMu30_up", &weightTrigOrMu30,"weightTrigOrMu30/F");
135     newtree->Branch("weightEleRecoAndId_up", &weightEleRecoAndId,"weightEleRecoAndId/F");
136     newtree->Branch("weightEleTrigJetMETPart_up", &weightEleTrigJetMETPart,"weightEleTrigJetMETPart/F");
137     newtree->Branch("weightEleTrigElePart_up", &weightEleTrigElePart,"weightEleTrigElePart/F");
138    
139     }
140    
141    
142    
143     // Float_t weightTrigUpdate;
144     // newtree->Branch("weightTrigUpdate" , &weightTrigUpdate , "weightTrigUpdate/F");
145    
146     for (Int_t i=0;i<nentries; i++) {
147     oldtree->GetEntry(i);
148    
149     std::vector<float> jet30eta;
150     std::vector<float> jet30pt;
151     for( int j = 0 ; j < nhJets; j++) if(hJet_pt[j]>30 ) { jet30eta.push_back(hJet_eta[j]); jet30pt.push_back(hJet_pt[j]); }
152     for( int j = 0 ; j < naJets; j++) if(aJet_pt[j]>30 ) { jet30eta.push_back(aJet_eta[j]); jet30pt.push_back(aJet_pt[j]); }
153    
154     if(Vtype == 0 ){
155     float cweightID = triggerWeight.scaleMuID(vLepton_pt[0],vLepton_eta[0]) * triggerWeight.scaleMuID(vLepton_pt[1],vLepton_eta[1]) ;
156     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLepton_pt[0],vLepton_eta[0]);
157     float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLepton_pt[1],vLepton_eta[1]);
158     float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
159     weightTrig = cweightID * cweightTrig;
160    
161     }
162     if( Vtype == 1 ){
163     std::vector<float> pt,eta;
164     pt.push_back(vLepton_pt[0]); eta.push_back(vLepton_eta[0]);
165     pt.push_back(vLepton_pt[1]); eta.push_back(vLepton_eta[1]);
166     weightEleRecoAndId=triggerWeight.scaleID95Ele(vLepton_pt[0],vLepton_eta[0]) * triggerWeight.scaleRecoEle(vLepton_pt[0],vLepton_eta[0]) *
167     triggerWeight.scaleID95Ele(vLepton_pt[1],vLepton_eta[1]) * triggerWeight.scaleRecoEle(vLepton_pt[1],vLepton_eta[1]);
168     weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
169     //REMOVE FOR "float" for newer ntupler and add branch
170     float weightEleTrigEleAugPart = triggerWeight.scaleDoubleEle17Ele8Aug(pt,eta);
171     weightTrig = (weightEleTrigElePart*1.14+weightEleTrigEleAugPart*0.98 )/2.12 * weightEleRecoAndId;
172    
173    
174    
175     }
176     if(Vtype == 2 ){
177     float cweightID = triggerWeight.scaleMuID(vLepton_pt[0],vLepton_eta[0]);
178     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLepton_pt[0],vLepton_eta[0]);
179     float cweightTrig = weightTrig1;
180     weightTrig = cweightID * cweightTrig;
181     float weightTrig1OrMu30 = triggerWeight.scaleMuOr30IsoHLT(vLepton_pt[0],vLepton_eta[0]);
182     weightTrigOrMu30 = cweightID*weightTrig1OrMu30;
183    
184     }
185     if( Vtype == 3 ){
186     weightTrigMay = triggerWeight.scaleSingleEleMay(vLepton_pt[0],vLepton_eta[0]);
187     weightTrigV4 = triggerWeight.scaleSingleEleV4(vLepton_pt[0],vLepton_eta[0]);
188     weightEleRecoAndId=triggerWeight.scaleID80Ele(vLepton_pt[0],vLepton_eta[0]) * triggerWeight.scaleRecoEle(vLepton_pt[0],vLepton_eta[0]);
189     weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30pt,jet30eta)*triggerWeight.scalePFMHTEle(MET.et);
190     weightEleTrigElePart= weightTrigV4; //this is for debugging only, checking only the V4 part
191    
192     weightTrigMay*=weightEleRecoAndId;
193     weightTrigV4*=weightEleRecoAndId;
194     weightTrigV4*=weightEleTrigJetMETPart;
195     // weightTrig = weightTrigMay * 0.187 + weightTrigV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
196     weightTrig = (weightTrigMay * 0.215 + weightTrigV4 * 1.915)/ 2.13; //FIXME: use proper lumi if we reload 2.fb
197    
198    
199     }
200     if( Vtype == 4 ){
201     float weightTrig1 = triggerWeight.scaleMetHLT(MET.et);
202     weightTrig = weightTrig1;
203     weightTrigMET = weightTrig1;
204    
205     }
206    
207    
208    
209    
210     newtree->Fill();
211     }
212     newtree->Print();
213     newtree->AutoSave();
214     delete oldfile;
215     delete newfile;
216     }