#include <HPDNoiseLibraryReader.h>

Public Member Functions
std::vector< std::pair < HcalDetId, const float * > >	getBiasedNoisyHcalDetIds ()
std::vector< std::pair < HcalDetId, const float * > >	getBiasedNoisyHcalDetIds (int timeSliceId)
double	getIonFeedbackNoise (HcalDetId id, double energy, double bias)
std::vector< std::pair < HcalDetId, const float * > >	getNoisyHcalDetIds (int timeSliceId)
std::vector< std::pair < HcalDetId, const float * > >	getNoisyHcalDetIds ()
	HPDNoiseLibraryReader (const edm::ParameterSet &)
	~HPDNoiseLibraryReader ()
Static Public Member Functions
static void	initializeServices ()
Public Attributes
std::vector< float >	theDischargeNoiseRate
std::string	theHPDName
std::vector< float >	theIonFeedbackFirstPeakRate
std::vector< float >	theIonFeedbackSecondPeakRate
std::vector< std::string >	theNames
std::vector< int >	theNoisyPhi
CLHEP::RandFlat *	theRandFlat
CLHEP::RandGaussQ *	theRandGaussQ
HPDNoiseReader *	theReader
Protected Member Functions
void	setRandomEngine ()
void	setRandomEngine (CLHEP::HepRandomEngine &engine)
Private Member Functions
void	clearPhi ()
void	fillRates ()
void	getBiasedNoisyPhis ()
HPDNoiseData *	getNoiseData (int iphi)
void	getNoisyPhis ()
bool	IsNoiseApplicable (int iphi)
std::string	itos (int i)
void	Rannor (double &a, double &b)
void	shuffleData (int timeSliceId, float *&data)

Detailed Description

Definition at line 44 of file HPDNoiseLibraryReader.h.

Constructor & Destructor Documentation

HPDNoiseLibraryReader::HPDNoiseLibraryReader ( const edm::ParameterSet & iConfig )

Definition at line 22 of file HPDNoiseLibraryReader.cc.

References HPDNoiseReader::allNames(), fillRates(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), setRandomEngine(), theHPDName, theNames, and theReader.

:theDischargeNoiseRate(0),
theIonFeedbackFirstPeakRate(0),
theIonFeedbackSecondPeakRate(0),
theNoisyPhi(0),
theRandFlat(0),
theRandGaussQ(0) {
    setRandomEngine();

    ParameterSet pSet = iConfig.getParameter < edm::ParameterSet > ("HPDNoiseLibrary");
    FileInPath filepath = pSet.getParameter < edm::FileInPath > ("FileName");

    theHPDName = pSet.getUntrackedParameter < string > ("HPDName", "HPD");
    string pName = filepath.fullPath();

    if (pName.find(".") == 0)
        pName.erase(0, 2);
    theReader = new HPDNoiseReader(pName);
    theNames = theReader->allNames();   // all 72x2 HPDs

    fillRates();
}

HPDNoiseLibraryReader::~HPDNoiseLibraryReader ( )

Definition at line 45 of file HPDNoiseLibraryReader.cc.

References theRandFlat, and theRandGaussQ.

                                              {
    if (theRandFlat)
        delete theRandFlat;

    if (theRandGaussQ)
        delete theRandGaussQ;
}

Member Function Documentation

void HPDNoiseLibraryReader::clearPhi ( ) [private]

Definition at line 393 of file HPDNoiseLibraryReader.cc.

References theNoisyPhi.

Referenced by getBiasedNoisyPhis(), and getNoisyPhis().

                                     {
    theNoisyPhi.clear();
}

void HPDNoiseLibraryReader::fillRates ( ) [private]

Definition at line 88 of file HPDNoiseLibraryReader.cc.

References dtNoiseDBValidation_cfg::cerr, HPDNoiseReader::dischargeRate(), HPDNoiseReader::getHandle(), i, HPDNoiseReader::ionFeedbackFirstPeakRate(), HPDNoiseReader::ionFeedbackSecondPeakRate(), theDischargeNoiseRate, theIonFeedbackFirstPeakRate, theIonFeedbackSecondPeakRate, theNames, theReader, and HPDNoiseReader::valid().

