11 #include <dqm_core/AlgorithmConfig.h>
14 #include <dqm_core/AlgorithmManager.h>
20 #include <TObjArray.h>
33 TH1::AddDirectory(
false);
50 const TObject&
object,
51 const dqm_core::AlgorithmConfig&
config ){
62 if(
object.
IsA()->InheritsFrom(
"TH1" ) ) {
65 throw dqm_core::BadConfig( ERS_HERE,
name,
"dimension != 2 " );
68 throw dqm_core::BadConfig( ERS_HERE,
name,
"does not inherit from TH1" );
84 catch( dqm_core::Exception & ex ) {
85 throw dqm_core::BadConfig( ERS_HERE,
name, ex.what(), ex );
102 if (
histogram->GetEntries() < minStat ) {
112 std::vector<int>
range;
116 catch( dqm_core::Exception & ex ) {
117 throw dqm_core::BadConfig( ERS_HERE,
name, ex.what(), ex );
128 if (maxPublish > 999) maxPublish=999;
139 bool findBinsOverAvg =
true;
140 bool findBinsUnderAvg =
true;
142 if (!lessThan && greaterThan) {
143 findBinsUnderAvg =
false;
145 else if (!greaterThan && lessThan) {
146 findBinsOverAvg =
false;
164 if (minBinsBeforeSkimming < 3) {
165 if ( !testConsistencyWithErrors ) minBinsBeforeSkimming = 50;
166 else minBinsBeforeSkimming = 8;
170 if (minBinsAfterSkimming < 2) {
171 if ( !testConsistencyWithErrors ) minBinsAfterSkimming = 8;
172 else minBinsAfterSkimming = 3;
199 config.getParameters(), 0.3);
219 std::vector<double> xBinCenters;
220 std::vector<double> yBinCenters;
221 double * xBinEdges =
new double[ixmax+1];
222 double * yBinEdges =
new double[iymax+1];
225 xBinCenters.push_back(
histogram->GetXaxis()->GetBinLowEdge(
range[0]) );
227 xBinCenters.push_back(
histogram->GetXaxis()->GetBinCenter(
i) );
231 xBinCenters.push_back( xBinEdges[ixmax] );
234 yBinCenters.push_back(
histogram->GetYaxis()->GetBinLowEdge(
range[2]) );
236 for (
int j =
range[2]; j <=
range[3]; j++ ) {
237 yBinCenters.push_back(
histogram->GetYaxis()->GetBinCenter(j) );
241 yBinCenters.push_back( yBinEdges[iymax] );
243 double * sBinEdges = xBinEdges;
245 sBinEdges = yBinEdges;
249 TH2F* inputBins =
new TH2F(
"inputBins",
histogram->GetTitle(), ixmax, xBinEdges, iymax, yBinEdges);
250 TH2F* binwiseStatus =
new TH2F(
"binewiseStatus",
"BinsDiffByStrips Status Results per Bin", ixmax, xBinEdges,
252 TH2F* binDeviations =
new TH2F(
"binDeviations",
"Input bin content normalized by skimmed average and variance per strip",
253 ixmax, xBinEdges, iymax, yBinEdges);
254 TH1F* stripAverages =
new TH1F(
"stripAverages",
"Average Value from Cells in Strip After Removal of Outliers",
256 TH1F* stripVariances =
new TH1F(
"stripVariances",
"Standard Variance of Cells in Strip After Removal of Outliers",
262 double maxDeviation = 0;
263 std::vector <tools::binContainer> yellowBins;
264 std::vector <tools::binContainer> redBins;
268 int nBinsUndefined = 0;
276 std::vector<std::vector<tools::binContainer> > strips;
278 std::vector<tools::binContainer> AllBinInOneStrip;
279 double emptyBinCounter=0;
283 for (
int is = 1; is <= ismax; is++ ) {
284 std::vector<tools::binContainer>
inputs;
292 for (
int in = 1; in <= inmax; in++ ) {
304 inputBins->SetBinContent(ix,iy,
value);
305 inputBins->SetBinError(ix,iy,
error);
306 if(
value==0) emptyBinCounter++;
308 AllBinInOneStrip.push_back(binContent_tmp);
309 if( (
value != ignoreValue) && (
error > minError) ){
312 inputs.push_back(binContent);
315 binwiseStatus->SetBinContent(ix,iy,dqm_core::Result::Disabled);
321 if( (!combineStripsBeforeSkimming) && ((
int)
inputs.size() < minBinsBeforeSkimming) ) {
322 nBinsUndefined +=
inputs.size();
324 else if(
inputs.size() != 0) {
330 if( combineStripsBeforeSkimming ) {
336 std::vector<double> inputValues;
337 std::vector<double> inputErrors;
338 std::vector<double> stripAveragesVector;
339 std::vector<double> stripErrors;
341 std::vector<tools::binContainer> deviations;
343 for( std::vector<std::vector<tools::binContainer> >::
iterator stripItr = strips.begin(); stripItr != strips.end(); ++stripItr ) {
345 if ( stripItr->empty() ) {
353 int nTot = stripItr->size();
358 tools::findOutliers(*stripItr, stripAvg, iScale, nIn, nIter, iVarExp , iterThresh, isbf, ibc);
366 if ( (nIn < minBinsAfterSkimming) ) {
367 for (std::vector<tools::binContainer>::const_iterator
it = stripItr->begin();
it != stripItr->end(); ++
it) {
369 deviations.push_back(deviationContainer);
375 double stripVariance = 0;
376 double stripVarianceError = 0;
377 double stripAvgError = 0;
380 if( !useMeanErrorForScale ) {
388 double sumSquaredDiffFromAvg = 0;
389 double sumCompensator = 0;
390 double err2Diff2Sum = 0;
396 for (std::vector<tools::binContainer>::const_iterator
it = stripItr->begin();
it != stripItr->end(); ++
it) {
409 double diffFromAvg =
it->value - stripAvg;
410 sumSquaredDiffFromAvg +=
std::pow( diffFromAvg, 2);
411 sumCompensator += ( diffFromAvg );
415 err2Sum += inputErr2;
416 err2Diff2Sum += inputErr2 *
std::pow( diffFromAvg, 2);
425 double countingVariance = 0;
426 double countingWeight = 0;
428 if( (nTot != 0) && (
range != 0) && ((nTot - nIn) != 0) ) {
429 double erfin = TMath::ErfInverse( (1.0 * nIn) / nTot );
433 countingVariance =
range / (erfin * 2 * std::sqrt(2));
439 double boundedVariance = 0;
440 double boundedWeight = 0;
441 double sumWeights = 0;
443 S2 = (sumSquaredDiffFromAvg - (
std::pow(sumCompensator,2)/nIn) )/(nIn - 1);
444 boundedVariance = std::sqrt(
S2);
445 double boundEffect = 1;
447 if( countingVariance != 0 ) {
448 U =
range / (2 * countingVariance);
450 else if( boundedVariance != 0 ) {
451 U =
range / (2 * boundedVariance);
454 boundEffect = TMath::Erf( U / 2 ) ;
456 if(boundEffect != 0) {
457 boundedVariance = boundedVariance / boundEffect;
459 double boundedErr2 = 0;
463 if( boundedErr2 != 0 ) {
464 boundedWeight = 1/boundedErr2;
466 sumWeights = ( countingWeight + boundedWeight );
468 if(sumWeights != 0) {
470 stripVariance = (countingVariance * countingWeight + boundedVariance * boundedWeight ) / sumWeights;
471 if(stripVariance != 0 ) {
472 boundEffect = TMath::Erf(
range / (4 * stripVariance) );
473 if( boundEffect != 0 ) {
474 boundedVariance = std::sqrt(
S2) / boundEffect;
482 else if ( iScale != 0 ){
484 boundedVariance = iScale;
485 boundedWeight =
std::pow(iScale, -2);
489 sumWeights = ( countingWeight + boundedWeight );
492 if(sumWeights != 0) {
493 stripVariance = (countingVariance * countingWeight + boundedVariance * boundedWeight ) / sumWeights;
494 stripVarianceError = 1./std::sqrt(sumWeights);
500 for (std::vector<tools::binContainer>::const_iterator
it = stripItr->begin();
it != stripItr->end(); ++
it) {
507 stripVariance = std::sqrt( err2Sum / nIn );
508 stripVarianceError = stripVariance / std::sqrt( nIn );
513 stripAvgError = std::sqrt(err2Sum)/nIn;
517 is = stripItr->front().iy;
520 is = stripItr->front().ix;
524 stripAverages->SetBinContent(is,stripAvg);
525 stripAverages->SetBinError(is,stripAvgError);
526 stripVariances->SetBinContent(is,stripVariance);
527 stripVariances->SetBinError(is,stripVarianceError);
530 if(doChiSquaredTest) {
531 for (std::vector<tools::binContainer>::const_iterator
it = stripItr->begin();
it != stripItr->end(); ++
it) {
532 inputValues.push_back(
it->value );
533 inputErrors.push_back( 0. );
534 stripAveragesVector.push_back( stripAvg );
535 stripErrors.push_back( stripVariance );
540 for (std::vector<tools::binContainer>::const_iterator
it = stripItr->begin();
it != stripItr->end(); ++
it) {
541 double deviation = 0;
542 double deviationError = -1;
543 double diffFromAvg =
it->value - stripAvg;
546 if(stripVariance > 0.00001){
547 deviation = diffFromAvg / stripVariance;
556 if( (
it->test==0) && (std::abs(
it->value)/std::abs(stripAvg)) > outstandingRatio ) {
561 if ( std::abs(diffFromAvg) < (absDiffGreenThresh + std::sqrt(
std::pow(
it->error,2) + (err2Sum/nIn))) ) {
566 deviations.push_back(deviationContainer);
568 int bin = binDeviations->GetBin(
it->ix,
it->iy);
569 binDeviations->SetBinContent(
bin,deviation);
570 binDeviations->SetBinError(
bin,deviationError);
576 std::vector<tools::binCluster>
clusters;
577 std::vector<tools::binCluster> redClusters;
578 std::vector<tools::binCluster> yellowClusters;
585 std::vector<std::vector<tools::binContainer*> > binMap =
makeBinMap(deviations, ixmax, iymax, topology);
588 if( (
it->value > seedThreshold +
it->error) && !
it->test ) {
604 double fabsDeviation = std::abs(
it->value);
605 double significanceBound = sigmaThresh *
it->error;
606 bool overAvg = (
it->value >= 0 );
607 if ( (findBinsOverAvg && overAvg) || (findBinsUnderAvg && !overAvg) ) {
608 if ( fabsDeviation > (gthreshold + significanceBound) ) {
609 if ( fabsDeviation > (rthreshold + significanceBound) ) {
611 if (publish || publishRed) {
612 redClusters.push_back(*
it);
616 nBinsYellow +=
it->n;
618 yellowClusters.push_back(*
it);
628 for (std::vector<tools::binContainer>::const_iterator
it = deviations.begin();
it != deviations.end(); ++
it) {
629 int bin = binwiseStatus->GetBin(
it->ix,
it->iy);
631 if (
it->error < 0 ) {
637 if (std::abs(
it->value) > std::abs(maxDeviation) ) {
639 maxDeviation =
it->value;
643 if (
it->test == 3 ) {
650 double fabsDeviation = std::abs(
it->value);
651 double significanceBound = sigmaThresh *
it->error;
654 if ( fabsDeviation < (gthreshold - significanceBound ) ) {
659 bool overAvg = (
it->value >= 0 );
660 if ( (findBinsOverAvg && overAvg) || (findBinsUnderAvg && !overAvg) ) {
661 if ( fabsDeviation > (gthreshold + significanceBound) ) {
662 if ( fabsDeviation > (rthreshold + significanceBound) ) {
667 if(
it->test != 10 ) {
669 if (publish || publishRed) {
670 redBins.push_back( *
it );
677 binwiseStatus->SetBinContent(
bin,dqm_core::Result::Yellow);
679 if(
it->test != 10 ) {
681 yellowBins.push_back( *
it );
694 double maxvalue_pre=-1;
695 std::vector<tools::binContainer> topBinEntries;
696 std::vector<tools::binContainer> topDeviations;
697 double emptyRatio_this = emptyBinCounter/(
histogram->GetNbinsX()*
histogram->GetNbinsY());
698 if (emptyRatio_this > emptyRatio) name_flag=1;
701 if(AllBinInOneStrip.size()>10 ) NTopEntries=10;
702 else NTopEntries = AllBinInOneStrip.size();
703 int AllBinInOneStrip_size = AllBinInOneStrip.size();
704 bool *bin_entries_status =
new bool [AllBinInOneStrip_size];
705 for(
int i=0;
i< AllBinInOneStrip_size;
i++) bin_entries_status[
i] =
true;
706 for(
int i=0;
i<NTopEntries;
i++){
710 for(std::vector<tools::binContainer>::const_iterator
it =AllBinInOneStrip.begin();
it!=AllBinInOneStrip.end();++
it){
711 int flag_my =
i==0 || ( bin_entries_status[
counter] &&
it->value <= maxvalue_pre);
712 if(std::abs(
it->value) >= std::abs(maxvalue) && flag_my) {
713 maxvalue =
it->value;
719 bin_entries_status[counter2]=0;
720 if (onebin_my.
x!= onebin_my_pre.
x||onebin_my.
y!=onebin_my_pre.
y) topBinEntries.push_back(onebin_my);
721 maxvalue_pre = maxvalue;
722 onebin_my_pre = onebin_my;
724 delete[] bin_entries_status;
728 if(deviations.size()>5) NTopdeviation = 5;
729 else NTopdeviation = deviations.size();
730 int deviations_size = deviations.size();
731 bool *bin_dev_status =
new bool [ deviations_size];
732 for(
int i=0;
i< deviations_size ;
i++) bin_dev_status[
i] =
true;
733 for(
int i=0;
i<NTopdeviation;
i++){
737 if(deviations.size()!=0){
738 for (std::vector<tools::binContainer>::const_iterator
it = deviations.begin();
it != deviations.end(); ++
it){
739 int flag_my =
i==0 || ( bin_dev_status[
counter] &&
it->value <= maxvalue_pre);
740 if(std::abs(
it->value) >= std::abs(maxvalue) && flag_my) {
741 maxvalue =
it->value;
747 bin_dev_status[counter2]=0;
748 if (onebin_my.
x!= onebin_my_pre.
x||onebin_my.
y!=onebin_my_pre.
y) {
749 topDeviations.push_back(onebin_my);
750 maxvalue_pre = maxvalue;
751 onebin_my_pre = onebin_my;
756 delete[] bin_dev_status;
758 if ( publish || publishRed) {
760 int objectsPublished = 0;
761 int clustersPublished = 0;
764 std::vector<tools::binCluster>::const_reverse_iterator rbcrbegin = redClusters.rbegin();
765 std::vector<tools::binCluster>::const_reverse_iterator rbcrend = redClusters.rend();
766 for (std::vector<tools::binCluster>::const_reverse_iterator
it = rbcrbegin;
767 it != rbcrend; ++
it) {
769 if (objectsPublished < maxPublish) {
771 sprintf(ctag,
"C%.3i-R-%.3iBins@ Eta=(%+.3f_to_%+.3f) Phi=(%+.3f_to_%+.3f) Center=(%+.3f,%+.3f) Radius=%+.3f",clustersPublished,
it->n,
772 xBinCenters[
it->ixmin],xBinCenters[
it->ixmax],yBinCenters[
it->iymin],yBinCenters[
it->iymax],
it->x,
it->y,
it->radius);
774 std::string
tag = ctag;
775 int sizeDiff = 30 -
tag.size();
777 tag.append(sizeDiff,
'_');
788 int binsPublished = 0;
791 for ( std::vector<tools::binContainer>::const_reverse_iterator rIter = redBins.rbegin();
792 rIter != redBins.rend();
794 if (objectsPublished < maxPublish) {
796 sprintf(ctag,
"%.3i-R-",binsPublished);
797 std::string
tag = ctag;
811 for ( std::vector<tools::binContainer>::const_reverse_iterator rIter = yellowBins.rbegin();
812 rIter != yellowBins.rend();
814 if (objectsPublished < maxPublish) {
816 sprintf(ctag,
"%.3i-Y-",binsPublished);
817 std::string
tag = ctag;
832 for (std::vector<tools::binCluster>::const_reverse_iterator
it = yellowClusters.rbegin();
833 it != yellowClusters.rend(); ++
it) {
834 if (objectsPublished < maxPublish) {
836 sprintf(ctag,
"C%.3i-Y-%.3iBins@ Eta=(%+.3f_to_%+.3f) Phi=(%+.3f_to_%+.3f) Center=(%+.3f,%+.3f) Radius=%+.3f",clustersPublished,
it->n,
837 xBinCenters[
it->ixmin],xBinCenters[
it->ixmax],yBinCenters[
it->iymin],yBinCenters[
it->iymax],
it->x,
it->y,
it->radius);
838 std::string
tag = ctag;
839 int sizeDiff = 30 -
tag.size();
841 tag.append(sizeDiff,
'_');
851 result->tags_[
"NBins_RED"] = nBinsRed;
852 result->tags_[
"NBins_YELLOW"] = nBinsYellow;
856 for(
unsigned int i=0;
i<topDeviations.size();
i++){
857 if(topDeviations[
i].
error !=-999.9 ) {
859 sprintf(
tmp,
"MaxDeviation%u-",
i);
860 std::string myString =
tmp;
862 result->tags_[myString] = topDeviations[
i].value;
867 if( maxDevBin.
error != -999.9 ) {
868 std::string devString =
"MaxDeviation-";
870 result->tags_[devString] = maxDeviation;
872 for(
unsigned int i=0;
i<topBinEntries.size();
i++){
874 sprintf(
tmp,
"LeadingBinContents%u-",
i);
875 std::string myString =
tmp;
877 result->tags_[myString] = topBinEntries[
i].value;
882 result->tags_[
"Algorithm--BinsDiffByStrips"] = 5;
887 if( andGroup != -99999 ) {
888 result->tags_[
"AndGroup"] = andGroup;
892 TObjArray * resultList =
new TObjArray(5,0);
896 resultList->Add(binwiseStatus);
897 resultList->Add(inputBins);
898 resultList->Add(binDeviations);
899 resultList->Add(stripAverages);
900 resultList->Add(stripVariances);
902 resultList->SetOwner(
true);
904 result->object_ = (boost::shared_ptr<TObject>)(TObject*)(resultList);
907 if( doChiSquaredTest ) {
910 result->tags_[
"SigmaChiSq"] = chiSquareResult.second;
912 if ( testConsistencyWithErrors ) {
914 if ( chiSquareResult.second <= gthreshold ) {
916 }
else if ( chiSquareResult.second < rthreshold ) {
917 result->status_ = dqm_core::Result::Yellow;
927 int nActiveBins = nBinsGreen + nBinsRed + nBinsYellow + nBinsUndefined;
928 if( nActiveBins != 0 ) {
929 if ( (nBinsRed >= nRedBinsToRedStatus) || ((nBinsRed * 1.0 / nActiveBins) > redFracToRedStatus) ) {
932 else if ( ((nBinsRed + nBinsYellow) >= nYellowBinsToYellowStatus)
933 || (((nBinsRed + nBinsYellow) / nActiveBins) > yellowFracToYellowStatus) ) {
934 result->status_ = dqm_core::Result::Yellow;
936 else if ( (nBinsGreen * 1.0 / nActiveBins ) >= greenFracToGreenStatus ) {
944 result->status_ = dqm_core::Result::Disabled;
954 out<<
"BinsDiffByStrips: Calculates Average bin value and variance for each strip of constant x and finds bins in that strip that our outliers from the mean\n"<<std::endl;
955 out<<
"Mandatory Green/Red Threshold: MaxDeviation: size of deviation from mean, in multiples of variance, for bin to be considered Green/Red. Overall result is result for worst bin"<<std::endl;
956 out<<
"Optional Parameter: MinStat: Minimum histogram statistics needed to return defined result"<<std::endl;
957 out<<
"Optional Parameter: ignoreval: value to be ignored for calculating average, variance, and identifying bad bins. Ignored bins are flagged as Disabled"<<std::endl;
958 out<<
"Optional Parameter: useStripsOfConstantY: compare bins in strips of constant Y rather than constant X, as is the defualt (set to 1 for true). Default = false"<<std::endl;
960 out<<
"Optional Parameter: xmin: minimum x range"<<std::endl;
961 out<<
"Optional Parameter: xmax: maximum x range"<<std::endl;
962 out<<
"Optional Parameter: ymin: minimum y range"<<std::endl;
963 out<<
"Optional Parameter: ymax: maximum y range"<<std::endl;
964 out<<
"Optional Parameter: GreaterThan: check for bins which are GreaterThan average (set to 1 for true). Default = true"<<std::endl;
965 out<<
"Optional Parameter: LessThan: check for bins which are LessThan average (set to 1 for true). Default = true"<<std::endl;
966 out<<
"Optional Parameter: PublishBins: Save all bin quality decisions in output histogram and save Red and Yellow bins in tag output (set to 1 for true). Default = false"<<std::endl;
967 out<<
"Optional Parameter: PublishRedBins: Save bins identified as Red in tag output(set to 1 for true). Default = false"<<std::endl;
968 out<<
"Optional Parameter: MaxPublish: Max number of bins to save in tag output. Default = 20"<<std::endl;
969 out<<
"Optional Parameter: MinBinsAfterSkimming : Minimum number of bins remaining after skimming for any bins to be flagged as Red/Yellow/Green. If threshold is not met, another attempt will be made at the iteration process with these few bin excluded in the first pass. Failing success, the entire strip will be flagged as Undefined. Default = 5"<<std::endl;
970 out<<
"Optional Parameter: MinAbsDiffFromAvg : Bins whose absolute difference from the mean is less than this amount, in a statistically significant fashion, will always be flagged as green."
971 <<
"Bin will be flagged as either Green or Undefined, depending on significance by which it passes this cut."<<std::endl;
972 out<<
"Optional Parameter: nIterations : Number of iterations to be used in calculating mean and variance while skimming outliers. Default = 10" <<std::endl;
973 out<<
"Optional Parameter: MaxNRetryIteration : Maxiumum number of times the iteration should be restarted should it fail (by leaving fewer than MinBinsAfterSkimming), reseeding with those cells that pass the iteration removed in the first pass. If MaxNRetryIteration = 0, will never retry. Default = 5. " <<std::endl;
974 out<<
"Optional Parameter: IterVariationExponent: exponent k used in calculating the generalized variance, iVar, where iVar = (sum[ (abs[ x - xbar ])^k ])^(1/k). If k = 2, for instance, iVar is equivalent to the standard variance. (If nIterations = 1, this parameter has no effect.) Default = -1.4235"<<std::endl;
975 out<<
"Optional Parameter: IterDeviationThresh : size of deviation from mean, in multiples of the generalized variance, for bin to be excluded during iterative calculation of mean and variance. (If nIterations = 1, this parameter has no effect.) Default = 1.3"<<std::endl;
977 out<<
"Optional Parameter: SigmaThresh : Minimum significance (distance from threshold in multiples of the error) for quality statement to be made about a bin). Default = 5"<<std::endl;
978 out<<
"Optional Parameter: FindOutliersUsingErrors : Switch to use outlier finding based on the pull of the mean of each bin given its error as compared to the error on the mean, rather than the default error independent aproach."<<std::endl;
979 out<<
"Optional Parameter: UseMeanErrorForScale : Use the average error of the non-outliers, rather than their estimated variance, as the scale against which deviations are measured."<<std::endl;
980 out<<
"Optional Parameter: DoChiSquaredTest : Calculate the Chi_Squared value for the distribution of all the bin values given the scale (estimated variance or, if UseMeanErrorForScale is used, mean error."<<std::endl;
981 out<<
"Optional Parameter: TestConsistencyWithErrors : Switch on FindOutliersUsingErrors, UseMeanErrorForScale, and DoChiSquaredTest and base the DQ decision on the Chi-Squared result"<<std::endl;
985 int where = name_tmp.find_last_of(
"_");
988 if(
name.compare(
"5Sigma")==0 || name_tmp.compare(
"m_EtavsPhi3")==0 || name_tmp.compare(
"m_EtavsPhi4")==0 ) name_flag=1;