// Dielectron SignalExt //
// //
/*
-Dielectron signal extraction class using functions as input.
-A function to describe the signal as well as one to describe the background
-has to be deployed by the user. Alternatively on of the default implementaions
-can be used.
+ Class used for extracting the signal from an invariant mass spectrum.
+ Used invariant mass spectra are provided via an array of histograms. There are serveral method
+ to estimate the background and to extract the raw yield from the background subtracted spectra.
+
+ Example usage:
+
+ AliDielectronSignalExt *sig = new AliDielectronSignalExt();
+
+
+ 1) invariant mass input spectra
+
+ 1.1) Assuming a AliDielectronCF container as data format (check class for more details)
+ AliDielectronCFdraw *cf = new AliDielectronCFdraw("path/to/the/output/file.root");
+ TObjArray *arrHists = cf->CollectMinvProj(cf->FindStep("Config"));
+
+ 1.2) Assuming a AliDielectronHF grid as data format (check class for more details)
+ AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
+ TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
+
+ 1.3) Assuming a single histograms
+ TObjArray *histoArray = new TObjArray();
+ arrHists->Add(signalPP); // add the spectrum histograms to the array
+ arrHists->Add(signalPM); // the order is important !!!
+ arrHists->Add(signalMM);
+
+
+ 2) background estimation
+
+ 2.1) set the method for the background estimation (methods can be found in AliDielectronSignalBase)
+ sig->SetMethod(AliDielectronSignalBase::kEventMixing);
+ 2.2) rebin the spectras if needed
+ // sig->SetRebin(2);
+ 2.3) normalize the backgound spectum to the odd-sign spectrum in the desired range(s)
+ sig->SetScaleRawToBackground(minScale, maxScale);
+ // sig->SetScaleRawToBackground(minScale, maxScale, minScale2, maxScale2);
+
+
+ 3) configure the signal extraction
+
+ 3.1) set the method for the signal extraction (methods can be found in AliDielectronSignalBase)
+ depending on the method serveral inputs are needed (e.g. MC shape, PDG code of the particle of interest)
+ // sig->SetParticleOfInterest(443); //default is jpsi
+ // sig->SetMCSignalShape(signalMC);
+ sig->SetIntegralRange(minInt, maxInt); // range for bin counting
+ sig->SetExtractionMethod(AliDielectronSignal::BinCounting); // this is the default
+
+
+ 4) start the processing
+
+ sig->Process(arrHists);
+ sig->Print(""); // print values and errors extracted
+
+
+ 5) access the spectra and values created
+
+ 5.1) standard spectras
+ TH1F *hsign = (TH1F*) sig->GetUnlikeSignHistogram(); // same as the input (rebinned)
+ TH1F *hbgrd = (TH1F*) sig->GetBackgroundHistogram(); // scaled input (rebinned)
+ TH1F *hextr = (TH1F*) sig->GetSignalHistogram(); // after backgound extraction (rebinned)
+ TObject *oPeak = (TObject*) sig->GetPeakShape(); // can be a TF1 or TH1 depending on the extraction method
+ TH1F *hrfac = (TH1F*) sig->GetRfactorHistogram(); // if like-sign correction was activated, o.w. 0x0
+ 5.2) access the extracted values and errors
+ sig->GetValues(); or GetErrors(); // yield extraction
*/
// //
#include <AliLog.h>
#include "AliDielectronSignalExt.h"
+#include "AliDielectron.h"
ClassImp(AliDielectronSignalExt)
switch ( fMethod ){
case kLikeSign :
case kLikeSignArithm :
+ case kLikeSignRcorr:
+ case kLikeSignArithmRcorr:
ProcessLS(arrhist); // process like-sign subtraction method
break;
//
// signal subtraction
//
- fHistDataPP = (TH1*)(arrhist->At(0))->Clone("histPP"); // ++ SE
- fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM"); // +- SE
- fHistDataMM = (TH1*)(arrhist->At(2))->Clone("histMM"); // -- SE
+ fHistDataPP = (TH1*)(arrhist->At(AliDielectron::kEv1PP))->Clone("histPP"); // ++ SE
+ fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
+ fHistDataMM = (TH1*)(arrhist->At(AliDielectron::kEv1MM))->Clone("histMM"); // -- SE
fHistDataPP->Sumw2();
fHistDataPM->Sumw2();
fHistDataMM->Sumw2();
fHistDataMM->Rebin(fRebin);
}
+ fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
+ fHistDataPM->GetXaxis()->GetNbins(),
+ fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
+ fHistRfactor->Sumw2();
+ fHistRfactor->SetDirectory(0);
+
fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
fHistDataPM->GetXaxis()->GetNbins(),
fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
fHistBackground->SetDirectory(0);
- // fill out background and subtracted histogram
+ // fill out background histogram
for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
- Float_t pm = fHistDataPM->GetBinContent(ibin);
Float_t pp = fHistDataPP->GetBinContent(ibin);
Float_t mm = fHistDataMM->GetBinContent(ibin);
- Float_t epm = fHistDataPM->GetBinError(ibin);
Float_t background = 2*TMath::Sqrt(pp*mm);
Float_t ebackground = TMath::Sqrt(mm+pp);
- if (fMethod==kLikeSignArithm){
+ if (fMethod==kLikeSignArithm || fMethod==kLikeSignArithmRcorr ){
//Arithmetic mean instead of geometric
background=(pp+mm);
ebackground=TMath::Sqrt(pp+mm);
if (TMath::Abs(ebackground)<1e-30) ebackground=1;
}
-// Float_t signal = pm - background;
-// Float_t error = TMath::Sqrt(epm*epm+mm+pp);
- fHistSignal->SetBinContent(ibin, pm);
- fHistSignal->SetBinError(ibin, epm);
fHistBackground->SetBinContent(ibin, background);
fHistBackground->SetBinError(ibin, ebackground);
}
+
+ //correct LS spectrum bin-by-bin with R factor obtained in mixed events
+ if(fMixingCorr || fMethod==kLikeSignRcorr || fMethod==kLikeSignArithmRcorr) {
+
+ TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP"); // ++ ME
+ TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM"); // -- ME
+
+ TH1* histMixPM = 0x0;
+ if (arrhist->At(AliDielectron::kEv1MEv2P)){
+ histMixPM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM"); // -+ ME
+ }
+
+ if (arrhist->At(AliDielectron::kEv1PEv2M)){
+ TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
+ if (!histMixPM) fHistDataME = (TH1*)htmp->Clone("mixPM"); // +- ME
+ else histMixPM->Add(htmp);
+ }
+
+ if (!histMixPM){
+ AliError("For R-factor correction the mixed event histograms are requires. No +- histogram found");
+ return;
+ }
+ histMixPM->Sumw2();
+
+ // rebin the histograms
+ if (fRebin>1) {
+ histMixPP->Rebin(fRebin);
+ histMixMM->Rebin(fRebin);
+ histMixPM->Rebin(fRebin);
+ }
+
+ // fill out rfactor histogram
+ for(Int_t ibin=1; ibin<=histMixPM->GetXaxis()->GetNbins(); ibin++) {
+ Float_t pp = histMixPP->GetBinContent(ibin);
+ Float_t ppe = histMixPP->GetBinError(ibin);
+ Float_t mm = histMixMM->GetBinContent(ibin);
+ Float_t mme = histMixMM->GetBinError(ibin);
+ Float_t pm = histMixPM->GetBinContent(ibin);
+ Float_t pme = histMixPM->GetBinError(ibin);
+
+ Float_t background = 2*TMath::Sqrt(pp*mm);
+ // do not use it since ME could be weighted err!=sqrt(entries)
+ // Float_t ebackground = TMath::Sqrt(mm+pp);
+ Float_t ebackground = TMath::Sqrt(mm/pp*ppe*ppe + pp/mm*mme*mme);
+ if (fMethod==kLikeSignArithm){
+ //Arithmetic mean instead of geometric
+ background=(pp+mm);
+ ebackground=TMath::Sqrt(ppe*ppe+mme*mme);
+ if (TMath::Abs(ebackground)<1e-30) ebackground=1;
+ }
+
+ Float_t rcon = 1.0;
+ Float_t rerr = 0.0;
+ if(background>0.0) {
+ rcon = pm/background;
+ rerr = TMath::Sqrt((1./background)*(1./background) * pme*pme +
+ (pm/(background*background))*(pm/(background*background)) * ebackground*ebackground);
+ }
+ fHistRfactor->SetBinContent(ibin, rcon);
+ fHistRfactor->SetBinError(ibin, rerr);
+ }
+
+ fHistBackground->Multiply(fHistRfactor);
+
+ if (histMixPP) delete histMixPP;
+ if (histMixMM) delete histMixMM;
+ if (histMixPM) delete histMixPM;
+ }
+
//scale histograms to match integral between fScaleMin and fScaleMax
// or if fScaleMax < fScaleMin use fScaleMin as scale factor
- if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
+ if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
+ else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
else if (fScaleMin>0.){
fScaleFactor=fScaleMin;
fHistBackground->Scale(fScaleFactor);
}
//subract background
+ fHistSignal->Add(fHistDataPM);
fHistSignal->Add(fHistBackground,-1);
- // signal
- fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
- fHistSignal->FindBin(fIntMax), fErrors(0));
// background
fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
fHistBackground->FindBin(fIntMax),
fErrors(1));
- printf("%f %f\n",fValues(0),fValues(1));
+
+ // signal depending on peak description method
+ DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
+ //printf("%f %f\n",fValues(0),fValues(1));
// S/B and significance
- SetSignificanceAndSOB();
+ // SetSignificanceAndSOB();
fProcessed = kTRUE;
}
void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
{
//
- // event mixing
+ // event mixing of +- and -+
//
- fHistDataPM = (TH1*)(arrhist->At(0))->Clone("histPMSE"); // +- SE
- fHistDataME = (TH1*)(arrhist->At(1))->Clone("histPMME"); // +- ME
- fHistDataPM->Sumw2();
- fHistDataME->Sumw2();
- fHistDataPM->SetDirectory(0);
- fHistDataME->SetDirectory(0);
- // rebin the histograms
- if (fRebin>1) {
- fHistDataPM->Rebin(fRebin);
- fHistDataME->Rebin(fRebin);
+ if (!arrhist->At(AliDielectron::kEv1PM) || !(arrhist->At(AliDielectron::kEv1MEv2P) || arrhist->At(AliDielectron::kEv1PEv2M)) ){
+ AliError("Either OS or mixed histogram missing");
+ return;
}
- fHistSignal = new TH1D("HistSignal", "Mixed events background substracted signal",
- fHistDataPM->GetXaxis()->GetNbins(),
- fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
- fHistSignal->SetDirectory(0);
- fHistBackground = new TH1D("HistBackground", "background contribution from mixed events",
- fHistDataPM->GetXaxis()->GetNbins(),
- fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
- fHistBackground->SetDirectory(0);
+ delete fHistDataPM; fHistDataPM=0x0;
+ delete fHistDataME; fHistDataME=0x0;
+ delete fHistBackground; fHistBackground=0x0;
+
+ fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
+ fHistDataPM->Sumw2();
+ fHistDataPM->SetDirectory(0x0);
- // fill out background and subtracted histogram
- for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
- Float_t pm = fHistDataPM->GetBinContent(ibin);
- Float_t epm = fHistDataPM->GetBinError(ibin);
- Float_t background = fHistDataME->GetBinContent(ibin);
- Float_t ebackground = fHistDataME->GetBinError(ibin);
+ if (arrhist->At(AliDielectron::kEv1MEv2P)){
+ fHistDataME = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME"); // -+ ME
+ }
+
+ if (arrhist->At(AliDielectron::kEv1PEv2M)){
+ TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M));
+ if (!fHistDataME) fHistDataME = (TH1*)htmp->Clone("histMPME"); // -+ ME
+ else fHistDataME->Add(htmp);
+ }
- fHistSignal->SetBinContent(ibin, pm);
- fHistSignal->SetBinError(ibin, epm);
- fHistBackground->SetBinContent(ibin, background);
- fHistBackground->SetBinError(ibin, ebackground);
+ fHistBackground = (TH1*)fHistDataME->Clone("ME_Background");
+ fHistBackground->SetDirectory(0x0);
+ fHistBackground->Sumw2();
+
+ // rebin the histograms
+ if (fRebin>1) {
+ fHistDataPM->Rebin(fRebin);
+ fHistDataME->Rebin(fRebin);
+ fHistBackground->Rebin(fRebin);
}
- if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
- else if (fScaleMin>0.){
+ //scale histograms to match integral between fScaleMin and fScaleMax
+ // or if fScaleMax < fScaleMin use fScaleMin as scale factor
+ if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
+ else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
+ else if (fScaleMin>0.){
fScaleFactor=fScaleMin;
fHistBackground->Scale(fScaleFactor);
}
-
- //subract background
- fHistSignal->Add(fHistBackground,-1);
-
- // signal
- fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
- fHistSignal->FindBin(fIntMax), fErrors(0));
+ fHistSignal=(TH1*)fHistDataPM->Clone("Signal");
+ fHistSignal->Sumw2();
+ // printf(" err: %f %f \n",fHistSignal->GetBinError(75),TMath::Sqrt(fHistSignal->GetBinContent(75)));
+ fHistSignal->Add(fHistBackground,-1.);
+ // printf(" err: %f %f \n",fHistSignal->GetBinError(75),TMath::Sqrt(fHistSignal->GetBinContent(75)));
+// // signal
+// fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
+// fHistSignal->FindBin(fIntMax), fErrors(0));
// background
fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
- fHistBackground->FindBin(fIntMax),
- fErrors(1));
- // S/B and significance
- SetSignificanceAndSOB();
+ fHistBackground->FindBin(fIntMax),
+ fErrors(1));
+
+ // signal depending on peak description method
+ DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
fProcessed = kTRUE;
}
//
// signal subtraction
//
- fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM"); // +- SE
- if (!fHistDataPM){
- AliError("Unlike sign histogram not available. Cannot extract the signal.");
+
+ if (!arrhist->At(AliDielectron::kEv1PM) || !arrhist->At(AliDielectron::kEv1PMRot) ){
+ AliError("Either OS or rotation histogram missing");
return;
}
+
+ fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE
fHistDataPM->Sumw2();
+ fHistDataPM->SetDirectory(0x0);
- fHistBackground=(TH1*)arrhist->At(10)->Clone("histRotation");
- if (!fHistBackground){
- AliError("Histgram from rotation not available. Cannot extract the signal.");
- delete fHistDataPM;
- fHistDataPM=0x0;
- return;
- }
+ fHistBackground = (TH1*)(arrhist->At(AliDielectron::kEv1PMRot))->Clone("histRotation");
fHistBackground->Sumw2();
+ fHistBackground->SetDirectory(0x0);
// rebin the histograms
if (fRebin>1) {
//scale histograms to match integral between fScaleMin and fScaleMax
// or if fScaleMax < fScaleMin use fScaleMin as scale factor
- if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
+ if (fScaleMax>fScaleMin && fScaleMax2>fScaleMin2) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax,fScaleMin2,fScaleMax2);
+ else if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
else if (fScaleMin>0.){
fScaleFactor=fScaleMin;
fHistBackground->Scale(fScaleFactor);
fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
fHistSignal->Add(fHistBackground,-1.);
+ fHistSignal->SetDirectory(0x0);
- // signal
- fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
- fHistSignal->FindBin(fIntMax), fErrors(0));
+ // // signal
+ // fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
+ // fHistSignal->FindBin(fIntMax), fErrors(0));
// background
fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
fHistBackground->FindBin(fIntMax),
fErrors(1));
- // S/B and significance
- SetSignificanceAndSOB();
+ // signal depending on peak description method
+ DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
fProcessed = kTRUE;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff