ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Anghel/macros/src/readTree_true_forMM.cc
Revision: 1.1
Committed: Fri Jan 29 23:08:33 2010 UTC (15 years, 3 months ago) by anghel
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
Revision 29Jan2010

File Contents

# User Rev Content
1 anghel 1.1 #include "myGlobalMuon.h"
2     #include "myGlobalMuon_LinkDef.h"
3     #include "myJet.h"
4     #include "myJet_LinkDef.h"
5     #include "myParticle.h"
6     #include "myParticle_LinkDef.h"
7     #include "myGenParticle.h"
8     #include "myGenParticle_LinkDef.h"
9     #include "myEvent.h"
10     #include "myEvent_LinkDef.h"
11    
12     #include "src/readTreeFunc.cc"
13     #include "interface/Utils.h"
14    
15     #include "TChain.h"
16     #include "TFile.h"
17     #include "TH1.h"
18     #include "TH1F.h"
19     #include "TH2.h"
20     #include "TH2D.h"
21     #include "TMath.h"
22     #include <iostream>
23     #include <vector>
24     #include <string>
25     #include <fstream>
26    
27     using namespace std;
28    
29     int main(int argc, char* argv[]) {
30    
31     string treeName;
32     string simType;
33     string flavorFilter;
34     string studyType;
35    
36     //arguments used when running
37     for(int i = 1; i != argc; ++i) {
38     string arg = argv[i];
39     if (arg == "-s") {treeName = string(argv[i+1]); }
40     else if (arg == "-m") {simType = string(argv[i+1]);}
41     else if (arg == "-f") {flavorFilter = string(argv[i+1]); }
42     else if (arg == "-p") {studyType = string(argv[i+1]); }
43     else if (arg == "-h") { WhichOption(); return 0; }
44     //else { std::cout << "\nUnknown parameter " << arg << std::endl; WhichOption(); return 0; }
45     }
46    
47     cout << treeName << " ,\t" << simType << " ,\t" << flavorFilter << " ,\t" << studyType << endl;
48    
49     //initialize variables to set up the names for the histograms for each jet multiplicity bin
50     string njets[5] = {"all","1jet","2jet","3jet","4jet"};
51     int njetbin = 5;
52    
53     //define the histograms that need to be read directly from the input file
54     TH1F *histo_weight;
55    
56     //define histograms put in the output file
57     TH1F* histo_metEt_loose[5];
58     TH1F* histo_dPhiMuMet_loose[5];
59     TH1F* histo_dPhiMuJet2_loose[5];
60     TH1F* histo_muPt_loose[5];
61     TH1F* histo_mud0d0sig_loose[5];
62     TH2D* histo_relIsoVSd0_loose[5];
63     TH2D* histo_relIsoVSmetEt_loose[5];
64     TH2D* histo_relIsoVSmuPt_loose[5];
65     TH2D* histo_relIsoVSmuEta_loose[5];
66     TH2D* histo_relIsoVSdPhimuMet_loose[5];
67    
68     for (int nj=0; nj != njetbin; nj++) {
69     histo_metEt_loose[nj] = new TH1F(("h_metEt_loose_"+njets[nj]).c_str(), "", 500, 0., 500.);
70     histo_dPhiMuMet_loose[nj] = new TH1F(("h_dPhi_MuandMet_loose_"+njets[nj]).c_str(), "", 40, 0., 4.); //dPhi(mu,MET)
71     histo_dPhiMuJet2_loose[nj] = new TH1F(("h_dPhi_Muonand2Jets_loose_"+njets[nj]).c_str(), "", 100, 0., 10.); //add dPhi(mu,jet1)+dPhi(mu,jet2)
72     histo_muPt_loose[nj] = new TH1F(("h_muPt_loose_"+njets[nj]).c_str(), "", 500, 0., 500.);
73     histo_mud0d0sig_loose[nj] = new TH1F(("h_d0d0sig_loose_"+njets[nj]).c_str(), "", 150, 0., 15.);
74    
75     histo_relIsoVSd0_loose[nj] = new TH2D(("h_muRelIso_vs_d0corr_loose_"+njets[nj]).c_str(), "", 10, 0., 1., 150, 0., 15.);
76     histo_relIsoVSmetEt_loose[nj] = new TH2D(("h_muRelIso_vs_metEt_loose_"+njets[nj]).c_str(), "", 10, 0., 1., 500, 0., 500.);
77     histo_relIsoVSmuPt_loose[nj] = new TH2D(("h_muRelIso_vs_muPt_loose_"+njets[nj]).c_str(), "", 10, 0., 1., 500, 0., 500.);
78     histo_relIsoVSmuEta_loose[nj] = new TH2D(("h_muRelIso_vs_muEta_loose_"+njets[nj]).c_str(), "", 10, 0., 1., 21, 0.0, 2.1);
79     histo_relIsoVSdPhimuMet_loose[nj] = new TH2D(("h_muRelIso_vs_dPhimuMet_loose_"+njets[nj]).c_str(), "", 10, 0., 1., 40, 0., 4.);
80     }
81    
82     TH1F histo_muIso_BeforeDR = TH1F("h_muIso_beforeDR", "", 10, 0., 1.);
83     TH1F histo_muIso_AfterDR = TH1F("h_muIso_afterDR", "", 10, 0., 1.);
84     TH1F histo_weightValue = TH1F("h_weight", "", 2, 0., 2.);
85    
86     //define the files location
87     string path = "/uscms_data/d2/omanghel/Analyses/PATtuples/";
88    
89     //input file
90     //TFile *infile;
91    
92     //load all the input files
93     string treeRealName;
94     if (simType == "FastSim") treeRealName = treeName+simType;
95     else treeRealName = treeName;
96     TChain *chain = new TChain(treeRealName.c_str(), "");
97     chain->Add((path+treeName+simType+"/"+treeName+simType+"_"+flavorFilter+"_"+studyType+"/"+treeName+simType+"_"+flavorFilter+"_*_treeP.root").c_str());
98     //chain->Add("/uscms_data/d2/omanghel/Analyses/PATtuples/TTjetsFastSim/TTjetsFastSimnone/TTjetsFastSim_100_treeP.root");
99     cout << (path+treeName+simType+"/"+treeName+simType+"_"+flavorFilter+"_"+studyType+"/*.root").c_str() << endl;
100    
101     //get the number of events
102     size_t nrEvents = chain->GetEntries();
103     cout <<"nrEvents = " << nrEvents <<endl;
104     //define the storing objects
105     vector<myJet> *jetPs_ = new vector<myJet>();
106     vector<myGlobalMuon> *muonPs_ = new vector<myGlobalMuon>();
107     vector<myParticle> *metPs_ = new vector<myParticle>();
108     vector<myGenParticle> *genPs_ = new vector<myGenParticle>();
109     vector<myEvent> *evtPs_ = new vector<myEvent>();
110    
111     //set the branches
112     chain->SetBranchStatus("*", 0);
113     chain->SetBranchStatus("jetParams*", 1);
114     chain->SetBranchAddress("jetParams", &jetPs_);
115    
116     chain->SetBranchStatus("muonParams*", 1);
117     chain->SetBranchAddress("muonParams", &muonPs_);
118    
119     chain->SetBranchStatus("metParams*", 1);
120     chain->SetBranchAddress("metParams", &metPs_);
121    
122     chain->SetBranchStatus("genParticleParams*", 1);
123     chain->SetBranchAddress("genParticleParams", &genPs_);
124    
125     chain->SetBranchStatus("eventParams*", 1);
126     chain->SetBranchAddress("eventParams", &evtPs_);
127    
128     typedef vector<myJet>::const_iterator myJetiter;
129     typedef vector<myGlobalMuon>::const_iterator myGlobalMuoniter;
130     typedef vector<myParticle>::const_iterator myParticleiter;
131     typedef vector<myGenParticle>::const_iterator myGenPiter;
132     typedef vector<myEvent>::const_iterator myEventiter;
133    
134     //variables
135     double weight = 0.0;
136     double dPhimujet = 0.0;
137     double dPhimujetSum = 0.0;
138     double dPhimumet = 0.0;
139     int isL = 0;
140     int isT = 0;
141     int NevLoose[5];
142     int NevTight[5];
143     string currentfileName;
144    
145     for(int i=0; i!=5; i++) {
146     NevLoose[i] = 0;
147     NevTight[i] = 0;
148     }
149    
150     //loop over the number of events
151     //nrEvents = 5000;
152     for(size_t n = 0; n != nrEvents; ++n) {
153     //cout <<"----------New Event---------------" << endl;
154    
155     isL = 0;
156     isT = 0;
157     TH1F* tmp_muIsoB;
158     TH1F* tmp_muIsoA;
159    
160     chain->GetEntry(n);
161    
162     //get input root file
163     //infile = TFile::Open(chain->GetName());
164     TFile *infile = chain->GetCurrentFile();
165     //cout << chain->GetName() << endl;
166     string tname = infile->GetName();
167     //chain->ls();
168     //cout <<infile->IsZombie() << endl;
169     //cout <<infile->IsOpen() << endl;
170     //cout << tname << endl;
171    
172     if (n==0) {
173     //get histograms already filled in the input file
174     //infile->cd();
175     histo_weight = (TH1F*) infile->Get("h_weightFactor");
176     weight = histo_weight->GetBinContent(1);
177     //cout << "weight factor = " << weight << endl;
178     currentfileName = tname;
179     }
180     else if (currentfileName != tname) {
181     currentfileName = tname;
182     infile->cd();
183     infile->cd("Muons");
184     gDirectory->GetObject("h_muIsoBeforeDR",tmp_muIsoB);
185     gDirectory->GetObject("h_muIsoAfterDR",tmp_muIsoA);
186     //TH1F* tmp_muIsoB = (TH1F*) infile->Get("h_muIsoBeforeDR");
187     //TH1F* tmp_muIsoA = (TH1F*) infile->Get("h_muIsoAfterDR");
188     histo_muIso_BeforeDR.Add(tmp_muIsoB);
189     histo_muIso_AfterDR.Add(tmp_muIsoA);
190     }
191     //delete infile;
192     //delete tmp_muIsoB;
193    
194     /***************************************************************************/
195     /*****************************read collections******************************/
196     /***************************************************************************/
197     /**************************************************************************/
198    
199     ////////////////////////////////////////////////////////////////////////////////
200     //////////////fill separately histograms as a function of the # of jets/////////
201     ////////////////////////////////////////////////////////////////////////////////
202     int njEv = jetPs_->end() - jetPs_->begin();
203     //cout << "njEv = " << njEv << endl;
204     if (njEv == 0) continue;
205     if (njEv >= 4) njEv = 4;
206    
207     for (myEventiter ev = evtPs_->begin(); ev != evtPs_->end(); ++ev) {
208     isL = ev->isLoose();
209     isT = ev->isTight();
210     }
211    
212     dPhimujet = 0.0;
213     dPhimujetSum = 0.0;
214     dPhimumet = 0.0;
215     bool isMetfilled = false;
216    
217     for (myGlobalMuoniter mu = muonPs_->begin(); mu != muonPs_->end(); ++mu) {
218     if (isL == 1) {
219     //fill histos for any # of jets
220     histo_muPt_loose[0]->Fill(mu->p4().Pt());
221     histo_mud0d0sig_loose[0]->Fill(fabs(mu->d0corr()));
222     histo_relIsoVSd0_loose[0]->Fill(mu->relIso(),fabs(mu->d0corr()));
223     histo_relIsoVSmuPt_loose[0]->Fill(mu->relIso(),mu->p4().Pt());
224     histo_relIsoVSmuEta_loose[0]->Fill(mu->relIso(),fabs(mu->p4().Eta()));
225     //fill histos separately for each jet bin
226     histo_muPt_loose[njEv]->Fill(mu->p4().Pt());
227     histo_mud0d0sig_loose[njEv]->Fill(fabs(mu->d0corr()));
228     histo_relIsoVSd0_loose[njEv]->Fill(mu->relIso(),fabs(mu->d0corr()));
229     histo_relIsoVSmuPt_loose[njEv]->Fill(mu->relIso(),mu->p4().Pt());
230     histo_relIsoVSmuEta_loose[njEv]->Fill(mu->relIso(),fabs(mu->p4().Eta()));
231     } // end if loose event
232    
233     for (myJetiter jp = jetPs_->begin(); jp != jetPs_->end(); ++jp) {
234     dPhimujet = fabs(getDeltaPhi(mu->p4(),jp->p4()));
235     dPhimujetSum += dPhimujet;
236     if ((jp - (jetPs_->begin())) == 1 ) {
237     if (isL == 1) {
238     histo_dPhiMuJet2_loose[0]->Fill(dPhimujetSum);
239     histo_dPhiMuJet2_loose[njEv]->Fill(dPhimujetSum);
240     }
241     break;
242     }
243     }//end loop over jets
244    
245     for (myParticleiter metp = metPs_->begin(); metp != metPs_->end(); ++metp) {
246     if (!isMetfilled) {
247     if (isL == 1) {
248     histo_metEt_loose[0]->Fill(metp->p4().Et());
249     histo_metEt_loose[njEv]->Fill(metp->p4().Et());
250     }
251     isMetfilled = true;
252     }
253    
254     dPhimumet = fabs(getDeltaPhi(mu->p4(),metp->p4()));
255    
256     if (isL == 1) {
257     histo_dPhiMuMet_loose[0]->Fill(dPhimumet);
258     histo_relIsoVSmetEt_loose[0]->Fill(mu->relIso(),metp->p4().Et());
259     histo_relIsoVSdPhimuMet_loose[0]->Fill(mu->relIso(),dPhimumet);
260     histo_dPhiMuMet_loose[njEv]->Fill(dPhimumet);
261     histo_relIsoVSmetEt_loose[njEv]->Fill(mu->relIso(),metp->p4().Et());
262     histo_relIsoVSdPhimuMet_loose[njEv]->Fill(mu->relIso(),dPhimumet);
263     }
264    
265     }//end loop over MET
266     }// end loop over muons
267    
268     if (isL == 1) {NevLoose[0]++; NevLoose[njEv]++;}
269     if (isT == 1) {NevTight[0]++; NevTight[njEv]++;}
270    
271     }//end loop # of events
272    
273     cout <<"weight - before filling ---------- = " << weight << endl;
274     histo_weightValue.SetBinContent(1,weight);
275    
276     for (int ij=0; ij!=njetbin; ij++){
277     cout << "The number of events in the loose sample (" << njets[ij] << ") = " << NevLoose[ij]<< endl;
278     cout << "The number of events in the tight sample (" << njets[ij] << ") = " << NevTight[ij] << endl;
279     }
280    
281     ofstream otxtfile;
282     string namefile = "inputMM_"+treeName+simType+"_"+flavorFilter+"_true";
283     //otxtfile.open((namefile+".txt").c_str(),ios_base::app);
284     otxtfile.open((namefile+".txt").c_str());
285     for (int i=0; i!=njetbin; i++) {
286     otxtfile << treeName+flavorFilter;
287     otxtfile << "\t" << weight;
288     otxtfile << "\t" << njets[i];
289     otxtfile << "\t" << NevLoose[i];
290     otxtfile << "\t" << NevTight[i];
291     otxtfile << "\n";
292     }
293     otxtfile << endl;
294    
295     //write the histograms into the root file
296     TFile *ofile;
297     ofile = new TFile((namefile+".root").c_str(), "UPDATE");
298     ofile->cd();
299     for (int i=0; i!=njetbin; i++) {
300     string muPtNameloose = histo_muPt_loose[i]->GetName();
301     histo_muPt_loose[i]->Write((muPtNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
302     string metEtNameloose = histo_metEt_loose[i]->GetName();
303     histo_metEt_loose[i]->Write((metEtNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
304     string dPhiMuJet2Nameloose = histo_dPhiMuJet2_loose[i]->GetName();
305     histo_dPhiMuJet2_loose[i]->Write((dPhiMuJet2Nameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
306     string dPhiMuMetNameloose = histo_dPhiMuMet_loose[i]->GetName();
307     histo_dPhiMuMet_loose[i]->Write((dPhiMuMetNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
308     string mud0d0sigNameloose = histo_mud0d0sig_loose[i]->GetName();
309     histo_mud0d0sig_loose[i]->Write((mud0d0sigNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
310     string relIsoVSd0Nameloose = histo_relIsoVSd0_loose[i]->GetName();
311     histo_relIsoVSd0_loose[i]->Write((relIsoVSd0Nameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
312     string relIsoVSmetEtNameloose = histo_relIsoVSmetEt_loose[i]->GetName();
313     histo_relIsoVSmetEt_loose[i]->Write((relIsoVSmetEtNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
314     string relIsoVSmuPtNameloose = histo_relIsoVSmuPt_loose[i]->GetName();
315     histo_relIsoVSmuPt_loose[i]->Write((relIsoVSmuPtNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
316     string relIsoVSmuEtaNameloose = histo_relIsoVSmuEta_loose[i]->GetName();
317     histo_relIsoVSmuEta_loose[i]->Write((relIsoVSmuEtaNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
318     string relIsoVSdPhimuMetNameloose = histo_relIsoVSdPhimuMet_loose[i]->GetName();
319     histo_relIsoVSdPhimuMet_loose[i]->Write((relIsoVSdPhimuMetNameloose+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
320    
321     if( i==0 ){
322     string muIso_BeforeDRName = histo_muIso_BeforeDR.GetName();
323     histo_muIso_BeforeDR.Write((muIso_BeforeDRName+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
324     string muIso_AfterDRName = histo_muIso_AfterDR.GetName();
325     histo_muIso_AfterDR.Write((muIso_AfterDRName+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
326    
327     string weightName = histo_weightValue.GetName();
328     histo_weightValue.Write((weightName+"_"+treeName+"_"+flavorFilter).c_str(), TObject::kOverwrite);
329     }
330     }
331    
332     ofile->Close();
333     return 0;
334     }