Fork me on GitHub

source: git/modules/Hector.cc@ caba091

ImprovedOutputFile Timing dual_readout llp
Last change on this file since caba091 was a8782e8, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago

replace examples by cards

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