ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/dasu/UltraFastSim/ZHAnalysis.cc

Comparing UserCode/dasu/UltraFastSim/ZHAnalysis.cc (file contents):
Revision 1.1 by dasu, Fri Feb 18 05:06:02 2011 UTC vs.
Revision 1.4 by brownson, Mon Feb 28 20:54:48 2011 UTC

#	Line 1 | Line 1
1		#include "ZHAnalysis.h"
2
3	–	#include <iostream>
4	–	#include <vector>
5	–	#include <string>
3		using namespace std;
4
5	<	#include "Pythia.h"
6	<	using namespace Pythia8;
7	<
8	<	#include "UltraFastSim.h"
12	<
13	<	#include "TFile.h"
14	<	#include "TH1F.h"
15	<
16	<	double CorrectPhi(double phi, double x, double y)
17	<	{
18	<	if(phi>0 && x<0 && y<0)phi+=M_PI;
19	<	if(phi<0 && x>0 && y<0)phi=(2*M_PI-fabs(phi));
20	<	if(phi<0 && x<0 && y>0)phi=(M_PI-fabs(phi));
21	<	//if(y<0)phi-=(2*M_PI);//without this you have 0<phi<2pi, but with this
22	<	return phi;//you get -pi<phi<pi
23	<	}
24	<
25	<	ZHAnalysis::ZHAnalysis(TFile *o, Pythia *p, UltraFastSim *u, bool v) : outFile(o), pythia(p), ufs(u), verbose_(v),
26	<	AnalyzedEvents(0),
27	<	LeptonPairPreSelection(0),LeptonPairSelection(0),
28	<	BJetPairPreSelection(0),BJetPairSelection(0),
29	<	ZptSelection(0),HptSelection(0),
30	<	DphiSelection(0),InvMassSelection(0)
31	<	{
32	<	outFile->cd();
33	<
34	<	hptmuons = new TH1F("Muons Pt","Muons Pt",300,0.,300.);
35	<	hetamuons = new TH1F("Muons Eta","Muons Eta",100,-5.,5.);
36	<	hphimuons = new TH1F("Muons Phi","Muons Phi",100,-7.,7.);
37	<
38	<	hetbjets = new TH1F("B-jets Et","B-jets Et",300,0.,300.);
39	<	hetabjets = new TH1F("B-jets Eta","B-jets Eta",100,-5.,5.);
40	<	hphibjets = new TH1F("B-jets Phi","B-jets Phi",100,-7.,7.);
41	<
42	<	hptz     = new TH1F("Z Pt","Z Pt",300,0.,300.);
43	<	hzinvmass= new TH1F("Z invmass","Z invmass",300,0.,300.);
44	<
45	<	hpth     = new TH1F("H Pt","H Pt",300,0.,300.);
46	<	hhinvmass= new TH1F("H invmass preselection","H invmass preselection",300,0.,300.);
47	<	hhinvmass2= new TH1F("H invmass","H invmass",300,0.,300.);
48	<
49	<	hphiz= new TH1F("Z Phi","Z Phi",100,-7.,7.);
50	<	hphih= new TH1F("H Phi","H Phi",100,-7.,7.);
51	<
52	<	hdphiZH= new TH1F("Dphi(Z,H)","Dphi(Z,H)",100,-7.,7.);
53	<
54	<
5	>	void ZHAnalysis::BookTree() {
6	>	myevent_ = new myevent;
7	>	t=new TTree("t","tree");
8	>	t->Branch("myevent","myevent",&myevent_,256000,1);
9		}
10
11
12	<	bool ZHAnalysis::run() {
12	>	bool ZHAnalysis::Run(UltraFastSim * ufs){
13
14	<	vector<TParticle> MyMuons = ufs->muonList();
15	<	vector<fastjet::PseudoJet> MyBJets = ufs->bJetList();
14	>	vector<TParticle>          GenTaus_ = ufs->genTauList();
15	>	vector<TParticle>          BQuarks_ = ufs->bQuarkList();
16	>	vector<TParticle>          CQuarks_ = ufs->cQuarkList();
17	>	vector<TParticle>          Photons_ = ufs->photonList();
18	>	vector<TParticle>          Electrons_ = ufs->electronList();
19	>	vector<TParticle>          Muons_ = ufs->muonList();
20	>	vector<TParticle>          MuonsStdGeom_ = ufs->muonListStdGeom();
21	>	vector<TParticle>          Taus_ = ufs->tauList();
22	>	vector<TParticle>          ChargedHadrons_ = ufs->chargedHadronList();
23	>	vector<TParticle>          NeutralHadrons_ = ufs->neutralHadronList();
24	>	vector<TLorentzVector>     Jets_ = ufs->jetList();
25	>	vector<TLorentzVector>     BJets_ = ufs->bJetList();
26	>	vector<TLorentzVector>     BJetsStdGeom_ = ufs->bJetListStdGeom();
27	>
28	>	myevent_->clear();
29	>
30	>	for(vector<TParticle>::const_iterator it=GenTaus_.begin();it!=GenTaus_.end();it++){
31	>	myobject obj;
32	>	obj.px=it->Px();
33	>	obj.py=it->Py();
34	>	obj.pz=it->Pz();
35	>	obj.E=it->Energy();
36	>	obj.pt=it->Pt();
37	>	obj.eta=it->Eta();
38	>	obj.phi=it->Phi();
39	>	(myevent_->GenTaus).push_back(obj);
40	>	}
41	>	for(vector<TParticle>::const_iterator it=BQuarks_.begin();it!=BQuarks_.end();it++){
42	>	myobject obj;
43	>	obj.px=it->Px();
44	>	obj.py=it->Py();
45	>	obj.pz=it->Pz();
46	>	obj.E=it->Energy();
47	>	obj.pt=it->Pt();
48	>	obj.eta=it->Eta();
49	>	obj.phi=it->Phi();
50	>	(myevent_->BQuarks).push_back(obj);
51	>	}
52	>	for(vector<TParticle>::const_iterator it=CQuarks_.begin();it!=CQuarks_.end();it++){
53	>	myobject obj;
54	>	obj.px=it->Px();
55	>	obj.py=it->Py();
56	>	obj.pz=it->Pz();
57	>	obj.E=it->Energy();
58	>	obj.pt=it->Pt();
59	>	obj.eta=it->Eta();
60	>	obj.phi=it->Phi();
61	>	(myevent_->CQuarks).push_back(obj);
62	>	}
63	>	for(vector<TParticle>::const_iterator it=Photons_.begin();it!=Photons_.end();it++){
64	>	myobject obj;
65	>	obj.px=it->Px();
66	>	obj.py=it->Py();
67	>	obj.pz=it->Pz();
68	>	obj.E=it->Energy();
69	>	obj.pt=it->Pt();
70	>	obj.eta=it->Eta();
71	>	obj.phi=it->Phi();
72	>	(myevent_->Photons).push_back(obj);
73	>	}
74	>	for(vector<TParticle>::const_iterator it=Electrons_.begin();it!=Electrons_.end();it++){
75	>	myobject obj;
76	>	obj.px=it->Px();
77	>	obj.py=it->Py();
78	>	obj.pz=it->Pz();
79	>	obj.E=it->Energy();
80	>	obj.pt=it->Pt();
81	>	obj.eta=it->Eta();
82	>	obj.phi=it->Phi();
83	>	(myevent_->Electrons).push_back(obj);
84	>	}
85	>	for(vector<TParticle>::const_iterator it=Muons_.begin();it!=Muons_.end();it++){
86	>	myobject obj;
87	>	obj.px=it->Px();
88	>	obj.py=it->Py();
89	>	obj.pz=it->Pz();
90	>	obj.E=it->Energy();
91	>	obj.pt=it->Pt();
92	>	obj.eta=it->Eta();
93	>	obj.phi=it->Phi();
94	>	(myevent_->Muons).push_back(obj);
95	>	}
96	>	for(vector<TParticle>::const_iterator it=MuonsStdGeom_.begin();it!=MuonsStdGeom_.end();it++){
97	>	myobject obj;
98	>	obj.px=it->Px();
99	>	obj.py=it->Py();
100	>	obj.pz=it->Pz();
101	>	obj.E=it->Energy();
102	>	obj.pt=it->Pt();
103	>	obj.eta=it->Eta();
104	>	obj.phi=it->Phi();
105	>	(myevent_->MuonsStdGeom).push_back(obj);
106	>	}
107	>	for(vector<TParticle>::const_iterator it=Taus_.begin();it!=Taus_.end();it++){
108	>	myobject obj;
109	>	obj.px=it->Px();
110	>	obj.py=it->Py();
111	>	obj.pz=it->Pz();
112	>	obj.E=it->Energy();
113	>	obj.pt=it->Pt();
114	>	obj.eta=it->Eta();
115	>	obj.phi=it->Phi();
116	>	(myevent_->Taus).push_back(obj);
117	>	}
118	>	for(vector<TParticle>::const_iterator it=ChargedHadrons_.begin();it!=ChargedHadrons_.end();it++){
119	>	myobject obj;
120	>	obj.px=it->Px();
121	>	obj.py=it->Py();
122	>	obj.pz=it->Pz();
123	>	obj.E=it->Energy();
124	>	obj.pt=it->Pt();
125	>	obj.eta=it->Eta();
126	>	obj.phi=it->Phi();
127	>	(myevent_->ChargedHadrons).push_back(obj);
128	>	}
129	>	for(vector<TParticle>::const_iterator it=NeutralHadrons_.begin();it!=NeutralHadrons_.end();it++){
130	>	myobject obj;
131	>	obj.px=it->Px();
132	>	obj.py=it->Py();
133	>	obj.pz=it->Pz();
134	>	obj.E=it->Energy();
135	>	obj.pt=it->Pt();
136	>	obj.eta=it->Eta();
137	>	obj.phi=it->Phi();
138	>	(myevent_->NeutralHadrons).push_back(obj);
139	>	}
140	>	for(vector<TLorentzVector>::const_iterator it=Jets_.begin();it!=Jets_.end();it++){
141	>	myobject obj;
142	>	obj.px=it->Px();
143	>	obj.py=it->Py();
144	>	obj.pz=it->Pz();
145	>	obj.E=it->E();
146	>	obj.pt=it->Pt();
147	>	obj.eta=it->Eta();
148	>	obj.phi=it->Phi();
149	>	(myevent_->Jets).push_back(obj);
150	>	}
151	>	for(vector<TLorentzVector>::const_iterator it=BJets_.begin();it!=BJets_.end();it++){
152	>	myobject obj;
153	>	obj.px=it->Px();
154	>	obj.py=it->Py();
155	>	obj.pz=it->Pz();
156	>	obj.E=it->E();
157	>	obj.pt=it->Pt();
158	>	obj.eta=it->Eta();
159	>	obj.phi=it->Phi();
160	>	(myevent_->BJets).push_back(obj);
161	>	}
162	>	for(vector<TLorentzVector>::const_iterator it=BJetsStdGeom_.begin();it!=BJetsStdGeom_.end();it++){
163	>	myobject obj;
164	>	obj.px=it->Px();
165	>	obj.py=it->Py();
166	>	obj.pz=it->Pz();
167	>	obj.E=it->E();
168	>	obj.pt=it->Pt();
169	>	obj.eta=it->Eta();
170	>	obj.phi=it->Phi();
171	>	(myevent_->BJetsStdGeom).push_back(obj);
172	>	}
173
174	<	for(vector<TParticle>::const_iterator itm=MyMuons.begin();itm!=MyMuons.end();itm++){
64	<	hptmuons->Fill(itm->Pt());
65	<	hetamuons->Fill(itm->Eta());
66	<	hphimuons->Fill(itm->Phi());
67	<	}
68	<
69	<	for(vector<fastjet::PseudoJet>::const_iterator itb=MyBJets.begin();itb!=MyBJets.end();itb++){
70	<	hetbjets->Fill(itb->perp());
71	<	hetabjets->Fill(itb->eta());
72	<	hphibjets->Fill(itb->phi());
73	<	}
74	<
75	<	//Event selection
76	<
77	<	double Zinvmass=0;
78	<	double Zpt=0;
79	<	double Hinvmass=0;
80	<	double Hpt=0;
81	<	double Hphi=0;
82	<	double Zphi=0;
83	<	double DphiZH=0;
84	<
85	<	AnalyzedEvents++;
86	<
87	<	if(MyMuons.size()==2){
88	<
89	<	LeptonPairPreSelection++;
90	<
91	<	if(MyMuons[0].Pt()>20. && MyMuons[1].Pt()>20.){
92	<
93	<	LeptonPairSelection++;
94	<
95	<	Zinvmass=(pow((MyMuons[0].Energy()+MyMuons[1].Energy()),2)-
96	<	pow((MyMuons[0].Px()+MyMuons[1].Px()),2)-
97	<	pow((MyMuons[0].Py()+MyMuons[1].Py()),2)-
98	<	pow((MyMuons[0].Pz()+MyMuons[1].Pz()),2));
99	<	if(Zinvmass>0)Zinvmass=sqrt(Zinvmass);
100	<	else Zinvmass=0;
101	<	Zpt=sqrt(pow((MyMuons[0].Px()+MyMuons[1].Px()),2)+pow((MyMuons[0].Py()+MyMuons[1].Py()),2));
102	<
103	<	double Zpx=0;
104	<	double Zpy=0;
105	<	Zpx=MyMuons[0].Px()+MyMuons[1].Px();
106	<	Zpy=MyMuons[0].Py()+MyMuons[1].Py();
107	<	Zphi=atan(Zpy/Zpx);
108	<	Zphi=CorrectPhi(Zphi,Zpx,Zpy);
109	<
110	<	hzinvmass->Fill(Zinvmass);
111	<	hptz->Fill(Zpt);
112	<	hphiz->Fill(Zphi);
113	<
114	<	if(MyBJets.size()==2){
115	<
116	<	BJetPairPreSelection++;
117	<
118	<	if(MyBJets[0].perp()>30. && MyBJets[1].perp()>30.){
119	<
120	<	BJetPairSelection++;
121	<
122	<	Hinvmass=(pow((MyBJets[0].E()+MyBJets[1].E()),2)-
123	<	pow((MyBJets[0].px()+MyBJets[1].px()),2)-
124	<	pow((MyBJets[0].py()+MyBJets[1].py()),2)-
125	<	pow((MyBJets[0].pz()+MyBJets[1].pz()),2));
126	<	if(Hinvmass>0)Hinvmass=sqrt(Hinvmass);
127	<	else Hinvmass=0;
128	<
129	<	Hpt=sqrt(pow((MyBJets[0].px()+MyBJets[1].px()),2)+pow((MyBJets[0].py()+MyBJets[1].py()),2));
130	<
131	<	double Hpx=0;
132	<	double Hpy=0;
133	<	Hpx=MyBJets[0].px()+MyBJets[1].px();
134	<	Hpy=MyBJets[0].py()+MyBJets[1].py();
135	<	Hphi=atan(Hpy/Hpx);
136	<	Hphi=CorrectPhi(Hphi,Hpx,Hpy);
137	<
138	<	hhinvmass->Fill(Hinvmass);
139	<	hpth->Fill(Hpt);
140	<	hphih->Fill(Hphi);
141	<
142	<	DphiZH=fabs(Hphi-Zphi);
143	<	if(DphiZH>M_PI)DphiZH = (2.*M_PI-DphiZH);
144	<	hdphiZH->Fill(DphiZH);
145	<
146	<	if(Zpt>150.){
147	<
148	<	ZptSelection++;
149	<
150	<	if(Hpt>150.){
151	<
152	<	HptSelection++;
153	<
154	<	hhinvmass2->Fill(Hinvmass);
155	<	if(Hinvmass<140. && Hinvmass>100){
156	<
157	<	InvMassSelection++;
158	<
159	<	if(DphiZH>2.5){
160	<
161	<	DphiSelection++;
162	<
163	<	}//dphi cut
164	<	}//inv mass cut
165	<	}//h pt cut
166	<	}//z pt cut
167	<	}//bjet et > 30 GeV
168	<	}//bjet pair
169	<	}//muon pt>20 GeV
170	<	}//muon.size==2
174	>	t->Fill();
175
176		return true;
177
178		}//Run
175	–
176	–	bool ZHAnalysis::end() {
177	–	outFile->cd();
178	–	return true;
179	–	}
180	–
181	–	ZHAnalysis::~ZHAnalysis() {
182	–	cout<<" ====================================== " <<endl;
183	–	cout<<" Analyzed Events               : "<<AnalyzedEvents<<endl;
184	–	cout<<" Lepton pair preselection eff.(%) : "<<100.*double(LeptonPairPreSelection)/AnalyzedEvents<<endl;
185	–	cout<<" Lepton pair    selection eff.(%) : "<<100.*double(LeptonPairSelection)/LeptonPairPreSelection<<endl;
186	–	cout<<" B-jet  pair preselection eff.(%) : "<<100.*double(BJetPairPreSelection)/LeptonPairSelection<<endl;
187	–	cout<<" B-jet  pair    selection eff.(%) : "<<100.*double(BJetPairSelection)/BJetPairPreSelection<<endl;
188	–	cout<<" Z pt                     eff.(%) : "<<100.*double(ZptSelection)/BJetPairSelection<<endl;
189	–	cout<<" H pt                     eff.(%) : "<<100.*double(HptSelection)/ZptSelection<<endl;
190	–	cout<<" H inv. mass              eff.(%) : "<<100.*double(InvMassSelection)/HptSelection<<endl;
191	–	cout<<" Dphi(H,Z)                eff.(%) : "<<100.*double(DphiSelection)/InvMassSelection<<endl;
192	–	cout<<" Total selection          eff.(%) : "<<100.*double(DphiSelection)/AnalyzedEvents<<endl;
193	–	cout<<" ====================================== " <<endl;
194	–	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines