search

cms-sw / cmssw

CMS Offline Software
http://cms-sw.github.io/
Apache License 2.0
1.08k stars 4.32k forks source link

134 PS weights instead of 46 in Pythia 8.3 #36705

Closed agrohsje closed 2 years ago

agrohsje commented 2 years ago

Dear all, producing GEN events in e.g. CMSSW_12_1_0 or CMSSW_11_3_4 (P8.3) I see that the weight vector per event, i.e. Events->Scan("GenEventInfoProduct_generatorGEN.obj.weights_") has 134 entries, while the header contains the correct info. I.e. LuminosityBlocks->Scan("GenLumiInfoHeader_generatorGEN.obj.weightNames_","","colsize=60") gives me 46 PS weights as expected. @SiewYan @mkirsano is one of you available to check what goes wrong? I didn't run P8 standalone yet.
I used this fragment for testing: https://www.desy.de/~agrohsje/hidden/PPD-Run3Autumn21wmLHEGS-00001-fragment.py Cheers, Alexander.

agrohsje commented 2 years ago

assign generators

cmsbuild commented 2 years ago

New categories assigned: generators

@mkirsano,@alberto-sanchez,@SiewYan,@GurpreetSinghChahal,@Saptaparna,@agrohsje you have been requested to review this Pull request/Issue and eventually sign? Thanks

cmsbuild commented 2 years ago

A new Issue was created by @agrohsje .

@Dr15Jones, @perrotta, @dpiparo, @makortel, @smuzaffar, @qliphy can you please review it and eventually sign/assign? Thanks.

cms-bot commands are listed here

agrohsje commented 2 years ago

FYI: @kdlong @smrenna

agrohsje commented 2 years ago

Discussed with @Dongwoon77 and @SanghyunKo as possible extension of relval to catch these changes in relval.

agrohsje commented 2 years ago

@smrenna @SiewYan @mkirsano problem seems to be on the cms side. i tested p8.303 standalone and ended up with the expected number of ps weights in hepmc file.

kdlong commented 2 years ago

@agrohsje, as for catching this in the future, there is a flag in the implementation of PR #32167 that would fail here because the number of weights don't match the weightNames meta info. I've flipped it to false now so it can pass the tests, but ideally it would be true (do fail) and would help us catch such things

SanghyunKo commented 2 years ago

I tested Pythia8Hadronizer.cc with 12_2_0 and had some observations:

  1. Indeed Pythia8 knows that it produces 45 PS weights & adds 45 weights here - fMasterGen->info.nWeights().size() = 45 (fMasterGen is std::unique_ptr<Pythia8::Pythia>)
  2. event()->weights().size() is 89 already at this point: 
    event() = std::make_unique<HepMC::GenEvent>();
bool py8hepmc = toHepMC.fill_next_event(*(fMasterGen.get()), event().get());

    so it seems that something is happening inside the HepMC::Pythia8ToHepMC

Saptaparna commented 2 years ago

Hi! Is anyone looking into HepMC::Pythia8ToHepMC for this PR?

smrenna commented 2 years ago

Okay, I confirm that something strange is happening in Pythia83. I should have a solution soon.

smrenna commented 2 years ago

Okay, I am not convinced that anything is incorrect in Pythia. It is just a matter of interpretation. The config files specifies 44 different variations. Each variation has its own group name. Each variation also has a key/value pair. In principle, one "variation" could be a combination of individual ones. It just so happens that CMS has groups of one key/value pair.

The weights dictionary has an entry for each key/value pair and for each group/value pair. The group/value pair could include the multiplication of several weights, but, again, this is not what CMS does. Therefore, at the end, you get a duplication of weights because each weight is also an individual group.

So, if this is not the information you want, then what do you want, and is it generalizable?

smrenna commented 2 years ago

