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