Ticket #406: TrackVertexSmearing.cc

/** \class TrackVertexSmearing
 *
 *  Performs track vertex smearing
 *
 *  $Date: 2014-16-03 14:57:44 +0100   
 *
 *
 *  \author F. Nesti,   IRB, Zagreb
 *
 */
 

#include "modules/TrackVertexSmearing.h"

#include "classes/DelphesClasses.h"
#include "classes/DelphesFactory.h"
#include "classes/DelphesFormula.h"

#include "ExRootAnalysis/ExRootResult.h"
#include "ExRootAnalysis/ExRootFilter.h"
#include "ExRootAnalysis/ExRootClassifier.h"

#include "TMath.h"
#include "TString.h"
#include "TFormula.h"
#include "TRandom3.h"
#include "TObjArray.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"

#include <algorithm> 
#include <stdexcept>
#include <iostream>
#include <sstream>

using namespace std;

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

TrackVertexSmearing::TrackVertexSmearing() :
  fFormula(0), fItInputArray(0)
{
  fFormula = new DelphesFormula;
}

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

TrackVertexSmearing::~TrackVertexSmearing()
{
  if(fFormula) delete fFormula;
}

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

void TrackVertexSmearing::Init()
{
  // read resolution formula

  fFormula->Compile(GetString("ResolutionFormula", "0.0"));

  // import input array

  fInputArray = ImportArray(GetString("InputArray", "TrackMerger/tracks"));
  fItInputArray = fInputArray->MakeIterator();

  // create output array

  fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
}

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

void TrackVertexSmearing::Finish()
{
  if(fItInputArray) delete fItInputArray;
}

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

void TrackVertexSmearing::Process()
{
  Candidate *candidate, *particle, *mother;
  Double_t x, y, z, t, sx, sy, sz;
  Double_t pt, eta;

  cout<<"New event"<<endl;


  fItInputArray->Reset();
  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
  {
  
    //take momentum 
    particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
  
    const TLorentzVector &candidateMomentum = particle->Momentum;
    //  const TLorentzVector &candidateMomentum = candidate->Momentum;
    const TLorentzVector &candidatePosition = particle->Position;
    
    eta = candidateMomentum.Eta();
    pt = candidateMomentum.Pt();
      
    // vertex before smearing
    x =  candidatePosition.X();
    y =  candidatePosition.Y();
    z =  candidatePosition.Z();
    t =  candidatePosition.T();
   
    // calculate smeared values   
    sx = gRandom->Gaus(0,fFormula->Eval(pt, eta));
    sy = gRandom->Gaus(0,fFormula->Eval(pt, eta));
    sz = gRandom->Gaus(0,fFormula->Eval(pt, eta)); 
       
    fprintf(stdout, "Smearing (%3d) (pt,eta)=(%6.3f, %6.3f), formula=%6.3f  by dx,dy,dz -> %6.4f  %6.4f  %6.4f\n", particle->PID, pt, eta, fFormula->Eval(pt, eta), sx, sy, sz);

    x += sx;
    y += sy;
    z += sz;
       

    //fill smeared values in candidate
  /* mother = candidate;
    
    candidate = static_cast<Candidate*>(candidate->Clone());
       
    candidate->Position.SetXYZT(x, y, z, t);
    
        cout << z  << endl ;
        cout << candidate->Position.Z()  << endl ;

    candidate->AddCandidate(mother);
    fOutputArray->Add(candidate);
 */
 
    
    candidate = static_cast<Candidate*>(candidate->Clone());
    mother = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
    mother = static_cast<Candidate*>(mother->Clone());
       
    mother->Position.SetXYZT(x, y, z, t);
    
        cout << z  << endl ;
        cout << mother->Position.Z()  << endl ;

    candidate->GetCandidates()->AddAt(mother, 0);
    
    fOutputArray->Add(candidate);
 
 
  }
}

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