diff --git a/ANALYSIS/ANALYSISalice/AliAnalysisTaskMCParticleFilter.cxx b/ANALYSIS/ANALYSISalice/AliAnalysisTaskMCParticleFilter.cxx index 82849a73a96..38141256336 100644 --- a/ANALYSIS/ANALYSISalice/AliAnalysisTaskMCParticleFilter.cxx +++ b/ANALYSIS/ANALYSISalice/AliAnalysisTaskMCParticleFilter.cxx @@ -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" @@ -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(mcEH); - if(!hiH){ - dpmH = dynamic_cast(mcEH); - if(!dpmH){ - tunedH = dynamic_cast(mcEH); - } + hiH = dynamic_cast(mcEH); + if(!hiH){ + dpmH = dynamic_cast(mcEH); + if(!dpmH){ + hmcH = dynamic_cast(mcEH); + if(!hmcH) { + tunedH = dynamic_cast(mcEH); + } } + } } if (hiH || dpmH) colG = dynamic_cast(mcEH); @@ -308,10 +313,11 @@ void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/) if (ccEH) { TList *genHeaders = ccEH->GetHeaders(); for (int imch=0; imchGetEntries(); imch++) { - if(!pyH)pyH = dynamic_cast(genHeaders->At(imch)); - if(!hiH)hiH = dynamic_cast(genHeaders->At(imch)); - if(!colG)colG = dynamic_cast(genHeaders->At(imch)); - if(!dpmH)dpmH = dynamic_cast(genHeaders->At(imch)); + if(!pyH)pyH = dynamic_cast(genHeaders->At(imch)); + if(!hiH)hiH = dynamic_cast(genHeaders->At(imch)); + if(!colG)colG = dynamic_cast(genHeaders->At(imch)); + if(!dpmH)dpmH = dynamic_cast(genHeaders->At(imch)); + if(!hmcH)hmcH = dynamic_cast(genHeaders->At(imch)); } } } @@ -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); @@ -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++){ diff --git a/EVGEN/AliGenReaderHepMC.cxx b/EVGEN/AliGenReaderHepMC.cxx index 97aba4b8265..fcb6b3f9aaf 100644 --- a/EVGEN/AliGenReaderHepMC.cxx +++ b/EVGEN/AliGenReaderHepMC.cxx @@ -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 @@ -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, @@ -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(); diff --git a/STEER/STEERBase/AliGenHepMCEventHeader.cxx b/STEER/STEERBase/AliGenHepMCEventHeader.cxx index 60f417c0fce..fcf467af15f 100644 --- a/STEER/STEERBase/AliGenHepMCEventHeader.cxx +++ b/STEER/STEERBase/AliGenHepMCEventHeader.cxx @@ -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 } @@ -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 } @@ -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), @@ -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 } @@ -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; +} diff --git a/STEER/STEERBase/AliGenHepMCEventHeader.h b/STEER/STEERBase/AliGenHepMCEventHeader.h index 4d5ca661f3d..467f57f2d37 100644 --- a/STEER/STEERBase/AliGenHepMCEventHeader.h +++ b/STEER/STEERBase/AliGenHepMCEventHeader.h @@ -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() {} @@ -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: @@ -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 }; diff --git a/TEvtGen/THepMCParser/THepMCParser.cxx b/TEvtGen/THepMCParser/THepMCParser.cxx index 75fd73f5e14..2e998b4792a 100644 --- a/TEvtGen/THepMCParser/THepMCParser.cxx +++ b/TEvtGen/THepMCParser/THepMCParser.cxx @@ -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 { @@ -632,13 +634,15 @@ void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, listheavy_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(); @@ -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 ""; } diff --git a/TEvtGen/THepMCParser/THepMCParser.h b/TEvtGen/THepMCParser/THepMCParser.h index 2aabfea1916..14074e5287b 100644 --- a/TEvtGen/THepMCParser/THepMCParser.h +++ b/TEvtGen/THepMCParser/THepMCParser.h @@ -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 @@ -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); };