MinL3AlgoUnivErr.h

Go to the documentation of this file.

#ifndef MinL3AlgoUnivErr_H
#define MinL3AlgoUnivErr_H

//=============================================================================

//=============================================================================

#include <vector>
#include <iostream>
#include <map>
#include <math.h>

//=============================================================================
template <class IDdet>
class MinL3AlgoUnivErr
{
public:
  typedef          std::map<IDdet,float>  IDmapF;
  typedef typename IDmapF::value_type     IDmapFvalue;
  typedef typename IDmapF::iterator       iter_IDmapF;
  typedef typename IDmapF::const_iterator citer_IDmapF;
  typedef          std::map<IDdet,int>    IDmapI;
  typedef typename IDmapI::value_type     IDmapIvalue;
  typedef typename IDmapI::iterator       iter_IDmapI;
  typedef typename IDmapI::const_iterator citer_IDmapI;

  //----------------------------------------------

  MinL3AlgoUnivErr(float kweight_ = 0.);


  //----------------------------------------------

  ~MinL3AlgoUnivErr();


  //----------------------------------------------

  IDmapF iterate( const std::vector<std::vector<float> > & eventMatrix,
                  const std::vector<std::vector<IDdet> > & idMatrix,
                  const std::vector<float>               & energyVector,
                  const int                              & nIter,
                  const bool                             & normalizeFlag
                                                                       = false,
                  const int                              & nSubsamples = 10
                );
  //----------------------------------------------

  float getError( IDdet id ) const;

  //----------------------------------------------

  IDmapF getError() const;

  //----------------------------------------------

  int numberOfWrongErrors() const;

  //----------------------------------------------

  float  getErrorQuality( IDdet id ) const;

  IDmapF getErrorQuality()           const;

  //----------------------------------------------

  float getMeanPartialSolution( IDdet id ) const;

  //----------------------------------------------

  IDmapF getMeanPartialSolution() const;

  //----------------------------------------------

  void addEvent( const std::vector<float> & myCluster,
                 const std::vector<IDdet> & idCluster,
                 const float              & energy
               );

  //----------------------------------------------

  std::vector<float>
  recalibrateEvent( const std::vector<float> & myCluster,
                    const std::vector<IDdet> & idCluster,
                    const IDmapF             & newCalibration,
                    const int                & isol = -1,
                    const int                & iter = -1
                ); 

  //----------------------------------------------

  IDmapF getSolution( const bool resetsolution=true );

  //----------------------------------------------

  void resetSolution(); 

  //----------------------------------------------

private:

  float   kweight;
  int     countEvents;
  IDmapF  wsum;
  IDmapF  Ewsum;
  IDmapI                   idToIndex;    // map: cell id -> index of cell info
                                         // in sumPartSolu... vectors
  std::vector<int>         sumPartSolu0; // number of partial solutions
  std::vector<float>       sumPartSolu1; // sum of partial solutions
  std::vector<float>       sumPartSolu2; // sum of squared partial solutions

  int nOfSubsamples ; //  a copy of an input argument of `iterate' function

  //----------------------------------------------
  // register a partial solution in data members
  void registerPartialSolution( const IDmapF & partialSolution );

  //----------------------------------------------

};                                 //end of class MinL3AlgoUnivErr prototype


//=============================================================================


template<class IDdet>
MinL3AlgoUnivErr<IDdet>::MinL3AlgoUnivErr(float kweight_)
  :kweight(kweight_), countEvents(0)
{
  std::cout << "MinL3AlgoUnivErr : L3 algo with a calculation"
            << " of stat.errors working..." << std::endl;
  resetSolution();
}

//=============================================================================

template<class IDdet>
MinL3AlgoUnivErr<IDdet>::~MinL3AlgoUnivErr()
{
}

//=============================================================================