Referenced by HPDNoiseLibraryReader().

                                      {
    for (size_t i = 0; i < theNames.size(); ++i) {
        HPDNoiseReader::Handle hpdObj = theReader->getHandle(theNames[i]);
        if (theReader->valid(hpdObj)) {
            theDischargeNoiseRate.push_back(theReader->dischargeRate(hpdObj));
            theIonFeedbackFirstPeakRate.push_back(theReader->ionFeedbackFirstPeakRate(hpdObj));
            theIonFeedbackSecondPeakRate.push_back(theReader->ionFeedbackSecondPeakRate(hpdObj));
        } else {
            std::cerr << "HPD Handle Object is not valid!" << endl;
        }
    }
}

vector< pair< HcalDetId, const float * > > HPDNoiseLibraryReader::getBiasedNoisyHcalDetIds ( )

Definition at line 236 of file HPDNoiseLibraryReader.cc.

References data, getBiasedNoisyPhis(), HPDNoiseData::getDataFrame(), HPDNoiseDataFrame::getFrame(), getNoiseData(), i, HPDNoiseDataFrame::id(), createXMLFile::iphi, query::result, HPDNoiseData::size(), and theNoisyPhi.

                                                                                           {

    vector < pair < HcalDetId, const float *> >result;

    // decide which phi are noisy by using noise rates 
    // and random numbers (to be called for each event)
    // at least one Phi is always noisy.
    getBiasedNoisyPhis();
    for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
        int iphi = theNoisyPhi[i];
        HPDNoiseData *data;

        data = getNoiseData(iphi);
        for (unsigned int i = 0; i < data->size(); ++i) {
            result.emplace_back(data->getDataFrame(i).id(), data->getDataFrame(i).getFrame());
        }
    }
    return result;
}

vector< pair< HcalDetId, const float * > > HPDNoiseLibraryReader::getBiasedNoisyHcalDetIds ( int timeSliceId )

Definition at line 213 of file HPDNoiseLibraryReader.cc.

References data, getBiasedNoisyPhis(), HPDNoiseData::getDataFrame(), HPDNoiseDataFrame::getFrame(), getNoiseData(), i, HPDNoiseDataFrame::id(), createXMLFile::iphi, query::result, shuffleData(), HPDNoiseData::size(), and theNoisyPhi.

                                                                                                          {

    vector < pair < HcalDetId, const float *> >result;

    // decide which phi are noisy by using noise rates 
    // and random numbers (to be called for each event)
    // at least one Phi is always noisy.
    getBiasedNoisyPhis();
    for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
        int iphi = theNoisyPhi[i];
        HPDNoiseData *data;

        data = getNoiseData(iphi);
        for (unsigned int i = 0; i < data->size(); ++i) {
            float* data_ = const_cast<float*>(data->getDataFrame(i).getFrame());
            shuffleData(timeSliceId, data_);
            const float* _data_ =const_cast<const float*>(data_);
            result.emplace_back(data->getDataFrame(i).id(), _data_);
        }
    }
    return result;
}

void HPDNoiseLibraryReader::getBiasedNoisyPhis ( ) [private]

Definition at line 153 of file HPDNoiseLibraryReader.cc.

References clearPhi(), i, theDischargeNoiseRate, theNoisyPhi, and theRandFlat.

Referenced by getBiasedNoisyHcalDetIds().

                                               {

    clearPhi();
    double rndm[144];

    theRandFlat->shootArray(144, rndm);
    for (int i = 0; i < 144; ++i) {
        if (rndm[i] < theDischargeNoiseRate[i]) {
            theNoisyPhi.push_back(i + 1);
        }
    }
    // make sure one HPD is always noisy
    if (theNoisyPhi.size() == 0) {
        int iPhi = (theRandFlat->fireInt(144)) + 1; // integer from interval [0-144[ + 1

        theNoisyPhi.push_back(iPhi);
    }
}

double HPDNoiseLibraryReader::getIonFeedbackNoise	(	HcalDetId	id,
		double	energy,
		double	bias
	)

