Fork me on GitHub

source: git/modules/TreeWriter.cc@ b4a3c55

Timing
Last change on this file since b4a3c55 was b4a3c55, checked in by Michele Selvaggi <michele.selvaggi@…>, 5 years ago

added energy loss module

  • Property mode set to 100644
File size: 27.3 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 entry->DeDx = candidate->DeDx; // in MeV/cm
349
350 entry->D0 = candidate->D0;
351 entry->ErrorD0 = candidate->ErrorD0;
352 entry->DZ = candidate->DZ;
353 entry->ErrorDZ = candidate->ErrorDZ;
354
355 entry->ErrorP = candidate->ErrorP;
356 entry->ErrorPT = candidate->ErrorPT;
357 entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
358 entry->ErrorPhi = candidate->ErrorPhi;
359
360 entry->Xd = candidate->Xd;
361 entry->Yd = candidate->Yd;
362 entry->Zd = candidate->Zd;
363 entry->Td = candidate->Td*1.0E-3/c_light;
364
365 if(candidate->ClusterIndex != -1)
366 {
367 entry->TOFreco = 1E-3*(candidate->Position.T() - candidate->InitialPosition.T())/c_light;
368 }
369 else
370 {
371 entry->TOFreco =-1E6;
372 }
373
374 const TLorentzVector &momentum = candidate->Momentum;
375
376 pt = momentum.Pt();
377 p = momentum.P();
378 phi = momentum.Phi();
379 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
380
381 cosTheta = TMath::Abs(momentum.CosTheta());
382 signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
383 eta = (cosTheta == 1.0 ? signz * 999.9 : momentum.Eta());
384 rapidity = (cosTheta == 1.0 ? signz * 999.9 : momentum.Rapidity());
385
386 entry->P = p;
387 entry->PT = pt;
388 entry->Eta = eta;
389 entry->Phi = phi;
390 entry->CtgTheta = ctgTheta;
391
392 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
393 entry->TOFgen = 1E-3*particle->L/(c_light*particle->Momentum.P()/particle->Momentum.E());
394
395 const TLorentzVector &initialPosition = particle->Position;
396
397 entry->X = initialPosition.X();
398 entry->Y = initialPosition.Y();
399 entry->Z = initialPosition.Z();
400 entry->T = initialPosition.T() * 1.0E-3 / c_light;
401
402 entry->Particle = particle;
403
404 entry->VertexIndex = candidate->ClusterIndex;
405 }
406}
407
408//------------------------------------------------------------------------------
409
410void TreeWriter::ProcessTowers(ExRootTreeBranch *branch, TObjArray *array)
411{
412 TIter iterator(array);
413 Candidate *candidate = 0;
414 Tower *entry = 0;
415 Double_t pt, signPz, cosTheta, eta, rapidity;
416 const Double_t c_light = 2.99792458E8;
417
418 // loop over all towers
419 iterator.Reset();
420 while((candidate = static_cast<Candidate *>(iterator.Next())))
421 {
422 const TLorentzVector &momentum = candidate->Momentum;
423 const TLorentzVector &position = candidate->Position;
424
425 pt = momentum.Pt();
426 cosTheta = TMath::Abs(momentum.CosTheta());
427 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
428 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
429 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
430
431 entry = static_cast<Tower *>(branch->NewEntry());
432
433 entry->SetBit(kIsReferenced);
434 entry->SetUniqueID(candidate->GetUniqueID());
435
436 entry->Eta = eta;
437 entry->Phi = momentum.Phi();
438 entry->ET = pt;
439 entry->E = momentum.E();
440 entry->Eem = candidate->Eem;
441 entry->Ehad = candidate->Ehad;
442 entry->Edges[0] = candidate->Edges[0];
443 entry->Edges[1] = candidate->Edges[1];
444 entry->Edges[2] = candidate->Edges[2];
445 entry->Edges[3] = candidate->Edges[3];
446
447 entry->T = position.T() * 1.0E-3 / c_light;
448 entry->NTimeHits = candidate->NTimeHits;
449
450 FillParticles(candidate, &entry->Particles);
451 }
452}
453
454//------------------------------------------------------------------------------
455
456void TreeWriter::ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array)
457{
458
459 TIter iterator(array);
460 Candidate *candidate = 0;
461 Candidate *particle = 0;
462 ParticleFlowCandidate *entry = 0;
463 Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi;
464 const Double_t c_light = 2.99792458E8;
465
466 // loop over all tracks
467 iterator.Reset();
468 while((candidate = static_cast<Candidate *>(iterator.Next())))
469 {
470 const TLorentzVector &position = candidate->Position;
471
472 cosTheta = TMath::Abs(position.CosTheta());
473 signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
474 eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
475 rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
476
477 entry = static_cast<ParticleFlowCandidate *>(branch->NewEntry());
478
479 entry->SetBit(kIsReferenced);
480 entry->SetUniqueID(candidate->GetUniqueID());
481
482 entry->PID = candidate->PID;
483
484 entry->Charge = candidate->Charge;
485
486 entry->EtaOuter = eta;
487 entry->PhiOuter = position.Phi();
488
489 entry->XOuter = position.X();
490 entry->YOuter = position.Y();
491 entry->ZOuter = position.Z();
492 entry->TOuter = position.T() * 1.0E-3 / c_light;
493
494 entry->L = candidate->L;
495 entry->DeDx = candidate->DeDx;
496
497 entry->D0 = candidate->D0;
498 entry->ErrorD0 = candidate->ErrorD0;
499 entry->DZ = candidate->DZ;
500 entry->ErrorDZ = candidate->ErrorDZ;
501
502 entry->ErrorP = candidate->ErrorP;
503 entry->ErrorPT = candidate->ErrorPT;
504 entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
505 entry->ErrorPhi = candidate->ErrorPhi;
506
507 entry->Xd = candidate->Xd;
508 entry->Yd = candidate->Yd;
509 entry->Zd = candidate->Zd;
510
511 const TLorentzVector &momentum = candidate->Momentum;
512
513 e = momentum.E();
514 pt = momentum.Pt();
515 p = momentum.P();
516 phi = momentum.Phi();
517 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
518
519 entry->E = e;
520 entry->P = p;
521 entry->PT = pt;
522 entry->Eta = eta;
523 entry->Phi = phi;
524 entry->CtgTheta = ctgTheta;
525
526 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
527 const TLorentzVector &initialPosition = particle->Position;
528
529 entry->X = initialPosition.X();
530 entry->Y = initialPosition.Y();
531 entry->Z = initialPosition.Z();
532 entry->T = initialPosition.T() * 1.0E-3 / c_light;
533
534 entry->VertexIndex = candidate->ClusterIndex;
535
536 entry->Eem = candidate->Eem;
537 entry->Ehad = candidate->Ehad;
538 entry->Edges[0] = candidate->Edges[0];
539 entry->Edges[1] = candidate->Edges[1];
540 entry->Edges[2] = candidate->Edges[2];
541 entry->Edges[3] = candidate->Edges[3];
542
543 entry->T = position.T() * 1.0E-3 / c_light;
544 entry->NTimeHits = candidate->NTimeHits;
545
546 FillParticles(candidate, &entry->Particles);
547
548 }
549}
550
551//------------------------------------------------------------------------------
552
553void TreeWriter::ProcessPhotons(ExRootTreeBranch *branch, TObjArray *array)
554{
555 TIter iterator(array);
556 Candidate *candidate = 0;
557 Photon *entry = 0;
558 Double_t pt, signPz, cosTheta, eta, rapidity;
559 const Double_t c_light = 2.99792458E8;
560
561 array->Sort();
562
563 // loop over all photons
564 iterator.Reset();
565 while((candidate = static_cast<Candidate *>(iterator.Next())))
566 {
567 TIter it1(candidate->GetCandidates());
568 const TLorentzVector &momentum = candidate->Momentum;
569 const TLorentzVector &position = candidate->Position;
570
571 pt = momentum.Pt();
572 cosTheta = TMath::Abs(momentum.CosTheta());
573 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
574 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
575 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
576
577 entry = static_cast<Photon *>(branch->NewEntry());
578
579 entry->Eta = eta;
580 entry->Phi = momentum.Phi();
581 entry->PT = pt;
582 entry->E = momentum.E();
583 entry->T = position.T() * 1.0E-3 / c_light;
584
585 // Isolation variables
586
587 entry->IsolationVar = candidate->IsolationVar;
588 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
589 entry->SumPtCharged = candidate->SumPtCharged;
590 entry->SumPtNeutral = candidate->SumPtNeutral;
591 entry->SumPtChargedPU = candidate->SumPtChargedPU;
592 entry->SumPt = candidate->SumPt;
593
594 entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad / candidate->Eem : 999.9;
595
596 // 1: prompt -- 2: non prompt -- 3: fake
597 entry->Status = candidate->Status;
598
599 FillParticles(candidate, &entry->Particles);
600 }
601}
602
603//------------------------------------------------------------------------------
604
605void TreeWriter::ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array)
606{
607 TIter iterator(array);
608 Candidate *candidate = 0;
609 Electron *entry = 0;
610 Double_t pt, signPz, cosTheta, eta, rapidity;
611 const Double_t c_light = 2.99792458E8;
612
613 array->Sort();
614
615 // loop over all electrons
616 iterator.Reset();
617 while((candidate = static_cast<Candidate *>(iterator.Next())))
618 {
619 const TLorentzVector &momentum = candidate->Momentum;
620 const TLorentzVector &position = candidate->Position;
621
622 pt = momentum.Pt();
623 cosTheta = TMath::Abs(momentum.CosTheta());
624 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
625 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
626 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
627
628 entry = static_cast<Electron *>(branch->NewEntry());
629
630 entry->Eta = eta;
631 entry->Phi = momentum.Phi();
632 entry->PT = pt;
633
634 entry->T = position.T() * 1.0E-3 / c_light;
635
636 // displacement
637 entry->D0 = candidate->D0;
638 entry->ErrorD0 = candidate->ErrorD0;
639 entry->DZ = candidate->DZ;
640 entry->ErrorDZ = candidate->ErrorDZ;
641
642 // Isolation variables
643 entry->IsolationVar = candidate->IsolationVar;
644 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
645 entry->SumPtCharged = candidate->SumPtCharged;
646 entry->SumPtNeutral = candidate->SumPtNeutral;
647 entry->SumPtChargedPU = candidate->SumPtChargedPU;
648 entry->SumPt = candidate->SumPt;
649
650 entry->Charge = candidate->Charge;
651
652 entry->EhadOverEem = 0.0;
653
654 entry->Particle = candidate->GetCandidates()->At(0);
655 }
656}
657
658//------------------------------------------------------------------------------
659
660void TreeWriter::ProcessMuons(ExRootTreeBranch *branch, TObjArray *array)
661{
662 TIter iterator(array);
663 Candidate *candidate = 0;
664 Muon *entry = 0;
665 Double_t pt, signPz, cosTheta, eta, rapidity;
666
667 const Double_t c_light = 2.99792458E8;
668
669 array->Sort();
670
671 // loop over all muons
672 iterator.Reset();
673 while((candidate = static_cast<Candidate *>(iterator.Next())))
674 {
675 const TLorentzVector &momentum = candidate->Momentum;
676 const TLorentzVector &position = candidate->Position;
677
678 pt = momentum.Pt();
679 cosTheta = TMath::Abs(momentum.CosTheta());
680 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
681 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
682 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
683
684 entry = static_cast<Muon *>(branch->NewEntry());
685
686 entry->SetBit(kIsReferenced);
687 entry->SetUniqueID(candidate->GetUniqueID());
688
689 entry->Eta = eta;
690 entry->Phi = momentum.Phi();
691 entry->PT = pt;
692
693 entry->T = position.T() * 1.0E-3 / c_light;
694
695 // displacement
696 entry->D0 = candidate->D0;
697 entry->ErrorD0 = candidate->ErrorD0;
698 entry->DZ = candidate->DZ;
699 entry->ErrorDZ = candidate->ErrorDZ;
700
701 // Isolation variables
702
703 entry->IsolationVar = candidate->IsolationVar;
704 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
705 entry->SumPtCharged = candidate->SumPtCharged;
706 entry->SumPtNeutral = candidate->SumPtNeutral;
707 entry->SumPtChargedPU = candidate->SumPtChargedPU;
708 entry->SumPt = candidate->SumPt;
709
710 entry->Charge = candidate->Charge;
711
712 entry->Particle = candidate->GetCandidates()->At(0);
713 }
714}
715
716//------------------------------------------------------------------------------
717
718void TreeWriter::ProcessJets(ExRootTreeBranch *branch, TObjArray *array)
719{
720 TIter iterator(array);
721 Candidate *candidate = 0, *constituent = 0;
722 Jet *entry = 0;
723 Double_t pt, signPz, cosTheta, eta, rapidity;
724 Double_t ecalEnergy, hcalEnergy;
725 const Double_t c_light = 2.99792458E8;
726 Int_t i;
727
728 array->Sort();
729
730 // loop over all jets
731 iterator.Reset();
732 while((candidate = static_cast<Candidate *>(iterator.Next())))
733 {
734 TIter itConstituents(candidate->GetCandidates());
735
736 const TLorentzVector &momentum = candidate->Momentum;
737 const TLorentzVector &position = candidate->Position;
738
739 pt = momentum.Pt();
740 cosTheta = TMath::Abs(momentum.CosTheta());
741 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
742 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
743 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
744
745 entry = static_cast<Jet *>(branch->NewEntry());
746
747 entry->Eta = eta;
748 entry->Phi = momentum.Phi();
749 entry->PT = pt;
750
751 entry->T = position.T() * 1.0E-3 / c_light;
752
753 entry->Mass = momentum.M();
754
755 entry->Area = candidate->Area;
756
757 entry->DeltaEta = candidate->DeltaEta;
758 entry->DeltaPhi = candidate->DeltaPhi;
759
760 entry->Flavor = candidate->Flavor;
761 entry->FlavorAlgo = candidate->FlavorAlgo;
762 entry->FlavorPhys = candidate->FlavorPhys;
763
764 entry->BTag = candidate->BTag;
765
766 entry->BTagAlgo = candidate->BTagAlgo;
767 entry->BTagPhys = candidate->BTagPhys;
768
769 entry->TauTag = candidate->TauTag;
770 entry->TauWeight = candidate->TauWeight;
771
772 entry->Charge = candidate->Charge;
773
774 itConstituents.Reset();
775 entry->Constituents.Clear();
776 ecalEnergy = 0.0;
777 hcalEnergy = 0.0;
778 while((constituent = static_cast<Candidate *>(itConstituents.Next())))
779 {
780 entry->Constituents.Add(constituent);
781 ecalEnergy += constituent->Eem;
782 hcalEnergy += constituent->Ehad;
783 }
784
785 entry->EhadOverEem = ecalEnergy > 0.0 ? hcalEnergy / ecalEnergy : 999.9;
786
787 //--- Pile-Up Jet ID variables ----
788
789 entry->NCharged = candidate->NCharged;
790 entry->NNeutrals = candidate->NNeutrals;
791
792 entry->NeutralEnergyFraction = candidate->NeutralEnergyFraction;
793 entry->ChargedEnergyFraction = candidate->ChargedEnergyFraction;
794 entry->Beta = candidate->Beta;
795 entry->BetaStar = candidate->BetaStar;
796 entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
797 entry->PTD = candidate->PTD;
798
799 //--- Sub-structure variables ----
800
801 entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
802 entry->NSubJetsPruned = candidate->NSubJetsPruned;
803 entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
804
805 entry->SoftDroppedJet = candidate->SoftDroppedJet;
806 entry->SoftDroppedSubJet1 = candidate->SoftDroppedSubJet1;
807 entry->SoftDroppedSubJet2 = candidate->SoftDroppedSubJet2;
808
809 for(i = 0; i < 5; i++)
810 {
811 entry->FracPt[i] = candidate->FracPt[i];
812 entry->Tau[i] = candidate->Tau[i];
813 entry->TrimmedP4[i] = candidate->TrimmedP4[i];
814 entry->PrunedP4[i] = candidate->PrunedP4[i];
815 entry->SoftDroppedP4[i] = candidate->SoftDroppedP4[i];
816 }
817
818 //--- exclusive clustering variables ---
819 entry->ExclYmerge23 = candidate->ExclYmerge23;
820 entry->ExclYmerge34 = candidate->ExclYmerge34;
821 entry->ExclYmerge45 = candidate->ExclYmerge45;
822 entry->ExclYmerge56 = candidate->ExclYmerge56;
823
824 FillParticles(candidate, &entry->Particles);
825 }
826}
827
828//------------------------------------------------------------------------------
829
830void TreeWriter::ProcessMissingET(ExRootTreeBranch *branch, TObjArray *array)
831{
832 Candidate *candidate = 0;
833 MissingET *entry = 0;
834
835 // get the first entry
836 if((candidate = static_cast<Candidate *>(array->At(0))))
837 {
838 const TLorentzVector &momentum = candidate->Momentum;
839
840 entry = static_cast<MissingET *>(branch->NewEntry());
841
842 entry->Eta = (-momentum).Eta();
843 entry->Phi = (-momentum).Phi();
844 entry->MET = momentum.Pt();
845 }
846}
847
848//------------------------------------------------------------------------------
849
850void TreeWriter::ProcessScalarHT(ExRootTreeBranch *branch, TObjArray *array)
851{
852 Candidate *candidate = 0;
853 ScalarHT *entry = 0;
854
855 // get the first entry
856 if((candidate = static_cast<Candidate *>(array->At(0))))
857 {
858 const TLorentzVector &momentum = candidate->Momentum;
859
860 entry = static_cast<ScalarHT *>(branch->NewEntry());
861
862 entry->HT = momentum.Pt();
863 }
864}
865
866//------------------------------------------------------------------------------
867
868void TreeWriter::ProcessRho(ExRootTreeBranch *branch, TObjArray *array)
869{
870 TIter iterator(array);
871 Candidate *candidate = 0;
872 Rho *entry = 0;
873
874 // loop over all rho
875 iterator.Reset();
876 while((candidate = static_cast<Candidate *>(iterator.Next())))
877 {
878 const TLorentzVector &momentum = candidate->Momentum;
879
880 entry = static_cast<Rho *>(branch->NewEntry());
881
882 entry->Rho = momentum.E();
883 entry->Edges[0] = candidate->Edges[0];
884 entry->Edges[1] = candidate->Edges[1];
885 }
886}
887
888//------------------------------------------------------------------------------
889
890void TreeWriter::ProcessWeight(ExRootTreeBranch *branch, TObjArray *array)
891{
892 Candidate *candidate = 0;
893 Weight *entry = 0;
894
895 // get the first entry
896 if((candidate = static_cast<Candidate *>(array->At(0))))
897 {
898 const TLorentzVector &momentum = candidate->Momentum;
899
900 entry = static_cast<Weight *>(branch->NewEntry());
901
902 entry->Weight = momentum.E();
903 }
904}
905
906//------------------------------------------------------------------------------
907
908void TreeWriter::ProcessHectorHit(ExRootTreeBranch *branch, TObjArray *array)
909{
910 TIter iterator(array);
911 Candidate *candidate = 0;
912 HectorHit *entry = 0;
913
914 // loop over all roman pot hits
915 iterator.Reset();
916 while((candidate = static_cast<Candidate *>(iterator.Next())))
917 {
918 const TLorentzVector &position = candidate->Position;
919 const TLorentzVector &momentum = candidate->Momentum;
920
921 entry = static_cast<HectorHit *>(branch->NewEntry());
922
923 entry->E = momentum.E();
924
925 entry->Tx = momentum.Px();
926 entry->Ty = momentum.Py();
927
928 entry->T = position.T();
929
930 entry->X = position.X();
931 entry->Y = position.Y();
932 entry->S = position.Z();
933
934 entry->Particle = candidate->GetCandidates()->At(0);
935 }
936}
937
938//------------------------------------------------------------------------------
939
940void TreeWriter::Process()
941{
942 TBranchMap::iterator itBranchMap;
943 ExRootTreeBranch *branch;
944 TProcessMethod method;
945 TObjArray *array;
946
947 for(itBranchMap = fBranchMap.begin(); itBranchMap != fBranchMap.end(); ++itBranchMap)
948 {
949 branch = itBranchMap->first;
950 method = itBranchMap->second.first;
951 array = itBranchMap->second.second;
952
953 (this->*method)(branch, array);
954 }
955}
956
957//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.