---

HaloDepositor.cxx

Go to the documentation of this file.

#include "FastShowerUtils/HaloDepositor.h"
#include "FastShowerUtils/IDeposits.h"
#include "FastShowerUtils/IConfigurer.h"
#include "FastShowerUtils/ParticleParameters.h"
#include "FastShowerUtils/Normalisations.h"
#include "FastShowerUtils/IFnOfParticleParameters.h"
#include "FastShowerUtils/LinearProcessor.h"
#include "FastShowerUtils/TriangleProcessor.h"
#include "FastShowerUtils/ProcessedFlat.h"
#include <boost/array.hpp>
#include <cmath>
#include <numeric>
#include <vector>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <algorithm>
#include <strstream>

namespace FastShower{
  //                                               //
  // Function objects used by STL algorithms       //
  //                                               //
  class CalcPeak{

  public:
    CalcPeak(const ParticleParameters& pp, 
             HaloDepositor::HaloFnsIter hfIter ):m_pp(pp), m_hfIter(hfIter){}
    double operator()(){
      m_sp.setBoundary( (*m_hfIter++)->value(m_pp) );
      if(m_sp.lower()){
        LinearProcessor lp(0.01, 0.0);
        return m_pf.sample(&lp);
      }else{
        TriangleProcessor tp(1.0, 0.01);
        return m_pf.sample(&tp);
      }
    }
  private:
    const ParticleParameters& m_pp;
    SplitDecision m_sp;
    ProcessedFlat m_pf;
    HaloDepositor::HaloFnsIter m_hfIter;
  };
  class ApplyNorm{
  public:
    ApplyNorm(double norm):m_norm(norm){}
    void operator()(double& proxDep){ 
      proxDep *= m_norm;
   }
  private:
    double m_norm;
  };
  class GreaterThanOne{
  public:
    bool operator()(double arg){return arg>1.0;}
  };
  class HaloProxToGeom{
  public:
    
    HaloProxToGeom(HaloDepositor::HaloDepIter pdIter, IDeposits& deposits):
      m_proxIter(pdIter), m_deposits(deposits){}
    
    void operator()(std::pair<int, int> geom){
      //cout<<"HaloDepositor HaloProxToGeom: [eta,phi]: ["
      //    <<geom.first<<","<<geom.second<<"]"<<endl;
      m_deposits.accept(geom.first, geom.second, *m_proxIter++);
    }
    
  private:
    HaloDepositor::HaloDepIter m_proxIter;
    IDeposits& m_deposits;
  };
  class PrintFnValue{

  public:
    PrintFnValue(const ParticleParameters& pp):m_pp(pp){}
    void operator()(IFnOfParticleParameters* f){
      cout<<( f->value(m_pp) )<<" ";
    }
  private:
    const ParticleParameters& m_pp;
  };
  
  template<class T>
  class PrintValues{
  public:
    void operator()(T iter, const T& end){
      cout<<setiosflags(ios::fixed);
      for(;iter!=end; ++iter){cout<<setw(4)<<setprecision(3)<<*iter<<" ";}
      cout<<endl;
    }
  };
        