Note, if you want only 46 weights, you could use this scheme: UncertaintyBands:List = {cmsVar isr:muRfac=0.707 fsr:muRfac=0.707 isr:muRfac=1.414 fsr:muRfac=1.414 isr:muRfac=0.5 fsr:muRfac=0.5 isr:muRfac=2.0 fsr:muRfac=2.0 isr:muRfac=0.25 fsr:muRfac=0.25 isr:muRfac=4.0 fsr:muRfac=4.0 fsr:G2GG:muRfac=0.5 fsr:G2GG:muRfac=2.0 fsr:G2QQ:muRfac=0.5 fsr:G2QQ:muRfac=2.0 fsr:Q2QG:muRfac=0.5 fsr:Q2QG:muRfac=2.0 fsr:X2XG:muRfac=0.5 fsr:X2XG:muRfac=2.0 fsr:G2GG:cNS=-2.0 fsr:G2GG:cNS=2.0 fsr:G2QQ:cNS=-2.0 fsr:G2QQ:cNS=2.0 fsr:Q2QG:cNS=-2.0 fsr:Q2QG:cNS=2.0 fsr:X2XG:cNS=-2.0 fsr:X2XG:cNS=2.0 isr:G2GG:muRfac=0.5 isr:G2GG:muRfac=2.0 isr:G2QQ:muRfac=0.5 isr:G2QQ:muRfac=2.0 isr:Q2QG:muRfac=0.5 isr:Q2QG:muRfac=2.0 isr:X2XG:muRfac=0.5 isr:X2XG:muRfac=2.0 isr:G2GG:cNS=-2.0 isr:G2GG:cNS=2.0 isr:G2QQ:cNS=-2.0 isr:G2QQ:cNS=2.0 isr:Q2QG:cNS=-2.0 isr:Q2QG:cNS=2.0 isr:X2XG:cNS=-2.0 isr:X2XG:cNS=2.0}

This will give you 46 weights -- the baseline, the group called cmsVar, and the 44 individual variations for you to use as you want.

smrenna commented 2 years ago

Also, note, this does not explain:

Events->Scan("GenEventInfoProduct_generator__GEN.obj.weights_")
has 134 entries

Except that 134 = 44 x 3 + 2.

smrenna commented 2 years ago

@agrohsje to help debug this, could you run using the variations defined as above (one big string) and see what are the number of weights in the GenEventInfoProduct? Also, can you list their names/values?

agrohsje commented 2 years ago

Will test.

agrohsje commented 2 years ago

@smrenna, that helped us going down to 90 (2*44+2), ideally we reduce the 2->1 and then we are good to go for run 3. :-) The weight names look good, i.e. I see LuminosityBlocks->Scan("GenLumiInfoHeader_generator_GEN.obj.weightNames","","colsize=30")

smrenna commented 2 years ago

So why are the weights getting doubled by the GenInfoProduct?

Saptaparna commented 2 years ago

@agrohsje Can you please point to a set up that allows me to reproduce the issue? I tried this: https://cms-pdmv.cern.ch/mcm/public/restapi/requests/get_test/EGM-Run3Winter22wmLHEGS-00003 but the weights don't seem to be filled at all.

agrohsje commented 2 years ago

Hi @Saptaparna . I can see 134 weights in your example. Maybe you had a typo in the root command?

agrohsje commented 2 years ago

