Changeset 405 in svn for trunk/src/VeryForward.cc
- Timestamp:
- May 20, 2009, 10:39:01 AM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/VeryForward.cc
r389 r405 31 31 32 32 #include "VeryForward.h" 33 #include "PdgParticle.h" 33 34 #include "H_RomanPot.h" 34 #include "PdgParticle.h" 35 35 36 #include <iostream> 37 #include <fstream> 36 38 #include<cmath> 37 39 38 40 using namespace std; 39 41 42 /* Notes on the correct initialisation for Hector 43 * -- these notes apply to the LHC beamlines 44 * 45 * beam1 : forward direction is for increasing 's' values 46 * beamline1 = new H_BeamLine(1,...); 47 * beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name); 48 * beam2 : forward direction is for decreasing 's' values 49 * beamline2 = new H_BeamLine(-1,...); 50 * beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name); 51 * 52 * relative_energy should be false -- kickers_on should be 1 53 * 54 */ 40 55 41 56 //------------------------------------------------------------------------------ 42 57 VeryForward::VeryForward() : 43 58 DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), 44 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine( 1,d_max)),45 rel ative_energy(true), // should always be true46 kickers _on(1) // should always be 159 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)), 60 rel_energy(true), // should always be true 61 kickers(1) // should always be 1 47 62 { 48 63 init(); //Initialisation of Hector … … 55 70 const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s); 56 71 beamline1 = new H_BeamLine(1,d_max); 57 beamline2 = new H_BeamLine( 1,d_max);72 beamline2 = new H_BeamLine(-1,d_max); 58 73 init(); //Initialisation of Hector 59 rel ative_energy = true; // should always be true60 kickers _on= 1; // should always be 174 rel_energy = true; // should always be true 75 kickers = 1; // should always be 1 61 76 } 62 77 63 78 VeryForward::VeryForward(const RESOLution * DetDatacard) : 64 79 DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), 65 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine( 1,d_max)),66 rel ative_energy(true), // should always be true67 kickers _on(1) // should always be 180 beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)), 81 rel_energy(true), // should always be true 82 kickers(1) // should always be 1 68 83 { 69 84 init(); //Initialisation of Hector … … 73 88 DET(new RESOLution(*(vf.DET))), d_max(vf.d_max), 74 89 beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))), 75 rel ative_energy(vf.relative_energy),76 kickers _on(vf.kickers_on) {90 rel_energy(vf.rel_energy), 91 kickers(vf.kickers) { 77 92 } 78 93 … … 83 98 beamline1 = new H_BeamLine(*(vf.beamline1)); 84 99 beamline2 = new H_BeamLine(*(vf.beamline2)); 85 rel ative_energy =vf.relative_energy;86 kickers _on = vf.kickers_on;100 rel_energy =vf.rel_energy; 101 kickers = vf.kickers; 87 102 return *this; 88 103 } … … 93 108 static unsigned int counter; 94 109 counter =0; 95 relative_energy = true; // should always be true 96 kickers_on = 1; // should always be 1 110 extern bool relative_energy; 111 extern int kickers_on; 112 relative_energy = rel_energy; // should always be true 113 kickers_on = kickers; // should always be 1 97 114 beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name); 98 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); 99 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); // RP 220m, 2mm, beam 1 100 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); // RP 420m, 4mm, beam 1 115 beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything 116 H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); 117 // RP 220m, 2mm, beam 1 118 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); 119 // RP 420m, 4mm, beam 1 120 //rp220_1->printProperties(); 121 //rp420_1->printProperties(); 101 122 beamline1->add(rp220_1); 102 123 beamline1->add(rp420_1); 103 124 104 125 beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name); 105 beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x); 106 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);// RP 220m, 2mm, beam 2 107 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);// RP 420m, 4mm, beam 2 126 beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything 127 H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6); 128 // RP 220m, 2mm, beam 2 129 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6); 130 // RP 420m, 4mm, beam 2 131 //rp220_2->printProperties(); 132 //rp420_2->printProperties(); 108 133 beamline2->add(rp220_2); 109 134 beamline2->add(rp420_2); 110 135 // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine 111 136 // do not put explicit delete 137 } 138 139 140 float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) { 141 // time of flight t is t = T + d/[ cos(theta) v ] 142 float cos_theta = 1; //very good approximation, if detector_etamin >3 143 if (detector_etamin<3) { // if smaller eta -> make the complete calculation 144 double tx = atan(particle->Px/particle->Pz); 145 double ty = atan(particle->Py/particle->Pz); 146 double theta = sqrt( pow(tx,2) + pow(ty,2) ); 147 // cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; 148 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 149 // eta = 2.6 <-> cos(theta) = 0.99 150 // eta = 3.0 <-> cos(theta) = 0.995 151 cos_theta = cos(theta); 152 } 153 // units from StdHEP : Z [mm] T[mm/c] 154 // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s] 155 double flight_distance = (detector_s - particle->Z*(1E-3))/cos_theta ; 156 // assumes highly relativistic particles 157 double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; 158 double timeS = gRandom->Gaus(flight_time,detector_t_resolution); 159 return timeS; 112 160 } 113 161 … … 125 173 elementZdc = (TRootZdcHits*) branchZDC->NewEntry(); 126 174 127 128 elementZdc->pid = particle->PID;129 130 131 // for compatibility with 'old' version132 175 TLorentzVector genMomentum; 133 176 genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E); 134 elementZdc->Set(genMomentum); 135 // ****************** 136 137 138 //particle->print(); 177 elementZdc->Set(genMomentum); // initialises the gen-level data 178 elementZdc->pid = particle->PID; 139 179 140 180 // 1) energy smearing … … 149 189 pow(DET->HAD_Szdc*sqrt(particle->E),2) )); 150 190 elementZdc->E = energyS; 151 152 153 // 2) time of flight t is t = T + d/[ cos(theta) v ] 154 float cos_theta = 1; //very good approximation, if eta_zdc >3 155 if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation 156 double tx = atan(particle->Px/particle->Pz); 157 double ty = atan(particle->Py/particle->Pz); 158 double theta = sqrt( pow(tx,2) + pow(ty,2) ); 159 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; 160 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 161 // eta = 2.6 <-> cos(theta) = 0.99 162 // eta = 3.0 <-> cos(theta) = 0.995 163 cos_theta = cos(theta); 164 } 165 // units from StdHEP : Z [mm] T[mm/c] 166 // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s] 167 double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ; 168 double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] 169 double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution); 170 elementZdc->T = timeS; 191 192 // 2) time of flight t is t = T + d/[ cos(theta) v ] + detector smearing on time 193 elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution); 171 194 172 195 // 3) side: which ZDC has been hit? … … 174 197 175 198 // 4) object nature : e.m. (photon) or had (neutron) ? 176 //elementZdc->hadronic_hit = (bool) (particle->PID==pN); 177 } 199 elementZdc->hadronic_hit = (bool) (particle->PID!=pGAMMA); 200 } // if neutrons or photons over E_threshold 201 178 202 } 179 203 … … 181 205 void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) 182 206 { 207 if(particle->Status != 1) return; // reject particles that are not final ones 208 extern bool relative_energy; 209 relative_energy = rel_energy; 210 extern int kickers_on; 211 kickers_on = kickers; 212 183 213 float charge = particle->Charge, mass = particle->M; 214 //float charge, mass, ctau; 215 //charge = mass = ctau = UNDEFINED; 184 216 if (mass<-999) { // unitialised! 185 217 PdgParticle pdg_part = DET->PDGtable[particle->PID]; 186 charge = pdg_part.charge(); 187 mass = pdg_part.mass(); 218 charge = pdg_part.charge(); // e+ 219 mass = pdg_part.mass(); // GeV 220 // ctau = pdg_part.ctau(); // m 221 // cout << "ctau = " << ctau << endl; 188 222 } 189 //if(particle->Charge!=1) return; // only particles with Q=+1 can hope to reach RP200/FP420 190 223 224 225 if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420 226 //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl; 191 227 TRootRomanPotHits* elementRP220; 192 228 //TRootForwardTaggerHits* elementFP420; … … 197 233 198 234 // to go faster, why not rejecting particles already going into the ZDC? 199 if( particle->PID == pP) 200 if( (particle->Status == 1) && (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) ) 235 236 // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ; 237 // so if RP/FP too far away, the particle must be a proton or a mu+ 238 if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return; 239 240 if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd ) 201 241 { 202 //cout << "VeryForward :: M = " << mass << "\t Q = " << charge << "\t\t " << particle->PID << endl;203 242 H_BeamParticle p1(mass,charge); 204 p1.smearAng(); p1.smearPos(); // vertex smearing 205 cout << "x = " << p1.getX() + DET->RP_cross_x 206 << " y= " << p1.getY() + DET->RP_cross_y 207 << " tx= " << p1.getTX() - kickers_on*DET->RP_cross_ang 208 << " ty=" << p1.getTY() << endl; 209 p1.setPosition(p1.getX()+DET->RP_cross_x,p1.getY()+DET->RP_cross_y,p1.getTX()-kickers_on*DET->RP_cross_ang,p1.getTY(),0); 210 //p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E); 243 double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad 244 double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad 245 //p1.smearAng(); p1.smearPos(); // vertex smearing do not put it here !!! 246 247 /* cout << "x = " << particle->X << " + " << p1.getX() << " + " << DET->RP_cross_x 248 << " y= " << particle->Y << " + " << p1.getY() << " + " << DET->RP_cross_y 249 << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x 250 << " ty= " << ty << " + " << p1.getTY() << " - " << kickers_on*DET->RP_cross_ang_y 251 << " z= " << particle->Z << endl;*/ 252 253 // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done 254 // all in micrometers or microradians 255 p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x, 256 (1E3)*particle->Y + p1.getY() + DET->RP_cross_y, 257 tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x, 258 ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y, 259 (1E3)*particle->Z); 211 260 p1.setE(particle->E); 212 261 … … 218 267 219 268 if(p1.stopped(beamline)) { 269 // roman pots at 220 m 220 270 if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") { 221 static unsigned int counter; 222 counter++;223 if (counter==1) {224 //p1.getPath(0,"p1path.txt");225 cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge226 << "\t" << mass << "=" << particle->M227 << "\t E=" << particle->E 228 << endl; 229 } 271 272 /*static unsigned int counter; 273 counter++; 274 if (counter==1) { 275 p1.getPath(0,"p1path.txt"); 276 cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge 277 << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl; 278 }*/ 279 230 280 p1.propagate(DET->RP_220_s); 231 281 elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry(); 282 283 // detector measurements 232 284 elementRP220->X = (1E-6)*p1.getX(); // [m] 233 285 elementRP220->Y = (1E-6)*p1.getY(); // [m] … … 235 287 elementRP220->Ty = (1E-6)*p1.getTY(); // [rad] 236 288 elementRP220->S = p1.getS(); // [m] 237 238 /* time of flight t is t = T + d/[ cos(theta) v ] 239 // nb: here we assume a straight path to the detector, which is not the case! 240 // this time estimate is always underestimated (while exact for the ZDC case) 241 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3 242 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation 243 double tx = atan(particle->Px/particle->Pz); 244 double ty = atan(particle->Py/particle->Pz); 245 double theta = sqrt( pow(tx,2) + pow(ty,2) ); 246 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; 247 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 248 // eta = 2.6 <-> cos(theta) = 0.99 249 // eta = 3.0 <-> cos(theta) = 0.995 250 cos_theta = cos(theta); 251 } 252 // units from StdHEP : Z [mm] T[mm/c] 253 // units from Delphes : p1.getS [m] speed_of_light [m/s] 254 //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ; 255 //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] 256 */ 257 elementRP220->E = p1.getE(); // not yet implemented 258 elementRP220->q2 = -1; // not yet implemented 259 elementRP220->side = sign(particle->Eta); 260 261 262 elementRP220->pid = particle->PID; 289 elementRP220->T = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution); 290 elementRP220->side = sign(particle->Eta); 263 291 264 } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") { 292 // reconstructed data 293 float sE = p1.getE(); // apply the smearing here!!! 294 elementRP220->E = sE; // not yet implemented 295 elementRP220->q2 = UNDEFINED; // not yet implemented 296 297 // generator level data 298 elementRP220->pid = particle->PID; 299 elementRP220->Set(genMomentum); 300 } // if RP220 301 302 // proton taggers at 420 m 303 else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") { 304 265 305 p1.propagate(DET->RP_420_s); 306 elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry(); 266 307 //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry(); 267 elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry(); 308 309 // detector measurements 268 310 elementFP420->X = (1E-6)*p1.getX(); // [m] 269 311 elementFP420->Y = (1E-6)*p1.getY(); // [m] … … 271 313 elementFP420->Ty = (1E-6)*p1.getTY(); // [rad] 272 314 elementFP420->S = p1.getS(); // [m] 273 274 // time of flight t is t = T + d/[ cos(theta) v ] 275 // nb: here we assume a straight path to the detector, which is not the case! 276 // this time estimate is always underestimated (while exact for the ZDC case) 277 float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3 278 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation 279 double tx = atan(particle->Px/particle->Pz); 280 double ty = atan(particle->Py/particle->Pz); 281 double theta = sqrt( pow(tx,2) + pow(ty,2) ); 282 //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; 283 // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 284 // eta = 2.6 <-> cos(theta) = 0.99 285 // eta = 3.0 <-> cos(theta) = 0.995 286 cos_theta = cos(theta); 287 } 288 // units from StdHEP : Z [mm] T[mm/c] 289 // units from Delphes : p1.getS [m] speed_of_light [m/s] 290 double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ; 291 elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] 315 elementFP420->T = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution); 316 elementFP420->side = sign(particle->Eta); 317 318 // reconstructed data 292 319 elementFP420->E = p1.getE(); // not yet implemented 293 elementFP420->q2 = -1; // not yet implemented 294 elementFP420->side = sign(particle->Eta); 295 296 297 elementFP420->pid = particle->PID; 298 } 299 300 } 301 // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100)) 302 // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl; 320 elementFP420->q2 = UNDEFINED; // not yet implemented 321 322 // generator level data 323 elementFP420->pid = particle->PID; 324 elementFP420->Set(genMomentum); 325 } // if FP420 326 327 } // if stopped 303 328 } // if forward proton 304 } 305 306 // Forward particles in CASTOR ? 307 // if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 308 // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) { 309 // 310 // 311 // } // CASTOR 312 // */ 313 329 330 }
Note:
See TracChangeset
for help on using the changeset viewer.