HPD Ion feedback simulation based on LED data. A simple simulation which uses gaussian fit to data. biased = false ==> HPD noise from Ion Feedback only, unbiased biased = true ==> HPD noise from Ion Feedback only, biased (rate is X times larger than nominal rate)

Definition at line 256 of file HPDNoiseLibraryReader.cc.

References a, b, Exception, createXMLFile::iphi, j, p4, pe2Charge, Rannor(), theIonFeedbackFirstPeakRate, theIonFeedbackSecondPeakRate, theRandFlat, and theRandGaussQ.

                                                                                          {

    // constants for simulation/parameterization
    double pe2Charge = 0.333333;    // fC/p.e.
    double GeVperfC = 0.177;    // in unit GeV/fC and this will be updated when it start reading from DB.
    double PedSigma = 0.8;
    double noise = 0.;          // fC

    int iphi = (id.ieta() > 0) ? (id.iphi()) : (id.iphi() + 72);
    double rateInTail = theIonFeedbackFirstPeakRate[iphi - 1];
    double rateInSecondTail = theIonFeedbackSecondPeakRate[iphi - 1];

    if (bias != 0.) {
        rateInTail = rateInTail * bias;
        rateInSecondTail = rateInSecondTail * bias;
    } else {
        edm::LogError("HPDNoise") << "HPDNoise: ion feedback error (biased or unbiased selection)." << bias << " failed";
        throw cms::Exception("Unknown", "biase selection ")
        << "Usage of " << bias << " fails\n";
    }
    double Charge = energy / GeVperfC;

    // three gauss fit is applied to data to get ion feedback distribution
    // the script is at neutralino: /home/tyetkin/work/hpd_noise/PlotFromPelin.C
    // a new fit woth double-sigmoids is under way.
    // parameters (in fC)
    // first gaussian
    // double p0 = 9.53192e+05;
    // double p1 = -3.13653e-01;
    // double p2 = 2.78350e+00;

    // second gaussian
    // double p3 = 2.41611e+03;
    double p4 = 2.06117e+01;
    double p5 = 1.09239e+01;

    // third gaussian
    // double p6 = 3.42793e+01;
    double p7 = 5.45548e+01;
    double p8 = 1.59696e+01;

    if (Charge > 3 * PedSigma) {    // 3 sigma away from pedestal mean
        int npe = int (Charge / pe2Charge);
        double a = 0.;
        double b = 0.;

        for (int j = 0; j < npe; ++j) {
            double probability = theRandFlat->shoot();

            if (probability < rateInTail) { // total tail
                if (probability < rateInSecondTail) {   // second tail
                    Rannor(a, b);
                    noise += b * p8 + p7;
                } else {
                    Rannor(a, b);
                    noise += b * p5 + p4;
                }
            }
        }
        // add pedestal 
        if (noise > 0.)
            noise += theRandGaussQ->fire(0, 2 * PedSigma);
    }
    return (noise * GeVperfC);  // returns noise in GeV.

}

HPDNoiseData * HPDNoiseLibraryReader::getNoiseData ( int iphi ) [private]

Definition at line 101 of file HPDNoiseLibraryReader.cc.

References dtNoiseDBValidation_cfg::cerr, data, HPDNoiseReader::getEntry(), HPDNoiseReader::getHandle(), IsNoiseApplicable(), itos(), mergeVDriftHistosByStation::name, AlCaHLTBitMon_QueryRunRegistry::string, theHPDName, theRandFlat, theReader, HPDNoiseReader::totalEntries(), and HPDNoiseReader::valid().

Referenced by getBiasedNoisyHcalDetIds(), and getNoisyHcalDetIds().

                                                          {


    HPDNoiseData *data = 0; //data->size() is checked wherever actually used  

    // make sure that iphi from HcalDetId is found noisy at first.
    // In other words, be sure that iphi is in the collection of
    // noisy Phis
    if (!(IsNoiseApplicable(iphi)))         
        return data;

    int zside = 1;

    if (iphi > 72) {
        iphi = iphi - 72;
        zside = -1;
    }
    std::string name;
    if (zside == 1) {
        name = "ZPlus" + theHPDName + itos(iphi);
    } else if (zside == -1) {
        name = "ZMinus" + theHPDName + itos(iphi);
    } else {
        cerr << " ZSide Calculation Error." << endl;
    }
    HPDNoiseReader::Handle hpdObj = theReader->getHandle(name);
    if (theReader->valid(hpdObj)) {
        // randomly select one entry from library for this HPD
        unsigned long entry = theRandFlat->fireInt(theReader->totalEntries(hpdObj));

        theReader->getEntry(hpdObj, entry, &data);
    } else {
        std::cerr << " HPD Name in the library is not valid." << std::endl;
    }
    return data;
}

