Fork me on GitHub

source: svn/trunk/src/VeryForward.cc@ 376

Last change on this file since 376 was 374, checked in by Xavier Rouby, 16 years ago

added Resolution terms for energy and timing for ZDC
File size: 12.8 KB
RevLine 
[260]1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
[374]6** \------------- arXiv:0903.2225v1 ------------/ **
[260]7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
[53]31
[219]32#include "VeryForward.h"
33#include "H_RomanPot.h"
[53]34#include <iostream>
35#include<cmath>
36
37using namespace std;
38
39
40//------------------------------------------------------------------------------
[374]41VeryForward::VeryForward() :
42 DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
43 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
44 relative_energy(true), // should always be true
45 kickers_on(1) // should always be 1
46 {
47 init(); //Initialisation of Hector
[219]48}
49
[374]50VeryForward::VeryForward(const string& DetDatacard) :
51 DET(new RESOLution())
52 {
[219]53 DET->ReadDataCard(DetDatacard);
[374]54 const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
55 beamline1 = new H_BeamLine(1,d_max);
56 beamline2 = new H_BeamLine(1,d_max);
57 init(); //Initialisation of Hector
58 relative_energy = true; // should always be true
59 kickers_on = 1; // should always be 1
[219]60}
61
[374]62VeryForward::VeryForward(const RESOLution * DetDatacard) :
63 DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
64 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
65 relative_energy(true), // should always be true
66 kickers_on(1) // should always be 1
67 {
68 init(); //Initialisation of Hector
[219]69}
70
[374]71VeryForward::VeryForward(const VeryForward& vf) :
72 DET(new RESOLution(*(vf.DET))), d_max(vf.d_max),
73 beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
74 relative_energy(vf.relative_energy),
75 kickers_on(vf.kickers_on) {
[219]76}
77
78VeryForward& VeryForward::operator=(const VeryForward& vf){
79 if (this==&vf) return *this;
80 DET = new RESOLution(*(vf.DET));
[374]81 d_max = vf.d_max;
[219]82 beamline1 = new H_BeamLine(*(vf.beamline1));
83 beamline2 = new H_BeamLine(*(vf.beamline2));
[374]84 relative_energy =vf.relative_energy;
85 kickers_on = vf.kickers_on;
[219]86 return *this;
87}
88
89
90void VeryForward::init() {
[53]91 //Initialisation of Hector
92 relative_energy = true; // should always be true
93 kickers_on = 1; // should always be 1
[257]94 beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
[252]95 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
[254]96 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*(1E6)); // RP 220m, 2mm, beam 1
97 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*(1E6)); // RP 420m, 4mm, beam 1
[242]98 beamline1->add(rp220_1);
[53]99 beamline1->add(rp420_1);
100
[257]101 beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
[252]102 beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x);
[254]103 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*(1E6));// RP 220m, 2mm, beam 2
104 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*(1E6));// RP 420m, 4mm, beam 2
[53]105 beamline2->add(rp220_2);
106 beamline2->add(rp420_2);
[242]107 // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
108 // do not put explicit delete
[53]109}
110
[242]111
[374]112void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
[53]113{
114 TRootZdcHits *elementZdc;
[374]115 float energy = particle->E;
116 //TLorentzVector genMomentum;
117
[53]118 // Zero degree calorimeter, for forward neutrons and photons
[374]119 if (particle->Status ==1 && ( (particle->PID==pN && energy>DET->ZDC_n_E) ||
120 (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
121 && fabs(particle->Eta) > DET->VFD_min_zdc ) {
[53]122 elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
[374]123
124 // 1) energy smearing
125 float energyS = -1.;
126 if (particle->PID == pGAMMA)
127 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
128 pow(DET->ELG_Czdc*particle->E,2) +
129 pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
130 else // smearing with hadronic resolution
131 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
132 pow(DET->HAD_Czdc*particle->E,2) +
133 pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
134 elementZdc->E = energyS;
135
136
137 // 2) time of flight t is t = T + d/[ cos(theta) v ]
138 float cos_theta = 1; //very good approximation, if eta_zdc >3
139 if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
140 double tx = atan(particle->Px/particle->Pz);
141 double ty = atan(particle->Py/particle->Pz);
142 double theta = sqrt( pow(tx,2) + pow(ty,2) );
143 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
144 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
145 // eta = 2.6 <-> cos(theta) = 0.99
146 // eta = 3.0 <-> cos(theta) = 0.995
147 cos_theta = cos(theta);
148 }
149 // units from StdHEP : Z [mm] T[mm/c]
150 // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
151 double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
152 double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
153 double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
154 elementZdc->T = timeS;
155
156 // 3) side: which ZDC has been hit?
[355]157 elementZdc->side = sign(particle->Eta);
[374]158
159 // 4) object nature : e.m. (photon) or had (neutron) ?
160 elementZdc->hadronic_hit = (bool) (particle->PID==pN);
[53]161 }
162
163}
164void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
165{
[355]166 int pid=particle->PID;
[53]167
168 TRootRomanPotHits* elementRP220;
[374]169 TRootForwardTaggerHits* elementFP420;
[53]170
171 TLorentzVector genMomentum;
172 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
173 // if forward proton
[100]174 if( (pid == pP) && (particle->Status == 1) && (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
[53]175 {
176 // !!!!!!!! put here particle->CHARGE and particle->MASS
177 H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
178 p1.smearAng();
179 p1.smearPos();
[254]180 p1.setPosition(p1.getX()+DET->RP_cross_x,p1.getY()+DET->RP_cross_y,p1.getTX()-1*kickers_on*DET->RP_cross_ang,p1.getTY(),0);
[53]181 p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
182
183 H_BeamLine *beamline;
184 if(genMomentum.Eta() >0) beamline = beamline1;
185 else beamline = beamline2;
186
187 p1.computePath(beamline,1);
188
189 if(p1.stopped(beamline)) {
190 if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
[100]191 p1.propagate(DET->RP_220_s);
[53]192 elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
193 elementRP220->X = (1E-6)*p1.getX(); // [m]
194 elementRP220->Y = (1E-6)*p1.getY(); // [m]
195 elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
196 elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
197 elementRP220->S = p1.getS(); // [m]
[374]198
199 /* time of flight t is t = T + d/[ cos(theta) v ]
200 // nb: here we assume a straight path to the detector, which is not the case!
201 // this time estimate is always underestimated (while exact for the ZDC case)
202 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3
203 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
204 double tx = atan(particle->Px/particle->Pz);
205 double ty = atan(particle->Py/particle->Pz);
206 double theta = sqrt( pow(tx,2) + pow(ty,2) );
207 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
208 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
209 // eta = 2.6 <-> cos(theta) = 0.99
210 // eta = 3.0 <-> cos(theta) = 0.995
211 cos_theta = cos(theta);
212 }
213 // units from StdHEP : Z [mm] T[mm/c]
214 // units from Delphes : p1.getS [m] speed_of_light [m/s]
215 //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
216 //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
217 */
218 elementRP220->E = p1.getE(); // not yet implemented
219 elementRP220->q2 = -1; // not yet implemented
220 elementRP220->side = sign(particle->Eta);
[53]221
222 } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
[100]223 p1.propagate(DET->RP_420_s);
[374]224 elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
[53]225 elementFP420->X = (1E-6)*p1.getX(); // [m]
226 elementFP420->Y = (1E-6)*p1.getY(); // [m]
227 elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
228 elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
[355]229 elementFP420->S = p1.getS(); // [m]
[374]230
231 // time of flight t is t = T + d/[ cos(theta) v ]
232 // nb: here we assume a straight path to the detector, which is not the case!
233 // this time estimate is always underestimated (while exact for the ZDC case)
234 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3
235 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
236 double tx = atan(particle->Px/particle->Pz);
237 double ty = atan(particle->Py/particle->Pz);
238 double theta = sqrt( pow(tx,2) + pow(ty,2) );
239 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
240 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
241 // eta = 2.6 <-> cos(theta) = 0.99
242 // eta = 3.0 <-> cos(theta) = 0.995
243 cos_theta = cos(theta);
244 }
245 // units from StdHEP : Z [mm] T[mm/c]
246 // units from Delphes : p1.getS [m] speed_of_light [m/s]
247 double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
248 elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
[53]249 elementFP420->E = p1.getE(); // not yet implemented
250 elementFP420->q2 = -1; // not yet implemented
[355]251 elementFP420->side = sign(particle->Eta);
[53]252 }
253
254 }
[355]255 // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
256 // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
[53]257 } // if forward proton
258}
259
260 // Forward particles in CASTOR ?
[355]261 // if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
262 // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
[53]263 //
264 //
[355]265 // } // CASTOR
266 // */
267
Note: See TracBrowser for help on using the repository browser.