template<class IDdet>
typename MinL3AlgoUnivErr<IDdet>::IDmapF
MinL3AlgoUnivErr<IDdet>::
iterate( const std::vector < std::vector<float> > & eventMatrix,
         const std::vector < std::vector<IDdet> > & idMatrix,
         const std::vector <float>                & energyVector,
         const int                                & nIter,
         const bool                               & normalizeFlag,
         const int                                & nSubsamples
       )
{
  // clear the data members which are filled inside the function
  //                               (in registerPartialSolution called from here)

  nOfSubsamples = nSubsamples; // keep the input argument for use in other
                               // functions

  idToIndex    .clear();
  sumPartSolu0 .clear();
  sumPartSolu1 .clear();
  sumPartSolu2 .clear();

  IDmapF totalSolution;

  // Loop over samples/solutions:
  //    isol = 0                 : all events with the solution stored in
  //                               totalSolution
  //    isol = 1,...,nSubsamples : partial solutions are found for sub-samples
  //                               with the info on the solutions stored in the
  //                               data members
  //                                    idToIndex, sumPartSolu...
  //                               in order to be able to estimate the stat.
  //                               errors later

  for (int isol = 0 ; isol <= nSubsamples; isol++) {

    IDmapF sampleSolution ;  // solution for the sample
    IDmapF iterSolution   ;  // intermediate solution after an iteration
    std::vector < std::vector<float> > myEventMatrix  ;
    std::vector < std::vector<IDdet> > myIdMatrix     ;
    std::vector <float>                myEnergyVector ;

    // Select the sample.
    // Fill myEventMatrix, myIdMatrix and myEnergyVector 
    // either with all evs or with independent event subsamples

    if (isol == 0)  // total solution
      {
        myEventMatrix  = eventMatrix  ;
        myIdMatrix     = idMatrix     ;
        myEnergyVector = energyVector ;
      }
    else            // partial solution # isol
      {
        // clear containers filled for the previous sample
        sampleSolution .clear() ;
        myEventMatrix  .clear() ;
        myIdMatrix     .clear() ;
        myEnergyVector .clear() ;

        for (int i = 0; i < static_cast<int>( eventMatrix.size() ); i++)
          {
            // select every nSubsamples'th event to the subsample
            if ( i % nSubsamples + 1  ==  isol )
              {
                myEventMatrix  .push_back (eventMatrix  [i]) ;
                myIdMatrix     .push_back (idMatrix     [i]) ;
                myEnergyVector .push_back (energyVector [i]) ;
              }
          }
      }

    int Nevents = myEventMatrix.size(); // Number of events to calibrate with
    countEvents = 0;

    int i;

    // Iterate the correction
    for (int iter=1 ; iter <= nIter ; iter++) 
      {

        // if normalization flag is set, normalize energies
        float sumOverEnergy;
        if (normalizeFlag)
          {
            float scale = 0.;
            
            for (i=0; i<Nevents; i++)
              {
                sumOverEnergy = 0.;
                for (unsigned j=0 ; j < myEventMatrix[i].size() ; j++)
                                       {sumOverEnergy += myEventMatrix[i][j];}
                sumOverEnergy /= myEnergyVector[i];
                scale += sumOverEnergy;
              }
            scale /= Nevents;
          
            for (i=0; i<Nevents; i++) {myEnergyVector[i] *= scale;}       
          } // end normalize energies

        // now the real work starts:
        for (int iEvt=0; iEvt < Nevents; iEvt++)
          {
            addEvent( myEventMatrix[iEvt], myIdMatrix[iEvt]
                                         , myEnergyVector[iEvt] );
          }
        iterSolution = getSolution();
        if (iterSolution.empty())
          { sampleSolution.clear();
            break;    // exit the iteration loop leaving sampleSolution empty 
          }

        // re-calibrate eventMatrix with solution
        for (int ievent = 0; ievent<Nevents; ievent++)
          {
            myEventMatrix[ievent] = recalibrateEvent (myEventMatrix[ievent],
                                                      myIdMatrix[ievent],
                                                      iterSolution,
                                                      isol,iter);
          }

        // save solution into the sampleSolution map
        for (iter_IDmapF i = iterSolution.begin(); i != iterSolution.end(); i++)
          {
            iter_IDmapF itotal = sampleSolution.find(i->first);
            if (itotal == sampleSolution.end())
              {
                sampleSolution.insert(IDmapFvalue(i->first,i->second));
              }
            else
              {
                itotal->second *= i->second;
              }
          }

        //      resetSolution(); // reset for new iteration,
        //               now: getSolution does it automatically if not vetoed
      } // end iterate correction
    
    if (isol == 0) // total solution
      {
        totalSolution = sampleSolution;
      }
    else           // partial solution => register it in sumPartSolu...
      {
        registerPartialSolution( sampleSolution );
      }

  }  // end of the loop over solutions/samples

  return totalSolution;
}


//=============================================================================

