source: trunk/modules/MadGraphMatchingTreeWriter.cc@ 4

Last change on this file since 4 was 2, checked in by Pavel Demin, 16 years ago

first commit
File size: 5.6 KB
Line 
1
2#include "modules/MadGraphMatchingTreeWriter.h"
3
4
5#include "ExRootAnalysis/ExRootResult.h"
6#include "ExRootAnalysis/ExRootClasses.h"
7#include "ExRootAnalysis/ExRootTreeBranch.h"
8
9#include "ExRootAnalysis/ExRootCandidate.h"
10
11#include "TClonesArray.h"
12
13#include "TH1.h"
14#include "TH2.h"
15#include "TString.h"
16#include "TCanvas.h"
17#include "TLorentzVector.h"
18
19#include <iostream>
20
21using namespace std;
22
23//------------------------------------------------------------------------------
24
25MadGraphMatchingTreeWriter::MadGraphMatchingTreeWriter()
26{
27}
28
29//------------------------------------------------------------------------------
30
31MadGraphMatchingTreeWriter::~MadGraphMatchingTreeWriter()
32{
33}
34
35//------------------------------------------------------------------------------
36
37void MadGraphMatchingTreeWriter::Init()
38{
39 fJetPTMin = GetDouble("JetPTMin", 20.0);
40 fJetEtaMax = GetDouble("JetEtaMax", 4.5);
41
42 // import array with output from filter/classifier/jetfinder modules
43
44 fInputArrayPartonJets = ImportArray(GetString("InputArrayPartonJets", "partonjetfinder/candidates"));
45 fItInputArrayPartonJets = fInputArrayPartonJets->MakeIterator();
46
47 fInputArrayHadronJets = ImportArray(GetString("InputArrayHadronJets", "hadronjetfinder/candidates"));
48 fItInputArrayHadronJets = fInputArrayHadronJets->MakeIterator();
49
50 fInputArrayMatching = ImportArray(GetString("InputArrayMatching", "partonjetfinder/matching"));
51 fItInputArrayMatching = fInputArrayMatching->MakeIterator();
52
53 fInputArrayPartons = ImportArray(GetString("InputArrayPartons", "initstateselection/candidates"));
54 fItInputArrayPartons = fInputArrayPartons->MakeIterator();
55
56 fBranchPartonJets = NewBranch("PartonJet", ExRootGenJet::Class());
57 fBranchHadronJets = NewBranch("HadronJet", ExRootGenJet::Class());
58 fBranchMatching = NewBranch("Match", ExRootMatching::Class());
59 fBranchPartons = NewBranch("Parton", ExRootGenParticle::Class());
60}
61
62//------------------------------------------------------------------------------
63
64void MadGraphMatchingTreeWriter::Finish()
65{
66 if(fItInputArrayPartonJets) delete fItInputArrayPartonJets;
67 if(fItInputArrayHadronJets) delete fItInputArrayHadronJets;
68 if(fItInputArrayMatching) delete fItInputArrayMatching;
69 if(fItInputArrayPartons) delete fItInputArrayPartons;
70}
71
72//------------------------------------------------------------------------------
73
74void MadGraphMatchingTreeWriter::Process()
75{
76 ExRootCandidate *candidate = 0;
77 ExRootGenParticle *entryParton = 0;
78 ExRootMatching *matching = 0, *entryMatching = 0;
79 ExRootGenJet *entryJet = 0;
80 Double_t pt, signPz, eta, rapidity;
81
82 // loop over all parton jets
83 fItInputArrayPartonJets->Reset();
84 while((candidate = static_cast<ExRootCandidate*>(fItInputArrayPartonJets->Next())))
85 {
86 const TLorentzVector &momentum = candidate->GetP4();
87
88 pt = momentum.Pt();
89 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
90 eta = (pt == 0.0 ? signPz*999.9 : momentum.Eta());
91 rapidity = (pt == 0.0 ? signPz*999.9 : momentum.Rapidity());
92
93 if(pt < fJetPTMin) continue;
94 if(TMath::Abs(eta) > fJetEtaMax) continue;
95
96 entryJet = static_cast<ExRootGenJet*>(fBranchPartonJets->NewEntry());
97
98 entryJet->E = momentum.E();
99 entryJet->Px = momentum.Px();
100 entryJet->Py = momentum.Py();
101 entryJet->Pz = momentum.Pz();
102
103 entryJet->Eta = eta;
104 entryJet->Phi = momentum.Phi();
105 entryJet->PT = pt;
106
107 entryJet->Rapidity = rapidity;
108
109 entryJet->Mass = momentum.M();
110 }
111
112 // loop over all hadron jets
113 fItInputArrayHadronJets->Reset();
114 while((candidate = static_cast<ExRootCandidate*>(fItInputArrayHadronJets->Next())))
115 {
116 const TLorentzVector &momentum = candidate->GetP4();
117
118 pt = momentum.Pt();
119 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
120 eta = (pt == 0.0 ? signPz*999.9 : momentum.Eta());
121 rapidity = (pt == 0.0 ? signPz*999.9 : momentum.Rapidity());
122
123 if(pt < fJetPTMin) continue;
124 if(TMath::Abs(eta) > fJetEtaMax) continue;
125
126 entryJet = static_cast<ExRootGenJet*>(fBranchHadronJets->NewEntry());
127
128 entryJet->E = momentum.E();
129 entryJet->Px = momentum.Px();
130 entryJet->Py = momentum.Py();
131 entryJet->Pz = momentum.Pz();
132
133 entryJet->Eta = eta;
134 entryJet->Phi = momentum.Phi();
135 entryJet->PT = pt;
136
137 entryJet->Rapidity = rapidity;
138
139 entryJet->Mass = momentum.M();
140 }
141
142 // loop over all matching
143 fItInputArrayMatching->Reset();
144 while((matching = static_cast<ExRootMatching*>(fItInputArrayMatching->Next())))
145 {
146 entryMatching = static_cast<ExRootMatching*>(fBranchMatching->NewEntry());
147
148 entryMatching->DMerge = matching->DMerge;
149 entryMatching->YMerge = matching->YMerge;
150 }
151
152 // loop over all partons
153 fItInputArrayPartons->Reset();
154 while((candidate = static_cast<ExRootCandidate*>(fItInputArrayPartons->Next())))
155 {
156 const TLorentzVector &momentum = candidate->GetP4();
157
158 entryParton = static_cast<ExRootGenParticle*>(fBranchPartons->NewEntry());
159
160 pt = momentum.Pt();
161 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
162 eta = (pt == 0.0 ? signPz*999.9 : momentum.Eta());
163 rapidity = (pt == 0.0 ? signPz*999.9 : momentum.Rapidity());
164
165 entryParton->PID = candidate->GetType()->PdgCode();
166
167 entryParton->E = momentum.E();
168 entryParton->Px = momentum.Px();
169 entryParton->Py = momentum.Py();
170 entryParton->Pz = momentum.Pz();
171
172 entryParton->Eta = eta;
173 entryParton->Phi = momentum.Phi();
174 entryParton->PT = pt;
175
176 entryParton->Rapidity = rapidity;
177 }
178}
179
180//------------------------------------------------------------------------------
181
Note: See TracBrowser for help on using the repository browser.