[53] | 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 |
|
---|
| 23 | using namespace std;
|
---|
| 24 |
|
---|
| 25 |
|
---|
| 26 | //------------------------------------------------------------------------------
|
---|
| 27 |
|
---|
| 28 | VeryForward::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 |
|
---|
| 55 | void 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 |
|
---|
| 83 | }
|
---|
| 84 |
|
---|
| 85 | }
|
---|
| 86 | void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
|
---|
| 87 | {
|
---|
| 88 | int pid=abs(particle->PID);
|
---|
| 89 | float eta=fabs(particle->Eta);
|
---|
| 90 |
|
---|
| 91 | TRootRomanPotHits* elementRP220;
|
---|
| 92 | TRootRomanPotHits* elementFP420;
|
---|
| 93 |
|
---|
| 94 | TLorentzVector genMomentum;
|
---|
| 95 | genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
|
---|
| 96 | // if forward proton
|
---|
| 97 | if( (pid == pP) && (particle->Status == 1) && (fabs(genMomentum.Eta()) > MAX_CALO_FWD) )
|
---|
| 98 | {
|
---|
| 99 | // !!!!!!!! put here particle->CHARGE and particle->MASS
|
---|
| 100 | H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
|
---|
| 101 | p1.smearAng();
|
---|
| 102 | p1.smearPos();
|
---|
| 103 | p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
|
---|
| 104 | p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
|
---|
| 105 |
|
---|
| 106 | H_BeamLine *beamline;
|
---|
| 107 | if(genMomentum.Eta() >0) beamline = beamline1;
|
---|
| 108 | else beamline = beamline2;
|
---|
| 109 |
|
---|
| 110 | p1.computePath(beamline,1);
|
---|
| 111 |
|
---|
| 112 | if(p1.stopped(beamline)) {
|
---|
| 113 | if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
|
---|
| 114 | p1.propagate(RP220_S);
|
---|
| 115 | elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
|
---|
| 116 | elementRP220->X = (1E-6)*p1.getX(); // [m]
|
---|
| 117 | elementRP220->Y = (1E-6)*p1.getY(); // [m]
|
---|
| 118 | elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
|
---|
| 119 | elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
|
---|
| 120 | elementRP220->S = p1.getS(); // [m]
|
---|
| 121 | elementRP220->T = -1; // not yet implemented
|
---|
| 122 | elementRP220->E = p1.getE(); // not yet implemented
|
---|
| 123 | elementRP220->q2 = -1; // not yet implemented
|
---|
| 124 | elementRP220->side = sign(eta);
|
---|
| 125 |
|
---|
| 126 | } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
|
---|
| 127 | p1.propagate(FP420_S);
|
---|
| 128 | elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
|
---|
| 129 | elementFP420->X = (1E-6)*p1.getX(); // [m]
|
---|
| 130 | elementFP420->Y = (1E-6)*p1.getY(); // [m]
|
---|
| 131 | elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
|
---|
| 132 | elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
|
---|
| 133 | elementFP420->S = p1.getS(); // [m]
|
---|
| 134 | elementFP420->T = -1; // not yet implemented
|
---|
| 135 | elementFP420->E = p1.getE(); // not yet implemented
|
---|
| 136 | elementFP420->q2 = -1; // not yet implemented
|
---|
| 137 | elementFP420->side = sign(eta);
|
---|
| 138 | }
|
---|
| 139 |
|
---|
| 140 | }
|
---|
| 141 | // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
|
---|
| 142 | // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
|
---|
| 143 | } // if forward proton
|
---|
| 144 | }
|
---|
| 145 |
|
---|
| 146 | // Forward particles in CASTOR ?
|
---|
| 147 | // /* if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
|
---|
| 148 | // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
|
---|
| 149 | //
|
---|
| 150 | //
|
---|
| 151 | // } // CASTOR
|
---|
| 152 | // */
|
---|
| 153 | //
|
---|