  //                                               //
  // End Function objects used by STL algorithms   //
  //                                               //
    
 
  HaloDepositor::HaloDepositor(const IConfigurer* c, const std::string& s):
    IDepositor(),DebugBase(s){

    cout<<"HaloDepositor: setting HaloOrder"<<endl;
    int order[14]={7, 6, 0, 8, 4, 1, 3, 13, 5, 10, 2, 11, 9, 12};
    for (std::size_t i=0; i<s_nHaloCells; ++i){
      m_order.push_back(order[i]);
    }
    assert(m_order.size() == s_nHaloCells);

    cout<<"HaloDepositor: getting HaloFns"<<endl;
    for (int i=1; i<=14; ++i){
      std::ostrstream no;
      no<<i<<ends;
      m_fns.push_back(c->findFnPP( text()+"Cell"+no.str() ));
    }
    assert(m_fns.size() == s_nHaloCells);

    cout<<"HaloDepositor: done with the configurer"<<endl;
    array<pair<int, int>, 14> q0, q1, q2, q3;
    // 
    q0[0] =pair<int, int>(  2,  0);
    q0[1] =pair<int, int>(  2,  1);
    q0[2] =pair<int, int>(  2,  2);
    q0[3] =pair<int, int>(  1,  2);
    q0[4] =pair<int, int>(  0,  2);
    q0[5] =pair<int, int>( -1,  2);
    q0[6] =pair<int, int>( -1,  1);
    q0[7] =pair<int, int>( -1,  0);
    q0[8] =pair<int, int>( -1, -1);
    q0[9] =pair<int, int>( -1, -2);
    q0[10]=pair<int, int>(  0, -2);
    q0[11]=pair<int, int>(  1, -2);
    q0[12]=pair<int, int>(  2, -2);
    q0[13]=pair<int, int>(  2, -1);
    assert(q0.size() == s_nHaloCells);
    m_mappers[0]=q0;
    // 
    q1[0] =pair<int, int>( -2,  0);
    q1[1] =pair<int, int>( -2,  1);
    q1[2] =pair<int, int>( -2,  2);
    q1[3] =pair<int, int>( -1,  2);
    q1[4] =pair<int, int>(  0,  2);
    q1[5] =pair<int, int>(  1,  2);
    q1[6] =pair<int, int>(  1,  1);
    q1[7] =pair<int, int>(  1,  0);
    q1[8] =pair<int, int>(  1, -1);
    q1[9] =pair<int, int>(  1, -2);
    q1[10]=pair<int, int>(  0, -2);
    q1[11]=pair<int, int>( -1, -2);
    q1[12]=pair<int, int>( -2, -2);
    q1[13]=pair<int, int>( -2, -1);
    assert(q1.size() == s_nHaloCells);
    m_mappers[1]=q1;
    //
    q2[0] =pair<int, int>( -2,  0);
    q2[1] =pair<int, int>( -2, -1);
    q2[2] =pair<int, int>( -2, -2);
    q2[3] =pair<int, int>( -1, -2);
    q2[4] =pair<int, int>(  0, -2);
    q2[5] =pair<int, int>(  1, -2);
    q2[6] =pair<int, int>(  1, -1);
    q2[7] =pair<int, int>(  1,  0);
    q2[8] =pair<int, int>(  1,  1);
    q2[9] =pair<int, int>(  1,  2);
    q2[10]=pair<int, int>(  0,  2);
    q2[11]=pair<int, int>( -1,  2);
    q2[12]=pair<int, int>( -2,  2);
    q2[13]=pair<int, int>( -2,  1);
    assert(q2.size() == s_nHaloCells);
    m_mappers[2]=q2;
    //
    q3[0] =pair<int, int>(  2,  0);
    q3[1] =pair<int, int>(  2, -1);
    q3[2] =pair<int, int>(  2, -2);
    q3[3] =pair<int, int>(  1, -2);
    q3[4] =pair<int, int>(  0, -2);
    q3[5] =pair<int, int>( -1, -2);
    q3[6] =pair<int, int>( -1, -1);
    q3[7] =pair<int, int>( -1,  0);
    q3[8] =pair<int, int>( -1,  1);
    q3[9] =pair<int, int>( -1,  2);
    q3[10]=pair<int, int>(  0,  2);
    q3[11]=pair<int, int>(  1,  2);
    q3[12]=pair<int, int>(  2,  2);
    q3[13]=pair<int, int>(  2,  1);
    assert(q3.size() == s_nHaloCells);
    m_mappers[3]=q3;
    //
  }  
  HaloDepositor::~HaloDepositor(){};
  void HaloDepositor::deposit(const ParticleParameters& pp,
                              const Normalisations& ns, 
                              IDeposits& geomDeposits) { 
    HaloDeposits proximDeps(s_nHaloCells);
    
    //Debug
    //PrintValues<vector<int>::iterator> ip;
    //PrintValues<vector<double>::iterator> dp;
    
    //cout<<"HaloDepositor: Shuffle order: "<<endl;
    //ip(m_order.begin(), m_order.end());

    //cout<<"HaloDepositor: Peak values: "<<endl;    
    //std::for_each(m_fns.begin(), m_fns.end(), PrintFnValue(pp));
    //cout<<endl;

    std::generate(proximDeps.begin(), 
                  proximDeps.end(), 
                  CalcPeak(pp, m_fns.begin() )
                  );


    //cout<<"HaloDepositor: ProxDepsafter Generation "<<endl;    
    //dp(proximDeps.begin(), proximDeps.end());
    //cout<<endl;

    double eDep = 0.0;
    eDep = accumulate(proximDeps.begin(), proximDeps.end(), eDep);

    //cout<<"HaloDepositor: eDep: "<<eDep<<endl;

    
    //    if(where != sumDeps.end()){
    if(eDep<=1.0){
      assert(eDep>0.001);
      HaloDepIter hdi = proximDeps.begin();
      for(; hdi != proximDeps.end(); *hdi++ *= 1/eDep ){}
      //cout<<"HaloDepositor: scale factor "<<1./eDep<<endl;
    }else{

      HaloDeposits shuffledDeps(s_nHaloCells);
      this->invShuffle(proximDeps, shuffledDeps);

      HaloDeposits sumDeps(s_nHaloCells);
      std::partial_sum(
                       shuffledDeps.begin(), 
                       shuffledDeps.end(), 
                       sumDeps.begin()
                       );

      //cout<<"HaloDepositor: SumDeps: "<<endl;
      //dp(sumDeps.begin(), sumDeps.end());
      
      HaloDepIter where = 
        std::find_if(sumDeps.begin(), sumDeps.end(), GreaterThanOne());
      
      HaloDepIter shdIter = shuffledDeps.begin()+(where-sumDeps.begin()); 
      
      *shdIter++ -= (*where)-1.; 
      for(; shdIter!=shuffledDeps.end(); *shdIter++=0.){}
      
      this->shuffle(shuffledDeps, proximDeps); // return to proximity scheme
    }
    
    //cout<<"HaloDepositor: Normalised Peaks: "<<endl;
    //dp(proximDeps.begin(), proximDeps.end());
      
    double depSum=0.;
    assert(abs(std::accumulate(proximDeps.begin(), 
                               proximDeps.end(), 
                               depSum)                
               )-1.0<0.001);
    
    //normalise
    std::for_each(proximDeps.begin(), proximDeps.end(), ApplyNorm(ns.halo()));
    
    // now map to geometery and send to deposits
    //cout<<"HaloDepositor mapping"<<endl;
    std::for_each(
                  m_mappers[pp.quadrant()].begin(),
                  m_mappers[pp.quadrant()].end(),
                  HaloProxToGeom( proximDeps.begin(), geomDeposits )
                  );
    //cout<<"HaloDepositor mapped"<<endl;
  }
  //
  void HaloDepositor::shuffle(const HaloDeposits& src, 
                              HaloDeposits& dest) const{
    HaloDepCIter sI;
    HaloArrCIter oI;
    for(sI = src.begin(), oI= m_order.begin(); sI != src.end(); ++sI, ++oI){
      dest[*oI] = *sI;
    }
  }
  //
  void HaloDepositor::invShuffle(const HaloDeposits& src, 
                                 HaloDeposits& dest) const{
    HaloDepIter  dI;
    HaloArrCIter oI;
    for(dI = dest.begin(), oI= m_order.begin(); dI != dest.end(); ++dI, ++oI){
      *dI = src[*oI];
    }
  }
  //
  IDepositor* HaloDepositor::clone() const{
    IDepositor* id = new HaloDepositor(*this);
    return id;
  }
  void HaloDepositor::components(IDebug::Cpts& v)const{
    v.clear();
    HaloFnsCIter itr;
    for (itr = m_fns.begin(); itr != m_fns.end(); ++itr) {
      v.push_back(*itr);
    }
  }
}//namespace

---