PWG3/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 Dielectron signal extraction class using functions as input.
  22
  23 A function to describe the signal as well as one to describe the background
  24 has to be deployed by the user. Alternatively on of the default implementaions
  25 can be used.
  26
  27 */
  28 //                                                                       //
  29 ///////////////////////////////////////////////////////////////////////////
  30 #include <TF1.h>
  31 #include <TH1.h>
  32 #include <TH2F.h>
  33 #include <TLatex.h>
  34 #include <TLegend.h>
  35 #include <TCanvas.h>
  36 #include <TMath.h>
  37 #include <TString.h>
  38 #include <TLine.h>
  39
  40 #include <AliLog.h>
  41
  42 #include "AliDielectronSignalExt.h"
  43
  44 ClassImp(AliDielectronSignalExt)
  45
  46 AliDielectronSignalExt::AliDielectronSignalExt() :
  47   AliDielectronSignalBase()
  48 {
  49   //
  50   // Default Constructor
  51   //
  52 }
  53
  54 //______________________________________________
  55 AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
  56   AliDielectronSignalBase(name, title)
  57 {
  58   //
  59   // Named Constructor
  60   //
  61 }
  62
  63 //______________________________________________
  64 AliDielectronSignalExt::~AliDielectronSignalExt()
  65 {
  66   //
  67   // Default Destructor
  68   //
  69 }
  70
  71 //______________________________________________
  72 void AliDielectronSignalExt::Process(TObjArray* const arrhist)
  73 {
  74   //
  75   // signal subtraction. support like-sign subtraction and event mixing method
  76   //
  77   switch ( fMethod ){
  78     case kLikeSign :
  79       ProcessLS(arrhist);    // process like-sign subtraction method
  80       break;
  81
  82     case kEventMixing :
  83       ProcessEM(arrhist);    // process event mixing method
  84       break;
  85
  86   case kRotation:
  87       ProcessRotation(arrhist);
  88       break;
  89
  90     default :
  91       AliWarning("Subtraction method not supported. Please check SetMethod() function.");
  92   }
  93 }
  94
  95 //______________________________________________
  96 void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
  97 {
  98   //
  99   // signal subtraction
 100   //
 101   fHistDataPP = (TH1*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 102   fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 103   fHistDataMM = (TH1*)(arrhist->At(2))->Clone("histMM");  // --    SE
 104   fHistDataPP->Sumw2();
 105   fHistDataPM->Sumw2();
 106   fHistDataMM->Sumw2();
 107   fHistDataPP->SetDirectory(0);
 108   fHistDataPM->SetDirectory(0);
 109   fHistDataMM->SetDirectory(0);
 110
 111   // rebin the histograms
 112   if (fRebin>1) {
 113     fHistDataPP->Rebin(fRebin);
 114     fHistDataPM->Rebin(fRebin);
 115     fHistDataMM->Rebin(fRebin);
 116   }
 117
 118   fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
 119                          fHistDataPM->GetXaxis()->GetNbins(),
 120                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 121   fHistSignal->SetDirectory(0);
 122   fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
 123                              fHistDataPM->GetXaxis()->GetNbins(),
 124                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 125   fHistBackground->SetDirectory(0);
 126
 127   // fill out background and subtracted histogram
 128   for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
 129     Float_t pm = fHistDataPM->GetBinContent(ibin);
 130     Float_t pp = fHistDataPP->GetBinContent(ibin);
 131     Float_t mm = fHistDataMM->GetBinContent(ibin);
 132     Float_t epm = fHistDataPM->GetBinError(ibin);
 133
 134     Float_t background = 2*TMath::Sqrt(pp*mm);
 135     Float_t ebackground = TMath::Sqrt(mm+pp);
 136 //     Float_t signal = pm - background;
 137 //     Float_t error = TMath::Sqrt(epm*epm+mm+pp);
 138
 139     fHistSignal->SetBinContent(ibin, pm);
 140     fHistSignal->SetBinError(ibin, epm);
 141     fHistBackground->SetBinContent(ibin, background);
 142     fHistBackground->SetBinError(ibin, ebackground);
 143   }
 144   //scale background histogram to match integral of the data histogram between fScaleMin and fScaleMax
 145   if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 146
 147   //subract background
 148   fHistSignal->Add(fHistBackground,-1);
 149
 150   // signal
 151   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 152                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 153   // background
 154   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 155                                                  fHistBackground->FindBin(fIntMax),
 156                                                  fErrors(1));
 157   // S/B and significance
 158   SetSignificanceAndSOB();
 159
 160   fProcessed = kTRUE;
 161 }
 162
 163 //______________________________________________
 164 void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
 165 {
 166   //
 167   // event mixing. not implemented yet.
 168   //
 169   printf("event mixing method is not implemented yet. Like-sign method will be used.\n");
 170   ProcessLS(arrhist);
 171 }
 172
 173 //______________________________________________
 174 void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
 175 {
 176   //
 177   // signal subtraction
 178   //
 179   fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 180   if (!fHistDataPM){
 181     AliError("Unlike sign histogram not available. Cannot extract the signal.");
 182     return;
 183   }
 184   fHistDataPM->Sumw2();
 185
 186   fHistBackground=(TH1*)arrhist->At(10)->Clone("histRotation");
 187   if (!fHistBackground){
 188     AliError("Histgram from rotation not available. Cannot extract the signal.");
 189     delete fHistDataPM;
 190     fHistDataPM=0x0;
 191     return;
 192   }
 193
 194   //scale histograms to match integral between fScaleMin and fScaleMax
 195   if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 196
 197   fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
 198   fHistSignal->Add(fHistBackground,-1.);
 199
 200     // signal
 201   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 202                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 203   // background
 204   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 205                                                  fHistBackground->FindBin(fIntMax),
 206                                                  fErrors(1));
 207   // S/B and significance
 208   SetSignificanceAndSOB();
 209
 210   fProcessed = kTRUE;
 211
 212 }
 213
 214 //______________________________________________
 215 void AliDielectronSignalExt::Draw(const Option_t* option)
 216 {
 217   //
 218   // Draw the fitted function
 219   //
 220   TString drawOpt(option);
 221   drawOpt.ToLower();
 222
 223   Float_t minY = 0.001;
 224   Float_t maxY = 1.2*fHistDataPM->GetMaximum();
 225   Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
 226   Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
 227   Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
 228   Float_t minMinY = fHistSignal->GetMinimum();
 229
 230   TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
 231   cSub->SetLeftMargin(0.15);
 232   cSub->SetRightMargin(0.0);
 233   cSub->SetTopMargin(0.002);
 234   cSub->SetBottomMargin(0.0);
 235   cSub->Divide(2,2,0.,0.);
 236   cSub->Draw();
 237
 238   TVirtualPad* pad = cSub->cd(1);
 239   pad->SetLeftMargin(0.15);
 240   pad->SetRightMargin(0.0);
 241   pad->SetTopMargin(0.005);
 242   pad->SetBottomMargin(0.0);
 243   TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
 244   range1->SetStats(kFALSE);
 245   range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 246   range1->GetYaxis()->CenterTitle();
 247   range1->GetYaxis()->SetLabelSize(0.05);
 248   range1->GetYaxis()->SetTitleSize(0.06);
 249   range1->GetYaxis()->SetTitleOffset(0.8);
 250   range1->Draw();
 251   fHistDataPM->SetLineColor(1);
 252   fHistDataPM->SetLineWidth(2);
 253   //  fHistDataPM->SetMarkerStyle(21);
 254   fHistDataPM->Draw("Psame");
 255   TLatex *latex = new TLatex();
 256   latex->SetNDC();
 257   latex->SetTextSize(0.05);
 258   latex->DrawLatex(0.2, 0.95, "Background un-substracted");
 259   TLine line;
 260   line.SetLineWidth(1);
 261   line.SetLineStyle(2);
 262   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 263   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 264
 265   pad = cSub->cd(2);
 266   pad->SetLeftMargin(0.);
 267   pad->SetRightMargin(0.005);
 268   pad->SetTopMargin(0.005);
 269   pad->SetBottomMargin(0.0);
 270   TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
 271   range2->SetStats(kFALSE);
 272   range2->Draw();
 273   fHistBackground->SetLineColor(4);
 274   fHistBackground->SetLineWidth(2);
 275   //  fHistBackground->SetMarkerColor(4);
 276   //  fHistBackground->SetMarkerStyle(6);
 277   fHistBackground->Draw("Psame");
 278   latex->DrawLatex(0.05, 0.95, "Like-sign background");
 279   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 280   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 281   TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
 282   legend->SetFillColor(0);
 283   legend->SetMargin(0.15);
 284   legend->AddEntry(fHistDataPM, "N_{+-}", "l");
 285   legend->AddEntry(fHistDataPP, "N_{++}", "l");
 286   legend->AddEntry(fHistDataMM, "N_{--}", "l");
 287   legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
 288   legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
 289   legend->Draw();
 290
 291
 292   pad = cSub->cd(3);
 293   pad->SetLeftMargin(0.15);
 294   pad->SetRightMargin(0.0);
 295   pad->SetTopMargin(0.0);
 296   pad->SetBottomMargin(0.15);
 297   TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
 298   range3->SetStats(kFALSE);
 299   range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 300   range3->GetYaxis()->CenterTitle();
 301   range3->GetYaxis()->SetLabelSize(0.05);
 302   range3->GetYaxis()->SetTitleSize(0.06);
 303   range3->GetYaxis()->SetTitleOffset(0.8);
 304   range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 305   range3->GetXaxis()->CenterTitle();
 306   range3->GetXaxis()->SetLabelSize(0.05);
 307   range3->GetXaxis()->SetTitleSize(0.06);
 308   range3->GetXaxis()->SetTitleOffset(1.0);
 309   range3->Draw();
 310   fHistDataPM->Draw("Psame");
 311   fHistDataPP->SetLineWidth(2);
 312   fHistDataPP->SetLineColor(6);
 313   fHistDataMM->SetLineWidth(2);
 314   fHistDataMM->SetLineColor(8);
 315   fHistDataPP->Draw("Psame");
 316   fHistDataMM->Draw("Psame");
 317   line.DrawLine(minX, 0.,maxX, 0.);
 318   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 319   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 320
 321   pad = cSub->cd(4);
 322   pad->SetLeftMargin(0.0);
 323   pad->SetRightMargin(0.005);
 324   pad->SetTopMargin(0.0);
 325   pad->SetBottomMargin(0.15);
 326   TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
 327   range4->SetStats(kFALSE);
 328   range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 329   range4->GetXaxis()->CenterTitle();
 330   range4->GetXaxis()->SetLabelSize(0.05);
 331   range4->GetXaxis()->SetTitleSize(0.06);
 332   range4->GetXaxis()->SetTitleOffset(1.0);
 333   range4->Draw();
 334   fHistSignal->SetLineWidth(2);
 335   fHistSignal->SetLineColor(2);
 336   fHistSignal->Draw("Psame");
 337   latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
 338   if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
 339   line.DrawLine(minX, 0.,maxX, 0.);
 340   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 341   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 342
 343   cSub->SaveAs(Form("%s_summary.png", fName.Data()));
 344 }
 345