EnergyLossPlain Class Reference

#include <EnergyLossPlain.h>

Public Member Functions
	EnergyLossPlain (const TrackerGeometry *theTracker_, double pixelToStripMultiplier, double pixelToStripExponent)
int	estimate (const Trajectory *trajectory, std::vector< std::pair< int, double > > &arithmeticMean, std::vector< std::pair< int, double > > &truncatedMean)
void	loadOptimalWeights ()
	~EnergyLossPlain ()

Private Member Functions
double	average (std::vector< std::pair< double, double > > &values)
double	expected (double Delta1, double Delta2)
double	logTruncate (std::vector< std::pair< double, double > > &values)
double	optimal (std::vector< std::pair< double, double > > &values)
void	process (LocalVector ldir, const SiPixelRecHit *recHit, std::vector< std::pair< double, double > > &values)
void	process (LocalVector ldir, const SiStripRecHit2D *recHit, std::vector< std::pair< double, double > > &values)
double	truncate (std::vector< std::pair< double, double > > &values)

Private Attributes
double	pixelToStripExponent
double	pixelToStripMultiplier
const TrackerGeometry *	theTracker

Static Private Attributes
static bool	isFirst = true
static float	optimalWeight [61][61]

Detailed Description

Definition at line 12 of file EnergyLossPlain.h.

Constructor & Destructor Documentation

EnergyLossPlain::EnergyLossPlain	(	const TrackerGeometry *	theTracker_,
		double	pixelToStripMultiplier,
		double	pixelToStripExponent
	)

Definition at line 52 of file EnergyLossPlain.cc.

References cuy::isFirst.

                                   : theTracker(theTracker_),
      pixelToStripMultiplier(pixelToStripMultiplier_),
      pixelToStripExponent  (pixelToStripExponent_  )
{
  // Load data
  if(isFirst == true)
  {
    loadOptimalWeights();
    isFirst = false;
  }
}

EnergyLossPlain::~EnergyLossPlain

(

)

Definition at line 67 of file EnergyLossPlain.cc.

{
}

Member Function Documentation

double EnergyLossPlain::average ( std::vector< std::pair< double, double > > &  values )

private

Definition at line 96 of file EnergyLossPlain.cc.

References i, and makeHLTPrescaleTable::values.

{
  int num = values.size();
  double sum[2] = {0.,0.};

  for(int i = 0; i < num; i++)
  {
    sum[1] += values[i].first;
    sum[0] += 1;
  }

  return sum[1] / sum[0];
}

int EnergyLossPlain::estimate	(	const Trajectory *	trajectory,
		std::vector< std::pair< int, double > > &	arithmeticMean,
		std::vector< std::pair< int, double > > &	truncatedMean
	)

Definition at line 297 of file EnergyLossPlain.cc.

References a, PDRates::average, TrackingRecHit::geographicalId(), i, TrackingRecHit::isValid(), m, Trajectory::measurements(), ProjectedSiStripRecHit2D::originalHit(), GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, LaserDQM_cfg::process, alignCSCRings::s, findQualityFiles::v, hltHiggsPostProcessors_cff::vstr, and detailsBasic3DVector::y.

Referenced by EnergyLossProducer::produce().

