00001
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
00013 #include "CLHEP/Vector/LorentzVector.h"
00014
00015
00016 #include "TH1.h"
00017 #include "TMath.h"
00018 #include "TLorentzVector.h"
00019 #include "TFile.h"
00020
00021
00022 #include "fastjet/PseudoJet.hh"
00023 #include "fastjet/ClusterSequence.hh"
00024 #include "fastjet/JetDefinition.hh"
00025 #include "fastjet/SISConePlugin.hh"
00026
00027
00028 #include "../include/ElasScatAnalysis.h"
00029
00030 using namespace std;
00031
00032
00033
00034
00035 ElasScatAnalysis::ElasScatAnalysis()
00036 {
00037 }
00038
00039 ElasScatAnalysis::~ElasScatAnalysis()
00040 {
00041 }
00042
00043
00044 int ElasScatAnalysis::Init(double tr_max_eta, double tr_min_pt)
00045 {
00046
00047 m_max_eta=tr_max_eta;
00048
00049
00050 m_min_pt=tr_min_pt;
00051
00052
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
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
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
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
00090
00091
00092
00093
00094
00095
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
00106 for ( HepMC::GenEvent::particle_const_iterator p = hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p )
00107 {
00108 int pid = (*p)->pdg_id();
00109
00110
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
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
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);
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);
00155
00156
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.;
00162
00163 double tau = log10(tmean);
00164 t_spect->Fill(tmean);
00165 tau_spect->Fill(tau);
00166 }
00167
00168 return true;
00169
00170 }