Monitoring laser eo p

EOverPCalibration/interface/ConvoluteTemplate.h

@@ -0,0 +1,14 @@
+#ifndef ConvoluteTemplate_h
+#define ConvoluteTemplate_h
+
+#include "histoFunc.h"
+
+#include "TH1.h"
+#include "TF1.h"
+#include "TVirtualFFT.h"
+
+
+TH1F* ConvoluteTemplate(const std::string& name, TH1F* h_template, TH1F* h_smearing,
+                        int nPoints, double min, double max);
+
+#endif

EOverPCalibration/interface/histoFunc.h

@@ -0,0 +1,87 @@
+#ifndef histoFunc_h
+#define histoFunc_h
+
+#include "TH1.h"
+
+
+
+
+
+
+class histoFunc
+{
+ public:
+
+
+  //! ctor
+  histoFunc(TH1F* histo)
+  {
+    histo_p = histo;
+  };
+
+
+  //! dtor
+  ~histoFunc()
+  {};
+
+
+  //! operator()
+  double operator()(double* x, double* par)
+  {
+    double xx = par[1] * (x[0] - par[2]);
+
+    double xMin = histo_p -> GetBinCenter(1);
+    double xMax = histo_p -> GetBinCenter(histo_p -> GetNbinsX());
+
+
+
+    if( (xx < xMin) || (xx >= xMax) )
+      return 1.e-10;
+
+    else
+    {  
+      int bin = histo_p -> FindBin(xx);
+      int bin1 = 0;
+      int bin2 = 0;
+
+      if(xx >= histo_p -> GetBinCenter(bin))
+      {
+        bin1 = bin;
+        bin2 = bin+1;
+      }
+
+      else
+      {
+        bin1 = bin-1;
+        bin2 = bin;
+      }
+
+
+      double x1 = histo_p -> GetBinCenter(bin1);
+      double y1 = histo_p -> GetBinContent(bin1);
+
+      double x2 = histo_p -> GetBinCenter(bin2);
+      double y2 = histo_p -> GetBinContent(bin2);
+
+      double m = 1. * (y2 - y1) / (x2 - x1);
+
+
+
+      if( (y1 + m * (xx - x1)) < 1.e-10)
+        return 1.e-10;
+
+
+      return par[0] * par[1] * (y1 + m * (xx - x1));
+    }
+
+    return 1.e-10;  
+  }
+
+
+
+ private:
+
+  TH1F* histo_p;
+};
+
+#endif

EOverPCalibration/interface/ntpleUtils2.h

@@ -0,0 +1,79 @@
+#ifndef ntupleUtils2_h
+#define ntupleUtils2_h
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <cmath>
+
+
+#include "TFile.h"
+#include "TChain.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TF1.h"
+
+#ifdef _MAKECINT_
+#pragma link C++ class vector<Sorter>+;
+#pragma link C++ class vector<SorterLC>+;
+#pragma link C++ class map<int,TH1F*>+;
+#pragma link C++ class map<string,TH1F*>+;
+#endif
+
+
+extern TH1F* templateHisto;
+extern TF1* templateFunc;
+
+extern std::vector<double>* mydata;
+
+void FitTemplate(const bool& draw = false);
+
+/*** double crystall ball ***/
+double crystalBallLowHigh_v2(double* x, double* par);
+
+/*** time sort a tree ***/
+struct Sorter
+{
+  int time;
+  int entry;
+
+  bool operator() (const Sorter& s1, const Sorter& s2)
+  {
+    return s1.time < s2.time;
+  }
+};
+
+/*** time sort a tree ***/
+struct myEvent
+{
+  int runId;
+  int timeStampHigh;
+
+  int region;
+
+  float scE;
+  float P;
+
+  float scLaserCorr;
+  float seedLaserAlpha;
+
+  bool operator() (const myEvent& s1, const myEvent& s2)
+  {
+    return s1.timeStampHigh < s2.timeStampHigh;
+  }
+};
+
+/*** time sort a tree ***/
+struct SorterLC
+{
+  float laserCorr;
+  int entry;
+
+  bool operator() (const SorterLC& s1, const SorterLC& s2)
+  {
+    return s1.laserCorr < s2.laserCorr;
+  }
+};
+
+
+#endif

