Fork me on GitHub

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

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

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