HepMCAnalysis: src/UEAnalysis.cc Source File

00001 #include <iostream>
00002 #include <stdio.h>
00003 #include "HepMC/GenEvent.h"
00004 #include "HepMC/IO_GenEvent.h"
00005 #include "HepMC/GenParticle.h"
00006 #include "HepMC/GenVertex.h"
00007 #include "HepMC/IO_AsciiParticles.h"
00008 #include "HepMC/SimpleVector.h"
00009 #include "HepPDT/ParticleData.hh"
00010 #include "CLHEP/Vector/LorentzVector.h"
00011 #define PI 3.141593
00012 using namespace std;
00013 
00014 // ROOT headers
00015 #include "TH1.h"
00016 #include "TH2.h"
00017 #include "TFile.h"
00018 #include "TMath.h"
00019 #include "TLorentzVector.h"
00020 
00021 #include "fastjet/PseudoJet.hh"
00022 #include "fastjet/ClusterSequence.hh"
00023 #include "fastjet/JetDefinition.hh"
00024 #include "fastjet/SISConePlugin.hh"
00025 
00026 //UE analysis header
00027 #include "../include/UEAnalysis.h"
00028 
00029 //**********************
00030 
00031 /// empty default constructor
00032 UEAnalysis::UEAnalysis() {
00033 
00034 }
00035 
00036 /// empty default destructor
00037 UEAnalysis::~UEAnalysis() {
00038 
00039 }
00040 
00041 int UEAnalysis::Init(double tr_max_eta, double tr_min_pt) {
00042   // specify default eta cut
00043   m_max_eta=tr_max_eta;
00044   
00045   // specify default pt cut
00046   m_min_pt=tr_min_pt;
00047   
00048   // default Output file name
00049   m_outputFileName="UEAnalysis.root";
00050   m_outputRootDir="UE";
00051   
00052   char name1[20],title1[20];
00053   char name2[20],title2[20];
00054   char name3[20],title3[20];
00055   //declaration of histograms
00056   for (int i=0; i<26; i++){
00057     //just give a different name to the histograms
00058     //for <Ncharged> 
00059     sprintf(name1,"NchargedToward_%2d",i);
00060     sprintf(name2,"NchargedTransverse_%2d",i);
00061     sprintf(name3,"NchargedAway_%2d",i);
00062     
00063     sprintf(title1,"NchargedToward_%2d",i);
00064     sprintf(title2,"NchargedTransverse_%2d",i);
00065     sprintf(title2,"NchargedAway_%2d",i);
00066     
00067     m_NchargedToward[i] = new TH1D(name1,title1, 20, 0, 100);
00068     m_NchargedTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
00069     m_NchargedAway[i] = new TH1D(name3,title3, 20, 0, 100);
00070 
00071    //for <Ptsum>
00072     sprintf(name1,"PtsumToward_%2d",i);
00073     sprintf(name2,"PtsumTransverse_%2d",i);
00074     sprintf(name3,"PtsumAway_%2d",i);
00075     
00076     sprintf(title1,"PtsumToward_%2d",i);
00077     sprintf(title2,"PtsumTransverse_%2d",i);
00078     sprintf(title2,"PtsumAway_%2d",i);
00079     
00080     m_PtsumToward[i] = new TH1D(name1,title1, 20, 0, 100);
00081     m_PtsumTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
00082     m_PtsumAway[i] = new TH1D(name3,title3, 20, 0, 100);
00083 
00084   }
00085 
00086 
00087   //setting up the variable bin size
00088   m_nbin_pT = 25;
00089   
00090   for (int i=0; i<11; i++) {
00091     m_nbinRange_pT[i] = i;
00092   }
00093   for (int i=11; i<21; i++) {
00094     m_nbinRange_pT[i] = 10.0+ 3.0*(i-10);
00095   }
00096   m_nbinRange_pT[21] = 50.; 
00097   m_nbinRange_pT[22] = 60.; 
00098   m_nbinRange_pT[23] = 70.; 
00099   m_nbinRange_pT[24] = 80.; 
00100   m_nbinRange_pT[25] = 120.; 
00101   
00102   //booking histograms
00103   m_NchargedMeanToward = baseAnalysis::initHistVariableBin( string("NchargedMeanToward"), string("NchargedMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00104   m_NchargedMeanTransverse = baseAnalysis::initHistVariableBin( string("NchargedMeanTransverse"), string("NchargedMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00105   m_NchargedMeanAway = baseAnalysis::initHistVariableBin( string("NchargedMeanAway"), string("NchargedMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
00106  
00107   m_PtsumMeanToward = baseAnalysis::initHistVariableBin( string("PtsumMeanToward"), string("PtsumMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00108   m_PtsumMeanTransverse = baseAnalysis::initHistVariableBin( string("PtsumMeanTransverse"), string("PtsumMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00109   m_PtsumMeanAway = baseAnalysis::initHistVariableBin( string("PtsumMeanAway"), string("PtsumMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
00110  
00111   m_Njet =  baseAnalysis::initHist( string("Njet"), string("Njet"), string("jet Number"), 16, -0.5, 15.5);
00112   m_Ptjet =  baseAnalysis::initHist( string("Ptjet"), string("jet p_{T}"), string("jet p_{T}"), 50, 0, 200);
00113   m_Ptjet_log =  baseAnalysis::initHist( string("Ptjet_logscale"), string("jet p_{T}"), string("jet p_{T} - logscale -"), 50, 0, 200);
00114   m_Ptleadingjet =  baseAnalysis::initHist( string("Ptleadingjet"), string("leading jet p_{T}"), string("leading jet p_{T}"), 50, 0, 200);
00115   m_Ptleadingjet_log =  baseAnalysis::initHist( string("Ptleadingjet_logscale"), string("leading jet p_{T}"), 
00116                                                 string("leading jet p_{T} - logscale -"), 50, 0, 200);
00117 
00118   m_charged_particle_pt=baseAnalysis::initHist("charged_particle_pt","p_{T} of charged particles","p_{T}",25,0,25.); 
00119 
00120   return true; 
00121 }
00122 
00123 
00124 int UEAnalysis::Process(HepMC::GenEvent* hepmcevt){ 
00125 
00126   // storage of input particles for jet
00127   double leadingJetPt = 0.0;
00128   CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00129   
00130   // check if jets were found
00131   // if(!m_inclusive_jets.size()) {
00132   //   std::cout<<"UEAnalysis::no jets found"<<std::endl;
00133   //   return 0;
00134   // }
00135   if(m_inclusive_jets.size()) {
00136     // use 4 vector of leading jet to define towards/away/transverse region
00137     lv_leadingJet.setPx(m_inclusive_jets[0].px());
00138     lv_leadingJet.setPy(m_inclusive_jets[0].py());
00139     lv_leadingJet.setPz(m_inclusive_jets[0].pz());
00140     lv_leadingJet.setE(m_inclusive_jets[0].e());
00141     
00142     leadingJetPt = m_inclusive_jets[0].perp(); 
00143     // fill jet histograms
00144     for(unsigned int  i=0; i<m_inclusive_jets.size(); i++) {
00145       m_Ptjet->Fill(m_inclusive_jets[i].perp());
00146       m_Ptjet_log->Fill(m_inclusive_jets[i].perp());
00147     }
00148     m_Njet->Fill(m_inclusive_jets.size()); 
00149     m_Ptleadingjet->Fill(m_inclusive_jets[0].perp());
00150     m_Ptleadingjet_log->Fill(m_inclusive_jets[0].perp());
00151   }
00152   int NchargedToward[25];
00153   int NchargedTransverse[25];
00154   int NchargedAway[25];
00155   
00156   double PtsumToward[25];
00157   double PtsumTransverse[25];
00158   double PtsumAway[25];
00159   
00160   for(int i=0; i<m_nbin_pT; i++){
00161     NchargedToward[i]=0;
00162     NchargedTransverse[i]=0;
00163     NchargedAway[i]=0;
00164     
00165     PtsumToward[i]=0;
00166     PtsumTransverse[i]=0;
00167     PtsumAway[i]=0;
00168   }
00169   
00170   // find the leadjet bin for this event
00171   int ileadjet = -1;
00172   //int ileadjet;
00173   for(int i=0; i<m_nbin_pT; i++) {
00174     if( leadingJetPt > m_nbinRange_pT[i] && leadingJetPt < m_nbinRange_pT[i+1])
00175       ileadjet = i;
00176   }
00177   
00178   // loop over particles in the event
00179   for ( HepMC::GenEvent::particle_const_iterator p =
00180           hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p ){
00181     
00182     // use the pdg id to identify the particles
00183     int pid = (*p)->pdg_id();
00184     
00185     // cut out the neutral particles (charge is not stored in HepMC)
00186     // just remove the neutrinos, gammas, neutral pions and
00187     // neutrons, K0s, ... from the final state
00188     if(!chargedParticle(pid)) continue;
00189     
00190     
00191     // check wether the track directly comes from a Z-Boson
00192     //if( trackfromPID(23,(*p)) ) continue;
00193     
00194     // check wether the track directly comes from a W-Boson
00195     //if( trackfromPID(24,(*p)) ) continue;
00196     //if( trackfromPID(-24,(*p)) ) continue;
00197 
00198     // check for final state particles, be careful, maybe there is another eta cut than 2.5 necessary
00199     if(!IsFinalStateParticle((*p))) continue;
00200 
00201     m_charged_particle_pt->Fill((*p)->momentum().perp());
00202     
00203     // to skip all the events with no jet
00204     if(ileadjet == -1) continue;
00205     //if(!ileadjet) continue;
00206     double pt = (*p)->momentum().perp();
00207     double px = (*p)->momentum().px();
00208     double py = (*p)->momentum().py();
00209     double pz = (*p)->momentum().pz();
00210     double pe = (*p)->momentum().e();
00211     
00212     CLHEP::HepLorentzVector lv_chargedParticle(px,py,pz,pe);
00213     
00214     //calculate for <Ncharged> vs pt of leading jet at toward region
00215     if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < PI/3.0 ) { 
00216       NchargedToward[ileadjet]++;       
00217       PtsumToward[ileadjet]+=pt;
00218     }  
00219     
00220     //calculate for <Ncharged> vs pt of leading jet at Transverse region
00221     if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= PI/3.0 &&
00222         fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < 2.0*PI/3.0 ) { 
00223       NchargedTransverse[ileadjet]++;
00224       PtsumTransverse[ileadjet]+=pt; 
00225     }
00226     //calculate for <Ncharged> vs pt of leading jet at Away region
00227     if( fabs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= 2.0*PI/3.0 ) { 
00228       NchargedAway[ileadjet]++;
00229       PtsumAway[ileadjet]+=pt;
00230     }  
00231       
00232     m_NchargedToward[ileadjet]->Fill(NchargedToward[ileadjet], 1.);
00233     m_NchargedTransverse[ileadjet]->Fill(NchargedTransverse[ileadjet], 1.);
00234     m_NchargedAway[ileadjet]->Fill(NchargedAway[ileadjet], 1.);
00235     
00236     m_PtsumToward[ileadjet]->Fill(PtsumToward[ileadjet], 1.);
00237     m_PtsumTransverse[ileadjet]->Fill(PtsumTransverse[ileadjet], 1.);
00238     m_PtsumAway[ileadjet]->Fill(PtsumAway[ileadjet], 1.);
00239     
00240   }
00241   
00242   return true;
00243 } //end of Process
00244 
00245 std::vector<TH1D*> UEAnalysis::averagedHistograms() {
00246 
00247   //Fill Histograms at Toward, Transverse and Away region
00248   for(int i=0; i<m_nbin_pT; i++){
00249     m_NchargedMeanToward->SetBinContent( i+1, m_NchargedToward[i]->GetMean() );  
00250     m_NchargedMeanToward->SetBinError( i+1, m_NchargedToward[i]->GetRMS() );  
00251     
00252     m_NchargedMeanTransverse->SetBinContent( i+1, m_NchargedTransverse[i]->GetMean() );  
00253     m_NchargedMeanTransverse->SetBinError( i+1, m_NchargedTransverse[i]->GetRMS() );  
00254     
00255     m_NchargedMeanAway->SetBinContent( i+1, m_NchargedAway[i]->GetMean() );  
00256     m_NchargedMeanAway->SetBinError( i+1, m_NchargedAway[i]->GetRMS() );  
00257     
00258     m_PtsumMeanToward->SetBinContent( i+1, m_PtsumToward[i]->GetMean() );  
00259     m_PtsumMeanToward->SetBinError( i+1, m_PtsumToward[i]->GetRMS() );  
00260     
00261     m_PtsumMeanTransverse->SetBinContent( i+1, m_PtsumTransverse[i]->GetMean() );  
00262     m_PtsumMeanTransverse->SetBinError( i+1, m_PtsumTransverse[i]->GetRMS() );  
00263     
00264     m_PtsumMeanAway->SetBinContent( i+1, m_PtsumAway[i]->GetMean() );  
00265     m_PtsumMeanAway->SetBinError( i+1, m_PtsumAway[i]->GetRMS() ); 
00266   }
00267   //m_histVector.insert(m_histVector.end(), m_histVectorVariableBin.begin(), m_histVectorVariableBin.end());
00268   //m_histVector.insert(m_histVector.begin()+2, m_histVectorVariableBin.begin(), m_histVectorVariableBin.end());
00269 
00270   return m_histVector;
00271 } //end of averagedHistograms
00272 
00273 int UEAnalysis::Finalize(TFile* output){
00274   
00275   averagedHistograms();  
00276 
00277   baseAnalysis::Finalize(output);
00278   
00279   return true;
00280 }