ElasScatAnalysis.cc

Go to the documentation of this file.

// Analysis example for elastic scattering with PYTHIA8 (provided by Hasko Stenzel)

#include <iostream>
#include <stdio.h>

#include "HepMC/GenEvent.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/IO_AsciiParticles.h"
#include "HepMC/SimpleVector.h"
//#include "HepPDT/ParticleData.hh"
#include "CLHEP/Vector/LorentzVector.h"

// ROOT headers
#include "TH1.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include "TFile.h"

// fastjet headers
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/SISConePlugin.hh"

// header for this analysis class
#include "../include/ElasScatAnalysis.h"

using namespace std;

//**********************

/// empty default constructor
ElasScatAnalysis::ElasScatAnalysis()
{
}
/// empty default destructor
ElasScatAnalysis::~ElasScatAnalysis()
{
}


int ElasScatAnalysis::Init(double tr_max_eta, double tr_min_pt)
{
  // specify default eta cut
  m_max_eta=tr_max_eta;
  
  // specify default pt cut
  m_min_pt=tr_min_pt;
  
  // default Output file name
  m_outputFileName="ElasScatAnalysis.root";
  m_outputRootDir="ElasScat";
 
  n_proton   = initHist("nproton", "nproton", "Number of Protons", 5, -0.5, 4.5);
  // some histograms for the outgoing protons
  p1_pt =  baseAnalysis::initHist( string("Ptp1"), string("p_{T} proton1"), string(" p_{T} [GeV]"), 100,0.0, 1.);
  p1_eta=baseAnalysis::initHist(string("eta1"),string("#eta 1st proton"),string("#eta"), 200, 9., 16.);
  p1_theta=baseAnalysis::initHist(string("theta1"),string("sacttering angle p1"),string("#theta [#mu rad]"), 200, 0.,100.);
  p1_phi=baseAnalysis::initHist(string("p1_phi"),string("#varphi p1"),string("#varphi"), 64, -3.2, 3.2);

  p2_pt =  baseAnalysis::initHist( string("Ptp2"), string("p_{T} proton2"), string(" p_{T} [GeV]"), 100, 0.0, 1.);
  p2_eta=baseAnalysis::initHist(string("eta2"),string("#eta 2nd proton"),string("#eta"), 200, -16.,-9.);
  p2_theta=baseAnalysis::initHist(string("theta2"),string("sacttering angle p2"),string("#theta [#mu rad]"), 200, -100.,0.);
  p2_phi=baseAnalysis::initHist(string("p2_phi"),string("#varphi p2"),string("#varphi"), 64, -3.2, 3.2);

  // two histograms for the t-spectrum

  t_spect =  baseAnalysis::initHist(string("tspect"), string("t-spectrum"), string("t"), 100, 0.0,0.1);
  tau_spect =  baseAnalysis::initHist(string("tauspect"), string("#tau-spectrum"), string("#tau=log_{10}|t|"), 100, -5.0,0.0);

  // beam proton histograms (to check beam class settings)
  pi1_en = baseAnalysis::initHist( string("pi1en"), string("Energy deviation"), string(" E-E_{0} [GeV]"), 100, -2.0, 2.);
  pi1_divx = baseAnalysis::initHist( string("pi1divx"), string("Divergence in x beam 1"), string("angle [#mu rad]"), 100, -2.0, 2.);
  pi1_divy = baseAnalysis::initHist( string("pi1divy"), string("Divergence in y beam 1"), string("angle [#mu rad]"), 100, -2.0, 2.);

  pi2_en = baseAnalysis::initHist( string("pi2en"), string("Energy deviation"), string(" E-E_{0} [GeV]"), 100, -2.0, 2.);
  pi2_divx = baseAnalysis::initHist( string("pi2divx"), string("Divergence in x beam 2"), string("angle [#mu rad]"), 100, -2.0, 2.);
  pi2_divy = baseAnalysis::initHist( string("pi2divy"), string("Divergence in y beam 2"), string("angle [#mu rad]"), 100, -2.0, 2.);


  return true; 
}


int ElasScatAnalysis::Process(HepMC::GenEvent* hepmcevt){ 

//   // print number of events
//   if (hepmcevt->event_number()%1 == 0) {
  //    std::cout << "Processing event " << hepmcevt->event_number() << std::endl;
  //   }
  
 
// elastic analysis
  int pcount=0;
  int pint=0;
  double const plhc=7000.;
  
  HepMC::GenParticle* pi1=0;
  HepMC::GenParticle* pi2=0;
  HepMC::GenParticle* p1=0;
  HepMC::GenParticle* p2=0;

  // loop over the final state particles 
  for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
    {
      int pid = (*p)->pdg_id();
      //      (*p)->print();
      //if( ((*p)->status() == -12))  {
      if( (hepmcevt->beam_particles().first->pdg_id()==2212 && !(*p)->production_vertex())
          || (hepmcevt->beam_particles().second->pdg_id()==2212 && !(*p)->production_vertex()) ) {
        if( pid == 2212){ 
          pint++;
          //(*p)->print();
          if(pint == 1){pi1=(*p);}
          if(pint == 2){pi2=(*p);}
        }
      }

      if( ((*p)->status()%1000 == 1) && (!(*p)->end_vertex()) )  {
        //if( (!(*p)->end_vertex()) )  {
        if( pid == 2212){ 
          pcount++;
          if(pcount == 1){p1=(*p);}
          if(pcount == 2){p2=(*p);}
          if(pcount > 2){cout << "Strange: More than two protons in this event!" << "\n";}
        }
      }
    }
    if(pint == 2){
      pi1_en->Fill((pi1->momentum().e())-plhc);
      double divx = (pi1->momentum().px())/(pi1->momentum().pz());
      double divy = (pi1->momentum().py())/(pi1->momentum().pz());
      pi1_divx->Fill(divx*1.E6);
      pi1_divy->Fill(divy*1.E6);

      pi2_en->Fill((pi2->momentum().e())-plhc);
      divx = (pi2->momentum().px())/(pi2->momentum().pz());
      divy = (pi2->momentum().py())/(pi2->momentum().pz());
      pi2_divx->Fill(divx*1.E6);
      pi2_divy->Fill(divy*1.E6);
    }
    
    if(pcount==2){
     n_proton->Fill(pcount);
     p1_pt->Fill(p1->momentum().perp());
     p1_eta->Fill(p1->momentum().eta());
     p1_phi->Fill(p1->momentum().phi());
     p1_theta->Fill((p1->momentum().theta())*1.E6); //micro-radians
     p2_pt->Fill(p2->momentum().perp());
     p2_eta->Fill(p2->momentum().eta());
     p2_phi->Fill(p2->momentum().phi());
     p2_theta->Fill((p2->momentum().theta()-M_PI)*1.E6); //micro-radians

     // t-calculation 
     double st1 =(p1->momentum().theta())*plhc;
     double st2 =(p2->momentum().theta()-M_PI)*plhc;
     double t1 = st1*st1; 
     double t2 = st2*st2; 
     double tmean = (t1+t2)/2.; // average t-value reconstructed from the scattering angles

     double tau = log10(tmean);
     t_spect->Fill(tmean); 
     tau_spect->Fill(tau); 
    }
  
  return true;

} //end of Process

---