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 #include "AliDielectron.h"
  44
  45 ClassImp(AliDielectronSignalExt)
  46
  47 AliDielectronSignalExt::AliDielectronSignalExt() :
  48   AliDielectronSignalBase()
  49 {
  50   //
  51   // Default Constructor
  52   //
  53 }
  54
  55 //______________________________________________
  56 AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
  57   AliDielectronSignalBase(name, title)
  58 {
  59   //
  60   // Named Constructor
  61   //
  62 }
  63
  64 //______________________________________________
  65 AliDielectronSignalExt::~AliDielectronSignalExt()
  66 {
  67   //
  68   // Default Destructor
  69   //
  70 }
  71
  72 //______________________________________________
  73 void AliDielectronSignalExt::Process(TObjArray* const arrhist)
  74 {
  75   //
  76   // signal subtraction. support like-sign subtraction and event mixing method
  77   //
  78   switch ( fMethod ){
  79     case kLikeSign :
  80     case kLikeSignArithm :
  81     case kLikeSignRcorr:
  82     case kLikeSignArithmRcorr:
  83       ProcessLS(arrhist);    // process like-sign subtraction method
  84       break;
  85
  86     case kEventMixing :
  87       ProcessEM(arrhist);    // process event mixing method
  88       break;
  89
  90   case kRotation:
  91       ProcessRotation(arrhist);
  92       break;
  93
  94     default :
  95       AliWarning("Subtraction method not supported. Please check SetMethod() function.");
  96   }
  97 }
  98
  99 //______________________________________________
 100 void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
 101 {
 102   //
 103   // signal subtraction
 104   //
 105   fHistDataPP = (TH1*)(arrhist->At(AliDielectron::kEv1PP))->Clone("histPP");  // ++    SE
 106   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 107   fHistDataMM = (TH1*)(arrhist->At(AliDielectron::kEv1MM))->Clone("histMM");  // --    SE
 108   fHistDataPP->Sumw2();
 109   fHistDataPM->Sumw2();
 110   fHistDataMM->Sumw2();
 111   fHistDataPP->SetDirectory(0);
 112   fHistDataPM->SetDirectory(0);
 113   fHistDataMM->SetDirectory(0);
 114
 115   // rebin the histograms
 116   if (fRebin>1) {
 117     fHistDataPP->Rebin(fRebin);
 118     fHistDataPM->Rebin(fRebin);
 119     fHistDataMM->Rebin(fRebin);
 120   }
 121
 122   fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
 123                           fHistDataPM->GetXaxis()->GetNbins(),
 124                           fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 125   fHistRfactor->Sumw2();
 126   fHistRfactor->SetDirectory(0);
 127
 128   fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
 129                          fHistDataPM->GetXaxis()->GetNbins(),
 130                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 131   fHistSignal->SetDirectory(0);
 132   fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
 133                              fHistDataPM->GetXaxis()->GetNbins(),
 134                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 135   fHistBackground->SetDirectory(0);
 136
 137   // fill out background histogram
 138   for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
 139     Float_t pp = fHistDataPP->GetBinContent(ibin);
 140     Float_t mm = fHistDataMM->GetBinContent(ibin);
 141
 142     Float_t background = 2*TMath::Sqrt(pp*mm);
 143     Float_t ebackground = TMath::Sqrt(mm+pp);
 144     if (fMethod==kLikeSignArithm || fMethod==kLikeSignArithmRcorr ){
 145       //Arithmetic mean instead of geometric
 146       background=(pp+mm);
 147       ebackground=TMath::Sqrt(pp+mm);
 148       if (TMath::Abs(ebackground)<1e-30) ebackground=1;
 149     }
 150
 151     fHistBackground->SetBinContent(ibin, background);
 152     fHistBackground->SetBinError(ibin, ebackground);
 153   }
 154
 155   //correct LS spectrum bin-by-bin with R factor obtained in mixed events
 156   if(fMixingCorr || fMethod==kLikeSignRcorr || fMethod==kLikeSignArithmRcorr) {
 157
 158     TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP");  // ++    ME
 159     TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM");  // --    ME
 160
 161     TH1* histMixPM = 0x0;
 162     if (arrhist->At(AliDielectron::kEv1MEv2P)){
 163       histMixPM   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM");  // -+    ME
 164     }
 165
 166     if (arrhist->At(AliDielectron::kEv1PEv2M)){
 167       TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
 168       if (!histMixPM) fHistDataME   = (TH1*)htmp->Clone("mixPM");                   // +-    ME
 169       else histMixPM->Add(htmp);
 170     }
 171
 172     if (!histMixPM){
 173       AliError("For R-factor correction the mixed event histograms are requires. No +- histogram found");
 174       return;
 175     }
 176     histMixPM->Sumw2();
 177
 178     // rebin the histograms
 179     if (fRebin>1) {
 180       histMixPP->Rebin(fRebin);
 181       histMixMM->Rebin(fRebin);
 182       histMixPM->Rebin(fRebin);
 183     }
 184
 185     // fill out rfactor histogram
 186     for(Int_t ibin=1; ibin<=histMixPM->GetXaxis()->GetNbins(); ibin++) {
 187       Float_t pp  = histMixPP->GetBinContent(ibin);
 188       Float_t ppe = histMixPP->GetBinError(ibin);
 189       Float_t mm  = histMixMM->GetBinContent(ibin);
 190       Float_t mme = histMixMM->GetBinError(ibin);
 191       Float_t pm  = histMixPM->GetBinContent(ibin);
 192       Float_t pme = histMixPM->GetBinError(ibin);
 193
 194       Float_t background = 2*TMath::Sqrt(pp*mm);
 195       // do not use it since ME could be weighted err!=sqrt(entries)
 196       //      Float_t ebackground = TMath::Sqrt(mm+pp);
 197       Float_t ebackground = TMath::Sqrt(mm/pp*ppe*ppe + pp/mm*mme*mme);
 198       if (fMethod==kLikeSignArithm){
 199         //Arithmetic mean instead of geometric
 200         background=(pp+mm);
 201         ebackground=TMath::Sqrt(ppe*ppe+mme*mme);
 202         if (TMath::Abs(ebackground)<1e-30) ebackground=1;
 203       }
 204
 205       Float_t rcon = 1.0;
 206       Float_t rerr = 0.0;
 207       if(background>0.0) {
 208         rcon = pm/background;
 209         rerr = TMath::Sqrt((1./background)*(1./background) * pme*pme +
 210                            (pm/(background*background))*(pm/(background*background)) * ebackground*ebackground);
 211       }
 212       fHistRfactor->SetBinContent(ibin, rcon);
 213       fHistRfactor->SetBinError(ibin, rerr);
 214     }
 215
 216     fHistBackground->Multiply(fHistRfactor);
 217
 218     if (histMixPP) delete histMixPP;
 219     if (histMixMM) delete histMixMM;
 220     if (histMixPM) delete histMixPM;
 221   }
 222
 223   //scale histograms to match integral between fScaleMin and fScaleMax
 224   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 225   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 226   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 227   else if (fScaleMin>0.){
 228     fScaleFactor=fScaleMin;
 229     fHistBackground->Scale(fScaleFactor);
 230   }
 231
 232   //subract background
 233   fHistSignal->Add(fHistDataPM);
 234   fHistSignal->Add(fHistBackground,-1);
 235
 236   // signal
 237   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 238                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 239   // background
 240   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 241                                                  fHistBackground->FindBin(fIntMax),
 242                                                  fErrors(1));
 243   //printf("%f  %f\n",fValues(0),fValues(1));
 244   // S/B and significance
 245   SetSignificanceAndSOB();
 246
 247   fProcessed = kTRUE;
 248 }
 249
 250 //______________________________________________
 251 void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
 252 {
 253   //
 254   // event mixing of +- and -+
 255   //
 256
 257   if (!arrhist->At(AliDielectron::kEv1PM) || !(arrhist->At(AliDielectron::kEv1MEv2P) || arrhist->At(AliDielectron::kEv1PEv2M)) ){
 258     AliError("Either OS or mixed histogram missing");
 259     return;
 260   }
 261
 262   delete fHistDataPM; fHistDataPM=0x0;
 263   delete fHistDataME; fHistDataME=0x0;
 264   delete fHistBackground; fHistBackground=0x0;
 265
 266   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 267   fHistDataPM->Sumw2();
 268   fHistDataPM->SetDirectory(0x0);
 269
 270   if (arrhist->At(AliDielectron::kEv1MEv2P)){
 271     fHistDataME   = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME");  // -+    ME
 272   }
 273
 274   if (arrhist->At(AliDielectron::kEv1PEv2M)){
 275     TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
 276     if (!fHistDataME) fHistDataME   = (TH1*)htmp->Clone("histMPME");  // -+    ME
 277     else fHistDataME->Add(htmp);
 278   }
 279
 280   fHistBackground = (TH1*)fHistDataME->Clone("ME_Background");
 281   fHistBackground->SetDirectory(0x0);
 282   fHistBackground->Sumw2();
 283
 284   // rebin the histograms
 285   if (fRebin>1) {
 286     fHistDataPM->Rebin(fRebin);
 287     fHistDataME->Rebin(fRebin);
 288     fHistBackground->Rebin(fRebin);
 289   }
 290
 291   //scale histograms to match integral between fScaleMin and fScaleMax
 292   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 293   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 294   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 295   else if (fScaleMin>0.){
 296     fScaleFactor=fScaleMin;
 297     fHistBackground->Scale(fScaleFactor);
 298   }
 299
 300   fHistSignal=(TH1*)fHistDataPM->Clone("Signal");
 301   fHistSignal->Add(fHistBackground,-1.);
 302
 303     // signal
 304   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 305                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 306   // background
 307   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 308                                                  fHistBackground->FindBin(fIntMax),
 309                                                  fErrors(1));
 310   // S/B and significance
 311   SetSignificanceAndSOB();
 312
 313   fProcessed = kTRUE;
 314 }
 315
 316 //______________________________________________
 317 void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
 318 {
 319   //
 320   // signal subtraction
 321   //
 322
 323   if (!arrhist->At(AliDielectron::kEv1PM) || !arrhist->At(AliDielectron::kEv1PMRot) ){
 324     AliError("Either OS or rotation histogram missing");
 325     return;
 326   }
 327
 328   fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM");  // +-    SE
 329   fHistDataPM->Sumw2();
 330   fHistDataPM->SetDirectory(0x0);
 331
 332   fHistBackground = (TH1*)(arrhist->At(AliDielectron::kEv1PMRot))->Clone("histRotation");
 333   fHistBackground->Sumw2();
 334   fHistBackground->SetDirectory(0x0);
 335
 336   // rebin the histograms
 337   if (fRebin>1) {
 338     fHistDataPM->Rebin(fRebin);
 339     fHistBackground->Rebin(fRebin);
 340   }
 341
 342   //scale histograms to match integral between fScaleMin and fScaleMax
 343   // or if fScaleMax <  fScaleMin use fScaleMin as scale factor
 344   if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
 345   else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
 346   else if (fScaleMin>0.){
 347     fScaleFactor=fScaleMin;
 348     fHistBackground->Scale(fScaleFactor);
 349   }
 350
 351   fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
 352   fHistSignal->Add(fHistBackground,-1.);
 353   fHistSignal->SetDirectory(0x0);
 354
 355     // signal
 356   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 357                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 358   // background
 359   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 360                                                  fHistBackground->FindBin(fIntMax),
 361                                                  fErrors(1));
 362   // S/B and significance
 363   SetSignificanceAndSOB();
 364
 365   fProcessed = kTRUE;
 366
 367 }
 368
 369 //______________________________________________
 370 void AliDielectronSignalExt::Draw(const Option_t* option)
 371 {
 372   //
 373   // Draw the fitted function
 374   //
 375   TString drawOpt(option);
 376   drawOpt.ToLower();
 377
 378   Float_t minY = 0.001;
 379   Float_t maxY = 1.2*fHistDataPM->GetMaximum();
 380   Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
 381   Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
 382   Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1));   // in MeV
 383   Float_t minMinY = fHistSignal->GetMinimum();
 384
 385   TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
 386   cSub->SetLeftMargin(0.15);
 387   cSub->SetRightMargin(0.0);
 388   cSub->SetTopMargin(0.002);
 389   cSub->SetBottomMargin(0.0);
 390   cSub->Divide(2,2,0.,0.);
 391   cSub->Draw();
 392
 393   TVirtualPad* pad = cSub->cd(1);
 394   pad->SetLeftMargin(0.15);
 395   pad->SetRightMargin(0.0);
 396   pad->SetTopMargin(0.005);
 397   pad->SetBottomMargin(0.0);
 398   TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
 399   range1->SetStats(kFALSE);
 400   range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 401   range1->GetYaxis()->CenterTitle();
 402   range1->GetYaxis()->SetLabelSize(0.05);
 403   range1->GetYaxis()->SetTitleSize(0.06);
 404   range1->GetYaxis()->SetTitleOffset(0.8);
 405   range1->Draw();
 406   fHistDataPM->SetLineColor(1);
 407   fHistDataPM->SetLineWidth(2);
 408   //  fHistDataPM->SetMarkerStyle(21);
 409   fHistDataPM->Draw("Psame");
 410   TLatex *latex = new TLatex();
 411   latex->SetNDC();
 412   latex->SetTextSize(0.05);
 413   latex->DrawLatex(0.2, 0.95, "Background un-substracted");
 414   TLine line;
 415   line.SetLineWidth(1);
 416   line.SetLineStyle(2);
 417   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 418   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 419
 420   pad = cSub->cd(2);
 421   pad->SetLeftMargin(0.);
 422   pad->SetRightMargin(0.005);
 423   pad->SetTopMargin(0.005);
 424   pad->SetBottomMargin(0.0);
 425   TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
 426   range2->SetStats(kFALSE);
 427   range2->Draw();
 428   fHistBackground->SetLineColor(4);
 429   fHistBackground->SetLineWidth(2);
 430   //  fHistBackground->SetMarkerColor(4);
 431   //  fHistBackground->SetMarkerStyle(6);
 432   fHistBackground->Draw("Psame");
 433   latex->DrawLatex(0.05, 0.95, "Like-sign background");
 434   line.DrawLine(fIntMin, minY, fIntMin, maxY);
 435   line.DrawLine(fIntMax, minY, fIntMax, maxY);
 436   TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
 437   legend->SetFillColor(0);
 438   legend->SetMargin(0.15);
 439   legend->AddEntry(fHistDataPM, "N_{+-}", "l");
 440   legend->AddEntry(fHistDataPP, "N_{++}", "l");
 441   legend->AddEntry(fHistDataMM, "N_{--}", "l");
 442   legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
 443   legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
 444   legend->Draw();
 445
 446
 447   pad = cSub->cd(3);
 448   pad->SetLeftMargin(0.15);
 449   pad->SetRightMargin(0.0);
 450   pad->SetTopMargin(0.0);
 451   pad->SetBottomMargin(0.15);
 452   TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
 453   range3->SetStats(kFALSE);
 454   range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
 455   range3->GetYaxis()->CenterTitle();
 456   range3->GetYaxis()->SetLabelSize(0.05);
 457   range3->GetYaxis()->SetTitleSize(0.06);
 458   range3->GetYaxis()->SetTitleOffset(0.8);
 459   range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 460   range3->GetXaxis()->CenterTitle();
 461   range3->GetXaxis()->SetLabelSize(0.05);
 462   range3->GetXaxis()->SetTitleSize(0.06);
 463   range3->GetXaxis()->SetTitleOffset(1.0);
 464   range3->Draw();
 465   fHistDataPM->Draw("Psame");
 466   fHistDataPP->SetLineWidth(2);
 467   fHistDataPP->SetLineColor(6);
 468   fHistDataMM->SetLineWidth(2);
 469   fHistDataMM->SetLineColor(8);
 470   fHistDataPP->Draw("Psame");
 471   fHistDataMM->Draw("Psame");
 472   line.DrawLine(minX, 0.,maxX, 0.);
 473   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 474   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 475
 476   pad = cSub->cd(4);
 477   pad->SetLeftMargin(0.0);
 478   pad->SetRightMargin(0.005);
 479   pad->SetTopMargin(0.0);
 480   pad->SetBottomMargin(0.15);
 481   TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
 482   range4->SetStats(kFALSE);
 483   range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
 484   range4->GetXaxis()->CenterTitle();
 485   range4->GetXaxis()->SetLabelSize(0.05);
 486   range4->GetXaxis()->SetTitleSize(0.06);
 487   range4->GetXaxis()->SetTitleOffset(1.0);
 488   range4->Draw();
 489   fHistSignal->SetLineWidth(2);
 490   fHistSignal->SetLineColor(2);
 491   fHistSignal->Draw("Psame");
 492   latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
 493   if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
 494   line.DrawLine(minX, 0.,maxX, 0.);
 495   line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
 496   line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
 497
 498   cSub->SaveAs(Form("%s_summary.png", fName.Data()));
 499 }
 500