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     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
 201   //scale histograms to match integral between fScaleMin and fScaleMax
 202   if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 203
 204   fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
 205   fHistSignal->Add(fHistBackground,-1.);
 206
 207     // signal
 208   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 209                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 210   // background
 211   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 212                                                  fHistBackground->FindBin(fIntMax),
 213                                                  fErrors(1));
 214   // S/B and significance
 215   SetSignificanceAndSOB();
 216
 217   fProcessed = kTRUE;
 218
 219 }
 220
 221 //______________________________________________
 222 void AliDielectronSignalExt::Draw(const Option_t* option)
 223 {
 224   //
 225   // Draw the fitted function
 226   //
 227   TString drawOpt(option);
 228   drawOpt.ToLower();
 229
 230   Float_t minY = 0.001;
 231   Float_t maxY = 1.2*fHistDataPM->GetMaximum();
 232   Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
 233   Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
 234   Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
 235   Float_t minMinY = fHistSignal->GetMinimum();
 236
 237   TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
 238   cSub->SetLeftMargin(0.15);
 239   cSub->SetRightMargin(0.0);
 240   cSub->SetTopMargin(0.002);
 241   cSub->SetBottomMargin(0.0);
 242   cSub->Divide(2,2,0.,0.);
 243   cSub->Draw();
 244
 245   TVirtualPad* pad = cSub->cd(1);
 246   pad->SetLeftMargin(0.15);
 247   pad->SetRightMargin(0.0);
 248   pad->SetTopMargin(0.005);
 249   pad->SetBottomMargin(0.0);
 250   TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
 251   range1->SetStats(kFALSE);
 252   range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 253   range1->GetYaxis()->CenterTitle();
 254   range1->GetYaxis()->SetLabelSize(0.05);
 255   range1->GetYaxis()->SetTitleSize(0.06);
 256   range1->GetYaxis()->SetTitleOffset(0.8);
 257   range1->Draw();
 258   fHistDataPM->SetLineColor(1);
 259   fHistDataPM->SetLineWidth(2);
 260   //  fHistDataPM->SetMarkerStyle(21);
 261   fHistDataPM->Draw("Psame");
 262   TLatex *latex = new TLatex();
 263   latex->SetNDC();
 264   latex->SetTextSize(0.05);
 265   latex->DrawLatex(0.2, 0.95, "Background un-substracted");
 266   TLine line;
 267   line.SetLineWidth(1);
 268   line.SetLineStyle(2);
 269   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 270   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 271
 272   pad = cSub->cd(2);
 273   pad->SetLeftMargin(0.);
 274   pad->SetRightMargin(0.005);
 275   pad->SetTopMargin(0.005);
 276   pad->SetBottomMargin(0.0);
 277   TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
 278   range2->SetStats(kFALSE);
 279   range2->Draw();
 280   fHistBackground->SetLineColor(4);
 281   fHistBackground->SetLineWidth(2);
 282   //  fHistBackground->SetMarkerColor(4);
 283   //  fHistBackground->SetMarkerStyle(6);
 284   fHistBackground->Draw("Psame");
 285   latex->DrawLatex(0.05, 0.95, "Like-sign background");
 286   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 287   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 288   TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
 289   legend->SetFillColor(0);
 290   legend->SetMargin(0.15);
 291   legend->AddEntry(fHistDataPM, "N_{+-}", "l");
 292   legend->AddEntry(fHistDataPP, "N_{++}", "l");
 293   legend->AddEntry(fHistDataMM, "N_{--}", "l");
 294   legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
 295   legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
 296   legend->Draw();
 297
 298
 299   pad = cSub->cd(3);
 300   pad->SetLeftMargin(0.15);
 301   pad->SetRightMargin(0.0);
 302   pad->SetTopMargin(0.0);
 303   pad->SetBottomMargin(0.15);
 304   TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
 305   range3->SetStats(kFALSE);
 306   range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 307   range3->GetYaxis()->CenterTitle();
 308   range3->GetYaxis()->SetLabelSize(0.05);
 309   range3->GetYaxis()->SetTitleSize(0.06);
 310   range3->GetYaxis()->SetTitleOffset(0.8);
 311   range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 312   range3->GetXaxis()->CenterTitle();
 313   range3->GetXaxis()->SetLabelSize(0.05);
 314   range3->GetXaxis()->SetTitleSize(0.06);
 315   range3->GetXaxis()->SetTitleOffset(1.0);
 316   range3->Draw();
 317   fHistDataPM->Draw("Psame");
 318   fHistDataPP->SetLineWidth(2);
 319   fHistDataPP->SetLineColor(6);
 320   fHistDataMM->SetLineWidth(2);
 321   fHistDataMM->SetLineColor(8);
 322   fHistDataPP->Draw("Psame");
 323   fHistDataMM->Draw("Psame");
 324   line.DrawLine(minX, 0.,maxX, 0.);
 325   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 326   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 327
 328   pad = cSub->cd(4);
 329   pad->SetLeftMargin(0.0);
 330   pad->SetRightMargin(0.005);
 331   pad->SetTopMargin(0.0);
 332   pad->SetBottomMargin(0.15);
 333   TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
 334   range4->SetStats(kFALSE);
 335   range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 336   range4->GetXaxis()->CenterTitle();
 337   range4->GetXaxis()->SetLabelSize(0.05);
 338   range4->GetXaxis()->SetTitleSize(0.06);
 339   range4->GetXaxis()->SetTitleOffset(1.0);
 340   range4->Draw();
 341   fHistSignal->SetLineWidth(2);
 342   fHistSignal->SetLineColor(2);
 343   fHistSignal->Draw("Psame");
 344   latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
 345   if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
 346   line.DrawLine(minX, 0.,maxX, 0.);
 347   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 348   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 349
 350   cSub->SaveAs(Form("%s_summary.png", fName.Data()));
 351 }
 352