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..956d87a6568 --- /dev/null +++ b/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 diff --git a/EOverPCalibration/interface/ntpleUtils2.h b/EOverPCalibration/interface/ntpleUtils2.h new file mode 100644 index 00000000000..d88d3838133 --- /dev/null +++ b/EOverPCalibration/interface/ntpleUtils2.h @@ -0,0 +1,79 @@ +#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..dce8053f1fb --- /dev/null +++ b/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; +} diff --git a/EOverPCalibration/src/ntpleUtils2.cc b/EOverPCalibration/src/ntpleUtils2.cc new file mode 100644 index 00000000000..486097c5d85 --- /dev/null +++ b/EOverPCalibration/src/ntpleUtils2.cc @@ -0,0 +1,83 @@ +#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..308a24e33ff --- /dev/null +++ b/EOverPCalibration/src/stabilityUtils.cc @@ -0,0 +1,146 @@ +#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: " < 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/bin/ZFitter.cpp b/ZFitter/bin/ZFitter.cpp index 47d6650d1b2..7e9efae9c7b 100644 --- a/ZFitter/bin/ZFitter.cpp +++ b/ZFitter/bin/ZFitter.cpp @@ -5,26 +5,49 @@ The aim of the program is to provide a common interface to all the Z fitting algorithms reading and combining in the proper way the configuration files. - \todo +\todo - remove commonCut from category name and add it in the ZFit_class in order to not repeate the cut - make alpha fitting more generic (look for alphaName) - Implement the iterative Et dependent scale corrections */ - -#include + #include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include -#include #include #include +#include #include #include #include +#include +#include #include #include +#include "TROOT.h" +#include "TSystem.h" +#include "TStyle.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TGraphErrors.h" +#include "TGraphAsymmErrors.h" +#include "TPaveStats.h" +#include "TLegend.h" +#include "TChain.h" +#include "TVirtualFitter.h" +#include "TLatex.h" + /// @cond SHOW /// \code @@ -45,6 +68,13 @@ configuration files. #include "../../EOverPCalibration/interface/CalibrationUtils.h" #include "../../EOverPCalibration/interface/FastCalibratorEB.h" #include "../../EOverPCalibration/interface/FastCalibratorEE.h" +//#include "setTDRStyle.h" +#include "../../EOverPCalibration/interface/ntpleUtils2.h" +#include "../../EOverPCalibration/interface/histoFunc.h" +#include "../../EOverPCalibration/interface/ConvoluteTemplate.h" +#include "../../EOverPCalibration/interface/stabilityUtils.h" +#include "../../EOverPCalibration/interface/geometryUtils.h" + /// \endcode /// @endcond @@ -168,7 +198,7 @@ void Dump(tag_chain_map_t& tagChainMap, TString tag="s", Long64_t firstentry=0){ /** * \param tagChainMap map of all the tags declared in the validation config file * \param tag name of the new \b tag created by the function, all the existent tags with name starting with -b tag are merged in the new \b tag + b tag are merged in the new \b tag * * A new tagChain with name=tag is added to the tagChainMap. All the tagChains with tag starting with \b tag are merged * After the merging the friend list is updated by \ref UpdateFriends @@ -182,10 +212,10 @@ void MergeSamples(tag_chain_map_t& tagChainMap, TString regionsFileNameTag, TStr for(tag_chain_map_t::const_iterator tag_chain_itr=tagChainMap.begin(); tag_chain_itr!=tagChainMap.end(); tag_chain_itr++){ - + // consider tags matching the tag input parameter if(tag_chain_itr->first.CompareTo(tag)==0 || !tag_chain_itr->first.Contains(tag)) continue; //do it for each sample - + // loop over all the trees for(chain_map_t::const_iterator chain_itr=tag_chain_itr->second.begin(); chain_itr!=tag_chain_itr->second.end(); @@ -248,7 +278,6 @@ TString energyBranchNameFromInvMassName(TString invMass_var){ return TString( (energyBranchNameFromInvMassName(std::string(invMass_var))).c_str()); } - int main(int argc, char **argv) { TStopwatch myClock; TStopwatch globalClock; @@ -294,6 +323,10 @@ int main(int argc, char **argv) { std::string minimType; std::vector branchList; + + + + //options for E/p std::string jsonFileName; std::string miscalibMap; @@ -320,7 +353,17 @@ int main(int argc, char **argv) { int nLoops; bool isDeadTriggerTower; std::string inputFileDeadXtal; - + std::string EBEE; + std::string EBEEpu; + int evtsPerPoint; + + int useRegression; + float yMIN; + float yMAX; + std::string dayMin; + std::string dayMax; + std::string dayZOOM; + std::string LUMI; //------------------------------ setting option categories po::options_description desc("Main options"); po::options_description outputOption("Output options"); @@ -329,6 +372,10 @@ int main(int argc, char **argv) { po::options_description smearerOption("Z smearer options"); po::options_description toyOption("toyMC options"); po::options_description EoverPOption("EoverP options"); + po::options_description laserMonitoringEPOption("laser monitoring with E/p options"); + po::options_description laserMonitoringEPvsPUOption("laser monitoring with E/p versus Pile Up options"); + + //po::options_description cmd_line_options; //cmd_line_options.add(desc).add(fitOption).add(smearOption); @@ -434,6 +481,27 @@ int main(int argc, char **argv) { ("constTermToy", po::value(&constTermToy)->default_value(0.01),"") ("eventsPerToy", po::value(&nEventsPerToy)->default_value(0),"=0: all events") ; + laserMonitoringEPOption.add_options() + ("laserMonitoringEP", "call the laser monitoring with E/p") + ("yMIN", po::value(&yMIN)->default_value(0.65),"y min") + ("yMAX", po::value(&yMAX)->default_value(1.05),"y max") + ("EBEE", po::value(&EBEE)->default_value("EB"),"barrel or endcap") + ("evtsPerPoint", po::value(&evtsPerPoint)->default_value(1000),"events per point") + ("useRegression", po::value(&useRegression)->default_value(0),"use regression") + ("dayMin", po::value(&dayMin)->default_value("1-7-2015"),"day min") + ("dayMax", po::value(&dayMax)->default_value("15-11-2015"),"day max") + ("dayZOOM", po::value(&dayZOOM)->default_value("10-8-2015"),"day ZOOM") + ("LUMI", po::value(&LUMI)->default_value("1.9"),"LUMI") + ; + laserMonitoringEPvsPUOption.add_options() + ("laserMonitoringEPvsPU", "call the laser monitoring with E/p") + // ("EBEEpu", po::value(&EBEEpu)->default_value("EB"),"barrel or endcap") + //("evtsPerPointpu", po::value(&evtsPerPoint)->default_value(1000),"events per point") + // ("useRegression", po::value(&useRegression)->default_value(0),"use regression") + // ("dayMin", po::value(&dayMin)->default_value("1-6-2015"),"day min") + // ("dayMax", po::value(&dayMax)->default_value("1-10-2015"),"day max") + // ("dayZOOM", po::value(&dayZOOM)->default_value("10-8-2015"),"day ZOOM") + ; EoverPOption.add_options() ("EOverPCalib", "call the E/p calibration") ("isBarrel", po::value(&isBarrel)->default_value(1),"1=barrel, 0=endcap") @@ -465,14 +533,16 @@ int main(int argc, char **argv) { ("isDeadTriggerTower", po::value(&isDeadTriggerTower)->default_value(false),"") ("inputFileDeadXtal", po::value(&inputFileDeadXtal)->default_value("NULL"),"") ; - + desc.add(inputOption); desc.add(outputOption); desc.add(fitterOption); desc.add(smearerOption); desc.add(toyOption); desc.add(EoverPOption); - + desc.add(laserMonitoringEPOption); + desc.add(laserMonitoringEPvsPUOption); + po::variables_map vm; // // po::store(po::parse_command_line(argc, argv, smearOption), vm); @@ -506,7 +576,8 @@ int main(int argc, char **argv) { if(!vm.count("regionsFile") && !vm.count("runDivide") && !vm.count("savePUTreeWeight") && - !vm.count("saveR9TreeWeight") && !vm.count("saveCorrEleTree") && !vm.count("EOverPCalib") + !vm.count("saveR9TreeWeight") && !vm.count("saveCorrEleTree") && !vm.count("EOverPCalib") && !vm.count("laserMonitoringEP") && + !vm.count("laserMonitoringEPvsPU") //&& !vm.count("saveRootMacro") ){ std::cerr << "[ERROR] Missing mandatory option \"regionsFile\"" << std::endl; @@ -536,19 +607,19 @@ int main(int argc, char **argv) { //============================== Check output folders bool checkDirectories=true; checkDirectories=checkDirectories && !system("[ -d "+TString(outDirFitResMC)+" ]"); - if(!checkDirectories && !vm.count("EOverPCalib")){ + if(!checkDirectories && !vm.count("EOverPCalib") && !vm.count("laserMonitoringEP") && !vm.count("laserMonitoringEPvsPU")){ std::cerr << "[ERROR] Directory " << outDirFitResMC << " not found" << std::endl; } checkDirectories=checkDirectories && !system("[ -d "+TString(outDirFitResData)+" ]"); - if(!checkDirectories && !vm.count("EOverPCalib")){ + if(!checkDirectories && !vm.count("EOverPCalib") && !vm.count("laserMonitoringEP") && !vm.count("laserMonitoringEPvsPU")){ std::cerr << "[ERROR] Directory " << outDirFitResData << " not found" << std::endl; } checkDirectories=checkDirectories && !system("[ -d "+TString(outDirImgMC)+" ]"); - if(!checkDirectories && !vm.count("EOverPCalib")){ + if(!checkDirectories && !vm.count("EOverPCalib") && !vm.count("laserMonitoringEP") && !vm.count("laserMonitoringEPvsPU")){ std::cerr << "[ERROR] Directory " << outDirImgMC << " not found" << std::endl; } checkDirectories=checkDirectories && !system("[ -d "+TString(outDirImgData)+" ]"); - if(!checkDirectories && !vm.count("EOverPCalib")){ + if(!checkDirectories && !vm.count("EOverPCalib") && !vm.count("laserMonitoringEP") && !vm.count("laserMonitoringEPvsPU")){ std::cerr << "[ERROR] Directory " << outDirImgData << " not found" << std::endl; } // checkDirectories=checkDirectories && !system("[ -d "+TString(outDirTable)+" ]"); @@ -562,6 +633,8 @@ int main(int argc, char **argv) { && !vm.count("saveR9TreeWeight") && !vm.count("saveRootMacro") && !vm.count("EOverPCalib") + && !vm.count("laserMonitoringEP") + && !vm.count("laserMonitoringEPvsPU") ) return 1; if(!dataPUFileName.empty()) dataPUFileNameVec.push_back(dataPUFileName.c_str()); @@ -628,6 +701,7 @@ int main(int argc, char **argv) { #endif } + //init chains and print for(tag_chain_map_t::const_iterator tag_chain_itr=tagChainMap.begin(); @@ -783,12 +857,13 @@ int main(int argc, char **argv) { } // end of data samples loop } // end of r9Weight - + //============================== // if(vm.count("dataPU")==0 && (tagChainMap["s"]).count("pileupHist")==0 && (tagChainMap["s"]).count("pileup")==0){ - if(vm.count("noPU")==0 && !vm.count("runToy")){ + if(vm.count("noPU")==0 && !vm.count("runToy") && !vm.count("laserMonitoringEP") && !vm.count("laserMonitoringEPvsPU") + && !vm.count("EOverPCalib")){ if(dataPUFileNameVec.empty() && (tagChainMap.count("s")!=0) && (tagChainMap["s"]).count("pileup")==0){ std::cerr << "[ERROR] Nor pileup mc tree configured in chain list file either dataPU histograms are not provided" << std::endl; return 1; @@ -991,6 +1066,8 @@ int main(int argc, char **argv) { } //end of sample loop } //end of branches loop + + //(tagChainMap["s"])["selected"]->GetEntries(); UpdateFriends(tagChainMap, regionsFileNameTag); @@ -1167,8 +1244,2494 @@ int main(int argc, char **argv) { mc = (tagChainMap["s"])["selected"]; } - if(vm.count("EOverPCalib") && vm.count("doEB")) { -///////// E/P calibration + + +//------------------------------ LASER MONITORING WITH E/P ------------------------------------------------------ + + if(vm.count("laserMonitoringEP")) { + + float timeLapse = 24.; // in hours + // int t1 = 1267401600; // 1 Mar 2010 + //int t2 = 1325289600; // 31 Dec 2011 + //int t1 = 1400000000; + //int t2 = 1600000000; + + // float yMIN = 0.65; + //float yMAX = 1.10; + + + // Set style options + setTDRStyle(); + gStyle->SetPadTickX(1); + gStyle->SetPadTickY(1); + gStyle->SetOptTitle(0); + gStyle->SetOptStat(1110); + gStyle->SetOptFit(1); + + // Set fitting options + TVirtualFitter::SetDefaultFitter("Fumili2"); + + + + //----------------- + // Input parameters + + + std::cout << "\n***************************************************************************************************************************" << std::endl; + + // std::string dayMin = ""; + //std::string dayMax = ""; + std::string dayMinLabel = ""; + std::string dayMaxLabel = ""; + std::string dayZOOMLabel =""; + float absEtaMin = -1.; + float absEtaMax = -1.; + int IetaMin = -1; + int IetaMax = -1; + int IphiMin = -1; + int IphiMax = -1; + + + + + int t1 = dateToInt(dayMin); + int t2 = dateToInt(dayMax); + int t3 = dateToInt(dayZOOM); + + /* + if(argc >= 5) + { + dayMin = std::string(argv[4])+" "+std::string(argv[5])+" "+std::string(argv[6]); + dayMax = std::string(argv[7])+" "+std::string(argv[8])+" "+std::string(argv[9]); + dayMinLabel = std::string(argv[4])+"_"+std::string(argv[5])+"_"+std::string(argv[6]); + dayMaxLabel = std::string(argv[7])+"_"+std::string(argv[8])+"_"+std::string(argv[9]); + + } + if(argc >= 11) + { + yMIN = atof(argv[10]); + yMAX = atof(argv[11]); + } + if(argc >= 13) + { + absEtaMin = atof(argv[12]); + absEtaMax = atof(argv[13]); + } + if(argc >= 15) + { + IetaMin = atoi(argv[14]); + IetaMax = atoi(argv[15]); + IphiMin = atoi(argv[16]); + IphiMax = atoi(argv[17]); + } + */ + + std::cout << "EBEE: " << EBEE << std::endl; + std::cout << "evtsPerPoint: " << evtsPerPoint << std::endl; + std::cout << "useRegression: " << useRegression << std::endl; + std::cout << "dayMin: " << dayMin << std::endl; + std::cout << "dayZOOM: " << dayZOOM << std::endl; + std::cout << "dayMax: " << dayMax << std::endl; + std::cout << "yMin: " << yMIN << std::endl; + std::cout << "yMax: " << yMAX << std::endl; + std::cout << "absEtaMin: " << absEtaMin << std::endl; + std::cout << "absEtaMax: " << absEtaMax << std::endl; + std::cout << "IetaMin: " << IetaMin << std::endl; + std::cout << "IetaMax: " << IetaMax << std::endl; + std::cout << "IphiMin: " << IphiMin << std::endl; + std::cout << "IphiMax: " << IphiMax << std::endl; + std::cout << "t1: " << t1 << std::endl; + std::cout << "t2: " << t2 << std::endl; + std::cout << "t3" << t3 << std::endl; + + std::string dayZOOM = ""; + std::string dayMin = ""; + std::string dayMax = ""; + + + //------------------- + // Define in/outfiles + + std::string folderName = std::string(EBEE) + "_" + dayMinLabel + "_" + dayMaxLabel; + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + char absEtaBuffer[50]; + sprintf(absEtaBuffer,"_%.2f-%.2f",absEtaMin,absEtaMax); + folderName += std::string(absEtaBuffer); + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + char absEtaBuffer[50]; + sprintf(absEtaBuffer,"_Ieta_%d-%d_Iphi_%d-%d",IetaMin,IetaMax,IphiMin,IphiMax); + folderName += std::string(absEtaBuffer); + } + + gSystem->mkdir(folderName.c_str()); + TFile* o = new TFile((folderName+"/"+folderName+"_histos.root").c_str(),"RECREATE"); + + + + // Get trees + std::cout << std::endl; + + /* + TChain* ntu_DA = new TChain("simpleNtupleEoverP/SimpleNtupleEoverP"); + FillChain(ntu_DA,"inputDATA.txt"); + std::cout << " DATA: " << std::setw(8) << ntu_DA->GetEntries() << " entries" << std::endl; + + TChain* ntu_MC = new TChain("simpleNtupleEoverP/SimpleNtupleEoverP"); + FillChain(ntu_MC,"inputMC.txt"); + std::cout << "REFERENCE: " << std::setw(8) << ntu_MC->GetEntries() << " entries" << std::endl; + */ + + if (data->GetEntries() == 0 || mc->GetEntries() == 0 ) + { + std::cout << "Error: At least one file is empty" << std::endl; + return -1; + } + + + + // Set branch addresses + int runNumber; + int runTime; + int nPU; + float avgLCSCEle[3], etaSCEle[3], phiSCEle[3], energySCEle[3], esEnergySCEle[3], pAtVtxGsfEle[3], energySCEle_corr[3]; + int seedXSCEle[3], seedYSCEle[3];//, seedZside; + // float seedLaserAlphaSCEle1; + + data->SetBranchStatus("*",0); + data->SetBranchStatus("runNumber",1); + data->SetBranchStatus("runTime",1); + data->SetBranchStatus("nPU",1); + data->SetBranchStatus("avgLCSCEle",1); + // data->SetBranchStatus("seedLaserAlphaSCEle1",1); + // data->SetBranchStatus("ele1_EOverP",1); + data->SetBranchStatus("etaSCEle",1); + data->SetBranchStatus("phiSCEle",1); + data->SetBranchStatus("energySCEle",1); + data->SetBranchStatus("energySCEle_corr",1); + data->SetBranchStatus("esEnergySCEle",1); + data->SetBranchStatus("pAtVtxGsfEle",1); + data->SetBranchStatus("seedXSCEle",1); + data->SetBranchStatus("seedYSCEle",1); + // data->SetBranchStatus("ele1_seedZside",1); + + data->SetBranchAddress("runNumber", &runNumber); + data->SetBranchAddress("runTime", &runTime); + data->SetBranchAddress("nPU", &nPU); + data->SetBranchAddress("avgLCSCEle", &avgLCSCEle[0]); + //data->SetBranchAddress("seedLaserAlphaSCEle1", &seedLaserAlphaSCEle1); + // data->SetBranchAddress("ele1_EOverP", &EoP); + data->SetBranchAddress("etaSCEle", &etaSCEle); + data->SetBranchAddress("phiSCEle", &phiSCEle); + if( useRegression < 1 ) + data->SetBranchAddress("energySCEle", &energySCEle); + else + data->SetBranchAddress("energySCEle_corr", &energySCEle_corr); + data->SetBranchAddress("energySCEle", &energySCEle); + data->SetBranchAddress("esEnergySCEle", &esEnergySCEle); + data->SetBranchAddress("pAtVtxGsfEle", &pAtVtxGsfEle); + data->SetBranchAddress("seedXSCEle", &seedXSCEle); + data->SetBranchAddress("seedYSCEle", &seedYSCEle); + // data->SetBranchAddress("ele1_seedZside", &seedZside); + + + mc->SetBranchStatus("*",0); + mc->SetBranchStatus("runNumber",1); + mc->SetBranchStatus("runTime",1); + mc->SetBranchStatus("nPU",1); + mc->SetBranchStatus("avgLCSCEle",1); + // mc->SetBranchStatus("seedLaserAlphaSCEle1",1); + // mc->SetBranchStatus("ele1_EOverP",1); + mc->SetBranchStatus("etaSCEle",1); + mc->SetBranchStatus("phiSCEle",1); + mc->SetBranchStatus("energySCEle",1); + mc->SetBranchStatus("energySCEle_corr",1); + mc->SetBranchStatus("esEnergySCEle",1); + mc->SetBranchStatus("pAtVtxGsfEle",1); + mc->SetBranchStatus("seedXSCEle",1); + mc->SetBranchStatus("seedYSCEle",1); + // mc->SetBranchStatus("ele1_seedZside",1); + + mc->SetBranchAddress("runNumber", &runNumber); + mc->SetBranchAddress("runTime", &runTime); + mc->SetBranchAddress("nPU", &nPU); + mc->SetBranchAddress("avgLCSCEle", &avgLCSCEle[0]); + //mc->SetBranchAddress("seedLaserAlphaSCEle1", &seedLaserAlphaSCEle1); + // mc->SetBranchAddress("ele1_EOverP", &EoP); + mc->SetBranchAddress("etaSCEle", &etaSCEle); + mc->SetBranchAddress("phiSCEle", &phiSCEle); + if( useRegression < 1 ) + mc->SetBranchAddress("energySCEle", &energySCEle); + else + mc->SetBranchAddress("energySCEle_corr", &energySCEle_corr); + mc->SetBranchAddress("energySCEle", &energySCEle); + mc->SetBranchAddress("esEnergySCEle", &esEnergySCEle); + mc->SetBranchAddress("pAtVtxGsfEle", &pAtVtxGsfEle); + mc->SetBranchAddress("seedXSCEle", &seedXSCEle); + mc->SetBranchAddress("seedYSCEle", &seedYSCEle); + // mc->SetBranchAddress("ele1_seedZside", &seedZside); + + + + + + + //-------------------------------------------------------- + // Define PU correction (to be used if useRegression == 0) + + // corr = p0 + p1 * nPU + float p0_EB; + float p1_EB; + float p0_EE; + float p1_EE; + + if( useRegression == 0 ) + { + //2012 EB + p0_EB = 0.9991; + p1_EB = 0.0001635; + //2012 EE + p0_EE = 0.9968; + p1_EE = 0.001046; + } + else + { + //2012 EB + p0_EB = 1.001; + p1_EB = -0.000143; + //2012 EE + p0_EE = 1.00327; + p1_EE = -0.000432; + } + + float p0 = -1.; + float p1 = -1.; + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + { + p0 = p0_EB; + p1 = p1_EB; + } + else + { + p0 = p0_EE; + p1 = p1_EE; + } + + //2015 + p0=1.; + p1=0.; + + + + + //--------------------------------- + // Build the reference distribution + + std::cout << std::endl; + std::cout << "***** Build reference for " << EBEE << " *****" << std::endl; + + TH1F* h_template = new TH1F("template", "", 2000, 0., 5.); + + for(int ientry = 0; ientry < mc->GetEntries(); ++ientry) + { + if( (ientry%10000 == 0) ) std::cout << "reading MC entry " << ientry << "\r" << std::endl;//std::flush; + mc->GetEntry(ientry); + + // selections + if( (strcmp(EBEE.c_str(),"EB") == 0) && (fabs(etaSCEle[0]) > 1.479) ) continue; // barrel + if( (strcmp(EBEE.c_str(),"EE") == 0) && (fabs(etaSCEle[0]) < 1.479 || fabs(etaSCEle[0])>2.5) ) continue; // endcap + + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + if( (fabs(etaSCEle[0]) < absEtaMin) || (fabs(etaSCEle[0]) > absEtaMax) ) continue; + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + if( (seedXSCEle[0] < IetaMin) || (seedXSCEle[0] > IetaMax) ) continue; + if( (seedYSCEle[0] < IphiMin) || (seedYSCEle[0] > IphiMax) ) continue; + } + + // PU correction + float PUCorr = (p0 + p1*nPU); + //std::cout << "p0: " << p0 << " p1: " << p1 << " nPU: " << nPU << std::endl; + + // fill the template histogram + h_template -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / PUCorr ); + } + + std::cout << "Reference built for " << EBEE << " - " << h_template->GetEntries() << " events" << std::endl; + + + + + + + //--------------------- + // Loop and sort events + + std::cout << std::endl; + std::cout << "***** Sort events and define bins *****" << std::endl; + + int nEntries = data -> GetEntriesFast(); + int nSavePts = 0; + std::vector isSavedEntries(nEntries); + std::vector sortedEntries; + std::vector timeStampFirst; + + for(int ientry = 0; ientry < nEntries; ++ientry) + { + data -> GetEntry(ientry); + isSavedEntries.at(ientry) = false; + + if( (ientry%10000 == 0) ) std::cout << "reading data entry " << ientry << "\r" << std::endl;//std::flush; + + // selections + if( (strcmp(EBEE.c_str(),"EB") == 0) && (fabs(etaSCEle[0]) > 1.479) ) continue; // barrel + if( (strcmp(EBEE.c_str(),"EE") == 0) && (fabs(etaSCEle[0]) < 1.479 || fabs(etaSCEle[0])>2.5) ) continue; // endcap + + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + if( (fabs(etaSCEle[0]) < absEtaMin) || (fabs(etaSCEle[0]) > absEtaMax) ) continue; + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + if( (seedXSCEle[0] < IetaMin) || (seedXSCEle[0] > IetaMax) ) continue; + if( (seedYSCEle[0] < IphiMin) || (seedYSCEle[0] > IphiMax) ) continue; + } + + if( runTime < t1 ) continue; + if( runTime > t2 ) continue; + + if( avgLCSCEle[0] <= 0. ) continue; + + isSavedEntries.at(ientry) = true; + + + // fill sorter + Sorter dummy; + dummy.time = runTime; + dummy.entry = ientry; + sortedEntries.push_back(dummy); + + ++nSavePts; + } + + // sort events + std::sort(sortedEntries.begin(),sortedEntries.end(),Sorter()); + std::cout << "Data sorted in " << EBEE << " - " << nSavePts << " events" << std::endl; + + std::map antiMap; + for(unsigned int iSaved = 0; iSaved < sortedEntries.size(); ++iSaved) + antiMap[sortedEntries.at(iSaved).entry] = iSaved; + + + //--------------------- + // Loop and define bins + + // "wide" bins - find events with time separation bigger than 1 day + int nWideBins = 1; + std::vector wideBinEntryMax; + int timeStampOld = -1; + + // TEventList* evlist=new TEventList("events"); + // data->Draw(">>events","","goff"); + // TEventList* evlist = (TEventList*) gDirectory->Get("events"); + + wideBinEntryMax.push_back(0); + for(int iSaved = 0; iSaved < nSavePts; ++iSaved) + { + if( iSaved%10000 == 0 ) std::cout << "reading saved entry " << iSaved << "\r" << std::endl;//std::flush; + data->GetEntry(sortedEntries[iSaved].entry); + // data->GetEntry(evlist->GetEntry(sortedEntries[iSaved].entry)); + + if( iSaved == 0 ) + { + timeStampOld = runTime; + continue; + } + + if( (runTime-timeStampOld)/3600. > timeLapse ) + { + ++nWideBins; + wideBinEntryMax.push_back(iSaved-1); + } + + timeStampOld = runTime; + } + std::cout << std::endl; + wideBinEntryMax.push_back(nSavePts); + + + // bins with approximatively evtsPerPoint events per bin + int nBins = 0; + std::vector binEntryMax; + + binEntryMax.push_back(0); + for(int wideBin = 0; wideBin < nWideBins; ++wideBin) + { + int nTempBins = std::max(1,int( (wideBinEntryMax.at(wideBin+1)-wideBinEntryMax.at(wideBin))/evtsPerPoint )); + int nTempBinEntries = int( (wideBinEntryMax.at(wideBin+1)-wideBinEntryMax.at(wideBin))/nTempBins ); + + for(int tempBin = 0; tempBin < nTempBins; ++tempBin) + { + ++nBins; + if( tempBin < nTempBins - 1 ) + binEntryMax.push_back( wideBinEntryMax.at(wideBin) + (tempBin+1)*nTempBinEntries ); + else + binEntryMax.push_back( wideBinEntryMax.at(wideBin+1) ); + } + } + + std::cout << "nBins = " << nBins << std::endl; + //for(int bin = 0; bin < nBins; ++bin) + // std::cout << "bin: " << bin + // << " entry min: " << setw(6) << binEntryMax.at(bin) + // << " entry max: " << setw(6) << binEntryMax.at(bin+1) + // << " events: " << setw(6) << binEntryMax.at(bin+1)-binEntryMax.at(bin) + // << std::endl; + + + + + + + //--------------------- + // histogram definition + + TH1F* h_scOccupancy_eta = new TH1F("h_scOccupancy_eta","", 298, -2.6, 2.6); + TH1F* h_scOccupancy_phi = new TH1F("h_scOccupancy_phi","", 363, -3.1765, 3.159); + SetHistoStyle(h_scOccupancy_eta); + SetHistoStyle(h_scOccupancy_phi); + + TH2F* h_seedOccupancy_EB = new TH2F("h_seedOccupancy_EB","", 171, -86., 85., 361, 0.,361.); + TH2F* h_seedOccupancy_EEp = new TH2F("h_seedOccupancy_EEp","", 101, 0.,101., 100, 0.,101.); + TH2F* h_seedOccupancy_EEm = new TH2F("h_seedOccupancy_EEm","", 101, 0.,101., 100, 0.,101.); + SetHistoStyle(h_seedOccupancy_EB); + SetHistoStyle(h_seedOccupancy_EEp); + SetHistoStyle(h_seedOccupancy_EEm); + + TH1F* h_EoP_spread = new TH1F("h_EoP_spread","",100,yMIN,yMAX); + TH1F* h_EoC_spread = new TH1F("h_EoC_spread","",100,yMIN,yMAX); + TH1F* h_EoP_spread_run = new TH1F("h_EoP_spread_run","",100,yMIN,yMAX); + TH1F* h_EoC_spread_run = new TH1F("h_EoC_spread_run","",100,yMIN,yMAX); + SetHistoStyle(h_EoP_spread,"EoP"); + SetHistoStyle(h_EoC_spread,"EoC"); + SetHistoStyle(h_EoP_spread_run,"EoP"); + SetHistoStyle(h_EoC_spread_run,"EoC"); + + TH1F* h_EoP_chi2 = new TH1F("h_EoP_chi2","",50,0.,5.); + TH1F* h_EoC_chi2 = new TH1F("h_EoC_chi2","",50,0.,5.); + SetHistoStyle(h_EoP_chi2,"EoP"); + SetHistoStyle(h_EoC_chi2,"EoC"); + + TH1F** h_EoP = new TH1F*[nBins]; + TH1F** h_EoC = new TH1F*[nBins]; + TH1F** h_Las = new TH1F*[nBins]; + TH1F** h_Tsp = new TH1F*[nBins]; + TH1F** h_Cvl = new TH1F*[nBins]; + + for(int i = 0; i < nBins; ++i) + { + char histoName[80]; + + sprintf(histoName, "EoP_%d", i); + h_EoP[i] = new TH1F(histoName, histoName, 2000, 0., 5.); + SetHistoStyle(h_EoP[i],"EoP"); + + sprintf(histoName, "EoC_%d", i); + h_EoC[i] = new TH1F(histoName, histoName, 2000, 0., 5.); + SetHistoStyle(h_EoC[i],"EoC"); + + sprintf(histoName, "Las_%d", i); + h_Las[i] = new TH1F(histoName, histoName, 500, 0.5, 1.5); + + sprintf(histoName, "Tsp_%d", i); + h_Tsp[i] = new TH1F(histoName, histoName, 500, 0.5, 1.5); + + } + + + // function definition + TF1** f_EoP = new TF1*[nBins]; + TF1** f_EoC = new TF1*[nBins]; + + + // graphs definition + TGraphAsymmErrors* g_fit = new TGraphAsymmErrors(); + TGraphAsymmErrors* g_c_fit = new TGraphAsymmErrors(); + + TGraphAsymmErrors* g_fit_run = new TGraphAsymmErrors(); + TGraphAsymmErrors* g_c_fit_run = new TGraphAsymmErrors(); + + TGraph* g_las = new TGraph(); + + TGraphErrors* g_LT = new TGraphErrors(); + + g_fit->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_fit->GetXaxis()->SetTimeDisplay(1); + g_c_fit->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_c_fit->GetXaxis()->SetTimeDisplay(1); + g_las->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_las->GetXaxis()->SetTimeDisplay(1); + g_LT->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_LT->GetXaxis()->SetTimeDisplay(1); + + + + + + + //------------------------------------ + // loop on the saved and sorted events + + std::cout << std::endl; + std::cout << "***** Fill and fit histograms *****" << std::endl; + + std::vector Entries(nBins); + std::vector AveTime(nBins); + std::vector MinTime(nBins); + std::vector MaxTime(nBins); + std::vector AveRun(nBins); + std::vector MinRun(nBins); + std::vector MaxRun(nBins); + std::vector AveLT(nBins); + std::vector AveLT2(nBins); + + int iSaved = -1; + for(int ientry = 0; ientry < nEntries; ++ientry) + { + if( (ientry%100000 == 0) ) std::cout << "reading entry " << ientry << "\r" << std::endl;//std::flush; + + if( isSavedEntries.at(ientry) == false ) continue; + + ++iSaved; + + int iSaved = antiMap[ientry]; + int bin = -1; + for(bin = 0; bin < nBins; ++bin) + if( iSaved >= binEntryMax.at(bin) && iSaved < binEntryMax.at(bin+1) ) + break; + + //std::cout << "bin = " << bin << " iSaved = "<< iSaved << std::endl; + data->GetEntry(ientry); + + + + Entries[bin] += 1; + + if( iSaved == binEntryMax.at(bin)+1 ) MinTime[bin] = runTime; + if( iSaved == binEntryMax.at(bin+1)-1 ) MaxTime[bin] = runTime; + AveTime[bin] += runTime; + + if( iSaved == binEntryMax.at(bin)+1 ) MinRun[bin] = runNumber; + if( iSaved == binEntryMax.at(bin+1)-1 ) MaxRun[bin] = runNumber; + AveRun[bin] += runNumber; + + // float LT = (-1. / seedLaserAlphaSCEle1 * log(avgLCSCEle[0])); + float LT = 1.; + AveLT[bin] += LT; + AveLT2[bin] += LT*LT; + + // PU correction + float PUCorr = (p0 + p1*nPU); + + // fill the histograms + (h_EoP[bin]) -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / avgLCSCEle[0] / PUCorr); + (h_EoC[bin]) -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / PUCorr ); + + (h_Las[bin]) -> Fill(avgLCSCEle[0]); + (h_Tsp[bin]) -> Fill(1./avgLCSCEle[0]); + + h_scOccupancy_eta -> Fill(etaSCEle[0]); + h_scOccupancy_phi -> Fill(phiSCEle[0]); + if(fabs(etaSCEle[0])<1.449) + h_seedOccupancy_EB -> Fill(seedXSCEle[0],seedYSCEle[0]); + else if(etaSCEle[0]>1.449) + h_seedOccupancy_EEp -> Fill(seedXSCEle[0],seedYSCEle[0]); + else if(etaSCEle[0]<-1.449) + h_seedOccupancy_EEm -> Fill(seedXSCEle[0],seedYSCEle[0]); + } + + for(int bin = 0; bin < nBins; ++bin) + { + AveTime[bin] = 1. * AveTime[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveRun[bin] = 1. * AveRun[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveLT[bin] = 1. * AveLT[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveLT2[bin] = 1. * AveLT2[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + //std::cout << date << " " << AveTime[i] << " " << MinTime[i] << " " << MaxTime[i] << std::endl; + } + + + + + + + int rebin = 2; + if( strcmp(EBEE.c_str(),"EE") == 0 ) rebin *= 2; + + h_template -> Rebin(rebin); + + + + float EoP_scale = 0.; + float EoP_err = 0.; + int EoP_nActiveBins = 0; + + float EoC_scale = 0.; + float EoC_err = 0.; + int EoC_nActiveBins = 0; + + float LCInv_scale = 0; + + std::vector validBins; + for(int i = 0; i < nBins; ++i) + { + bool isValid = true; + + h_EoP[i] -> Rebin(rebin); + h_EoC[i] -> Rebin(rebin); + + + + //------------------------------------ + // Fill the graph for uncorrected data + + // define the fitting function + // N.B. [0] * ( [1] * f( [1]*(x-[2]) ) ) + + //o -> cd(); + char convolutionName[50]; + sprintf(convolutionName,"h_convolution_%d",i); + //h_Cvl[i] = ConvoluteTemplate(std::string(convolutionName),h_template,h_Las[i],32768,-5.,5.); + h_Cvl[i] = MellinConvolution(std::string(convolutionName),h_template,h_Tsp[i]); + + histoFunc* templateHistoFunc = new histoFunc(h_template); + histoFunc* templateConvolutedHistoFunc = new histoFunc(h_Cvl[i]); + char funcName[50]; + + sprintf(funcName,"f_EoP_%d",i); + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + f_EoP[i] = new TF1(funcName, templateConvolutedHistoFunc, 0.8*(h_Tsp[i]->GetMean()), 1.4*(h_Tsp[i]->GetMean()), 3, "histoFunc"); + else + f_EoP[i] = new TF1(funcName, templateConvolutedHistoFunc, 0.75*(h_Tsp[i]->GetMean()), 1.5*(h_Tsp[i]->GetMean()), 3, "histoFunc"); + + f_EoP[i] -> SetParName(0,"Norm"); + f_EoP[i] -> SetParName(1,"Scale factor"); + f_EoP[i] -> SetLineWidth(1); + f_EoP[i] -> SetNpx(10000); + + double xNorm = h_EoP[i]->GetEntries()/h_template->GetEntries() * + h_EoP[i]->GetBinWidth(1)/h_template->GetBinWidth(1); + + f_EoP[i] -> FixParameter(0, xNorm); + f_EoP[i] -> SetParameter(1, 1.); + f_EoP[i] -> FixParameter(2, 0.); + f_EoP[i] -> SetLineColor(kRed+2); + + int fStatus = 0; + int nTrials = 0; + TFitResultPtr rp; + + rp = h_EoP[i] -> Fit(funcName, "ERLS+"); + while( (fStatus != 0) && (nTrials < 10) ) + { + rp = h_EoP[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + if(fStatus == 0) break; + ++nTrials; + } + + // fill the graph + double eee = f_EoP[i]->GetParError(1); + //float k = f_EoP[i]->GetParameter(1); + float k = f_EoP[i]->GetParameter(1) / h_Tsp[i]->GetMean(); //needed when using mellin's convolution + + /* + std::cout << i <<"--nocorr---- "<< 1./k << std::endl; + std::cout <<" condizione 1: " << h_EoP[i]->GetEntries() << " fStatus: " << fStatus << " eee: " << eee << "con eee che ci piace essere maggiore di : " << 0.05*h_template->GetRMS()/sqrt(evtsPerPoint) << std::endl ; + getchar(); + */ + + + // if( (h_EoP[i]->GetEntries() > 3) && (fStatus == 0) && (eee > 0.05*h_template->GetRMS()/sqrt(evtsPerPoint)) ) + if( (h_EoP[i]->GetEntries() > 500) && (fStatus == 0) ) + { + float date = (float)AveTime[i]; + float dLow = (float)(AveTime[i]-MinTime[i]); + float dHig = (float)(MaxTime[i]-AveTime[i]); + float run = (float)AveRun[i]; + float rLow = (float)(AveRun[i]-MinRun[i]); + float rHig = (float)(MaxRun[i]-AveRun[i]); + + g_fit -> SetPoint(i, date , 1./k); + g_fit -> SetPointError(i, dLow , dHig, eee/k/k, eee/k/k); + + g_fit_run -> SetPoint(i, run , 1./k); + g_fit_run -> SetPointError(i, rLow , rHig, eee/k/k, eee/k/k); + + std::cout <<"************-------------------*****************" < Fill(f_EoP[i]->GetChisquare()/f_EoP[i]->GetNDF()); + + EoP_scale += 1./k; + EoP_err += eee/k/k; + ++EoP_nActiveBins; + } + else + { + std::cout << "Fitting uncorrected time bin: " << i << " Fail status: " << fStatus << " sigma: " << eee << std::endl; + isValid = false; + } + + //---------------------------------- + // Fill the graph for corrected data + + // define the fitting function + // N.B. [0] * ( [1] * f( [1]*(x-[2]) ) ) + + sprintf(funcName,"f_EoC_%d",i); + if( strcmp(EBEE.c_str(),"EB") == 0 ) + f_EoC[i] = new TF1(funcName, templateHistoFunc, 0.8, 1.4, 3, "histoFunc"); + else + f_EoC[i] = new TF1(funcName, templateHistoFunc, 0.75, 1.5, 3, "histoFunc"); + f_EoC[i] -> SetParName(0,"Norm"); + f_EoC[i] -> SetParName(1,"Scale factor"); + f_EoC[i] -> SetLineWidth(1); + f_EoC[i] -> SetNpx(10000); + + xNorm = h_EoC[i]->GetEntries()/h_template->GetEntries() * + h_EoC[i]->GetBinWidth(1)/h_template->GetBinWidth(1); + + f_EoC[i] -> FixParameter(0, xNorm); + f_EoC[i] -> SetParameter(1, 0.99); + f_EoC[i] -> FixParameter(2, 0.); + f_EoC[i] -> SetLineColor(kGreen+2); + + + rp = h_EoC[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + nTrials = 0; + while( (fStatus != 0) && (nTrials < 10) ) + { + rp = h_EoC[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + if(fStatus == 0) break; + ++nTrials; + } + + + // fill the graph + k = f_EoC[i]->GetParameter(1); + eee = f_EoC[i]->GetParError(1); + + /* std::cout << i <<"--corr---- "<< 1./k << std::endl; + std::cout <<" condizione 1: " << h_EoP[i]->GetEntries() << " fStatus: " << fStatus << " eee: " << eee << "con eee che ci piace essere maggiore di : " << 0.05*h_template->GetRMS()/sqrt(evtsPerPoint) << std::endl ; + getchar(); + */ + + + if( (h_EoC[i]->GetEntries() > 500) && (fStatus == 0) ) + // if( (h_EoC[i]->GetEntries() > 10) && (fStatus == 0) && (eee > 0.05*h_template->GetRMS()/sqrt(evtsPerPoint)) ) + { + float date = (float)AveTime[i]; + float dLow = (float)(AveTime[i]-MinTime[i]); + float dHig = (float)(MaxTime[i]-AveTime[i]); + float run = (float)AveRun[i]; + float rLow = (float)(AveRun[i]-MinRun[i]); + float rHig = (float)(MaxRun[i]-AveRun[i]); + + g_c_fit -> SetPoint(i, date , 1./k); + g_c_fit -> SetPointError(i, dLow , dHig , eee/k/k, eee/k/k); + + g_c_fit_run -> SetPoint(i, run , 1./k); + g_c_fit_run -> SetPointError(i, rLow , rHig, eee/k/k, eee/k/k); + std::cout <<"************-------------------*****************" < Fill(f_EoC[i]->GetChisquare()/f_EoP[i]->GetNDF()); + + EoC_scale += 1./k; + EoC_err += eee/k/k; + ++EoC_nActiveBins; + } + else + { + std::cout << "Fitting corrected time bin: " << i << " Fail status: " << fStatus << " sigma: " << eee << std::endl; + isValid = false; + } + + if( isValid == true ) validBins.push_back(i); + } + + EoP_scale /= EoP_nActiveBins; + EoP_err /= EoP_nActiveBins; + + EoC_scale /= EoC_nActiveBins; + EoC_err /= EoC_nActiveBins; + + + + + + + //---------------------------------------- + // Fill the graph for avg laser correction + + //fede + for(unsigned int itr = 0; itr < validBins.size(); ++itr) + { + int i = validBins.at(itr); + g_las -> SetPoint(itr, (float)AveTime[i], h_Tsp[i]->GetMean()); + g_LT -> SetPoint(itr, (float)AveTime[i], AveLT[i] ); + g_LT -> SetPointError(itr, 0., sqrt(AveLT2[i]-AveLT[i]*AveLT[i]) / sqrt(Entries[i]) ); + + LCInv_scale += h_Tsp[i]->GetMean(); + } + + LCInv_scale /= validBins.size(); + + + + + + + //--------------- + // Rescale graphs + + float yscale = 1.; + //float yscale = 1./EoC_scale; + + for(unsigned int itr = 0; itr < validBins.size(); ++itr) + { + double x,y; + g_fit -> GetPoint(itr,x,y); + g_fit -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoP_spread -> Fill(y*yscale); + + g_c_fit -> GetPoint(itr,x,y); + g_c_fit -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoC_spread -> Fill(y*yscale); + + g_fit_run -> GetPoint(itr,x,y); + g_fit_run -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoP_spread_run -> Fill(y*yscale); + + g_c_fit_run -> GetPoint(itr,x,y); + g_c_fit_run -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoC_spread_run -> Fill(y*yscale); + + g_las -> GetPoint(itr,x,y); + g_las -> SetPoint(itr,x,y*yscale*EoP_scale/LCInv_scale); + } + TF1 EoC_pol0("EoC_pol0","pol0",t1,t2); + EoC_pol0.SetLineColor(kGreen+2); + EoC_pol0.SetLineWidth(2); + EoC_pol0.SetLineStyle(2); + g_c_fit -> Fit("EoC_pol0","QNR"); + + + + + + + + + //---------------------------- + // Print out global quantities + + std::cout << std::endl; + std::cout << "***** Mean scales and errors *****" << std::endl; + std::cout << std::fixed; + std::cout << std::setprecision(4); + std::cout << "Mean EoP scale: " << std::setw(6) << EoP_scale << " mean EoP error: " << std::setw(8) << EoP_err << std::endl; + std::cout << "Mean EoC scale: " << std::setw(6) << EoC_scale << " mean EoC error: " << std::setw(8) << EoC_err << std::endl; + std::cout << "Mean 1/LC scale: " << std::setw(6) << LCInv_scale << std::endl; + + + + + + + //----------------- + // Occupancy plots + //----------------- + + TCanvas* c_scOccupancy = new TCanvas("c_scOccupancy","SC occupancy",0,0,1000,500); + c_scOccupancy -> Divide(2,1); + + c_scOccupancy -> cd(1); + h_scOccupancy_eta -> GetXaxis() -> SetTitle("sc #eta"); + h_scOccupancy_eta -> GetYaxis() -> SetTitle("events"); + h_scOccupancy_eta -> Draw(); + + c_scOccupancy -> cd(2); + h_scOccupancy_phi -> GetXaxis() -> SetTitle("sc #phi"); + h_scOccupancy_phi -> GetYaxis() -> SetTitle("events"); + h_scOccupancy_phi -> Draw(); + + TCanvas* c_seedOccupancy = new TCanvas("c_seedOccupancy","seed occupancy",0,0,1500,500); + c_seedOccupancy -> Divide(3,1); + + c_seedOccupancy -> cd(1); + h_seedOccupancy_EB -> GetXaxis() -> SetTitle("seed i#eta"); + h_seedOccupancy_EB -> GetYaxis() -> SetTitle("seed i#phi"); + h_seedOccupancy_EB -> Draw("COLZ"); + + c_seedOccupancy -> cd(2); + h_seedOccupancy_EEp -> GetXaxis() -> SetTitle("seed ix"); + h_seedOccupancy_EEp -> GetYaxis() -> SetTitle("seed iy"); + h_seedOccupancy_EEp -> Draw("COLZ"); + + c_seedOccupancy -> cd(3); + h_seedOccupancy_EEm -> GetXaxis() -> SetTitle("seed ix"); + h_seedOccupancy_EEm -> GetYaxis() -> SetTitle("seed iy"); + h_seedOccupancy_EEm -> Draw("COLZ"); + + + + //----------- + // Chi2 plots + //----------- + + TCanvas* c_chi2 = new TCanvas("c_chi2","fit chi2",0,0,500,500); + c_chi2 -> cd(); + + h_EoC_chi2 -> GetXaxis() -> SetTitle("#chi^{2}/N_{dof}"); + h_EoC_chi2 -> Draw(""); + gPad -> Update(); + + TPaveStats* s_EoC = new TPaveStats; + s_EoC = (TPaveStats*)(h_EoC_chi2->GetListOfFunctions()->FindObject("stats")); + s_EoC -> SetStatFormat("1.4g"); + s_EoC -> SetTextColor(kGreen+2); + s_EoC->SetY1NDC(0.59); + s_EoC->SetY2NDC(0.79); + s_EoC -> Draw("sames"); + gPad -> Update(); + + h_EoP_chi2 -> GetXaxis() -> SetTitle("#chi^{2}/N_{dof}"); + h_EoP_chi2 -> Draw("sames"); + gPad -> Update(); + + TPaveStats* s_EoP = new TPaveStats; + s_EoP = (TPaveStats*)(h_EoP_chi2->GetListOfFunctions()->FindObject("stats")); + s_EoP -> SetStatFormat("1.4g"); + s_EoP -> SetTextColor(kRed+2); + s_EoP->SetY1NDC(0.79); + s_EoP->SetY2NDC(0.99); + s_EoP -> Draw("sames"); + gPad -> Update(); + + //ciao + //------------------- + // RMS vs Num evts -BARREL + //------------------- + double x[13]={2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.,24.,30.}; + double y[13]={0.001049,0.001114,0.0009367,0.0008480,0.0007669,0.0007892,0.0006699,0.0006473,0.0006235,0.0005903,0.0005815,0.0005459,0.0005506}; + + TCanvas* RMSeb = new TCanvas("plot", "plot",0,0,500,500); + TGraph* gRMSeb = new TGraph(13,x,y); + + gRMSeb->Draw("APC"); + gRMSeb -> SetMarkerColor(38); + gRMSeb -> SetLineColor(38); + gRMSeb->GetXaxis()->SetTitle("Number of Events - Barrel"); + gRMSeb->GetYaxis()->SetTitle("RMS"); + + RMSeb -> Print((folderName+"/"+folderName+"_RMSeb"+".png").c_str(),"png"); + RMSeb -> Print((folderName+"/"+folderName+"_RMSeb"+".pdf").c_str(),"pdf"); + + //------------------- + // RMS vs Num evts -ENDCAP + //------------------- + + double xx[11]={2.,4.,6.,8.,10.,12.,14.,16.,18.,20.,22.}; + double yy[11]={0.007234,0.005759,0.004174,0.004255,0.003833,0.004037,0.003912,0.004251,0.003598,0.004067,0.004138}; + + TCanvas* RMSee = new TCanvas("plot", "plot",0,0,500,500); + TGraph* gRMSee = new TGraph(11,xx,yy); + + gRMSee->Draw("APC"); + gRMSee -> SetMarkerColor(38); + gRMSee -> SetLineColor(38); + gRMSee->GetXaxis()->SetTitle("Number of Events - Endcap"); + gRMSee->GetYaxis()->SetTitle("RMS"); + + + RMSee -> Print((folderName+"/"+folderName+"_RMSee"+".png").c_str(),"png"); + RMSee -> Print((folderName+"/"+folderName+"_RMSee"+".pdf").c_str(),"pdf"); + + + //ciao + //------------------- + // histos + //------------------- + + + /* + for ( int i = 0; i < nBins; ++i) + { + + TCanvas* histoEoP = new TCanvas("histo","histo",0,0,500,500); + histoEoP -> cd(); + + h_EoP[i] -> Draw(); + f_EoP[i] -> SetLineWidth(2); + f_EoP[i] -> SetLineColor(4); + f_EoP[i] -> Draw("same"); + // histoEoP -> Update(); + + histoEoP -> Print((folderName+"/"+folderName+"_histoEoP"+std::to_string(i)+".png").c_str(),"png"); + histoEoP -> Print((folderName+"/"+folderName+"_histoEoP"+std::to_string(i)+".pdf").c_str(),"pdf"); + + + TCanvas* histoEoC = new TCanvas("histo","histo",0,0,500,500); + histoEoC -> cd(); + + h_EoC[i] -> Draw(); + f_EoC[i] -> SetLineWidth(2); + f_EoC[i] -> SetLineColor(4); + f_EoC[i] -> Draw("same"); + //histoEoC -> Update(); + + histoEoC -> Print((folderName+"/"+folderName+"_histoEoC"+to_string(i)+".png").c_str(),"png"); + histoEoC -> Print((folderName+"/"+folderName+"_histoEoC"+to_string(i)+".pdf").c_str(),"pdf"); + + } + */ + + //------------------- + // Final Plot vs date + //------------------- + + TCanvas* cplot = new TCanvas("cplot", "history plot vs date",100,100,1000,500); + cplot->cd(); + + TPad *cLeft = new TPad("pad_0","pad_0",0.00,0.00,0.75,1.00); + TPad *cRight = new TPad("pad_1","pad_1",0.75,0.00,1.00,1.00); + + cLeft->SetLeftMargin(0.15); + cLeft->SetRightMargin(0.025); + cRight->SetLeftMargin(0.025); + + cLeft->Draw(); + cRight->Draw(); + + float tYoffset = 1.0; + float labSize = 0.05; + float labSize2 = 0.06; + + cLeft->cd(); + + cLeft->SetGridx(); + cLeft->SetGridy(); + + TH1F *hPad = (TH1F*)gPad->DrawFrame(t1,0.9,t2,1.05); + + //hPad->GetXaxis()->SetLimits(t3,t2); + hPad->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + hPad->GetXaxis()->SetTimeDisplay(1); + hPad->GetXaxis() -> SetRangeUser(MinTime[0]-43200,MaxTime[nBins-1]+43200); + hPad->GetXaxis()->SetTitle("date (day/month)"); + //ciao + //hPad->GetXaxis()->SetLabelSize(0.025); + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + else + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + hPad->GetYaxis()->SetTitleOffset(tYoffset); + hPad->GetXaxis()->SetLabelSize(0.03); + hPad->GetXaxis()->SetTitleSize(labSize2); + hPad->GetYaxis()->SetLabelSize(labSize); + hPad->GetYaxis()->SetTitleSize(labSize2); + hPad -> SetMinimum(yMIN); + hPad -> SetMaximum(yMAX); + + // draw history plot + g_fit -> SetMarkerStyle(24); + g_fit -> SetMarkerSize(0.7); + g_fit -> SetMarkerColor(kRed+2); + g_fit -> SetLineColor(kRed+2); + g_fit -> Draw("P"); + g_c_fit -> SetMarkerStyle(20); + g_c_fit -> SetMarkerColor(kGreen+2); + g_c_fit -> SetLineColor(kGreen+2); + g_c_fit -> SetMarkerSize(0.7); + g_c_fit -> Draw("EP,same"); + g_las -> SetLineColor(kAzure-2); + g_las -> SetLineWidth(2); + g_las -> Draw("L,same"); + + TLegend* legend = new TLegend(0.60,0.78,0.90,0.94); + legend -> SetLineColor(kWhite); + legend -> SetLineWidth(0); + legend -> SetFillColor(kWhite); + legend -> SetFillStyle(0); + legend -> SetTextFont(42); + legend -> SetTextSize(0.04); + legend -> AddEntry(g_c_fit,"with LM correction","PL"); + legend -> AddEntry(g_fit, "without LM correction","PL"); + legend -> AddEntry(g_las, "1 / LM correction","L"); + legend -> Draw("same"); + + char latexBuffer[250]; + + sprintf(latexBuffer,"CMS 2015 Preliminary"); + TLatex* latex = new TLatex(0.18,0.89,latexBuffer); + latex -> SetNDC(); + latex -> SetTextFont(62); + latex -> SetTextSize(0.05); + latex -> Draw("same"); + + //sprintf(latexBuffer,"#sqrt{s} = 8 TeV L = 3.95 fb^{-1}"); + sprintf(latexBuffer,"#sqrt{s} = 13 TeV, L =%s fb^{-1} ", LUMI.c_str()); + + // sprintf(latexBuffer, LUMI.c_str()); + TLatex* latex2 = new TLatex(0.18,0.84,latexBuffer); + latex2 -> SetNDC(); + latex2 -> SetTextFont(42); + latex2 -> SetTextSize(0.05); + latex2 -> Draw("same"); + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + sprintf(latexBuffer,"ECAL Barrel"); + else + sprintf(latexBuffer,"ECAL Endcap"); + TLatex* latex3 = new TLatex(0.18,0.19,latexBuffer); + latex3 -> SetNDC(); + latex3 -> SetTextFont(42); + latex3 -> SetTextSize(0.05); + latex3 -> Draw("same"); + + //sprintf(latexBuffer,"%.2E events",1.*nSavePts); + //TLatex* latex4 = new TLatex(0.18,0.24,latexBuffer); + //latex4 -> SetNDC(); + //latex4 -> SetTextFont(42); + //latex4 -> SetTextSize(0.04); + //latex4 -> Draw("same"); + // + //sprintf(latexBuffer,"%d events/bin - %d bins",evtsPerPoint,nBins); + //TLatex* latex5 = new TLatex(0.18,0.19,latexBuffer); + //latex5 -> SetNDC(); + //latex5 -> SetTextFont(42); + //latex5 -> SetTextSize(0.04); + //latex5 -> Draw("same"); + + + cRight -> cd(); + + TPaveStats* s_EoP_spread = new TPaveStats(); + TPaveStats* s_EoC_spread = new TPaveStats(); + + + h_EoC_spread -> SetFillStyle(3001); + h_EoC_spread -> SetFillColor(kGreen+2); + h_EoC_spread->GetYaxis()->SetLabelSize(0.09); + h_EoC_spread->GetYaxis()->SetLabelOffset(-0.03); + h_EoC_spread->GetYaxis()->SetTitleSize(0.08); + h_EoC_spread->GetYaxis()->SetNdivisions(505); + h_EoC_spread->GetXaxis()->SetLabelOffset(1000); + + h_EoC_spread -> Draw("hbar"); + gPad -> Update(); + + s_EoC_spread = (TPaveStats*)(h_EoC_spread->GetListOfFunctions()->FindObject("stats")); + s_EoC_spread -> SetStatFormat("1.4g"); + s_EoC_spread->SetX1NDC(0.06); //new x start position + s_EoC_spread->SetX2NDC(0.71); //new x end position + s_EoC_spread->SetY1NDC(0.43); //new x start position + s_EoC_spread->SetY2NDC(0.34); //new x end position + s_EoC_spread -> SetOptStat(1100); + s_EoC_spread ->SetTextColor(kGreen+2); + s_EoC_spread ->SetTextSize(0.08); + s_EoC_spread -> Draw("sames"); + + + h_EoP_spread -> SetFillStyle(3001); + h_EoP_spread -> SetFillColor(kRed+2); + h_EoP_spread -> Draw("hbarsames"); + gPad -> Update(); + s_EoP_spread = (TPaveStats*)(h_EoP_spread->GetListOfFunctions()->FindObject("stats")); + s_EoP_spread -> SetStatFormat("1.4g"); + s_EoP_spread->SetX1NDC(0.06); //new x start position + s_EoP_spread->SetX2NDC(0.71); //new x end position + s_EoP_spread->SetY1NDC(0.33); //new x start position + s_EoP_spread->SetY2NDC(0.24); //new x end position + s_EoP_spread ->SetOptStat(1100); + s_EoP_spread ->SetTextColor(kRed+2); + s_EoP_spread ->SetTextSize(0.08); + s_EoP_spread -> Draw("sames"); + + /* + h_EoP_spread -> SetFillStyle(3001); + h_EoP_spread -> SetFillColor(kRed+2); + h_EoP_spread -> Draw("hbarsame"); + gPad -> Update(); + */ + + + + //------------------ + // Final plot vs run + //------------------ + + TCanvas* cplot_run = new TCanvas("cplot_run", "history plot vs run",100,100,1000,500); + cplot_run->cd(); + + cLeft = new TPad("pad_0_run","pad_0_run",0.00,0.00,0.75,1.00); + cRight = new TPad("pad_1_run","pad_1_run",0.75,0.00,1.00,1.00); + + cLeft->SetLeftMargin(0.15); + cLeft->SetRightMargin(0.025); + cRight->SetLeftMargin(0.025); + + cLeft->Draw(); + cRight->Draw(); + + tYoffset = 1.5; + labSize = 0.04; + labSize2 = 0.07; + + cLeft->cd(); + + cLeft->SetGridx(); + cLeft->SetGridy(); + + hPad = (TH1F*)gPad->DrawFrame(MinRun[0]-1000,0.9,MaxRun[nBins-1]+1000,1.05); + hPad->GetXaxis()->SetTitle("run"); + if( strcmp(EBEE.c_str(),"EB") == 0 ) + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + else + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + hPad->GetYaxis()->SetTitleOffset(tYoffset); + hPad->GetXaxis()->SetLabelSize(labSize); + hPad->GetXaxis()->SetTitleSize(labSize); + hPad->GetYaxis()->SetLabelSize(labSize); + hPad->GetYaxis()->SetTitleSize(labSize); + hPad -> SetMinimum(yMIN); + hPad -> SetMaximum(yMAX); + + // draw history plot + g_fit_run -> SetMarkerStyle(20); + g_fit_run -> SetMarkerSize(0.7); + g_fit_run -> SetMarkerColor(kRed+2); + g_fit_run -> SetLineColor(kRed+2); + g_fit_run -> Draw("P"); + g_c_fit_run -> SetMarkerStyle(20); + g_c_fit_run -> SetMarkerColor(kGreen+2); + g_c_fit_run -> SetLineColor(kGreen+2); + g_c_fit_run -> SetMarkerSize(0.7); + g_c_fit_run -> Draw("P,same"); + + + cRight -> cd(); + + s_EoP_spread = new TPaveStats(); + s_EoC_spread = new TPaveStats(); + + + h_EoC_spread_run -> SetFillStyle(3001); + h_EoC_spread_run -> SetFillColor(kGreen+2); + h_EoC_spread_run->GetYaxis()->SetLabelSize(labSize2); + h_EoC_spread_run->GetYaxis()->SetTitleSize(labSize2); + h_EoC_spread_run->GetYaxis()->SetNdivisions(505); + h_EoC_spread_run->GetYaxis()->SetLabelOffset(-0.02); + h_EoC_spread_run->GetXaxis()->SetLabelOffset(1000); + + h_EoC_spread_run -> Draw("hbar"); + gPad -> Update(); + + s_EoC_spread = (TPaveStats*)(h_EoC_spread_run->GetListOfFunctions()->FindObject("stats")); + s_EoC_spread ->SetTextColor(kGreen+2); + s_EoC_spread ->SetTextSize(0.06); + s_EoC_spread->SetX1NDC(0.49); //new x start position + s_EoC_spread->SetX2NDC(0.99); //new x end position + s_EoC_spread->SetY1NDC(0.475); //new x start position + s_EoC_spread->SetY2NDC(0.590); //new x end position + s_EoC_spread -> SetOptStat(1100); + s_EoC_spread -> Draw("sames"); + + h_EoP_spread_run -> SetFillStyle(3001); + h_EoP_spread_run -> SetFillColor(kRed+2); + h_EoP_spread_run -> Draw("hbarsames"); + gPad -> Update(); + + s_EoP_spread = (TPaveStats*)(h_EoP_spread_run->GetListOfFunctions()->FindObject("stats")); + s_EoP_spread->SetX1NDC(0.49); //new x start position + s_EoP_spread->SetX2NDC(0.99); //new x end position + s_EoP_spread->SetY1NDC(0.350); //new x start position + s_EoP_spread->SetY2NDC(0.475); //new x end position + s_EoP_spread ->SetOptStat(1100); + s_EoP_spread ->SetTextColor(kRed+2); + s_EoP_spread ->SetTextSize(0.06); + s_EoP_spread -> Draw("sames"); + + + + + c_chi2 -> Print((folderName+"/"+folderName+"_fitChi2.png").c_str(),"png"); + c_scOccupancy -> Print((folderName+"/"+folderName+"_scOccupancy.png").c_str(),"png"); + c_seedOccupancy -> Print((folderName+"/"+folderName+"_seedOccupancy.png").c_str(),"png"); + cplot -> Print((folderName+"/"+folderName+"_history_vsTime.png").c_str(),"png"); + cplot_run -> Print((folderName+"/"+folderName+"_history_vsRun.png").c_str(),"png"); + + c_chi2 -> Print((folderName+"/"+folderName+"_fitChi2.pdf").c_str(),"pdf"); + c_scOccupancy -> Print((folderName+"/"+folderName+"_scOccupancy.pdf").c_str(),"pdf"); + c_seedOccupancy -> Print((folderName+"/"+folderName+"_seedOccupancy.pdf").c_str(),"pdf"); + cplot -> Print((folderName+"/"+folderName+"_history_vsTime.pdf").c_str(),"pdf"); + cplot_run -> Print((folderName+"/"+folderName+"_history_vsRun.pdf").c_str(),"pdf"); + + cplot -> SaveAs((folderName+"/"+folderName+"_history_vsTime.C").c_str()); + cplot_run -> SaveAs((folderName+"/"+folderName+"_history_vsRun.C").c_str()); + + cplot -> SaveAs((folderName+"/"+folderName+"_history_vsTime.root").c_str()); + cplot_run -> SaveAs((folderName+"/"+folderName+"_history_vsRun.root").c_str()); + + + + o -> cd(); + + h_template -> Write(); + + h_scOccupancy_eta -> Write(); + h_scOccupancy_phi -> Write(); + h_seedOccupancy_EB -> Write(); + h_seedOccupancy_EEp -> Write(); + h_seedOccupancy_EEm -> Write(); + + g_fit -> Write("g_fit"); + g_c_fit -> Write("g_c_fit"); + g_fit_run -> Write("g_fit_run"); + g_c_fit_run -> Write("g_c_fit_run"); + g_las -> Write("g_las"); + g_LT -> Write("g_LT"); + + h_EoP_chi2 -> Write(); + h_EoC_chi2 -> Write(); + + //ciao + + for(int i = 0; i < nBins; ++i) + { + gStyle->SetOptFit(1111); + + h_EoP[i] -> Write(); + h_EoC[i] -> Write(); + f_EoP[i] -> Write(); + f_EoC[i] -> Write(); + // h_Tsp[i] -> Write(); + // + // h_Cvl[i] -> Write(); + } + + o -> Close(); + + return 0; + } + //------------------------------ LASER MONITORING WITH E/P versus PILE UP ------------------------------------------------------ + + if(vm.count("laserMonitoringEPvsPU")) { + + + //float yMIN = 0.95; + //float yMAX = 1.05; + + + // Set style options + setTDRStyle(); + gStyle->SetPadTickX(1); + gStyle->SetPadTickY(1); + gStyle->SetOptTitle(0); + gStyle->SetOptStat(1110); + gStyle->SetOptFit(1); + + // Set fitting options + TVirtualFitter::SetDefaultFitter("Fumili2"); + + + + //----------------- + // Input parameters + + + std::cout << "\n***************************************************************************************************************************" << std::endl; + + std::string dayMinLabel = ""; + std::string dayMaxLabel = ""; + std::string dayZOOMLabel =""; + float absEtaMin = -1.; + float absEtaMax = -1.; + int IetaMin = -1; + int IetaMax = -1; + int IphiMin = -1; + int IphiMax = -1; + + + + + int t1 = 0; + int t2 = 60; + + std::cout << "EBEE: " << EBEE << std::endl; + std::cout << "evtsPerPoint: " << evtsPerPoint << std::endl; + std::cout << "useRegression: " << useRegression << std::endl; + std::cout << "dayMin: " << dayMin << std::endl; + std::cout << "dayZOOM: " << dayZOOM << std::endl; + std::cout << "dayMax: " << dayMax << std::endl; + std::cout << "yMin: " << yMIN << std::endl; + std::cout << "yMax: " << yMAX << std::endl; + std::cout << "absEtaMin: " << absEtaMin << std::endl; + std::cout << "absEtaMax: " << absEtaMax << std::endl; + std::cout << "IetaMin: " << IetaMin << std::endl; + std::cout << "IetaMax: " << IetaMax << std::endl; + std::cout << "IphiMin: " << IphiMin << std::endl; + std::cout << "IphiMax: " << IphiMax << std::endl; + std::cout << "t1: " << t1 << std::endl; + std::cout << "t2: " << t2 << std::endl; + + + std::string dayZOOM = ""; + std::string dayMin = ""; + std::string dayMax = ""; + + + //------------------- + // Define in/outfiles + + std::string folderName = std::string(EBEE) + "_" + dayMinLabel + "_" + dayMaxLabel; + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + char absEtaBuffer[50]; + sprintf(absEtaBuffer,"_%.2f-%.2f",absEtaMin,absEtaMax); + folderName += std::string(absEtaBuffer); + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + char absEtaBuffer[50]; + sprintf(absEtaBuffer,"_Ieta_%d-%d_Iphi_%d-%d",IetaMin,IetaMax,IphiMin,IphiMax); + folderName += std::string(absEtaBuffer); + } + + gSystem->mkdir(folderName.c_str()); + TFile* o = new TFile((folderName+"/"+folderName+"_histos.root").c_str(),"RECREATE"); + + + + // Get trees + std::cout << std::endl; + + + if (data->GetEntries() == 0 || mc->GetEntries() == 0 ) + { + std::cout << "Error: At least one file is empty" << std::endl; + return -1; + } + + + + // Set branch addresses + int runNumber; + int nPV; + int nPU; + float avgLCSCEle[3], seedLaserAlphaSCEle1, etaSCEle[3], phiSCEle[3], energySCEle[3], esEnergySCEle[3], pAtVtxGsfEle[3], energySCEle_corr[3]; + int seedXSCEle[3], seedYSCEle[3];//, seedZside; + + data->SetBranchStatus("*",0); + data->SetBranchStatus("runNumber",1); + data->SetBranchStatus("nPV",1); + data->SetBranchStatus("nPU",1); + data->SetBranchStatus("avgLCSCEle",1); + data->SetBranchStatus("seedLaserAlphaSCEle1",1); + // data->SetBranchStatus("ele1_EOverP",1); + data->SetBranchStatus("etaSCEle",1); + data->SetBranchStatus("phiSCEle",1); + data->SetBranchStatus("energySCEle",1); + data->SetBranchStatus("energySCEle_corr",1); + data->SetBranchStatus("esEnergySCEle",1); + data->SetBranchStatus("pAtVtxGsfEle",1); + data->SetBranchStatus("seedXSCEle",1); + data->SetBranchStatus("seedYSCEle",1); + // data->SetBranchStatus("ele1_seedZside",1); + + data->SetBranchAddress("runNumber", &runNumber); + data->SetBranchAddress("nPV", &nPV); + data->SetBranchAddress("nPU", &nPU); + data->SetBranchAddress("avgLCSCEle", &avgLCSCEle[0]); + data->SetBranchAddress("seedLaserAlphaSCEle1", &seedLaserAlphaSCEle1); + // data->SetBranchAddress("ele1_EOverP", &EoP); + data->SetBranchAddress("etaSCEle", &etaSCEle); + data->SetBranchAddress("phiSCEle", &phiSCEle); + if( useRegression < 1 ) + data->SetBranchAddress("energySCEle", &energySCEle); + else + data->SetBranchAddress("energySCEle_corr", &energySCEle_corr); + data->SetBranchAddress("energySCEle", &energySCEle); + data->SetBranchAddress("esEnergySCEle", &esEnergySCEle); + data->SetBranchAddress("pAtVtxGsfEle", &pAtVtxGsfEle); + data->SetBranchAddress("seedXSCEle", &seedXSCEle); + data->SetBranchAddress("seedYSCEle", &seedYSCEle); + // data->SetBranchAddress("ele1_seedZside", &seedZside); + + + mc->SetBranchStatus("*",0); + mc->SetBranchStatus("runNumber",1); + mc->SetBranchStatus("nPV",1); + mc->SetBranchStatus("nPU",1); + mc->SetBranchStatus("avgLCSCEle",1); + mc->SetBranchStatus("seedLaserAlphaSCEle1",1); + // mc->SetBranchStatus("ele1_EOverP",1); + mc->SetBranchStatus("etaSCEle",1); + mc->SetBranchStatus("phiSCEle",1); + mc->SetBranchStatus("energySCEle",1); + mc->SetBranchStatus("energySCEle_corr",1); + mc->SetBranchStatus("esEnergySCEle",1); + mc->SetBranchStatus("pAtVtxGsfEle",1); + mc->SetBranchStatus("seedXSCEle",1); + mc->SetBranchStatus("seedYSCEle",1); + // mc->SetBranchStatus("ele1_seedZside",1); + + mc->SetBranchAddress("runNumber", &runNumber); + mc->SetBranchAddress("nPV", &nPV); + mc->SetBranchAddress("nPU", &nPU); + mc->SetBranchAddress("avgLCSCEle", &avgLCSCEle[0]); + mc->SetBranchAddress("seedLaserAlphaSCEle1", &seedLaserAlphaSCEle1); + // mc->SetBranchAddress("ele1_EOverP", &EoP); + mc->SetBranchAddress("etaSCEle", &etaSCEle); + mc->SetBranchAddress("phiSCEle", &phiSCEle); + if( useRegression < 1 ) + mc->SetBranchAddress("energySCEle", &energySCEle); + else + mc->SetBranchAddress("energySCEle_corr", &energySCEle_corr); + mc->SetBranchAddress("energySCEle", &energySCEle); + mc->SetBranchAddress("esEnergySCEle", &esEnergySCEle); + mc->SetBranchAddress("pAtVtxGsfEle", &pAtVtxGsfEle); + mc->SetBranchAddress("seedXSCEle", &seedXSCEle); + mc->SetBranchAddress("seedYSCEle", &seedYSCEle); + // mc->SetBranchAddress("ele1_seedZside", &seedZside); + + + + + + + //-------------------------------------------------------- + // Define PU correction (to be used if useRegression == 0) + + // corr = p0 + p1 * nPU + float p0_EB; + float p1_EB; + float p0_EE; + float p1_EE; + + if( useRegression == 0 ) + { + //2012 EB + p0_EB = 0.9991; + p1_EB = 0.0001635; + //2012 EE + p0_EE = 0.9968; + p1_EE = 0.001046; + } + else + { + //2012 EB + p0_EB = 1.001; + p1_EB = -0.000143; + //2012 EE + p0_EE = 1.00327; + p1_EE = -0.000432; + } + + float p0 = -1.; + float p1 = -1.; + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + { + p0 = p0_EB; + p1 = p1_EB; + } + else + { + p0 = p0_EE; + p1 = p1_EE; + } + + //2015 + p0=1.; + p1=0.; + + + + + //--------------------------------- + // Build the reference distribution + + std::cout << std::endl; + std::cout << "***** Build reference for " << EBEE << " *****" << std::endl; + + TH1F* h_template = new TH1F("template", "", 2000, 0., 5.); + + for(int ientry = 0; ientry < mc->GetEntries(); ++ientry) + { + if( (ientry%100000 == 0) ) std::cout << "reading MC entry " << ientry << "\r" << std::endl;//std::flush; + mc->GetEntry(ientry); + + // selections + if( (strcmp(EBEE.c_str(),"EB") == 0) && (fabs(etaSCEle[0]) > 1.479) ) continue; // barrel + if( (strcmp(EBEE.c_str(),"EE") == 0) && (fabs(etaSCEle[0]) < 1.479 || fabs(etaSCEle[0])>2.5) ) continue; // endcap + + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + if( (fabs(etaSCEle[0]) < absEtaMin) || (fabs(etaSCEle[0]) > absEtaMax) ) continue; + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + if( (seedXSCEle[0] < IetaMin) || (seedXSCEle[0] > IetaMax) ) continue; + if( (seedYSCEle[0] < IphiMin) || (seedYSCEle[0] > IphiMax) ) continue; + } + + // PU correction + float PUCorr = (p0 + p1*nPU); + //std::cout << "p0: " << p0 << " p1: " << p1 << " nPU: " << nPU << std::endl; + + // fill the template histogram + h_template -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / PUCorr ); + } + + std::cout << "Reference built for " << EBEE << " - " << h_template->GetEntries() << " events" << std::endl; + + + + + + + //--------------------- + // Loop and sort events + + std::cout << std::endl; + std::cout << "***** Sort events and define bins *****" << std::endl; + + int nEntries = data -> GetEntriesFast(); + int nSavePts = 0; + std::vector isSavedEntries(nEntries); + std::vector sortedEntries; + std::vector timeStampFirst; + + for(int ientry = 0; ientry < nEntries; ++ientry) + { + data -> GetEntry(ientry); + isSavedEntries.at(ientry) = false; + + // selections + if( (strcmp(EBEE.c_str(),"EB") == 0) && (fabs(etaSCEle[0]) > 1.479) ) continue; // barrel + if( (strcmp(EBEE.c_str(),"EE") == 0) && (fabs(etaSCEle[0]) < 1.479 || fabs(etaSCEle[0])>2.5) ) continue; // endcap + + if( (absEtaMin != -1.) && (absEtaMax != -1.) ) + { + if( (fabs(etaSCEle[0]) < absEtaMin) || (fabs(etaSCEle[0]) > absEtaMax) ) continue; + } + + if( (IetaMin != -1.) && (IetaMax != -1.) && (IphiMin != -1.) && (IphiMax != -1.) ) + { + if( (seedXSCEle[0] < IetaMin) || (seedXSCEle[0] > IetaMax) ) continue; + if( (seedYSCEle[0] < IphiMin) || (seedYSCEle[0] > IphiMax) ) continue; + } + + if( nPV < t1 ) continue; + if( nPV > t2 ) continue; + + if( avgLCSCEle[0] <= 0. ) continue; + + isSavedEntries.at(ientry) = true; + + + // fill sorter + Sorter dummy; + dummy.time = nPV; + dummy.entry = ientry; + sortedEntries.push_back(dummy); + + ++nSavePts; + } + + // sort events + std::sort(sortedEntries.begin(),sortedEntries.end(),Sorter()); + std::cout << "Data sorted in " << EBEE << " - " << nSavePts << " events" << std::endl; + + std::map antiMap; + for(unsigned int iSaved = 0; iSaved < sortedEntries.size(); ++iSaved) + antiMap[sortedEntries.at(iSaved).entry] = iSaved; + + + //--------------------- + // Loop and define bins + + // "wide" bins - find events with time separation bigger than 1 day + int nWideBins = 1; + std::vector wideBinEntryMax; + //int timeStampOld = -1; + + wideBinEntryMax.push_back(0); + + for(int iSaved = 0; iSaved < nSavePts; ++iSaved) + { + /* if( iSaved%100000 == 0 ) std::cout << "reading saved entry " << iSaved << "\r" << std::flush; + data->GetEntry(sortedEntries[iSaved].entry); + + if( iSaved == 0 ) + { + timeStampOld = nPV; + continue; + } + + if( (nPV-timeStampOld)/3600. > timeLapse ) + { + ++nWideBins; + wideBinEntryMax.push_back(iSaved-1); + } + + + timeStampOld = nPV; + */ + } + + + std::cout << std::endl; + wideBinEntryMax.push_back(nSavePts); + + // bins with approximatively evtsPerPoint events per bin + int nBins = 0; + std::vector binEntryMax; + + binEntryMax.push_back(0); + for(int wideBin = 0; wideBin < nWideBins; ++wideBin) + { + int nTempBins = std::max(1,int( (wideBinEntryMax.at(wideBin+1)-wideBinEntryMax.at(wideBin))/evtsPerPoint )); + int nTempBinEntries = int( (wideBinEntryMax.at(wideBin+1)-wideBinEntryMax.at(wideBin))/nTempBins ); + + for(int tempBin = 0; tempBin < nTempBins; ++tempBin) + { + ++nBins; + if( tempBin < nTempBins - 1 ) + binEntryMax.push_back( wideBinEntryMax.at(wideBin) + (tempBin+1)*nTempBinEntries ); + else + binEntryMax.push_back( wideBinEntryMax.at(wideBin+1) ); + } + } + + // std::cout << "nBins = " << nBins << std::endl; + //for(int bin = 0; bin < nBins; ++bin) + // std::cout << "bin: " << bin + // << " entry min: " << setw(6) << binEntryMax.at(bin) + // << " entry max: " << setw(6) << binEntryMax.at(bin+1) + // << " events: " << setw(6) << binEntryMax.at(bin+1)-binEntryMax.at(bin) + // << std::endl; + + + + + + + //--------------------- + // histogram definition + + TH1F* h_scOccupancy_eta = new TH1F("h_scOccupancy_eta","", 298, -2.6, 2.6); + TH1F* h_scOccupancy_phi = new TH1F("h_scOccupancy_phi","", 363, -3.1765, 3.159); + SetHistoStyle(h_scOccupancy_eta); + SetHistoStyle(h_scOccupancy_phi); + + TH2F* h_seedOccupancy_EB = new TH2F("h_seedOccupancy_EB","", 171, -86., 85., 361, 0.,361.); + TH2F* h_seedOccupancy_EEp = new TH2F("h_seedOccupancy_EEp","", 101, 0.,101., 100, 0.,101.); + TH2F* h_seedOccupancy_EEm = new TH2F("h_seedOccupancy_EEm","", 101, 0.,101., 100, 0.,101.); + SetHistoStyle(h_seedOccupancy_EB); + SetHistoStyle(h_seedOccupancy_EEp); + SetHistoStyle(h_seedOccupancy_EEm); + + TH1F* h_EoP_spread = new TH1F("h_EoP_spread","",100,yMIN,yMAX); + TH1F* h_EoC_spread = new TH1F("h_EoC_spread","",100,yMIN,yMAX); + TH1F* h_EoP_spread_run = new TH1F("h_EoP_spread_run","",100,yMIN,yMAX); + TH1F* h_EoC_spread_run = new TH1F("h_EoC_spread_run","",100,yMIN,yMAX); + SetHistoStyle(h_EoP_spread,"EoP"); + SetHistoStyle(h_EoC_spread,"EoC"); + SetHistoStyle(h_EoP_spread_run,"EoP"); + SetHistoStyle(h_EoC_spread_run,"EoC"); + + TH1F* h_EoP_chi2 = new TH1F("h_EoP_chi2","",50,0.,5.); + TH1F* h_EoC_chi2 = new TH1F("h_EoC_chi2","",50,0.,5.); + SetHistoStyle(h_EoP_chi2,"EoP"); + SetHistoStyle(h_EoC_chi2,"EoC"); + + TH1F** h_EoP = new TH1F*[nBins]; + TH1F** h_EoC = new TH1F*[nBins]; + TH1F** h_Las = new TH1F*[nBins]; + TH1F** h_Tsp = new TH1F*[nBins]; + TH1F** h_Cvl = new TH1F*[nBins]; + + for(int i = 0; i < nBins; ++i) + { + char histoName[80]; + + sprintf(histoName, "EoP_%d", i); + h_EoP[i] = new TH1F(histoName, histoName, 2000, 0., 5.); + SetHistoStyle(h_EoP[i],"EoP"); + + sprintf(histoName, "EoC_%d", i); + h_EoC[i] = new TH1F(histoName, histoName, 2000, 0., 5.); + SetHistoStyle(h_EoC[i],"EoC"); + + sprintf(histoName, "Las_%d", i); + h_Las[i] = new TH1F(histoName, histoName, 500, 0.5, 1.5); + + sprintf(histoName, "Tsp_%d", i); + h_Tsp[i] = new TH1F(histoName, histoName, 500, 0.5, 1.5); + + } + + + // function definition + TF1** f_EoP = new TF1*[nBins]; + TF1** f_EoC = new TF1*[nBins]; + + + // graphs definition + TGraphAsymmErrors* g_fit = new TGraphAsymmErrors(); + TGraphAsymmErrors* g_c_fit = new TGraphAsymmErrors(); + + TGraphAsymmErrors* g_fit_run = new TGraphAsymmErrors(); + TGraphAsymmErrors* g_c_fit_run = new TGraphAsymmErrors(); + + TGraph* g_las = new TGraph(); + + TGraphErrors* g_LT = new TGraphErrors(); + + g_fit->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_fit->GetXaxis()->SetTimeDisplay(1); + g_c_fit->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_c_fit->GetXaxis()->SetTimeDisplay(1); + g_las->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_las->GetXaxis()->SetTimeDisplay(1); + g_LT->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + g_LT->GetXaxis()->SetTimeDisplay(1); + + + + + + + //------------------------------------ + // loop on the saved and sorted events + + std::cout << std::endl; + std::cout << "***** Fill and fit histograms *****" << std::endl; + + std::vector Entries(nBins); + std::vector AveTime(nBins); + std::vector MinTime(nBins); + std::vector MaxTime(nBins); + std::vector AveRun(nBins); + std::vector MinRun(nBins); + std::vector MaxRun(nBins); + std::vector AveLT(nBins); + std::vector AveLT2(nBins); + + int iSaved = -1; + for(int ientry = 0; ientry < nEntries; ++ientry) + { + if( (ientry%100000 == 0) ) std::cout << "reading entry " << ientry << "\r" << std::endl;//std::flush; + + if( isSavedEntries.at(ientry) == false ) continue; + + ++iSaved; + + int iSaved = antiMap[ientry]; + int bin = -1; + for(bin = 0; bin < nBins; ++bin) + if( iSaved >= binEntryMax.at(bin) && iSaved < binEntryMax.at(bin+1) ) + break; + + //std::cout << "bin = " << bin << " iSaved = "<< iSaved << std::endl; + data->GetEntry(ientry); + + + + Entries[bin] += 1; + + if( iSaved == binEntryMax.at(bin)+1 ) MinTime[bin] = nPV; + if( iSaved == binEntryMax.at(bin+1)-1 ) MaxTime[bin] = nPV; + AveTime[bin] += nPV; + + if( iSaved == binEntryMax.at(bin)+1 ) MinRun[bin] = runNumber; + if( iSaved == binEntryMax.at(bin+1)-1 ) MaxRun[bin] = runNumber; + AveRun[bin] += runNumber; + + float LT = (-1. / seedLaserAlphaSCEle1 * log(avgLCSCEle[0])); + AveLT[bin] += LT; + AveLT2[bin] += LT*LT; + + // PU correction + float PUCorr = (p0 + p1*nPU); + + // fill the histograms + (h_EoP[bin]) -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / avgLCSCEle[0] / PUCorr); + (h_EoC[bin]) -> Fill( (energySCEle[0]-esEnergySCEle[0])/(pAtVtxGsfEle[0]-esEnergySCEle[0]) / PUCorr ); + + (h_Las[bin]) -> Fill(avgLCSCEle[0]); + (h_Tsp[bin]) -> Fill(1./avgLCSCEle[0]); + + h_scOccupancy_eta -> Fill(etaSCEle[0]); + h_scOccupancy_phi -> Fill(phiSCEle[0]); + if(fabs(etaSCEle[0])<1.449) + h_seedOccupancy_EB -> Fill(seedXSCEle[0],seedYSCEle[0]); + else if(etaSCEle[0]>1.449) + h_seedOccupancy_EEp -> Fill(seedXSCEle[0],seedYSCEle[0]); + else if(etaSCEle[0]<-1.449) + h_seedOccupancy_EEm -> Fill(seedXSCEle[0],seedYSCEle[0]); + } + + for(int bin = 0; bin < nBins; ++bin) + { + AveTime[bin] = 1. * AveTime[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveRun[bin] = 1. * AveRun[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveLT[bin] = 1. * AveLT[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + AveLT2[bin] = 1. * AveLT2[bin] / (binEntryMax.at(bin+1)-binEntryMax.at(bin)); + //std::cout << date << " " << AveTime[i] << " " << MinTime[i] << " " << MaxTime[i] << std::endl; + } + + + + + + + int rebin = 2; + if( strcmp(EBEE.c_str(),"EE") == 0 ) rebin *= 2; + + h_template -> Rebin(rebin); + + + + float EoP_scale = 0.; + float EoP_err = 0.; + int EoP_nActiveBins = 0; + + float EoC_scale = 0.; + float EoC_err = 0.; + int EoC_nActiveBins = 0; + + float LCInv_scale = 0; + + std::vector validBins; + for(int i = 0; i < nBins; ++i) + { + bool isValid = true; + + h_EoP[i] -> Rebin(rebin); + h_EoC[i] -> Rebin(rebin); + + + + //------------------------------------ + // Fill the graph for uncorrected data + + // define the fitting function + // N.B. [0] * ( [1] * f( [1]*(x-[2]) ) ) + + //o -> cd(); + char convolutionName[50]; + sprintf(convolutionName,"h_convolution_%d",i); + //h_Cvl[i] = ConvoluteTemplate(std::string(convolutionName),h_template,h_Las[i],32768,-5.,5.); + h_Cvl[i] = MellinConvolution(std::string(convolutionName),h_template,h_Tsp[i]); + + histoFunc* templateHistoFunc = new histoFunc(h_template); + histoFunc* templateConvolutedHistoFunc = new histoFunc(h_Cvl[i]); + char funcName[50]; + + sprintf(funcName,"f_EoP_%d",i); + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + f_EoP[i] = new TF1(funcName, templateConvolutedHistoFunc, 0.8*(h_Tsp[i]->GetMean()), 1.4*(h_Tsp[i]->GetMean()), 3, "histoFunc"); + else + f_EoP[i] = new TF1(funcName, templateConvolutedHistoFunc, 0.75*(h_Tsp[i]->GetMean()), 1.5*(h_Tsp[i]->GetMean()), 3, "histoFunc"); + + f_EoP[i] -> SetParName(0,"Norm"); + f_EoP[i] -> SetParName(1,"Scale factor"); + f_EoP[i] -> SetLineWidth(1); + f_EoP[i] -> SetNpx(10000); + + double xNorm = h_EoP[i]->GetEntries()/h_template->GetEntries() * + h_EoP[i]->GetBinWidth(1)/h_template->GetBinWidth(1); + + f_EoP[i] -> FixParameter(0, xNorm); + f_EoP[i] -> SetParameter(1, 1.); + f_EoP[i] -> FixParameter(2, 0.); + f_EoP[i] -> SetLineColor(kRed+2); + + int fStatus = 0; + int nTrials = 0; + TFitResultPtr rp; + + rp = h_EoP[i] -> Fit(funcName, "ERLS+"); + while( (fStatus != 0) && (nTrials < 10) ) + { + rp = h_EoP[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + if(fStatus == 0) break; + ++nTrials; + } + + // fill the graph + double eee = f_EoP[i]->GetParError(1); + //float k = f_EoP[i]->GetParameter(1); + float k = f_EoP[i]->GetParameter(1) / h_Tsp[i]->GetMean(); //needed when using mellin's convolution + // std::cout << "noCORR" << std::endl; + // std::cout << "eee: " << f_EoP[i]->GetParError(1) << "k: " << f_EoP[i]->GetParameter(1) / h_Tsp[i]->GetMean() << std::endl; + + /* + std::cout << i <<"--nocorr---- "<< 1./k << std::endl; + std::cout <<" condizione 1: " << h_EoP[i]->GetEntries() << " fStatus: " << fStatus << " eee: " << eee << "con eee che ci piace essere maggiore di : " << 0.05*h_template->GetRMS()/sqrt(evtsPerPoint) << std::endl ; + getchar(); + */ + // if( (h_EoP[i]->GetEntries() > 3) && (fStatus == 0) && (eee > 0.05*h_template->GetRMS()/sqrt(evtsPerPoint)) ) + + if( (h_EoP[i]->GetEntries() > 500) && (fStatus == 0) ) + { + float date = (float)AveTime[i]; + float dLow = (float)(AveTime[i]-MinTime[i]); + float dHig = (float)(MaxTime[i]-AveTime[i]); + float run = (float)AveRun[i]; + float rLow = (float)(AveRun[i]-MinRun[i]); + float rHig = (float)(MaxRun[i]-AveRun[i]); + g_fit -> SetPoint(i, date , 1./k); + g_fit -> SetPointError(i, dLow , dHig, eee/k/k, eee/k/k); + g_fit_run -> SetPoint(i, run , 1./k); + g_fit_run -> SetPointError(i, rLow , rHig, eee/k/k, eee/k/k); + + std::cout <<"************-------------------*****************" < Fill(f_EoP[i]->GetChisquare()/f_EoP[i]->GetNDF()); + + EoP_scale += 1./k; + EoP_err += eee/k/k; + ++EoP_nActiveBins; + } + else + { + std::cout << "Fitting uncorrected time bin: " << i << " Fail status: " << fStatus << " sigma: " << eee << std::endl; + isValid = false; + } + + //---------------------------------- + // Fill the graph for corrected data + + // define the fitting function + // N.B. [0] * ( [1] * f( [1]*(x-[2]) ) ) + + sprintf(funcName,"f_EoC_%d",i); + if( strcmp(EBEE.c_str(),"EB") == 0 ) + f_EoC[i] = new TF1(funcName, templateHistoFunc, 0.8, 1.4, 3, "histoFunc"); + else + f_EoC[i] = new TF1(funcName, templateHistoFunc, 0.75, 1.5, 3, "histoFunc"); + f_EoC[i] -> SetParName(0,"Norm"); + f_EoC[i] -> SetParName(1,"Scale factor"); + f_EoC[i] -> SetLineWidth(1); + f_EoC[i] -> SetNpx(10000); + + xNorm = h_EoC[i]->GetEntries()/h_template->GetEntries() * + h_EoC[i]->GetBinWidth(1)/h_template->GetBinWidth(1); + + f_EoC[i] -> FixParameter(0, xNorm); + f_EoC[i] -> SetParameter(1, 0.99); + f_EoC[i] -> FixParameter(2, 0.); + f_EoC[i] -> SetLineColor(kGreen+2); + + + rp = h_EoC[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + nTrials = 0; + while( (fStatus != 0) && (nTrials < 10) ) + { + rp = h_EoC[i] -> Fit(funcName, "ERLS+"); + fStatus = rp; + if(fStatus == 0) break; + ++nTrials; + } + + + // fill the graph + k = f_EoC[i]->GetParameter(1); + eee = f_EoC[i]->GetParError(1); + //std::cout << "CORR" << std::endl; + //std::cout << "eee: " << f_EoP[i]->GetParError(1) << "k: " << f_EoP[i]->GetParameter(1) / h_Tsp[i]->GetMean() << std::endl; + //getchar(); + /* std::cout << i <<"--corr---- "<< 1./k << std::endl; + std::cout <<" condizione 1: " << h_EoP[i]->GetEntries() << " fStatus: " << fStatus << " eee: " << eee << "con eee che ci piace essere maggiore di : " << 0.05*h_template->GetRMS()/sqrt(evtsPerPoint) << std::endl ; + getchar(); + */ + + + if( (h_EoC[i]->GetEntries() > 500) && (fStatus == 0) ) + // if( (h_EoC[i]->GetEntries() > 10) && (fStatus == 0) && (eee > 0.05*h_template->GetRMS()/sqrt(evtsPerPoint)) ) + { + float date = (float)AveTime[i]; + float dLow = (float)(AveTime[i]-MinTime[i]); + float dHig = (float)(MaxTime[i]-AveTime[i]); + float run = (float)AveRun[i]; + float rLow = (float)(AveRun[i]-MinRun[i]); + float rHig = (float)(MaxRun[i]-AveRun[i]); + + g_c_fit -> SetPoint(i, date , 1./k); + g_c_fit -> SetPointError(i, dLow , dHig , eee/k/k, eee/k/k); + + g_c_fit_run -> SetPoint(i, run , 1./k); + g_c_fit_run -> SetPointError(i, rLow , rHig, eee/k/k, eee/k/k); + std::cout <<"************-------------------*****************" < Fill(f_EoC[i]->GetChisquare()/f_EoP[i]->GetNDF()); + + EoC_scale += 1./k; + EoC_err += eee/k/k; + ++EoC_nActiveBins; + } + else + { + std::cout << "Fitting corrected time bin: " << i << " Fail status: " << fStatus << " sigma: " << eee << std::endl; + isValid = false; + } + + if( isValid == true ) validBins.push_back(i); + } + + EoP_scale /= EoP_nActiveBins; + EoP_err /= EoP_nActiveBins; + + EoC_scale /= EoC_nActiveBins; + EoC_err /= EoC_nActiveBins; + + + + + + + //---------------------------------------- + // Fill the graph for avg laser correction + + + for(unsigned int itr = 0; itr < validBins.size(); ++itr) + { + //float k0 = f_EoP[0]->GetParameter(1) / h_Tsp[0]->GetMean(); + int i = validBins.at(itr); + // g_las -> SetPoint(itr, (float)AveTime[i], (h_Tsp[i]->GetMean())+((1/k0)-(h_Tsp[0]->GetMean())) ); + g_las -> SetPoint(itr, (float)AveTime[i], h_Tsp[i]->GetMean() ); + + + //g_las -> SetPointffa(itr, (float)AveTime[i], h_Tsp[i]->GetMean()); + g_LT -> SetPoint(itr, (float)AveTime[i], AveLT[i] ); + g_LT -> SetPointError(itr, 0., sqrt(AveLT2[i]-AveLT[i]*AveLT[i]) / sqrt(Entries[i]) ); + + LCInv_scale += h_Tsp[i]->GetMean(); + } + + LCInv_scale /= validBins.size(); + + + + + + + //--------------- + // Rescale graphs + + float yscale = 1.; + //float yscale = 1./EoC_scale; + + for(unsigned int itr = 0; itr < validBins.size(); ++itr) + { + double x,y; + g_fit -> GetPoint(itr,x,y); + g_fit -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoP_spread -> Fill(y*yscale); + + g_c_fit -> GetPoint(itr,x,y); + g_c_fit -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoC_spread -> Fill(y*yscale); + + g_fit_run -> GetPoint(itr,x,y); + g_fit_run -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoP_spread_run -> Fill(y*yscale); + + g_c_fit_run -> GetPoint(itr,x,y); + g_c_fit_run -> SetPoint(itr,x,y*yscale); + if ( (x > t1) && (x < t2) ) h_EoC_spread_run -> Fill(y*yscale); + + g_las -> GetPoint(itr,x,y); + g_las -> SetPoint(itr,x,y*yscale*EoP_scale/LCInv_scale); + } + //ciao + TF1 EoC_pol0("EoC_pol0","pol0",t1,t2); + EoC_pol0.SetLineColor(kGreen+2); + EoC_pol0.SetLineWidth(2); + EoC_pol0.SetLineStyle(2); + g_c_fit -> Fit("EoC_pol0","QNR"); + + + + + + + + + //---------------------------- + // Print out global quantities + + std::cout << std::endl; + std::cout << "***** Mean scales and errors *****" << std::endl; + std::cout << std::fixed; + std::cout << std::setprecision(4); + std::cout << "Mean EoP scale: " << std::setw(6) << EoP_scale << " mean EoP error: " << std::setw(8) << EoP_err << std::endl; + std::cout << "Mean EoC scale: " << std::setw(6) << EoC_scale << " mean EoC error: " << std::setw(8) << EoC_err << std::endl; + std::cout << "Mean 1/LC scale: " << std::setw(6) << LCInv_scale << std::endl; + + + + + + + //------------------- + // Final Plot vs Vertex + //------------------- + + TCanvas* cplot = new TCanvas("cplot", "history plot vs Vertex",100,100,1000,500); + cplot->cd(); + + TPad *cLeft = new TPad("pad_0","pad_0",0.00,0.00,0.75,1.00); + TPad *cRight = new TPad("pad_1","pad_1",0.75,0.00,1.00,1.00); + + cLeft->SetLeftMargin(0.15); + cLeft->SetRightMargin(0.025); + cRight->SetLeftMargin(0.025); + + cLeft->Draw(); + cRight->Draw(); + + float tYoffset = 1.0; + float labSize = 0.05; + float labSize2 = 0.06; + + cLeft->cd(); + + cLeft->SetGridx(); + cLeft->SetGridy(); + + TH1F *hPad = (TH1F*)gPad->DrawFrame(t1,0.9,t2,1.05); + + hPad->GetXaxis()->SetLimits(0,46); + //hPad->GetXaxis()->SetTimeFormat("%d/%m%F1970-01-01 00:00:00"); + //hPad->GetXaxis()->SetTimeDisplay(1); + //hPad->GetXaxis() -> SetRangeUser(MinTime[0]-43200,MaxTime[nBins-1]+43200); + hPad->GetXaxis()->SetTitle(" Number of Vertices"); + hPad->GetXaxis()->SetTitleOffset(0.8); + + //ciao + + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + else + hPad->GetYaxis()->SetTitle("Relative E/p scale"); + hPad->GetYaxis()->SetTitleOffset(tYoffset); + hPad->GetXaxis()->SetLabelSize(labSize); + hPad->GetXaxis()->SetTitleSize(labSize2); + hPad->GetYaxis()->SetLabelSize(labSize); + hPad->GetYaxis()->SetTitleSize(labSize2); + hPad -> SetMinimum(yMIN); + hPad -> SetMaximum(yMAX); + + // draw history plot + g_fit -> SetMarkerStyle(24); + g_fit -> SetMarkerSize(0.7); + g_fit -> SetMarkerColor(kRed+2); + g_fit -> SetLineColor(kRed+2); + //g_fit -> Draw("P"); + g_c_fit -> SetMarkerStyle(20); + g_c_fit -> SetMarkerColor(kGreen+2); + g_c_fit -> SetLineColor(kGreen+2); + g_c_fit -> SetMarkerSize(0.7); + g_c_fit -> Draw("EP"); + g_las -> SetLineColor(kAzure-2); + g_las -> SetLineWidth(2); + //g_las -> Draw("L,same"); + + TLegend* legend = new TLegend(0.60,0.78,0.90,0.94); + legend -> SetLineColor(kWhite); + legend -> SetLineWidth(0); + legend -> SetFillColor(kWhite); + legend -> SetFillStyle(0); + legend -> SetTextFont(42); + legend -> SetTextSize(0.04); + legend -> AddEntry(g_c_fit,"with LM correction","PL"); + //legend -> AddEntry(g_fit, "without LM correction","PL"); + //legend -> AddEntry(g_las, "1 / LM correction","L"); + legend -> Draw("same"); + + char latexBuffer[250]; + + sprintf(latexBuffer,"CMS 2015 Preliminary"); + TLatex* latex = new TLatex(0.18,0.89,latexBuffer); + latex -> SetNDC(); + latex -> SetTextFont(62); + latex -> SetTextSize(0.05); + latex -> Draw("same"); + + //sprintf(latexBuffer,"#sqrt{s} = 8 TeV L = 3.95 fb^{-1}"); + sprintf(latexBuffer,"#sqrt{s} = 13 TeV"); + TLatex* latex2 = new TLatex(0.18,0.84,latexBuffer); + latex2 -> SetNDC(); + latex2 -> SetTextFont(42); + latex2 -> SetTextSize(0.05); + latex2 -> Draw("same"); + + if( strcmp(EBEE.c_str(),"EB") == 0 ) + sprintf(latexBuffer,"ECAL Barrel"); + else + sprintf(latexBuffer,"ECAL Endcap"); + TLatex* latex3 = new TLatex(0.18,0.19,latexBuffer); + latex3 -> SetNDC(); + latex3 -> SetTextFont(42); + latex3 -> SetTextSize(0.05); + latex3 -> Draw("same"); + + //sprintf(latexBuffer,"%.2E events",1.*nSavePts); + //TLatex* latex4 = new TLatex(0.18,0.24,latexBuffer); + //latex4 -> SetNDC(); + //latex4 -> SetTextFont(42); + //latex4 -> SetTextSize(0.04); + //latex4 -> Draw("same"); + // + //sprintf(latexBuffer,"%d events/bin - %d bins",evtsPerPoint,nBins); + //TLatex* latex5 = new TLatex(0.18,0.19,latexBuffer); + //latex5 -> SetNDC(); + //latex5 -> SetTextFont(42); + //latex5 -> SetTextSize(0.04); + //latex5 -> Draw("same"); + + + cRight -> cd(); + + TPaveStats* s_EoP_spread = new TPaveStats(); + TPaveStats* s_EoC_spread = new TPaveStats(); + + + h_EoC_spread -> SetFillStyle(3001); + h_EoC_spread -> SetFillColor(kGreen+2); + h_EoC_spread->GetYaxis()->SetLabelSize(0.09); + h_EoC_spread->GetYaxis()->SetLabelOffset(-0.03); + h_EoC_spread->GetYaxis()->SetTitleSize(0.08); + h_EoC_spread->GetYaxis()->SetNdivisions(505); + h_EoC_spread->GetXaxis()->SetLabelOffset(1000); + + h_EoC_spread -> Draw("hbar"); + gPad -> Update(); + + s_EoC_spread = (TPaveStats*)(h_EoC_spread->GetListOfFunctions()->FindObject("stats")); + s_EoC_spread -> SetStatFormat("1.4g"); + s_EoC_spread->SetX1NDC(0.06); //new x start position + s_EoC_spread->SetX2NDC(0.71); //new x end position + s_EoC_spread->SetY1NDC(0.43); //new x start position + s_EoC_spread->SetY2NDC(0.34); //new x end position + s_EoC_spread -> SetOptStat(1100); + s_EoC_spread ->SetTextColor(kGreen+2); + s_EoC_spread ->SetTextSize(0.08); + s_EoC_spread -> Draw("sames"); + + + h_EoP_spread -> SetFillStyle(3001); + h_EoP_spread -> SetFillColor(kRed+2); + h_EoP_spread -> Draw("hbarsames"); + gPad -> Update(); + s_EoP_spread = (TPaveStats*)(h_EoP_spread->GetListOfFunctions()->FindObject("stats")); + s_EoP_spread -> SetStatFormat("1.4g"); + s_EoP_spread->SetX1NDC(0.06); //new x start position + s_EoP_spread->SetX2NDC(0.71); //new x end position + s_EoP_spread->SetY1NDC(0.33); //new x start position + s_EoP_spread->SetY2NDC(0.24); //new x end position + s_EoP_spread ->SetOptStat(1100); + s_EoP_spread ->SetTextColor(kRed+2); + s_EoP_spread ->SetTextSize(0.08); + s_EoP_spread -> Draw("sames"); + + /* + h_EoP_spread -> SetFillStyle(3001); + h_EoP_spread -> SetFillColor(kRed+2); + h_EoP_spread -> Draw("hbarsame"); + gPad -> Update(); + */ + + + + + + cplot -> Print((folderName+"/"+folderName+"_history_vsVertex.png").c_str(),"png"); + + cplot -> Print((folderName+"/"+folderName+"_history_vsVertex.pdf").c_str(),"pdf"); + + cplot -> SaveAs((folderName+"/"+folderName+"_history_vsVertex.C").c_str()); + + + + o -> cd(); + + h_template -> Write(); + + h_scOccupancy_eta -> Write(); + h_scOccupancy_phi -> Write(); + h_seedOccupancy_EB -> Write(); + h_seedOccupancy_EEp -> Write(); + h_seedOccupancy_EEm -> Write(); + + g_fit -> Write("g_fit"); + g_c_fit -> Write("g_c_fit"); + g_fit_run -> Write("g_fit_run"); + g_c_fit_run -> Write("g_c_fit_run"); + g_las -> Write("g_las"); + g_LT -> Write("g_LT"); + + h_EoP_chi2 -> Write(); + h_EoC_chi2 -> Write(); + + //ciao + + for(int i = 0; i < nBins; ++i) + { + gStyle->SetOptFit(1111); + + h_EoP[i] -> Write(); + h_EoC[i] -> Write(); + f_EoP[i] -> Write(); + f_EoC[i] -> Write(); + // h_Tsp[i] -> Write(); + // + // h_Cvl[i] -> Write(); + } + + o -> Close(); + + return 0; + } + ///////////--------------------------- E/P calibration ---------------------------------------------------------------------- + + if(vm.count("EOverPCalib") && vm.count("doEB")) { std::cout<<"---- START E/P CALIBRATION: BARREL ----"<1) selection_data = cutter.GetCut(category, false,0,true); selection_data.Print(); - return NULL; + // return NULL; selection_data+=selection; TCut selection_MC=""; if(category.Sizeof()>1) selection_MC = cutter.GetCut(category, false,0); selection_MC+=selection; + std::cout<<"qui"<GetXaxis()->SetBinLabel(index,(*hlt_itr).c_str()); } } else { + + std::cout<<"qui"<Draw(branchNameData+">>data_hist"+binning, selection_data); if(mc!=NULL){ if(usePU) mc->Draw(branchNameMC+">>mc_hist"+binning, selection_MC *"puWeight"); else mc->Draw(branchNameMC+">>mc_hist"+binning, selection_MC); + std::cout<<"qui"<SaveAs("tmp/d_hist.root"); s->SaveAs("tmp/s_hist.root"); + std::cout<<"qui"<GetBinWidth(2), yLabelUnit.Data()); float max = 1.1 * std::max( @@ -570,7 +579,7 @@ TCanvas *PlotDataMCMC(TChain *data, TChain *mc, TChain *mc2, TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchname, TString binning, TString category, TString selection, - TString dataLabel, std::vector mcLabel_vec, TString xLabel, TString yLabelUnit, + TString dataLabel, std::vector mcLabel_vec, TString xLabel, TString yLabelUnit, TString outputPath, TString label4Print, bool logy=false, bool usePU=true, bool ratio=true,bool smear=false, bool scale=false, bool useR9Weight=false, TString pdfIndex=""){ TStopwatch watch; watch.Start(); @@ -586,7 +595,7 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn TPad * pad1 = new TPad("pad1", "pad1",0.01,0.13,0.75,1.); TPad * pad2 = new TPad("pad2", "pad2",0.01,0.001,0.75,0.2); TPad * pad3 = new TPad("pad3", "pad3",0.68,0.001,1.,0.2); - + pad1->SetRightMargin(0.1); pad1->SetLogy(); pad1->Draw(); @@ -606,7 +615,7 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn pad3->cd(); pad1->cd(); - + TString branchNameData=branchname; TString branchNameMC=branchname; @@ -623,6 +632,7 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn // data->SetBranchStatus("invMass_SC_regrCorrSemiParV5_ele", 1); if(branchNameData.Contains("energySCEle_regrCorrSemiParV5_pho")) cutter.energyBranchName="energySCEle_regrCorrSemiParV5_pho"; else if(branchNameData.Contains("energySCEle_regrCorrSemiParV5_ele")) cutter.energyBranchName="energySCEle_regrCorrSemiParV5_ele"; + else if (branchNameData.Contains("energySCEle")) cutter.energyBranchName="energySCEle"; TCut selection_data=""; if(category.Sizeof()>1) selection_data = cutter.GetCut(category, false,0,scale); @@ -688,7 +698,7 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn c->Clear(); - TLegend *leg = new TLegend(0.65,0.8,1,1); + TLegend *leg = new TLegend(0.6,0.75,0.9,0.9); leg->SetBorderSize(1); leg->SetFillColor(0); leg->SetTextSize(0.04); @@ -698,6 +708,8 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn TH1F *d = (TH1F *) gROOT->FindObject("data_hist"); if(dataLabel !="") leg->AddEntry(d,dataLabel,"p"); + d->SetStats(0); + d->SetTitle(""); d->SetMarkerStyle(20); d->SetMarkerSize(1); @@ -730,6 +742,8 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn for(int i=0; i < nHist; i++){ TString mcHistName; mcHistName+=i; mcHistName+="_hist"; TH1F *s = (TH1F *) gROOT->FindObject(mcHistName); + s->SetStats(0); + s->SetTitle(""); if(s==NULL) continue; std::cout << "nEvents signal: " << s->Integral() << "\t" << s->GetEntries() << std::endl; if(d->Integral()==0 && s->Integral()==0){ @@ -791,13 +805,19 @@ TCanvas *PlotDataMCs(TChain *data, std::vector mc_vec, TString branchn if(mcLabel_vec.size()!=0) leg->Draw(); TPaveText *pv = new TPaveText(0.25,0.95,0.65,1,"NDC"); - pv->AddText("CMS Preliminary 2012"); + pv->AddText("CMS Preliminary 2015"); pv->SetFillColor(0); pv->SetBorderSize(0); pv->Draw(); watch.Stop(); watch.Print(); + + c->SaveAs(outputPath+label4Print+".png","png"); + c->SaveAs(outputPath+label4Print+".pdf","pdf"); + c->SaveAs(outputPath+label4Print+".eps","eps"); + c->SaveAs(outputPath+label4Print+".C","C"); + return c; } diff --git a/ZFitter/macro/standardDataMC.C b/ZFitter/macro/standardDataMC.C index 7f3b63cbddb..c749021e868 100644 --- a/ZFitter/macro/standardDataMC.C +++ b/ZFitter/macro/standardDataMC.C @@ -1,29 +1,31 @@ { - std::vector mcLabel_vec; - mcLabel_vec.push_back("Madgraph"); - mcLabel_vec.push_back("Powheg"); - mcLabel_vec.push_back("Sherpa"); - - - // pileup - c = PlotDataMC(data, signal, "nPV", "(30,0,30)", ""+commonCut, "", dataLabel+*l_itr, mcLabel, "nVtx", ""); - c->SaveAs(outputPath+"nPV.eps"); - delete c; - -// c = PlotDataMCs(data, MakeChainVector(signalA,signalB,signalC), "etaEle", "(100,-2.5,2.5)", "eleID_loose-trigger-noPF-Et_25", "", "", mcLabel_vec, "#eta", "", false, true, true,false,false,false); -// c->SaveAs(outputPath+"etaEle-Et_25.eps"); -// //c->SaveAs(outputPath+"etaEle-Et_25.png"); -// c->SaveAs(outputPath+"etaEle-Et_25.C"); -// delete c; - - -} - - + int i=0; + TCanvas *c[200]; + TString dataLabel="data"; + std::vector mcLabel_vec; + mcLabel_vec.push_back("MC"); + // mcLabel_vec.push_back("Powheg"); + // mcLabel_vec.push_back("Sherpa"); + + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "invMass_SC", "(30,60,120)", "EB", "", dataLabel,mcLabel_vec, "invMass_SC", "",outputPath,"invMass_SC",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "esEnergySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0,0.2)", "EB", "", dataLabel,mcLabel_vec, "ES energy fraction: ES/(rawSC+ES)", "",outputPath,"esEnergyFraction",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EB-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES})", "",outputPath,"energyCorrections_EB",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EE-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "",outputPath,"energyCorrections_EE",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "nPV", "(30,0,30)", "eleID_7-Et_25", "", dataLabel,mcLabel_vec, "nVtx ", "",outputPath,"nPV",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "rho", "(60,-1,29)", "eleID_7-Et_25", "", dataLabel,mcLabel_vec, "#rho ", "",outputPath,"rho",false,false)); i++; + // kinematic variables + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "energySCEle", "(100,0,200)", "eleID_7", "", dataLabel,mcLabel_vec, "energy SC [GeV] ", "",outputPath,"energySCEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "etaEle", "(100,-2.5,2.5)", "eleID_7-Et_25", "", dataLabel,mcLabel_vec, "#eta ", "",outputPath,"etaEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "phiEle", "(100,-3.14,3.14)", "eleID_7-Et_25", "", dataLabel,mcLabel_vec, "#phi ", "",outputPath,"phiEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "R9Ele", "(100,0.3,1.4)", "EB-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EB",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "R9Ele", "(100,0.3,1.4)", "EB-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EB-log",true,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "R9Ele", "(100,0.3,1.4)", "EE-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EE",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(data,MakeChainVector(signalA), "R9Ele", "(100,0.3,1.4)", "EE-eleID_7-Et_25", "", dataLabel,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EE-log",true,false)); i++; +/* TString dataLabel="22Jan ", mcLabel="Simulation"; std::vector runPeriods, labels; @@ -41,12 +43,13 @@ p_itr != runPeriods.end(); p_itr++,l_itr++){ int index=p_itr-runPeriods.begin(); + } TString commonCut="eleID_WP80PU-Et_25-noPF-trigger"+*p_itr; //TString outputPath="tmp/"; //test/dato/rereco/rereco29Jun-RUN2011/rereco29Jun-RUN2011/WP80_PU/img/"; //TString outputPath="test/dato/moriond2013/WP80_PU/img/"; outputPath=outPath+"/"+commonCut+"/"; - + system(("mkdir -p "+outputPath).Data()); c = PlotDataMC(data, signal, "nHitsSCEle", "(100,0,100)", "EB-"+commonCut,"", dataLabel+*l_itr, mcLabel, "nHits SC", ""); @@ -142,26 +145,6 @@ c->SaveAs(outputPath+"R9Ele_EE-log.eps"); // Checking fit selection -} -} - - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-"+commonCut,"", dataLabel+*l_itr, mcLabel, "Reference Region map", "",1); - c->SaveAs(outputPath+"EERefRegMap.eps"); - delete c; -} - - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(171,-85,86,361,0,361)", "EBRefReg-gold-"+commonCut,"", dataLabel+*l_itr, mcLabel, "Reference Region map", "",1); - c->SaveAs(outputPath+"EBRefRegMap.eps"); - delete c; - - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-"+commonCut,"", dataLabel+*l_itr, mcLabel, "Reference Region map", "",1); - c->SaveAs(outputPath+"EERefRegMap.eps"); - delete c; - - //PlotDataMC(data, signal, "energySCEle", "(100,0,200)", "", dataLabel+*l_itr, "DY Summer 12", "energy SC [GeV]", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)<1.5", dataLabel+*l_itr, "DY Summer 12", "R9", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)>1.5", dataLabel+*l_itr, "DY Summer 12", "R9", ""); -// PlotDataMC(data, signal, "R9Ele:etaEle", "(100,-3,3,100,0,1.1)", "", dataLabel+*l_itr, "DY Summer 12", "#eta", "R9"); -} +*/ } diff --git a/ZFitter/macro/standardMCMC.C b/ZFitter/macro/standardMCMC.C index 3bfa6f2c2c3..76c1a60cb69 100644 --- a/ZFitter/macro/standardMCMC.C +++ b/ZFitter/macro/standardMCMC.C @@ -1,206 +1,101 @@ { - //TString outputPath="tmp/"; //test/dato/rereco/rereco29Jun-RUN2011/rereco29Jun-RUN2011/WP80_PU/img/"; - //TString outputPath="test/dato/moriond2013/WP80_PU/img/"; - TString mc1Label="standard Madgraph"; - TString mc2Label="runDep Madgraph"; - - // regression comparison - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_ele/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCEle_stdSC-EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_pho/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Photon (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCPho_stdSC-EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Egamma) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrEleEgamma_stdSC-EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_ele/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCEle_stdSC-EE.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_pho/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Photon (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCPho_stdSC-EE.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Egamma) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrEleEgamma_stdSC-EE.eps"); - delete c; - - - // cluster corrections and ES energy fraction - c = PlotDataMC(data, signal, "esEnergySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0,0.2)", "absEta_1.7_2.5-eleID_7-Et_25","", mc1Label,mc2Label, "ES energy fraction: ES/(rawSC+ES) ", "", false, false); - c->SaveAs(outputPath+"/esEnergyFraction.eps"); - delete c; - - c = PlotDataMC(data, signal, "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false); - c->SaveAs(outputPath+"/energyCorrections_EB.eps"); - delete c; + int i=0; - c = PlotDataMC(data, signal, "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false); - c->SaveAs(outputPath+"/energyCorrections_EE.eps"); - delete c; - - // pileup - c = PlotDataMC(data, signal, "nPV", "(30,0,30)", "eleID_7-Et_25", "", mc1Label,mc2Label, "nVtx", "", false, false); - c->SaveAs(outputPath+"nPV.eps"); - delete c; - - c = PlotDataMC(data, signal, "rho", "(60,-1,29)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#rho", "", false, false); - c->SaveAs(outputPath+"rho.eps"); - delete c; + TCanvas *c[200]; + //TString outputPath="tmp/"; //test/dato/rereco/rereco29Jun-RUN2011/rereco29Jun-RUN2011/WP80_PU/img/"; + //TString outputPath="test/dato/moriond2013/WP80_PU/img/"; + TString mc1Label="Asympt50ns"; + TString mc2Label="StartupFlat10To50"; + + std::vector mcLabel_vec; // create a vector of labels for the MC samples + mcLabel_vec.push_back("StartupFlat10To50"); + //mcLabel_vec.push_back("Powheg"); + //mcLabel_vec.push_back("Sherpa"); + + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "invMass_SC", "(30,60,120)", "EB", "", mc1Label,mcLabel_vec, "invMass_SC", "",outputPath,"invMass_SC",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "esEnergySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0,0.2)", "EB", "", mc1Label,mcLabel_vec, "ES energy fraction: ES/(rawSC+ES)", "",outputPath,"esEnergyFraction",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EB-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES})", "",outputPath,"energyCorrections_EB",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EE-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "",outputPath,"energyCorrections_EE",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "nPV", "(30,0,30)", "eleID_7-Et_25", "", mc1Label,mcLabel_vec, "nVtx ", "",outputPath,"nPV",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "rho", "(60,-1,29)", "eleID_7-Et_25", "", mc1Label,mcLabel_vec, "#rho ", "",outputPath,"rho",false,false)); i++; // kinematic variables - c = PlotDataMC(data, signal, "energySCEle", "(100,0,200)", "eleID_7", "", mc1Label,mc2Label, "energy SC [GeV]", "", false, false); - c->SaveAs(outputPath+"energySCEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "etaEle", "(100,-2.5,2.5)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#eta", "", false, false); - c->SaveAs(outputPath+"etaEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "phiEle", "(100,-5,5)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#phi", "", false, false); - c->SaveAs(outputPath+"phiEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "", false, false); - c->SaveAs(outputPath+"R9Ele_EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "",true, false); + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "energySCEle", "(100,0,200)", "eleID_7", "", mc1Label,mcLabel_vec, "energy SC [GeV] ", "",outputPath,"energySCEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "etaEle", "(100,-2.5,2.5)", "eleID_7-Et_25", "", mc1Label,mcLabel_vec, "#eta ", "",outputPath,"etaEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "phiEle", "(100,-3.14,3.14)", "eleID_7-Et_25", "", mc1Label,mcLabel_vec, "#phi ", "",outputPath,"phiEle",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "R9Ele", "(100,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EB",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "R9Ele", "(100,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EB-log",true,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "R9Ele", "(100,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EE",false,false)); i++; + c[i] = new TCanvas (PlotDataMCs(signalA,MakeChainVector(signalB), "R9Ele", "(100,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9} ", "",outputPath,"R9Ele_EE-log",true,false)); i++; + +/* + + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele", "(110,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9}", "",true, false)); c->SaveAs(outputPath+"R9Ele_EB-log.eps"); delete c; - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9}", "", false, false)); c->SaveAs(outputPath+"R9Ele_EE.eps"); delete c; - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "",true, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mcLabel_vec, "R_{9}", "",true, false)); c->SaveAs(outputPath+"R9Ele_EE-log.eps"); // Checking fit selection -} - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); + c = new TCanvas( PlotDataMC2D(signalA, MakeChainVector(signalB), "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Reference Region map", "",1, false, false)); c->SaveAs(outputPath+"EERefRegMap.eps"); delete c; -} - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(171,-85,86,361,0,361)", "EBRefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); + + c = new TCanvas( PlotDataMC2D(signalA, MakeChainVector(signalB), "seedYSCEle:seedXSCEle", "(171,-85,86,361,0,361)", "EBRefReg-gold-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Reference Region map", "",1, false, false)); c->SaveAs(outputPath+"EBRefRegMap.eps"); delete c; - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); + c = new TCanvas( PlotDataMC2D(signalA, MakeChainVector(signalB), "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Reference Region map", "",1, false, false)); c->SaveAs(outputPath+"EERefRegMap.eps"); delete c; - //PlotDataMC(data, signal, "energySCEle", "(100,0,200)", "", mc1Label, "DY Summer 12", "energy SC [GeV]", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)<1.5", mc1Label, "DY Summer 12", "R9", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)>1.5", mc1Label, "DY Summer 12", "R9", ""); -// PlotDataMC(data, signal, "R9Ele:etaEle", "(100,-3,3,100,0,1.1)", "", mc1Label, "DY Summer 12", "#eta", "R9"); -} + //PlotDataMCs(signalA, MakeChainVector(signalB), "energySCEle", "(100,0,200)", "", mc1Label, "DY Summer 12", "energy SC [GeV]", ""); + //PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele", "(100,0,1.1)", "abs(etaEle)<1.5", mc1Label, "DY Summer 12", "R9", ""); + //PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele", "(100,0,1.1)", "abs(etaEle)>1.5", mc1Label, "DY Summer 12", "R9", ""); +// PlotDataMCs(signalA, MakeChainVector(signalB), "R9Ele:etaEle", "(100,-3,3,100,0,1.1)", "", mc1Label, "DY Summer 12", "#eta", "R9"); -{ //TString outputPath="tmp/"; //test/dato/rereco/rereco29Jun-RUN2011/rereco29Jun-RUN2011/WP80_PU/img/"; //TString outputPath="test/dato/moriond2013/WP80_PU/img/"; TString mc1Label="standard Madgraph"; - TString mc2Label="runDep Madgraph"; + TString mcLabel_vec="runDep Madgraph"; - // regression comparison - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_ele/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCEle_stdSC-EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_pho/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Photon (Hgg) regression/std. SC -1 ", "", false, false); - c->SaveAs(outputPath+"/regrSCPho_stdSC-EB.eps"); - delete c; - c = PlotDataMC(data, signal, "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Egamma) regression/std. SC -1 ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.2,0.2)", "EB-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Electron (Egamma) regression/std. SC -1 ", "", false, false)); c->SaveAs(outputPath+"/regrEleEgamma_stdSC-EB.eps"); delete c; - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_ele/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Hgg) regression/std. SC -1 ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "invMass_SC_regrCorr_ele/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Electron (Hgg) regression/std. SC -1 ", "", false, false)); c->SaveAs(outputPath+"/regrSCEle_stdSC-EE.eps"); delete c; - c = PlotDataMC(data, signal, "invMass_SC_regrCorr_pho/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Photon (Hgg) regression/std. SC -1 ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "invMass_SC_regrCorr_pho/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Photon (Hgg) regression/std. SC -1 ", "", false, false)); c->SaveAs(outputPath+"/regrSCPho_stdSC-EE.eps"); delete c; - c = PlotDataMC(data, signal, "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Electron (Egamma) regression/std. SC -1 ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "invMass_regrCorr_egamma/invMass_SC -1", "(500,-0.1,0.1)", "EE-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Electron (Egamma) regression/std. SC -1 ", "", false, false)); c->SaveAs(outputPath+"/regrEleEgamma_stdSC-EE.eps"); delete c; // cluster corrections and ES energy fraction - c = PlotDataMC(data, signal, "esEnergySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0,0.2)", "absEta_1.7_2.5-eleID_7-Et_25","", mc1Label,mc2Label, "ES energy fraction: ES/(rawSC+ES) ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "esEnergySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0,0.2)", "absEta_1.7_2.5-eleID_7-Et_25","", mc1Label,mcLabel_vec, "ES energy fraction: ES/(rawSC+ES) ", "", false, false)); c->SaveAs(outputPath+"/esEnergyFraction.eps"); delete c; - c = PlotDataMC(data, signal, "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EB-eleID_7-Et_25","", mc1Label,mc2Label, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EB-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false)); c->SaveAs(outputPath+"/energyCorrections_EB.eps"); delete c; - c = PlotDataMC(data, signal, "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EE-eleID_7-Et_25","", mc1Label,mc2Label, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false); + c = new TCanvas( PlotDataMCs(signalA, MakeChainVector(signalB), "energySCEle/(rawEnergySCEle+esEnergySCEle)", "(100,0.99,1.15)", "EE-eleID_7-Et_25","", mc1Label,mcLabel_vec, "Energy corrections F: E_{SC}/(E_{rawSC}+E_{ES}) ", "", false, false)); c->SaveAs(outputPath+"/energyCorrections_EE.eps"); delete c; - // pileup - c = PlotDataMC(data, signal, "nPV", "(30,0,30)", "eleID_7-Et_25", "", mc1Label,mc2Label, "nVtx", "", false, false); - c->SaveAs(outputPath+"nPV.eps"); - delete c; - - c = PlotDataMC(data, signal, "rho", "(60,-1,29)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#rho", "", false, false); - c->SaveAs(outputPath+"rho.eps"); - delete c; - - // kinematic variables - c = PlotDataMC(data, signal, "energySCEle", "(100,0,200)", "eleID_7", "", mc1Label,mc2Label, "energy SC [GeV]", "", false, false); - c->SaveAs(outputPath+"energySCEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "etaEle", "(100,-2.5,2.5)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#eta", "", false, false); - c->SaveAs(outputPath+"etaEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "phiEle", "(100,-5,5)", "eleID_7-Et_25", "", mc1Label,mc2Label, "#phi", "", false, false); - c->SaveAs(outputPath+"phiEle.eps"); - delete c; - - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "", false, false); - c->SaveAs(outputPath+"R9Ele_EB.eps"); - delete c; - - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EB-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "",true, false); - c->SaveAs(outputPath+"R9Ele_EB-log.eps"); - delete c; - - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "", false, false); - c->SaveAs(outputPath+"R9Ele_EE.eps"); - delete c; - c = PlotDataMC(data, signal, "R9Ele", "(110,0.3,1.4)", "EE-eleID_7-Et_25", "", mc1Label,mc2Label, "R_{9}", "",true, false); - c->SaveAs(outputPath+"R9Ele_EE-log.eps"); - - // Checking fit selection -} - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); - c->SaveAs(outputPath+"EERefRegMap.eps"); - delete c; - -} - - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(171,-85,86,361,0,361)", "EBRefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); - c->SaveAs(outputPath+"EBRefRegMap.eps"); - delete c; - - c = PlotDataMC2D(data, signal, "seedYSCEle:seedXSCEle", "(100,0,100,100,0,100)", "EERefReg-gold-eleID_7-Et_25","", mc1Label,mc2Label, "Reference Region map", "",1, false, false); - c->SaveAs(outputPath+"EERefRegMap.eps"); - delete c; +*/ - //PlotDataMC(data, signal, "energySCEle", "(100,0,200)", "", mc1Label, "DY Summer 12", "energy SC [GeV]", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)<1.5", mc1Label, "DY Summer 12", "R9", ""); - //PlotDataMC(data, signal, "R9Ele", "(100,0,1.1)", "abs(etaEle)>1.5", mc1Label, "DY Summer 12", "R9", ""); -// PlotDataMC(data, signal, "R9Ele:etaEle", "(100,-3,3,100,0,1.1)", "", mc1Label, "DY Summer 12", "#eta", "R9"); } - diff --git a/ZFitter/script/GenRootChain.sh b/ZFitter/script/GenRootChain.sh index 1207b813728..08ebc5084bd 100755 --- a/ZFitter/script/GenRootChain.sh +++ b/ZFitter/script/GenRootChain.sh @@ -2,11 +2,11 @@ #tag_name="" -commonCut=Et_25-trigger-noPF -selection=loose -invMass_var=invMass_SC_regrCorr_ele -configFile=data/validation/monitoring_2012_53X.dat -regionsFile=data/regions/scaleStep2smearing_9.dat +commonCut=Et_25 +selection= +invMass_var=invMass_SC +configFile=data/validation/monitoring_2015.dat +regionsFile=data/regions/validation.dat runRangesFile=data/runRanges/monitoring.dat baseDir=test @@ -139,7 +139,7 @@ fi # saving the root files with the chains rm tmp/*_chain.root -./bin/ZFitter.exe --saveRootMacro -f ${configFile} --regionsFile=${regionsFile} ${noPU} ${addBranchList} ${corrEleFile} ${corrEleType} ${fitterOptions} || exit 1 +ZFitter.exe --saveRootMacro -f ${configFile} --regionsFile=${regionsFile} ${noPU} ${addBranchList} ${corrEleFile} ${corrEleType} ${fitterOptions} || exit 1 # adding all the chains in one file for file in tmp/s[0-9]*_selected_chain.root tmp/d_selected_chain.root tmp/s_selected_chain.root @@ -156,20 +156,20 @@ done cat > tmp/load.C <ProcessLine(".L macro/PlotDataMC.C+"); - gROOT->ProcessLine(".L src/setTDRStyle.C"); +// gROOT->ProcessLine(".L src/setTDRStyle.C"); gROOT->ProcessLine(".L macro/addChainWithFriends.C+"); - setTDRStyle(); +// setTDRStyle(); - TChain *data = (TChain *) _file0->Get("selected"); - TChain *signalA = (TChain *) _file1->Get("selected"); - TChain *signalB = (TChain *) _file2->Get("selected"); + TChain *signalA = (TChain *) _file0->Get("selected"); + TChain *signalB = (TChain *) _file1->Get("selected"); + TChain *data = (TChain *) _file2->Get("selected"); //TChain *signalC = (TChain *) _file3->Get("selected"); - ReassociateFriends(_file0, data); - ReassociateFriends(_file1, signalA); - ReassociateFriends(_file2, signalB); + ReassociateFriends(_file0, signalA); + ReassociateFriends(_file1, signalB); + ReassociateFriends(_file2, data); // ReassociateFriends(_file3, signalC); TDirectory *curDir = new TDirectory("curDir",""); @@ -183,7 +183,9 @@ cat > tmp/load.C < @@ -14,7 +16,7 @@ ElectronCategory_class::~ElectronCategory_class(){ ElectronCategory_class::ElectronCategory_class(bool isRooFit_, bool roofitNameAsNtuple_): _isRooFit(isRooFit_), _roofitNameAsNtuple(roofitNameAsNtuple_), - energyBranchName("energySCEle_regrCorrSemiParV5_pho"), + energyBranchName("energySCEle"), _corrEle(false){ return; @@ -91,6 +93,8 @@ std::set ElectronCategory_class::GetCutSet(TString region){ TCut cut_string; cut_string.Clear(); + eleIDMap eleID_map; + std::set cutSet; // events %2 == 0 are dropped TCut odd_cut = "eventNumber%2 !=0"; @@ -695,6 +699,8 @@ std::set ElectronCategory_class::GetCutSet(TString region){ TObjString *Objstring1 = (TObjString *) splitted->At(1); TString string1 = Objstring1->GetString(); + string1.Form("%d",eleID_map.eleIDmap[string1.Data()]); + /* if(string1=="loose") string1="2"; else if(string1=="medium") string1="6"; else if(string1=="tight") string1="14"; @@ -706,9 +712,10 @@ std::set ElectronCategory_class::GetCutSet(TString region){ else if(string1=="loose50nsRun2") string1="1024"; else if(string1=="medium50nsRun2") string1="3072"; else if(string1=="tight50nsRun2") string1="7168"; + */ - TCut cutEle1("(eleID_ele1 & "+string1+")=="+string1); - TCut cutEle2("(eleID_ele2 & "+string1+")=="+string1); + TCut cutEle1("(eleID[0] & "+string1+")=="+string1); + TCut cutEle2("(eleID[1] & "+string1+")=="+string1); cut_string+=cutEle1 && cutEle2; cutSet.insert(TString(cutEle1 && cutEle2)); 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 + + + + + + diff --git a/ZNtupleDumper/interface/eleIDMap.h b/ZNtupleDumper/interface/eleIDMap.h new file mode 100644 index 00000000000..1fe28f546bc --- /dev/null +++ b/ZNtupleDumper/interface/eleIDMap.h @@ -0,0 +1,33 @@ +#ifndef __eleIDMap__ +#define __eleIDMap__ + +#include +#include + +class eleIDMap { + + public: + std::map eleIDmap; + + eleIDMap() { + + eleIDmap["fiducial"] =0x0001; + eleIDmap["loose"] =0x0002; + eleIDmap["medium"] =0x0004; + eleIDmap["tight"] =0x0008; + eleIDmap["WP90PU"] =0x0010; + + eleIDmap["WP80PU"] =0x0020; + eleIDmap["WP70PU"] =0x0040; + eleIDmap["loose25nsRun2"] =0x0080; + eleIDmap["medium25nsRun2"] =0x0100; + eleIDmap["tight25nsRun2"] =0x0200; + eleIDmap["loose50nsRun2"] =0x0400; + eleIDmap["medium50nsRun2"] =0x0800; + eleIDmap["tight50nsRun2"] =0x1000; + } + +}; + +#endif + diff --git a/ZNtupleDumper/plugins/ZNtupleDumper.cc b/ZNtupleDumper/plugins/ZNtupleDumper.cc index fd60cea3663..91449aab884 100644 --- a/ZNtupleDumper/plugins/ZNtupleDumper.cc +++ b/ZNtupleDumper/plugins/ZNtupleDumper.cc @@ -31,6 +31,7 @@ #include #include #include +#include // root include files #include @@ -127,6 +128,8 @@ // alcaSkimPaths #include "DataFormats/Provenance/interface/ParameterSetID.h" +#include "Calibration/ZNtupleDumper/interface/eleIDMap.h" + //#define DEBUG //////////////////////////////////////////////// @@ -237,7 +240,7 @@ class ZNtupleDumper : public edm::EDAnalyzer { Int_t nPU[5]; //[nBX] ///< number of PU (filled only for MC) // selection - Int_t eleID[3]; ///< bit mask for eleID: 1=fiducial, 2=loose, 6=medium, 14=tight, 16=WP90PU, 48=WP80PU, 112=WP70PU, 128=loose25nsRun2, 384=medium25nsRun2, 896=tight25nsRun2, 1024=loose50nsRun2, 3072=medium50nsRun2, 7168=tight50nsRun2. Selection from https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification#Electron_ID_Working_Points + UInt_t eleID[3]; ///< bit mask for eleID: 1=fiducial, 2=loose, 6=medium, 14=tight, 16=WP90PU, 48=WP80PU, 112=WP70PU, 128=loose25nsRun2, 384=medium25nsRun2, 896=tight25nsRun2, 1024=loose50nsRun2, 3072=medium50nsRun2, 7168=tight50nsRun2. Selection from https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification#Electron_ID_Working_Points Int_t chargeEle[3]; ///< -100: no electron, 0: SC or photon, -1 or +1:electron or muon Float_t etaSCEle[3], phiSCEle[3]; ///< phi of the SC @@ -315,7 +318,7 @@ class ZNtupleDumper : public edm::EDAnalyzer { std::vector energyRecHitSCEle[3]; std::vector LCRecHitSCEle[3]; std::vector ICRecHitSCEle[3]; - std::vector AlphaRecHitSCEle[3]; + float seedLaserAlphaSCEle[3]; //============================== //============================== check ele-id and iso @@ -401,6 +404,7 @@ class ZNtupleDumper : public edm::EDAnalyzer { EcalClusterLazyTools *clustertools; // EcalClusterLocal _ecalLocal; + const EcalLaserAlphaMap* theEcalLaserAlphaMap; std::set alcaSkimPathIndexes; edm::ParameterSetID alcaSkimPathID; @@ -512,6 +516,7 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe pEvent = &iEvent; pSetup = &iSetup; + // filling infos runNumber, eventNumber, lumi if( !iEvent.isRealData() ){ @@ -595,6 +600,11 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe clustertools = new EcalClusterLazyTools (iEvent, iSetup, recHitCollectionEBTAG, recHitCollectionEETAG); + edm::ESHandle theEcalLaserAlphas; + iSetup.get().get(theEcalLaserAlphas); + theEcalLaserAlphaMap = theEcalLaserAlphas.product(); + + //------------------------------ electrons if (eventType==ZMMG) { //------------------------------ muons @@ -665,9 +675,9 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe for( pat::ElectronCollection::const_iterator eleIter1 = electronsHandle->begin(); eleIter1 != electronsHandle->end(); eleIter1++){ - if( eleIter1->electronID("tight") ) ++nWP70; - else if( eleIter1->electronID("medium") ) ++nMedium; - else if( eleIter1->electronID("loose") ) ++nWP90; + if( eleIter1->electronID("tight50nsRun2") ) ++nWP70; + else if( eleIter1->electronID("medium50nsRun2") ) ++nMedium; + else if( eleIter1->electronID("loose50nsRun2") ) ++nWP90; } } bool doFill=false; @@ -712,9 +722,10 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe // if(!eleIter1->ecalDriven()){ //to make alcareco/alcarereco ntuples coeherent // continue; // } + if(! (eleIter1->electronID("loose50nsRun2")) ) continue; if(eventType==WENU){ - if(! eleIter1->electronID("tight") ) continue; + if(! (eleIter1->electronID("tight50nsRun2")) ) continue; if( nWP70 != 1 || nWP90 > 0 ) continue; //to be a Wenu event request only 1 ele WP70 in the event // MET/MT selection @@ -746,7 +757,7 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe // } // should exit when eleIter1 == end-1 - //if(! eleIter2->electronID("loose") ) continue; + if(! (eleIter2->electronID("loose50nsRun2")) ) continue; // float mass=(eleIter1->p4()+eleIter2->p4()).mass(); @@ -756,14 +767,20 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe double t2=TMath::Exp(-eleIter2->eta()); double t2q = t2*t2; + + float mass=-99; + + if(eleIter1->parentSuperCluster().isNonnull() && eleIter2->parentSuperCluster().isNonnull()) { + double angle=1- ( (1-t1q)*(1-t2q)+4*t1*t2*cos(eleIter1->phi()-eleIter2->phi()))/( (1+t1q)*(1+t2q) ); - float mass = sqrt(2*eleIter1->parentSuperCluster()->energy()*eleIter2->parentSuperCluster()->energy() *angle); //use mustache SC, in order to have the same number of events between alcareco and alcarereco ntuples + mass = sqrt(2*eleIter1->parentSuperCluster()->energy()*eleIter2->parentSuperCluster()->energy() *angle); //use mustache SC, in order to have the same number of events between alcareco and alcarereco ntuples // std::cout<<" ele1 SC: "<superCluster()->energy()<<" ele1 SC must: "<parentSuperCluster()->energy()<<" eta1: "<eta()<<" phi1: "<phi()<superCluster()->energy()<<" ele2 SC must: "<parentSuperCluster()->energy()<<" eta2: "<eta()<<" phi2: "<phi()<<"mass: "< 125)) continue; doFill=true; @@ -857,14 +874,14 @@ void ZNtupleDumper::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe PatEle1++){ // consider electrons passing at least the loose identification - if(!PatEle1->electronID("loose") ) continue; + if(!PatEle1->electronID("loose50nsRun2") ) continue; // take the highest pt tight electrons if it exists (the collection is ordered in pt) // consider only the electrons passing the tightest electron identification if(nWP70>0){ // if there are tight electrons, consider only those - if(!PatEle1->electronID("tight") ) continue; + if(!PatEle1->electronID("tight50nsRun2") ) continue; }else if(nMedium>0){ // if there are only medium electrons, consider only those - if(!PatEle1->electronID("medium") ) continue; + if(!PatEle1->electronID("medium50nsRun2") ) continue; } // if(!PatEle1->ecalDriven()){ //to make alcareco/alcarereco ntuples coeherent @@ -1000,6 +1017,7 @@ void ZNtupleDumper::endJob() if(tree->GetEntries()>0){ tree->BuildIndex("runNumber","eventNumber"); if(doEleIDTree) eleIDTree->BuildIndex("runNumber","eventNumber"); + if(doExtraCalibTree) extraCalibTree->BuildIndex("runNumber","eventNumber"); } // save the tree into the file tree_file->cd(); @@ -1007,13 +1025,11 @@ void ZNtupleDumper::endJob() tree_file->Close(); if(doExtraCalibTree){ - extraCalibTree->BuildIndex("runNumber","eventNumber"); extraCalibTreeFile->cd(); extraCalibTree->Write(); extraCalibTreeFile->Close(); } if(doEleIDTree){ - eleIDTree->BuildIndex("runNumber","eventNumber"); eleIDTreeFile->cd(); eleIDTree->Write(); eleIDTreeFile->Close(); @@ -1087,7 +1103,7 @@ void ZNtupleDumper::InitNewTree(){ tree->Branch("nPV", &nPV, "nPV/I"); - tree->Branch("eleID",eleID, "eleID[3]/I"); + tree->Branch("eleID",eleID, "eleID[3]/i"); // tree->Branch("nBCSCEle", nBCSCEle, "nBCSCEle[3]/I"); tree->Branch("chargeEle", chargeEle, "chargeEle[3]/I"); //[nEle] @@ -1357,7 +1373,8 @@ void ZNtupleDumper::TreeSetSingleElectronVar(const pat::Electron& electron1, int } energySCEle[index] = electron1.superCluster()->energy(); - energySCEle_must[index] = electron1.parentSuperCluster()->energy(); + if(electron1.parentSuperCluster().isNonnull()) + energySCEle_must[index] = electron1.parentSuperCluster()->energy(); rawEnergySCEle[index] = electron1.superCluster()->rawEnergy(); esEnergySCEle[index] = electron1.superCluster()->preshowerEnergy(); #ifndef CMSSW42X @@ -1386,21 +1403,13 @@ void ZNtupleDumper::TreeSetSingleElectronVar(const pat::Electron& electron1, int // R9Ele[index] = R9Ele[index]*1.0086-0.0007; // } - // make it a function - eleID[index] = ((bool) electron1.electronID("fiducial")) << 0; - eleID[index] += ((bool) electron1.electronID("loose")) << 1; - eleID[index] += ((bool) electron1.electronID("medium")) << 2; - eleID[index] += ((bool) electron1.electronID("tight")) << 3; - eleID[index] += ((bool) electron1.electronID("WP90PU")) << 4; - eleID[index] += ((bool) electron1.electronID("WP80PU")) << 5; - eleID[index] += ((bool) electron1.electronID("WP70PU")) << 6; - //LUCA: need to decide if to put the run2 selection here (bits 7-12) or not. Also, need to modify the ElectronCategory_class according to this.. - eleID[index] += ((bool) electron1.electronID("loose25nsRun2")) << 7; - eleID[index] += ((bool) electron1.electronID("medium25nsRun2")) << 8; - eleID[index] += ((bool) electron1.electronID("tight25nsRun2")) << 9; - eleID[index] += ((bool) electron1.electronID("loose50nsRun2")) << 10; - eleID[index] += ((bool) electron1.electronID("medium50nsRun2")) << 11; - eleID[index] += ((bool) electron1.electronID("tight50nsRun2")) << 12; + eleIDMap eleID_map; + + eleID[index]=0; + for (std::map::iterator it=eleID_map.eleIDmap.begin(); it!=eleID_map.eleIDmap.end(); ++it) { + if ((bool) electron1.electronID(it->first)) + eleID[index] |= it->second; + } classificationEle[index] = electron1.classification(); @@ -1645,7 +1654,8 @@ void ZNtupleDumper::TreeSetSinglePhotonVar(const pat::Photon& photon, int index) } energySCEle[index] = photon.superCluster()->energy(); - energySCEle_must[index] = photon.parentSuperCluster()->energy(); + if(photon.parentSuperCluster().isNonnull()) + energySCEle_must[index] = photon.parentSuperCluster()->energy(); rawEnergySCEle[index] = photon.superCluster()->rawEnergy(); esEnergySCEle[index] = photon.superCluster()->preshowerEnergy(); // energySCEle_corr[index] = photon.scEcalEnergy(); //but, I don't think this is the correct energy.. @@ -1866,8 +1876,10 @@ void ZNtupleDumper::InitExtraCalibTree(){ extraCalibTree->Branch("LCRecHitSCEle2", &(LCRecHitSCEle[1])); extraCalibTree->Branch("ICRecHitSCEle1", &(ICRecHitSCEle[0])); extraCalibTree->Branch("ICRecHitSCEle2", &(ICRecHitSCEle[1])); - extraCalibTree->Branch("AlphaRecHitSCEle1", &(AlphaRecHitSCEle[0])); - extraCalibTree->Branch("AlphaRecHitSCEle2", &(AlphaRecHitSCEle[1])); + // extraCalibTree->Branch("AlphaRecHitSCEle1", &(AlphaRecHitSCEle[0])); + //extraCalibTree->Branch("AlphaRecHitSCEle2", &(AlphaRecHitSCEle[1])); + extraCalibTree->Branch("seedLaserAlphaSCEle1", &(seedLaserAlphaSCEle[0])); + extraCalibTree->Branch("seedLaserAlphaSCEle2", &(seedLaserAlphaSCEle[1])); return; } @@ -1897,7 +1909,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Electron& electron1, int ind energyRecHitSCEle[-index].clear(); LCRecHitSCEle[-index].clear(); ICRecHitSCEle[-index].clear(); - AlphaRecHitSCEle[-index].clear(); + seedLaserAlphaSCEle[-index] =-99.; return; } @@ -1909,7 +1921,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Electron& electron1, int ind energyRecHitSCEle[index].clear(); LCRecHitSCEle[index].clear(); ICRecHitSCEle[index].clear(); - AlphaRecHitSCEle[index].clear(); + seedLaserAlphaSCEle[index] = -99.; // EcalIntercalibConstantMap icMap = icHandle->get() std::vector< std::pair > hitsAndFractions_ele1 = electron1.superCluster()->hitsAndFractions(); @@ -1965,6 +1977,15 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Electron& electron1, int ind ICRecHitSCEle[index].push_back(icalconst); } + DetId seedDetId = electron1.seed()->seed(); + if(seedDetId.null()){ + seedDetId = findSCseed(*(electron1.superCluster())); + } + + EcalLaserAlphaMap::const_iterator italpha = theEcalLaserAlphaMap->find(seedDetId); + if( italpha != theEcalLaserAlphaMap->end() ) + seedLaserAlphaSCEle[index] = (*italpha); + return; } @@ -1979,7 +2000,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const reco::SuperCluster& electron1, in energyRecHitSCEle[-index].clear(); LCRecHitSCEle[-index].clear(); ICRecHitSCEle[-index].clear(); - AlphaRecHitSCEle[-index].clear(); + seedLaserAlphaSCEle[-index] =-99.; return; } @@ -1991,7 +2012,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const reco::SuperCluster& electron1, in energyRecHitSCEle[index].clear(); LCRecHitSCEle[index].clear(); ICRecHitSCEle[index].clear(); - AlphaRecHitSCEle[index].clear(); + seedLaserAlphaSCEle[index] =-99.; std::vector< std::pair > hitsAndFractions_ele1 = electron1.hitsAndFractions(); nHitsSCEle[index] = hitsAndFractions_ele1.size(); @@ -2071,6 +2092,15 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const reco::SuperCluster& electron1, in ICRecHitSCEle[index].push_back(icalconst); } + + DetId seedDetId = electron1.seed()->seed(); + if(seedDetId.null()){ + seedDetId = findSCseed(electron1); + } + + EcalLaserAlphaMap::const_iterator italpha = theEcalLaserAlphaMap->find(seedDetId); + if( italpha != theEcalLaserAlphaMap->end() ) + seedLaserAlphaSCEle[index] = (*italpha); return; } @@ -2096,7 +2126,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Photon& photon, int index){ energyRecHitSCEle[-index].clear(); LCRecHitSCEle[-index].clear(); ICRecHitSCEle[-index].clear(); - AlphaRecHitSCEle[-index].clear(); + seedLaserAlphaSCEle[-index] =-99.; return; } @@ -2108,7 +2138,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Photon& photon, int index){ energyRecHitSCEle[index].clear(); LCRecHitSCEle[index].clear(); ICRecHitSCEle[index].clear(); - AlphaRecHitSCEle[index].clear(); + seedLaserAlphaSCEle[index] =-99.; // EcalIntercalibConstantMap icMap = icHandle->get() std::vector< std::pair > hitsAndFractions_ele1 = photon.superCluster()->hitsAndFractions(); @@ -2164,6 +2194,15 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Photon& photon, int index){ ICRecHitSCEle[index].push_back(icalconst); } + DetId seedDetId = photon.seed()->seed(); + if(seedDetId.null()){ + seedDetId = findSCseed(*(photon.superCluster())); + } + + EcalLaserAlphaMap::const_iterator italpha = theEcalLaserAlphaMap->find(seedDetId); + if( italpha != theEcalLaserAlphaMap->end() ) + seedLaserAlphaSCEle[index] = (*italpha); + return; } @@ -2178,7 +2217,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Muon& muon1, int index){ energyRecHitSCEle[-index].clear(); LCRecHitSCEle[-index].clear(); ICRecHitSCEle[-index].clear(); - AlphaRecHitSCEle[-index].clear(); + seedLaserAlphaSCEle[-index] =-99.; return; } @@ -2190,7 +2229,7 @@ void ZNtupleDumper::TreeSetExtraCalibVar(const pat::Muon& muon1, int index){ energyRecHitSCEle[index].clear(); LCRecHitSCEle[index].clear(); ICRecHitSCEle[index].clear(); - AlphaRecHitSCEle[index].clear(); + seedLaserAlphaSCEle[index] =-99.; return; } @@ -2264,9 +2303,9 @@ void ZNtupleDumper::TreeSetEleIDVar(const pat::Electron& electron1, int index){ // dzvtx[index] = electron1.gsfTrack()->dz(); // } - eleIDloose[index] = electron1.electronID("loose"); - eleIDmedium[index] = electron1.electronID("medium"); - eleIDtight[index] = electron1.electronID("tight"); + eleIDloose[index] = electron1.electronID("loose50nsRun2"); + eleIDmedium[index] = electron1.electronID("medium50nsRun2"); + eleIDtight[index] = electron1.electronID("tight50nsRun2"); return; }