[Rivet] Rivet analysis ATLAS_2010_S8914702

[Andy Buckley]

On 20/07/16 17:24, Mike Hance wrote:
> 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.

Thanks -- there were about 6 analyses with the same copy & paste typo, 
so this has really allowed me to clean up the analysis code repo!


>>       // 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 for explaining. And here's the elegant way (IMHO):

ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);

Not bad! I've made this change in my 25x branch, Holger; it'll be pushed 
soon.

Andy

-- 
Dr Andy Buckley, Lecturer / Royal Society University Research Fellow
Particle Physics Expt Group, University of Glasgow