/*********************************************************************** ** ** ** /----------------------------------------------\ ** ** | Delphes, a framework for the fast simulation | ** ** | of a generic collider experiment | ** ** \------------- arXiv:0903.2225v1 ------------/ ** ** ** ** ** ** This package uses: ** ** ------------------ ** ** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] ** ** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] ** ** FROG: [hep-ex/0901.2718v1] ** ** ** ** ------------------------------------------------------------------ ** ** ** ** Main authors: ** ** ------------- ** ** ** ** Severine Ovyn Xavier Rouby ** ** severine.ovyn@uclouvain.be xavier.rouby@cern ** ** ** ** Center for Particle Physics and Phenomenology (CP3) ** ** Universite catholique de Louvain (UCL) ** ** Louvain-la-Neuve, Belgium ** ** ** ** Copyright (C) 2008-2009, ** ** All rights reserved. ** ** ** ***********************************************************************/ #include "VeryForward.h" #include "H_RomanPot.h" #include #include using namespace std; //------------------------------------------------------------------------------ VeryForward::VeryForward() : DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)), relative_energy(true), // should always be true kickers_on(1) // should always be 1 { init(); //Initialisation of Hector } VeryForward::VeryForward(const string& DetDatacard) : DET(new RESOLution()) { DET->ReadDataCard(DetDatacard); const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s); beamline1 = new H_BeamLine(1,d_max); beamline2 = new H_BeamLine(1,d_max); init(); //Initialisation of Hector relative_energy = true; // should always be true kickers_on = 1; // should always be 1 } VeryForward::VeryForward(const RESOLution * DetDatacard) : DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)), relative_energy(true), // should always be true kickers_on(1) // should always be 1 { init(); //Initialisation of Hector } VeryForward::VeryForward(const VeryForward& vf) : DET(new RESOLution(*(vf.DET))), d_max(vf.d_max), beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))), relative_energy(vf.relative_energy), kickers_on(vf.kickers_on) { } VeryForward& VeryForward::operator=(const VeryForward& vf){ if (this==&vf) return *this; DET = new RESOLution(*(vf.DET)); d_max = vf.d_max; beamline1 = new H_BeamLine(*(vf.beamline1)); beamline2 = new H_BeamLine(*(vf.beamline2)); relative_energy =vf.relative_energy; kickers_on = vf.kickers_on; return *this; } void VeryForward::init() { //Initialisation of Hector relative_energy = true; // should always be true kickers_on = 1; // should always be 1 beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name); beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*(1E6)); // RP 220m, 2mm, beam 1 H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*(1E6)); // RP 420m, 4mm, beam 1 beamline1->add(rp220_1); beamline1->add(rp420_1); beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name); beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x); H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*(1E6));// RP 220m, 2mm, beam 2 H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*(1E6));// RP 420m, 4mm, beam 2 beamline2->add(rp220_2); beamline2->add(rp420_2); // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine // do not put explicit delete } void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle) { TRootZdcHits *elementZdc; float energy = particle->E; //TLorentzVector genMomentum; // Zero degree calorimeter, for forward neutrons and photons if (particle->Status ==1 && ( (particle->PID==pN && energy>DET->ZDC_n_E) || (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) ) && fabs(particle->Eta) > DET->VFD_min_zdc ) { elementZdc = (TRootZdcHits*) branchZDC->NewEntry(); // 1) energy smearing float energyS = -1.; if (particle->PID == pGAMMA) energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) + pow(DET->ELG_Czdc*particle->E,2) + pow(DET->ELG_Szdc*sqrt(particle->E),2) )); else // smearing with hadronic resolution energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) + pow(DET->HAD_Czdc*particle->E,2) + pow(DET->HAD_Szdc*sqrt(particle->E),2) )); elementZdc->E = energyS; // 2) time of flight t is t = T + d/[ cos(theta) v ] float cos_theta = 1; //very good approximation, if eta_zdc >3 if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation double tx = atan(particle->Px/particle->Pz); double ty = atan(particle->Py/particle->Pz); double theta = sqrt( pow(tx,2) + pow(ty,2) ); //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 // eta = 2.6 <-> cos(theta) = 0.99 // eta = 3.0 <-> cos(theta) = 0.995 cos_theta = cos(theta); } // units from StdHEP : Z [mm] T[mm/c] // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s] double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ; double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution); elementZdc->T = timeS; // 3) side: which ZDC has been hit? elementZdc->side = sign(particle->Eta); // 4) object nature : e.m. (photon) or had (neutron) ? elementZdc->hadronic_hit = (bool) (particle->PID==pN); } } void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) { int pid=particle->PID; TRootRomanPotHits* elementRP220; TRootForwardTaggerHits* elementFP420; TLorentzVector genMomentum; genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E); // if forward proton if( (pid == pP) && (particle->Status == 1) && (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) ) { // !!!!!!!! put here particle->CHARGE and particle->MASS H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS p1.smearAng(); p1.smearPos(); 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); p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E); H_BeamLine *beamline; if(genMomentum.Eta() >0) beamline = beamline1; else beamline = beamline2; p1.computePath(beamline,1); if(p1.stopped(beamline)) { if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") { p1.propagate(DET->RP_220_s); elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry(); elementRP220->X = (1E-6)*p1.getX(); // [m] elementRP220->Y = (1E-6)*p1.getY(); // [m] elementRP220->Tx = (1E-6)*p1.getTX(); // [rad] elementRP220->Ty = (1E-6)*p1.getTY(); // [rad] elementRP220->S = p1.getS(); // [m] /* time of flight t is t = T + d/[ cos(theta) v ] // nb: here we assume a straight path to the detector, which is not the case! // this time estimate is always underestimated (while exact for the ZDC case) float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation double tx = atan(particle->Px/particle->Pz); double ty = atan(particle->Py/particle->Pz); double theta = sqrt( pow(tx,2) + pow(ty,2) ); //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 // eta = 2.6 <-> cos(theta) = 0.99 // eta = 3.0 <-> cos(theta) = 0.995 cos_theta = cos(theta); } // units from StdHEP : Z [mm] T[mm/c] // units from Delphes : p1.getS [m] speed_of_light [m/s] //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ; //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] */ elementRP220->E = p1.getE(); // not yet implemented elementRP220->q2 = -1; // not yet implemented elementRP220->side = sign(particle->Eta); } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") { p1.propagate(DET->RP_420_s); elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry(); elementFP420->X = (1E-6)*p1.getX(); // [m] elementFP420->Y = (1E-6)*p1.getY(); // [m] elementFP420->Tx = (1E-6)*p1.getTX(); // [rad] elementFP420->Ty = (1E-6)*p1.getTY(); // [rad] elementFP420->S = p1.getS(); // [m] // time of flight t is t = T + d/[ cos(theta) v ] // nb: here we assume a straight path to the detector, which is not the case! // this time estimate is always underestimated (while exact for the ZDC case) float cos_theta = 1; //very good approximation, if CEN_max_calo_fwd >3 if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation double tx = atan(particle->Px/particle->Pz); double ty = atan(particle->Py/particle->Pz); double theta = sqrt( pow(tx,2) + pow(ty,2) ); //cout << "tx = " << tx << " ty = " << ty << " theta = " << theta << " cos(theta) = " << cos(theta) << endl; // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1 // eta = 2.6 <-> cos(theta) = 0.99 // eta = 3.0 <-> cos(theta) = 0.995 cos_theta = cos(theta); } // units from StdHEP : Z [mm] T[mm/c] // units from Delphes : p1.getS [m] speed_of_light [m/s] double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ; elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s] elementFP420->E = p1.getE(); // not yet implemented elementFP420->q2 = -1; // not yet implemented elementFP420->side = sign(particle->Eta); } } // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100)) // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl; } // if forward proton } // Forward particles in CASTOR ? // if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) // && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) { // // // } // CASTOR // */