{
  // (dE/dx, dx)
  std::vector<pair<double,double> > vpix,vstr; 

  for(std::vector<TrajectoryMeasurement>::const_iterator
        meas = trajectory->measurements().begin();
        meas!= trajectory->measurements().end(); meas++)
  {
    const TrackingRecHit* recHit = meas->recHit()->hit();
    DetId id = recHit->geographicalId();

    if(recHit->isValid())
    {
      LocalVector ldir = meas->updatedState().localDirection();
  
      if(theTracker->idToDet(id)->subDetector() ==
           GeomDetEnumerators::PixelBarrel ||
         theTracker->idToDet(id)->subDetector() ==
           GeomDetEnumerators::PixelEndcap)
      { 
        // Pixel
        const SiPixelRecHit* pixelRecHit =
          dynamic_cast<const SiPixelRecHit *>(recHit);
  
        if(pixelRecHit != 0)
          process(ldir,pixelRecHit , vpix);
      }
      else
      {
        // Strip
        const SiStripMatchedRecHit2D* stripMatchedRecHit =
          dynamic_cast<const SiStripMatchedRecHit2D *>(recHit);
        const ProjectedSiStripRecHit2D* stripProjectedRecHit =
          dynamic_cast<const ProjectedSiStripRecHit2D *>(recHit);
        const SiStripRecHit2D* stripRecHit =
          dynamic_cast<const SiStripRecHit2D *>(recHit);

        std::pair<double,double> v;
  
        if(stripMatchedRecHit != 0)
        {
          auto m = stripMatchedRecHit->monoHit();
          auto s = stripMatchedRecHit->stereoHit();
          process(ldir,&m, vstr);
          process(ldir,&s, vstr);
        } 

        if(stripProjectedRecHit != 0)
          process(ldir,&(stripProjectedRecHit->originalHit()), vstr);
  
        if(stripRecHit != 0)
          process(ldir,stripRecHit, vstr);
      }
    }
  }

  // Transform
  std::vector<pair<double,double> > vall;

  for(unsigned int i = 0; i < vpix.size(); i++)
  {
    float a = 0.795;
    float s = 10.1;

    std::pair<double,double> str(vpix[i]);

    double y = str.first / a / s;
    if(y > 0.9999) y =  0.9999;
    if(y <-0.9999) y = -0.9999;

    str.first = s * atanh(y);

    vall.push_back(str);
  }

  for(unsigned int i = 0; i < vstr.size(); i++) vall.push_back(vstr[i]);

  // Arithmetic mean
  arithmeticMean.push_back(std::pair<int,double>(vpix.size(), average(vpix)));
  arithmeticMean.push_back(std::pair<int,double>(vstr.size(), average(vstr)));
  arithmeticMean.push_back(std::pair<int,double>(vall.size(), average(vall)));

  // Wighted mean
  truncatedMean.push_back(std::pair<int,double>(vpix.size(), optimal(vpix)));
  truncatedMean.push_back(std::pair<int,double>(vstr.size(), optimal(vstr)));
  truncatedMean.push_back(std::pair<int,double>(vall.size(), optimal(vall)));

  return vall.size();
}

double EnergyLossPlain::expected ( double  Delta1, double  Delta2  )

private

Definition at line 205 of file EnergyLossPlain.cc.

References create_public_lumi_plots::log.

{
  return log(Delta2/Delta1) / (1/Delta1 - 1/Delta2);
//    return 1 + (Delta2*log(Delta1) - Delta1*log(Delta2)) / (Delta2 - Delta1);
}

void EnergyLossPlain::loadOptimalWeights

(

)

Definition at line 73 of file EnergyLossPlain.cc.

References edm::FileInPath::fullPath(), i, LogTrace, n, optimalWeight, and w().

{
  edm::FileInPath
    fileInPath("QCDAnalysis/ChargedHadronSpectra/data/energyWeights.dat");
  ifstream inFile(fileInPath.fullPath().c_str());

  while(inFile.eof() == false)
  {
    int i; float w; int n;
    inFile >> i;
    inFile >> w;
    inFile >> n;

    EnergyLossPlain::optimalWeight[n][i] = w;
  }

  inFile.close();

  LogTrace("MinBiasTracking")
    << " [EnergyLossEstimator] optimal weights loaded";
}

double EnergyLossPlain::logTruncate ( std::vector< std::pair< double, double > > &  values )

private

Definition at line 111 of file EnergyLossPlain.cc.

References i, funct::pow(), python.multivaluedict::sort(), findQualityFiles::v, makeHLTPrescaleTable::values, and histoStyle::weight.

{
  std::vector<double> values;
  for(std::vector<pair<double,double> >::iterator
      v = values_.begin(); v!= values_.end(); v++)
    values.push_back((*v).first);

  sort(values.begin(), values.end());

  int num = values.size();
  double sum[2] = {0.,1.};

  for(int i = 0; i < (num+1)/2; i++)
  { 
    double weight = 1.;

    if(num%2 == 1 && i == (num-1)/2) weight = 1./2;
  
    sum[1] *= pow(values[i], weight); 
    sum[0] += weight;

  }
  
  return pow(sum[1], 1./sum[0]);
} 

