Page MenuHomeHEPForge
Files F19252041

plotBKstarGamma.C
No OneTemporary
Actions

Size
6 KB
Referenced Files
None
Subscribers
None

plotBKstarGamma.C
View Options

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
// Macro to create validation plots for testing FSR for
// B0 -> K*0' gamma, K*0' -> K+ pi- nGamma, where n = 0,1,2,...
// Uses the ROOT files created by genRootDecayChain program
#include <string>
using std::string;
void initHist( TH1D* theHist )
{
theHist->SetDirectory( 0 );
theHist->SetTitleOffset( 1.1, "Y" );
theHist->SetLineWidth( 2 );
TAxis* xAxis = theHist->GetXaxis();
xAxis->SetLabelSize( 0.05 );
xAxis->SetTitleSize( 0.05 );
TAxis* yAxis = theHist->GetYaxis();
yAxis->SetLabelSize( 0.05 );
yAxis->SetTitleSize( 0.05 );
}
void plotBKstarGamma( string inFileName = "BKstarGamma.root" )
{
TFile* theFile = new TFile( inFileName.c_str(), "read" );
TTree* theTree = dynamic_cast<TTree*>( theFile->Get( "Data" ) );
// Create plots of K*0' mass, momentum of 1st and FSR gammas,
// mtm of K+, pi-, and the multiplicity of photons
TH1D* KstMHist = new TH1D( "KstMHist", "", 100, 0.0, 5.0 );
KstMHist->SetXTitle( "K^{*0'} mass (GeV/c^{2})" );
KstMHist->SetYTitle( "Frequency/50 (MeV/c^{2})" );
initHist( KstMHist );
TH1D* gamPHist = new TH1D( "gamPHist", "", 100, 0.0, 3.0 );
gamPHist->SetXTitle( "1st #gamma mtm (GeV/c)" );
gamPHist->SetYTitle( "Frequency/30 (MeV/c^{2})" );
initHist( gamPHist );
TH1D* FSRPHist = new TH1D( "FSRPHist", "", 100, 0.0, 2.0 );
FSRPHist->SetXTitle( "FSR #gamma mtm (GeV/c)" );
FSRPHist->SetYTitle( "Frequency/20 (MeV/c)" );
initHist( FSRPHist );
TH1D* KPHist = new TH1D( "KPHist", "", 100, 0.0, 3.0 );
KPHist->SetXTitle( "K^{+} mtm (GeV/c)" );
KPHist->SetYTitle( "Frequency/30 (MeV/c)" );
initHist( KPHist );
TH1D* piPHist = new TH1D( "piPHist", "", 100, 0.0, 3.0 );
piPHist->SetXTitle( "#pi^{-} mtm (GeV/c)" );
piPHist->SetYTitle( "Frequency/30 (MeV/c)" );
initHist( piPHist );
TH1D* gamNHist = new TH1D( "gamNHist", "", 6, 0, 6 );
gamNHist->SetXTitle( "Number of photons" );
gamNHist->SetYTitle( "Fraction of events" );
initHist( gamNHist );
// Setup the TTree variables
int eventId( 0 ), PDGId( 0 ), pVtx( 0 ), daug( 0 );
double pL( 0.0 ), m( 0.0 );
theTree->SetBranchStatus( "*", 0 );
theTree->SetBranchStatus( "eventId", 1 );
theTree->SetBranchStatus( "PDGId", 1 );
theTree->SetBranchStatus( "daug", 1 );
theTree->SetBranchStatus( "pVtx", 1 );
theTree->SetBranchStatus( "pL", 1 );
theTree->SetBranchStatus( "m", 1 );
theTree->SetBranchAddress( "eventId", &eventId );
theTree->SetBranchAddress( "PDGId", &PDGId );
theTree->SetBranchAddress( "pVtx", &pVtx );
theTree->SetBranchAddress( "daug", &daug );
theTree->SetBranchAddress( "pL", &pL );
theTree->SetBranchAddress( "m", &m );
int nEntries = theTree->GetEntries();
int i( 0 );
int oldEvent( -1 );
int nPhotons( 0 );
int nEvents = theTree->GetMaximum( "eventId" ) + 1;
for ( i = 0; i < nEntries; i++ ) {
theTree->GetEntry( i );
if ( oldEvent != eventId ) {
// We have a new event
// Fill the photon multiplicity histogram
gamNHist->Fill( nPhotons * 1.0 );
// Reset photon counter
nPhotons = 0;
// Set the old event index to the current one
oldEvent = eventId;
}
if ( PDGId == 100313 ) {
// K*0' mass distribution
KstMHist->Fill( m );
}
if ( PDGId == 22 ) {
if ( pVtx == 0 && daug == 2 ) {
// 1st photon momentum
gamPHist->Fill( pL );
} else {
// Other FSR photon momentum
FSRPHist->Fill( pL );
}
// Increment the photon counter
nPhotons++;
}
if ( PDGId == 321 ) {
KPHist->Fill( pL );
}
if ( PDGId == -211 ) {
piPHist->Fill( pL );
}
}
// Create the plots
gROOT->SetStyle( "Plain" );
gStyle->SetOptStat( 0 );
TCanvas* theCanvas = new TCanvas( "theCanvas", "", 1200, 1100 );
theCanvas->UseCurrentStyle();
theCanvas->Divide( 2, 3 );
// Same normalisation number for all plots (except FSR momenta)
double scale = 1.0 / ( nEvents * 1.0 );
// K*0' mass
theCanvas->cd( 1 );
gPad->SetLogy( 0 );
KstMHist->Scale( scale );
KstMHist->SetMaximum( 0.16 );
KstMHist->Draw();
// 1st photon momentum
theCanvas->cd( 2 );
gPad->SetLogy( 1 );
gamPHist->Scale( scale );
gamPHist->SetMaximum( 0.5 );
gamPHist->Draw();
// K+ momentum
theCanvas->cd( 3 );
KPHist->Scale( scale );
KPHist->SetMaximum( 0.025 );
KPHist->Draw();
// FSR photon momenta
theCanvas->cd( 4 );
gPad->SetLogy( 1 );
// Scale is equal to the area of the histogram, since there may
// be events with more than one FSR photon
double scale1( 1.0 );
double FSRIntegral = FSRPHist->Integral();
if ( FSRIntegral > 0.0 ) {
scale1 = 1.0 / FSRIntegral;
}
FSRPHist->Scale( scale1 );
FSRPHist->SetMaximum( 1.0 );
FSRPHist->Draw();
// pi- momentum
theCanvas->cd( 5 );
piPHist->Scale( scale );
piPHist->SetMaximum( 0.025 );
piPHist->Draw();
// Photon multiplicity
theCanvas->cd( 6 );
gamNHist->Scale( scale );
gamNHist->SetMaximum( 1.0 );
gamNHist->Draw();
theCanvas->cd( 1 );
TLatex latex;
latex.SetNDC();
latex.SetTextSize( 0.045 );
latex.DrawLatex(
0.1, 0.95,
"B0 #rightarrow K^{*0'} #gamma, K^{*0'} #rightarrow K^{+} #pi^{-} with FSR" );
theCanvas->Print( "plotBKstarGamma.eps" );
//theCanvas->Print("plotBKstarGamma.png");
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:13 AM (11 h, 55 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6561970
Default Alt Text
plotBKstarGamma.C (6 KB)

Event Timeline

HEPForge.org · .