diff --git a/EOverPCalibration/interface/ConvoluteTemplate.h b/EOverPCalibration/interface/ConvoluteTemplate.h new file mode 100644 index 00000000000..745b2d8945a --- /dev/null +++ b/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 diff --git a/EOverPCalibration/interface/histoFunc.h b/EOverPCalibration/interface/histoFunc.h new file mode 100644 index 00000000000..58347c6fc3b --- /dev/null +++ b/EOverPCalibration/interface/histoFunc.h @@ -0,0 +1,84 @@ +#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 diff --git a/EOverPCalibration/interface/ntpleUtils2.h b/EOverPCalibration/interface/ntpleUtils2.h new file mode 100644 index 00000000000..5c5da549b8d --- /dev/null +++ b/EOverPCalibration/interface/ntpleUtils2.h @@ -0,0 +1,76 @@ +#ifndef ntupleUtils2_h +#define ntupleUtils2_h + +#include +#include +#include +#include + + +#include "TFile.h" +#include "TChain.h" +#include "TCanvas.h" +#include "TH1F.h" +#include "TF1.h" + +#ifdef _MAKECINT_ +#pragma link C++ class vector+; +#pragma link C++ class vector+; +#pragma link C++ class map+; +#pragma link C++ class map+; +#endif + + +extern TH1F* templateHisto; +extern TF1* templateFunc; + +extern std::vector* 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 diff --git a/EOverPCalibration/interface/stabilityUtils.h b/EOverPCalibration/interface/stabilityUtils.h new file mode 100644 index 00000000000..717a323d688 --- /dev/null +++ b/EOverPCalibration/interface/stabilityUtils.h @@ -0,0 +1,24 @@ +#ifndef stabilityUtils_h +#define stabilityUtils_h + +#include "histoFunc.h" +#include +#include +#include + +#include +#include + +#include "TH1F.h" +#include "TF1.h" + +#include +#include + +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 diff --git a/EOverPCalibration/src/ConvoluteTemplate.cc b/EOverPCalibration/src/ConvoluteTemplate.cc new file mode 100644 index 00000000000..edc70fa39dd --- /dev/null +++ b/EOverPCalibration/src/ConvoluteTemplate.cc @@ -0,0 +1,153 @@ +#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; +} diff --git a/EOverPCalibration/src/ntpleUtils2.cc b/EOverPCalibration/src/ntpleUtils2.cc new file mode 100644 index 00000000000..dcc76e2529c --- /dev/null +++ b/EOverPCalibration/src/ntpleUtils2.cc @@ -0,0 +1,77 @@ +#include "../interface/ntpleUtils2.h" + +TH1F* templateHisto; +TF1* templateFunc; + +std::vector* mydata; + +/*** fit the template ***/ +void FitTemplate(const bool& draw) +{ + TCanvas* c_template; + if( draw ) { + c_template = new TCanvas("c_template", "template"); + c_template -> cd(); + c_template -> SetGridx(); + c_template -> SetGridy(); + + templateHisto -> Scale(1. / templateHisto->GetEntries()); + templateHisto -> SetFillColor(kCyan + 2); + templateHisto -> Draw(); + } + + templateFunc = new TF1("templateFunc", crystalBallLowHigh_v2, 0., 10., 8); + templateFunc -> SetNpx(10000); + templateFunc -> SetLineWidth(2); + templateFunc -> SetLineColor(kRed); + + templateFunc -> SetParameters(templateHisto->GetMaximum(), 1., 0.05, 1., 2., 2., 1., 1.); + templateFunc -> FixParameter(7, 1.); + templateFunc -> SetParLimits(3, 0., 10.); + templateFunc -> SetParLimits(5, 0., 10.); + + templateFunc -> SetParName(0, "N"); + templateFunc -> SetParName(1, "#mu"); + templateFunc -> SetParName(2, "#sigma"); + templateFunc -> SetParName(3, "#alpha_{high}"); + templateFunc -> SetParName(4, "n_{high}"); + templateFunc -> SetParName(5, "#alpha_{low}"); + templateFunc -> SetParName(6, "n_{low}"); + + templateHisto -> Fit("templateFunc", "NQR+", "", 0.5, 3.); + + if( draw ) + templateFunc -> Draw("same"); +} + + +/*** double crystall ball ***/ +double crystalBallLowHigh_v2(double* x, double* par) +{ + //[0] = N + //[1] = mean + //[2] = sigma + //[3] = alpha + //[4] = n + //[5] = alpha2 + //[6] = n2 + //[7] = scale + double xx = x[0] * par[7]; + double mean = par[1]; + double sigma = par[2]; + double alpha = par[3]; + double n = par[4]; + double alpha2 = par[5]; + double n2 = par[6]; + if( (xx - mean) / sigma > fabs(alpha) ) { + double A = pow(n / fabs(alpha), n) * exp(-0.5 * alpha * alpha); + double B = n / fabs(alpha) - fabs(alpha); + return par[0] * par[7] * A * pow(B + (xx - mean) / sigma, -1.*n); + } else if( (xx - mean) / sigma < -1.*fabs(alpha2) ) { + double A = pow(n2 / fabs(alpha2), n2) * exp(-0.5 * alpha2 * alpha2); + double B = n2 / fabs(alpha2) - fabs(alpha2); + return par[0] * par[7] * A * pow(B - (xx - mean) / sigma, -1.*n2); + } else { + return par[0] * par[7] * exp(-1. * (xx - mean) * (xx - mean) / (2 * sigma * sigma) ); + } +} diff --git a/EOverPCalibration/src/stabilityUtils.cc b/EOverPCalibration/src/stabilityUtils.cc new file mode 100644 index 00000000000..d1ce9bcffca --- /dev/null +++ b/EOverPCalibration/src/stabilityUtils.cc @@ -0,0 +1,141 @@ +#include "../interface/stabilityUtils.h" + +using namespace std; + +int dateToInt(const std::string& date) +{ + int day; + int month; + int year; + int tempo; + + TString date2 = date; + + TObjArray *splitted = date2.Tokenize("-"); + + TObjString *Objstring1 = (TObjString *) splitted->At(0); + TObjString *Objstring2 = (TObjString *) splitted->At(1); + TObjString *Objstring3 = (TObjString *) splitted->At(2); + + + TString string1 = Objstring1->GetString(); + TString string2 = Objstring2->GetString(); + TString string3 = Objstring3->GetString(); + + + std::cout << "day: " << string1 << std::endl; + std::cout << "month: " << string2 << std::endl; + std::cout << "year: " << string3 << std::endl; + + + std::string stringday = string1.Data(); + std::string stringmonth = string2.Data(); + std::string stringyear = string3.Data(); + + day = atoi(stringday.c_str()); + month = atoi(stringmonth.c_str()); + year = atoi(stringyear.c_str()); + /* + istringstream buffer1(stringday); + buffer1 >> day; + istringstream buffer2(stringmonth); + buffer2 >> month; + istringstream buffer3(stringyear); + buffer3 >> year; + */ + + + std::cout << "day: " << day << std::endl; + std::cout << "month: " << month << std::endl; + std::cout << "year: " << year << std::endl; + delete splitted; + + + tm time; + time.tm_sec = 0; + time.tm_min = 0; + time.tm_hour = 0; + time.tm_mday = day; + time.tm_mon = month - 1; + time.tm_year = year - 1900; + + + tempo = timegm(&time); + + std::cout << "tempo: " << tempo << std::endl; + + + return mktime(&time); + +} + + + + +void SetHistoStyle(TH1* h, const std::string& label) +{ + h -> SetLineWidth(2); + + h -> GetXaxis() -> SetTitleSize(0.04); + h -> GetXaxis() -> SetLabelSize(0.04); + + h -> GetYaxis() -> SetTitleSize(0.04); + h -> GetYaxis() -> SetLabelSize(0.04); + + h -> GetYaxis() -> SetTitleOffset(1.5); + + if( label == "EoP") { + h -> SetMarkerColor(kRed + 2); + h -> SetLineColor(kRed + 2); + h -> SetFillColor(kRed + 2); + h -> SetFillStyle(3004); + h -> SetMarkerStyle(7); + } + + if( label == "EoC") { + h -> SetMarkerColor(kGreen + 2); + h -> SetLineColor(kGreen + 2); + h -> SetFillColor(kGreen + 2); + h -> SetFillStyle(3004); + h -> SetMarkerStyle(7); + } +} + + + + + + +TH1F* MellinConvolution(const std::string& name, TH1F* h_template, TH1F* h_Las) +{ + histoFunc* templateHistoFunc = new histoFunc(h_template); + + TF1* f_template = new TF1("f_template", templateHistoFunc, 0.8 * (h_Las->GetMean()), 1.4 * (h_Las->GetMean()), 3, "histoFunc"); + + f_template -> SetNpx(10000); + f_template -> FixParameter(0, 1.); + f_template -> SetParameter(1, 1.); + f_template -> FixParameter(2, 0.); + + TH1F* h_convolutedTemplate = (TH1F*)( h_template->Clone(name.c_str()) ); + h_convolutedTemplate -> Reset(); + + for(int bin = 1; bin <= h_Las->GetNbinsX(); ++bin) { + float scale = h_Las->GetBinCenter(bin); + float weight = h_Las->GetBinContent(bin); + + if( weight > 0. ) { + f_template -> SetParameter(0, weight); + f_template -> SetParameter(1, 1. / scale); + + for(int bin2 = 1; bin2 <= h_convolutedTemplate->GetNbinsX(); ++bin2) { + float binCenter = h_convolutedTemplate -> GetBinCenter(bin2); + h_convolutedTemplate -> Fill(binCenter, f_template->Eval(binCenter)); + } + } + } + + h_convolutedTemplate -> Scale(h_template->Integral() / h_convolutedTemplate->Integral()); + + return h_convolutedTemplate; +} diff --git a/ZFitter/submit_monitoring_jobs.py b/ZFitter/submit_monitoring_jobs.py new file mode 100644 index 00000000000..e04b4abb2be --- /dev/null +++ b/ZFitter/submit_monitoring_jobs.py @@ -0,0 +1,108 @@ +#! /usr/bin/env python +import os +import glob +import math +from array import array +import sys +import time +import subprocess + +from optparse import OptionParser + +parser = OptionParser() + +parser.add_option('--generateOnly', action='store_true', dest='generateOnly', default=False, help='generate jobs only, without submitting them') + +(options, args) = parser.parse_args() + +localReco = "multifit" #type of local reco: "multifit" or "weights" + +currentDir = os.getcwd(); +CMSSWDir = currentDir+"/../"; + +fn = "Job_monitoring_"+localReco+"/Job_"+"EB"; +outScript = open(fn+".sh","w"); +command = "ZFitter.exe -f EoverPcalibration_batch_"+localReco+"_hadd.dat --evtsPerPoint 5000 --laserMonitoringEP --EBEE EB --yMIN 0.85 --LUMI 2.5" +print command; +outScript.write('#!/bin/bash'); +outScript.write("\n"+'cd '+CMSSWDir); +outScript.write("\n"+'eval `scram runtime -sh`'); +outScript.write("\n"+'cd -'); +outScript.write("\necho $PWD"); +outScript.write("\nls"); + +outScript.write("\necho \"copia1\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/"+localReco+".root ./") +outScript.write("\necho \"copia2\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/"+localReco+"-extraCalibTree.root ./") + +outScript.write("\necho \"copia13\" "); +# outScript.write("\ncp -v /afs/cern.ch/user/l/lbrianza/work/public/ntupleEoP/* .") +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/momentumCalibration2015_EB_pTk.root ./") +outScript.write("\necho \"copia14\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/momentumCalibration2015_EB_scE.root ./") +outScript.write("\necho \"copia15\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/EoverPcalibration_batch_"+localReco+"_hadd.dat ./") +outScript.write("\necho \"fine copia\" "); + +outScript.write("\nls") +outScript.write("\necho \"eseguo: "+command+"\" ") +outScript.write("\n"+command); +outScript.write("\nls") +outScript.write("\ncp -v -r EB__/ "+currentDir+"/EB_"+localReco+"/") +outScript.close(); +os.system("chmod 777 "+currentDir+"/"+fn+".sh"); +command2 = "bsub -q cmscaf1nd -cwd "+currentDir+" "+currentDir+"/"+fn+".sh"; +if not options.generateOnly: + os.system(command2); +print command2 + + + + + + + +fn = "Job_monitoring_"+localReco+"/Job_"+"EE"; +outScript = open(fn+".sh","w"); +command = "ZFitter.exe -f EoverPcalibration_batch_"+localReco+"_hadd.dat --evtsPerPoint 5000 --laserMonitoringEP --EBEE EE --yMIN 0.65 --yMAX 1.15 --LUMI 2.5" +print command; +outScript.write('#!/bin/bash'); +outScript.write("\n"+'cd '+CMSSWDir); +outScript.write("\n"+'eval `scram runtime -sh`'); +outScript.write("\n"+'cd -'); +outScript.write("\necho $PWD"); +outScript.write("\nls"); + + +outScript.write("\necho \"copia1\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/"+localReco+".root ./") +outScript.write("\necho \"copia2\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/"+localReco+"-extraCalibTree.root ./") + +outScript.write("\necho \"copia13\" "); +# outScript.write("\ncp -v /afs/cern.ch/user/l/lbrianza/work/public/ntupleEoP/* .") +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/momentumCalibration2015_EE_pTk.root ./") +outScript.write("\necho \"copia14\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/momentumCalibration2015_EE_scE.root ./") +outScript.write("\necho \"copia15\" "); +outScript.write("\ncmsStage /store/group/dpg_ecal/alca_ecalcalib/ecalMIBI/lbrianza/ntupleEoP/EoverPcalibration_batch_"+localReco+"_hadd.dat ./") +outScript.write("\necho \"fine copia\" "); + +outScript.write("\nls") +outScript.write("\necho \"eseguo: "+command+"\" ") +outScript.write("\n"+command); +outScript.write("\nls") +outScript.write("\ncp -v -r EE__/ "+currentDir+"/EE_"+localReco+"/") +outScript.close(); +os.system("chmod 777 "+currentDir+"/"+fn+".sh"); +command2 = "bsub -q cmscaf1nd -cwd "+currentDir+" "+currentDir+"/"+fn+".sh"; +if not options.generateOnly: + os.system(command2); +print command2 + + + + + +