| 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 |
|
|---|
| 42 | using 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 | //------------------------------------------------------------------------------
|
|---|
| 59 | VeryForward::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 |
|
|---|
| 68 | VeryForward::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 |
|
|---|
| 80 | VeryForward::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 |
|
|---|
| 89 | VeryForward::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 |
|
|---|
| 96 | VeryForward& 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 |
|
|---|
| 108 | void 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 |
|
|---|
| 127 | beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
|
|---|
| 128 | beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
|
|---|
| 129 | H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);
|
|---|
| 130 | // RP 220m, 2mm, beam 2
|
|---|
| 131 | H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);
|
|---|
| 132 | // RP 420m, 4mm, beam 2
|
|---|
| 133 | //rp220_2->printProperties();
|
|---|
| 134 | //rp420_2->printProperties();
|
|---|
| 135 | beamline2->add(rp220_2);
|
|---|
| 136 | beamline2->add(rp420_2);
|
|---|
| 137 | // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
|
|---|
| 138 | // do not put explicit delete
|
|---|
| 139 | }
|
|---|
| 140 |
|
|---|
| 141 |
|
|---|
| 142 | float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) {
|
|---|
| 143 | // time of flight t is t = T + d/[ cos(theta) v ]
|
|---|
| 144 | float cos_theta = 1; //very good approximation, if detector_etamin >3
|
|---|
| 145 | if (detector_etamin<3) { // if smaller eta -> make the complete calculation
|
|---|
| 146 | double tx = atan(particle->Px/particle->Pz);
|
|---|
| 147 | double ty = atan(particle->Py/particle->Pz);
|
|---|
| 148 | double theta = sqrt( pow(tx,2) + pow(ty,2) );
|
|---|
| 149 | // cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl;
|
|---|
| 150 | // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
|
|---|
| 151 | // eta = 2.6 <-> cos(theta) = 0.99
|
|---|
| 152 | // eta = 3.0 <-> cos(theta) = 0.995
|
|---|
| 153 | cos_theta = cos(theta);
|
|---|
| 154 | }
|
|---|
| 155 | // units from StdHEP : Z [mm] T[mm/c]
|
|---|
| 156 | // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
|
|---|
| 157 | double flight_distance = (detector_s - particle->Z*(1E-3))/cos_theta ;
|
|---|
| 158 | // assumes highly relativistic particles
|
|---|
| 159 | double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light;
|
|---|
| 160 | double timeS = gRandom->Gaus(flight_time,detector_t_resolution);
|
|---|
| 161 | return timeS;
|
|---|
| 162 | }
|
|---|
| 163 |
|
|---|
| 164 |
|
|---|
| 165 |
|
|---|
| 166 | void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
|
|---|
| 167 | {
|
|---|
| 168 | TRootZdcHits *elementZdc;
|
|---|
| 169 | float energy = particle->E;
|
|---|
| 170 |
|
|---|
| 171 | // Zero degree calorimeter, for forward neutrons and photons
|
|---|
| 172 | if (particle->Status ==1 && ( (particle->PID==pN && energy>DET->ZDC_n_E) ||
|
|---|
| 173 | (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
|
|---|
| 174 | && fabs(particle->Eta) > DET->VFD_min_zdc ) {
|
|---|
| 175 | elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
|
|---|
| 176 |
|
|---|
| 177 | TLorentzVector genMomentum;
|
|---|
| 178 | genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
|
|---|
| 179 | elementZdc->Set(genMomentum); // initialises the gen-level data
|
|---|
| 180 | elementZdc->pid = particle->PID;
|
|---|
| 181 |
|
|---|
| 182 | // 1) energy smearing
|
|---|
| 183 | float energyS = -1.;
|
|---|
| 184 | if (particle->PID == pGAMMA)
|
|---|
| 185 | energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
|
|---|
| 186 | pow(DET->ELG_Czdc*particle->E,2) +
|
|---|
| 187 | pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
|
|---|
| 188 | else // smearing with hadronic resolution
|
|---|
| 189 | energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
|
|---|
| 190 | pow(DET->HAD_Czdc*particle->E,2) +
|
|---|
| 191 | pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
|
|---|
| 192 | elementZdc->E = energyS;
|
|---|
| 193 |
|
|---|
| 194 | // 2) time of flight t is t = T + d/[ cos(theta) v ] + detector smearing on time
|
|---|
| 195 | elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution);
|
|---|
| 196 |
|
|---|
| 197 | // 3) side: which ZDC has been hit?
|
|---|
| 198 | elementZdc->side = sign(particle->Eta);
|
|---|
| 199 |
|
|---|
| 200 | // 4) object nature : e.m. (photon) or had (neutron) ?
|
|---|
| 201 | elementZdc->hadronic_hit = (bool) (particle->PID!=pGAMMA);
|
|---|
| 202 | } // if neutrons or photons over E_threshold
|
|---|
| 203 |
|
|---|
| 204 | }
|
|---|
| 205 |
|
|---|
| 206 |
|
|---|
| 207 | void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
|
|---|
| 208 | {
|
|---|
| 209 | if(particle->Status != 1) return; // reject particles that are not final ones
|
|---|
| 210 | extern bool relative_energy;
|
|---|
| 211 | relative_energy = rel_energy;
|
|---|
| 212 | extern int kickers_on;
|
|---|
| 213 | kickers_on = kickers;
|
|---|
| 214 |
|
|---|
| 215 | float charge = particle->Charge, mass = particle->M;
|
|---|
| 216 | //float charge, mass, ctau;
|
|---|
| 217 | //charge = mass = ctau = UNDEFINED;
|
|---|
| 218 | if (mass<-999) { // unitialised!
|
|---|
| 219 | PdgParticle pdg_part = DET->PDGtable[particle->PID];
|
|---|
| 220 | charge = pdg_part.charge(); // e+
|
|---|
| 221 | mass = pdg_part.mass(); // GeV
|
|---|
| 222 | // ctau = pdg_part.ctau(); // m
|
|---|
| 223 | // cout << "ctau = " << ctau << endl;
|
|---|
| 224 | }
|
|---|
| 225 |
|
|---|
| 226 |
|
|---|
| 227 | if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420
|
|---|
| 228 | //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl;
|
|---|
| 229 | TRootRomanPotHits* elementRP220;
|
|---|
| 230 | //TRootForwardTaggerHits* elementFP420;
|
|---|
| 231 | TRootRomanPotHits* elementFP420;
|
|---|
| 232 |
|
|---|
| 233 | TLorentzVector genMomentum;
|
|---|
| 234 | genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
|
|---|
| 235 |
|
|---|
| 236 | // to go faster, why not rejecting particles already going into the ZDC?
|
|---|
| 237 |
|
|---|
| 238 | // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ;
|
|---|
| 239 | // so if RP/FP too far away, the particle must be a proton or a mu+
|
|---|
| 240 | if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return;
|
|---|
| 241 |
|
|---|
| 242 | if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd )
|
|---|
| 243 | {
|
|---|
| 244 | H_BeamParticle p1(mass,charge);
|
|---|
| 245 | double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad
|
|---|
| 246 | double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad
|
|---|
| 247 | //p1.smearAng(); p1.smearPos(); // vertex smearing do not put it here !!!
|
|---|
| 248 |
|
|---|
| 249 | /* cout << "x = " << particle->X << " + " << p1.getX() << " + " << DET->RP_cross_x
|
|---|
| 250 | << " y= " << particle->Y << " + " << p1.getY() << " + " << DET->RP_cross_y
|
|---|
| 251 | << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x
|
|---|
| 252 | << " ty= " << ty << " + " << p1.getTY() << " - " << kickers_on*DET->RP_cross_ang_y
|
|---|
| 253 | << " z= " << particle->Z << endl;*/
|
|---|
| 254 |
|
|---|
| 255 | // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done
|
|---|
| 256 | // all in micrometers or microradians
|
|---|
| 257 | p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x,
|
|---|
| 258 | (1E3)*particle->Y + p1.getY() + DET->RP_cross_y,
|
|---|
| 259 | tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x,
|
|---|
| 260 | ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y,
|
|---|
| 261 | 0*(1E3)*particle->Z);
|
|---|
| 262 | p1.setE(particle->E);
|
|---|
| 263 |
|
|---|
| 264 | H_BeamLine *beamline;
|
|---|
| 265 | if(genMomentum.Eta() >0) beamline = beamline1;
|
|---|
| 266 | else beamline = beamline2;
|
|---|
| 267 |
|
|---|
| 268 | p1.computePath(beamline,1);
|
|---|
| 269 |
|
|---|
| 270 | if(p1.stopped(beamline)) {
|
|---|
| 271 | // roman pots at 220 m
|
|---|
| 272 | if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
|
|---|
| 273 |
|
|---|
| 274 | /*static unsigned int counter;
|
|---|
| 275 | counter++;
|
|---|
| 276 | if (counter==1) {
|
|---|
| 277 | p1.getPath(0,"p1path.txt");
|
|---|
| 278 | cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge
|
|---|
| 279 | << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl;
|
|---|
| 280 | }*/
|
|---|
| 281 |
|
|---|
| 282 | p1.propagate(DET->RP_220_s);
|
|---|
| 283 | elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
|
|---|
| 284 |
|
|---|
| 285 | // detector measurements
|
|---|
| 286 | elementRP220->X = (1E-6)*p1.getX(); // [m]
|
|---|
| 287 | elementRP220->Y = (1E-6)*p1.getY(); // [m]
|
|---|
| 288 | elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
|
|---|
| 289 | elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
|
|---|
| 290 | elementRP220->S = p1.getS(); // [m]
|
|---|
| 291 | elementRP220->T = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution);
|
|---|
| 292 | elementRP220->side = sign(particle->Eta);
|
|---|
| 293 |
|
|---|
| 294 | // reconstructed data
|
|---|
| 295 | float sE = p1.getE(); // apply the smearing here!!!
|
|---|
| 296 | elementRP220->E = sE; // not yet implemented
|
|---|
| 297 | elementRP220->q2 = UNDEFINED; // not yet implemented
|
|---|
| 298 |
|
|---|
| 299 | // generator level data
|
|---|
| 300 | elementRP220->pid = particle->PID;
|
|---|
| 301 | elementRP220->Set(genMomentum);
|
|---|
| 302 | } // if RP220
|
|---|
| 303 |
|
|---|
| 304 | // proton taggers at 420 m
|
|---|
| 305 | else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
|
|---|
| 306 |
|
|---|
| 307 | p1.propagate(DET->RP_420_s);
|
|---|
| 308 | elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
|
|---|
| 309 | //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
|
|---|
| 310 |
|
|---|
| 311 | // detector measurements
|
|---|
| 312 | elementFP420->X = (1E-6)*p1.getX(); // [m]
|
|---|
| 313 | elementFP420->Y = (1E-6)*p1.getY(); // [m]
|
|---|
| 314 | elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
|
|---|
| 315 | elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
|
|---|
| 316 | elementFP420->S = p1.getS(); // [m]
|
|---|
| 317 | elementFP420->T = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution);
|
|---|
| 318 | elementFP420->side = sign(particle->Eta);
|
|---|
| 319 |
|
|---|
| 320 | // reconstructed data
|
|---|
| 321 | elementFP420->E = p1.getE(); // not yet implemented
|
|---|
| 322 | elementFP420->q2 = UNDEFINED; // not yet implemented
|
|---|
| 323 |
|
|---|
| 324 | // generator level data
|
|---|
| 325 | elementFP420->pid = particle->PID;
|
|---|
| 326 | elementFP420->Set(genMomentum);
|
|---|
| 327 | } // if FP420
|
|---|
| 328 |
|
|---|
| 329 | } // if stopped
|
|---|
| 330 | } // if forward proton
|
|---|
| 331 |
|
|---|
| 332 | }
|
|---|
