Fork me on GitHub

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

Last change on this file since 1078 was 512, checked in by Xavier Rouby, 15 years ago

new commented lines for the misalignment studies
File size: 15.4 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** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34#include "VeryForward.h"
35#include "PdgParticle.h"
36#include "H_RomanPot.h"
37
38#include <iostream>
39#include <fstream>
40#include<cmath>
41
42using namespace std;
43
44/* Notes on the correct initialisation for Hector
45 * -- these notes apply to the LHC beamlines
46 *
47 * beam1 : forward direction is for increasing 's' values
48 * beamline1 = new H_BeamLine(1,...);
49 * beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
50 * beam2 : forward direction is for decreasing 's' values
51 * beamline2 = new H_BeamLine(-1,...);
52 * beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
53 *
54 * relative_energy should be false -- kickers_on should be 1
55 *
56 */
57
58//------------------------------------------------------------------------------
59VeryForward::VeryForward() :
60 DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
61 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
62 rel_energy(true), // should always be true
63 kickers(1) // should always be 1
64 {
65 init(); //Initialisation of Hector
66}
67
68VeryForward::VeryForward(const string& DetDatacard) :
69 DET(new RESOLution())
70 {
71 DET->ReadDataCard(DetDatacard);
72 const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
73 beamline1 = new H_BeamLine(1,d_max);
74 beamline2 = new H_BeamLine(-1,d_max);
75 init(); //Initialisation of Hector
76 rel_energy = true; // should always be true
77 kickers = 1; // should always be 1
78}
79
80VeryForward::VeryForward(const RESOLution * DetDatacard) :
81 DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
82 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
83 rel_energy(true), // should always be true
84 kickers(1) // should always be 1
85 {
86 init(); //Initialisation of Hector
87}
88
89VeryForward::VeryForward(const VeryForward& vf) :
90 DET(new RESOLution(*(vf.DET))), d_max(vf.d_max),
91 beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
92 rel_energy(vf.rel_energy),
93 kickers(vf.kickers) {
94}
95
96VeryForward& VeryForward::operator=(const VeryForward& vf){
97 if (this==&vf) return *this;
98 DET = new RESOLution(*(vf.DET));
99 d_max = vf.d_max;
100 beamline1 = new H_BeamLine(*(vf.beamline1));
101 beamline2 = new H_BeamLine(*(vf.beamline2));
102 rel_energy =vf.rel_energy;
103 kickers = vf.kickers;
104 return *this;
105}
106
107
108void VeryForward::init() {
109 //Initialisation of Hector
110 static unsigned int counter;
111 counter =0;
112 extern bool relative_energy;
113 extern int kickers_on;
114 relative_energy = rel_energy; // should always be true
115 kickers_on = kickers; // should always be 1
116 beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
117 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
118 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6);
119 // RP 220m, 2mm, beam 1
120 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6);
121 // RP 420m, 4mm, beam 1
122 //rp220_1->printProperties();
123 //rp420_1->printProperties();
124 beamline1->add(rp220_1);
125 beamline1->add(rp420_1);
126 //beamline1->alignElement("\"MQXA.1R5\"",+0.0005,0);
127 //beamline1->alignElement("\"MQM.9R5.B1\"",-0.0001,0);
128 //beamline1->alignElement("\"MQML.5R5.B1\"",+0.0001,0);
129
130 beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
131 beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
132 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);
133 // RP 220m, 2mm, beam 2
134 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);
135 // RP 420m, 4mm, beam 2
136 //rp220_2->printProperties();
137 //rp420_2->printProperties();
138 beamline2->add(rp220_2);
139 beamline2->add(rp420_2);
140 //beamline2->alignElement("\"MQXA.1L5\"",+0.0005,0);
141 // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
142 // do not put explicit delete
143}
144
145
146float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) {
147 // time of flight t is t = T + d/[ cos(theta) v ]
148 float cos_theta = 1; //very good approximation, if detector_etamin >3
149 if (detector_etamin<3) { // if smaller eta -> make the complete calculation
150 double tx = atan(particle->Px/particle->Pz);
151 double ty = atan(particle->Py/particle->Pz);
152 double theta = sqrt( pow(tx,2) + pow(ty,2) );
153 // cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
154 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
155 // eta = 2.6 <-> cos(theta) = 0.99
156 // eta = 3.0 <-> cos(theta) = 0.995
157 cos_theta = cos(theta);
158 }
159 // units from StdHEP : Z [mm] T[mm/c]
160 // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
161 double flight_distance = (detector_s - particle->Z*(1E-3))/cos_theta ;
162 // assumes highly relativistic particles
163 double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light;
164 double timeS = gRandom->Gaus(flight_time,detector_t_resolution);
165 return timeS;
166}
167
168
169
170void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
171{
172 TRootZdcHits *elementZdc;
173 float energy = particle->E;
174
175 // Zero degree calorimeter, for forward neutrons and photons
176 if (particle->Status ==1 && ( (particle->PID==pN && energy>DET->ZDC_n_E) ||
177 (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
178 && fabs(particle->Eta) > DET->VFD_min_zdc ) {
179 elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
180
181 TLorentzVector genMomentum;
182 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
183 elementZdc->Set(genMomentum); // initialises the gen-level data
184 elementZdc->pid = particle->PID;
185
186 // 1) energy smearing
187 float energyS = -1.;
188 if (particle->PID == pGAMMA)
189 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
190 pow(DET->ELG_Czdc*particle->E,2) +
191 pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
192 else // smearing with hadronic resolution
193 energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
194 pow(DET->HAD_Czdc*particle->E,2) +
195 pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
196 elementZdc->E = energyS;
197
198 // 2) time of flight t is t = T + d/[ cos(theta) v ] + detector smearing on time
199 elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution);
200
201 // 3) side: which ZDC has been hit?
202 elementZdc->side = sign(particle->Eta);
203
204 // 4) object nature : e.m. (photon) or had (neutron) ?
205 elementZdc->hadronic_hit = (bool) (particle->PID!=pGAMMA);
206 } // if neutrons or photons over E_threshold
207
208}
209
210
211void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
212{
213 if(particle->Status != 1) return; // reject particles that are not final ones
214 extern bool relative_energy;
215 relative_energy = rel_energy;
216 extern int kickers_on;
217 kickers_on = kickers;
218
219 float charge = particle->Charge, mass = particle->M;
220 //float charge, mass, ctau;
221 //charge = mass = ctau = UNDEFINED;
222 if (mass<-999) { // unitialised!
223 PdgParticle pdg_part = DET->PDGtable[particle->PID];
224 charge = pdg_part.charge(); // e+
225 mass = pdg_part.mass(); // GeV
226 // ctau = pdg_part.ctau(); // m
227 // cout << "ctau = " << ctau << endl;
228 }
229
230
231 if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420
232 //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl;
233 TRootRomanPotHits* elementRP220;
234 //TRootForwardTaggerHits* elementFP420;
235 TRootRomanPotHits* elementFP420;
236
237 TLorentzVector genMomentum;
238 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
239
240 // to go faster, why not rejecting particles already going into the ZDC?
241
242 // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ;
243 // so if RP/FP too far away, the particle must be a proton or a mu+
244 if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return;
245 if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd )
246 {
247 H_BeamParticle p1(mass,charge);
248 double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad
249 double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad
250 //p1.smearAng(); p1.smearPos(); // vertex smearing do not put it here !!!
251
252 /* cout << "x = " << particle->X << " + " << p1.getX() << " + " << DET->RP_cross_x
253 << " y= " << particle->Y << " + " << p1.getY() << " + " << DET->RP_cross_y
254 << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x
255 << " ty= " << ty << " + " << p1.getTY() << " - " << kickers_on*DET->RP_cross_ang_y
256 << " z= " << particle->Z << endl;*/
257
258 // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done
259 // all in micrometers or microradians
260
261 // for checking purposes
262 /*p1.smearAng(); p1.smearPos(); tx=0; ty=0; // to be removed !!!! test only !!!!
263 ofstream xdist("xdist",ios::ios_base::app);
264 xdist << endl;
265 xdist << (1E3)*particle->X + p1.getX() + DET->RP_cross_x << " = " << (1E3)*particle->X << " + " << p1.getX() << " + " << DET->RP_cross_x << endl;
266 xdist << (1E3)*particle->Y + p1.getY() + DET->RP_cross_y << " = " << (1E3)*particle->Y << " + " << p1.getY() << " + " << DET->RP_cross_y << endl;
267 xdist << tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x << " = " << tx << " + " << p1.getTX() << " - " << kickers_on << "*" << DET->RP_cross_ang_x<< endl;
268 xdist << ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y << " = " << ty << " + " << p1.getTY() << " - " << kickers_on << "*" << DET->RP_cross_ang_y<< endl;
269 xdist.close();
270 */
271
272
273 p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x,
274 (1E3)*particle->Y + p1.getY() + DET->RP_cross_y,
275 tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x,
276 ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y,
277 0*(1E3)*particle->Z);
278 p1.setE(particle->E);
279
280 H_BeamLine *beamline;
281 if(genMomentum.Eta() >0) beamline = beamline1;
282 else beamline = beamline2;
283
284 p1.computePath(beamline,1);
285
286 if(p1.stopped(beamline)) {
287 // roman pots at 220 m
288 if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
289
290 /*static unsigned int counter;
291 counter++;
292 if (counter==1) {
293 p1.getPath(0,"p1path.txt");
294 cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge
295 << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl;
296 }*/
297
298 p1.propagate(DET->RP_220_s);
299 elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
300
301 // detector measurements
302 elementRP220->X = (1E-6)*p1.getX(); // [m]
303 elementRP220->Y = (1E-6)*p1.getY(); // [m]
304 elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
305 elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
306 elementRP220->S = p1.getS(); // [m]
307 elementRP220->T = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution);
308 elementRP220->side = sign(particle->Eta);
309
310 // reconstructed data
311 float sE = p1.getE(); // apply the smearing here!!!
312 elementRP220->E = sE; // not yet implemented
313 elementRP220->q2 = UNDEFINED; // not yet implemented
314
315 // generator level data
316 elementRP220->pid = particle->PID;
317 elementRP220->Set(genMomentum);
318 } // if RP220
319
320 // proton taggers at 420 m
321 else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
322
323 p1.propagate(DET->RP_420_s);
324 elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
325 //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
326
327 // detector measurements
328 elementFP420->X = (1E-6)*p1.getX(); // [m]
329 elementFP420->Y = (1E-6)*p1.getY(); // [m]
330 elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
331 elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
332 elementFP420->S = p1.getS(); // [m]
333 elementFP420->T = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution);
334 elementFP420->side = sign(particle->Eta);
335
336 // reconstructed data
337 elementFP420->E = p1.getE(); // not yet implemented
338 elementFP420->q2 = UNDEFINED; // not yet implemented
339
340 // generator level data
341 elementFP420->pid = particle->PID;
342 elementFP420->Set(genMomentum);
343 } // if FP420
344
345 } // if stopped
346 } // if forward proton
347
348}
Note: See TracBrowser for help on using the repository browser.