Fork me on GitHub

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

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

not final
File size: 14.2 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
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***********************************************************************/
31
32#include "VeryForward.h"
33#include "H_RomanPot.h"
34#include "PdgParticle.h"
35#include <iostream>
36#include<cmath>
37
38using namespace std;
39
40
41//------------------------------------------------------------------------------
42VeryForward::VeryForward() :
43 DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
44 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
45 relative_energy(true), // should always be true
46 kickers_on(1) // should always be 1
47 {
48 init(); //Initialisation of Hector
49}
50
51VeryForward::VeryForward(const string& DetDatacard) :
52 DET(new RESOLution())
53 {
54 DET->ReadDataCard(DetDatacard);
55 const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
56 beamline1 = new H_BeamLine(1,d_max);
57 beamline2 = new H_BeamLine(1,d_max);
58 init(); //Initialisation of Hector
59 relative_energy = true; // should always be true
60 kickers_on = 1; // should always be 1
61}
62
63VeryForward::VeryForward(const RESOLution * DetDatacard) :
64 DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
65 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
66 relative_energy(true), // should always be true
67 kickers_on(1) // should always be 1
68 {
69 init(); //Initialisation of Hector
70}
71
72VeryForward::VeryForward(const VeryForward& vf) :
73 DET(new RESOLution(*(vf.DET))), d_max(vf.d_max),
74 beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
75 relative_energy(vf.relative_energy),
76 kickers_on(vf.kickers_on) {
77}
78
79VeryForward& VeryForward::operator=(const VeryForward& vf){
80 if (this==&vf) return *this;
81 DET = new RESOLution(*(vf.DET));
82 d_max = vf.d_max;
83 beamline1 = new H_BeamLine(*(vf.beamline1));
84 beamline2 = new H_BeamLine(*(vf.beamline2));
85 relative_energy =vf.relative_energy;
86 kickers_on = vf.kickers_on;
87 return *this;
88}
89
90
91void VeryForward::init() {
92 //Initialisation of Hector
93 static unsigned int counter;
94 counter =0;
95 relative_energy = true; // should always be true
96 kickers_on = 1; // should always be 1
97 beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
98 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
99 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); // RP 220m, 2mm, beam 1
100 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); // RP 420m, 4mm, beam 1
101 beamline1->add(rp220_1);
102 beamline1->add(rp420_1);
103
104 beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
105 beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x);
106 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);// RP 220m, 2mm, beam 2
107 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);// RP 420m, 4mm, beam 2
108 beamline2->add(rp220_2);
109 beamline2->add(rp420_2);
110 // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
111 // do not put explicit delete
112}
113
114
115
116void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
117{
118 TRootZdcHits *elementZdc;
119 float energy = particle->E;
120
121 // Zero degree calorimeter, for forward neutrons and photons
122 if (particle->Status ==1 && ( (particle->PID==pN && energy>DET->ZDC_n_E) ||
123 (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
124 && fabs(particle->Eta) > DET->VFD_min_zdc ) {
125 elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
126
127
128 elementZdc->pid = particle->PID;
129
130
131 // for compatibility with 'old' version
132 TLorentzVector genMomentum;
133 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
134 elementZdc->Set(genMomentum);
135 // ******************
136
137
138 //particle->print();
139
140 // 1) energy smearing
141 float energyS = -1.;
142 if (particle->PID == pGAMMA)
143 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
144 pow(DET->ELG_Czdc*particle->E,2) +
145 pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
146 else // smearing with hadronic resolution
147 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
148 pow(DET->HAD_Czdc*particle->E,2) +
149 pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
150 elementZdc->E = energyS;
151
152
153 // 2) time of flight t is t = T + d/[ cos(theta) v ]
154 float cos_theta = 1; //very good approximation, if eta_zdc >3
155 if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
156 double tx = atan(particle->Px/particle->Pz);
157 double ty = atan(particle->Py/particle->Pz);
158 double theta = sqrt( pow(tx,2) + pow(ty,2) );
159 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
160 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
161 // eta = 2.6 <-> cos(theta) = 0.99
162 // eta = 3.0 <-> cos(theta) = 0.995
163 cos_theta = cos(theta);
164 }
165 // units from StdHEP : Z [mm] T[mm/c]
166 // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
167 double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
168 double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
169 double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
170 elementZdc->T = timeS;
171
172 // 3) side: which ZDC has been hit?
173 elementZdc->side = sign(particle->Eta);
174
175 // 4) object nature : e.m. (photon) or had (neutron) ?
176 //elementZdc->hadronic_hit = (bool) (particle->PID==pN);
177 }
178}
179
180
181void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
182{
183 float charge = particle->Charge, mass = particle->M;
184 if (mass<-999) { // unitialised!
185 PdgParticle pdg_part = DET->PDGtable[particle->PID];
186 charge = pdg_part.charge();
187 mass = pdg_part.mass();
188 }
189 //if(particle->Charge!=1) return; // only particles with Q=+1 can hope to reach RP200/FP420
190
191 TRootRomanPotHits* elementRP220;
192 //TRootForwardTaggerHits* elementFP420;
193 TRootRomanPotHits* elementFP420;
194
195 TLorentzVector genMomentum;
196 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
197
198 // to go faster, why not rejecting particles already going into the ZDC?
199 if( particle->PID == pP)
200 if( (particle->Status == 1) && (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
201 {
202 //cout << "VeryForward :: M = " << mass << "\t Q = " << charge << "\t\t " << particle->PID << endl;
203 H_BeamParticle p1(mass,charge);
204 p1.smearAng(); p1.smearPos(); // vertex smearing
205 cout << "x = " << p1.getX() + DET->RP_cross_x
206 << " y= " << p1.getY() + DET->RP_cross_y
207 << " tx= " << p1.getTX() - kickers_on*DET->RP_cross_ang
208 << " ty=" << p1.getTY() << endl;
209 p1.setPosition(p1.getX()+DET->RP_cross_x,p1.getY()+DET->RP_cross_y,p1.getTX()-kickers_on*DET->RP_cross_ang,p1.getTY(),0);
210 //p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
211 p1.setE(particle->E);
212
213 H_BeamLine *beamline;
214 if(genMomentum.Eta() >0) beamline = beamline1;
215 else beamline = beamline2;
216
217 p1.computePath(beamline,1);
218
219 if(p1.stopped(beamline)) {
220 if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
221 static unsigned int counter;
222 counter++;
223 if (counter==1) {
224 p1.getPath(0,"p1path.txt");
225 cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge
226 << "\t" << mass << "=" << particle->M
227 << "\t E=" << particle->E
228 << endl;
229 }
230 p1.propagate(DET->RP_220_s);
231 elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
232 elementRP220->X = (1E-6)*p1.getX(); // [m]
233 elementRP220->Y = (1E-6)*p1.getY(); // [m]
234 elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
235 elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
236 elementRP220->S = p1.getS(); // [m]
237
238 /* time of flight t is t = T + d/[ cos(theta) v ]
239 // nb: here we assume a straight path to the detector, which is not the case!
240 // this time estimate is always underestimated (while exact for the ZDC case)
241 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3
242 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
243 double tx = atan(particle->Px/particle->Pz);
244 double ty = atan(particle->Py/particle->Pz);
245 double theta = sqrt( pow(tx,2) + pow(ty,2) );
246 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
247 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
248 // eta = 2.6 <-> cos(theta) = 0.99
249 // eta = 3.0 <-> cos(theta) = 0.995
250 cos_theta = cos(theta);
251 }
252 // units from StdHEP : Z [mm] T[mm/c]
253 // units from Delphes : p1.getS [m] speed_of_light [m/s]
254 //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
255 //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
256 */
257 elementRP220->E = p1.getE(); // not yet implemented
258 elementRP220->q2 = -1; // not yet implemented
259 elementRP220->side = sign(particle->Eta);
260
261
262 elementRP220->pid = particle->PID;
263
264 } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
265 p1.propagate(DET->RP_420_s);
266 //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
267 elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
268 elementFP420->X = (1E-6)*p1.getX(); // [m]
269 elementFP420->Y = (1E-6)*p1.getY(); // [m]
270 elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
271 elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
272 elementFP420->S = p1.getS(); // [m]
273
274 // time of flight t is t = T + d/[ cos(theta) v ]
275 // nb: here we assume a straight path to the detector, which is not the case!
276 // this time estimate is always underestimated (while exact for the ZDC case)
277 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3
278 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
279 double tx = atan(particle->Px/particle->Pz);
280 double ty = atan(particle->Py/particle->Pz);
281 double theta = sqrt( pow(tx,2) + pow(ty,2) );
282 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
283 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
284 // eta = 2.6 <-> cos(theta) = 0.99
285 // eta = 3.0 <-> cos(theta) = 0.995
286 cos_theta = cos(theta);
287 }
288 // units from StdHEP : Z [mm] T[mm/c]
289 // units from Delphes : p1.getS [m] speed_of_light [m/s]
290 double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
291 elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
292 elementFP420->E = p1.getE(); // not yet implemented
293 elementFP420->q2 = -1; // not yet implemented
294 elementFP420->side = sign(particle->Eta);
295
296
297 elementFP420->pid = particle->PID;
298 }
299
300 }
301 // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
302 // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
303 } // if forward proton
304}
305
306 // Forward particles in CASTOR ?
307 // if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
308 // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
309 //
310 //
311 // } // CASTOR
312 // */
313
Note: See TracBrowser for help on using the repository browser.