@smrenna : that is the question. I used the cms list and specification in P8.303: UncertaintyBands:List = {isrRedHi isr:muRfac=0.707,fsrRedHi fsr:muRfac=0.707,isrRedLo isr:muRfac=1.414,fsrRedLo fsr:muRfac=1.414,isrDefHi isr:muRfac=0.5,fsrDefHi fsr:muRfac=0.5,isrDefLo isr:muRfac=2.0,fsrDefLo fsr:muRfac=2.0,isrConHi isr:muRfac=0.25,fsrConHi fsr:muRfac=0.25,isrConLo isr:muRfac=4.0,fsrConLo fsr:muRfac=4.0,fsr_G2GG_muR_dn fsr:G2GG:muRfac=0.5,fsr_G2GG_muR_up fsr:G2GG:muRfac=2.0,fsr_G2QQ_muR_dn fsr:G2QQ:muRfac=0.5,fsr_G2QQ_muR_up fsr:G2QQ:muRfac=2.0,fsr_Q2QG_muR_dn fsr:Q2QG:muRfac=0.5,fsr_Q2QG_muR_up fsr:Q2QG:muRfac=2.0,fsr_X2XG_muR_dn fsr:X2XG:muRfac=0.5,fsr_X2XG_muR_up fsr:X2XG:muRfac=2.0,fsr_G2GG_cNS_dn fsr:G2GG:cNS=-2.0,fsr_G2GG_cNS_up fsr:G2GG:cNS=2.0,fsr_G2QQ_cNS_dn fsr:G2QQ:cNS=-2.0,fsr_G2QQ_cNS_up fsr:G2QQ:cNS=2.0,fsr_Q2QG_cNS_dn fsr:Q2QG:cNS=-2.0,fsr_Q2QG_cNS_up fsr:Q2QG:cNS=2.0,fsr_X2XG_cNS_dn fsr:X2XG:cNS=-2.0,fsr_X2XG_cNS_up fsr:X2XG:cNS=2.0,isr_G2GG_muR_dn isr:G2GG:muRfac=0.5,isr_G2GG_muR_up isr:G2GG:muRfac=2.0,isr_G2QQ_muR_dn isr:G2QQ:muRfac=0.5,isr_G2QQ_muR_up isr:G2QQ:muRfac=2.0,isr_Q2QG_muR_dn isr:Q2QG:muRfac=0.5,isr_Q2QG_muR_up isr:Q2QG:muRfac=2.0,isr_X2XG_muR_dn isr:X2XG:muRfac=0.5,isr_X2XG_muR_up isr:X2XG:muRfac=2.0,isr_G2GG_cNS_dn isr:G2GG:cNS=-2.0,isr_G2GG_cNS_up isr:G2GG:cNS=2.0,isr_G2QQ_cNS_dn isr:G2QQ:cNS=-2.0,isr_G2QQ_cNS_up isr:G2QQ:cNS=2.0,isr_Q2QG_cNS_dn isr:Q2QG:cNS=-2.0,isr_Q2QG_cNS_up isr:Q2QG:cNS=2.0,isr_X2XG_cNS_dn isr:X2XG:cNS=-2.0,isr_X2XG_cNS_up isr:X2XG:cNS=2.0} and my hepmc file has the right number of events. Different from what I see when running our hadronizer.

Saptaparna commented 2 years ago

Very strange! I generated the file by using the test command from McM and then did this: Events->Scan("GenEventInfoProduct_generator__GEN.obj.weights_"). However, clicking through the branch, I can see the weights. So, at least the file was generated correctly. Mystery solved: GenEventInfoProduct_generator__SIM.obj.weights_ since I obviously ran SIM as the command was from McM.

agrohsje commented 2 years ago

You need to replace GEN with SIM. Your example is wmLHEGS.

Saptaparna commented 2 years ago

I know. I just did that. I am using this example to set up other diagnostic tests, hence the questions.

Saptaparna commented 2 years ago

One more question: how do I associate a weight with an "id"? Is there an easy way to do that? Otherwise, I'll write a quick edanalyzer.

SanghyunKo commented 2 years ago

Hi all, as @agrohsje already confirmed, I also see 46 weights in the HepMC file after adopting @smrenna's comment. The reason of having 91 weights in GenEventInfoProduct is that Pythia8.3 is directly adding PS weights to HepMC::GenEvent at the following:

event() = std::make_unique<HepMC::GenEvent>();
bool py8hepmc = toHepMC.fill_next_event(*(fMasterGen.get()), event().get());
// event()->weights().size() = 46

while in 8.2, the weights weren't added there so in Pythia8Hadronizer.cc we have an additional block to add PS weights manually:

// fill shower weights
// http://home.thep.lu.se/~torbjorn/pythia82html/Variations.html
if (fMasterGen->info.nWeights() > 1) {
  for (int i = 0; i < fMasterGen->info.nWeights(); ++i) {
    double wgt = fMasterGen->info.weight(i);
    event()->weights().push_back(wgt);
  }
}

and removing this section for Pythia8.3 gives 46 weights as intended.

The only remaining issue is that the baseline is not equal to the nominal. In my trial of Pythia8 DYToLL WeakSingleBoson:ffbar2gmZ (obviously not a CKKW, UNLOPS, .etc), both weights()[0] and fMasterGen.get()->info.mergingWeightNLO() give 1.0, but not for the weights()[1].

