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