double EnergyLossPlain::optimal ( std::vector< std::pair< double, double > > &  values )

private

Definition at line 164 of file EnergyLossPlain.cc.

References first, i, create_public_lumi_plots::log, edm::second(), python.multivaluedict::sort(), findQualityFiles::v, makeHLTPrescaleTable::values, and histoStyle::weight.

{
  std::vector<double> values;
  for(std::vector<pair<double,double> >::iterator
      v = values_.begin(); v!= values_.end(); v++)
    values.push_back((*v).first);

  int num = values.size();
  sort(values.begin(), values.end());

  // First guess
  double sum = 0.;
  for(int i = 0; i < num; i++)
  {
    double weight = optimalWeight[num][i];
    sum += weight * values[i];
  }

  // Correct every deposit with log(path)
  for(int i = 0; i < num; i++)
    values_[i].first -= 0.178*log(values_[i].second/0.03) * 0.38 * sum;

  // Sort again 
  values.clear();
  for(std::vector<pair<double,double> >::iterator
      v = values_.begin(); v!= values_.end(); v++)
    values.push_back((*v).first);
  sort(values.begin(), values.end()); 

  // Corrected estimate
  sum = 0.;
  for(int i = 0; i < num; i++)
  {
    double weight = optimalWeight[num][i];
    sum += weight * values[i];
  }

  return sum;
}

void EnergyLossPlain::process ( LocalVector  ldir, const SiPixelRecHit *  recHit, std::vector< std::pair< double, double > > &  values  )

private

Definition at line 213 of file EnergyLossPlain.cc.

References Surface::bounds(), PV3DBase< T, PVType, FrameType >::mag(), PixelTopology::ncolumns(), PixelTopology::nrows(), PixelGeomDetUnit::specificTopology(), GeomDet::surface(), Bounds::thickness(), makeHLTPrescaleTable::values, x, and PV3DBase< T, PVType, FrameType >::z().

Referenced by ConfigBuilder.ConfigBuilder::addExtraStream(), ConfigBuilder.ConfigBuilder::completeInputCommand(), ConfigBuilder.ConfigBuilder::doNotInlineEventContent(), ConfigBuilder.ConfigBuilder.PrintAllModules::leave(), ConfigBuilder.ConfigBuilder::prepare_FASTSIM(), ConfigBuilder.ConfigBuilder::prepare_HLT(), ConfigBuilder.ConfigBuilder::prepare_LHE(), ConfigBuilder.ConfigBuilder::prepare_VALIDATION(), ConfigBuilder.ConfigBuilder::renameHLTprocessInSequence(), ConfigBuilder.ConfigBuilder::renameInputTagsInSequence(), and ConfigBuilder.ConfigBuilder::scheduleSequence().

{
  DetId id = recHit->geographicalId();
  const PixelGeomDetUnit* pixelDet =
    dynamic_cast<const PixelGeomDetUnit*> (theTracker->idToDet(id));

  // Check if cluster does not touch the edge
  if(recHit->cluster()->minPixelRow() == 0 ||
     recHit->cluster()->maxPixelRow() ==
       pixelDet->specificTopology().nrows() - 1 ||
     recHit->cluster()->minPixelCol() == 0 ||
     recHit->cluster()->maxPixelCol() == 
       pixelDet->specificTopology().ncolumns() - 1)
    return;

  double MeVperElec = 3.61e-6;

  // Collect adc
  double elec = recHit->cluster()->charge();
  double Delta = elec * MeVperElec;

  // Length
  double x = ldir.mag()/fabsf(ldir.z()) *
             pixelDet->surface().bounds().thickness();

  double pix = Delta/x;

  // MeV/cm, only if not low deposit
  if(pix > 1.5 * 0.795)
    values.push_back(std::pair<double,double>(pix, x)); 
}

void EnergyLossPlain::process ( LocalVector  ldir, const SiStripRecHit2D *  recHit, std::vector< std::pair< double, double > > &  values  )

private

Definition at line 247 of file EnergyLossPlain.cc.