smrenna commented 2 years ago

Great work! I am on spring break but can look at this one when I am back.

Get Outlook for iOShttps://aka.ms/o0ukef

From: Sanghyun Ko @.> Sent: Wednesday, March 30, 2022 5:56:02 AM To: cms-sw/cmssw @.> Cc: Stephen Mrenna @.>; Mention @.> Subject: Re: [cms-sw/cmssw] 134 PS weights instead of 46 in Pythia 8.3 (Issue #36705)

Hi all, as @agrohsjehttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_agrohsje&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=P6TR6eMSzDEOLYeqW0ECsQ&m=4cT3681usRNZRcKGQCmPxFUnmCS3EU2bEufOI_DrtD6ntlw1_L5QXnRinte3DTqa&s=lqYAJq8oZcyPR83NHzIIDOYqnqGmwLUMuNLzADLa-kA&e= already confirmed, I also see 46 weights in the HepMC file after adopting @smrennahttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_smrenna&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=P6TR6eMSzDEOLYeqW0ECsQ&m=4cT3681usRNZRcKGQCmPxFUnmCS3EU2bEufOI_DrtD6ntlw1_L5QXnRinte3DTqa&s=RQglgOVH_C0oxj6BJyFW06rbKOAOU0V7aRO8QggvH0s&e='s comment. The reason of having 91 weights in GenEventInfoProduct is that Pythia8.3 is directly adding PS weights to HepMC::GenEvent at the following:

event() = std::make_unique(); bool py8hepmc = toHepMC.fill_next_event(*(fMasterGen.get()), event().get()); // event()->weights().size() = 46

while in 8.2, the weights weren't added there so in Pythia8Hadronizer.cc we have an additional block to add PS weights manually:

// fill shower weights // http://home.thep.lu.se/~torbjorn/pythia82html/Variations.htmlhttps://urldefense.proofpoint.com/v2/url?u=http-3A__home.thep.lu.se_-7Etorbjorn_pythia82html_Variations.html&d=DwQCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=P6TR6eMSzDEOLYeqW0ECsQ&m=4cT3681usRNZRcKGQCmPxFUnmCS3EU2bEufOI_DrtD6ntlw1_L5QXnRinte3DTqa&s=dJheY1gTxD6B-XABfvD-UeFYqe8z5eC96UcKyKyUz5c&e= if (fMasterGen->info.nWeights() > 1) { for (int i = 0; i < fMasterGen->info.nWeights(); ++i) { double wgt = fMasterGen->info.weight(i); event()->weights().push_back(wgt); } }

and removing this section for Pythia8.3 gives 46 weights as intended.

The only remaining issue is that the baseline is not equal to the nominal. In my trial of Pythia8 DYToLL WeakSingleBoson:ffbar2gmZ (obviously not a CKKW, UNLOPS, .etc), both weights()[0] and fMasterGen.get()->info.mergingWeightNLO() give 1.0, but not for the weights()[1].

— Reply to this email directly, view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_cms-2Dsw_cmssw_issues_36705-23issuecomment-2D1082878574&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=P6TR6eMSzDEOLYeqW0ECsQ&m=4cT3681usRNZRcKGQCmPxFUnmCS3EU2bEufOI_DrtD6ntlw1_L5QXnRinte3DTqa&s=8Vr9fgd8iNCZCn9D6e3ldt3gmVmh9WzenaRNZ79XbGg&e=, or unsubscribehttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_ABLSEJHQJLB55JDCYXFNQG3VCQQLFANCNFSM5L6RUXGQ&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=P6TR6eMSzDEOLYeqW0ECsQ&m=4cT3681usRNZRcKGQCmPxFUnmCS3EU2bEufOI_DrtD6ntlw1_L5QXnRinte3DTqa&s=b84zrqYzWFjltaxaWWUk4Yd1MFan6GHkrJTqNag8xfs&e=. You are receiving this because you were mentioned.Message ID: @.***>

SanghyunKo commented 2 years ago

@Saptaparna @GurpreetSinghChahal @menglu21 @smrenna Hi all, I just want to remind you guys that the deadline for the last open pre-release for 12_4_X (for Run 3 MC production) is approaching, and I believe we don't want to waste our storage for duplicated weights for O(10B) events.

Can we simply remove the code section that I mentioned above + fix the default weight list in the config file? and how should we proceed with the discrepancy between the nominal & baseline weights even when we have no merging feature?

Saptaparna commented 2 years ago

I think that should be fine as long as @smrenna agrees.

agrohsje commented 2 years ago

As time is indeed an issue we should clarify responsibilities. @SiewYan @mkirsano is one of you available to implement the fix as proposed by @SanghyunKo. let me know and i can help if you have no time. @SanghyunKo should we validate the weight impact against P8.240 once we have the fix ready?

SanghyunKo commented 2 years ago

@agrohsje The PR will be relatively simple, so I can create it right away (we just need a quick sign-off before the release). For the validation, fortunately we now have #36994 so the change should be already visible at the Relval stage. If we observe unexpected changes there (these should be the duplication of weights, not more than that), then we can launch a quick private validation with a simple process + PS weights.

SanghyunKo commented 2 years ago

FYI @shimashimarin

smrenna commented 2 years ago

The only remaining issue is that the baseline is not equal to the nominal. In my trial of Pythia8 DYToLL WeakSingleBoson:ffbar2gmZ (obviously not a CKKW, UNLOPS, .etc), both weights()[0] and fMasterGen.get()->info.mergingWeightNLO() give 1.0, but not for the weights()[1].

I don't understand what is meant by baseline and nominal.

SanghyunKo commented 2 years ago

Hi @smrenna, the 'nominal' weight here is the 0th weight calculated in the CMS code, which is merging weight * aMCatNLO norm factor (if there exists, otherwise it is 1.0), while 'baseline' is the 0th weight given by the Pythia8 (and this is 1st weight in the CMS code since we push back everything from Pythia8 after the nominal weight).

Pythia8::Info->mergingWeightNLO() * Pythia8::amcnlo_unitarised_interface->getNormFactor(); // nominal
Pythia8::Info->weight(0); // baseline

The thing is that the baseline weight from Pythia8.24x was always 1.0, but it seems that in 8.306 it's not 1.0, but some floating value fluctuating near 1.0. So I wonder what was the change for the 'baseline' (or 0th weight in Pythia) from 8.3.

smrenna commented 2 years ago

This is just rounding error coming from the multiplication of two 1.0s. In 8.2, the value of weight(0) was just stored internally:

double Info::weight(int iWeight) const {
   double wt = (iWeight <= 0 || iWeight >= int(weightSave.size()))
   ? weightSave[0] : weightSave[iWeight];
   return (abs(lhaStrategySave) == 4) ? CONVERTMB2PB * wt : wt;
}

In 8.3, there is an additional multiplication:

double Info::weight(int iWeight) const {
  double wt = weightContainerPtr->weightNominal;
  if (iWeight < 0 ||
    iWeight >= int(weightContainerPtr->weightsShowerPtr->getWeightsSize()))
     return wt;
  else
    return wt * weightContainerPtr->weightsShowerPtr->getWeightsValue(iWeight);
 }
SanghyunKo commented 2 years ago

ahhhh now I got this. There are two ways of retrieving weights from Pythia8. The first one is the one we used in the CMS with 8.24x (case a)

for (int i = 0; i < Pythia8::Pythia8->info.nWeights(); ++i) {
  double wgt = Pythia8::Pythia8->info.weight(i);
  HepMC::GenEvent()->weights().push_back(wgt);
  // std::cout << "a) name = " << Pythia8::Pythia8->info.weightLabel(i)
  //           <<  " value = " << Pythia8::Pythia8->info.weight(i) << std::endl;
}

and the second is the one newly added in Pythia8ToHepMC::fill_next_event (case b) in 8.30x (doing both a & b gave duplication)

for (int iweight = 0; iweight < Pythia8::Info->numberOfWeights(); ++iweight) {
  std::string name  = Pythia8::Info->weightNameByIndex(iweight);
  double value      = Pythia8::Info->weightValueByIndex(iweight);
  HepMC::GenEvent()->weights()[name] = value;
  // std::cout << "b) name = " << name << " value = " << value << std::endl;
}

and indeed both method should give exactly the same result (with slightly different weight names):

a) name = fsr:murfac=0.707 value = 1.028
a) name = Baseline value = 1.000
a) name = fsr:murfac=1.414 value = 0.970
a) name = fsr:murfac=0.5 value = 1.053
a) name = fsr:murfac=2.0 value = 0.942
a) name = fsr:murfac=0.25 value = 0.905
a) name = fsr:murfac=4.0 value = 0.900
a) name = fsr:g2gg:murfac=0.5 value = 1.168
a) name = fsr:g2gg:murfac=2.0 value = 0.895
a) name = fsr:g2qq:murfac=0.5 value = 0.993
a) name = fsr:g2qq:murfac=2.0 value = 1.004
a) name = fsr:q2qg:murfac=0.5 value = 0.909
a) name = fsr:q2qg:murfac=2.0 value = 1.048
a) name = fsr:x2xg:murfac=0.5 value = 1.000
a) name = fsr:x2xg:murfac=2.0 value = 1.000
a) name = fsr:g2gg:cns=-2.0 value = 1.000
a) name = fsr:g2gg:cns=2.0 value = 1.000
a) name = fsr:g2qq:cns=-2.0 value = 1.000
a) name = fsr:g2qq:cns=2.0 value = 1.000
a) name = fsr:q2qg:cns=-2.0 value = 1.000
a) name = fsr:q2qg:cns=2.0 value = 1.000
a) name = fsr:x2xg:cns=-2.0 value = 1.000
a) name = fsr:x2xg:cns=2.0 value = 1.000
a) name = isr:murfac=0.707 value = 1.150
a) name = isr:murfac=1.414 value = 0.860
a) name = isr:murfac=0.5 value = 1.258
a) name = isr:murfac=2.0 value = 0.745
a) name = isr:murfac=0.25 value = 1.421
a) name = isr:murfac=4.0 value = 0.554
a) name = isr:g2gg:murfac=0.5 value = 1.215
a) name = isr:g2gg:murfac=2.0 value = 0.822
a) name = isr:g2qq:murfac=0.5 value = 1.298
a) name = isr:g2qq:murfac=2.0 value = 0.794
a) name = isr:q2qg:murfac=0.5 value = 0.797
a) name = isr:q2qg:murfac=2.0 value = 1.142
a) name = isr:x2xg:murfac=0.5 value = 1.000
a) name = isr:x2xg:murfac=2.0 value = 1.000
a) name = isr:g2gg:cns=-2.0 value = 1.004
a) name = isr:g2gg:cns=2.0 value = 0.996
a) name = isr:g2qq:cns=-2.0 value = 1.000
a) name = isr:g2qq:cns=2.0 value = 1.000
a) name = isr:q2qg:cns=-2.0 value = 1.001
a) name = isr:q2qg:cns=2.0 value = 0.999
a) name = isr:x2xg:cns=-2.0 value = 1.000
a) name = isr:x2xg:cns=2.0 value = 1.000

