MomentumScaleCalibration/test/Probs.C

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