test/Macros/BackgroundCheck.C

0001 /**
0002  * This macro is used to compare the background function from MuScleFit to the effective background distribution seen in input.
0003  * It needs the result of MuScleFit on three different runs:
0004  * - all events together
0005  * - background only
0006  * - resonance only
0007  * The first run should be on all events together. Once the parameters for the scale fit are found, they can be introduced as
0008  * bias for the following two separate runs so as to get the same corrected muons.
0009  * The following two runs should be made without any fit, they only need to produce the reconstructed mass distributions.
0010  *
0011  * This macro takes the reconstructed mass distributions from the last two runs and shows them on the same plot as the all-events
0012  * distributions. On the same canvas it also shows the background function with the fitted parameters.
0013  */
0014
0015 #include <iostream>
0016 #include "TFile.h"
0017 #include "TH1F.h"
0018 #include "TCanvas.h"
0019 #include "TF1.h"
0020 #include "TLegend.h"
0021 #include <vector>
0022 #include "TString.h"
0023 #include <sstream>
0024
0025
0026 // Linear function as in Functions.h.
0027 // Note that if the type1 background function in Functions.h changes this should be updated too.
0028 class Linear
0029 {
0030  public:
0031   Linear(const double & lowerLimit, const double & upperLimit) : lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
0032
0033   double operator()( const double * x, const double * parval ) const
0034   {
0035     double a = parval[1];
0036     double b = parval[2];
0037
0038     double norm = -(a*lowerLimit_ + b*lowerLimit_*lowerLimit_/2.);
0039
0040     if( -a/b > upperLimit_ ) norm += a*upperLimit_ + b*upperLimit_*upperLimit_/2.;
0041     else norm += -a*a/(2*b);
0042
0043     if( x[0] < -a/b && norm != 0 ) return (a + b*x[0])/norm;
0044     else return 0;
0045   }
0046 protected:
0047   double lowerLimit_, upperLimit_;
0048 };
0049
0050 class Atan
0051 {
0052 public:
0053   Atan(const double & lowerLimit, const double & upperLimit) : lowerLimit_(lowerLimit), upperLimit_(upperLimit) {}
0054
0055   double operator()( const double * x, const double * parval ) const
0056   {
0057     double value = atan(parval[1]*(x[0] + parval[2]));
0058     double upperX = parval[1]*(upperLimit_ + parval[2]);
0059     double lowerX = parval[1]*(lowerLimit_ + parval[2]);
0060     double norm = (upperX*atan(upperX) - 0.5*log(1+upperX*upperX) - (lowerX*atan(lowerX) - 0.5*log(1+lowerX*lowerX)))/parval[1];
0061
0062     if( norm != 0 ) return value/norm;
0063     return 0.;
0064   }
0065 protected:
0066   double lowerLimit_, upperLimit_;
0067 };
0068
0069 TH1F * subRangeHisto( const double * resMass, const double * resHalfWidth,
0070                       const int iRes, TH1F * histo, const TString & name )
0071 {
0072   std::stringstream ss;
0073   ss << iRes;
0074
0075   TH1F* newHisto = (TH1F*)histo->Clone(name+ss.str());
0076   newHisto->SetAxisRange( resMass[iRes] - resHalfWidth[iRes], resMass[iRes] + resHalfWidth[iRes] );
0077   // Workaround for when we fit all the Upsilon resonances together
0078   if( iRes == 3 ) {
0079     newHisto->SetAxisRange( resMass[iRes] - resHalfWidth[iRes], resMass[1] + resHalfWidth[iRes] );
0080   }
0081   return newHisto;
0082 }
0083
0084 TH1F * buildHistogram(const double * ResMass, const double * ResHalfWidth, const int xBins, const double & deltaX, const double & xMin, const double & xMax,
0085                       const int ires, const double & Bgrp1, const double & a, const double & leftWindowBorder, const double & rightWindowBorder,
0086               const TH1F* allHisto, const int backgroundType, const double & b = 0);
0087
0088 void BackgroundCheck()
0089 {
0090   TFile * allFile = new TFile("0_MuScleFit.root", "READ");
0091   TFile * resonanceFile = new TFile("0_MuScleFit.root", "READ");
0092   TFile * backgroundFile = new TFile("0_MuScleFit.root", "READ");
0093
0094   // TH1F * allHisto = (TH1F*)allFile->Get("hRecBestRes_Mass");
0095   TH1F * allHisto = (TH1F*)allFile->Get("hRecBestResAllEvents_Mass");
0096   TH1F * resonanceHisto = (TH1F*)resonanceFile->Get("hRecBestRes_Mass");
0097   TH1F * backgroundHisto = (TH1F*)backgroundFile->Get("hRecBestRes_Mass");
0098
0099   int rebinCount = 4;
0100
0101   allHisto->Rebin(rebinCount);
0102   resonanceHisto->Rebin(rebinCount);
0103   backgroundHisto->Rebin(rebinCount);
0104
0105   double xMin = allHisto->GetXaxis()->GetXmin();
0106   double xMax = allHisto->GetXaxis()->GetXmax();
0107   double deltaX = xMax - xMin;
0108
0109   // The fitted background function gives the background fraction (is normalized).
0110   // We multiply it by the integral to get the background value.
0111   int xBins = allHisto->GetNbinsX();
0112
0113   double OriginalResMass[] = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.096916};
0114   double ResMass[] = {91.1876, 0., 0., (10.3552 + 10.0233 + 9.4603)/3., (3.68609 + 3.096916)/2., (3.68609 + 3.096916)/2.};
0115   double ResHalfWidth[] = {20., 0.5, 0.5, 0.5, 0.2, 0.2};
0116   TString ResName[] = {"Z", "Upsilon3S", "Upsilon2S", "Upsilon1S", "Psi2S", "J/Psi"};
0117
0118   std::vector<int> ires;
0119   std::vector<double> Bgrp1;
0120   std::vector<double> a;
0121   std::vector<double> leftWindowBorder;
0122   std::vector<double> rightWindowBorder;
0123   int backgroundType = 0;
0124
0125   // IMPORTANT: parameters to change
0126   // -------------------------------
0127   ires.push_back(5);
0128
0129   // Exponential
0130   backgroundType = 0;
0131   Bgrp1.push_back(0.730566);
0132   a.push_back(0.826053);
0133   a.push_back(0);
0134   leftWindowBorder.push_back(2.5);
0135   rightWindowBorder.push_back(5.);
0136
0137   // // Linear
0138   // backgroundType.push_back(1);
0139   // Bgrp1.push_back(0.833);
0140   // a.push_back(0.00016);
0141   // a.push_back(-1.27329e-11);
0142   // leftWindowBorder.push_back(1);
0143   // rightWindowBorder.push_back(1);
0144
0145   // // Atan
0146   // backgroundType.push_back(2);
0147   // Bgrp1.push_back(0.000547261);
0148   // a.push_back(10.0145);
0149   // a.push_back(-0.320232);
0150   // leftWindowBorder.push_back(1);
0151   // rightWindowBorder.push_back(1);
0152
0153   // -------------------------------
0154
0155   // Create histograms for the background functions
0156   std::vector<TH1F*> backgroundFunctionHisto;
0157   for( unsigned int i=0; i<ires.size(); ++i ) {
0158     backgroundFunctionHisto.push_back( buildHistogram(ResMass, ResHalfWidth, xBins, deltaX, xMin, xMax,
0159                                                       ires[i], Bgrp1[i], a[i], leftWindowBorder[i], rightWindowBorder[i],
0160                               allHisto, backgroundType, a[i+1]) );
0161   }
0162
0163   TLegend * legend = new TLegend( 0.55, 0.65, 0.76, 0.82 );
0164   TCanvas * canvas = new TCanvas("ResMassCanvas", "ResMassCanvas", 1000, 800);
0165   canvas->cd();
0166   legend->AddEntry(allHisto, "All events");
0167   allHisto->SetLineWidth(2);
0168   allHisto->Draw();
0169
0170   resonanceHisto->SetLineColor(kRed);
0171   resonanceHisto->SetLineWidth(2);
0172   resonanceHisto->Draw("same");
0173   legend->AddEntry(resonanceHisto, "Resonance events");
0174
0175   backgroundHisto->SetLineWidth(2);
0176   backgroundHisto->SetLineColor(kBlue);
0177   backgroundHisto->Draw("same");
0178   legend->AddEntry(backgroundHisto, "Background events");
0179
0180   for( unsigned int i=0; i<ires.size(); ++i ) {
0181     backgroundFunctionHisto[i]->SetLineWidth(2);
0182     if( i == 1 ) backgroundFunctionHisto[i]->SetLineColor(kRed);
0183     else backgroundFunctionHisto[i]->SetLineColor(kGreen);
0184
0185     backgroundFunctionHisto[i]->Draw("same");
0186
0187     TH1F * whiteHisto = subRangeHisto(OriginalResMass, ResHalfWidth, ires[i], backgroundFunctionHisto[i], "whiteHisto");
0188     whiteHisto->SetLineWidth(2);
0189     whiteHisto->SetLineColor(kWhite);
0190     whiteHisto->Draw("same");
0191
0192     TH1F * dashedHisto = subRangeHisto(OriginalResMass, ResHalfWidth, ires[i], backgroundFunctionHisto[i], "dashedHisto");
0193     dashedHisto->SetLineWidth(2);
0194     dashedHisto->SetLineStyle(7);
0195     dashedHisto->Draw("same");
0196
0197     legend->AddEntry(backgroundFunctionHisto[i], "Background function for "+ResName[ires[i]]);
0198   }
0199
0200   legend->Draw("same");
0201
0202   canvas->GetPad(0)->SetLogy(true);
0203
0204   canvas->Print("BackgroundCheck.pdf");
0205 }
0206
0207 TH1F * buildHistogram(const double * ResMass, const double * ResHalfWidth, const int xBins,
0208               const double & deltaX, const double & xMin, const double & xMax,
0209                       const int ires, const double & Bgrp1, const double & a,
0210               const double & leftWindowBorder, const double & rightWindowBorder,
0211               const TH1F* allHisto, const int backgroundType, const double & b)
0212 {
0213   // For J/Psi exclude the Upsilon from the background normalization as the bin is not used by the fit.
0214   double lowWindowValue = leftWindowBorder;
0215   double upWindowValue = rightWindowBorder;
0216
0217   int lowBin = int((lowWindowValue)*xBins/deltaX);
0218   int upBin = int((upWindowValue)*xBins/deltaX);
0219
0220   std::cout << "lowBin = " << lowBin << ", upBin = " << upBin << std::endl;
0221   std::cout << "lowWindowValue = " << lowWindowValue << ", upWindowValue = " << upWindowValue << std::endl;
0222
0223   double xWidth = deltaX/xBins;
0224
0225
0226   TF1 * backgroundFunction = 0;
0227   TH1F * backgroundFunctionHisto = 0;
0228
0229   if( backgroundType == 0 ) {
0230     // Exponential
0231     // -----------
0232     // backgroundFunction = new TF1("backgroundFunction", "[0]*([1]*exp(-[1]*x))", xMin, xMax );
0233     std::stringstream ssUp;
0234     std::stringstream ssDown;
0235     ssUp << upWindowValue;
0236     ssDown << lowWindowValue;
0237     std::string functionString("[0]*(-[1]*exp(-[1]*x)/(exp(-[1]*("+ssUp.str()+")) - exp(-[1]*("+ssDown.str()+")) ))");
0238
0239     backgroundFunction = new TF1("backgroundFunction", functionString.c_str(), xMin, xMax );
0240     backgroundFunction->SetParameter(0, 1);
0241     backgroundFunction->SetParameter(1, a);
0242     backgroundFunctionHisto = new TH1F("backgroundFunctionHisto", "backgroundFunctionHisto", xBins, xMin, xMax);
0243     for( int xBin = 0; xBin < xBins; ++xBin ) {
0244       backgroundFunctionHisto->SetBinContent(xBin+1, backgroundFunction->Integral(xBin*xWidth, (xBin+1)*xWidth));
0245     }
0246   }
0247   else if( backgroundType == 2 ) {
0248
0249     std::stringstream ssUp;
0250     std::stringstream ssDown;
0251     ssUp << upWindowValue;
0252     ssDown << lowWindowValue;
0253
0254     Atan * atanFunction = new Atan(lowWindowValue, upWindowValue);
0255     backgroundFunction = new TF1("backgroundFunction", atanFunction, 0, 1, 3, "atanFunction");
0256
0257     backgroundFunction->SetParameter(0, 1);
0258     backgroundFunction->SetParameter(1, a);
0259     backgroundFunction->SetParameter(2, b);
0260     backgroundFunctionHisto = new TH1F("backgroundFunctionHisto", "backgroundFunctionHisto", xBins, xMin, xMax);
0261     for( int xBin = 0; xBin < xBins; ++xBin ) {
0262       backgroundFunctionHisto->SetBinContent(xBin+1, backgroundFunction->Integral(xBin*xWidth, (xBin+1)*xWidth));
0263     }
0264   }
0265   else {
0266     // Linear
0267     // ------
0268     Linear * linearFunction = new Linear(lowWindowValue, upWindowValue);
0269     backgroundFunction = new TF1("backgroundFunction", linearFunction, 0, 1, 3, "linearFunction");
0270     backgroundFunction->SetParameter(0, 1);
0271     backgroundFunction->SetParameter(1, a);
0272     backgroundFunction->SetParameter(2, b);
0273     backgroundFunctionHisto = new TH1F("backgroundFunctionHisto", "backgroundFunctionHisto", xBins, xMin, xMax);
0274     for( int xBin = 0; xBin < xBins; ++xBin ) {
0275       backgroundFunctionHisto->SetBinContent(xBin+1, backgroundFunction->Integral(xBin*xWidth, (xBin+1)*xWidth));
0276     }
0277   }
0278
0279   // The integral of the background function is 1
0280   // The function must be rescaled multiplying it to the number of events and k = Bgrp1
0281   double totEvents = allHisto->Integral(lowBin, upBin);
0282   backgroundFunctionHisto->Scale(Bgrp1*totEvents);
0283   backgroundFunction->SetParameter(0, Bgrp1*totEvents);
0284
0285   std::cout << "Total events in the background window = " << totEvents << std::endl;
0286   std::cout << "FunctionHisto integral = " << backgroundFunctionHisto->Integral(lowBin, upBin) << std::endl;
0287   std::cout << "Function integral = " << backgroundFunction->Integral(lowWindowValue, upWindowValue) << std::endl;
0288
0289   backgroundFunctionHisto->SetAxisRange(lowWindowValue, upWindowValue);
0290
0291   return backgroundFunctionHisto;
0292 }