b) name = Weight value = 1.000
b) name = AUX_fsr:murfac=0.707 value = 1.028
b) name = AUX_fsr:murfac=1.414 value = 0.970
b) name = AUX_fsr:murfac=0.5 value = 1.053
b) name = AUX_fsr:murfac=2.0 value = 0.942
b) name = AUX_fsr:murfac=0.25 value = 0.905
b) name = AUX_fsr:murfac=4.0 value = 0.900
b) name = AUX_fsr:g2gg:murfac=0.5 value = 1.168
b) name = AUX_fsr:g2gg:murfac=2.0 value = 0.895
b) name = AUX_fsr:g2qq:murfac=0.5 value = 0.993
b) name = AUX_fsr:g2qq:murfac=2.0 value = 1.004
b) name = AUX_fsr:q2qg:murfac=0.5 value = 0.909
b) name = AUX_fsr:q2qg:murfac=2.0 value = 1.048
b) name = AUX_fsr:x2xg:murfac=0.5 value = 1.000
b) name = AUX_fsr:x2xg:murfac=2.0 value = 1.000
b) name = AUX_fsr:g2gg:cns=-2.0 value = 1.000
b) name = AUX_fsr:g2gg:cns=2.0 value = 1.000
b) name = AUX_fsr:g2qq:cns=-2.0 value = 1.000
b) name = AUX_fsr:g2qq:cns=2.0 value = 1.000
b) name = AUX_fsr:q2qg:cns=-2.0 value = 1.000
b) name = AUX_fsr:q2qg:cns=2.0 value = 1.000
b) name = AUX_fsr:x2xg:cns=-2.0 value = 1.000
b) name = AUX_fsr:x2xg:cns=2.0 value = 1.000
b) name = AUX_isr:murfac=0.707 value = 1.150
b) name = AUX_isr:murfac=1.414 value = 0.860
b) name = AUX_isr:murfac=0.5 value = 1.258
b) name = AUX_isr:murfac=2.0 value = 0.745
b) name = AUX_isr:murfac=0.25 value = 1.421
b) name = AUX_isr:murfac=4.0 value = 0.554
b) name = AUX_isr:g2gg:murfac=0.5 value = 1.215
b) name = AUX_isr:g2gg:murfac=2.0 value = 0.822
b) name = AUX_isr:g2qq:murfac=0.5 value = 1.298
b) name = AUX_isr:g2qq:murfac=2.0 value = 0.794
b) name = AUX_isr:q2qg:murfac=0.5 value = 0.797
b) name = AUX_isr:q2qg:murfac=2.0 value = 1.142
b) name = AUX_isr:x2xg:murfac=0.5 value = 1.000
b) name = AUX_isr:x2xg:murfac=2.0 value = 1.000
b) name = AUX_isr:g2gg:cns=-2.0 value = 1.004
b) name = AUX_isr:g2gg:cns=2.0 value = 0.996
b) name = AUX_isr:g2qq:cns=-2.0 value = 1.000
b) name = AUX_isr:g2qq:cns=2.0 value = 1.000
b) name = AUX_isr:q2qg:cns=-2.0 value = 1.001
b) name = AUX_isr:q2qg:cns=2.0 value = 0.999
b) name = AUX_isr:x2xg:cns=-2.0 value = 1.000
b) name = AUX_isr:x2xg:cns=2.0 value = 1.000
b) name = AUX_cmsvar value = 0.547