template<class IDdet>
void MinL3AlgoUnivErr<IDdet>::
addEvent( const std::vector<float> & myCluster,
          const std::vector<IDdet> & idCluster,
          const float              & energy
        )
{
  countEvents++;

  float w, invsumXmatrix;
  float eventw;

  // Loop over the crystal matrix to find the sum
  float sumXmatrix=0.;
  for (unsigned i=0; i<myCluster.size(); i++) { sumXmatrix += myCluster[i]; }
      
  // event weighting
  eventw = 1 - fabs(1 - sumXmatrix/energy);
  eventw = pow(eventw,kweight);
      
  if (sumXmatrix != 0.)
    {
      invsumXmatrix = 1/sumXmatrix;
      // Loop over the crystal matrix (3x3,5x5,7x7) again
      // and calculate the weights for each xtal
      for (unsigned i=0; i<myCluster.size(); i++) 
        {               
          w = myCluster[i] * invsumXmatrix;

          // include the weights into wsum, Ewsum
          iter_IDmapF iwsum = wsum.find(idCluster[i]);
          if (iwsum == wsum.end())
                               wsum.insert(IDmapFvalue(idCluster[i],w*eventw));
          else iwsum->second += w*eventw;

          iter_IDmapF iEwsum = Ewsum.find(idCluster[i]);
          if (iEwsum == Ewsum.end())
                Ewsum.insert(IDmapFvalue (idCluster[i],
                                          (w*eventw * energy * invsumXmatrix)
                                         ));
          else iEwsum->second += (w*eventw * energy * invsumXmatrix);
        }
    }
  //  else {std::cout << " Debug: dropping null event: " << countEvents << std::endl;}
}

//=============================================================================

template<class IDdet>
typename MinL3AlgoUnivErr<IDdet>::IDmapF
MinL3AlgoUnivErr<IDdet>::
getSolution( const bool resetsolution )
{
  IDmapF solution;

  for (iter_IDmapF i = wsum.begin(); i != wsum.end(); i++)
    {
      iter_IDmapF iEwsum = Ewsum.find(i->first);
      float myValue = 1;
      if (i->second != 0) myValue = iEwsum->second / i->second;

      solution.insert( IDmapFvalue(i->first,myValue) );
    }
  
  if (resetsolution) resetSolution();

  return solution;
}

//=============================================================================

template<class IDdet>
void MinL3AlgoUnivErr<IDdet>::resetSolution()
{
  wsum.clear();
  Ewsum.clear();
}

//=============================================================================

template<class IDdet>
std::vector<float>
MinL3AlgoUnivErr<IDdet>::
recalibrateEvent( const std::vector <float> & myCluster,
                  const std::vector <IDdet> & idCluster,
                  const IDmapF              & newCalibration,
                  const int                 & isol,     // for a printout only
                  const int                 & iter      // for a printout only
                )
{
  std::vector<float> newCluster(myCluster);

  for (unsigned i=0; i<myCluster.size(); i++) 
    {
      citer_IDmapF icalib = newCalibration.find(idCluster[i]);
      if (icalib != newCalibration.end())
        {
          newCluster[i] *= icalib->second;
        }
      else
        {
          std::cout << "No calibration available for this element."
                    << std::endl;
          std::cout  << "   isol = " << isol
                     << "   iter = " << iter
                     << "   idCluster[i] = " << idCluster[i] << "\n";
        }

    }

  return newCluster;
}

//=============================================================================

template<class IDdet>
void
MinL3AlgoUnivErr<IDdet>::
registerPartialSolution( const IDmapF & partialSolution )

{
  int lastIndex = sumPartSolu0.size() - 1; // index of the last element
                                           // of the parallel vectors

  for (citer_IDmapF cell  = partialSolution.begin() ;
                    cell != partialSolution.  end() ; ++cell)
    {
      IDdet id   = cell->first ;
      float corr = cell->second;

      // where is the cell in another map?
      iter_IDmapI cell2 =  idToIndex.find( id );

      if ( cell2 == idToIndex.end() )
        {
          // the cell is met for the first time in patial solutions
          // => insert the info to the end of the vectors

          sumPartSolu0.push_back( 1 );
          sumPartSolu1.push_back( corr );
          sumPartSolu2.push_back( corr * corr );
          idToIndex.insert( IDmapIvalue( id, ++lastIndex ) );
        }
      else
        {
          // add the info to the already registered cell

          int index = cell2 -> second;
          sumPartSolu0[ index ] +=           1;
          sumPartSolu1[ index ] +=        corr;
          sumPartSolu2[ index ] += corr * corr;
        }

    }

}

