EcalTiming: interface/EcalCrystalTimingCalibration.h Source File

00001 #include "EcalTiming/EcalTiming/interface/EcalTimingEvent.h"
00002 #include <vector>
00003 #include <TTree.h>
00004 
00005 /** \class EcalCrystalTimingCalibration EcalCrystalTimingCalibration.h EcalTiming/EcalTiming/interface/EcalCrystalTimingCalibration.h
00006  *
00007  * Description: add a description here
00008  * This class contains all the timing information for a single crystal
00009  */
00010 
00011 
00012 class EcalCrystalTimingCalibration
00013 {
00014 public:
00015 //  DetId _detId; ///< detId of the channel
00016 
00017 
00018 private:
00019 
00020     float _sum; ///< scalar sum of the time of each timingEvent
00021     float _sum2; ///< scalar sum of the square of the time of each timingEvent
00022     unsigned long int _num; ///< number of timingEvents;
00023 
00024     float _sumE; ///< scalar sum of the energy of each timingEvent: needed for average energy
00025 
00026     std::map<float, float> _sumWithinNSigma, _sum2WithinNSigma, _sum3WithinNSigma, _sumEWithinNSigma; ///< variables for calculation of mean, stdDev within n-times the origina stdDev (to remove tails)
00027     std::map<float, unsigned int> _numWithinNSigma; ///< variables for calculation of mean, stdDev within n-times the origina stdDev (to remove tails)
00028 
00029     std::vector<EcalTimingEvent> timingEvents; ///< vector containing  all the events for this crystal
00030     std::vector<EcalTimingEvent>::iterator maxChi2Itr;
00031 
00032 
00033 public:
00034     /// default constructor
00035     EcalCrystalTimingCalibration(bool weightMean = true) :
00036         //_detId(),
00037         _sum(0), _sum2(0), _num(0), _sumE(0)
00038         //totalChi2(-1),
00039         //useWeightedMean(weightMean)
00040     {
00041     }
00042 
00043     inline unsigned int num() const
00044     {
00045         return _num;
00046     };
00047     inline float mean() const
00048     {
00049         return _sum / _num;
00050     }; ///< average time (mean of the time distribution)
00051     inline float stdDev() const  ///< standard deviation of the time distribution
00052     {
00053         float mean_ = mean();
00054         return sqrt(_sum2 / _num - mean_ * mean_);
00055     };
00056     inline float meanError() const
00057     {
00058         return stdDev() / sqrt(_num);
00059     };
00060     inline float meanE() const
00061     {
00062         return _sumE / _num;
00063     }; ///< average Energy (mean of the Energy  distribution)
00064     /* bool operator<( EcalCrystalTimingCalibration& rhs) */
00065     /* { */
00066     /*  if(_detId < rhs._detId) return true; */
00067     /*  else return false; */
00068     /* } */
00069     //float totalChi2;
00070 
00071     float getMeanWithinNSigma(float sigma, float maxRange); ///< returns the mean time within abs(mean+ n * stdDev) to reject tails
00072     float getStdDevWithinNSigma(float sigma, float maxRange); ///< returns the stdDev calculated within abs(mean+ n * stdDev) to reject tails
00073     float getSkewnessWithinNSigma(float sigma, float maxRange); ///< returns the skewness calculated within abs(mean+ n * stdDev) to reject tails
00074 
00075     friend std::ostream& operator<< (std::ostream& os, const EcalCrystalTimingCalibration& s)
00076     {
00077         os << s.mean() << "\t" << s.stdDev() << "\t" << s.num();
00078         return os;
00079     }
00080 
00081     /// add new event for this crystal
00082     inline bool add(EcalTimingEvent te_)
00083     {
00084         return insertEvent(te_);
00085     }
00086     inline void clear()
00087     {
00088         _sum = 0.0f;
00089         _sum2 = 0.0f;
00090         _num = 0;
00091         _sumE = 0.0f;
00092         _sumWithinNSigma.clear();
00093         _sum2WithinNSigma.clear();
00094         _sum3WithinNSigma.clear();
00095         _numWithinNSigma.clear();
00096 
00097         timingEvents.clear();
00098     }
00099 
00100     void dumpToTree(TTree *tree, int ix_, int iy_, int iz_); ///< dump the full set of events in a TTree:  need an empty tree
00101 
00102     /// checks if the time measurement is stable changing the min energy threshold
00103     bool isStableInEnergy(float min, float max, float step);
00104 
00105 private:
00106     void calcAllWithinNSigma(float n_sigma, float maxRange = 10); ///< calculate sum, sum2, sum3, n for time if time within n x stdDev and store the result
00107     // since the values are stored, the calculation is done only once with only one loop over the events
00108 
00109     /// \todo weighted average by timeError
00110     bool insertEvent(EcalTimingEvent te_)
00111     {
00112         if(te_.timeError() > 0 && te_.timeError() < 1000 && te_.timeError() < 3) { //exclude values with wrong timeError estimation
00113             _sum += te_.time();
00114             _sum2 += te_.time() * te_.time();
00115             _sumE += te_.energy();
00116             _num++;
00117             timingEvents.push_back(te_);
00118             //updateChi2();
00119             return true;
00120         } else {
00121             return false;
00122         }
00123     }
00124 
00125     /* int filterOutliers(float threshold = 0.5) */
00126     /* { */
00127     /*  int numPointsErased = 0; */
00128 
00129     /*  while(timingEvents.size() > 4) { */
00130     /*      updateChi2(); */
00131     /*      float oldMean = mean; */
00132     /*      // Erase largest chi2 event */
00133     /*      EcalTimingEvent toRemove = *maxChi2Itr; */
00134     /*      timingEvents.erase(maxChi2Itr); */
00135     /*      //Calculate new mean/error */
00136     /*      updateChi2(); */
00137 
00138     /*      //Compare to old mean and break if |(newMean-oldMean)| < newSigma */
00139     /*      //TODO: study acceptance threshold */
00140     /*      if(fabs(mean - oldMean) < threshold * meanE) { */
00141     /*          insertEvent(toRemove); */
00142     /*          break; */
00143     /*      } else { */
00144     /*          numPointsErased++; */
00145     /*      } */
00146     /*  } */
00147 
00148     /*  return numPointsErased; */
00149     /* } */
00150 
00151 private:
00152     /* bool useWeightedMean; */
00153 
00154     /* //calculate chi2: assume a gaussian distribution */
00155     /* void updateChi2() // update individual, total, maxChi2s */
00156     /* { */
00157     /*  if(useWeightedMean) { */
00158     /*      updateChi2Weighted(); */
00159     /*  } else { */
00160     /*      updateChi2Unweighted(); */
00161     /*  } */
00162     /* } */
00163 
00164     /* void updateChi2Weighted() */
00165     /* { */
00166     /*  updateMeanWeighted(); */
00167     /*  float chi2 = 0; */
00168     /*  maxChi2Itr = timingEvents.begin(); */
00169 
00170     /*  for(std::vector<EcalTimingEvent>::iterator itr = timingEvents.begin(); */
00171     /*          itr != timingEvents.end(); ++itr) { */
00172     /*      float singleChi = (itr->time - mean) / itr->sigmaTime; */
00173     /*      itr->chi2 = singleChi * singleChi; */
00174     /*      chi2 += singleChi * singleChi; */
00175 
00176     /*      if(itr->chi2 > maxChi2Itr->chi2) { */
00177     /*          maxChi2Itr = itr; */
00178     /*      } */
00179     /*  } */
00180 
00181     /*  totalChi2 = chi2; */
00182     /* } */
00183 
00184     /* void updateChi2Unweighted() */
00185     /* { */
00186     /*  updateMeanUnweighted(); */
00187     /*  float chi2 = 0; */
00188     /*  maxChi2Itr = timingEvents.begin(); */
00189 
00190     /*  for(std::vector<EcalTimingEvent>::iterator itr = timingEvents.begin(); */
00191     /*          itr != timingEvents.end(); ++itr) { */
00192     /*      float singleChi = (itr->time - mean); */
00193     /*      itr->chi2 = singleChi * singleChi; */
00194     /*      chi2 += singleChi * singleChi; */
00195 
00196     /*      if(itr->chi2 > maxChi2Itr->chi2) { */
00197     /*          maxChi2Itr = itr; */
00198     /*      } */
00199     /*  } */
00200 
00201     /*  totalChi2 = chi2; */
00202     /* } */
00203 
00204     /* void updateMeanWeighted() */
00205     /* { */
00206     /*  float meanTmp = 0; */
00207     /*  float mean2Tmp = 0; */
00208     /*  float sigmaTmp = 0; */
00209 
00210     /*  for(std::vector<EcalTimingEvent>::const_iterator itr = timingEvents.begin(); */
00211     /*          itr != timingEvents.end(); ++itr) { */
00212     /*      float sigmaT2 = itr->sigmaTime; */
00213     /*      sigmaT2 *= sigmaT2; */
00214     /*      sigmaTmp += 1 / (sigmaT2); */
00215     /*      meanTmp += (itr->time) / (sigmaT2); */
00216     /*      mean2Tmp += ((itr->time) * (itr->time)) / (sigmaT2); */
00217     /*  } */
00218 
00219     /*  meanE = sqrt(1 / sigmaTmp); */
00220     /*  mean = meanTmp / sigmaTmp; */
00221     /*  rms = sqrt(mean2Tmp / sigmaTmp); */
00222     /*  stdDev = sqrt(rms * rms - mean * mean); */
00223     /* } */
00224 
00225     /* void updateMeanUnweighted() */
00226     /* { */
00227     /*  float meanTmp = 0; */
00228     /*  float mean2Tmp = 0; */
00229 
00230     /*  for(std::vector<EcalTimingEvent>::const_iterator itr = timingEvents.begin(); */
00231     /*          itr != timingEvents.end(); ++itr) { */
00232     /*      meanTmp += itr->time; */
00233     /*      mean2Tmp += (itr->time) * (itr->time); */
00234     /*  } */
00235 
00236     /*  mean = meanTmp / timingEvents.size(); */
00237     /*  rms = sqrt(mean2Tmp / timingEvents.size()); */
00238     /*  stdDev = sqrt(rms * rms - mean * mean); */
00239     /*  meanE = stdDev / sqrt(timingEvents.size()); // stdDev/sqrt(n) */
00240     /* } */
00241 
00242 };