1 |
//--------------------------------------------------------------------------------------------------
|
2 |
// $Id: MCParticle.h,v 1.25 2012/03/29 23:41:55 paus Exp $
|
3 |
//
|
4 |
// MCParticle
|
5 |
//
|
6 |
// Stores MC information for both generated and simulated particles.
|
7 |
//
|
8 |
// Authors: C.Loizides, J.Bendavid, C.Paus
|
9 |
//--------------------------------------------------------------------------------------------------
|
10 |
|
11 |
#ifndef MITANA_DATATREE_MCPARTICLE_H
|
12 |
#define MITANA_DATATREE_MCPARTICLE_H
|
13 |
|
14 |
#include <TDatabasePDG.h>
|
15 |
#include <TParticlePDG.h>
|
16 |
#include "MitCommon/DataFormats/interface/Vect3.h"
|
17 |
#include "MitCommon/DataFormats/interface/Vect4M.h"
|
18 |
#include "MitAna/DataCont/interface/Ref.h"
|
19 |
#include "MitAna/DataTree/interface/CompositeParticle.h"
|
20 |
|
21 |
namespace mithep
|
22 |
{
|
23 |
class MCParticle : public CompositeParticle
|
24 |
{
|
25 |
public:
|
26 |
enum EPartType {
|
27 |
kUnknown=0,
|
28 |
kUp=1, kDown=2, kStrange=3, kCharm=4, kBottom=5, kTop=6,
|
29 |
kEl=11, kMu=13, kTau=15,
|
30 |
kElNu=12, kMuNu=14, kTauNu=16,
|
31 |
kGlu=21, kGamma=22, kZ=23, kW=24, kH=25,
|
32 |
kZp=32, kZpp=33, kWp=34, kH0=35, kA0=36, kHp=37,
|
33 |
kPi0=111, kEta=221, kKLong=130, kKShort=310, kK0=311,
|
34 |
kD0=413, kB0=511, kB0Bar=513, kJPsi=443, kUpsilon=553,
|
35 |
kProton=2122, kNeutron=2122, kDeltaPlusPlus = 2224,
|
36 |
kDeltaPlus = 2214, kDelta0 = 2114, kDeltaMinus=1114
|
37 |
};
|
38 |
|
39 |
MCParticle() : fPdgId(0), fStatus(0), fIsGenerated(kFALSE), fIsSimulated(kFALSE) {}
|
40 |
MCParticle(Double_t px, Double_t py, Double_t pz, Double_t e, Int_t id, Int_t s) :
|
41 |
fPdgId(id), fStatus(s), fMom(FourVector(px,py,pz,e)),
|
42 |
fIsGenerated(kFALSE), fIsSimulated(kFALSE) {}
|
43 |
|
44 |
Int_t AbsPdgId() const { return (fPdgId<0 ? -fPdgId:fPdgId); }
|
45 |
void AddDaughter(const MCParticle *p) { fDaughters.Add(p); }
|
46 |
const ThreeVector DecayVertex() const { return fDecayVertex.V(); }
|
47 |
const MCParticle *Daughter(UInt_t i) const;
|
48 |
const MCParticle *DistinctMother() const;
|
49 |
using CompositeParticle::HasDaughter;
|
50 |
const MCParticle *FindDaughter(Int_t pid,
|
51 |
Bool_t checkCharge=kFALSE, const MCParticle *start=0) const;
|
52 |
const MCParticle *FindMother(Int_t pid, Bool_t checkCharge=kFALSE) const;
|
53 |
Bool_t HasDaughter(Int_t pid, Bool_t checkCharge=kFALSE) const;
|
54 |
Bool_t HasMother() const { return fMother.IsValid(); }
|
55 |
Bool_t HasMother(const MCParticle *m) const;
|
56 |
Bool_t HasMother(Int_t pid, Bool_t checkCharge=kFALSE) const;
|
57 |
Bool_t Is(Int_t pid, Bool_t checkCharge=kFALSE) const;
|
58 |
Bool_t IsCharged() const { return !IsNeutral(); }
|
59 |
Bool_t IsGenerated() const { return fIsGenerated; }
|
60 |
Bool_t IsGluon() const { return fPdgId == kGlu; }
|
61 |
Bool_t IsLepton() const;
|
62 |
Bool_t IsNeutral() const;
|
63 |
Bool_t IsNeutrino() const;
|
64 |
Bool_t IsNot(Int_t pid, Bool_t checkCharge=kFALSE) const;
|
65 |
Bool_t IsParton() const { return (IsGluon() || IsQuark()); }
|
66 |
Bool_t IsQuark() const { return (AbsPdgId()>0 && AbsPdgId()<7); }
|
67 |
Bool_t IsSimulated() const { return fIsSimulated; }
|
68 |
const MCParticle *Mother() const { return fMother.Obj(); }
|
69 |
EObjType ObjType() const { return kMCParticle; }
|
70 |
void SetIsGenerated(Bool_t t=kTRUE) { fIsGenerated = t; }
|
71 |
void SetIsSimulated(Bool_t t=kTRUE) { fIsSimulated = t; }
|
72 |
TParticlePDG *PdgEntry() const;
|
73 |
Int_t PdgId() const { return fPdgId; }
|
74 |
Double_t PdgMass() const;
|
75 |
void SetMom(Double_t px, Double_t py, Double_t pz, Double_t e);
|
76 |
void SetMother(const MCParticle *p) { fMother = p; }
|
77 |
void SetStatus(Int_t s) { fStatus = s; }
|
78 |
void SetVertex(Double_t x, Double_t y, Double_t z);
|
79 |
void SetPdgId(Int_t s) { fPdgId = s; }
|
80 |
Int_t Status() const { return fStatus; }
|
81 |
void Print(Option_t *opt="") const;
|
82 |
|
83 |
// Some structural tools
|
84 |
void MarkMother(UInt_t i=1) const;
|
85 |
|
86 |
protected:
|
87 |
Double_t GetCharge() const;
|
88 |
void GetMom() const;
|
89 |
|
90 |
Int_t fPdgId; //pdg identifier
|
91 |
Short_t fStatus; //status flag of generator or simulation
|
92 |
Vect4M fMom; //four momentum vector
|
93 |
Vect3 fDecayVertex; //gen decay vertex
|
94 |
Ref<MCParticle> fMother; //reference to mother
|
95 |
Bool_t fIsGenerated; //=true if generated particle
|
96 |
Bool_t fIsSimulated; //=true if simulated particle
|
97 |
|
98 |
ClassDef(MCParticle,1) // Generated particle class
|
99 |
};
|
100 |
}
|
101 |
|
102 |
//--------------------------------------------------------------------------------------------------
|
103 |
inline void mithep::MCParticle::MarkMother(UInt_t ib) const
|
104 |
{
|
105 |
// mark myself
|
106 |
mithep::DataObject::Mark(ib);
|
107 |
// mark only the mother reference (as opposed to standard Mark() function which marks
|
108 |
// only daughters, keep them seperate to avoid circular references
|
109 |
if (fMother.IsValid())
|
110 |
fMother.Obj()->Mark(ib);
|
111 |
|
112 |
}
|
113 |
|
114 |
//--------------------------------------------------------------------------------------------------
|
115 |
inline const mithep::MCParticle *mithep::MCParticle::Daughter(UInt_t i) const
|
116 |
{
|
117 |
// Return daughter corresponding to given index.
|
118 |
|
119 |
return static_cast<const MCParticle*>(fDaughters.At(i));
|
120 |
}
|
121 |
|
122 |
//--------------------------------------------------------------------------------------------------
|
123 |
inline const mithep::MCParticle *mithep::MCParticle::DistinctMother() const
|
124 |
{
|
125 |
// Return mother, walking up the tree until a particle with a different pdg from this one
|
126 |
// is found.
|
127 |
|
128 |
const mithep::MCParticle *mother = Mother();
|
129 |
|
130 |
if (!mother)
|
131 |
return 0;
|
132 |
|
133 |
while (mother && mother->PdgId()==fPdgId)
|
134 |
mother = mother->Mother();
|
135 |
|
136 |
return mother;
|
137 |
}
|
138 |
|
139 |
//--------------------------------------------------------------------------------------------------
|
140 |
inline Double_t mithep::MCParticle::GetCharge() const
|
141 |
{
|
142 |
// Get charge from pdg lookup.
|
143 |
|
144 |
TParticlePDG *pdgEntry = PdgEntry();
|
145 |
if (pdgEntry)
|
146 |
return pdgEntry->Charge()/3.0;
|
147 |
else {
|
148 |
Error("Charge", "Pdg code for %i not found in table, returning charge=-99.0", fPdgId);
|
149 |
return -99.0;
|
150 |
}
|
151 |
}
|
152 |
|
153 |
//--------------------------------------------------------------------------------------------------
|
154 |
inline void mithep::MCParticle::GetMom() const
|
155 |
{
|
156 |
// Get momentum values from stored values.
|
157 |
|
158 |
fCachedMom.SetCoordinates(fMom.Pt(), fMom.Eta(), fMom.Phi(), fMom.M());
|
159 |
}
|
160 |
|
161 |
//--------------------------------------------------------------------------------------------------
|
162 |
inline Bool_t mithep::MCParticle::HasDaughter(Int_t pid, Bool_t checkCharge) const
|
163 |
{
|
164 |
// Return true if a particle with given pdg code is found amoung daughters.
|
165 |
// If checkCharge is false then just the type of particle is checked
|
166 |
// (ie particle and anti-particle).
|
167 |
|
168 |
if (checkCharge) {
|
169 |
for (UInt_t i=0; i<NDaughters(); ++i)
|
170 |
if (Daughter(i)->PdgId()==pid)
|
171 |
return kTRUE;
|
172 |
} else {
|
173 |
Int_t apid = pid>0?pid:-pid;
|
174 |
for (UInt_t i=0; i<NDaughters(); ++i)
|
175 |
if (Daughter(i)->AbsPdgId()==apid)
|
176 |
return kTRUE;
|
177 |
}
|
178 |
return kFALSE;
|
179 |
}
|
180 |
|
181 |
//--------------------------------------------------------------------------------------------------
|
182 |
inline Bool_t mithep::MCParticle::HasMother(const MCParticle *m) const
|
183 |
{
|
184 |
// Return true if the given particle is among mothers. (Note the comparison
|
185 |
// is made on pointers and therefore will fail if you work on copies.)
|
186 |
|
187 |
if (!m)
|
188 |
return kFALSE;
|
189 |
|
190 |
const mithep::MCParticle *mother = Mother();
|
191 |
while (mother && mother!=m)
|
192 |
mother = mother->Mother();
|
193 |
|
194 |
if (mother)
|
195 |
return kTRUE;
|
196 |
return kFALSE;
|
197 |
}
|
198 |
|
199 |
//--------------------------------------------------------------------------------------------------
|
200 |
inline Bool_t mithep::MCParticle::HasMother(Int_t pid, Bool_t checkCharge) const
|
201 |
{
|
202 |
// Return true if a particle with given pdg code is found amoung mothers.
|
203 |
// If checkCharge is false then just the type of particle is checked
|
204 |
// (ie particle and anti-particle).
|
205 |
|
206 |
const mithep::MCParticle *mother = Mother();
|
207 |
|
208 |
if (checkCharge) {
|
209 |
while (mother && mother->PdgId()!=pid)
|
210 |
mother = mother->Mother();
|
211 |
} else {
|
212 |
Int_t apid = pid>0?pid:-pid;
|
213 |
while (mother && mother->AbsPdgId()!=apid)
|
214 |
mother = mother->Mother();
|
215 |
}
|
216 |
|
217 |
if (mother)
|
218 |
return kTRUE;
|
219 |
return kFALSE;
|
220 |
}
|
221 |
|
222 |
//--------------------------------------------------------------------------------------------------
|
223 |
inline Bool_t mithep::MCParticle::Is(Int_t pid, Bool_t checkCharge) const
|
224 |
{
|
225 |
// Return true if particle is of given type. If checkCharge is false then just the type of
|
226 |
// particle is checked (ie particle and anti-particle).
|
227 |
|
228 |
if (checkCharge)
|
229 |
return (PdgId() == pid);
|
230 |
|
231 |
Int_t apid = pid>0?pid:-pid;
|
232 |
return (AbsPdgId() == apid);
|
233 |
}
|
234 |
|
235 |
//--------------------------------------------------------------------------------------------------
|
236 |
inline Bool_t mithep::MCParticle::IsNot(Int_t pid, Bool_t checkCharge) const
|
237 |
{
|
238 |
// Return true if particle is not of given type. If checkCharge is false then just the type of
|
239 |
// particle is checked (ie particle and anti-particle).
|
240 |
|
241 |
return !Is(pid, checkCharge);
|
242 |
}
|
243 |
|
244 |
//--------------------------------------------------------------------------------------------------
|
245 |
inline Bool_t mithep::MCParticle::IsLepton() const
|
246 |
{
|
247 |
// Return true if particle is a lepton.
|
248 |
|
249 |
Int_t ap = AbsPdgId();
|
250 |
if ((ap==kEl) ||
|
251 |
(ap==kMu) ||
|
252 |
(ap==kTau))
|
253 |
return kTRUE;
|
254 |
|
255 |
return kFALSE;
|
256 |
}
|
257 |
|
258 |
//--------------------------------------------------------------------------------------------------
|
259 |
inline Bool_t mithep::MCParticle::IsNeutral() const
|
260 |
{
|
261 |
// Return true if particle is a neutral.
|
262 |
|
263 |
return (Charge()!=0);
|
264 |
}
|
265 |
|
266 |
//--------------------------------------------------------------------------------------------------
|
267 |
inline Bool_t mithep::MCParticle::IsNeutrino() const
|
268 |
{
|
269 |
// Return true if particle is a neutrino.
|
270 |
|
271 |
Int_t ap = AbsPdgId();
|
272 |
if ((ap==kElNu) ||
|
273 |
(ap==kMuNu) ||
|
274 |
(ap==kTauNu))
|
275 |
return kTRUE;
|
276 |
|
277 |
return kFALSE;
|
278 |
}
|
279 |
|
280 |
//--------------------------------------------------------------------------------------------------
|
281 |
inline TParticlePDG *mithep::MCParticle::PdgEntry() const
|
282 |
{
|
283 |
// Return entry to pdg database.
|
284 |
|
285 |
return TDatabasePDG::Instance()->GetParticle(fPdgId);
|
286 |
}
|
287 |
|
288 |
//--------------------------------------------------------------------------------------------------
|
289 |
inline Double_t mithep::MCParticle::PdgMass() const
|
290 |
{
|
291 |
// Return mass obtained from pdg database.
|
292 |
|
293 |
TParticlePDG *pdg = PdgEntry();
|
294 |
if (pdg)
|
295 |
return pdg->Mass();
|
296 |
return 0;
|
297 |
}
|
298 |
|
299 |
//--------------------------------------------------------------------------------------------------
|
300 |
inline void mithep::MCParticle::SetMom(Double_t px, Double_t py, Double_t pz, Double_t e)
|
301 |
{
|
302 |
// Set four vector.
|
303 |
|
304 |
fMom.SetXYZT(px, py, pz, e);
|
305 |
ClearMom();
|
306 |
}
|
307 |
|
308 |
//--------------------------------------------------------------------------------------------------
|
309 |
inline void mithep::MCParticle::SetVertex(Double_t x, Double_t y, Double_t z)
|
310 |
{
|
311 |
// Set decay vertex.
|
312 |
|
313 |
fDecayVertex.SetXYZ(x,y,z);
|
314 |
}
|
315 |
#endif
|