MomentumScaleCalibration/test/Probs_SM.C

0001 /**
0002  * This macro produces the probability distributions that can be used with standalone muons.
0003  */
0004
0005 #include <TH1.h>
0006 #include <TH2.h>
0007 #include <TH3.h>
0008 #include <TH2D.h>
0009 #include <TF1.h>
0010 #include <TProfile.h>
0011 #include <TFile.h>
0012 #include <TTree.h>
0013 #include <TBranch.h>
0014 #include <TLeaf.h>
0015 #include <TCanvas.h>
0016 #include <TPostScript.h>
0017 #include <TLine.h>
0018 #include <TPaveText.h>
0019 #include <TPad.h>
0020 #include <TStyle.h>
0021 #include <TMath.h>
0022 #include <TChain.h>
0023 #include <TProfile2D.h>
0024 #include <TLegend.h>
0025 #include <TROOT.h>
0026 #include <TMath.h>
0027
0028 #include <iostream>
0029 #include <algorithm>
0030
0031 //TF1 * GL = new TF1 ("GL",
0032 //      "0.5/TMath::PI*[0]/(pow(x-[1],2)+pow(0.5*[0],2))*exp(-0.5*pow((x-[2])/[3],2))/([3]*sqrt(2*TMath::PI))",
0033 TF1 * GL = new TF1 ("GL",
0034             "0.5/TMath::Pi()*[0]/((x-[1])*(x-[1])+(0.5*[0])*(0.5*[0]))*exp(-0.5*(x-[2])*(x-[2])/([3]*[3]))/([3]*sqrt(2*TMath::Pi()))",
0035             0, 1000);
0036 TF1 * L = new TF1 ("L", "0.5/TMath::Pi()*[0]/(pow(x-[1],2)+pow(0.5*[0],2))", 0, 1000);
0037 void Probs (int nbins) {
0038
0039   double Gamma[6] = {2.4952, 0.0000934, 0.000337, 0.000054, 0.000032, 0.000020};
0040   double Mass[6] = {90.986, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969};
0041   double ResHalfWidth[6] = {50., 2.5, 2.5, 2.5, 1.5, 1.5};
0042   double ResMaxSigma[6] = {100., 10., 10., 10., 5., 5.};
0043
0044   double Xmin[6];
0045   double Xmax[6];
0046   double Ymin[6];
0047   double Ymax[6];
0048   for (int i=0; i<6; i++) {
0049     Xmin[i] = Mass[i]-ResHalfWidth[i];
0050     Xmax[i] = Mass[i]+ResHalfWidth[i];
0051     Ymin[i] = 0.;
0052     Ymax[i] = ResMaxSigma[i];
0053   }
0054
0055   TH2D * I[6];
0056   char name[20];
0057   Int_t np = 20000;
0058   double * x = new double[np];
0059   double * w = new double[np];
0060
0061   for (int i=0; i<6; i++) {
0062     sprintf (name, "GL%d", i);
0063     I[i] = new TH2D (name, "Gaussian x Lorentz", nbins+1, Xmin[i], Xmax[i], nbins+1, Ymin[i], Ymax[i]);
0064
0065     std::cout << "Processing resonance " << i << std::endl;
0066     for (int ix=0; ix<=nbins/2; ix++) {
0067       double mass = Xmin[i]+(Xmax[i]-Xmin[i])*((double)ix)/(double)nbins;
0068       double sigma;
0069       double P;
0070       for (int iy=0; iy<=nbins; iy++) {
0071     sigma = Ymin[i]+(Ymax[i]-Ymin[i])*((double)iy)/(double)nbins;
0072     P = 0.;
0073     if (iy==0) {
0074       np = 10000;
0075       sigma = 0.1*(ResMaxSigma[i])/(double)nbins; //////////////////////////// NNBB approximation
0076       GL->SetParameters (Gamma[i], Mass[i], mass, sigma);
0077       GL->CalcGaussLegendreSamplingPoints (np, x, w, 0.1e-18);
0078       P = GL->IntegralFast (np, x, w, Mass[i]-10*Gamma[i], Mass[i]+10*Gamma[i]);
0079       std::cout << "For Resonance #" << i << ": mass = " << mass << ", sigma = " << sigma
0080            << ", P = " << P << std::endl;
0081     } else if (iy<10) {
0082       np = 2000;
0083       GL->SetParameters (Gamma[i], Mass[i], mass, sigma);
0084       GL->CalcGaussLegendreSamplingPoints (np, x, w, 0.1e-16);
0085       P = GL->IntegralFast (np, x, w, Mass[i]-10*Gamma[i], Mass[i]+10*Gamma[i]);
0086       // P = GL->Integral(Mass[i]-10*Gamma[i], Mass[i]+10*Gamma[i]);
0087       std::cout << "For Resonance #" << i << ": mass = " << mass << ", sigma = " << sigma
0088            << ", P = " << P << std::endl;
0089     } else {
0090       GL->SetParameters (Gamma[i], Mass[i], mass, sigma);
0091       P = GL->Integral(Mass[i]-10*Gamma[i], Mass[i]+10*Gamma[i]);
0092     }
0093
0094     I[i]->SetBinContent(ix+1, iy+1, P);
0095     if (ix<nbins/2) I[i]->SetBinContent(nbins-ix+1, iy+1, P);
0096       }
0097     }
0098   }
0099
0100   sprintf (name, "Probs_SM_%d.root", nbins);
0101   TFile * Out = new TFile (name, "RECREATE");
0102   Out->cd();
0103   for (int i=0; i<6; i++) {
0104     I[i]->Write();
0105   }
0106   Out->Close();
0107
0108 }
0109