Fork me on GitHub

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

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

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