vector< pair< HcalDetId, const float * > > HPDNoiseLibraryReader::getNoisyHcalDetIds ( )

Definition at line 172 of file HPDNoiseLibraryReader.cc.

References data, HPDNoiseData::getDataFrame(), HPDNoiseDataFrame::getFrame(), getNoiseData(), getNoisyPhis(), i, HPDNoiseDataFrame::id(), createXMLFile::iphi, query::result, HPDNoiseData::size(), and theNoisyPhi.

Referenced by HPDNoiseGenerator::fillNoiseSignals().

                                                                                  {

    vector <pair< HcalDetId, const float *> >result;

    // decide which phi are noisy by using noise rates 
    // and random numbers (to be called for each event)
    getNoisyPhis();
    for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
        int iphi = theNoisyPhi[i];
        HPDNoiseData *data;

        data = getNoiseData(iphi);
        for (unsigned int i = 0; i < data->size(); ++i) {

            result.emplace_back(data->getDataFrame(i).id(), data->getDataFrame(i).getFrame());
        }
    }
    return result;
}

vector< pair< HcalDetId, const float * > > HPDNoiseLibraryReader::getNoisyHcalDetIds ( int timeSliceId )

Definition at line 192 of file HPDNoiseLibraryReader.cc.

References data, HPDNoiseData::getDataFrame(), HPDNoiseDataFrame::getFrame(), getNoiseData(), getNoisyPhis(), i, HPDNoiseDataFrame::id(), createXMLFile::iphi, query::result, shuffleData(), HPDNoiseData::size(), and theNoisyPhi.

{
    vector <pair< HcalDetId, const float *> >result;
    // decide which phi are noisy by using noise rates 
    // and random numbers (to be called for each event)
    getNoisyPhis();
    for (int i = 0; i < int (theNoisyPhi.size()); ++i) {
        int iphi = theNoisyPhi[i];
        HPDNoiseData *data;

        data = getNoiseData(iphi);
        for (unsigned int i = 0; i < data->size(); ++i) {
            float* data_ = const_cast<float*>(data->getDataFrame(i).getFrame());
            shuffleData(timeSliceId, data_);
            const float* _data_ =const_cast<const float*>(data_);
            result.emplace_back(data->getDataFrame(i).id(), _data_);
        }
    }
    return result;

}

void HPDNoiseLibraryReader::getNoisyPhis ( ) [private]

Definition at line 139 of file HPDNoiseLibraryReader.cc.

References clearPhi(), i, theDischargeNoiseRate, theNoisyPhi, and theRandFlat.

Referenced by getNoisyHcalDetIds().

                                         {

    clearPhi();
    double rndm[144];

    theRandFlat->shootArray(144, rndm);

    for (int i = 0; i < 144; ++i) {
        if (rndm[i] < theDischargeNoiseRate[i]) {
            theNoisyPhi.push_back(i + 1);
        }
    }
}

void HPDNoiseLibraryReader::initializeServices ( ) [static]

Definition at line 53 of file HPDNoiseLibraryReader.cc.

References HDQMDatabaseProducer::config, edmplugin::standard::config(), edmplugin::PluginManager::configure(), edm::ServiceRegistry::createServicesFromConfig(), edmplugin::PluginManager::isAvailable(), and AlCaHLTBitMon_QueryRunRegistry::string.

                                               {
    if (not edmplugin::PluginManager::isAvailable()) {
        edmplugin::PluginManager::configure(edmplugin::standard::config());
    }

    std::string config =
      "process CorrNoise = {"
      "service = RandomNumberGeneratorService" "{" "untracked uint32 sourceSeed = 123456789" "}" "}";

    // create the services
    edm::ServiceToken tempToken = edm::ServiceRegistry::createServicesFromConfig(config);

    // make the services available
    edm::ServiceRegistry::Operate operate(tempToken);
}

