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 //starting functions
00032 UEAnalysis::UEAnalysis()
00033 {
00034 
00035 }
00036 
00037 UEAnalysis::~UEAnalysis()
00038 {
00039   while (!m_histVector.empty())
00040     {
00041       delete m_histVector.back();
00042       m_histVector.pop_back();
00043     }
00044   
00045   while (!m_histVectorVariableBin.empty())
00046     {
00047       delete m_histVectorVariableBin.back();
00048       m_histVectorVariableBin.pop_back();
00049     }
00050 }
00051 
00052 int UEAnalysis::Init(double tr_max_eta, double tr_min_pt)
00053 {
00054   // specify default eta cut
00055   m_max_eta=tr_max_eta;
00056   
00057   // specify default pt cut
00058   m_min_pt=tr_min_pt;
00059   
00060   // default Output file name
00061   m_outputFileName="UEAnalysis.root";
00062   m_outputRootDir="UE";
00063 
00064   char name1[20],title1[20];
00065   char name2[20],title2[20];
00066   char name3[20],title3[20];
00067   //declaration of histograms
00068   for (int i=0; i<26; i++){
00069     //just give a different name to the histograms
00070     //for <Ncharged> 
00071     sprintf(name1,"NchargedToward_%2d",i);
00072     sprintf(name2,"NchargedTransverse_%2d",i);
00073     sprintf(name3,"NchargedAway_%2d",i);
00074     
00075     sprintf(title1,"NchargedToward_%2d",i);
00076     sprintf(title2,"NchargedTransverse_%2d",i);
00077     sprintf(title2,"NchargedAway_%2d",i);
00078     
00079     m_NchargedToward[i] = new TH1D(name1,title1, 20, 0, 100);
00080     m_NchargedTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
00081     m_NchargedAway[i] = new TH1D(name3,title3, 20, 0, 100);
00082 
00083    //for <Ptsum>
00084     sprintf(name1,"PtsumToward_%2d",i);
00085     sprintf(name2,"PtsumTransverse_%2d",i);
00086     sprintf(name3,"PtsumAway_%2d",i);
00087     
00088     sprintf(title1,"PtsumToward_%2d",i);
00089     sprintf(title2,"PtsumTransverse_%2d",i);
00090     sprintf(title2,"PtsumAway_%2d",i);
00091     
00092     m_PtsumToward[i] = new TH1D(name1,title1, 20, 0, 100);
00093     m_PtsumTransverse[i] = new TH1D(name2,title2, 20, 0, 100);
00094     m_PtsumAway[i] = new TH1D(name3,title3, 20, 0, 100);
00095 
00096   }
00097 
00098 
00099   //setting up the variable bin size
00100   m_nbin_pT = 25;
00101   
00102   for (int i=0; i<11; i++)
00103     {
00104       m_nbinRange_pT[i] = i;
00105     }
00106   for (int i=11; i<21; i++)
00107     {
00108       m_nbinRange_pT[i] = 10.0+ 3.0*(i-10);
00109     }
00110   m_nbinRange_pT[21] = 50.; 
00111   m_nbinRange_pT[22] = 60.; 
00112   m_nbinRange_pT[23] = 70.; 
00113   m_nbinRange_pT[24] = 80.; 
00114   m_nbinRange_pT[25] = 120.; 
00115   
00116   //booking histograms
00117   m_NchargedMeanToward = baseAnalysis::initHistVariableBin( string("NchargedMeanToward"), string("NchargedMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00118   m_NchargedMeanTransverse = baseAnalysis::initHistVariableBin( string("NchargedMeanTransverse"), string("NchargedMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00119   m_NchargedMeanAway = baseAnalysis::initHistVariableBin( string("NchargedMeanAway"), string("NchargedMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
00120  
00121   m_PtsumMeanToward = baseAnalysis::initHistVariableBin( string("PtsumMeanToward"), string("PtsumMeanToward"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00122   m_PtsumMeanTransverse = baseAnalysis::initHistVariableBin( string("PtsumMeanTransverse"), string("PtsumMeanTransverse"), string("Pt_leadingJet"), m_nbin_pT, m_nbinRange_pT);
00123   m_PtsumMeanAway = baseAnalysis::initHistVariableBin( string("PtsumMeanAway"), string("PtsumMeanAway"), string("Pt_leadingJet"),m_nbin_pT, m_nbinRange_pT);
00124  
00125   m_Njet =  baseAnalysis::initHist( string("Njet"), string("Njet"), string("jet Number"), 16, -0.5, 15.5);
00126   m_Ptjet =  baseAnalysis::initHist( string("Ptjet"), string("Ptjet"), string("jet Pt"), 100, 0, 300);
00127   m_Ptjet_log =  baseAnalysis::initHist( string("Ptjet_logscale"), string("Ptjet"), string("jet Pt - logscale -"), 100, 0, 300);
00128   m_Ptleadingjet =  baseAnalysis::initHist( string("Ptleadingjet"), string("Ptleadingjet"), string("leadingjet Pt"), 100, 0, 300);
00129   m_Ptleadingjet_log =  baseAnalysis::initHist( string("Ptleadingjet_logscale"), string("Ptleadingjet"), 
00130                                                 string("leadingjet Pt - logscale -"), 100, 0, 300);
00131 
00132   m_charged_particle_pt=baseAnalysis::initHist("charged_particle_pt","p_{T} of charged particles","p_{T}",25,0,25.); 
00133   
00134   for(unsigned int hist=0; hist< m_histVector.size(); hist++) {
00135     std::ostringstream counter_s;
00136     counter_s << hist;
00137     std::string histogram_name=m_outputRootDir;
00138     histogram_name+="_";
00139     histogram_name+=counter_s.str();
00140     histogram_name+="_";
00141     histogram_name+=m_histVector[hist]->GetName();
00142     m_histVector[hist]->SetName(histogram_name.c_str());
00143   }
00144   return SUCCESS; 
00145 }
00146 
00147 
00148 int UEAnalysis::Process(HepMC::GenEvent* hepmcevt){ 
00149   
00150   // storage of input particles for jet
00151   double leadingJetPt = 0.0;
00152   CLHEP::HepLorentzVector lv_leadingJet(0,0,0,0);
00153   
00154   // if jet finder wasn't run before do it now:
00155   if(!m_inclusive_jets.size()) {
00156     int test = FindJet(hepmcevt);
00157     if(test<=0){
00158       cout <<"no jets found - return failure" << endl;
00159       //return FAILURE;
00160     }
00161   }
00162 
00163   // return if no jet was found
00164   if( ! (m_inclusive_jets.size() && m_clust_seq)) {
00165     cout << "UEAnalysis::ProcessEvent() no jet found - return " << endl;
00166     //return FAILURE;
00167     return 0;
00168   }
00169   
00170   // use 4 vector of leading jet to define towards/away/transverse region
00171   lv_leadingJet.setPx(m_inclusive_jets[0].px());
00172   lv_leadingJet.setPy(m_inclusive_jets[0].py());
00173   lv_leadingJet.setPz(m_inclusive_jets[0].pz());
00174   lv_leadingJet.setE(m_inclusive_jets[0].e());
00175   
00176   leadingJetPt = m_inclusive_jets[0].perp(); 
00177   
00178   // fill jet histograms
00179   for(unsigned int  i=0; i<m_inclusive_jets.size(); i++) {
00180     m_Ptjet->Fill(m_inclusive_jets[i].perp());
00181     m_Ptjet_log->Fill(m_inclusive_jets[i].perp());
00182   }
00183   m_Njet->Fill(m_inclusive_jets.size()); 
00184   m_Ptleadingjet->Fill(m_inclusive_jets[0].perp());
00185   m_Ptleadingjet_log->Fill(m_inclusive_jets[0].perp());
00186   
00187   int NchargedToward[25];
00188   int NchargedTransverse[25];
00189   int NchargedAway[25];
00190 
00191   double PtsumToward[25];
00192   double PtsumTransverse[25];
00193   double PtsumAway[25];
00194   
00195   for(int i=0; i<m_nbin_pT; i++){
00196     NchargedToward[i]=0;
00197     NchargedTransverse[i]=0;
00198     NchargedAway[i]=0;
00199 
00200     PtsumToward[i]=0;
00201     PtsumTransverse[i]=0;
00202     PtsumAway[i]=0;
00203   }
00204   
00205   // find the leadjet bin for this event
00206   int ileadjet = -1;
00207   //int ileadjet;
00208   for(int i=0; i<m_nbin_pT; i++) {
00209     if( leadingJetPt > m_nbinRange_pT[i] && leadingJetPt < m_nbinRange_pT[i+1])
00210       ileadjet = i;
00211   }
00212   
00213   // loop over particles in the event
00214   for ( HepMC::GenEvent::particle_const_iterator p =
00215           hepmcevt->particles_begin(); p != hepmcevt->particles_end(); ++p ){
00216     
00217     // use the pdg id to identify the particles
00218     int pid = (*p)->pdg_id();
00219     
00220     // now collect the tracks for the Jet algorithm
00221     // first of all we need stable particle
00222     if (!((*p)->status() == 1)) continue;
00223     
00224     // cut out the neutral particles (charge is not stored in HepMC)
00225     // just remove the neutrinos, gammas, neutral pions and
00226     // neutrons, K0s, ... from the final state
00227     if(chargedParticle(pid)==NOTCHARGED) continue;
00228     
00229     // They need to have a production Vertex --> we will not select incoming Protons from the Beam
00230     if (!(*p)->production_vertex()) continue;
00231     // And since they are stable they should not have a decay vertex.
00232     if ((*p)->end_vertex()) continue;
00233     
00234     // check wether the track directly comes from a Z-Boson
00235     if( trackfromPID(23,(*p)) ) continue;
00236     
00237     // check wether the track directly comes from a W-Boson
00238     if( trackfromPID(24,(*p)) ) continue;
00239     if( trackfromPID(-24,(*p)) ) continue;
00240     
00241     // set the pt cut for charged particles 
00242     double pt = (*p)->momentum().perp(); 
00243     if(pt < m_min_pt) continue;
00244     
00245     // set the eta cut for charged particles 
00246     double eta = (*p)->momentum().eta(); 
00247     if(abs(eta) > m_max_eta) continue;
00248     
00249     m_charged_particle_pt->Fill((*p)->momentum().perp());
00250     
00251     // to skip all the events with no jet
00252     if(ileadjet == -1) continue;
00253     //if(!ileadjet) continue;
00254     
00255     double px = (*p)->momentum().px();
00256     double py = (*p)->momentum().py();
00257     double pz = (*p)->momentum().pz();
00258     double pe = (*p)->momentum().e();
00259     
00260     CLHEP::HepLorentzVector lv_chargedParticle(px,py,pz,pe);
00261     
00262     //calculate for <Ncharged> vs pt of leading jet at toward region
00263     if( abs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < PI/3.0 )
00264      { NchargedToward[ileadjet]++;      
00265        PtsumToward[ileadjet]+=pt;
00266      }  
00267     
00268     //calculate for <Ncharged> vs pt of leading jet at Transverse region
00269     if( abs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= PI/3.0 &&
00270         abs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) < 2.0*PI/3.0 )
00271      { NchargedTransverse[ileadjet]++;
00272        PtsumTransverse[ileadjet]+=pt; 
00273      }
00274     //calculate for <Ncharged> vs pt of leading jet at Away region
00275     if( abs(lv_chargedParticle.vect().deltaPhi(lv_leadingJet.vect())) >= 2.0*PI/3.0 ) 
00276      { NchargedAway[ileadjet]++;
00277        PtsumAway[ileadjet]+=pt;
00278      }  
00279     
00280     m_NchargedToward[ileadjet]->Fill(NchargedToward[ileadjet], 1.);
00281     m_NchargedTransverse[ileadjet]->Fill(NchargedTransverse[ileadjet], 1.);
00282     m_NchargedAway[ileadjet]->Fill(NchargedAway[ileadjet], 1.);
00283 
00284     m_PtsumToward[ileadjet]->Fill(PtsumToward[ileadjet], 1.);
00285     m_PtsumTransverse[ileadjet]->Fill(PtsumTransverse[ileadjet], 1.);
00286     m_PtsumAway[ileadjet]->Fill(PtsumAway[ileadjet], 1.);
00287     
00288   }
00289   return SUCCESS;
00290 } //end of Process
00291 
00292 
00293 int UEAnalysis::finalize(TFile* output){
00294   
00295   //Fill Histograms at Toward, Transverse and Away region
00296   for(int i=0; i<m_nbin_pT; i++){
00297     m_NchargedMeanToward->SetBinContent( i+1, m_NchargedToward[i]->GetMean() );  
00298     //  m_NchargedMeanToward->SetBinError( i+1, m_NchargedToward[i]->GetRMS() );  
00299     
00300     m_NchargedMeanTransverse->SetBinContent( i+1, m_NchargedTransverse[i]->GetMean() );  
00301     //   m_NchargedMeanTransverse->SetBinError( i+1, m_NchargedTransverse[i]->GetRMS() );  
00302     
00303     m_NchargedMeanAway->SetBinContent( i+1, m_NchargedAway[i]->GetMean() );  
00304     //   m_NchargedMeanAway->SetBinError( i+1, m_NchargedAway[i]->GetRMS() );  
00305 
00306 
00307     m_PtsumMeanToward->SetBinContent( i+1, m_PtsumToward[i]->GetMean() );  
00308     //  m_PtsumMeanToward->SetBinError( i+1, m_PtsumToward[i]->GetRMS() );  
00309     
00310     m_PtsumMeanTransverse->SetBinContent( i+1, m_PtsumTransverse[i]->GetMean() );  
00311     //   m_PtsumMeanTransverse->SetBinError( i+1, m_PtsumTransverse[i]->GetRMS() );  
00312     
00313     m_PtsumMeanAway->SetBinContent( i+1, m_PtsumAway[i]->GetMean() );  
00314     //   m_PtsumMeanAway->SetBinError( i+1, m_PtsumAway[i]->GetRMS() ); 
00315 
00316   }
00317   
00318   baseAnalysis::Finalize(output);
00319   
00320   return SUCCESS;
00321 }