Fork me on GitHub

source: git/modules/TreeWriter.cc @ 4d7014e

ImprovedOutputFileTiming
Last change on this file since 4d7014e was 4d7014e, checked in by Michele Selvaggi <michele.selvaggi@…>, 12 months ago

added PFcandidate class

  • Property mode set to 100644
File size: 26.7 KB
Line 
1/*
2 *  Delphes: a framework for fast simulation of a generic collider experiment
3 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
4 *
5 *  This program is free software: you can redistribute it and/or modify
6 *  it under the terms of the GNU General Public License as published by
7 *  the Free Software Foundation, either version 3 of the License, or
8 *  (at your option) any later version.
9 *
10 *  This program is distributed in the hope that it will be useful,
11 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 *  GNU General Public License for more details.
14 *
15 *  You should have received a copy of the GNU General Public License
16 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/** \class TreeWriter
20 *
21 *  Fills ROOT tree branches.
22 *
23 *  \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "modules/TreeWriter.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesFormula.h"
32
33#include "ExRootAnalysis/ExRootClassifier.h"
34#include "ExRootAnalysis/ExRootFilter.h"
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootTreeBranch.h"
37
38#include "TDatabasePDG.h"
39#include "TFormula.h"
40#include "TLorentzVector.h"
41#include "TMath.h"
42#include "TObjArray.h"
43#include "TROOT.h"
44#include "TRandom3.h"
45#include "TString.h"
46
47#include <algorithm>
48#include <iostream>
49#include <sstream>
50#include <stdexcept>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56TreeWriter::TreeWriter()
57{
58}
59
60//------------------------------------------------------------------------------
61
62TreeWriter::~TreeWriter()
63{
64}
65
66//------------------------------------------------------------------------------
67
68void TreeWriter::Init()
69{
70  fClassMap[GenParticle::Class()] = &TreeWriter::ProcessParticles;
71  fClassMap[Vertex::Class()] = &TreeWriter::ProcessVertices;
72  fClassMap[Track::Class()] = &TreeWriter::ProcessTracks;
73  fClassMap[Tower::Class()] = &TreeWriter::ProcessTowers;
74  fClassMap[ParticleFlowCandidate::Class()] = &TreeWriter::ProcessParticleFlowCandidates;
75  fClassMap[Photon::Class()] = &TreeWriter::ProcessPhotons;
76  fClassMap[Electron::Class()] = &TreeWriter::ProcessElectrons;
77  fClassMap[Muon::Class()] = &TreeWriter::ProcessMuons;
78  fClassMap[Jet::Class()] = &TreeWriter::ProcessJets;
79  fClassMap[MissingET::Class()] = &TreeWriter::ProcessMissingET;
80  fClassMap[ScalarHT::Class()] = &TreeWriter::ProcessScalarHT;
81  fClassMap[Rho::Class()] = &TreeWriter::ProcessRho;
82  fClassMap[Weight::Class()] = &TreeWriter::ProcessWeight;
83  fClassMap[HectorHit::Class()] = &TreeWriter::ProcessHectorHit;
84
85  TBranchMap::iterator itBranchMap;
86  map<TClass *, TProcessMethod>::iterator itClassMap;
87
88  // read branch configuration and
89  // import array with output from filter/classifier/jetfinder modules
90
91  ExRootConfParam param = GetParam("Branch");
92  Long_t i, size;
93  TString branchName, branchClassName, branchInputArray;
94  TClass *branchClass;
95  TObjArray *array;
96  ExRootTreeBranch *branch;
97
98  size = param.GetSize();
99  for(i = 0; i < size / 3; ++i)
100  {
101    branchInputArray = param[i * 3].GetString();
102    branchName = param[i * 3 + 1].GetString();
103    branchClassName = param[i * 3 + 2].GetString();
104
105    branchClass = gROOT->GetClass(branchClassName);
106
107    if(!branchClass)
108    {
109      cout << "** ERROR: cannot find class '" << branchClassName << "'" << endl;
110      continue;
111    }
112
113    itClassMap = fClassMap.find(branchClass);
114    if(itClassMap == fClassMap.end())
115    {
116      cout << "** ERROR: cannot create branch for class '" << branchClassName << "'" << endl;
117      continue;
118    }
119
120    array = ImportArray(branchInputArray);
121    branch = NewBranch(branchName, branchClass);
122
123    fBranchMap.insert(make_pair(branch, make_pair(itClassMap->second, array)));
124  }
125}
126
127//------------------------------------------------------------------------------
128
129void TreeWriter::Finish()
130{
131}
132
133//------------------------------------------------------------------------------
134
135void TreeWriter::FillParticles(Candidate *candidate, TRefArray *array)
136{
137  TIter it1(candidate->GetCandidates());
138  it1.Reset();
139  array->Clear();
140 
141  while((candidate = static_cast<Candidate *>(it1.Next())))
142  {
143    TIter it2(candidate->GetCandidates());
144
145    // particle
146    if(candidate->GetCandidates()->GetEntriesFast() == 0)
147    {
148      array->Add(candidate);
149      continue;
150    }
151
152    // track
153    candidate = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
154    if(candidate->GetCandidates()->GetEntriesFast() == 0)
155    {
156      array->Add(candidate);
157      continue;
158    }
159
160    // tower
161    it2.Reset();
162    while((candidate = static_cast<Candidate *>(it2.Next())))
163    {
164      array->Add(candidate->GetCandidates()->At(0));
165    }
166  }
167}
168
169//------------------------------------------------------------------------------
170
171void TreeWriter::ProcessParticles(ExRootTreeBranch *branch, TObjArray *array)
172{
173  TIter iterator(array);
174  Candidate *candidate = 0;
175  GenParticle *entry = 0;
176  Double_t pt, signPz, cosTheta, eta, rapidity;
177
178  const Double_t c_light = 2.99792458E8;
179
180  // loop over all particles
181  iterator.Reset();
182  while((candidate = static_cast<Candidate *>(iterator.Next())))
183  {
184    const TLorentzVector &momentum = candidate->Momentum;
185    const TLorentzVector &position = candidate->Position;
186
187    entry = static_cast<GenParticle *>(branch->NewEntry());
188
189    entry->SetBit(kIsReferenced);
190    entry->SetUniqueID(candidate->GetUniqueID());
191
192    pt = momentum.Pt();
193    cosTheta = TMath::Abs(momentum.CosTheta());
194    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
195    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
196    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
197
198    entry->PID = candidate->PID;
199
200    entry->Status = candidate->Status;
201    entry->IsPU = candidate->IsPU;
202
203    entry->M1 = candidate->M1;
204    entry->M2 = candidate->M2;
205
206    entry->D1 = candidate->D1;
207    entry->D2 = candidate->D2;
208
209    entry->Charge = candidate->Charge;
210    entry->Mass = candidate->Mass;
211
212    entry->E = momentum.E();
213    entry->Px = momentum.Px();
214    entry->Py = momentum.Py();
215    entry->Pz = momentum.Pz();
216
217    entry->D0 = candidate->D0;
218    entry->DZ = candidate->DZ;
219    entry->P = candidate->P;
220    entry->PT = candidate->PT;
221    entry->CtgTheta = candidate->CtgTheta;
222    entry->Phi = candidate->Phi;
223
224    entry->Eta = eta;
225    entry->Phi = momentum.Phi();
226    entry->PT = pt;
227
228    entry->Rapidity = rapidity;
229
230    entry->X = position.X();
231    entry->Y = position.Y();
232    entry->Z = position.Z();
233    entry->T = position.T() * 1.0E-3 / c_light;
234  }
235}
236
237//------------------------------------------------------------------------------
238
239void TreeWriter::ProcessVertices(ExRootTreeBranch *branch, TObjArray *array)
240{
241  TIter iterator(array);
242  Candidate *candidate = 0, *constituent = 0;
243  Vertex *entry = 0;
244
245  const Double_t c_light = 2.99792458E8;
246
247  Double_t x, y, z, t, xError, yError, zError, tError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
248  UInt_t index, ndf;
249
250  CompBase *compare = Candidate::fgCompare;
251  Candidate::fgCompare = CompSumPT2<Candidate>::Instance();
252  array->Sort();
253  Candidate::fgCompare = compare;
254
255  // loop over all vertices
256  iterator.Reset();
257  while((candidate = static_cast<Candidate *>(iterator.Next())))
258  {
259
260    index = candidate->ClusterIndex;
261    ndf = candidate->ClusterNDF;
262    sigma = candidate->ClusterSigma;
263    sumPT2 = candidate->SumPT2;
264    btvSumPT2 = candidate->BTVSumPT2;
265    genDeltaZ = candidate->GenDeltaZ;
266    genSumPT2 = candidate->GenSumPT2;
267
268    x = candidate->Position.X();
269    y = candidate->Position.Y();
270    z = candidate->Position.Z();
271    t = candidate->Position.T() * 1.0E-3 / c_light;
272
273    xError = candidate->PositionError.X();
274    yError = candidate->PositionError.Y();
275    zError = candidate->PositionError.Z();
276    tError = candidate->PositionError.T() * 1.0E-3 / c_light;
277
278    entry = static_cast<Vertex *>(branch->NewEntry());
279
280    entry->Index = index;
281    entry->NDF = ndf;
282    entry->Sigma = sigma;
283    entry->SumPT2 = sumPT2;
284    entry->BTVSumPT2 = btvSumPT2;
285    entry->GenDeltaZ = genDeltaZ;
286    entry->GenSumPT2 = genSumPT2;
287
288    entry->X = x;
289    entry->Y = y;
290    entry->Z = z;
291    entry->T = t;
292
293    entry->ErrorX = xError;
294    entry->ErrorY = yError;
295    entry->ErrorZ = zError;
296    entry->ErrorT = tError;
297
298    TIter itConstituents(candidate->GetCandidates());
299    itConstituents.Reset();
300    entry->Constituents.Clear();
301    while((constituent = static_cast<Candidate *>(itConstituents.Next())))
302    {
303      entry->Constituents.Add(constituent);
304    }
305  }
306}
307
308//------------------------------------------------------------------------------
309
310void TreeWriter::ProcessTracks(ExRootTreeBranch *branch, TObjArray *array)
311{
312  TIter iterator(array);
313  Candidate *candidate = 0;
314  Candidate *particle = 0;
315  Track *entry = 0;
316  Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
317  const Double_t c_light = 2.99792458E8;
318
319  // loop over all tracks
320  iterator.Reset();
321  while((candidate = static_cast<Candidate *>(iterator.Next())))
322  {
323    const TLorentzVector &position = candidate->Position;
324
325    cosTheta = TMath::Abs(position.CosTheta());
326    signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
327    eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
328    rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
329
330    entry = static_cast<Track *>(branch->NewEntry());
331
332    entry->SetBit(kIsReferenced);
333    entry->SetUniqueID(candidate->GetUniqueID());
334
335    entry->PID = candidate->PID;
336
337    entry->Charge = candidate->Charge;
338
339    entry->EtaOuter = eta;
340    entry->PhiOuter = position.Phi();
341
342    entry->XOuter = position.X();
343    entry->YOuter = position.Y();
344    entry->ZOuter = position.Z();
345    entry->TOuter = position.T() * 1.0E-3 / c_light;
346
347    entry->L = candidate->L;
348
349    entry->D0 = candidate->D0;
350    entry->ErrorD0 = candidate->ErrorD0;
351    entry->DZ = candidate->DZ;
352    entry->ErrorDZ = candidate->ErrorDZ;
353
354    entry->ErrorP = candidate->ErrorP;
355    entry->ErrorPT = candidate->ErrorPT;
356    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
357    entry->ErrorPhi = candidate->ErrorPhi;
358
359    entry->Xd = candidate->Xd;
360    entry->Yd = candidate->Yd;
361    entry->Zd = candidate->Zd;
362
363    const TLorentzVector &momentum = candidate->Momentum;
364
365    pt = momentum.Pt();
366    p = momentum.P();
367    phi = momentum.Phi();
368    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
369
370    cosTheta = TMath::Abs(momentum.CosTheta());
371    signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
372    eta = (cosTheta == 1.0 ? signz * 999.9 : momentum.Eta());
373    rapidity = (cosTheta == 1.0 ? signz * 999.9 : momentum.Rapidity());
374
375    entry->P = p;
376    entry->PT = pt;
377    entry->Eta = eta;
378    entry->Phi = phi;
379    entry->CtgTheta = ctgTheta;
380
381    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
382    const TLorentzVector &initialPosition = particle->Position;
383
384    entry->X = initialPosition.X();
385    entry->Y = initialPosition.Y();
386    entry->Z = initialPosition.Z();
387    entry->T = initialPosition.T() * 1.0E-3 / c_light;
388
389    entry->Particle = particle;
390
391    entry->VertexIndex = candidate->ClusterIndex;
392  }
393}
394
395//------------------------------------------------------------------------------
396
397void TreeWriter::ProcessTowers(ExRootTreeBranch *branch, TObjArray *array)
398{
399  TIter iterator(array);
400  Candidate *candidate = 0;
401  Tower *entry = 0;
402  Double_t pt, signPz, cosTheta, eta, rapidity;
403  const Double_t c_light = 2.99792458E8;
404
405  // loop over all towers
406  iterator.Reset();
407  while((candidate = static_cast<Candidate *>(iterator.Next())))
408  {
409    const TLorentzVector &momentum = candidate->Momentum;
410    const TLorentzVector &position = candidate->Position;
411
412    pt = momentum.Pt();
413    cosTheta = TMath::Abs(momentum.CosTheta());
414    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
415    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
416    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
417
418    entry = static_cast<Tower *>(branch->NewEntry());
419
420    entry->SetBit(kIsReferenced);
421    entry->SetUniqueID(candidate->GetUniqueID());
422
423    entry->Eta = eta;
424    entry->Phi = momentum.Phi();
425    entry->ET = pt;
426    entry->E = momentum.E();
427    entry->Eem = candidate->Eem;
428    entry->Ehad = candidate->Ehad;
429    entry->Edges[0] = candidate->Edges[0];
430    entry->Edges[1] = candidate->Edges[1];
431    entry->Edges[2] = candidate->Edges[2];
432    entry->Edges[3] = candidate->Edges[3];
433
434    entry->T = position.T() * 1.0E-3 / c_light;
435    entry->NTimeHits = candidate->NTimeHits;
436
437    FillParticles(candidate, &entry->Particles);
438  }
439}
440
441//------------------------------------------------------------------------------
442
443void TreeWriter::ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array)
444{
445
446  TIter iterator(array);
447  Candidate *candidate = 0;
448  Candidate *particle = 0;
449  ParticleFlowCandidate *entry = 0;
450  Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
451  const Double_t c_light = 2.99792458E8;
452
453  // loop over all tracks
454  iterator.Reset();
455  while((candidate = static_cast<Candidate *>(iterator.Next())))
456  {
457    const TLorentzVector &position = candidate->Position;
458
459    cosTheta = TMath::Abs(position.CosTheta());
460    signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
461    eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
462    rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
463
464    entry = static_cast<ParticleFlowCandidate *>(branch->NewEntry());
465
466    entry->SetBit(kIsReferenced);
467    entry->SetUniqueID(candidate->GetUniqueID());
468
469    entry->PID = candidate->PID;
470
471    entry->Charge = candidate->Charge;
472
473    entry->EtaOuter = eta;
474    entry->PhiOuter = position.Phi();
475
476    entry->XOuter = position.X();
477    entry->YOuter = position.Y();
478    entry->ZOuter = position.Z();
479    entry->TOuter = position.T() * 1.0E-3 / c_light;
480
481    entry->L = candidate->L;
482
483    entry->D0 = candidate->D0;
484    entry->ErrorD0 = candidate->ErrorD0;
485    entry->DZ = candidate->DZ;
486    entry->ErrorDZ = candidate->ErrorDZ;
487
488    entry->ErrorP = candidate->ErrorP;
489    entry->ErrorPT = candidate->ErrorPT;
490    entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
491    entry->ErrorPhi = candidate->ErrorPhi;
492
493    entry->Xd = candidate->Xd;
494    entry->Yd = candidate->Yd;
495    entry->Zd = candidate->Zd;
496
497    const TLorentzVector &momentum = candidate->Momentum;
498
499    e = momentum.E();
500    pt = momentum.Pt();
501    p = momentum.P();
502    phi = momentum.Phi();
503    ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
504
505    entry->E = e;
506    entry->P = p;
507    entry->PT = pt;
508    entry->Eta = eta;
509    entry->Phi = phi;
510    entry->CtgTheta = ctgTheta;
511
512    particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
513    const TLorentzVector &initialPosition = particle->Position;
514
515    entry->X = initialPosition.X();
516    entry->Y = initialPosition.Y();
517    entry->Z = initialPosition.Z();
518    entry->T = initialPosition.T() * 1.0E-3 / c_light;
519
520    entry->VertexIndex = candidate->ClusterIndex;
521
522    entry->Eem = candidate->Eem;
523    entry->Ehad = candidate->Ehad;
524    entry->Edges[0] = candidate->Edges[0];
525    entry->Edges[1] = candidate->Edges[1];
526    entry->Edges[2] = candidate->Edges[2];
527    entry->Edges[3] = candidate->Edges[3];
528
529    entry->T = position.T() * 1.0E-3 / c_light;
530    entry->NTimeHits = candidate->NTimeHits;
531
532    FillParticles(candidate, &entry->Particles);
533
534  }
535}
536
537//------------------------------------------------------------------------------
538
539void TreeWriter::ProcessPhotons(ExRootTreeBranch *branch, TObjArray *array)
540{
541  TIter iterator(array);
542  Candidate *candidate = 0;
543  Photon *entry = 0;
544  Double_t pt, signPz, cosTheta, eta, rapidity;
545  const Double_t c_light = 2.99792458E8;
546
547  array->Sort();
548
549  // loop over all photons
550  iterator.Reset();
551  while((candidate = static_cast<Candidate *>(iterator.Next())))
552  {
553    TIter it1(candidate->GetCandidates());
554    const TLorentzVector &momentum = candidate->Momentum;
555    const TLorentzVector &position = candidate->Position;
556
557    pt = momentum.Pt();
558    cosTheta = TMath::Abs(momentum.CosTheta());
559    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
560    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
561    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
562
563    entry = static_cast<Photon *>(branch->NewEntry());
564
565    entry->Eta = eta;
566    entry->Phi = momentum.Phi();
567    entry->PT = pt;
568    entry->E = momentum.E();
569    entry->T = position.T() * 1.0E-3 / c_light;
570
571    // Isolation variables
572
573    entry->IsolationVar = candidate->IsolationVar;
574    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
575    entry->SumPtCharged = candidate->SumPtCharged;
576    entry->SumPtNeutral = candidate->SumPtNeutral;
577    entry->SumPtChargedPU = candidate->SumPtChargedPU;
578    entry->SumPt = candidate->SumPt;
579
580    entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad / candidate->Eem : 999.9;
581
582    // 1: prompt -- 2: non prompt -- 3: fake
583    entry->Status = candidate->Status;
584
585    FillParticles(candidate, &entry->Particles);
586  }
587}
588
589//------------------------------------------------------------------------------
590
591void TreeWriter::ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array)
592{
593  TIter iterator(array);
594  Candidate *candidate = 0;
595  Electron *entry = 0;
596  Double_t pt, signPz, cosTheta, eta, rapidity;
597  const Double_t c_light = 2.99792458E8;
598
599  array->Sort();
600
601  // loop over all electrons
602  iterator.Reset();
603  while((candidate = static_cast<Candidate *>(iterator.Next())))
604  {
605    const TLorentzVector &momentum = candidate->Momentum;
606    const TLorentzVector &position = candidate->Position;
607
608    pt = momentum.Pt();
609    cosTheta = TMath::Abs(momentum.CosTheta());
610    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
611    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
612    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
613
614    entry = static_cast<Electron *>(branch->NewEntry());
615
616    entry->Eta = eta;
617    entry->Phi = momentum.Phi();
618    entry->PT = pt;
619
620    entry->T = position.T() * 1.0E-3 / c_light;
621
622    // displacement
623    entry->D0 = candidate->D0;
624    entry->ErrorD0 = candidate->ErrorD0;
625    entry->DZ = candidate->DZ;
626    entry->ErrorDZ = candidate->ErrorDZ;
627
628    // Isolation variables
629    entry->IsolationVar = candidate->IsolationVar;
630    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
631    entry->SumPtCharged = candidate->SumPtCharged;
632    entry->SumPtNeutral = candidate->SumPtNeutral;
633    entry->SumPtChargedPU = candidate->SumPtChargedPU;
634    entry->SumPt = candidate->SumPt;
635
636    entry->Charge = candidate->Charge;
637
638    entry->EhadOverEem = 0.0;
639
640    entry->Particle = candidate->GetCandidates()->At(0);
641  }
642}
643
644//------------------------------------------------------------------------------
645
646void TreeWriter::ProcessMuons(ExRootTreeBranch *branch, TObjArray *array)
647{
648  TIter iterator(array);
649  Candidate *candidate = 0;
650  Muon *entry = 0;
651  Double_t pt, signPz, cosTheta, eta, rapidity;
652
653  const Double_t c_light = 2.99792458E8;
654
655  array->Sort();
656
657  // loop over all muons
658  iterator.Reset();
659  while((candidate = static_cast<Candidate *>(iterator.Next())))
660  {
661    const TLorentzVector &momentum = candidate->Momentum;
662    const TLorentzVector &position = candidate->Position;
663
664    pt = momentum.Pt();
665    cosTheta = TMath::Abs(momentum.CosTheta());
666    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
667    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
668    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
669
670    entry = static_cast<Muon *>(branch->NewEntry());
671
672    entry->SetBit(kIsReferenced);
673    entry->SetUniqueID(candidate->GetUniqueID());
674
675    entry->Eta = eta;
676    entry->Phi = momentum.Phi();
677    entry->PT = pt;
678
679    entry->T = position.T() * 1.0E-3 / c_light;
680
681    // displacement
682    entry->D0 = candidate->D0;
683    entry->ErrorD0 = candidate->ErrorD0;
684    entry->DZ = candidate->DZ;
685    entry->ErrorDZ = candidate->ErrorDZ;
686
687    // Isolation variables
688
689    entry->IsolationVar = candidate->IsolationVar;
690    entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
691    entry->SumPtCharged = candidate->SumPtCharged;
692    entry->SumPtNeutral = candidate->SumPtNeutral;
693    entry->SumPtChargedPU = candidate->SumPtChargedPU;
694    entry->SumPt = candidate->SumPt;
695
696    entry->Charge = candidate->Charge;
697
698    entry->Particle = candidate->GetCandidates()->At(0);
699  }
700}
701
702//------------------------------------------------------------------------------
703
704void TreeWriter::ProcessJets(ExRootTreeBranch *branch, TObjArray *array)
705{
706  TIter iterator(array);
707  Candidate *candidate = 0, *constituent = 0;
708  Jet *entry = 0;
709  Double_t pt, signPz, cosTheta, eta, rapidity;
710  Double_t ecalEnergy, hcalEnergy;
711  const Double_t c_light = 2.99792458E8;
712  Int_t i;
713
714  array->Sort();
715
716  // loop over all jets
717  iterator.Reset();
718  while((candidate = static_cast<Candidate *>(iterator.Next())))
719  {
720    TIter itConstituents(candidate->GetCandidates());
721
722    const TLorentzVector &momentum = candidate->Momentum;
723    const TLorentzVector &position = candidate->Position;
724
725    pt = momentum.Pt();
726    cosTheta = TMath::Abs(momentum.CosTheta());
727    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
728    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
729    rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
730
731    entry = static_cast<Jet *>(branch->NewEntry());
732
733    entry->Eta = eta;
734    entry->Phi = momentum.Phi();
735    entry->PT = pt;
736
737    entry->T = position.T() * 1.0E-3 / c_light;
738
739    entry->Mass = momentum.M();
740
741    entry->Area = candidate->Area;
742
743    entry->DeltaEta = candidate->DeltaEta;
744    entry->DeltaPhi = candidate->DeltaPhi;
745
746    entry->Flavor = candidate->Flavor;
747    entry->FlavorAlgo = candidate->FlavorAlgo;
748    entry->FlavorPhys = candidate->FlavorPhys;
749
750    entry->BTag = candidate->BTag;
751
752    entry->BTagAlgo = candidate->BTagAlgo;
753    entry->BTagPhys = candidate->BTagPhys;
754
755    entry->TauTag = candidate->TauTag;
756    entry->TauWeight = candidate->TauWeight;
757
758    entry->Charge = candidate->Charge;
759
760    itConstituents.Reset();
761    entry->Constituents.Clear();
762    ecalEnergy = 0.0;
763    hcalEnergy = 0.0;
764    while((constituent = static_cast<Candidate *>(itConstituents.Next())))
765    {
766      entry->Constituents.Add(constituent);
767      ecalEnergy += constituent->Eem;
768      hcalEnergy += constituent->Ehad;
769    }
770
771    entry->EhadOverEem = ecalEnergy > 0.0 ? hcalEnergy / ecalEnergy : 999.9;
772
773    //---   Pile-Up Jet ID variables ----
774
775    entry->NCharged = candidate->NCharged;
776    entry->NNeutrals = candidate->NNeutrals;
777    entry->Beta = candidate->Beta;
778    entry->BetaStar = candidate->BetaStar;
779    entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
780    entry->PTD = candidate->PTD;
781
782    //--- Sub-structure variables ----
783
784    entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
785    entry->NSubJetsPruned = candidate->NSubJetsPruned;
786    entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
787
788    entry->SoftDroppedJet = candidate->SoftDroppedJet;
789    entry->SoftDroppedSubJet1 = candidate->SoftDroppedSubJet1;
790    entry->SoftDroppedSubJet2 = candidate->SoftDroppedSubJet2;
791
792    for(i = 0; i < 5; i++)
793    {
794      entry->FracPt[i] = candidate->FracPt[i];
795      entry->Tau[i] = candidate->Tau[i];
796      entry->TrimmedP4[i] = candidate->TrimmedP4[i];
797      entry->PrunedP4[i] = candidate->PrunedP4[i];
798      entry->SoftDroppedP4[i] = candidate->SoftDroppedP4[i];
799    }
800
801    //--- exclusive clustering variables ---
802    entry->ExclYmerge23 = candidate->ExclYmerge23;
803    entry->ExclYmerge34 = candidate->ExclYmerge34;
804    entry->ExclYmerge45 = candidate->ExclYmerge45;
805    entry->ExclYmerge56 = candidate->ExclYmerge56;
806
807    FillParticles(candidate, &entry->Particles);
808  }
809}
810
811//------------------------------------------------------------------------------
812
813void TreeWriter::ProcessMissingET(ExRootTreeBranch *branch, TObjArray *array)
814{
815  Candidate *candidate = 0;
816  MissingET *entry = 0;
817
818  // get the first entry
819  if((candidate = static_cast<Candidate *>(array->At(0))))
820  {
821    const TLorentzVector &momentum = candidate->Momentum;
822
823    entry = static_cast<MissingET *>(branch->NewEntry());
824
825    entry->Eta = (-momentum).Eta();
826    entry->Phi = (-momentum).Phi();
827    entry->MET = momentum.Pt();
828  }
829}
830
831//------------------------------------------------------------------------------
832
833void TreeWriter::ProcessScalarHT(ExRootTreeBranch *branch, TObjArray *array)
834{
835  Candidate *candidate = 0;
836  ScalarHT *entry = 0;
837
838  // get the first entry
839  if((candidate = static_cast<Candidate *>(array->At(0))))
840  {
841    const TLorentzVector &momentum = candidate->Momentum;
842
843    entry = static_cast<ScalarHT *>(branch->NewEntry());
844
845    entry->HT = momentum.Pt();
846  }
847}
848
849//------------------------------------------------------------------------------
850
851void TreeWriter::ProcessRho(ExRootTreeBranch *branch, TObjArray *array)
852{
853  TIter iterator(array);
854  Candidate *candidate = 0;
855  Rho *entry = 0;
856
857  // loop over all rho
858  iterator.Reset();
859  while((candidate = static_cast<Candidate *>(iterator.Next())))
860  {
861    const TLorentzVector &momentum = candidate->Momentum;
862
863    entry = static_cast<Rho *>(branch->NewEntry());
864
865    entry->Rho = momentum.E();
866    entry->Edges[0] = candidate->Edges[0];
867    entry->Edges[1] = candidate->Edges[1];
868  }
869}
870
871//------------------------------------------------------------------------------
872
873void TreeWriter::ProcessWeight(ExRootTreeBranch *branch, TObjArray *array)
874{
875  Candidate *candidate = 0;
876  Weight *entry = 0;
877
878  // get the first entry
879  if((candidate = static_cast<Candidate *>(array->At(0))))
880  {
881    const TLorentzVector &momentum = candidate->Momentum;
882
883    entry = static_cast<Weight *>(branch->NewEntry());
884
885    entry->Weight = momentum.E();
886  }
887}
888
889//------------------------------------------------------------------------------
890
891void TreeWriter::ProcessHectorHit(ExRootTreeBranch *branch, TObjArray *array)
892{
893  TIter iterator(array);
894  Candidate *candidate = 0;
895  HectorHit *entry = 0;
896
897  // loop over all roman pot hits
898  iterator.Reset();
899  while((candidate = static_cast<Candidate *>(iterator.Next())))
900  {
901    const TLorentzVector &position = candidate->Position;
902    const TLorentzVector &momentum = candidate->Momentum;
903
904    entry = static_cast<HectorHit *>(branch->NewEntry());
905
906    entry->E = momentum.E();
907
908    entry->Tx = momentum.Px();
909    entry->Ty = momentum.Py();
910
911    entry->T = position.T();
912
913    entry->X = position.X();
914    entry->Y = position.Y();
915    entry->S = position.Z();
916
917    entry->Particle = candidate->GetCandidates()->At(0);
918  }
919}
920
921//------------------------------------------------------------------------------
922
923void TreeWriter::Process()
924{
925  TBranchMap::iterator itBranchMap;
926  ExRootTreeBranch *branch;
927  TProcessMethod method;
928  TObjArray *array;
929
930  for(itBranchMap = fBranchMap.begin(); itBranchMap != fBranchMap.end(); ++itBranchMap)
931  {
932    branch = itBranchMap->first;
933    method = itBranchMap->second.first;
934    array = itBranchMap->second.second;
935
936    (this->*method)(branch, array);
937  }
938}
939
940//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.