ElasScatAnalysis.cc

Go to the documentation of this file.
00001 // Analysis example for elastic scattering with PYTHIA8 (provided by Hasko Stenzel)
00002 
00003 #include <iostream>
00004 #include <stdio.h>
00005 
00006 #include "HepMC/GenEvent.h"
00007 #include "HepMC/IO_GenEvent.h"
00008 #include "HepMC/GenParticle.h"
00009 #include "HepMC/GenVertex.h"
00010 #include "HepMC/IO_AsciiParticles.h"
00011 #include "HepMC/SimpleVector.h"
00012 //#include "HepPDT/ParticleData.hh"
00013 #include "CLHEP/Vector/LorentzVector.h"
00014 
00015 // ROOT headers
00016 #include "TH1.h"
00017 #include "TMath.h"
00018 #include "TLorentzVector.h"
00019 #include "TFile.h"
00020 
00021 // fastjet headers
00022 #include "fastjet/PseudoJet.hh"
00023 #include "fastjet/ClusterSequence.hh"
00024 #include "fastjet/JetDefinition.hh"
00025 #include "fastjet/SISConePlugin.hh"
00026 
00027 // header for this analysis class
00028 #include "../include/ElasScatAnalysis.h"
00029 
00030 using namespace std;
00031 
00032 //**********************
00033 
00034 /// empty default constructor
00035 ElasScatAnalysis::ElasScatAnalysis()
00036 {
00037 }
00038 /// empty default destructor
00039 ElasScatAnalysis::~ElasScatAnalysis()
00040 {
00041 }
00042 
00043 
00044 int ElasScatAnalysis::Init(double tr_max_eta, double tr_min_pt)
00045 {
00046   // specify default eta cut
00047   m_max_eta=tr_max_eta;
00048   
00049   // specify default pt cut
00050   m_min_pt=tr_min_pt;
00051   
00052   // default Output file name
00053   m_outputFileName="ElasScatAnalysis.root";
00054   m_outputRootDir="ElasScat";
00055  
00056   n_proton   = initHist("nproton", "nproton", "Number of Protons", 5, -0.5, 4.5);
00057   // some histograms for the outgoing protons
00058   p1_pt =  baseAnalysis::initHist( string("Ptp1"), string("p_{T} proton1"), string(" p_{T} [GeV]"), 100,0.0, 1.);
00059   p1_eta=baseAnalysis::initHist(string("eta1"),string("#eta 1st proton"),string("#eta"), 200, 9., 16.);
00060   p1_theta=baseAnalysis::initHist(string("theta1"),string("sacttering angle p1"),string("#theta [#mu rad]"), 200, 0.,100.);
00061   p1_phi=baseAnalysis::initHist(string("p1_phi"),string("#varphi p1"),string("#varphi"), 64, -3.2, 3.2);
00062 
00063   p2_pt =  baseAnalysis::initHist( string("Ptp2"), string("p_{T} proton2"), string(" p_{T} [GeV]"), 100, 0.0, 1.);
00064   p2_eta=baseAnalysis::initHist(string("eta2"),string("#eta 2nd proton"),string("#eta"), 200, -16.,-9.);
00065   p2_theta=baseAnalysis::initHist(string("theta2"),string("sacttering angle p2"),string("#theta [#mu rad]"), 200, -100.,0.);
00066   p2_phi=baseAnalysis::initHist(string("p2_phi"),string("#varphi p2"),string("#varphi"), 64, -3.2, 3.2);
00067 
00068   // two histograms for the t-spectrum
00069 
00070   t_spect =  baseAnalysis::initHist(string("tspect"), string("t-spectrum"), string("t"), 100, 0.0,0.1);
00071   tau_spect =  baseAnalysis::initHist(string("tauspect"), string("#tau-spectrum"), string("#tau=log_{10}|t|"), 100, -5.0,0.0);
00072 
00073   // beam proton histograms (to check beam class settings)
00074   pi1_en = baseAnalysis::initHist( string("pi1en"), string("Energy deviation"), string(" E-E_{0} [GeV]"), 100, -2.0, 2.);
00075   pi1_divx = baseAnalysis::initHist( string("pi1divx"), string("Divergence in x beam 1"), string("angle [#mu rad]"), 100, -2.0, 2.);
00076   pi1_divy = baseAnalysis::initHist( string("pi1divy"), string("Divergence in y beam 1"), string("angle [#mu rad]"), 100, -2.0, 2.);
00077 
00078   pi2_en = baseAnalysis::initHist( string("pi2en"), string("Energy deviation"), string(" E-E_{0} [GeV]"), 100, -2.0, 2.);
00079   pi2_divx = baseAnalysis::initHist( string("pi2divx"), string("Divergence in x beam 2"), string("angle [#mu rad]"), 100, -2.0, 2.);
00080   pi2_divy = baseAnalysis::initHist( string("pi2divy"), string("Divergence in y beam 2"), string("angle [#mu rad]"), 100, -2.0, 2.);
00081 
00082 
00083   return true; 
00084 }
00085 
00086 
00087 int ElasScatAnalysis::Process(HepMC::GenEvent* hepmcevt){ 
00088 
00089 //   // print number of events
00090 //   if (hepmcevt->event_number()%1 == 0) {
00091   //    std::cout << "Processing event " << hepmcevt->event_number() << std::endl;
00092   //   }
00093   
00094  
00095 // elastic analysis
00096   int pcount=0;
00097   int pint=0;
00098   double const plhc=7000.;
00099   
00100   HepMC::GenParticle* pi1=0;
00101   HepMC::GenParticle* pi2=0;
00102   HepMC::GenParticle* p1=0;
00103   HepMC::GenParticle* p2=0;
00104 
00105   // loop over the final state particles 
00106   for ( HepMC::GenEvent::particle_const_iterator p =  hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00107     {
00108       int pid = (*p)->pdg_id();
00109       //      (*p)->print();
00110       //if( ((*p)->status() == -12))  {
00111       if( (hepmcevt->beam_particles().first->pdg_id()==2212 && !(*p)->production_vertex())
00112           || (hepmcevt->beam_particles().second->pdg_id()==2212 && !(*p)->production_vertex()) ) {
00113         if( pid == 2212){ 
00114           pint++;
00115           //(*p)->print();
00116           if(pint == 1){pi1=(*p);}
00117           if(pint == 2){pi2=(*p);}
00118         }
00119       }
00120 
00121       if( ((*p)->status()%1000 == 1) && (!(*p)->end_vertex()) )  {
00122         //if( (!(*p)->end_vertex()) )  {
00123         if( pid == 2212){ 
00124           pcount++;
00125           if(pcount == 1){p1=(*p);}
00126           if(pcount == 2){p2=(*p);}
00127           if(pcount > 2){cout << "Strange: More than two protons in this event!" << "\n";}
00128         }
00129       }
00130     }
00131     if(pint == 2){
00132       pi1_en->Fill((pi1->momentum().e())-plhc);
00133       double divx = (pi1->momentum().px())/(pi1->momentum().pz());
00134       double divy = (pi1->momentum().py())/(pi1->momentum().pz());
00135       pi1_divx->Fill(divx*1.E6);
00136       pi1_divy->Fill(divy*1.E6);
00137 
00138       pi2_en->Fill((pi2->momentum().e())-plhc);
00139       divx = (pi2->momentum().px())/(pi2->momentum().pz());
00140       divy = (pi2->momentum().py())/(pi2->momentum().pz());
00141       pi2_divx->Fill(divx*1.E6);
00142       pi2_divy->Fill(divy*1.E6);
00143     }
00144     
00145     if(pcount==2){
00146      n_proton->Fill(pcount);
00147      p1_pt->Fill(p1->momentum().perp());
00148      p1_eta->Fill(p1->momentum().eta());
00149      p1_phi->Fill(p1->momentum().phi());
00150      p1_theta->Fill((p1->momentum().theta())*1.E6); //micro-radians
00151      p2_pt->Fill(p2->momentum().perp());
00152      p2_eta->Fill(p2->momentum().eta());
00153      p2_phi->Fill(p2->momentum().phi());
00154      p2_theta->Fill((p2->momentum().theta()-M_PI)*1.E6); //micro-radians
00155 
00156      // t-calculation 
00157      double st1 =(p1->momentum().theta())*plhc;
00158      double st2 =(p2->momentum().theta()-M_PI)*plhc;
00159      double t1 = st1*st1; 
00160      double t2 = st2*st2; 
00161      double tmean = (t1+t2)/2.; // average t-value reconstructed from the scattering angles
00162 
00163      double tau = log10(tmean);
00164      t_spect->Fill(tmean); 
00165      tau_spect->Fill(tau); 
00166     }
00167   
00168   return true;
00169 
00170 } //end of Process

Generated on Wed Aug 31 09:44:48 2011 for HepMCAnalysis by  doxygen 1.4.7