However, the problem occurs when saving weights to HepMC: the former is a simple pushback, so it's straightforward. But in case of the latter, the container is std::map<std::string,double> - the problem is that the std::map is an ordered container so it'll sort them out by it's own rule std::string::operator<. So the order of (key,val) becomes

(AUX_cmsvar,0.547)
(AUX_fsr:g2gg:cns=-2.0,1.000)
(AUX_fsr:g2gg:cns=2.0,1.000)
(AUX_fsr:g2gg:murfac=0.5,1.168)
(AUX_fsr:g2gg:murfac=2.0,0.895)
(AUX_fsr:g2qq:cns=-2.0,1.000)
(AUX_fsr:g2qq:cns=2.0,1.000)
(AUX_fsr:g2qq:murfac=0.5,0.993)
(AUX_fsr:g2qq:murfac=2.0,1.004)
(AUX_fsr:murfac=0.25,0.905)
(AUX_fsr:murfac=0.5,1.053)
(AUX_fsr:murfac=0.707,1.028)
(AUX_fsr:murfac=1.414,0.970)
(AUX_fsr:murfac=2.0,0.942)
(AUX_fsr:murfac=4.0,0.900)
(AUX_fsr:q2qg:cns=-2.0,1.000)
(AUX_fsr:q2qg:cns=2.0,1.000)
(AUX_fsr:q2qg:murfac=0.5,0.909)
(AUX_fsr:q2qg:murfac=2.0,1.048)
(AUX_fsr:x2xg:cns=-2.0,1.000)
(AUX_fsr:x2xg:cns=2.0,1.000)
(AUX_fsr:x2xg:murfac=0.5,1.000)
(AUX_fsr:x2xg:murfac=2.0,1.000)
(AUX_isr:g2gg:cns=-2.0,1.004)
(AUX_isr:g2gg:cns=2.0,0.996)
(AUX_isr:g2gg:murfac=0.5,1.215)
(AUX_isr:g2gg:murfac=2.0,0.822)
(AUX_isr:g2qq:cns=-2.0,1.000)
(AUX_isr:g2qq:cns=2.0,1.000)
(AUX_isr:g2qq:murfac=0.5,1.298)
(AUX_isr:g2qq:murfac=2.0,0.794)
(AUX_isr:murfac=0.25,1.421)
(AUX_isr:murfac=0.5,1.258)
(AUX_isr:murfac=0.707,1.150)
(AUX_isr:murfac=1.414,0.860)
(AUX_isr:murfac=2.0,0.745)
(AUX_isr:murfac=4.0,0.554)
(AUX_isr:q2qg:cns=-2.0,1.001)
(AUX_isr:q2qg:cns=2.0,0.999)
(AUX_isr:q2qg:murfac=0.5,0.797)
(AUX_isr:q2qg:murfac=2.0,1.142)
(AUX_isr:x2xg:cns=-2.0,1.000)
(AUX_isr:x2xg:cns=2.0,1.000)
(AUX_isr:x2xg:murfac=0.5,1.000)
(AUX_isr:x2xg:murfac=2.0,1.000)
(Weight,1.000)

