Fork me on GitHub

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

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

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