//=============================================================================

template<class IDdet>
float
MinL3AlgoUnivErr<IDdet>::
getError( IDdet id ) const

{
  float error;
  citer_IDmapI cell =  idToIndex.find( id );
  if ( cell == idToIndex.end() ) error = -2.;     // no info for the cell
  else
    {
      int i = cell->second;
      int n = sumPartSolu0[ i ];
      if (n <= 1 ) error = -1.;     // 1 entry => error estimate impossible
      else
        {
          float meanX  = sumPartSolu1[ i ] / n;
          float meanX2 = sumPartSolu2[ i ] / n;

          error = sqrt( fabs(meanX2 - meanX * meanX) / (n - 1.)  ) ;
        }
        
    }
  return error;
}

//=============================================================================

template<class IDdet>
typename MinL3AlgoUnivErr<IDdet>::IDmapF
MinL3AlgoUnivErr<IDdet>::
getError() const

{
  IDmapF errors;
  float  error;

  for (citer_IDmapI cell  = idToIndex.begin();
                    cell != idToIndex.  end(); ++cell )
    {
      int i = cell->second;
      int n = sumPartSolu0[ i ];
      if (n <= 1 ) error = -1.;     // 1 entry => error estimate impossible
      else
        {
          float meanX  = sumPartSolu1[ i ] / n;
          float meanX2 = sumPartSolu2[ i ] / n;

          error = sqrt( fabs(meanX2 - meanX * meanX) / (n - 1)  ) ;
        }

      errors.insert( IDmapFvalue( cell->first , error ) );
    }
  return errors;
}
//=============================================================================

template<class IDdet>
float
MinL3AlgoUnivErr<IDdet>::
getErrorQuality( IDdet id ) const

{
  float fraction;
  citer_IDmapI cell =  idToIndex.find( id );
  if ( cell == idToIndex.end() ) fraction = -1.;   // no info for the cell
                                                   // => return a special value
  else
    {
      int i = cell->second;
      int n = sumPartSolu0[ i ];
      if (n < nOfSubsamples ) fraction = float( n ) / nOfSubsamples; 
      else                    fraction = 1.;
    }
  return fraction;
}

//=============================================================================

template<class IDdet>
typename MinL3AlgoUnivErr<IDdet>::IDmapF
MinL3AlgoUnivErr<IDdet>::
getErrorQuality()           const

{
  IDmapF fractions;
  float  fraction;

  for (citer_IDmapI cell  = idToIndex.begin();
                    cell != idToIndex.  end(); ++cell )
    {
      int i = cell->second;
      int n = sumPartSolu0[ i ];

      if (n < nOfSubsamples ) fraction = float( n ) / nOfSubsamples; 
      else                    fraction = 1.;

      fractions.insert( IDmapFvalue( cell->first , fraction ) );
    }

  return fractions;
}


//=============================================================================

template<class IDdet>
int
MinL3AlgoUnivErr<IDdet>::
numberOfWrongErrors() const

{
  int nWrong = 0 ;

  for (citer_IDmapI cell  = idToIndex.begin();
                    cell != idToIndex.  end(); ++cell )
    {
      int i = cell->second;
      int n = sumPartSolu0[ i ];

      if (n < nOfSubsamples ) nWrong++;
    }

  return nWrong;
}

//=============================================================================

template<class IDdet>
float
MinL3AlgoUnivErr<IDdet>::
getMeanPartialSolution( IDdet id ) const

{
  float meanX;
  citer_IDmapI cell =  idToIndex.find( id );
  if ( cell == idToIndex.end() ) meanX = 999.;     // no info for the cell
  else
    {
      int i = cell->second;
      int   n      = sumPartSolu0[ i ];
      float meanX  = sumPartSolu1[ i ] / n;
    }
  return meanX;

}

//=============================================================================

template<class IDdet>
typename MinL3AlgoUnivErr<IDdet>::IDmapF
MinL3AlgoUnivErr<IDdet>::
getMeanPartialSolution () const

{
  IDmapF solution;

  for (citer_IDmapI cell  = idToIndex.begin();
                    cell != idToIndex.  end(); ++cell )
    {
      int i = cell->second;
      int   n      = sumPartSolu0[ i ];
      float meanX  = sumPartSolu1[ i ] / n;
      solution.insert( IDmapFvalue( cell->first , meanX ) );
    }
  return solution;
}

//=============================================================================

#endif // MinL3AlgoUnivErr_H