Skip to content

Commit

Permalink
Adding fields for cross section, pthard and ntrials to hepmc event he…
Browse files Browse the repository at this point in the history 
…ader

Fields are mainly of relevance for pp generators (Herwig,
SHERPA, ...)

Propagate cross section, number of trials and pt-hard from
HepMC header to AliAODMCHeader

Including a few indentation fix exporting tab to space
  • Loading branch information
@mfasDa @chiarazampolli
mfasDa authored and chiarazampolli committed Nov 12, 2020
1 parent 86b5e22 commit 46ecad4
Show file tree
Hide file tree
Showing 6 changed files with 106 additions and 23 deletions.
32 changes: 22 additions & 10 deletions ANALYSIS/ANALYSISalice/AliAnalysisTaskMCParticleFilter.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
#include "AliGenPythiaEventHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliGenEventHeaderTunedPbPb.h"
#include "AliGenHepMCEventHeader.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"
#include "AliAODPid.h"
Expand Down Expand Up @@ -283,17 +284,21 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
AliCollisionGeometry *colG = 0;
AliGenDPMjetEventHeader *dpmH = 0;
AliGenEventHeaderTunedPbPb *tunedH = 0;
AliGenHepMCEventHeader *hmcH = 0;

// it can be only one save some casts
// assuming PYTHIA and HIJING are the most likely ones...
if(!pyH){
hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
if(!hiH){
dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
if(!dpmH){
tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
}
hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
if(!hiH){
dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
if(!dpmH){
hmcH = dynamic_cast<AliGenHepMCEventHeader*>(mcEH);
if(!hmcH) {
tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
}
}
}
}

if (hiH || dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);
Expand All @@ -308,10 +313,11 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
if (ccEH) {
TList *genHeaders = ccEH->GetHeaders();
for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
if(!hmcH)hmcH = dynamic_cast<AliGenHepMCEventHeader*>(genHeaders->At(imch));
}
}
}
Expand All @@ -320,6 +326,7 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
if(hiH)ntrials = hiH->Trials();
if(dpmH)ntrials = dpmH->Trials();
if(pyH)ntrials = pyH->Trials();
if(hmcH)ntrials = hmcH->ntrials();
if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);


Expand All @@ -334,6 +341,11 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
fAODMcHeader->SetCrossSection(tunedH->GetCentrality());
}

if(hmcH) {
fAODMcHeader->SetCrossSection(hmcH->sigma_gen());
fAODMcHeader->SetPtHard(hmcH->pthard());
}


