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

/***********************************************************************
**                                                                    **
**          /----------------------------------------------\          **
**         |  Delphes, a framework for the fast simulation  |         **
**         |       of a  generic collider experiment        |         **
**          \-------------  arXiv:0903.2225v1  ------------/          **
**                                                                    **
**                                                                    **
**   This package uses:                                               **
**   ------------------                                               **
**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
**   FROG: [hep-ex/0901.2718v1]                                       **
**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
**                                                                    **
** ------------------------------------------------------------------ **
**                                                                    **
**   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 "PdgParticle.h"
#include "H_RomanPot.h"

#include <iostream>
#include <fstream>
#include<cmath>

using namespace std;

/* Notes on the correct initialisation for Hector
 * -- these notes apply to the LHC beamlines
 *
 *  beam1 : forward direction is for increasing 's' values
 *    beamline1 = new H_BeamLine(1,...);
 *    beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
 *  beam2 : forward direction is for decreasing 's' values
 *    beamline2 = new H_BeamLine(-1,...);
 *    beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
 *
 *  relative_energy should be false -- kickers_on should be 1
 *
 */

//------------------------------------------------------------------------------
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)),
   rel_energy(true), // should always be true
   kickers(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
   rel_energy = true; // should always be true
   kickers = 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)),
        rel_energy(true), // should always be true
        kickers(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))),
   rel_energy(vf.rel_energy),
   kickers(vf.kickers) {
}

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));
   rel_energy =vf.rel_energy;
   kickers = vf.kickers;
   return *this;
}


void VeryForward::init() { 
  //Initialisation of Hector
  static unsigned int counter;
  counter =0;
  extern bool relative_energy;
  extern int kickers_on;
  relative_energy = rel_energy; // should always be true
  kickers_on = kickers;         // should always be 1
  beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                                
  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
  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  
  //rp220_1->printProperties();
  //rp420_1->printProperties();
  beamline1->add(rp220_1); 
  beamline1->add(rp420_1);
  //beamline1->alignElement("\"MQXA.1R5\"",+0.0005,0);
  //beamline1->alignElement("\"MQM.9R5.B1\"",-0.0001,0);
  //beamline1->alignElement("\"MQML.5R5.B1\"",+0.0001,0);
  
  beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);                              
  beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
  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     
  //rp220_2->printProperties();
  //rp420_2->printProperties();
  beamline2->add(rp220_2);
  beamline2->add(rp420_2); 
  //beamline2->alignElement("\"MQXA.1L5\"",+0.0005,0);
  // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
  // do not put explicit delete
}


float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) {
    // time of flight t is t = T + d/[ cos(theta) v ]
    float cos_theta = 1;     //very good approximation, if detector_etamin >3
    if (detector_etamin<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 = (detector_s - particle->Z*(1E-3))/cos_theta ;
    // assumes highly relativistic particles
    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; 
    double timeS = gRandom->Gaus(flight_time,detector_t_resolution);
    return timeS;
}



void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
{
  TRootZdcHits *elementZdc;
  float energy = particle->E;

  // 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();

    TLorentzVector genMomentum;
    genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    elementZdc->Set(genMomentum);  // initialises the gen-level data
    elementZdc->pid = particle->PID;
 
    // 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 ] + detector smearing on time
    elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution);

    // 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!=pGAMMA);
  } // if neutrons or photons over E_threshold

}


void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) 
{
  if(particle->Status != 1) return; // reject particles that are not final ones
  extern bool relative_energy;
  relative_energy = rel_energy;
  extern int kickers_on;
  kickers_on = kickers;

  float charge = particle->Charge, mass = particle->M;
  //float charge, mass, ctau;
  //charge = mass = ctau = UNDEFINED;
  if (mass<-999) { // unitialised!
          PdgParticle pdg_part = DET->PDGtable[particle->PID];
          charge  = pdg_part.charge();  // e+
          mass    = pdg_part.mass();    // GeV
  //        ctau    = pdg_part.ctau();    // m
  // cout << "ctau = " << ctau << endl;
  }


  if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420
  //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl;
  TRootRomanPotHits* elementRP220;
  //TRootForwardTaggerHits* elementFP420;
  TRootRomanPotHits* elementFP420;
  
  TLorentzVector genMomentum;
  genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
  
  // to go faster, why not rejecting particles already going into the ZDC?
  
  // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ; 
  // so if RP/FP too far away, the particle must be a proton or a mu+
  if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return;
  if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd )
    {
      H_BeamParticle p1(mass,charge); 
      double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad
      double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad
      //p1.smearAng();   p1.smearPos(); // vertex smearing do not put it here !!!

      /* cout << "x = "  << particle->X << " + " << p1.getX()   << " + " <<  DET->RP_cross_x 
             << " y= "  << particle->Y << " + " << p1.getY()   << " + " << DET->RP_cross_y 
             << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x 
             << " ty= " << ty << " + " << p1.getTY()  << " - " << kickers_on*DET->RP_cross_ang_y  
             << " z= "  << particle->Z << endl;*/

      // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done
      // all in micrometers or microradians

      // for checking purposes
      /*p1.smearAng();   p1.smearPos(); tx=0; ty=0;      // to be removed !!!! test only !!!!
      ofstream xdist("xdist",ios::ios_base::app);
      xdist << endl;
      xdist << (1E3)*particle->X + p1.getX() + DET->RP_cross_x << " = " << (1E3)*particle->X << " + " << p1.getX() << " + " << DET->RP_cross_x << endl;
      xdist << (1E3)*particle->Y + p1.getY() + DET->RP_cross_y << " = " << (1E3)*particle->Y << " + " << p1.getY() << " + " << DET->RP_cross_y << endl;
      xdist << tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x << " = " << tx << " + " << p1.getTX() << " - " << kickers_on << "*" << DET->RP_cross_ang_x<< endl;
      xdist << ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y << " = " << ty << " + " << p1.getTY() << " - " << kickers_on << "*" << DET->RP_cross_ang_y<< endl;
      xdist.close();
      */

      
      p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x, 
                     (1E3)*particle->Y + p1.getY() + DET->RP_cross_y, 
                     tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x,
                     ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y,
                     0*(1E3)*particle->Z);
      p1.setE(particle->E);
       
      H_BeamLine *beamline;
      if(genMomentum.Eta() >0) beamline = beamline1;
      else beamline = beamline2;
      
      p1.computePath(beamline,1);
      
      if(p1.stopped(beamline)) {
        // roman pots at 220 m
        if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {

        /*static unsigned int counter;
        counter++;
        if (counter==1) {
           p1.getPath(0,"p1path.txt");
           cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge 
                << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl;
          }*/

          p1.propagate(DET->RP_220_s);
          elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
         
          // detector measurements
          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 = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution);
          elementRP220->side = sign(particle->Eta);
          
          // reconstructed data  
          float sE = p1.getE(); // apply the smearing here!!!
          elementRP220->E = sE;           // not yet implemented
          elementRP220->q2 = UNDEFINED;   // not yet implemented

          // generator level data
          elementRP220->pid = particle->PID;
          elementRP220->Set(genMomentum);
        } // if RP220

        // proton taggers at 420 m
        else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {

          p1.propagate(DET->RP_420_s);
          elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
          //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();

          // detector measurements
          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 = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution);
          elementFP420->side = sign(particle->Eta);

          // reconstructed data
          elementFP420->E = p1.getE();          // not yet implemented
          elementFP420->q2 = UNDEFINED;         // not yet implemented

          // generator level data
          elementFP420->pid = particle->PID;
          elementFP420->Set(genMomentum);
        }  // if FP420

      } // if stopped
    } // if forward proton

}