EOverPCalibration/interface/stabilityUtils.h

@@ -0,0 +1,24 @@
+#ifndef stabilityUtils_h
+#define stabilityUtils_h
+
+#include "histoFunc.h"
+#include <TString.h>
+#include <TObjArray.h>
+#include <TObjString.h>
+
+#include <sstream>
+#include <ctime>
+
+#include "TH1F.h"
+#include "TF1.h"
+
+#include <iostream>
+#include <string>
+
+int dateToInt(const std::string& date);
+
+void SetHistoStyle(TH1* h, const std::string& label = "");
+
+TH1F* MellinConvolution(const std::string& name, TH1F* h_template, TH1F* h_Las);
+
+#endif

EOverPCalibration/src/ConvoluteTemplate.cc

@@ -0,0 +1,157 @@
+#include "../interface/ConvoluteTemplate.h"
+
+
+
+TH1F* ConvoluteTemplate(const std::string& name, TH1F* h_template, TH1F* h_smearing,
+                        int nPoints, double min, double max)
+{ 
+  double width = 1.*(max-min)/nPoints;
+
+  double* FFT_re_func1 = new double[nPoints];
+  double* FFT_im_func1 = new double[nPoints];
+  double* FFT_re_func2 = new double[nPoints];
+  double* FFT_im_func2 = new double[nPoints];
+  double* FFT_re_convolution = new double[nPoints];
+  double* FFT_im_convolution = new double[nPoints];
+
+
+
+  //-----------------------------
+  // define the initial functions
+
+  char funcName[50];
+
+
+  sprintf(funcName,"f_template_temp");
+  histoFunc* hf_template = new histoFunc(h_template);
+  TF1* f_template = new TF1(funcName, hf_template, min, max, 3, "histoFunc");
+
+  f_template -> FixParameter(0,1.);
+  f_template -> FixParameter(1,1.);
+  f_template -> FixParameter(2,0.);
+
+
+  sprintf(funcName,"f_smearing_temp");  
+  histoFunc* hf_smearing = new histoFunc(h_smearing);
+  TF1* f_smearing = new TF1(funcName, hf_smearing, min, max, 3, "histoFunc");
+  f_smearing -> FixParameter(0,1.);
+  f_smearing -> FixParameter(1,1.);
+  f_smearing -> FixParameter(2,-1.*h_smearing->GetMean());
+
+
+  // use a gaussian as smearing function
+  //TF1* f_smearing = new TF1("f_smearing","1./([0]*sqrt(2.*3.14159))*exp(-1.*x*x/(2.*[0]*[0]))",min,max);
+  //f_smearing -> FixParameter(0,0.05);  
+
+
+
+  //-----------------------------------------
+  // define the histograms to contain the FFT
+  TH1D* h_func1     = new TH1D("h_func1",    "",nPoints,min,max);
+  TH1D* h_func1_FFT = new TH1D("h_func1_FFT","",nPoints,min,max);
+
+  TH1D* h_func2     = new TH1D("h_func2",    "",nPoints,min,max);
+  TH1D* h_func2_FFT = new TH1D("h_func2_FFT","",nPoints,min,max);
+
+  for(int bin = 1; bin <= nPoints; ++bin)
+  {
+    h_func1 -> SetBinContent(bin,f_template->Eval(min+width*(bin-1)));
+    h_func2 -> SetBinContent(bin,f_smearing->Eval(min+width*(bin-1)));
+  }
+
+  //h_func1 -> Write();
+  //h_func2 -> Write();
+
+
+
+  //-----------
+  // do the FFT
+
+  h_func1 -> FFT(h_func1_FFT,"MAG R2C M");
+  TVirtualFFT* FFT_func1 = TVirtualFFT::GetCurrentTransform();
+  FFT_func1 -> GetPointsComplex(FFT_re_func1,FFT_im_func1);
+  FFT_func1 -> Delete();
+
+  h_func2 -> FFT(h_func2_FFT,  "MAG R2C M");
+  TVirtualFFT* FFT_func2 = TVirtualFFT::GetCurrentTransform();
+  FFT_func2 -> GetPointsComplex(FFT_re_func2,FFT_im_func2);
+  FFT_func2 -> Delete();
+
+
+
+  //---------------------------------
+  // convolution in the Fourier space
+
+  for(int bin = 1; bin <= nPoints; ++bin)
+    {
+      FFT_re_convolution[bin-1] = FFT_re_func1[bin-1]*FFT_re_func2[bin-1] - FFT_im_func1[bin-1]*FFT_im_func2[bin-1];
+      FFT_im_convolution[bin-1] = FFT_re_func1[bin-1]*FFT_im_func2[bin-1] + FFT_im_func1[bin-1]*FFT_re_func2[bin-1];
+    }
+
+
+
+  //----------------
+  // do the anti FFT
+
+  TVirtualFFT* AFFT_convolution = TVirtualFFT::FFT(1, &nPoints, "C2R M K");
+  AFFT_convolution -> SetPointsComplex(FFT_re_convolution,FFT_im_convolution);
+  AFFT_convolution -> Transform();  
+
+  TH1* h_convolution_temp = NULL;
+  h_convolution_temp = TH1::TransformHisto(AFFT_convolution,h_convolution_temp,"Re");
+  h_convolution_temp -> SetName("h_convolution_temp");
+  h_convolution_temp -> Scale((max-min)/nPoints/nPoints);
+  AFFT_convolution -> Delete();
+
+  TH1D* h_convolution = new TH1D("h_convolution","",nPoints,min,max);
+  for(int bin = 1; bin <= nPoints; ++bin)
+  {
+    h_convolution -> Fill((max-min)/nPoints*(bin-1),h_convolution_temp->GetBinContent(bin));
+  }
+
+  histoFunc* convolutionHistoFunc = new histoFunc((TH1F*)(h_convolution));
+  TF1* f_convolution = new TF1("f_convolution", convolutionHistoFunc, min, max, 3, "histoFunc");
+  f_convolution -> FixParameter(0, 1);
+  f_convolution -> FixParameter(1, 1);
+  f_convolution -> FixParameter(2, 0.);  
+
+
+
+  int nBins = h_template->GetNbinsX();
+  float xMin = h_template->GetBinLowEdge(1);
+  float xMax = h_template->GetBinLowEdge(nBins) + h_template->GetBinWidth(nBins);
+  float xWidth = (xMax-xMin)/nBins;
+  TH1F* h_final = new TH1F(name.c_str(),"",nBins,xMin,xMax);
+
+  for(int bin = 1; bin <=nBins; ++bin)
+  {
+    float xMinBin = h_final->GetBinLowEdge(bin);
+    float xMaxBin = h_final->GetBinLowEdge(bin) + xWidth;
+
+    h_final -> SetBinContent(bin,f_convolution->Integral(xMinBin,xMaxBin)/xWidth);
+  }
+
+  h_final -> Scale(h_template->Integral()/h_final->Integral());
+  //h_final -> Write(); 
+
+
+
+  f_template -> Delete();
+  f_smearing -> Delete();
+
+  h_func1_FFT -> Delete();
+  h_func1 -> Delete(); 
+
+  h_func2_FFT -> Delete();
+  h_func2 -> Delete();
+
+  h_convolution_temp -> Delete();
+  h_convolution      -> Delete();
+
+  f_convolution -> Delete();
+  delete convolutionHistoFunc;
+
+
+
+  return h_final;
+}