References begin, Surface::bounds(), delta, end, i, PV3DBase< T, PVType, FrameType >::mag(), StripTopology::nstrips(), StripGeomDetUnit::specificTopology(), GeomDet::surface(), Bounds::thickness(), Tmax, makeHLTPrescaleTable::values, x, and PV3DBase< T, PVType, FrameType >::z().

Referenced by ConfigBuilder.ConfigBuilder::addExtraStream(), ConfigBuilder.ConfigBuilder::completeInputCommand(), ConfigBuilder.ConfigBuilder::doNotInlineEventContent(), ConfigBuilder.ConfigBuilder.PrintAllModules::leave(), ConfigBuilder.ConfigBuilder::prepare_FASTSIM(), ConfigBuilder.ConfigBuilder::prepare_HLT(), ConfigBuilder.ConfigBuilder::prepare_LHE(), ConfigBuilder.ConfigBuilder::prepare_VALIDATION(), ConfigBuilder.ConfigBuilder::renameHLTprocessInSequence(), ConfigBuilder.ConfigBuilder::renameInputTagsInSequence(), and ConfigBuilder.ConfigBuilder::scheduleSequence().

{
  DetId id = recHit->geographicalId();
  const StripGeomDetUnit* stripDet =
    dynamic_cast<const StripGeomDetUnit*> (theTracker->idToDet(id)); 

  // Check if cluster does not touch the edge
  if(recHit->cluster()->firstStrip() == 0 ||
     int(recHit->cluster()->firstStrip() +
         recHit->cluster()->amplitudes().size()) ==
       stripDet->specificTopology().nstrips())
    return;

  double MeVperAdc = 313 * 3.61e-6;

  // Collect adc
  double Delta = 0;
//  for(std::vector<uint16_t>::const_iterator
  for(std::vector<uint8_t>::const_iterator
    i = (recHit->cluster()->amplitudes()).begin();
    i!= (recHit->cluster()->amplitudes()).end(); i++)
  {
    // MeV
    double delta = (*i + 0.5) * MeVperAdc;

    if(*i >= 254)
    {
      double delta1, delta2; 
      if(*i == 254) { delta1 = 254 * MeVperAdc; delta2 = 512 * MeVperAdc; }
               else { delta1 = 512 * MeVperAdc; delta2 = Tmax; }

      delta = expected(delta1,delta2);
    }

    Delta += delta;
  }

  // Length
  double x = ldir.mag()/fabsf(ldir.z()) *
             stripDet->surface().bounds().thickness();

  double str = Delta/x;

  // MeV/cm, only if not low deposit
  if(str > 1.5)
    values.push_back(std::pair<double,double>(str, x));
}

double EnergyLossPlain::truncate ( std::vector< std::pair< double, double > > &  values )

private

Definition at line 138 of file EnergyLossPlain.cc.

References i, python.multivaluedict::sort(), findQualityFiles::v, makeHLTPrescaleTable::values, and histoStyle::weight.

{
  std::vector<double> values;
  for(std::vector<pair<double,double> >::iterator 
      v = values_.begin(); v!= values_.end(); v++)
    values.push_back((*v).first);

  sort(values.begin(), values.end());

  int num = values.size();
  double sum[2] = {0.,0.};

  for(int i = 0; i < (num+1)/2; i++)
  {
    double weight = 1.;

    if(num%2 == 1 && i == (num-1)/2) weight = 1./2;

    sum[1] += weight * values[i];
    sum[0] += weight;
  }

  return sum[1] / sum[0];
}

Member Data Documentation

bool EnergyLossPlain::isFirst = true

staticprivate

Definition at line 40 of file EnergyLossPlain.h.

float EnergyLossPlain::optimalWeight

staticprivate

Definition at line 41 of file EnergyLossPlain.h.

Referenced by loadOptimalWeights().

double EnergyLossPlain::pixelToStripExponent

private

Definition at line 38 of file EnergyLossPlain.h.

double EnergyLossPlain::pixelToStripMultiplier

private

Definition at line 38 of file EnergyLossPlain.h.

const TrackerGeometry* EnergyLossPlain::theTracker

private

Definition at line 37 of file EnergyLossPlain.h.