[Rivet] Rivet analysis ATLAS_2010_S8914702

[Mike Hance]

Dear Holger,

Thanks for your mail!  Some comments inline below:

On 07/20/2016 02:10 AM, Holger Schulz wrote:
> Hi Michael,
>
> during testing we found that your analysis plugin crashes for certain events
> when trying to calculate the median pTDensity.
>
> Could I please ask you to have a look at these code snippets.
>
>
>       // Get the jet pT densities
>       vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
>       FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
>       const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = fastjets.clusterSeqArea();
>       for (const Jet& jet : fastjets.jets()) {
>         const double area = clust_seq_area->area(jet); //< Implicit call to pseudojet()
>         */// @todo Should be 1e-4?*
>         if (area > *10e-4* && jet.abseta() < _eta_bins_areaoffset.back()) {
>           ptDensities.at(getEtaBin(jet.abseta(), true)) += jet.pT()/area;
>         }
>       }

The above change should be fine, and I agree "10e-4" looks like a typo.  The area should always be 
greater than 0.001 or 0.0001 though, so I don't think this should cause any problems either way.

>
>       // Now compute the median energy densities
>       vector<double> ptDensity;
>       for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
>         *if (ptDensities[b].size() >0 )* ptDensity +=
> median(ptDensities[b]);
>       }
>
> The red bit was added by us to prevent calculating the median of an empty vector --- do you remember
> what the logic in your analysis code was for cases like that?

The red bit isn't safe.  The code is not super-elegant as written, but the "ptDensity" vector needs 
to have an entry for each "b" in the loop above.

One option is to change the code to be something like:

       // Now compute the median energy densities
       vector<double> ptDensity;
       for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
         if (ptDensities[b].size() >0 ) ptDensity += median(ptDensities[b]);
         else                           ptDensity += 0;
       }

Again, not super-elegant, but I don't have time (or a local setup) to test any changes, so I'm 
reluctant to do any extensive cleanup.  Let me know what you think.

Thanks!

-Mike


>
> If you find the time, could you please comment on the blue bit as well?
>
> You can find the version of the code in the current release here:
>
>     https://rivet.hepforge.org/hg/rivet/file/0175188f89df/src/Analyses/ATLAS_2010_S8914702.cc
>
>
> Cheers,
> Holger
>