PWGDQ/dielectron/AliDielectronSignalExt.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                                                                       //
  18 //                      Dielectron SignalExt                             //
  19 //                                                                       //
  20 /*
  21
  22   Class used for extracting the signal from an invariant mass spectrum.
  23   Used invariant mass spectra are provided via an array of histograms. There are serveral method
  24   to estimate the background and to extract the raw yield from the background subtracted spectra.
  25
  26   Example usage:
  27
  28   AliDielectronSignalExt *sig = new AliDielectronSignalExt();
  29
  30
  31   1) invariant mass input spectra
  32
  33   1.1) Assuming a AliDielectronCF container as data format (check class for more details)
  34   AliDielectronCFdraw *cf = new AliDielectronCFdraw("path/to/the/output/file.root");
  35   TObjArray *arrHists = cf->CollectMinvProj(cf->FindStep("Config"));
  36
  37   1.2) Assuming a AliDielectronHF grid as data format (check class for more details)
  38   AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
  39   TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
  40
  41   1.3) Assuming a single histograms
  42   TObjArray *histoArray = new TObjArray();
  43   arrHists->Add(signalPP);            // add the spectrum histograms to the array
  44   arrHists->Add(signalPM);            // the order is important !!!
  45   arrHists->Add(signalMM);
  46
  47
  48   2) background estimation
  49
  50   2.1) set the method for the background estimation (methods can be found in AliDielectronSignalBase)
  51   sig->SetMethod(AliDielectronSignalBase::kEventMixing);
  52   2.2) rebin the spectras if needed
  53   //  sig->SetRebin(2);
  54   2.3) normalize the backgound spectum to the odd-sign spectrum in the desired range(s)
  55   sig->SetScaleRawToBackground(minScale, maxScale);
  56   //  sig->SetScaleRawToBackground(minScale, maxScale, minScale2, maxScale2);
  57
  58
  59   3) configure the signal extraction
  60
  61   3.1) set the method for the signal extraction (methods can be found in AliDielectronSignalBase)
  62   depending on the method serveral inputs are needed (e.g. MC shape, PDG code of the particle of interest)
  63   //  sig->SetParticleOfInterest(443); //default is jpsi
  64   //  sig->SetMCSignalShape(signalMC);
  65   sig->SetIntegralRange(minInt, maxInt);  // range for bin counting
  66   sig->SetExtractionMethod(AliDielectronSignal::BinCounting); // this is the default
  67
  68
  69   4) start the processing
  70
  71   sig->Process(arrHists);
  72   sig->Print(""); // print values and errors extracted
  73
  74
  75   5) access the spectra and values created
  76
  77   5.1) standard spectras
  78   TH1F *hsign = (TH1F*) sig->GetUnlikeSignHistogram();  // same as the input (rebinned)
  79   TH1F *hbgrd = (TH1F*) sig->GetBackgroundHistogram();  // scaled input      (rebinned)
  80   TH1F *hextr = (TH1F*) sig->GetSignalHistogram();      // after backgound extraction (rebinned)
  81   TObject *oPeak = (TObject*) sig->GetPeakShape();      // can be a TF1 or TH1 depending on the extraction method
  82   TH1F *hrfac = (TH1F*) sig->GetRfactorHistogram();     // if like-sign correction was activated, o.w. 0x0
  83   5.2) access the extracted values and errors
  84   sig->GetValues();     or GetErrors();                 // yield extraction
  85
  86 */
  87 //                                                                       //
  88 ///////////////////////////////////////////////////////////////////////////
  89 #include <TF1.h>
  90 #include <TH1.h>
  91 #include <TH2F.h>
  92 #include <TLatex.h>
  93 #include <TLegend.h>
  94 #include <TCanvas.h>
  95 #include <TMath.h>
  96 #include <TString.h>
  97 #include <TLine.h>
  98
  99 #include <AliLog.h>
 100
 101 #include "AliDielectronSignalExt.h"
 102 #include "AliDielectron.h"
 103
 104 ClassImp(AliDielectronSignalExt)
 105
 106 AliDielectronSignalExt::AliDielectronSignalExt() :
 107   AliDielectronSignalBase()
 108 {
 109   //
 110   // Default Constructor
 111   //
 112 }
 113
 114 //______________________________________________
 115 AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
 116   AliDielectronSignalBase(name, title)
 117 {
 118   //
 119   // Named Constructor
 120   //
 121 }
 122
 123 //______________________________________________
 124 AliDielectronSignalExt::~AliDielectronSignalExt()
 125 {
 126   //
 127   // Default Destructor
 128   //
 129 }
 130
 131 //______________________________________________
 132 void AliDielectronSignalExt::Process(TObjArray* const arrhist)
 133 {
 134   //
 135   // signal subtraction. support like-sign subtraction and event mixing method
 136   //
 137   switch ( fMethod ){
 138     case kLikeSign :
 139     case kLikeSignArithm :
 140     case kLikeSignRcorr:
 141     case kLikeSignArithmRcorr:
 142       ProcessLS(arrhist);    // process like-sign subtraction method
 143       break;
 144
 145     case kEventMixing :
 146       ProcessEM(arrhist);    // process event mixing method
 147       break;
 148
 149   case kRotation:
 150       ProcessRotation(arrhist);
 151       break;
 152
 153     default :
 154       AliWarning("Subtraction method not supported. Please check SetMethod() function.");
 155   }
 156 }
 157
 158 //______________________________________________
 159 void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
 160 {
 161   //
 162   // signal subtraction
 163   //
 164   fHistDataPP = (TH1*)(arrhist->At(AliDielectron::kEv1PP))->Clone("histPP");  // ++    SE
 165   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 166   fHistDataMM = (TH1*)(arrhist->At(AliDielectron::kEv1MM))->Clone("histMM");  // --    SE
 167   fHistDataPP->Sumw2();
 168   fHistDataPM->Sumw2();
 169   fHistDataMM->Sumw2();
 170   fHistDataPP->SetDirectory(0);
 171   fHistDataPM->SetDirectory(0);
 172   fHistDataMM->SetDirectory(0);
 173
 174   // rebin the histograms
 175   if (fRebin>1) {
 176     fHistDataPP->Rebin(fRebin);
 177     fHistDataPM->Rebin(fRebin);
 178     fHistDataMM->Rebin(fRebin);
 179   }
 180
 181   fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
 182                           fHistDataPM->GetXaxis()->GetNbins(),
 183                           fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 184   fHistRfactor->Sumw2();
 185   fHistRfactor->SetDirectory(0);
 186
 187   fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
 188                          fHistDataPM->GetXaxis()->GetNbins(),
 189                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 190   fHistSignal->SetDirectory(0);
 191   fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
 192                              fHistDataPM->GetXaxis()->GetNbins(),
 193                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 194   fHistBackground->SetDirectory(0);
 195
 196   // fill out background histogram
 197   for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
 198     Float_t pp = fHistDataPP->GetBinContent(ibin);
 199     Float_t mm = fHistDataMM->GetBinContent(ibin);
 200
 201     Float_t background = 2*TMath::Sqrt(pp*mm);
 202     Float_t ebackground = TMath::Sqrt(mm+pp);
 203     if (fMethod==kLikeSignArithm || fMethod==kLikeSignArithmRcorr ){
 204       //Arithmetic mean instead of geometric
 205       background=(pp+mm);
 206       ebackground=TMath::Sqrt(pp+mm);
 207       if (TMath::Abs(ebackground)<1e-30) ebackground=1;
 208     }
 209
 210     fHistBackground->SetBinContent(ibin, background);
 211     fHistBackground->SetBinError(ibin, ebackground);
 212   }
 213
 214   //correct LS spectrum bin-by-bin with R factor obtained in mixed events
 215   if(fMixingCorr || fMethod==kLikeSignRcorr || fMethod==kLikeSignArithmRcorr) {
 216
 217     TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP");  // ++    ME
 218     TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM");  // --    ME
 219
 220     TH1* histMixPM = 0x0;
 221     if (arrhist->At(AliDielectron::kEv1MEv2P)){
 222       histMixPM   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM");  // -+    ME
 223     }
 224
 225     if (arrhist->At(AliDielectron::kEv1PEv2M)){
 226       TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
 227       if (!histMixPM) fHistDataME   = (TH1*)htmp->Clone("mixPM");                   // +-    ME
 228       else histMixPM->Add(htmp);
 229     }
 230
 231     if (!histMixPM){
 232       AliError("For R-factor correction the mixed event histograms are requires. No +- histogram found");
 233       return;
 234     }
 235     histMixPM->Sumw2();
 236
 237     // rebin the histograms
 238     if (fRebin>1) {
 239       histMixPP->Rebin(fRebin);
 240       histMixMM->Rebin(fRebin);
 241       histMixPM->Rebin(fRebin);
 242     }
 243
 244     // fill out rfactor histogram
 245     for(Int_t ibin=1; ibin<=histMixPM->GetXaxis()->GetNbins(); ibin++) {
 246       Float_t pp  = histMixPP->GetBinContent(ibin);
 247       Float_t ppe = histMixPP->GetBinError(ibin);
 248       Float_t mm  = histMixMM->GetBinContent(ibin);
 249       Float_t mme = histMixMM->GetBinError(ibin);
 250       Float_t pm  = histMixPM->GetBinContent(ibin);
 251       Float_t pme = histMixPM->GetBinError(ibin);
 252
 253       Float_t background = 2*TMath::Sqrt(pp*mm);
 254       // do not use it since ME could be weighted err!=sqrt(entries)
 255       //      Float_t ebackground = TMath::Sqrt(mm+pp);
 256       Float_t ebackground = TMath::Sqrt(mm/pp*ppe*ppe + pp/mm*mme*mme);
 257       if (fMethod==kLikeSignArithm){
 258         //Arithmetic mean instead of geometric
 259         background=(pp+mm);
 260         ebackground=TMath::Sqrt(ppe*ppe+mme*mme);
 261         if (TMath::Abs(ebackground)<1e-30) ebackground=1;
 262       }
 263
 264       Float_t rcon = 1.0;
 265       Float_t rerr = 0.0;
 266       if(background>0.0) {
 267         rcon = pm/background;
 268         rerr = TMath::Sqrt((1./background)*(1./background) * pme*pme +
 269                            (pm/(background*background))*(pm/(background*background)) * ebackground*ebackground);
 270       }
 271       fHistRfactor->SetBinContent(ibin, rcon);
 272       fHistRfactor->SetBinError(ibin, rerr);
 273     }
 274
 275     fHistBackground->Multiply(fHistRfactor);
 276
 277     if (histMixPP) delete histMixPP;
 278     if (histMixMM) delete histMixMM;
 279     if (histMixPM) delete histMixPM;
 280   }
 281
 282   //scale histograms to match integral between fScaleMin and fScaleMax
 283   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 284   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 285   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 286   else if (fScaleMin>0.){
 287     fScaleFactor=fScaleMin;
 288     fHistBackground->Scale(fScaleFactor);
 289   }
 290
 291   //subract background
 292   fHistSignal->Add(fHistDataPM);
 293   fHistSignal->Add(fHistBackground,-1);
 294
 295   // background
 296   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 297                                                  fHistBackground->FindBin(fIntMax),
 298                                                  fErrors(1));
 299
 300   // signal depending on peak description method
 301   DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
 302   //printf("%f  %f\n",fValues(0),fValues(1));
 303   // S/B and significance
 304   //  SetSignificanceAndSOB();
 305
 306   fProcessed = kTRUE;
 307 }
 308
 309 //______________________________________________
 310 void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
 311 {
 312   //
 313   // event mixing of +- and -+
 314   //
 315
 316   if (!arrhist->At(AliDielectron::kEv1PM) || !(arrhist->At(AliDielectron::kEv1MEv2P) || arrhist->At(AliDielectron::kEv1PEv2M)) ){
 317     AliError("Either OS or mixed histogram missing");
 318     return;
 319   }
 320
 321   delete fHistDataPM; fHistDataPM=0x0;
 322   delete fHistDataME; fHistDataME=0x0;
 323   delete fHistBackground; fHistBackground=0x0;
 324
 325   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 326   fHistDataPM->Sumw2();
 327   fHistDataPM->SetDirectory(0x0);
 328
 329   if (arrhist->At(AliDielectron::kEv1MEv2P)){
 330     fHistDataME   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME");  // -+    ME
 331   }
 332
 333   if (arrhist->At(AliDielectron::kEv1PEv2M)){
 334     TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
 335     if (!fHistDataME) fHistDataME   = (TH1*)htmp->Clone("histMPME");  // -+    ME
 336     else fHistDataME->Add(htmp);
 337   }
 338
 339   fHistBackground = (TH1*)fHistDataME->Clone("ME_Background");
 340   fHistBackground->SetDirectory(0x0);
 341   fHistBackground->Sumw2();
 342
 343   // rebin the histograms
 344   if (fRebin>1) {
 345     fHistDataPM->Rebin(fRebin);
 346     fHistDataME->Rebin(fRebin);
 347     fHistBackground->Rebin(fRebin);
 348   }
 349
 350   //scale histograms to match integral between fScaleMin and fScaleMax
 351   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 352   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 353   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 354   else if (fScaleMin>0.){
 355     fScaleFactor=fScaleMin;
 356     fHistBackground->Scale(fScaleFactor);
 357   }
 358   fHistSignal=(TH1*)fHistDataPM->Clone("Signal");
 359   fHistSignal->Sumw2();
 360   //  printf(" err: %f %f \n",fHistSignal->GetBinError(75),TMath::Sqrt(fHistSignal->GetBinContent(75)));
 361   fHistSignal->Add(fHistBackground,-1.);
 362   //  printf(" err: %f %f \n",fHistSignal->GetBinError(75),TMath::Sqrt(fHistSignal->GetBinContent(75)));
 363 //     // signal
 364 //   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 365 //                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 366   // background
 367   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 368                                                  fHistBackground->FindBin(fIntMax),
 369                                                  fErrors(1));
 370
 371   // signal depending on peak description method
 372   DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
 373
 374   fProcessed = kTRUE;
 375 }
 376
 377 //______________________________________________
 378 void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
 379 {
 380   //
 381   // signal subtraction
 382   //
 383
 384   if (!arrhist->At(AliDielectron::kEv1PM) || !arrhist->At(AliDielectron::kEv1PMRot) ){
 385     AliError("Either OS or rotation histogram missing");
 386     return;
 387   }
 388
 389   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 390   fHistDataPM->Sumw2();
 391   fHistDataPM->SetDirectory(0x0);
 392
 393   fHistBackground = (TH1*)(arrhist->At(AliDielectron::kEv1PMRot))->Clone("histRotation");
 394   fHistBackground->Sumw2();
 395   fHistBackground->SetDirectory(0x0);
 396
 397   // rebin the histograms
 398   if (fRebin>1) {
 399     fHistDataPM->Rebin(fRebin);
 400     fHistBackground->Rebin(fRebin);
 401   }
 402
 403   //scale histograms to match integral between fScaleMin and fScaleMax
 404   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 405   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 406   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 407   else if (fScaleMin>0.){
 408     fScaleFactor=fScaleMin;
 409     fHistBackground->Scale(fScaleFactor);
 410   }
 411
 412   fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
 413   fHistSignal->Add(fHistBackground,-1.);
 414   fHistSignal->SetDirectory(0x0);
 415
 416   //     // signal
 417   //   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 418   //                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 419   // background
 420   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 421                                                  fHistBackground->FindBin(fIntMax),
 422                                                  fErrors(1));
 423   // signal depending on peak description method
 424   DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
 425
 426   fProcessed = kTRUE;
 427
 428 }
 429
 430 //______________________________________________
 431 void AliDielectronSignalExt::Draw(const Option_t* option)
 432 {
 433   //
 434   // Draw the fitted function
 435   //
 436   TString drawOpt(option);
 437   drawOpt.ToLower();
 438
 439   Float_t minY = 0.001;
 440   Float_t maxY = 1.2*fHistDataPM->GetMaximum();
 441   Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
 442   Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
 443   Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
 444   Float_t minMinY = fHistSignal->GetMinimum();
 445
 446   TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
 447   cSub->SetLeftMargin(0.15);
 448   cSub->SetRightMargin(0.0);
 449   cSub->SetTopMargin(0.002);
 450   cSub->SetBottomMargin(0.0);
 451   cSub->Divide(2,2,0.,0.);
 452   cSub->Draw();
 453
 454   TVirtualPad* pad = cSub->cd(1);
 455   pad->SetLeftMargin(0.15);
 456   pad->SetRightMargin(0.0);
 457   pad->SetTopMargin(0.005);
 458   pad->SetBottomMargin(0.0);
 459   TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
 460   range1->SetStats(kFALSE);
 461   range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 462   range1->GetYaxis()->CenterTitle();
 463   range1->GetYaxis()->SetLabelSize(0.05);
 464   range1->GetYaxis()->SetTitleSize(0.06);
 465   range1->GetYaxis()->SetTitleOffset(0.8);
 466   range1->Draw();
 467   fHistDataPM->SetLineColor(1);
 468   fHistDataPM->SetLineWidth(2);
 469   //  fHistDataPM->SetMarkerStyle(21);
 470   fHistDataPM->Draw("Psame");
 471   TLatex *latex = new TLatex();
 472   latex->SetNDC();
 473   latex->SetTextSize(0.05);
 474   latex->DrawLatex(0.2, 0.95, "Background un-substracted");
 475   TLine line;
 476   line.SetLineWidth(1);
 477   line.SetLineStyle(2);
 478   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 479   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 480
 481   pad = cSub->cd(2);
 482   pad->SetLeftMargin(0.);
 483   pad->SetRightMargin(0.005);
 484   pad->SetTopMargin(0.005);
 485   pad->SetBottomMargin(0.0);
 486   TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
 487   range2->SetStats(kFALSE);
 488   range2->Draw();
 489   fHistBackground->SetLineColor(4);
 490   fHistBackground->SetLineWidth(2);
 491   //  fHistBackground->SetMarkerColor(4);
 492   //  fHistBackground->SetMarkerStyle(6);
 493   fHistBackground->Draw("Psame");
 494   latex->DrawLatex(0.05, 0.95, "Like-sign background");
 495   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 496   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 497   TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
 498   legend->SetFillColor(0);
 499   legend->SetMargin(0.15);
 500   legend->AddEntry(fHistDataPM, "N_{+-}", "l");
 501   legend->AddEntry(fHistDataPP, "N_{++}", "l");
 502   legend->AddEntry(fHistDataMM, "N_{--}", "l");
 503   legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
 504   legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
 505   legend->Draw();
 506
 507
 508   pad = cSub->cd(3);
 509   pad->SetLeftMargin(0.15);
 510   pad->SetRightMargin(0.0);
 511   pad->SetTopMargin(0.0);
 512   pad->SetBottomMargin(0.15);
 513   TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
 514   range3->SetStats(kFALSE);
 515   range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 516   range3->GetYaxis()->CenterTitle();
 517   range3->GetYaxis()->SetLabelSize(0.05);
 518   range3->GetYaxis()->SetTitleSize(0.06);
 519   range3->GetYaxis()->SetTitleOffset(0.8);
 520   range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 521   range3->GetXaxis()->CenterTitle();
 522   range3->GetXaxis()->SetLabelSize(0.05);
 523   range3->GetXaxis()->SetTitleSize(0.06);
 524   range3->GetXaxis()->SetTitleOffset(1.0);
 525   range3->Draw();
 526   fHistDataPM->Draw("Psame");
 527   fHistDataPP->SetLineWidth(2);
 528   fHistDataPP->SetLineColor(6);
 529   fHistDataMM->SetLineWidth(2);
 530   fHistDataMM->SetLineColor(8);
 531   fHistDataPP->Draw("Psame");
 532   fHistDataMM->Draw("Psame");
 533   line.DrawLine(minX, 0.,maxX, 0.);
 534   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 535   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 536
 537   pad = cSub->cd(4);
 538   pad->SetLeftMargin(0.0);
 539   pad->SetRightMargin(0.005);
 540   pad->SetTopMargin(0.0);
 541   pad->SetBottomMargin(0.15);
 542   TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
 543   range4->SetStats(kFALSE);
 544   range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 545   range4->GetXaxis()->CenterTitle();
 546   range4->GetXaxis()->SetLabelSize(0.05);
 547   range4->GetXaxis()->SetTitleSize(0.06);
 548   range4->GetXaxis()->SetTitleOffset(1.0);
 549   range4->Draw();
 550   fHistSignal->SetLineWidth(2);
 551   fHistSignal->SetLineColor(2);
 552   fHistSignal->Draw("Psame");
 553   latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
 554   if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
 555   line.DrawLine(minX, 0.,maxX, 0.);
 556   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 557   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 558
 559   cSub->SaveAs(Form("%s_summary.png", fName.Data()));
 560 }
 561