Fork me on GitHub

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

Last change on this file since 256 was 254, checked in by severine ovyn, 16 years ago

add shifts X,Y,angl
File size: 7.4 KB
Line 
1 /*
2 * ---- Delphes ----
3 * A Fast Simulator for general purpose LHC detector
4 * S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5 *
6 * Center for Particle Physics and Phenomenology (CP3)
7 * Universite Catholique de Louvain (UCL)
8 * Louvain-la-Neuve, Belgium
9 * */
10
11#include "VeryForward.h"
12#include "H_RomanPot.h"
13#include <iostream>
14#include<cmath>
15
16using namespace std;
17
18
19//------------------------------------------------------------------------------
20
21VeryForward::VeryForward() {
22 DET = new RESOLution();
23 beamline1 = new H_BeamLine(1,500.);
24 beamline2 = new H_BeamLine(1,500.);
25 init();
26 //Initialisation of Hector
27 relative_energy = true; // should always be true
28 kickers_on = 1; // should always be 1
29
30}
31
32VeryForward::VeryForward(const string& DetDatacard) {
33 DET = new RESOLution();
34 DET->ReadDataCard(DetDatacard);
35 beamline1 = new H_BeamLine(1,500.);
36 beamline2 = new H_BeamLine(1,500.);
37 init();
38 //Initialisation of Hector
39 relative_energy = true; // should always be true
40 kickers_on = 1; // should always be 1
41
42}
43
44VeryForward::VeryForward(const RESOLution * DetDatacard) {
45 DET = new RESOLution(*DetDatacard);
46 beamline2 = new H_BeamLine(1,500.);
47 beamline1 = new H_BeamLine(1,500.);
48
49 init();
50 //Initialisation of Hector
51 relative_energy = true; // should always be true
52 kickers_on = 1; // should always be 1
53
54}
55
56VeryForward::VeryForward(const VeryForward& vf) {
57 DET = new RESOLution(*(vf.DET));
58 beamline1 = new H_BeamLine(*(vf.beamline1));
59 beamline2 = new H_BeamLine(*(vf.beamline2));
60}
61
62VeryForward& VeryForward::operator=(const VeryForward& vf){
63 if (this==&vf) return *this;
64 DET = new RESOLution(*(vf.DET));
65 beamline1 = new H_BeamLine(*(vf.beamline1));
66 beamline2 = new H_BeamLine(*(vf.beamline2));
67 return *this;
68}
69
70
71void VeryForward::init() {
72 //Initialisation of Hector
73 relative_energy = true; // should always be true
74 kickers_on = 1; // should always be 1
75 // user should provide : (1) optics file for each beamline, and IPname,
76 // and offset data (s,x) for optical elements
77 beamline1->fill(DET->RP_beam1Card,1,"IP5");
78 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
79 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*(1E6)); // RP 220m, 2mm, beam 1
80 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*(1E6)); // RP 420m, 4mm, beam 1
81 beamline1->add(rp220_1);
82 beamline1->add(rp420_1);
83
84 beamline2->fill(DET->RP_beam2Card,-1,"IP5");
85 beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x);
86 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*(1E6));// RP 220m, 2mm, beam 2
87 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*(1E6));// RP 420m, 4mm, beam 2
88 beamline2->add(rp220_2);
89 beamline2->add(rp420_2);
90 // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
91 // do not put explicit delete
92}
93
94
95void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
96{
97 int pid=abs(particle->PID);
98 float eta=fabs(particle->Eta);
99
100
101 TRootZdcHits *elementZdc;
102 TLorentzVector genMomentum;
103 // Zero degree calorimeter, for forward neutrons and photons
104 if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->VFD_min_zdc ) {
105 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
106 // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
107 // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
108 // !!!!!!!!! pas de smearing ! on garde trop d'info !
109 elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
110 elementZdc->Set(genMomentum);
111
112 // time of flight t is t = T + d/[ cos(theta) v ]
113 //double tx = acos(particle->Px/particle->Pz);
114 //double ty = acos(particle->Py/particle->Pz);
115 //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
116 //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
117 double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-6));
118 // assumes also that the emission angle is so small that 1/(cos theta) = 1
119 elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
120 elementZdc->side = sign(eta);
121
122
123 }
124
125}
126void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
127{
128 int pid=abs(particle->PID);
129 float eta=fabs(particle->Eta);
130
131 TRootRomanPotHits* elementRP220;
132 TRootRomanPotHits* elementFP420;
133
134 TLorentzVector genMomentum;
135 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
136 // if forward proton
137 if( (pid == pP) && (particle->Status == 1) && (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
138 {
139 // !!!!!!!! put here particle->CHARGE and particle->MASS
140 H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
141 p1.smearAng();
142 p1.smearPos();
143 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);
144 p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
145
146 H_BeamLine *beamline;
147 if(genMomentum.Eta() >0) beamline = beamline1;
148 else beamline = beamline2;
149
150 p1.computePath(beamline,1);
151
152 if(p1.stopped(beamline)) {
153 if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
154 p1.propagate(DET->RP_220_s);
155 elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
156 elementRP220->X = (1E-6)*p1.getX(); // [m]
157 elementRP220->Y = (1E-6)*p1.getY(); // [m]
158 elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
159 elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
160 elementRP220->S = p1.getS(); // [m]
161 elementRP220->T = -1; // not yet implemented
162 elementRP220->E = p1.getE(); // not yet implemented
163 elementRP220->q2 = -1; // not yet implemented
164 elementRP220->side = sign(eta);
165
166 } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
167 p1.propagate(DET->RP_420_s);
168 elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
169 elementFP420->X = (1E-6)*p1.getX(); // [m]
170 elementFP420->Y = (1E-6)*p1.getY(); // [m]
171 elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
172 elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
173 elementFP420->S = p1.getS(); // [m]
174 elementFP420->T = -1; // not yet implemented
175 elementFP420->E = p1.getE(); // not yet implemented
176 elementFP420->q2 = -1; // not yet implemented
177 elementFP420->side = sign(eta);
178 }
179
180 }
181 // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
182 // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
183 } // if forward proton
184}
185
186 // Forward particles in CASTOR ?
187 // /* if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
188 // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
189 //
190 //
191 // } // CASTOR
192 // */
193 //
Note: See TracBrowser for help on using the repository browser.