Fork me on GitHub

source: git/modules/Hector.cc@ ededa33

ImprovedOutputFile Timing
Last change on this file since ededa33 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 5.1 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 Hector
20 *
21 * Propagates candidates using Hector library.
22 *
23 * \author P. Demin - UCL, Louvain-la-Neuve
24 *
25 */
26
27#include "modules/Hector.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
37#include "TDatabasePDG.h"
38#include "TFormula.h"
39#include "TLorentzVector.h"
40#include "TMath.h"
41#include "TObjArray.h"
42#include "TRandom3.h"
43#include "TString.h"
44
45#include <algorithm>
46#include <iostream>
47#include <sstream>
48#include <stdexcept>
49
50#include "Hector/H_BeamLine.h"
51#include "Hector/H_BeamParticle.h"
52#include "Hector/H_RecRPObject.h"
53
54using namespace std;
55
56//------------------------------------------------------------------------------
57
58Hector::Hector() :
59 fBeamLine(0), fItInputArray(0)
60{
61}
62
63//------------------------------------------------------------------------------
64
65Hector::~Hector()
66{
67}
68
69//------------------------------------------------------------------------------
70
71void Hector::Init()
72{
73 // read Hector parameters
74
75 fDirection = GetInt("Direction", 1);
76 fBeamLineLength = GetDouble("BeamLineLength", 430.0);
77 fDistance = GetDouble("Distance", 420.0);
78 fOffsetX = GetDouble("OffsetX", 0.0);
79 fOffsetS = GetDouble("OffsetS", 120.0);
80 fSigmaE = GetDouble("SigmaE", 0.0);
81 fSigmaX = GetDouble("SigmaX", 0.0);
82 fSigmaY = GetDouble("SigmaY", 0.0);
83 fSigmaT = GetDouble("SigmaT", 0.0);
84 fEtaMin = GetDouble("EtaMin", 5.0);
85
86 fBeamLine = new H_BeamLine(fDirection, fBeamLineLength + 0.1);
87 fBeamLine->fill(GetString("BeamLineFile", "cards/LHCB1IR5_5TeV.tfs"), fDirection, GetString("IPName", "IP5"));
88 fBeamLine->offsetElements(fOffsetS, fOffsetX);
89 fBeamLine->calcMatrix();
90
91 // import input array
92
93 fInputArray = ImportArray(GetString("InputArray", "ParticlePropagator/stableParticles"));
94 fItInputArray = fInputArray->MakeIterator();
95
96 // create output array
97
98 fOutputArray = ExportArray(GetString("OutputArray", "hits"));
99}
100
101//------------------------------------------------------------------------------
102
103void Hector::Finish()
104{
105 if(fItInputArray) delete fItInputArray;
106 if(fBeamLine) delete fBeamLine;
107}
108
109//------------------------------------------------------------------------------
110
111void Hector::Process()
112{
113 Candidate *candidate, *mother;
114 Double_t pz;
115 Double_t x, y, z, tx, ty, theta;
116 Double_t distance, time;
117
118 const Double_t c_light = 2.99792458E8;
119
120 fItInputArray->Reset();
121 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
122 {
123 const TLorentzVector &candidatePosition = candidate->Position;
124 const TLorentzVector &candidateMomentum = candidate->Momentum;
125 pz = candidateMomentum.Pz();
126
127 if(TMath::Abs(candidateMomentum.Eta()) <= fEtaMin || TMath::Sign(pz, Double_t(fDirection)) != pz) continue;
128
129 x = 1.0E3 * candidatePosition.X();
130 y = 1.0E3 * candidatePosition.Y();
131 z = 1.0E-3 * candidatePosition.Z();
132
133 // tx = 1.0E6 * TMath::ATan(candidateMomentum.Px()/pz);
134 // ty = 1.0E6 * TMath::ATan(candidateMomentum.Py()/pz);
135
136 tx = 0.0;
137 ty = 0.0;
138
139 theta = TMath::Hypot(TMath::ATan(candidateMomentum.Px() / pz), TMath::ATan(candidateMomentum.Py() / pz));
140 distance = (fDistance - 1.0E-3 * candidatePosition.Z()) / TMath::Cos(theta);
141 time = gRandom->Gaus((distance + 1.0E-3 * candidatePosition.T()) / c_light, fSigmaT);
142
143 H_BeamParticle particle(candidate->Mass, candidate->Charge);
144 // particle.set4Momentum(candidateMomentum);
145 particle.set4Momentum(candidateMomentum.Px(), candidateMomentum.Py(),
146 candidateMomentum.Pz(), candidateMomentum.E());
147 particle.setPosition(x, y, tx, ty, z);
148
149 particle.smearAng(fSigmaX, fSigmaY, gRandom);
150 particle.smearE(fSigmaE, gRandom);
151
152 particle.computePath(fBeamLine);
153
154 if(particle.stopped(fBeamLine)) continue;
155
156 particle.propagate(fDistance);
157
158 mother = candidate;
159 candidate = static_cast<Candidate *>(candidate->Clone());
160 candidate->Position.SetXYZT(particle.getX(), particle.getY(), particle.getS(), time);
161 candidate->Momentum.SetPxPyPzE(particle.getTX(), particle.getTY(), 0.0, particle.getE());
162 candidate->AddCandidate(mother);
163
164 fOutputArray->Add(candidate);
165 }
166}
167
168//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.