Int_t j=0;
for (Int_t ip = 0; ip < np; ip++){
Expand Down
10 changes: 8 additions & 2 deletions EVGEN/AliGenReaderHepMC.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ void AliGenReaderHepMC::Init()

Int_t AliGenReaderHepMC::NextEvent()
{
const Double_t kPbToMb = 1e-9; // Cross section units: pb to mb;
// Clean memory
if (fGenEvent) delete fGenEvent;
// Read the next event
Expand All @@ -59,7 +60,8 @@ Int_t AliGenReaderHepMC::NextEvent()
fParticleIterator->Reset();
THepMCParser::HeavyIonHeader_t heavyIonHeader;
THepMCParser::PdfHeader_t pdfHeader;
THepMCParser::ParseGenEvent2HeaderStructs(fGenEvent,heavyIonHeader,pdfHeader,true,true);
THepMCParser::CrossSectionHeader_t csHeader;
THepMCParser::ParseGenEvent2HeaderStructs(fGenEvent,heavyIonHeader,pdfHeader,csHeader,true,true,true);
fGenEventHeader = new AliGenHepMCEventHeader(
heavyIonHeader.Ncoll_hard,
heavyIonHeader.Npart_proj,
Expand All @@ -82,7 +84,11 @@ Int_t AliGenReaderHepMC::NextEvent()
pdfHeader.x2,
pdfHeader.scalePDF,
pdfHeader.pdf1,
pdfHeader.pdf2
pdfHeader.pdf2,
csHeader.sigma_gen * kPbToMb,
csHeader.sigma_err * kPbToMb,
csHeader.pt_hard,
csHeader.ntrials
);
// propagate the event weight from HepMC to the event header
HepMC::WeightContainer weights = fGenEvent->weights();
Expand Down
31 changes: 27 additions & 4 deletions STEER/STEERBase/AliGenHepMCEventHeader.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,11 @@ AliGenHepMCEventHeader::AliGenHepMCEventHeader():
fx2(0.0),
fscalePDF(0.0),
fpdf1(0.0),
fpdf2(0.0)
fpdf2(0.0),
fSigmaGen(0.0),
fSigmaErr(0.0),
fPtHard(-1.0),
fNtrials(0.0)
{
// Default Constructor
}
Expand Down Expand Up @@ -69,7 +73,11 @@ AliGenHepMCEventHeader::AliGenHepMCEventHeader(const char* name):
fx2(0.0),
fscalePDF(0.0),
fpdf1(0.0),
fpdf2(0.0)
fpdf2(0.0),
fSigmaGen(0.0),
fSigmaErr(0.0),
fPtHard(-1.0),
fNtrials(0.0)
{
// Constructor
}
Expand All @@ -96,7 +104,11 @@ AliGenHepMCEventHeader::AliGenHepMCEventHeader(
Double_t var_x2, // fraction of beam momentum carried by second parton ("target side")
Double_t var_scalePDF, // Q-scale used in evaluation of PDF's (in GeV)
Double_t var_pdf1, // PDF (id1, x1, Q) - x*f(x)
Double_t var_pdf2 // PDF (id2, x2, Q) - x*f(x)
Double_t var_pdf2, // PDF (id2, x2, Q) - x*f(x)
Float_t sigma_gen, // Cross section (for pp collisions)
Float_t sigma_err, // Error of the cross section (for pp collisions, if supported by generator)
Float_t pthard, // pt-hard (event scale, for pp collisions)
Int_t ntrials // number of trials (in case of jet-jet productions, for pp collisions)
):
fNcoll_hard(var_Ncoll_hard),
fNpart_proj(var_Npart_proj),
Expand All @@ -119,7 +131,11 @@ AliGenHepMCEventHeader::AliGenHepMCEventHeader(
fx2(var_x2),
fscalePDF(var_scalePDF),
fpdf1(var_pdf1),
fpdf2(var_pdf2)
fpdf2(var_pdf2),
fSigmaGen(sigma_gen),
fSigmaErr(sigma_err),
fPtHard(pthard),
fNtrials(ntrials)
{
// The Constructor
}
Expand Down Expand Up @@ -151,3 +167,10 @@ Bool_t AliGenHepMCEventHeader::PDFValid() {
fpdf1 != 0.0 ||
fpdf2 != 0.0;
}

Bool_t AliGenHepMCEventHeader::CrossSectionValid() {
return fSigmaGen != 0. ||
fSigmaErr != 0. ||
fPtHard != 0. ||
fNtrials != 0;
}
20 changes: 18 additions & 2 deletions STEER/STEERBase/AliGenHepMCEventHeader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,11 @@ class AliGenHepMCEventHeader : public AliGenEventHeader
Double_t x2, // fraction of beam momentum carried by second parton ("target side")
Double_t scalePDF, // Q-scale used in evaluation of PDF's (in GeV)
Double_t pdf1, // PDF (id1, x1, Q) - x*f(x)
Double_t pdf2 // PDF (id2, x2, Q) - x*f(x)
Double_t pdf2, // PDF (id2, x2, Q) - x*f(x)
Float_t sigma_gen, // Cross section (for pp collisions)
Float_t sigma_err, // Error of the cross section (for pp collisions)
Float_t pthard, // pt-hard (event scale, for pp collisions)
Int_t ntrials // number of trials (in case of jet-jet productions, for pp collisions)
);
virtual ~AliGenHepMCEventHeader() {}

Expand Down Expand Up @@ -65,9 +69,15 @@ class AliGenHepMCEventHeader : public AliGenEventHeader
Double_t pdf1() const {return fpdf1;} // PDF (id1, x1, Q) - x*f(x)
Double_t pdf2() const {return fpdf2;} // PDF (id2, x2, Q) - x*f(x)

Float_t sigma_gen() const { return fSigmaGen; }
Float_t sigma_err() const { return fSigmaErr; }
Float_t pthard() const { return fPtHard; }
Int_t ntrials() const { return fNtrials; }

// convenience functions to check if the headers are containing information
Bool_t HeavyIonInfoValid();
Bool_t PDFValid();
Bool_t CrossSectionValid();

protected:

Expand Down Expand Up @@ -95,7 +105,13 @@ class AliGenHepMCEventHeader : public AliGenEventHeader
Double_t fpdf1; // PDF (id1, x1, Q) - x*f(x)
Double_t fpdf2; // PDF (id2, x2, Q) - x*f(x)

ClassDef(AliGenHepMCEventHeader, 2) // Event header for HepMC event
// For pp generators
Double_t fSigmaGen; ///< Cross section
Double_t fSigmaErr; ///< Error on the cross section
Double_t fPtHard; ///< Pt of the hard process
Int_t fNtrials; ///< Number of trials

ClassDef(AliGenHepMCEventHeader, 3) // Event header for HepMC event
};


Expand Down
27 changes: 23 additions & 4 deletions TEvtGen/THepMCParser/THepMCParser.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,12 @@ void THepMCParser::init(HepMC::IO_BaseClass * events)
fTree->Branch("HeavyIonInfo", &heavyIonHeader, THepMCParser::fgHeavyIonHeaderBranchString);
THepMCParser::PdfHeader_t pdfHeader;
fTree->Branch("PdfInfo", &pdfHeader, THepMCParser::fgPdfHeaderBranchString);
THepMCParser::CrossSectionHeader_t csHeader;
fTree->Branch("CrossSectionInfo", &csHeader, THepMCParser::fgCsHeaderBranchString);
HepMC::GenEvent* evt = 0; // nullptr
while ((evt = events->read_next_event())) {
string errMsg1 = ParseGenEvent2TCloneArray(evt,array);
string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader);
string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader,csHeader);
if (errMsg1.length() == 0 && errMsg2.length() == 0) {
fTree->Fill();
} else {
Expand Down Expand Up @@ -632,13 +634,15 @@ void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, list<HepMC::GenVerte

const char * THepMCParser::fgHeavyIonHeaderBranchString = "Ncoll_hard/I:Npart_proj:Npart_targ:Ncoll:spectator_neutrons:spectator_protons:N_Nwounded_collisions:Nwounded_N_collisions:Nwounded_Nwounded_collisions:impact_parameter/F:event_plane_angle:eccentricity:sigma_inel_NN";
const char * THepMCParser::fgPdfHeaderBranchString = "id1/I:id2:pdf_id1:pdf_id2:x1/D:x2:scalePDF:pdf1:pdf2";
const char * THepMCParser::fgCsHeaderBranchString = "sigma_gen/D:sigma_err,pt_hard,ntrials";

string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf)
string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, CrossSectionHeader_t &csheader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf, bool fillZeroOnMissingCs)
{
HepMC::HeavyIon * heavyIon = genEvent->heavy_ion();
HepMC::PdfInfo * pdfInfo = genEvent->pdf_info();
if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf)) {
return "HeavyIonInfo and/or PdfInfo not defined for this event, did you read it with IO_GenEvent?";
HepMC::GenCrossSection *csInfo = genEvent->cross_section();
if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf) || (!csInfo && !fillZeroOnMissingCs)) {
return "HeavyIonInfo and/or PdfInfo and/or CrossSectionInfo not defined for this event, did you read it with IO_GenEvent?";
}
if (heavyIon) {
heavyIonHeader.Ncoll_hard = heavyIon->Ncoll_hard();
Expand Down Expand Up @@ -690,6 +694,21 @@ string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, Hea
pdfHeader.pdf1 = 0.0;
pdfHeader.pdf2 = 0.0;
}
if(csInfo) {
csheader.sigma_gen = csInfo->cross_section();
csheader.sigma_err = csInfo->cross_section_error();
csheader.pt_hard = genEvent->event_scale();
if(genEvent->weights().has_key("NTrial")) {
csheader.ntrials = genEvent->weights()["NTrial"];
} else {
csheader.ntrials = 1;
}
} else {
csheader.sigma_gen = 0;
csheader.sigma_err = 0;
csheader.pt_hard = 0;
csheader.ntrials = 0;
}
return "";
}

