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

 /*
 *   ---- Delphes ----
 *     A Fast Simulator for general purpose LHC detector
 *       S. Ovyn ~~~~ severine.ovyn@uclouvain.be
 *
 *     Center for Particle Physics and Phenomenology (CP3)
 *     Universite Catholique de Louvain (UCL)
 *     Louvain-la-Neuve, Belgium
 *             */

#include "VeryForward.h"
#include "H_RomanPot.h"
#include <iostream>
#include<cmath>

using namespace std;


//------------------------------------------------------------------------------

VeryForward::VeryForward() {
   DET = new RESOLution();
   init();
}

VeryForward::VeryForward(const string& DetDatacard) {
   DET = new RESOLution();
   DET->ReadDataCard(DetDatacard);
   init();
}

VeryForward::VeryForward(const RESOLution * DetDatacard) {
   DET = new RESOLution(*DetDatacard);
   init();
}

VeryForward::VeryForward(const VeryForward& vf) {
   DET = new RESOLution(*(vf.DET));
   beamline1 = new H_BeamLine(*(vf.beamline1));
   beamline2 = new H_BeamLine(*(vf.beamline2));
}

VeryForward& VeryForward::operator=(const VeryForward& vf){
   if (this==&vf) return *this;
   DET = new RESOLution(*(vf.DET));
   beamline1 = new H_BeamLine(*(vf.beamline1));
   beamline2 = new H_BeamLine(*(vf.beamline2));
   return *this;
}


void VeryForward::init() { 
  //Initialisation of Hector
  relative_energy = true; // should always be true
  kickers_on = 1;         // should always be 1
  
  // user should provide : (1) optics file for each beamline, and IPname,
  // and offset data (s,x) for optical elements
  beamline1 = new H_BeamLine(1,500.);
  beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
  beamline1->offsetElements(120,-0.097);
  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",220,2000); // RP 220m, 2mm, beam 1
  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",420,4000); // RP 420m, 4mm, beam 1
  beamline1->add(rp220_1);
  beamline1->add(rp420_1);
  
  beamline2 = new H_BeamLine(1,500.);
  beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
  beamline2->offsetElements(120,+0.097);
  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",220,2000);// RP 220m, 2mm, beam 2
  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",420,4000);// RP 420m, 4mm, beam 2
  beamline2->add(rp220_2);
  beamline2->add(rp420_2); 
  
}

void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
{
  int pid=abs(particle->PID);
  float eta=fabs(particle->Eta);
  
  
  TRootZdcHits *elementZdc;
  TLorentzVector genMomentum; 
  // Zero degree calorimeter, for forward neutrons and photons
  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->VFD_min_zdc ) {
    genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
    // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
    // !!!!!!!!! pas de smearing ! on garde trop d'info !
    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
    elementZdc->Set(genMomentum);
    
    // time of flight t is t = T + d/[ cos(theta) v ]
    //double tx = acos(particle->Px/particle->Pz);
    //double ty = acos(particle->Py/particle->Pz);
    //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
    //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-6));
    // assumes also that the emission angle is so small that 1/(cos theta) = 1
    elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    elementZdc->side = sign(eta);
    
    
  }
  
}
void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) 
{
  int pid=abs(particle->PID);
  float eta=fabs(particle->Eta);
  
  TRootRomanPotHits* elementRP220;
  TRootRomanPotHits* 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()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,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]
          elementRP220->T = -1;                 // not yet implemented
          elementRP220->E = p1.getE();  // not yet implemented
          elementRP220->q2 = -1;                // not yet implemented
          elementRP220->side = sign(eta);
          
        } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
          p1.propagate(DET->RP_420_s);
          elementFP420 = (TRootRomanPotHits*) 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]
          elementFP420->T = -1;                       // not yet implemented
          elementFP420->E = p1.getE();          // not yet implemented
          elementFP420->q2 = -1;                // not yet implemented
          elementFP420->side = sign(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
    //                                                         */
    //