Fork me on GitHub

source: git/modules/TreeWriter.cc@ 56fb0be

Last change on this file since 56fb0be was b750b0a, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Fix CovarianceMatrix scaling error. Vertexing improvements.

  • Property mode set to 100644
File size: 29.2 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 param = GetParam("Info");
127 TString infoName;
128 Double_t infoValue;
129
130 size = param.GetSize();
131 for(i = 0; i < size / 2; ++i)
132 {
133 infoName = param[i * 2].GetString();
134 infoValue = param[i * 2 + 1].GetDouble();
135
136 AddInfo(infoName, infoValue);
137 }
138}
139
140//------------------------------------------------------------------------------
141
142void TreeWriter::Finish()
143{
144}
145
146//------------------------------------------------------------------------------
147
148void TreeWriter::FillParticles(Candidate *candidate, TRefArray *array)
149{
150 TIter it1(candidate->GetCandidates());
151 it1.Reset();
152 array->Clear();
153
154 while((candidate = static_cast<Candidate *>(it1.Next())))
155 {
156 TIter it2(candidate->GetCandidates());
157
158 // particle
159 if(candidate->GetCandidates()->GetEntriesFast() == 0)
160 {
161 array->Add(candidate);
162 continue;
163 }
164
165 // track
166 candidate = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
167 if(candidate->GetCandidates()->GetEntriesFast() == 0)
168 {
169 array->Add(candidate);
170 continue;
171 }
172
173 // tower
174 it2.Reset();
175 while((candidate = static_cast<Candidate *>(it2.Next())))
176 {
177 array->Add(candidate->GetCandidates()->At(0));
178 }
179 }
180}
181
182//------------------------------------------------------------------------------
183
184void TreeWriter::ProcessParticles(ExRootTreeBranch *branch, TObjArray *array)
185{
186 TIter iterator(array);
187 Candidate *candidate = 0;
188 GenParticle *entry = 0;
189 Double_t pt, signPz, cosTheta, eta, rapidity;
190
191 const Double_t c_light = 2.99792458E8;
192
193 // loop over all particles
194 iterator.Reset();
195 while((candidate = static_cast<Candidate *>(iterator.Next())))
196 {
197 const TLorentzVector &momentum = candidate->Momentum;
198 const TLorentzVector &position = candidate->Position;
199
200 entry = static_cast<GenParticle *>(branch->NewEntry());
201
202 entry->SetBit(kIsReferenced);
203 entry->SetUniqueID(candidate->GetUniqueID());
204
205 pt = momentum.Pt();
206 cosTheta = TMath::Abs(momentum.CosTheta());
207 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
208 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
209 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
210
211 entry->PID = candidate->PID;
212
213 entry->Status = candidate->Status;
214 entry->IsPU = candidate->IsPU;
215
216 entry->M1 = candidate->M1;
217 entry->M2 = candidate->M2;
218
219 entry->D1 = candidate->D1;
220 entry->D2 = candidate->D2;
221
222 entry->Charge = candidate->Charge;
223 entry->Mass = candidate->Mass;
224
225 entry->E = momentum.E();
226 entry->Px = momentum.Px();
227 entry->Py = momentum.Py();
228 entry->Pz = momentum.Pz();
229
230 entry->Eta = eta;
231 entry->Phi = momentum.Phi();
232 entry->PT = pt;
233
234 entry->Rapidity = rapidity;
235
236 entry->X = position.X();
237 entry->Y = position.Y();
238 entry->Z = position.Z();
239 entry->T = position.T() * 1.0E-3 / c_light;
240 }
241}
242
243//------------------------------------------------------------------------------
244
245void TreeWriter::ProcessVertices(ExRootTreeBranch *branch, TObjArray *array)
246{
247 TIter iterator(array);
248 Candidate *candidate = 0, *constituent = 0;
249 Vertex *entry = 0;
250
251 const Double_t c_light = 2.99792458E8;
252
253 Double_t x, y, z, t, xError, yError, zError, tError, sigma, sumPT2, btvSumPT2, genDeltaZ, genSumPT2;
254 UInt_t index, ndf;
255
256 CompBase *compare = Candidate::fgCompare;
257 Candidate::fgCompare = CompSumPT2<Candidate>::Instance();
258 array->Sort();
259 Candidate::fgCompare = compare;
260
261 // loop over all vertices
262 iterator.Reset();
263 while((candidate = static_cast<Candidate *>(iterator.Next())))
264 {
265
266 index = candidate->ClusterIndex;
267 ndf = candidate->ClusterNDF;
268 sigma = candidate->ClusterSigma;
269 sumPT2 = candidate->SumPT2;
270 btvSumPT2 = candidate->BTVSumPT2;
271 genDeltaZ = candidate->GenDeltaZ;
272 genSumPT2 = candidate->GenSumPT2;
273
274 x = candidate->Position.X();
275 y = candidate->Position.Y();
276 z = candidate->Position.Z();
277 t = candidate->Position.T() * 1.0E-3 / c_light;
278
279 xError = candidate->PositionError.X();
280 yError = candidate->PositionError.Y();
281 zError = candidate->PositionError.Z();
282 tError = candidate->PositionError.T() * 1.0E-3 / c_light;
283
284 entry = static_cast<Vertex *>(branch->NewEntry());
285
286 entry->Index = index;
287 entry->NDF = ndf;
288 entry->Sigma = sigma;
289 entry->SumPT2 = sumPT2;
290 entry->BTVSumPT2 = btvSumPT2;
291 entry->GenDeltaZ = genDeltaZ;
292 entry->GenSumPT2 = genSumPT2;
293
294 entry->X = x;
295 entry->Y = y;
296 entry->Z = z;
297 entry->T = t;
298
299 entry->ErrorX = xError;
300 entry->ErrorY = yError;
301 entry->ErrorZ = zError;
302 entry->ErrorT = tError;
303
304 TIter itConstituents(candidate->GetCandidates());
305 itConstituents.Reset();
306 entry->Constituents.Clear();
307 while((constituent = static_cast<Candidate *>(itConstituents.Next())))
308 {
309 entry->Constituents.Add(constituent);
310 }
311 }
312}
313
314//------------------------------------------------------------------------------
315
316void TreeWriter::ProcessTracks(ExRootTreeBranch *branch, TObjArray *array)
317{
318 TIter iterator(array);
319 Candidate *candidate = 0;
320 Candidate *particle = 0;
321 Track *entry = 0;
322 Double_t pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
323 const Double_t c_light = 2.99792458E8;
324
325 // loop over all tracks
326 iterator.Reset();
327 while((candidate = static_cast<Candidate *>(iterator.Next())))
328 {
329 const TLorentzVector &position = candidate->Position;
330
331 cosTheta = TMath::Abs(position.CosTheta());
332 signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
333 eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
334 rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
335
336 entry = static_cast<Track *>(branch->NewEntry());
337
338 entry->SetBit(kIsReferenced);
339 entry->SetUniqueID(candidate->GetUniqueID());
340
341 entry->PID = candidate->PID;
342
343 entry->Charge = candidate->Charge;
344
345 entry->EtaOuter = eta;
346 entry->PhiOuter = position.Phi();
347
348 entry->XOuter = position.X();
349 entry->YOuter = position.Y();
350 entry->ZOuter = position.Z();
351 entry->TOuter = position.T() * 1.0E-3 / c_light;
352
353 entry->L = candidate->L;
354
355 entry->D0 = candidate->D0;
356 entry->DZ = candidate->DZ;
357 entry->Nclusters = candidate->Nclusters;
358 entry->dNdx = candidate->dNdx;
359
360 entry->ErrorP = candidate->ErrorP;
361 entry->ErrorPT = candidate->ErrorPT;
362
363 // diagonal covariance matrix terms
364 entry->ErrorD0 = candidate->ErrorD0;
365 entry->ErrorC = candidate->ErrorC;
366 entry->ErrorPhi = candidate->ErrorPhi;
367 entry->ErrorDZ = candidate->ErrorDZ;
368 entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
369
370 // add some offdiagonal covariance matrix elements
371 entry->ErrorD0Phi = candidate->TrackCovariance(0,1)*1.e3;
372 entry->ErrorD0C = candidate->TrackCovariance(0,2);
373 entry->ErrorD0DZ = candidate->TrackCovariance(0,3)*1.e6;
374 entry->ErrorD0CtgTheta = candidate->TrackCovariance(0,4)*1.e3;
375 entry->ErrorPhiC = candidate->TrackCovariance(1,2)*1.e-3;
376 entry->ErrorPhiDZ = candidate->TrackCovariance(1,3)*1.e3;
377 entry->ErrorPhiCtgTheta = candidate->TrackCovariance(1,4);
378 entry->ErrorCDZ = candidate->TrackCovariance(2,3);
379 entry->ErrorCCtgTheta = candidate->TrackCovariance(2,4)*1.e-3;
380 entry->ErrorDZCtgTheta = candidate->TrackCovariance(3,4)*1.e3;
381
382 entry->Xd = candidate->Xd;
383 entry->Yd = candidate->Yd;
384 entry->Zd = candidate->Zd;
385
386 const TLorentzVector &momentum = candidate->Momentum;
387
388 pt = momentum.Pt();
389 p = momentum.P();
390 phi = momentum.Phi();
391 m = momentum.M();
392 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
393
394 cosTheta = TMath::Abs(momentum.CosTheta());
395 signz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
396 eta = (cosTheta == 1.0 ? signz * 999.9 : momentum.Eta());
397 rapidity = (cosTheta == 1.0 ? signz * 999.9 : momentum.Rapidity());
398
399 entry->P = p;
400 entry->PT = pt;
401 entry->Eta = eta;
402 entry->Phi = phi;
403 entry->CtgTheta = ctgTheta;
404 entry->C = candidate->C;
405 entry->Mass = m;
406
407 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
408 //const TLorentzVector &initialPosition = particle->Position;
409 const TLorentzVector &initialPosition = candidate->InitialPosition;
410
411 entry->X = initialPosition.X();
412 entry->Y = initialPosition.Y();
413 entry->Z = initialPosition.Z();
414 entry->T = initialPosition.T() * 1.0E-3 / c_light;
415 entry->ErrorT =candidate-> ErrorT * 1.0E-3 / c_light;
416
417 entry->Particle = particle;
418
419 entry->VertexIndex = candidate->ClusterIndex;
420 }
421}
422
423//------------------------------------------------------------------------------
424
425void TreeWriter::ProcessTowers(ExRootTreeBranch *branch, TObjArray *array)
426{
427 TIter iterator(array);
428 Candidate *candidate = 0;
429 Tower *entry = 0;
430 Double_t pt, signPz, cosTheta, eta, rapidity;
431 const Double_t c_light = 2.99792458E8;
432
433 // loop over all towers
434 iterator.Reset();
435 while((candidate = static_cast<Candidate *>(iterator.Next())))
436 {
437 const TLorentzVector &momentum = candidate->Momentum;
438 const TLorentzVector &position = candidate->Position;
439
440 pt = momentum.Pt();
441 cosTheta = TMath::Abs(momentum.CosTheta());
442 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
443 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
444 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
445
446 entry = static_cast<Tower *>(branch->NewEntry());
447
448 entry->SetBit(kIsReferenced);
449 entry->SetUniqueID(candidate->GetUniqueID());
450
451 entry->Eta = eta;
452 entry->Phi = momentum.Phi();
453 entry->ET = pt;
454 entry->E = momentum.E();
455 entry->Eem = candidate->Eem;
456 entry->Ehad = candidate->Ehad;
457 entry->Etrk = candidate->Etrk;
458 entry->Edges[0] = candidate->Edges[0];
459 entry->Edges[1] = candidate->Edges[1];
460 entry->Edges[2] = candidate->Edges[2];
461 entry->Edges[3] = candidate->Edges[3];
462
463 entry->T = position.T() * 1.0E-3 / c_light;
464 entry->NTimeHits = candidate->NTimeHits;
465
466 FillParticles(candidate, &entry->Particles);
467 }
468}
469
470//------------------------------------------------------------------------------
471
472void TreeWriter::ProcessParticleFlowCandidates(ExRootTreeBranch *branch, TObjArray *array)
473{
474
475 TIter iterator(array);
476 Candidate *candidate = 0;
477 Candidate *particle = 0;
478 ParticleFlowCandidate *entry = 0;
479 Double_t e, pt, signz, cosTheta, eta, rapidity, p, ctgTheta, phi, m;
480 const Double_t c_light = 2.99792458E8;
481
482 // loop over all tracks
483 iterator.Reset();
484 while((candidate = static_cast<Candidate *>(iterator.Next())))
485 {
486 const TLorentzVector &position = candidate->Position;
487
488 cosTheta = TMath::Abs(position.CosTheta());
489 signz = (position.Pz() >= 0.0) ? 1.0 : -1.0;
490 eta = (cosTheta == 1.0 ? signz * 999.9 : position.Eta());
491 rapidity = (cosTheta == 1.0 ? signz * 999.9 : position.Rapidity());
492
493 entry = static_cast<ParticleFlowCandidate *>(branch->NewEntry());
494
495 entry->SetBit(kIsReferenced);
496 entry->SetUniqueID(candidate->GetUniqueID());
497
498 entry->PID = candidate->PID;
499
500 entry->Charge = candidate->Charge;
501
502 entry->EtaOuter = eta;
503 entry->PhiOuter = position.Phi();
504
505 entry->XOuter = position.X();
506 entry->YOuter = position.Y();
507 entry->ZOuter = position.Z();
508 entry->TOuter = position.T() * 1.0E-3 / c_light;
509
510 entry->L = candidate->L;
511
512 entry->D0 = candidate->D0;
513 entry->DZ = candidate->DZ;
514 entry->Nclusters = candidate->Nclusters;
515 entry->dNdx = candidate->dNdx;
516
517 entry->ErrorP = candidate->ErrorP;
518 entry->ErrorPT = candidate->ErrorPT;
519 entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
520
521
522 // diagonal covariance matrix terms
523
524 entry->ErrorD0 = candidate->ErrorD0;
525 entry->ErrorC = candidate->ErrorC;
526 entry->ErrorPhi = candidate->ErrorPhi;
527 entry->ErrorDZ = candidate->ErrorDZ;
528 entry->ErrorCtgTheta = candidate->ErrorCtgTheta;
529
530 // add some offdiagonal covariance matrix elements
531 entry->ErrorD0Phi = candidate->TrackCovariance(0,1);
532 entry->ErrorD0C = candidate->TrackCovariance(0,2);
533 entry->ErrorD0DZ = candidate->TrackCovariance(0,3);
534 entry->ErrorD0CtgTheta = candidate->TrackCovariance(0,4);
535 entry->ErrorPhiC = candidate->TrackCovariance(1,2);
536 entry->ErrorPhiDZ = candidate->TrackCovariance(1,3);
537 entry->ErrorPhiCtgTheta = candidate->TrackCovariance(1,4);
538 entry->ErrorCDZ = candidate->TrackCovariance(2,3);
539 entry->ErrorCCtgTheta = candidate->TrackCovariance(2,4);
540 entry->ErrorDZCtgTheta = candidate->TrackCovariance(3,4);
541
542 entry->Xd = candidate->Xd;
543 entry->Yd = candidate->Yd;
544 entry->Zd = candidate->Zd;
545
546 const TLorentzVector &momentum = candidate->Momentum;
547
548 e = momentum.E();
549 pt = momentum.Pt();
550 p = momentum.P();
551 phi = momentum.Phi();
552 m = momentum.M();
553 ctgTheta = (TMath::Tan(momentum.Theta()) != 0) ? 1 / TMath::Tan(momentum.Theta()) : 1e10;
554
555 entry->E = e;
556 entry->P = p;
557 entry->PT = pt;
558 entry->Eta = eta;
559 entry->Phi = phi;
560 entry->CtgTheta = ctgTheta;
561 entry->C = candidate->C;
562 entry->Mass = m;
563
564 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
565 //const TLorentzVector &initialPosition = particle->Position;
566 const TLorentzVector &initialPosition = candidate->InitialPosition;
567
568 entry->X = initialPosition.X();
569 entry->Y = initialPosition.Y();
570 entry->Z = initialPosition.Z();
571 entry->T = initialPosition.T() * 1.0E-3 / c_light;
572 entry->ErrorT = candidate-> ErrorT * 1.0E-3 / c_light;
573
574 entry->VertexIndex = candidate->ClusterIndex;
575
576 entry->Eem = candidate->Eem;
577 entry->Ehad = candidate->Ehad;
578 entry->Etrk = candidate->Etrk;
579 entry->Edges[0] = candidate->Edges[0];
580 entry->Edges[1] = candidate->Edges[1];
581 entry->Edges[2] = candidate->Edges[2];
582 entry->Edges[3] = candidate->Edges[3];
583
584 //entry->T = position.T() * 1.0E-3 / c_light;
585 entry->NTimeHits = candidate->NTimeHits;
586
587 FillParticles(candidate, &entry->Particles);
588
589 }
590}
591
592//------------------------------------------------------------------------------
593
594void TreeWriter::ProcessPhotons(ExRootTreeBranch *branch, TObjArray *array)
595{
596 TIter iterator(array);
597 Candidate *candidate = 0;
598 Photon *entry = 0;
599 Double_t pt, signPz, cosTheta, eta, rapidity;
600 const Double_t c_light = 2.99792458E8;
601
602 array->Sort();
603
604 // loop over all photons
605 iterator.Reset();
606 while((candidate = static_cast<Candidate *>(iterator.Next())))
607 {
608 TIter it1(candidate->GetCandidates());
609 const TLorentzVector &momentum = candidate->Momentum;
610 const TLorentzVector &position = candidate->Position;
611
612 pt = momentum.Pt();
613 cosTheta = TMath::Abs(momentum.CosTheta());
614 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
615 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
616 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
617
618 entry = static_cast<Photon *>(branch->NewEntry());
619
620 entry->Eta = eta;
621 entry->Phi = momentum.Phi();
622 entry->PT = pt;
623 entry->E = momentum.E();
624 entry->T = position.T() * 1.0E-3 / c_light;
625
626 // Isolation variables
627
628 entry->IsolationVar = candidate->IsolationVar;
629 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
630 entry->SumPtCharged = candidate->SumPtCharged;
631 entry->SumPtNeutral = candidate->SumPtNeutral;
632 entry->SumPtChargedPU = candidate->SumPtChargedPU;
633 entry->SumPt = candidate->SumPt;
634
635 entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad / candidate->Eem : 999.9;
636
637 // 1: prompt -- 2: non prompt -- 3: fake
638 entry->Status = candidate->Status;
639
640 FillParticles(candidate, &entry->Particles);
641 }
642}
643
644//------------------------------------------------------------------------------
645
646void TreeWriter::ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array)
647{
648 TIter iterator(array);
649 Candidate *candidate = 0;
650 Electron *entry = 0;
651 Double_t pt, signPz, cosTheta, eta, rapidity;
652 const Double_t c_light = 2.99792458E8;
653
654 array->Sort();
655
656 // loop over all electrons
657 iterator.Reset();
658 while((candidate = static_cast<Candidate *>(iterator.Next())))
659 {
660 const TLorentzVector &momentum = candidate->Momentum;
661 const TLorentzVector &position = candidate->Position;
662
663 pt = momentum.Pt();
664 cosTheta = TMath::Abs(momentum.CosTheta());
665 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
666 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
667 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
668
669 entry = static_cast<Electron *>(branch->NewEntry());
670
671 entry->Eta = eta;
672 entry->Phi = momentum.Phi();
673 entry->PT = pt;
674
675 entry->T = position.T() * 1.0E-3 / c_light;
676
677 // displacement
678 entry->D0 = candidate->D0;
679 entry->ErrorD0 = candidate->ErrorD0;
680 entry->DZ = candidate->DZ;
681 entry->ErrorDZ = candidate->ErrorDZ;
682
683 // Isolation variables
684 entry->IsolationVar = candidate->IsolationVar;
685 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
686 entry->SumPtCharged = candidate->SumPtCharged;
687 entry->SumPtNeutral = candidate->SumPtNeutral;
688 entry->SumPtChargedPU = candidate->SumPtChargedPU;
689 entry->SumPt = candidate->SumPt;
690
691 entry->Charge = candidate->Charge;
692
693 entry->EhadOverEem = 0.0;
694
695 entry->Particle = candidate->GetCandidates()->At(0);
696 }
697}
698
699//------------------------------------------------------------------------------
700
701void TreeWriter::ProcessMuons(ExRootTreeBranch *branch, TObjArray *array)
702{
703 TIter iterator(array);
704 Candidate *candidate = 0;
705 Muon *entry = 0;
706 Double_t pt, signPz, cosTheta, eta, rapidity;
707
708 const Double_t c_light = 2.99792458E8;
709
710 array->Sort();
711
712 // loop over all muons
713 iterator.Reset();
714 while((candidate = static_cast<Candidate *>(iterator.Next())))
715 {
716 const TLorentzVector &momentum = candidate->Momentum;
717 const TLorentzVector &position = candidate->Position;
718
719 pt = momentum.Pt();
720 cosTheta = TMath::Abs(momentum.CosTheta());
721 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
722 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
723 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
724
725 entry = static_cast<Muon *>(branch->NewEntry());
726
727 entry->SetBit(kIsReferenced);
728 entry->SetUniqueID(candidate->GetUniqueID());
729
730 entry->Eta = eta;
731 entry->Phi = momentum.Phi();
732 entry->PT = pt;
733
734 entry->T = position.T() * 1.0E-3 / c_light;
735
736 // displacement
737 entry->D0 = candidate->D0;
738 entry->ErrorD0 = candidate->ErrorD0;
739 entry->DZ = candidate->DZ;
740 entry->ErrorDZ = candidate->ErrorDZ;
741
742 // Isolation variables
743
744 entry->IsolationVar = candidate->IsolationVar;
745 entry->IsolationVarRhoCorr = candidate->IsolationVarRhoCorr;
746 entry->SumPtCharged = candidate->SumPtCharged;
747 entry->SumPtNeutral = candidate->SumPtNeutral;
748 entry->SumPtChargedPU = candidate->SumPtChargedPU;
749 entry->SumPt = candidate->SumPt;
750
751 entry->Charge = candidate->Charge;
752
753 entry->Particle = candidate->GetCandidates()->At(0);
754 }
755}
756
757//------------------------------------------------------------------------------
758
759void TreeWriter::ProcessJets(ExRootTreeBranch *branch, TObjArray *array)
760{
761 TIter iterator(array);
762 Candidate *candidate = 0, *constituent = 0;
763 Jet *entry = 0;
764 Double_t pt, signPz, cosTheta, eta, rapidity;
765 Double_t ecalEnergy, hcalEnergy;
766 const Double_t c_light = 2.99792458E8;
767 Int_t i;
768
769 array->Sort();
770
771 // loop over all jets
772 iterator.Reset();
773 while((candidate = static_cast<Candidate *>(iterator.Next())))
774 {
775 TIter itConstituents(candidate->GetCandidates());
776
777 const TLorentzVector &momentum = candidate->Momentum;
778 const TLorentzVector &position = candidate->Position;
779
780 pt = momentum.Pt();
781 cosTheta = TMath::Abs(momentum.CosTheta());
782 signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
783 eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
784 rapidity = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Rapidity());
785
786 entry = static_cast<Jet *>(branch->NewEntry());
787
788 entry->Eta = eta;
789 entry->Phi = momentum.Phi();
790 entry->PT = pt;
791
792 entry->T = position.T() * 1.0E-3 / c_light;
793
794 entry->Mass = momentum.M();
795
796 entry->Area = candidate->Area;
797
798 entry->DeltaEta = candidate->DeltaEta;
799 entry->DeltaPhi = candidate->DeltaPhi;
800
801 entry->Flavor = candidate->Flavor;
802 entry->FlavorAlgo = candidate->FlavorAlgo;
803 entry->FlavorPhys = candidate->FlavorPhys;
804
805 entry->BTag = candidate->BTag;
806
807 entry->BTagAlgo = candidate->BTagAlgo;
808 entry->BTagPhys = candidate->BTagPhys;
809
810 entry->TauTag = candidate->TauTag;
811 entry->TauWeight = candidate->TauWeight;
812
813 entry->Charge = candidate->Charge;
814
815 itConstituents.Reset();
816 entry->Constituents.Clear();
817 ecalEnergy = 0.0;
818 hcalEnergy = 0.0;
819 while((constituent = static_cast<Candidate *>(itConstituents.Next())))
820 {
821 entry->Constituents.Add(constituent);
822 ecalEnergy += constituent->Eem;
823 hcalEnergy += constituent->Ehad;
824 }
825
826 entry->EhadOverEem = ecalEnergy > 0.0 ? hcalEnergy / ecalEnergy : 999.9;
827
828 //--- Pile-Up Jet ID variables ----
829
830 entry->NCharged = candidate->NCharged;
831 entry->NNeutrals = candidate->NNeutrals;
832
833 entry->NeutralEnergyFraction = candidate->NeutralEnergyFraction;
834 entry->ChargedEnergyFraction = candidate->ChargedEnergyFraction;
835 entry->Beta = candidate->Beta;
836 entry->BetaStar = candidate->BetaStar;
837 entry->MeanSqDeltaR = candidate->MeanSqDeltaR;
838 entry->PTD = candidate->PTD;
839
840 //--- Sub-structure variables ----
841
842 entry->NSubJetsTrimmed = candidate->NSubJetsTrimmed;
843 entry->NSubJetsPruned = candidate->NSubJetsPruned;
844 entry->NSubJetsSoftDropped = candidate->NSubJetsSoftDropped;
845
846 entry->SoftDroppedJet = candidate->SoftDroppedJet;
847 entry->SoftDroppedSubJet1 = candidate->SoftDroppedSubJet1;
848 entry->SoftDroppedSubJet2 = candidate->SoftDroppedSubJet2;
849
850 for(i = 0; i < 5; i++)
851 {
852 entry->FracPt[i] = candidate->FracPt[i];
853 entry->Tau[i] = candidate->Tau[i];
854 entry->TrimmedP4[i] = candidate->TrimmedP4[i];
855 entry->PrunedP4[i] = candidate->PrunedP4[i];
856 entry->SoftDroppedP4[i] = candidate->SoftDroppedP4[i];
857 }
858
859 //--- exclusive clustering variables ---
860 entry->ExclYmerge23 = candidate->ExclYmerge23;
861 entry->ExclYmerge34 = candidate->ExclYmerge34;
862 entry->ExclYmerge45 = candidate->ExclYmerge45;
863 entry->ExclYmerge56 = candidate->ExclYmerge56;
864
865 FillParticles(candidate, &entry->Particles);
866 }
867}
868
869//------------------------------------------------------------------------------
870
871void TreeWriter::ProcessMissingET(ExRootTreeBranch *branch, TObjArray *array)
872{
873 Candidate *candidate = 0;
874 MissingET *entry = 0;
875
876 // get the first entry
877 if((candidate = static_cast<Candidate *>(array->At(0))))
878 {
879 const TLorentzVector &momentum = candidate->Momentum;
880
881 entry = static_cast<MissingET *>(branch->NewEntry());
882
883 entry->Eta = (-momentum).Eta();
884 entry->Phi = (-momentum).Phi();
885 entry->MET = momentum.Pt();
886 }
887}
888
889//------------------------------------------------------------------------------
890
891void TreeWriter::ProcessScalarHT(ExRootTreeBranch *branch, TObjArray *array)
892{
893 Candidate *candidate = 0;
894 ScalarHT *entry = 0;
895
896 // get the first entry
897 if((candidate = static_cast<Candidate *>(array->At(0))))
898 {
899 const TLorentzVector &momentum = candidate->Momentum;
900
901 entry = static_cast<ScalarHT *>(branch->NewEntry());
902
903 entry->HT = momentum.Pt();
904 }
905}
906
907//------------------------------------------------------------------------------
908
909void TreeWriter::ProcessRho(ExRootTreeBranch *branch, TObjArray *array)
910{
911 TIter iterator(array);
912 Candidate *candidate = 0;
913 Rho *entry = 0;
914
915 // loop over all rho
916 iterator.Reset();
917 while((candidate = static_cast<Candidate *>(iterator.Next())))
918 {
919 const TLorentzVector &momentum = candidate->Momentum;
920
921 entry = static_cast<Rho *>(branch->NewEntry());
922
923 entry->Rho = momentum.E();
924 entry->Edges[0] = candidate->Edges[0];
925 entry->Edges[1] = candidate->Edges[1];
926 }
927}
928
929//------------------------------------------------------------------------------
930
931void TreeWriter::ProcessWeight(ExRootTreeBranch *branch, TObjArray *array)
932{
933 Candidate *candidate = 0;
934 Weight *entry = 0;
935
936 // get the first entry
937 if((candidate = static_cast<Candidate *>(array->At(0))))
938 {
939 const TLorentzVector &momentum = candidate->Momentum;
940
941 entry = static_cast<Weight *>(branch->NewEntry());
942
943 entry->Weight = momentum.E();
944 }
945}
946
947//------------------------------------------------------------------------------
948
949void TreeWriter::ProcessHectorHit(ExRootTreeBranch *branch, TObjArray *array)
950{
951 TIter iterator(array);
952 Candidate *candidate = 0;
953 HectorHit *entry = 0;
954
955 // loop over all roman pot hits
956 iterator.Reset();
957 while((candidate = static_cast<Candidate *>(iterator.Next())))
958 {
959 const TLorentzVector &position = candidate->Position;
960 const TLorentzVector &momentum = candidate->Momentum;
961
962 entry = static_cast<HectorHit *>(branch->NewEntry());
963
964 entry->E = momentum.E();
965
966 entry->Tx = momentum.Px();
967 entry->Ty = momentum.Py();
968
969 entry->T = position.T();
970
971 entry->X = position.X();
972 entry->Y = position.Y();
973 entry->S = position.Z();
974
975 entry->Particle = candidate->GetCandidates()->At(0);
976 }
977}
978
979//------------------------------------------------------------------------------
980
981void TreeWriter::Process()
982{
983 TBranchMap::iterator itBranchMap;
984 ExRootTreeBranch *branch;
985 TProcessMethod method;
986 TObjArray *array;
987
988 for(itBranchMap = fBranchMap.begin(); itBranchMap != fBranchMap.end(); ++itBranchMap)
989 {
990 branch = itBranchMap->first;
991 method = itBranchMap->second.first;
992 array = itBranchMap->second.second;
993
994 (this->*method)(branch, array);
995 }
996}
997
998//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.