Expand Down
9 changes: 8 additions & 1 deletion TEvtGen/THepMCParser/THepMCParser.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,13 @@ class THepMCParser {
Double_t pdf1; // PDF (id1, x1, Q) - x*f(x)
Double_t pdf2; // PDF (id2, x2, Q) - x*f(x)
};
static const char * fgCsHeaderBranchString; // "sigma_gen/D:sigma_err,pt_hard,ntrials";
struct CrossSectionHeader_t {
Double_t sigma_gen; ///< gen. cross section
Double_t sigma_err; ///< error of the gen. cross section
Double_t pt_hard; ///< pt of the hard process
Int_t ntrials; ///< number of trials
};

// Default constructor/destructor stuff, don't inherit from this class unless you handle the tree pointer
inline THepMCParser() : fTree(0) {;} // nullptr in c++11
Expand Down Expand Up @@ -114,7 +121,7 @@ class THepMCParser {
// heavy ions or parton distribution functions available. This function will pull them
// out if available.
// The caller must supply allocated structures which will then be filled.
static std::string ParseGenEvent2HeaderStructs(HepMC::GenEvent *, HeavyIonHeader_t &, PdfHeader_t &, bool fillZeroOnMissingHeavyIon = true, bool fillZeroOnMissingPdf = true);
static std::string ParseGenEvent2HeaderStructs(HepMC::GenEvent *, HeavyIonHeader_t &, PdfHeader_t &, CrossSectionHeader_t &, bool fillZeroOnMissingHeavyIon = true, bool fillZeroOnMissingPdf = true, bool fillZeroOnMissingCs = true);


};
Expand Down

0 comments on commit 46ecad4

Please sign in to comment.