and this was the reason why the 1st element of the weight container was not 1.0.

So I suggest the following - to keep desired order of the HepMC::WeightContainer (as this becomes the weight order of CMS EDM as it is), we keep pursuing the method a) and turn off b) by Pythia8ToHepMC::set_store_weights(false) to prevent duplication of weights.

SanghyunKo commented 2 years ago

another aspect (I think this is the main reason) - Pythia8 didn't push any PS weights in 8.24x but only 0th weight (different from the baseline, CMS pushbacks the baseline & PS weights manually afterwards)

// Store cross-section information in pb and event weight. The latter is
// usually dimensionless, but in units of pb for Les Houches strategies +-4.
if (m_store_xsec && pyinfo != 0) {
  HepMC::GenCrossSection xsec;
  xsec.set_cross_section( pyinfo->sigmaGen() * 1e9,
    pyinfo->sigmaErr() * 1e9);
  evt->set_cross_section(xsec);
  evt->weights().push_back( pyinfo->weight() );
}

while in 8.30x it directly pushbacks PS weights without pushing Pythia8::Info->weight() beforehand

// Store cross-section information in pb and event weight. The latter is
// usually dimensionless, but in units of pb for Les Houches strategies +-4.
if (m_store_xsec && pyinfo != 0) {
  HepMC::GenCrossSection xsec;
  xsec.set_cross_section( pyinfo->sigmaGen() * 1e9,
    pyinfo->sigmaErr() * 1e9);
  evt->set_cross_section(xsec);
  // If multiweights with possibly different xsec, overwrite central value
  std::vector<double> xsecVec = pyinfo->weightContainerPtr->getTotalXsec();
  if (xsecVec.size() > 0) {
    xsec.set_cross_section(xsecVec[0]*1e9);
    evt->set_cross_section(xsec);
  }
}

if (m_store_weights && pyinfo != 0) {
  for (int iweight = 0; iweight < pyinfo->numberOfWeights();
    ++iweight) {
    std::string name  = pyinfo->weightNameByIndex(iweight);
    double value      = pyinfo->weightValueByIndex(iweight);
    evt->weights()[name] = value;
  }
}

so pushing Pythia8::Info->weight() manually could leave the weight container structure as same as before. This could be also related to the cms-sw/cmsdist#6688

SanghyunKo commented 2 years ago

since there is no follow-up anymore, I've just created #38012

SanghyunKo commented 2 years ago

Hi all, @Dongwoon77 had checked some basic ISR up/down & FSR up/down variables with DY01234JetsToLL MLM + Pythia8 PS weights between P8240 and P8306 with the patch I made.

https://dongwoon.web.cern.ch/p240vsp306/plots_comp/Generator_GenWeight.html

You can navigate to the leadJetPt & leadLepPt ISR & FSR up/down plots, whereas they represent leading jet & lepton Pt above 20 GeV weighted to corresponding PS weights. They have good agreements between the two versions after the fix.