bool HPDNoiseLibraryReader::IsNoiseApplicable ( int iphi ) [private]

Definition at line 323 of file HPDNoiseLibraryReader.cc.

References spr::find(), createXMLFile::iphi, and theNoisyPhi.

Referenced by getNoiseData().

                                                      {

    bool isAccepted = false;
    vector < int >::iterator phi_iter;

    phi_iter = find(theNoisyPhi.begin(), theNoisyPhi.end(), iphi);
    if (phi_iter != theNoisyPhi.end()) {
        isAccepted = true;
    }
    return isAccepted;
}

string HPDNoiseLibraryReader::itos ( int i ) [private]

Definition at line 386 of file HPDNoiseLibraryReader.cc.

References i, and alignCSCRings::s.

Referenced by getNoiseData().

                                        {
    stringstream s;

    s << i;
    return s.str();
}

void HPDNoiseLibraryReader::Rannor	(	double &	a,
		double &	b
	)		`[private]`

Definition at line 373 of file HPDNoiseLibraryReader.cc.

References alignCSCRings::r, theRandFlat, x, detailsBasic3DVector::y, and z.

Referenced by getIonFeedbackNoise().

                                                       {
    double r,
      x,
      y,
      z;

    y = theRandFlat->shoot();
    z = theRandFlat->shoot();
    x = z * 6.28318530717958623;
    r = TMath::Sqrt(-2 * TMath::Log(y));
    a = r * TMath::Sin(x);
    b = r * TMath::Cos(x);
}

void HPDNoiseLibraryReader::setRandomEngine ( ) [protected]

Definition at line 68 of file HPDNoiseLibraryReader.cc.

References Exception, and edm::Service< T >::isAvailable().

Referenced by HPDNoiseLibraryReader().

                                            {
    edm::Service < edm::RandomNumberGenerator > rng;
    if (!rng.isAvailable()) {
        throw cms::Exception("Configuration") << "HcalHPDNoiseLibrary requires the RandomNumberGeneratorService\n"
          "which is not present in the configuration file.  You must add the service\n"
          "in the configuration file or remove the modules that require it.";
    }
    setRandomEngine(rng->getEngine());
}

void HPDNoiseLibraryReader::setRandomEngine ( CLHEP::HepRandomEngine & engine ) [protected]

Definition at line 77 of file HPDNoiseLibraryReader.cc.

References theRandFlat, and theRandGaussQ.

                                                                         {
    if (theRandGaussQ)
        delete theRandGaussQ;

    if (theRandFlat)
        delete theRandFlat;

    theRandGaussQ = new CLHEP::RandGaussQ(engine);
    theRandFlat = new CLHEP::RandFlat(engine);
}

void HPDNoiseLibraryReader::shuffleData	(	int	timeSliceId,
		float *&	data
	)		`[private]`

Definition at line 334 of file HPDNoiseLibraryReader.cc.

References createTree::dd, i, gen::k, and max().

Referenced by getBiasedNoisyHcalDetIds(), and getNoisyHcalDetIds().

{
   if(timeSliceId == -1 || (timeSliceId>=10)) return;
   //make a local copy of input data
   float Data[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
   for(int i=0;i<10;++i){
       Data[i] = data[i];
   }
   int ts_max = -1;
   float max = -999.;
   for(int i=0;i<10;++i){
       if(Data[i]>max){
           max = data[i];
           ts_max = i;
       }
   }
   if((ts_max == -1)){//couldn't find ts_max, return the same value.
       return;
   }else{
       // always shift the noise to the right by putting zeroes to the previous slices.
       // the noise is pedestal subtracted. 0 value is acceptable.
       int k = -1;
       for(int i=0;i<10;++i){
           data[i] = 0.;
           int newIdx = timeSliceId+k;
           float dd = 0.;
           if(newIdx < 10){
               data[newIdx] = Data[ts_max+k];
               dd = Data[ts_max+k];
               i = newIdx;
           }
           data[i] = dd;
           ++k;
       }
                                                                                                           
   }
}