Fork me on GitHub

source: git/modules/Hector.cc@ 239e1d0

ImprovedOutputFile Timing dual_readout llp 3.3.3pre11
Last change on this file since 239e1d0 was cab38f6, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

remove svn tags and fix formatting

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