// Author: Gero Flucke <mailto:flucke@mail.desy.de>
//____________________________________
// GFDstarHistsAnalysis
// Author: Gero Flucke
// Date: May 14th, 2002
// last update: $Date: 2005/11/16 00:39:01 $
// by: $Author: flucke $
//
#include <iostream>
#include <TROOT.h>
#include <TFile.h>
#include <TH1.h>
#include <TLegend.h>
#include <TH2.h>
#include <TH3.h>
#include <TProfile.h>
#include <TF1.h>
#include <TGraphAsymmErrors.h>
#include <TString.h>
#include <TObjString.h>
#include <TAxis.h>
#include <TObjArray.h>
#include "TDatabasePDG.h"
#include "TCanvas.h"
#include "H1PhysUtils/H1DedxNLH.h"
#include "GFAnalysis/GFDstarHistsAnalysis.h"
#include "GFAnalysis/GFAxisVariables.h"
#include "GFAnalysis/GFMultiDmFitter.h"
#include "GFAnalysis/GFDstarCount.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFHistManager.h"
#include "GFUtils/GFMath.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFDstarRunPeriods.h"
#include "/afs/desy.de/user/f/flucke/h1/root/DstarDmFitter/DstarDmFitter.h"
#include "/afs/desy.de/user/f/flucke/h1/root/DstarDmFitter/D0Fitter.h"
using namespace std;
ClassImp(GFDstarHistsAnalysis)
// Bool_t GFDstarHistsAnalysis::fgNoDelete = kFALSE;
const char* GFDstarHistsAnalysis::kDefDir = "/user/flucke/dstar/histFiles/";
//"/data/usr/flucke/dstar/data/oo/histFiles/";
const char* GFDstarHistsAnalysis::kDefName = "DstarAnalysisHists";
//__________________________________________
GFDstarHistsAnalysis::GFDstarHistsAnalysis(const char* fileIn)//,const char*fileOut)
{
this->Init(fileIn);
}
//__________________________________________
GFDstarHistsAnalysis::GFDstarHistsAnalysis(const char* prefix, const char* period)
{
TString fileInString(kDefDir);
(((fileInString += prefix) += kDefName) += period) += ".root";
this->Init(fileInString.Data());
}
//__________________________________________
void GFDstarHistsAnalysis::Init(const char* fileIn)
{
fFileIn = (TFile*)gROOT->GetListOfFiles()->FindObject(fileIn);
if (!fFileIn) {
fFileIn = TFile::Open(fileIn);
}
if(!fFileIn){
this->Error("Init", "Could not open file %s", fileIn);
} else {
cout << "input file: " << fileIn << endl;
}
this->SetHistPointerNull();
fTrigger = 83;
// fFitMode = 1011;
fFitModeD0 = 11;
fD0Fitter = NULL;
fMultiDmFitter = new GFMultiDmFitter;
fMultiDmFitter->SetFitModeBins(1011);//fFitMode);
// removal of fit option 'M' for improve as this often leads to endless loops
TString fitOption = fMultiDmFitter->GetDstarDmFitter()->GetFitOption();
fitOption.ReplaceAll("M", NULL);
fMultiDmFitter->GetDstarDmFitter()->SetFitOption(fitOption);
fHistManager = new GFHistManager;
fTmpGraph = NULL;
fVariables = new TObjArray;
fVariables->Add(new TObjString("pt"));
fVariables->Add(new TObjString("eta"));
fVariables->Add(new TObjString("phi"));
fVariables->Add(new TObjString("ptSumEt"));
fVariables->Add(new TObjString("wGammaP"));
fVariables->Add(new TObjString("y"));
fVariables->Add(new TObjString("nJet"));
// fVariables->Add(new TObjString("ptJet1"));
fVariables->Add(new TObjString("ptJet2"));
fVariables->Add(new TObjString("etaJet1"));
fVariables->Add(new TObjString("etaJet2"));
fVariables->Add(new TObjString("xGamma"));
}
//____________________________________
GFDstarHistsAnalysis::~GFDstarHistsAnalysis()
{
if(fFileIn) {
// fFileIn->Close();
// delete fFileIn; //FIXME: memleak???, but maybe not the only one opening this file...
fFileIn = NULL;
}
// if(fDstarFitter) delete fDstarFitter;
if(fD0Fitter) delete fD0Fitter;
if(fHistManager) delete fHistManager;
delete fMultiDmFitter;
if(fHistSignalBack) delete fHistSignalBack;
if(fHist2D) delete fHist2D;
if(fHist1D) delete fHist1D;
if(fHistNumDstar) delete fHistNumDstar;
if(fVariables){
fVariables->Delete();
delete fVariables;
}
}
//____________________________________
void GFDstarHistsAnalysis::SetTrigger(Int_t trigger)
{
fTrigger = trigger;
}
//____________________________________
void GFDstarHistsAnalysis::SetFitModeBins(Int_t mode)
{
fMultiDmFitter->SetFitModeBins(mode);//fFitMode);
}
//____________________________________
Int_t GFDstarHistsAnalysis::GetFitModeBins() const
{
return fMultiDmFitter->GetFitModeBins();
}
//____________________________________
void GFDstarHistsAnalysis::Draw(TH1* (GFDstarHistsAnalysis::*funcPtr) (const char*,Bool_t, const char* ), Bool_t perBinW, const char *dir)
{
const Bool_t wasBatch = this->SetBatch();
GFHistArray hists;
for(Int_t i = 0; i < fVariables->GetEntriesFast(); ++i){
hists.Add((this->*funcPtr)(fVariables->At(i)->GetName(), perBinW, dir));
}
this->SetBatch(wasBatch);
fHistManager->Clear();
fHistManager->AddHists(&hists);
fHistManager->Draw();
}
//____________________________________
void GFDstarHistsAnalysis::DrawDmAll(Bool_t perBinW)
{
this->Draw(&GFDstarHistsAnalysis::DrawDm, perBinW);
}
//____________________________________
TArrayD GFDstarHistsAnalysis::DrawDm(Bool_t both)
{
// the returned array has lenght 2:
// array[0] = numDstar
// array[1] = numDstarError
// (Both values for histogram with cuts appropriate to trigger setting!!)
// fDstarFitter->ResetFitParameters();
fHistManager->Clear();
TH1* dmHist = this->GetHistPtCut();
// fDstarFitter->FitAll(dmHist);
// TArrayD array(2);
// array[0] = fDstarFitter->GetNDstar();
// array[1] = fDstarFitter->GetNDstarErr();
TArrayD array = fMultiDmFitter->FitAll(dmHist);
fHistManager->AddHist(dmHist);
if(both){
TH1* hist2 = this->GetTrigger() == 83 ? this->GetHistPt1Cut()
: this->GetHistPt2Cut();
// fDstarFitter->FitAll(hist2);
fMultiDmFitter->FitAll(hist2);
fHistManager->AddHist(hist2);
}
fHistManager->SetHistsOption("e");
fHistManager->SetHistsXTitle("m(K#pi#pi_{s}) - m(K#pi) [GeV]");
fHistManager->Draw();
return array;
}
//____________________________________
TArrayD GFDstarHistsAnalysis::DrawD0(Bool_t both)
{
// the returned array has lenght 2:
// array[0] = numD0
// array[1] = numD0Error
// (Both values for histogram with cuts appropriate to trigger setting!!)
// fDstarFitter->ResetFitParameters();
fHistManager->Clear();
TH1* d0Hist = this->GetHistD0PtCut();
if(!fD0Fitter) fD0Fitter = new D0Fitter;
fD0Fitter->Fit(0,d0Hist, kTRUE);
TArrayD array(2);
array[0] = fD0Fitter->GetND0();
array[1] = fD0Fitter->GetND0Err();
fHistManager->AddHist(d0Hist);
if(both){
TH1* hist2 = this->GetTrigger() == 83 ? this->GetHistD0Pt1Cut()
: this->GetHistD0Pt2Cut();
fD0Fitter->Fit(0,hist2, kTRUE);
fHistManager->AddHist(hist2);
}
fHistManager->SetHistsOption("e");
fHistManager->SetHistsXTitle("m(K#pi) [GeV]");
fHistManager->Draw();
return array;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawCheckD0(Int_t mode)
{
// mode 0: L4 triggered by 'not class 15'
// 1: just omit L4 cut
// 2: omit L4 cut, but apply weights
const char* name = NULL;
switch(mode){
case 0:
name = "D0not15L4";
break;
case 1:
name = "D0noL4cut";
break;
case 2:
name = "D0L4weight";
break;
}
TH1* hD0 = this->GetHist(name, "D0Checks");
if(hD0){
if(!fD0Fitter) fD0Fitter = new D0Fitter;
fD0Fitter->Fit(0, hD0);
fHistManager->Clear();
fHistManager->AddHist(hD0);
fHistManager->SetHistsOption("e");
fHistManager->SetHistsXTitle("m(K#pi) [GeV]");
fHistManager->Draw();
}
return hD0;
}
//____________________________________
TArrayD GFDstarHistsAnalysis::DrawDmDsJet(const char* forwBack, UInt_t dmRefFlag)
{
// if(dmRefFlag != 0): First fit a reference hist with all params free, then fit with fFitMode
// >=2 [default]: D* inclusive
// 1 : D* + 1 jet inclusive
// fDstarFitter->ResetFitParameters();
fHistManager->Clear();
TArrayD result(4);
TH1* dmHist = this->GetHistDs1Jet(forwBack);
if(!dmHist) return result;
if(dmRefFlag && (dmRefFlag != 1 || (forwBack && strlen(forwBack)))){
if(!fHistManager->IsBatch()){
fMultiDmFitter->FitAll(dmHist);
fHistManager->AddHist((TH1*) dmHist->Clone());
}
TH1* histRef = this->GetHistDs1JetRef(forwBack, dmRefFlag);
fMultiDmFitter->FitAll(histRef);
fHistManager->AddHist(histRef);
result = fMultiDmFitter->FitSingle(dmHist, NULL);
} else {
result = fMultiDmFitter->FitAll(dmHist);
}
fHistManager->AddHist(dmHist, 1);
fHistManager->SetHistsOption("e");
fHistManager->SetHistsXTitle("m(K#pi#pi_{s}) - m(K#pi) [GeV]");
fHistManager->Draw();
return result;
}
//____________________________________
TArrayD GFDstarHistsAnalysis::DrawDmDsDiJet(UInt_t dmRefFlag)
{
// First fit a reference hist with all params free, then fit with fFitMode
// if dmRefFla == 2 take 'D* inclusive' as reference
// fDstarFitter->ResetFitParameters();
fHistManager->Clear();
TArrayD result(4);
TH1 *dmHist = this->GetHist("dstarDiJet", "dsDiJets");
if (!dmHist) return result;
TH1* dmHistRef = NULL;
if (dmRefFlag == 2) dmHistRef = this->GetHistPtCut();
if (dmHistRef) {
fMultiDmFitter->FitAll(dmHistRef);
fHistManager->AddHist(dmHistRef);
result = fMultiDmFitter->FitSingle(dmHist, NULL);
} else {
result = fMultiDmFitter->FitAll(dmHist);
}
fHistManager->AddHist(dmHist, 1);
fHistManager->SetHistsOption("e");
fHistManager->SetHistsXTitle("m(K#pi#pi_{s}) - m(K#pi) [GeV]");
fHistManager->Draw();
return result;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDm(const char* var, Bool_t perBinW, const char *dir)
{
// if(perBinW) result is divided through bin width)
fHistManager->Clear();
fHist2D = this->GetHist2D(var, dir);
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable(GFAxisVariables::GetAxisLabel(var), NULL, NULL, perBinW);
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmPt()//Int_t fitMode)
{
fHistManager->Clear();
fHistManager->SetNumHistsY(3);
fHistManager->SetNumHistsX(3);
fHist2D = this->GetHist2D("pt");
fHist1D = this->GetHistPtCut();
TH1* h = this->DrawDmVariable("p_{t} [GeV]");
fHistManager->SetLogY(3);
return h;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmPhi()
{
fHistManager->Clear();
fHistManager->SetNumHistsY(5);
fHist2D = this->GetHist2D("phi");
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("#phi");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmEta()
{
fHistManager->Clear();
fHistManager->SetNumHistsY(3);
fHist2D = this->GetHist2D("eta");
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("#eta");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmWgammaP()
{
fHistManager->Clear();
// fHistManager->SetNumHistsY(3);
fHist2D = this->GetHist2D("wGammaP");
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("W_{#gammap} [GeV]");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDm(const char* var1, const char* var2,
Int_t firstBin, Int_t lastBin, Bool_t perBinW)
{
// if(perBinW) result is divided through bin width)
fHistManager->Clear();
Bool_t flip = kFALSE;
TH3 *h3D = static_cast<TH3*>(this->GetHist2Diff(flip, var1, var2, "Hist2D", NULL));
if(!h3D) return NULL;
TAxis *axis = (flip ? h3D->GetYaxis() : h3D->GetZaxis());
if(!GFHistManip::BinsAreSubRange(axis, firstBin, lastBin)){
// could be extended to tolerate under/overflow?
this->Error("DrawDm", "problem with bin range %d to %d", firstBin, lastBin);
delete h3D;
return NULL;
}
axis->SetRange(firstBin, lastBin);
TString opt(flip?"zx":"yx");
if(h3D->GetSumw2N()) opt += 'e';// weighted 3D hist needs to propagate weights
fHist2D = static_cast<TH2*>(h3D->Project3D(opt));
this->AddOutOfRange(h3D, firstBin, lastBin, flip, fHist2D);
fHist2D->SetTitle(Form("%.2f #leq %s < %.2f", axis->GetBinLowEdge(firstBin),
GFAxisVariables::GetAxisLabel(var2), axis->GetBinUpEdge(lastBin)));
fHist1D = this->GetHistPtCut();
delete h3D;
TH2 *h2 = fHist2D;
TH1 *h = this->DrawDmVariable(GFAxisVariables::GetAxisLabel(var1),NULL,NULL, perBinW);
h->SetTitle(Form("%s, %s", h->GetTitle(), h2->GetTitle()));
cout << "| " << h2->GetTitle()
<< "n------------------------------------------------------------------------------"
<< endl;
return h;
}
//____________________________________
void GFDstarHistsAnalysis::AddOutOfRange(TH3 *h3D, Int_t firstBin, Int_t lastBin,
Bool_t flip, TH2* h2D) const
{
// add everything which is not in firstBin->lastBin range into
// y-underflow of h2D
if (!h3D || !h2D || firstBin > lastBin) {
this->Error("AddOutOfRange", "problem with input: %p, %p, %d %d",
h3D, h2D, firstBin, lastBin);
return;
}
TAxis *axis = (flip ? h3D->GetYaxis() : h3D->GetZaxis());
const Bool_t err = h3D->GetSumw2N();
const TString opt(Form("x%s", err?"e":""));//weighted 3D hist: propagate weights
TH1 *h1D1 = NULL;
TH1 *h1D2 = NULL;
// if (0 < firstBin) { // set range with 0 does not work...
// axis->SetRange(0, firstBin-1);
if (1 < firstBin) {
axis->SetRange(1, firstBin-1);
h1D1 = h3D->Project3D(opt);
}
// const Int_t overFlowBin = axis->GetNbins() + 1; // setRange(overflow) ...
// if (lastBin < overFlowBin) { /// .. does not work
// axis->SetRange(lastBin + 1, overFlowBin);
const Int_t lastAxBin = axis->GetNbins();
if (lastAxBin > lastBin) {
axis->SetRange(lastBin + 1, lastAxBin);
h1D2 = h3D->Project3D(opt+'2');
}
const Int_t nBinsX = h3D->GetXaxis()->GetNbins();
for (Int_t iBin = 0; iBin <= nBinsX + 1; ++iBin) {
Double_t cont = 0.;
Double_t err2 = 0.;
if (h1D1) {
cont += h1D1->GetBinContent(iBin);
if (err) err2 += h1D1->GetBinError(iBin)*h1D1->GetBinError(iBin);
}
if (h1D2) {
cont += h1D2->GetBinContent(iBin);
if (err) err2 += h1D2->GetBinError(iBin)*h1D2->GetBinError(iBin);
}
// HACK: since over-/underflow does not work for SetRange, add them by hand:
TAxis *othAxis = (flip ? h3D->GetZaxis() : h3D->GetYaxis());
for (Int_t i = 0; i <= othAxis->GetNbins() + 1; ++i) {
cont += (flip ? h3D->GetBinContent(iBin, 0, i) : h3D->GetBinContent(iBin, i, 0));
cont += (flip ? h3D->GetBinContent(iBin, lastAxBin+1, i) :
h3D->GetBinContent(iBin, i, lastAxBin+1));
if (err) {
Double_t tmpE = (flip ? h3D->GetBinError(iBin, 0, i) : h3D->GetBinError(iBin, i, 0));
err2 += tmpE*tmpE;
tmpE = (flip ? h3D->GetBinError(iBin, lastAxBin+1, i) :
h3D->GetBinError(iBin, i, lastAxBin+1));
err2 += tmpE*tmpE;
}
}
// end HACK
cont += h2D->GetBinContent(iBin, 0);
h2D->SetBinContent(iBin, 0, cont);
if (err) {
err2 += h2D->GetBinError(iBin, 0)*h2D->GetBinError(iBin, 0);
h2D->SetBinError(iBin, 0, TMath::Sqrt(err2));
}
}
delete h1D1;
delete h1D2;
}
//____________________________________
void GFDstarHistsAnalysis::AddOutOfRange(const TH2 *h2D, Int_t firstBin, Int_t lastBin,
Bool_t flip, TH1 *h1D) const
{
// add everything which is not in firstBin->lastBin range into
// underflow of h1D
if (!h2D || !h1D || firstBin > lastBin) {
this->Error("AddOutOfRange", "problem with input: %p, %p, %d %d",
h2D, h1D, firstBin, lastBin);
return;
}
const TAxis *axis = (flip ? h2D->GetXaxis() : h2D->GetYaxis());
const Bool_t err = h2D->GetSumw2N();
Double_t all = 0.;
Double_t err2 = 0.;
const Int_t overFlowBin = axis->GetNbins() + 1;
for (Int_t iBin = 0; iBin <= overFlowBin; ++iBin) {
if (iBin >= firstBin && iBin <= lastBin) continue;
const TAxis *othAxis = (flip ? h2D->GetYaxis() : h2D->GetXaxis());
for (Int_t iOthBin = 0; iOthBin <= othAxis->GetNbins() + 1; ++iOthBin) {
all += (flip ? h2D->GetBinContent(iBin, iOthBin) : h2D->GetBinContent(iOthBin, iBin));
if (err) {
const Double_t e = (flip ? h2D->GetBinError(iBin, iOthBin) :
h2D->GetBinError(iOthBin, iBin));
err2 += e*e;
}
}
}
h1D->SetBinContent(0, all + h1D->GetBinContent(0) + h1D->GetBinContent(overFlowBin));
if (err) {
err2 += h1D->GetBinError(0)*h1D->GetBinError(0);
err2 += h1D->GetBinError(overFlowBin)*h1D->GetBinError(overFlowBin);
h1D->SetBinError(0, TMath::Sqrt(err2));
}
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmDs1Jet(const char* varJetDs, const char* forwBack,
Bool_t perBinW, UInt_t dmRefFlag)
// Bool_t dmRefIsAllDsJet)
{
// varJetDs: Pt/Eta + Jet/Ds, Deta, Dphi, forwBack:'Forw', 'Back' or empty ''
// // dmRefIsAllDsJet: if true,
// dmRefFlag gives reference hist for determining fixed parameters of fit in bins
// >=2 [default]: D* inclusive
// 1 : D* + 1 jet inclusive
// 0 : determined by 'forwBack'
fHistManager->Clear();
fHist2D = (TH2*) this->GetHist(Form("dstar1Jet%s%s", varJetDs, forwBack) , "dstarJets");
fHist1D = this->GetHistDs1JetRef(forwBack, dmRefFlag);
if(!fHist2D || !fHist1D){
this->Error("DrawDmDs1Jet", "missed a histogram");
return NULL;
}
TH1* h = this->DrawDmVariable(GFAxisVariables::GetAxisLabel(varJetDs), NULL, NULL, perBinW);
h->SetName(Form("%s%s", forwBack, h->GetName()));
return h;
}
//____________________________________
void GFDstarHistsAnalysis::DrawDmDs1JetAll(const char* forwBack)
{
fHistManager->Clear();
GFHistManager* realMan = fHistManager;
GFHistManager tmp;
tmp.SetBatch(kTRUE);
fHistManager = &tmp;
TObjArray var;
var.Add(new TObjString("EtaDs"));
var.Add(new TObjString("PtDs"));
var.Add(new TObjString("EtaJet"));
var.Add(new TObjString("PtJet"));
var.Add(new TObjString("Dphi"));
var.Add(new TObjString("Deta"));
var.Add(new TObjString("Dpt"));
for(Int_t i = 0; i < var.GetEntriesFast(); ++i){
TH1* h = this->DrawDmDs1Jet(var[i]->GetName(), forwBack);
realMan->AddHist(h, 0);
realMan->AddLayer(fHistManager, 0);
fHistManager->Clear();
}
var.Delete();
fHistManager = realMan;
fHistManager->Draw();
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmDsDiJet(const char* varDiJet,
Bool_t perBinW, UInt_t dmRefFlag)
{
// varJetDs: Pt/Eta + DsJet/OthJet/Ds, Deta, Dphi, PtDijet, CosTh, xGam,
// DeltaR, DeltaRNoDjetDs
// dmRefFlag gives reference hist for determining fixed parameters of fit in bins
// >=2 [default]: D* inclusive
// 1 : D* dijet inclusive
// 0 : just the projection of the 2D-hist
fHistManager->Clear();
fHist2D = (TH2*) this->GetHist(Form("dstarDiJet%s", varDiJet) ,"dsDiJets");
switch(dmRefFlag){
case 0:
fHist1D = NULL;
break;
case 1:
fHist1D = this->GetHist("dstarDiJet" ,"dsDiJets");
break;
default:
this->Warning("DrawDmDsDiJet", "Don't know dmRefFlag %d, use 2",dmRefFlag);
case 2:
fHist1D = this->GetHistPtCut();
}
if(!fHist2D || (!fHist1D && dmRefFlag != 0)){
this->Error("DrawDmDsDiJet", "missed a histogram");
return NULL;
}
return this->DrawDmVariable(GFAxisVariables::GetAxisLabel(varDiJet), NULL, NULL, perBinW);
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmBgFinder()
{
fHistManager->Clear();
fHist2D = static_cast<TH2*>(this->GetHist("bgBits", "bg"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("BG-finder");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmTriggers(Bool_t allHists)
{
// selections:
fHistManager->Clear();
fHistManager->SetBatch();
// fHistManager->SetNumHistsX(4);
// fHistManager->SetNumHistsY(6);
TH1* hSave = this->GetHistPtCut();
GFHistManager* origMan = fHistManager;
GFHistManager tmp; tmp.SetBatch();
fHistManager = &tmp;
/*
fHist1D = hSave; // ???
fHist2D = this->GetHist2DTriggers("L4v");
TH1* hL4v = this->DrawDmVariable("L4vST");
if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
fHist1D = hSave; // ???
fHist2D = this->GetHist2DTriggers("L1ac");
TH1* hL1ac = this->DrawDmVariable("L1acST");
if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
fHist1D = hSave; // ??? SetHistPointerNull
fHist2D = this->GetHist2DTriggers("L1L4");
TH1* hL1L4 = this->DrawDmVariable("L1acL4vST");
if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
const Int_t oldLast = origMan->GetNumLayers();
origMan->AddHist(hL1L4, oldLast, "L1 ac && L4 v");
origMan->AddHistSame(hL4v, oldLast, 0, "L4 verified");
origMan->AddHistSame(hL1ac, oldLast, 0, "L1 actual");
origMan->AddHistsOption("HIST", oldLast);
origMan->AddLegend(oldLast, 0, Form("full D* selection: S%d", fTrigger)); // false?
origMan->SetHistsXTitle("ST", oldLast);
*/
if(fTrigger != -1){
TString notTrig("Ntrig");
fHist1D = hSave; // ???
fHist2D = this->GetHist2DTriggers("L4v"+notTrig);
TH1* hL4vNtrig = this->DrawDmVariable("L4vNtrST");
// if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
fHist2D = this->GetHist2DTriggers("L1ac"+notTrig);
fHist1D = hSave; // ???
TH1* hL1acNtrig = this->DrawDmVariable("L1acNtrST");
// if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
fHist2D = this->GetHist2DTriggers("L1L4"+notTrig);
fHist1D = hSave; // ??? SetHistPointerNull
TH1* hL1L4Ntrig = this->DrawDmVariable("L1L4NtrST");
if(allHists) origMan->AddLayers(&tmp);
tmp.Clear();
const Int_t oldLast = origMan->GetNumLayers();
origMan->AddHist(hL1L4Ntrig, oldLast, "L1 ac && L4 v");
origMan->AddHistSame(hL4vNtrig, oldLast, 0, "L4 verified");
origMan->AddHistSame(hL1acNtrig, oldLast, 0, "L1 actual");
origMan->AddHistsOption("HIST", oldLast);
origMan->AddLegend(oldLast, 0, Form("visible range: S%d", fTrigger)); // false?
origMan->SetHistsXTitle("ST", oldLast);
}
fHistManager = origMan;
fHistManager->SetBatch(kFALSE);
fHistManager->Draw();
// return hL1L4;
return NULL;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmDoubleTrigPt()
{
fHistManager->Clear();
fHist2D = this->GetHist2DDoubleTrigPt();
fHist1D = fHist2D->ProjectionX();//this->GetHistPtCut();
fHist1D->Draw();
return NULL;
fHistManager->AddHist(fHist1D,3);
return this->DrawDmVariable("subtrigger");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmEventClasses(Int_t mode)
{
// mode = 0: all events/classes (default)
// = 1: weight has to be == 1
// else: weight has to be > 1
fHistManager->Clear();
fHistManager->SetNumHistsX(4);
fHistManager->SetNumHistsY(4);
TH3* hist = static_cast<TH3*>(this->GetHist("eventClass"));
switch(mode){
case 0:
break;
case 1:
hist->GetZaxis()->SetRange(2,2);
break;
default:
TAxis* z = hist->GetZaxis();
z->SetRange(3,z->GetLast()); //?
}
fHist2D = static_cast<TH2*>(hist->Project3D("yx"));
delete hist; hist = NULL;// get rid of axis range for next ?
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("event classes");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmEventClassesNo15()
{
fHistManager->Clear();
fHistManager->SetNumHistsX(4);
fHistManager->SetNumHistsY(4);
fHist2D = static_cast<TH2*>(this->GetHist("eventClassNo15"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("event classes (no 15)");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmEventClassesNoL4Dstar()
{
fHistManager->Clear();
fHistManager->SetNumHistsX(4);
fHistManager->SetNumHistsY(4);
fHist2D = static_cast<TH2*>(this->GetHist("eventClassNotL4"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("event classes (no L4 D*)");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmWeight()
{
fHistManager->Clear();
fHistManager->SetNumHistsX(5);
fHistManager->SetNumHistsY(5);
// TH3* hist = static_cast<TH3*>(this->GetHist("eventClass")); // double
// fHist2D = static_cast<TH2*>(hist->Project3D("zx")); // counting!
fHist2D = static_cast<TH2*>(this->GetHist("eventWeight"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("weights");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmL4OpenCharm()
{
fHistManager->Clear();
fHistManager->SetNumHistsX(4);
fHistManager->SetNumHistsY(4);
fHist2D = static_cast<TH2*>(this->GetHist("openCharmL4tgP"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("L4 open charm finder");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmPtDsOverSumEt(const char* histName)
{
// PtDstarOverSumETDm, sumETDm
fHistManager->Clear();
fHistManager->SetNumHistsY(3);
fHist2D = //this->GetHist2DPtDsOverSumEt();
static_cast<TH2*>(this->GetHist(histName));
//static_cast<TH2*>(this->GetHist("PtDstarOverSumETDm"));
//static_cast<TH2*>(this->GetHist("sumETDm"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable(histName);//"p_{t}(D*) / #sum(E_{t}^{#theta > 10})");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmTest()
{
fHistManager->Clear();
fHistManager->SetNumHistsY(3);
fHistManager->SetNumHistsX(4);
TFile* f = TFile::Open("phtaghists.root");
fHist2D = (TH2*) f->Get("photdm");
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable("E_{#gamma tag}");
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDm2D(const char* variablePair)
{
// yeta, wGammaPeta
TH3* h3 = this->GetHist3D(variablePair);
// Int_t lastBin = h3->GetZaxis()->FindBin(kStartDmBackgr);
// Int_t firstBin = h3->GetZaxis()->FindBin(kStartDmSignal);
// cout << firstBin << " " << lastBin << endl;
// cout << kStartDmBackgr << " " << kStartDmSignal << endl;
// h3->GetXaxis()->SetRange(8, 14);//firstBin, lastBin-1);
TH2* h2 = static_cast<TH2*>(h3->Project3D("yz")); // its a TH2*!
h2->SetOption("BOX");
GFHistArray y;
for(Int_t i = 1; i <= h2->GetNbinsY(); ++i){
TString name(variablePair);
TH1* hist1D = h2->ProjectionX((name+="y")+=i, i, i);
y.Add(hist1D);
}
GFHistArray x;
for(Int_t i = 1; i <= h2->GetNbinsX(); ++i){
TString name(variablePair);
TH1* hist1D = h2->ProjectionY((name+="x")+=i, i, i);
x.Add(hist1D);
}
fHistManager->Clear();
fHistManager->AddHist(h2);
fHistManager->AddHists(&y,1);
fHistManager->AddHists(&x,2);
// fHistManager->AddHist(h3);
fHistManager->Draw();
return h2;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawD0(const char* variable)
{
fHistManager->Clear();
fHist2D = this->GetHist2DD0(variable);
fHist1D = this->GetHistD0PtCut();
return this->DrawD0Variable(GFAxisVariables::GetAxisLabel(variable));
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawD0NoL4Bias(const char* variable)
{
if(fTrigger != -1) {
this->Error("DrawD0Check", "need trigger -1. not %d", fTrigger);
return NULL;
}
const TString hName(Form("%sHist2DnoL4biasD0", variable));
fHist2D = static_cast<TH2*>(this->GetHist(hName.Data(), "D0Checks"));
if(!fHist2D) return NULL;
fHist1D = fHist2D->ProjectionX("_px", -1, -1, "e");
fHistManager->Clear();
return this->DrawD0Variable(GFAxisVariables::GetAxisLabel(variable));
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawD0Variable(const char* varName)
{
TString name(varName);
fHistNumDstar = this->CreateHistYBins(fHist2D);
// fHistSignalBack = this->CreateHistYBins(fHist2D);
fHistNumDstar->SetXTitle(varName);
fHistNumDstar->SetTitle("N(D^{0}) / #Delta"+name);
this->MakeStringName(name+=fTrigger); // changes 'name'!
fHistNumDstar->SetName("numD0"+name);
// fHistSignalBack->SetXTitle(varName);
// fHistSignalBack->SetTitle("S/B");
// fHistSignalBack->SetName("signalBack"+name);
if(!fD0Fitter) fD0Fitter = new D0Fitter;
fD0Fitter->Fit(0, fHist1D, kTRUE);//true: determine new start values
GFHistArray* hists1D = this->CreateHistArray1D(fHist2D, varName);
this->FitD0Hists(hists1D);
GFHistManip::CorrectForBinWidth(fHistNumDstar);
fHist2D->SetYTitle(varName);
fHistManager->AddHists(hists1D,0);
fHistManager->AddHist(fHist1D,1);
fHistManager->SetHistsOption("e");
fHistManager->AddHist(fHist2D,2);
fHistManager->SetHistsOption("BOX",2);
fHistManager->SetHistsXTitle("m(K#pi) [GeV]");
// fHistManager->AddHist(fHistSignalBack,2);
fHistManager->AddHist(fHistNumDstar,3);
fHistManager->Draw();
// fHistManager->WriteHistos(fFileOut);
delete hists1D;
TH1* result = fHistNumDstar;
this->SetHistPointerNull();
TString title("N(D0) in bins of ");
this->PrintTable(result, title += varName);
return result;
}
//________________________________________________
TH1* GFDstarHistsAnalysis::DrawTrack(const char* varTrackName)
{
//var+track: PtK (Pt, Theta, Phi, Length, Start, Dca, Nhit (Dz0?), + K, Pi Pis)
fHist2D = this->GetHist2DTrack(varTrackName);
fHist1D = this->GetHistPtCut();
return this->DrawTrackReally(GFAxisVariables::GetAxisLabel(varTrackName));
}
//____________________________________
void GFDstarHistsAnalysis::DrawDedx(Bool_t lhCuts)
{
// lhCuts = true (default): draw seperate colours for entries that passed or failed
// the LH-cut on dE/dx
fHistManager->Clear();
TH1* histK = this->GetHist("trackHistDedxK", "tracks");
TH1* histPi = this->GetHist("trackHistDedxPi", "tracks");
TH1* histPis = this->GetHist("trackHistDedxPis", "tracks");
fHistManager->AddHist(histK);
fHistManager->AddHist(histPi);
fHistManager->AddHist(histPis);
TH1* histAll = this->GetHist("trackHistDedxAll", "tracks");
fHistManager->AddHist(histAll, 2);
TH1* histPisSR = this->GetHist("trackHistDedxPis2", "tracks");
fHistManager->AddHist(histPisSR, 1);
fHistManager->SetNumHistsXY(1,3);
fHistManager->SetHistsXTitle("p [GeV]");
fHistManager->SetHistsYTitle("dE/dx [mip]");
histK->GetXaxis()->SetMoreLogLabels();
histK->GetXaxis()->SetNoExponent();
histPi->GetXaxis()->SetMoreLogLabels();
histPi->GetXaxis()->SetNoExponent();
histPis->GetXaxis()->SetMoreLogLabels();
histPis->GetXaxis()->SetNoExponent();
histAll->GetXaxis()->SetMoreLogLabels();
histAll->GetXaxis()->SetNoExponent();
histPisSR->GetXaxis()->SetMoreLogLabels();
histPisSR->GetXaxis()->SetNoExponent();
if(lhCuts){
histK->SetMarkerColor(kWhite);
histPi->SetMarkerColor(kWhite);
histPis->SetMarkerColor(kWhite);
histK = this->GetHist("trackHistDedxLHKK", "tracks");
histPi = this->GetHist("trackHistDedxLHPiPi", "tracks");
histPis = this->GetHist("trackHistDedxLHPiPis", "tracks");
fHistManager->AddHistSame(histK, 0, 0);
fHistManager->AddHistSame(histPi, 0, 1);
fHistManager->AddHistSame(histPis, 0, 2);
TH1* histK2 = this->GetHist("trackHistDedxNoLHKK", "tracks");
TH1* histPi2 = this->GetHist("trackHistDedxNoLHPiPi", "tracks");
TH1* histPis2 = this->GetHist("trackHistDedxNoLHPiPis", "tracks");
histK2->SetMarkerColor(kRed);
histPi2->SetMarkerColor(kRed);
histPis2->SetMarkerColor(kRed);
fHistManager->AddHistSame(histK2, 0, 0);
fHistManager->AddHistSame(histPi2, 0, 1);
fHistManager->AddHistSame(histPis2, 0, 2);
}
TLegend* legends[5];
for(Int_t i = 0; i < 3; ++i){
legends[i] = fHistManager->AddLegend(0, i, NULL, kFALSE);
}
legends[3] = fHistManager->AddLegend(2, 0, NULL, kFALSE);
legends[4] = fHistManager->AddLegend(1, 0, NULL, kFALSE);
fHistManager->Draw();
TCanvas* can = fHistManager->GetCanvas(0);
if(!can) return;
TF1* dedxK = H1DedxNLH::CreateDedxParam(H1DedxNLH::kKaon, histPis->GetXaxis()->GetXmin(),
histK->GetXaxis()->GetXmax());
// dedxK->SetLineColor(kRed);
TF1* dedxPi = H1DedxNLH::CreateDedxParam(H1DedxNLH::kPion, histPis->GetXaxis()->GetXmin(),
histK->GetXaxis()->GetXmax());
// dedxPi->SetLineColor(kBlue);
TF1* dedxProton = H1DedxNLH::CreateDedxParam(H1DedxNLH::kProton, histPis->GetXaxis()
->GetXmin(), histK->GetXaxis()->GetXmax());
// dedxProton->SetLineColor(kGreen);
TList funcs; funcs.Add(dedxPi); funcs.Add(dedxK); funcs.Add(dedxProton);
TF1* dedxPisCut = H1DedxNLH::CreateDedxParam(H1DedxNLH::kPion, histPis->GetXaxis()->GetXmin(),
histK->GetXaxis()->GetXmax(), 1.);
dedxPisCut->SetLineWidth(1);
dedxPisCut->SetName("dedxUpperBoundPis");
funcs.Add(dedxPisCut);
TIter next(&funcs); Int_t style = 1;
while(TF1* f = static_cast<TF1*>(next())){
f->SetFillStyle(0); // why needed for 2nd object...?
f->SetLineStyle(style++);
f->SetLineColor(kBlue);
}
funcs.Remove(dedxPisCut);
for(Int_t i = 1; i < 6 ; ++i){
if(i == 5) fHistManager->GetCanvas(2)->cd(1);
else if(i == 4) fHistManager->GetCanvas(1)->cd(1);
else can->cd(i);
gPad->SetLogx();
funcs.ForEach(TF1,DrawCopy)("SAME");
// dedxK->DrawCopy("SAME");
// dedxPi->DrawCopy("SAME");
// dedxProton->DrawCopy("SAME");
legends[i-1]->AddEntry(dedxPi, "pion", "l");
legends[i-1]->AddEntry(dedxK, "kaon", "l");
legends[i-1]->AddEntry(dedxProton, "proton", "l");
}
fHistManager->GetCanvas(1)->cd(1);
dedxPisCut->DrawCopy("SAME");
}
//____________________________________
void GFDstarHistsAnalysis::DrawDedxDiff()
{
fHistManager->Clear();
GFHistArray all;
GFHistArray not1; // not OK
all.Add(this->GetHist("trackHistDedxOKmExpK","tracks"));
not1.Add(this->GetHist("trackHistDedxNotOKmExpK", "tracks"));
all.Add(this->GetHist("trackHistDedxOKmExpPi","tracks"));
not1.Add(this->GetHist("trackHistDedxNotOKmExpPi", "tracks"));
all.Add(this->GetHist("trackHistDedxOKmExpPis","tracks"));
not1.Add(this->GetHist("trackHistDedxNotOKmExpPis", "tracks"));
for(Int_t i = 0; i < all.GetEntriesFast(); ++i){
all[i]->Add(not1[i]);
TH1* hNot1D = static_cast<TH2*>(not1[i])->ProjectionY();
TH1* hAll1D = static_cast<TH2*>(all[i])->ProjectionY();
fHistManager->AddHist(hAll1D, 0, "all tracks", "l");
fHistManager->AddHistSame(hNot1D, 0, i, "rejected tracks", "lf");
hNot1D->SetFillStyle(3004);
hNot1D->SetFillColor(kRed);
Int_t step = 1, count = 1;
for(Int_t j = 1; j <= not1[i]->GetNbinsX(); ++count, j += step){
TH1* hN =static_cast<TH2*>(not1[i])->ProjectionY(Form("histNot%d%d",i,count), j,j+step-1);
TH1* hA =static_cast<TH2*>(all[i])->ProjectionY(Form("histAll%d%d",i,count), j, j+step-1);
fHistManager->AddHist(hA, count, "all tracks", "l");
fHistManager->AddHistSame(hN, count, i, "rejected tracks", "lf");
hA->SetTitle(Form("(dE/dx - param) for %s: %.2f <= p < %.2f",
(i == 0 ? "K" : (i == 1 ? "#pi" : "#pi_{s}")),
not1[i]->GetBinLowEdge(j), not1[i]->GetBinLowEdge(j+step)));
hN->SetFillStyle(3004);
hN->SetFillColor(kRed);
++step;
}
}
fHistManager->SetNumHistsXY(2,2);
fHistManager->Draw();
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDedxLH(const char* track, Int_t pBinL, Int_t pBinU)
{
// // K, Pi, Pis
TH3* h = static_cast<TH3*>(this->GetHist(Form("trackHistLH%svsP", track), "tracks"));
if(!h) return NULL;
if(pBinL > pBinU || pBinL < 1 || pBinU > h->GetNbinsY()) {
this->Error("DrawDedxLH", "bins must be in range 1...%d and pBinL <= pBinU",h->GetNbinsY());
return NULL;
}
h->GetYaxis()->SetRange(pBinL, pBinU);
fHist2D = static_cast<TH2*>(h->Project3D((h->GetSumw2N() ? "zxe" : "zx"))); // "e" for errors
TString txt = fHist2D->GetName();
GFHistManip::MakeUniqueName(txt);
fHist2D->SetName(txt);
txt += Form("%s, %.1f #leq p < %.1f", h->GetTitle(), h->GetYaxis()->GetBinLowEdge(pBinL),
h->GetYaxis()->GetBinLowEdge(pBinU+1));
fHist2D->SetTitle(txt);
fHist1D = this->GetHistPtCut();
TH1* result = this->DrawTrackReally(GFAxisVariables::GetAxisLabel(Form("LH%s", track)));
txt = Form("%s, %.1f #leq p(%s) < %.1f", result->GetTitle(),
h->GetYaxis()->GetBinLowEdge(pBinL), track, h->GetYaxis()->GetBinLowEdge(pBinU+1));
result->SetTitle(txt);
return result;
}
//____________________________________
void GFDstarHistsAnalysis::DrawYresolution(const char* var, Bool_t mcYtag)
{
//
TString varUp(var);
TString start(mcYtag ? "yjbResYtag" : "yjbRes");
fHist2D = this->GetHist2D(start + varUp);
if(!fHist2D) {
this->Error("DrawYresolution", "Couldn't get 2D hist!");
return;
}
TH1* histInVarBinsMean = this->CreateHistYBins(fHist2D, ", mean");
TH1* histInVarBinsSigma = this->CreateHistYBins(fHist2D, ", #sigma");
GFHistArray* hists1D = this->CreateHistArray1D(fHist2D,
GFAxisVariables::GetAxisLabel(varUp));
this->FitGaus(hists1D, histInVarBinsMean, histInVarBinsSigma);
TProfile* prof = fHist2D->ProfileY(TString(fHist2D->GetName()) += "Prof");
// prof->SetXTitle(var);
// prof->SetYTitle("(y_{jb} - y_{gen/tag})/y_{gen/tag}");
prof->SetYTitle("mean #pm RMS");
fHist1D = this->GetHist(start + "Hist");
TH2* h = this->GetHist2D(start);
h->SetXTitle("y_{jb}");
h->SetYTitle("y_{gen/tag}");
fHistManager->Clear();
fHistManager->AddHist(fHist2D);
fHistManager->SetHistsYTitle(GFAxisVariables::GetAxisLabel(varUp));
fHistManager->AddHist(fHist1D,3);
fHistManager->AddHists(hists1D, 1);
fHistManager->SetHistsXTitle("(y_{jb} - y_{gen/tag})/y_{gen/tag}");
fHistManager->AddHist(prof, 2);
fHistManager->AddHist(histInVarBinsMean, 2);
fHistManager->AddHist(histInVarBinsSigma, 2);
fHistManager->SetHistsXTitle(var, 2);
fHistManager->AddHist(h, 0);
fHistManager->AddHistsOption("BOX", 0);
fHistManager->SetNumHistsX(5);
fHistManager->SetNumHistsY(5);
fHistManager->Draw();
fHist1D = fHist2D = NULL;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawJetComposition(const char* component)
{
//HadCl,EmCl,Track + 1/2/3/.../Pt1/Pt2/...
fHistManager->Clear();
TString name(Form("jetCompEt%s", component));
fHist2D = static_cast<TH2*>(this->GetHist(name.Data(), "jetComp"));
if(!fHist2D) return NULL;
fHist1D = this->GetHistPtCut();
// new:
TH1* hComp = fMultiDmFitter->Fit(fHist2D, fHist1D);
this->MakeStringName(name+=fTrigger); // changes 'name'!
hComp->SetName(name);
// GFHistManip::Normalise(hComp);
TH1* hNorm = GFHistManip::CreateProjection(fHist2D, "hNorm", -1, -1);
TArrayD nDstarErr = fMultiDmFitter->FitSingle(hNorm, NULL);
if(nDstarErr[0]) GFHistManip::Scale(hComp, 1./nDstarErr[0]);
else this->Error("DrawJetComposition", "no D*!");
delete hNorm; hNorm = NULL;
TString title(hComp->GetTitle());
const Ssiz_t startBin = title.Index(",");
if(startBin != kNPOS) title.Remove(0, startBin+2);
else title = "all jets";
if(strncmp(component, "Track", 5) == 0){
// title+=";p_{t}^{track}/p_{t}^{jet}";
title+=";p_{t}(track)/p_{t}";
} else if(strncmp(component, "EmCl", 4) == 0){
// title+=";p_{t}^{em}/p_{t}^{jet}";
title+=";p_{t}(em)/p_{t}";
} else if(strncmp(component, "HadCl", 5) == 0){
// title+=";p_{t}^{had}/p_{t}^{jet}";
title+=";p_{t}(had)/p_{t}";
}
hComp->SetTitle(title+=";n_{jet}/N");
fHistManager->AddHists(fMultiDmFitter->GetHists1D(),0);
fHistManager->AddHist(fHist1D,1);
fHistManager->SetHistsOption("e");
fHistManager->AddHist(fHist2D,2);
fHistManager->SetHistsOption("BOX",2);
fHistManager->AddHist(fMultiDmFitter->GetHistSignalBack(), 2);
fHistManager->AddHist(hComp,3);
fHistManager->Draw();
this->SetHistPointerNull();
// old:
// TH1* hComp = this->DrawDmVariable(name.Data(), NULL, NULL, kFALSE);// not per bin width
// hComp->SetTitle(fHist2D->GetTitle());
// GFHistManip::Normalise(hComp);
// both:
return hComp;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawJetCompSRMean(const char* component,
const char* pteta)
{
//HadCl,EmCl,Track; Eta/Pt
GFHistArray hists;
TString xyLabel;
TH2* hBins = NULL;
TString name;
if(strcmp(pteta, "Eta") == 0){
xyLabel = ";#eta(jet)";
hBins = static_cast<TH2*>(this->GetHist("dstar1JetEtaJet", "dstarJets"));
name = "jetCompEt%sSR1D%d";
} else if(strcmp(pteta, "Pt") == 0){
xyLabel = ";p_{t}(jet)";
name = "jetCompEt%sSR1DPt%d";
hBins = static_cast<TH2*>(this->GetHist("dstar1JetPtJet", "dstarJets"));
}
if(!hBins) return NULL;
for(Int_t i = 1; i <= hBins->GetNbinsY(); ++i){
hists.Add(this->GetHist(Form(name.Data(), component, i), "jetComp"));
}
if(hists.GetEntriesFast() != hBins->GetNbinsY()) return NULL;
name = Form("jetCompEt%sSR1D", component);
TH1* hOverall = this->GetHist(name.Data(), "jetComp");
if(!hOverall) {
hists.Delete();
return NULL;
}
this->MakeStringName(name+=fTrigger); // changes 'name'!
hOverall->SetName(name);
GFHistManip::Normalise(hOverall);
TH1* hResult = GFHistManip::CreateHistYBins(hBins);
name = Form("jetCompEt%sSRMean", component);
this->MakeStringName(name+=fTrigger); // changes 'name'!
hResult->SetName(name);
for(Int_t i = 1; hBins && i <= hBins->GetNbinsY(); ++i){
hResult->SetBinContent(i, hists.At(i-1)->GetMean());
// Stat_t s[11];
// hists.At(i-1)->GetStats(s);//0: S(w), 2: S(wx), 3: S(wx^2)
// Stat_t var = (s[3]-s[2]*s[2]/s[0])/s[0];// cf. Blobel's book, page 135
// if(var >= 0.) hResult->SetBinError(i, TMath::Sqrt(var));
// else this->Error("DrawJetCompSRMean", "Problem of error of mean(%d)-%f",
// i, hists.At(i-1)->GetMean());
}
if(strncmp(component, "Track", 5) == 0){
xyLabel+=";<p_{t}(track)/p_{t}>";
} else if(strncmp(component, "EmCl", 4) == 0){
xyLabel+=";<p_{t}(em)/p_{t}>";
} else if(strncmp(component, "HadCl", 5) == 0){
xyLabel+=";<p_{t}(had)/p_{t}>";
}
hResult->SetTitle(Form("%s fraction%s", component, xyLabel.Data()));
fHistManager->Clear();
fHistManager->AddHists(&hists);
fHistManager->AddHist(hOverall, 1);
fHistManager->AddHist(hResult, 2);
fHistManager->Draw();
return hResult;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawJetProfile(const char* jetPhiEta)
{
// jet1Phi, jet2Eta + ' '/Pt1/Pt2/.../Eta1/Eta2/...
TString name(jetPhiEta);
const Bool_t isPhi = name.Contains("phi", TString::kIgnoreCase);
name += "Prof";
fHistManager->Clear();
TH2* hist2D = this->GetHist2D(name, "jetProf");
if(!hist2D) return NULL;
fHist2D = GFHistManip::CreateHistFlipXY(hist2D);
delete hist2D; hist2D = NULL;
fHist1D = this->GetHistPtCut();
TH1* hNorm=this->GetHist(Form("jet%cNorm%sHist2D",*(name.Data()+3),name.Data()+7),"jetProf");
if(!hNorm){
this->Warning("DrawJetProfile", "correct normalisation not found!");
hNorm = this->GetHistDs1Jet("");
}
if(!hNorm || !fHist1D) return NULL;
TH1* profile = fMultiDmFitter->Fit(fHist2D, fHist1D);
this->MakeStringName(name+=fTrigger); // changes 'name'!
profile->SetName(name);
GFHistManip::CorrectForBinWidth(profile);
if(isPhi){
profile->SetXTitle("#Delta#phi");
// profile->SetYTitle("1/N #upoint p_{t}/d(#Delta#phi) [GeV]");
profile->SetYTitle("<p_{t}>/#Delta(#Delta#phi) [GeV]");
} else {
profile->SetXTitle("#eta_{HFS}-#eta_{jet}");
profile->SetYTitle("<p_{t}>/#Delta(#eta_{HFS}-#eta_{jet}) [GeV]");
}
profile->SetMarkerStyle(8);
// TH1* hNorm = GFHistManip::CreateProjection(fHist2D, name + "AllProjDm", -1, -1);
TArrayD nDstarErr = fMultiDmFitter->FitSingle(hNorm, NULL);
if(nDstarErr[0]) GFHistManip::Scale(profile, 1./nDstarErr[0]);
else this->Error("DrawJetProfile", "no D*!");
fHistManager->AddHists(fMultiDmFitter->GetHists1D(),0);
fHistManager->AddHist(fHist1D,1);
fHistManager->AddHist(hNorm, 1);
fHistManager->SetHistsOption("e");
fHistManager->AddHist(fHist2D,2);
fHistManager->SetHistsOption("BOX",2);
fHistManager->AddHist(fMultiDmFitter->GetHistSignalBack(), 2);
fHistManager->AddHist(profile,3);
fHistManager->Draw();
this->SetHistPointerNull();
// TString title("N(D*)/binwidth in bins of ");
// this->PrintTable(histNumDstar, title += varName);
return profile;
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmJetMult(const char* histName, Bool_t perBinW)
{
// ptDstarNoJet, dstarInJetPt, dstarInJetEta, nNonDstarJet, nJet
// if(perBinW) result is divided through bin width)
fHistManager->Clear();
fHist2D = static_cast<TH2*>(this->GetHist(histName, "jetMul"));
fHist1D = this->GetHistPtCut();
return this->DrawDmVariable(histName, NULL, NULL, perBinW);
}
//____________________________________
TH1* GFDstarHistsAnalysis::CreateL4EffSR(const char* var, Int_t dstar,Int_t fitFlag)
{
// L4 eff in signal region, if(dstar): test D* finder (in gamma p), else just L4 class
// var: pt, eta, phi, ptSumEt, wGammaP, y, nTrackEv, nTrack
const Int_t oldFlag = fMultiDmFitter->GetFitOrNotFlag();
fMultiDmFitter->SetFitOrNotFlag(fitFlag); // 2: count signal region
TH1* h = this->CreateL4Eff(var, dstar);
fMultiDmFitter->SetFitOrNotFlag(oldFlag);
return h;
}
//____________________________________
TH1* GFDstarHistsAnalysis::CreateL4Eff(const char* var, Int_t dstar)
{
// L4 eff in signal region, if(dstar): test D* finder (in gamma p), else just L4 class
// var: pt, eta, phi, ptSumEt, wGammaP, y, nTrackEv, nTrack
TString name(var);
fHistManager->Clear();
// const Double_t oldStart = fMultiDmFitter->GetStartDmSignal();
// const Double_t oldEnd = fMultiDmFitter->GetEndDmSignal();
// fMultiDmFitter->SetStartDmSignal(0.150501);
// fMultiDmFitter->SetEndDmSignal(0.156499);
TH1* histL4Yes = NULL;
TH1* histL4No = NULL;
TH1* hRef = NULL;
TH1* hRefFit = this->GetHistPtCut();
if(name == "nTrackEv"){
TH1* histYes = this->GetHist(dstar ? "nTrackEvHistL4DsTrig" : "nTrackEvHistL4Trig");
TH1* histNo = this->GetHist(dstar ? "nTrackEvHistNoL4DsTrig" : "nTrackEvHistNoL4Trig");
if(!histYes || ! histNo) return NULL;
name += "L4EffSR";
this->MakeStringName(name+=fTrigger);
// this->DeleteObject(name+"Yes");
histL4Yes = static_cast<TH1*>(histYes->Clone(name+"Yes"));
// this->DeleteObject(name+"No");
histL4No = static_cast<TH1*>(histNo ->Clone(name+"No" ));
hRef = static_cast<TH1*>(histL4Yes->Clone(name+"NumAll"));
hRef->Add(histL4No);
} else {
TH2* histL4Yes2D = static_cast<TH2*>
(this->GetHist(name + (dstar ? "HistL4DsTrig" : "HistL4Trig")));
TH2* histL4No2D = static_cast<TH2*>
(this->GetHist(name + (dstar ? "HistNoL4DsTrig" : "HistNoL4Trig")));
if(!histL4Yes2D || ! histL4No2D) return NULL;
fHistManager->AddHist(histL4Yes2D, 0);
fHistManager->AddHist(histL4No2D, 0);
name += "L4EffSR";
this->MakeStringName(name+=fTrigger);
histL4Yes = fMultiDmFitter->Fit(histL4Yes2D, hRefFit);// NULL, NULL, NULL);
fHistManager->AddHists(fMultiDmFitter->GetHists1D(), 3);
// fMultiDmFitter->GetHists1D()->Delete();
delete fMultiDmFitter->GetHistSignalBack();
if(fMultiDmFitter->GetFitOrNotFlag() == 0){
TH2* hRef2D = static_cast<TH2*>(histL4Yes2D->Clone(name+"2DNumAll"));
hRef2D->Add(histL4No2D);
hRef = fMultiDmFitter->Fit(hRef2D, hRefFit);//, NULL, NULL);
} else {
histL4No = fMultiDmFitter->Fit(histL4No2D);// NULL, NULL, NULL);
hRef = static_cast<TH1*>(histL4Yes->Clone(name+"NumAll"));
hRef->Add(histL4No);
}
fHistManager->AddHists(fMultiDmFitter->GetHists1D(), 4);
// fMultiDmFitter->GetHists1D()->Delete();
delete fMultiDmFitter->GetHistSignalBack();
}
// fMultiDmFitter->SetStartDmSignal(oldStart);
// fMultiDmFitter->SetEndDmSignal(oldEnd);
// TH1* hRef = static_cast<TH1*>(histL4Yes->Clone(name+"NumAll"));
// hRef->Add(histL4No);
TGraphAsymmErrors* grPtr = NULL;
TH1* effHist = NULL;
if(fMultiDmFitter->GetFitOrNotFlag() == 2 || fMultiDmFitter->GetFitOrNotFlag() == 3){//count!
effHist = GFHistManip::CreateRatioBinomErrors(histL4Yes, hRef, NULL, NULL, NULL, &grPtr);
if(!effHist) this->Warning("CreateL4Eff", "Cannot do correct errors from binomial!");
}
if(effHist){
effHist->SetName(name);
} else {
effHist = GFHistManip::CreateRatioGassnerErrors(histL4Yes, hRef);
this->Info("CreateL4Eff", "changed to Blobel errors");
effHist = GFHistManip::CreateRatioBlobelErrors(histL4Yes, hRef);
}
fHistManager->AddHist(histL4Yes, 1);
if(histL4No) fHistManager->AddHist(histL4No, 1);
fHistManager->AddHist(hRef, 1);
fHistManager->AddHist(effHist, 2);
if(grPtr) {
fHistManager->AddObject(grPtr, 2, 0, "P");
// for(Int_t i = 1; i <= effHist->GetNbinsX(); ++i) effHist->SetBinError(i, 1.e-7);
// effHist->SetOption("HIST");
effHist->TAttLine::Copy(*grPtr);
effHist->TAttMarker::Copy(*grPtr);
// } else {
// effHist->SetOption("E");
}
fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var));
TString title("#epsilon_{L4}, ");
title += (dstar ? "(D*)" : "(class)");
effHist->SetTitle(title);
effHist->SetYTitle("#epsilon_{L4}");
fHistManager->Draw();
((title += " from ")+=fPeriod) += " in bins of ";
this->PrintTable(effHist, title += var);
this->TmpGraph(grPtr);
return effHist;
}
//____________________________________
TArrayD GFDstarHistsAnalysis::TotalL4EffSR(Int_t dstar)
{
TH2* histL4Yes2D = static_cast<TH2*>
(this->GetHist(dstar ? "etaHistL4DsTrig" : "etaHistL4Trig"));
TH2* histL4No2D = static_cast<TH2*>
(this->GetHist(dstar ? "etaHistNoL4DsTrig" : "etaHistNoL4Trig"));
Int_t first = histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgStartDmSignal);
if(histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgStartDmSignal + 0.0001)> first)
++first;
Int_t last = histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgEndDmSignal);
if(histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgEndDmSignal - 0.0001) < last) --last;
TString name(Form("%s_py", histL4Yes2D->GetName()));
GFHistManip::MakeUniqueName(name);
TH1* histL4Yes = histL4Yes2D->ProjectionY(name, first, last, "E");
name = Form("%s_py", histL4No2D->GetName());
GFHistManip::MakeUniqueName(name);
TH1* histL4No = histL4No2D->ProjectionY(name, first, last, "E");
Double_t sumYes = 0.;
Double_t sumErr2All = 0.;
Double_t sumNo = 0.;
// Double_t sumErr2No = 0.;
for(Int_t i = 1; i <= histL4Yes->GetNbinsX(); ++i){
sumYes += histL4Yes->GetBinContent(i);
const Double_t errY = histL4Yes->GetBinError(i);
// sumErr2Yes += (errY*errY);
sumErr2All += (errY*errY);
sumNo += histL4No->GetBinContent(i);
const Double_t errN = histL4No->GetBinError(i);
// sumErr2No += (errN*errN);
sumErr2All += (errN*errN);
}
const Double_t sumAll = sumYes + sumNo;
// const Double_t errAll = TMath::Sqrt(sumErr2Yes + sumErr2No);
// const Double_t errYes = TMath::Sqrt(sumErr2Yes);
// effective number of entries:
const Double_t nEffAll = sumErr2All ? sumAll*sumAll/sumErr2All : 0.;
Double_t nFillAll = 0.; // number of filling actions
if(TMath::Abs(nEffAll - sumAll) < 1.e-5){// no weights
nFillAll = sumAll;
} else {
//overestimation, but cannot recover how many we really had in projected range
nFillAll = histL4No2D->GetEntries() + histL4Yes2D->GetEntries();
}
TArrayD result = GFMath::Efficiency(sumYes, sumAll, nEffAll, nFillAll);
TString title("#epsilon_{L4}, SR ");
title += (dstar ? "(D*)" : "(class)");
(title += " from ")+=fPeriod;
this->PrintPlusMinus(result, title);
return result;
}
//____________________________________
TH1* GFDstarHistsAnalysis::CreateL4EffDsClassSR(const char* var)
{
// effciency of L4-D*-finder if L4-class already requested
// var: pt, eta, phi, ptSumEt, wGammaP, y
TString name(var);
fHistManager->Clear();
TH2* histL42D = static_cast<TH2*>(this->GetHist(name + "HistL4Trig"));
TH2* histL4Yes2D = static_cast<TH2*>(this->GetHist(name + "HistL4TrigNotDs"));
histL4Yes2D->Add(histL42D, histL4Yes2D, 1, -1);
fHistManager->AddHist(histL4Yes2D, 0);
fHistManager->AddHist(histL42D, 0);
name += "L4EffDsClassSR";
this->MakeStringName(name+=fTrigger);
Int_t first = histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgStartDmSignal);
if(histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgStartDmSignal + 0.0001)> first)
++first;
Int_t last = histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgEndDmSignal);
if(histL4Yes2D->GetXaxis()->FindBin(GFDstarRunPeriods::fgEndDmSignal - 0.0001) < last) --last;
// this->DeleteObject(name+"Yes");
TH1* histL4Yes = histL4Yes2D->ProjectionY(name+"Yes", first, last, "E");
// this->DeleteObject(name+"All");
TH1* histL4All = histL42D->ProjectionY(name+"All", first, last, "E");
TH1* effHist = static_cast<TH1*>(histL4All->Clone(name));
// if there are other values than '1.',
// ==> look in html-docu of TH1::Divide(...) in 3.02_07 and 3.03_?
effHist->Divide(histL4Yes, effHist, 1., 1., "B");
fHistManager->AddHist(histL4Yes, 1);
fHistManager->AddHist(histL4All, 1);
fHistManager->AddHist(effHist, 2);
fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var));
TString title("#epsilon_{L4}: D* if class, SR ");
effHist->SetTitle(title);
effHist->SetOption("E");
fHistManager->Draw();
((title += " from ")+=fPeriod) += " in bins of ";
this->PrintTable(effHist, title += var);
return effHist;
}
//____________________________________
TH1* GFDstarHistsAnalysis::CreateL4EffFit(const char* var, Int_t fitFlag,
Bool_t dstar)
{
// L4 eff via fit,
// if(dstar): test D* finder (in gamma p), else just L4 class
// var: pt, eta, phi, ptSumEt, wGammaP, y (not nTrackEv!)
TString name(var);
fHistManager->Clear();
TH2* histL4Yes2D = static_cast<TH2*>
(this->GetHist(name+ (dstar ? "HistL4DsTrig" : "HistL4Trig")));
TH2* histL4No2D = static_cast<TH2*>
(this->GetHist(name+ (dstar ? "HistNoL4DsTrig" : "HistNoL4Trig")));
TH2* histL4All2D = static_cast<TH2*>
(histL4No2D->Clone(TString(histL4No2D->GetName()).ReplaceAll("No", "All")));
histL4All2D->Add(histL4Yes2D);
histL4All2D->SetTitle("#epsilon(L4) reference");
fHistManager->AddHist(histL4Yes2D, 0);
fHistManager->AddHist(histL4No2D, 0);
fHistManager->AddHist(histL4All2D, 0);
TH1* histAll = this->GetHistPtCut();
fHistManager->SetHistsYTitle(GFAxisVariables::GetAxisLabel(var));
name += "L4EffFit";
this->MakeStringName(name+=fTrigger);
GFMultiDmFitter fi;
fi.SetFitOrNotFlag(fitFlag);
// fi.SetFitModeBins(1111);
TH1* histL4All = fi.Fit(histL4All2D, histAll);
fHistManager->AddHists(fi.GetHists1D(),1);
TH1* histL4Yes = fi.Fit(histL4Yes2D, histAll);
fHistManager->AddHists(fi.GetHists1D(),2);
fHistManager->SetHistsXTitle("#Delta m");
fHistManager->AddHist(histL4All, 3);
fHistManager->AddHist(histL4Yes, 3);
fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var), 3);
TH1* effHist = static_cast<TH1*>(histL4Yes->Clone(name));
// if there are other values than '1.',
// ==> look in html-docu of TH1::Divide(...) in 3.02_07 and 3.03_?
effHist->Divide(histL4Yes, histL4All, 1., 1., "B");
fHistManager->AddHist(effHist, 4);
fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var), 4);
TString title("#epsilon_{L4}, Fit ");
title += (dstar ? "(D*)" : "(class)");
effHist->SetTitle(title);
effHist->SetOption("E");
fHistManager->Draw();
((title += " from ")+=fPeriod) += " in bins of ";
this->PrintTable(effHist, title += var);
return effHist;
}
//____________________________________
GFHistArray* GFDstarHistsAnalysis::CreateHistArray1D(TH2* hist2D,
const char* varName,
Bool_t noOverUnderFlow,
Bool_t projectX) const
{
TString name(varName);
this->MakeStringName(name+=fTrigger);
// // some clean up, works only if the hists get a name in the 'correct' way...
// Int_t nBins = projectX ? hist2D->GetNbinsY() : hist2D->GetNbinsX();
// for(Long_t i = noOverUnderFlow ? 1 : 0; i <= (noOverUnderFlow ? nBins : nBins+1); ++i){
// this->DeleteObject(name + i);
// }
return GFHistManip::CreateHistArray1D(hist2D, varName, name, noOverUnderFlow, projectX);
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHist(const Char_t* name, const char* dirName)
{
// returns histogram with name 'name' from fFileIn, looking into directory 'dirName'
// (but according to fTrigger it may add sth. to the name...)
// name of actual period is added to hist's name
TString histName(name);
switch(fTrigger){
case -1:
break;
case 83:
histName += "S83";
break;
case 84:
histName += "S84";
break;
default:
this->Warning("GetHist","Trigger %d not supported, taking 'all trigger hist'!",
fTrigger);
}
TDirectory* dir = dirName && fFileIn->Get(dirName)
? static_cast<TDirectory*>(fFileIn->Get(dirName)) : fFileIn;
TH1* hist = static_cast<TH1*>(dir->Get(histName));
if(!hist) {
this->Error("GetHist","No histogram '%s' found%s.", histName.Data(),
(dirName && fFileIn->Get(dirName) ? Form(" in directory '%s'", dirName) : ""));
if(dir != fFileIn){
hist = static_cast<TH1*>(fFileIn->Get(histName));
if(hist) this->Error("GetHist", "But found in file!");
}
}
if(hist){
TString name = hist->GetName();
this->MakeStringName(name);
hist->SetName(name);
}
return hist;
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistPt1Cut()
{
return this->GetHist("dmHist1");
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistPt2Cut()
{
return this->GetHist("dmHist2");
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistPtCut()
{
// gives histPt1/Cut or ptHist2Cut corresponding to trigger!
return this->GetTrigger() == 83 ? this->GetHistPt2Cut()
: this->GetHistPt1Cut();
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistDs1Jet(const char* forwback)
{
// delta m hist for D* and (other) jet: forwback = 'Forw', 'Back' or empty
return this->GetHist(Form("dstar1Jet%s",forwback), "dstarJets");
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistDs1JetRef(const char* forwBack, UInt_t refFlag)
{
// if(dmRefFlag != 0): First fit a reference hist with all params free, then fit with fFitMode
// >=2 [default]: D* inclusive
// 1 : D* + 1 jet inclusive
switch(refFlag){
case 0:
return this->GetHistDs1Jet(forwBack);
case 1:
return this->GetHistDs1Jet("");
case 2:
default:
return this->GetHistPtCut();
}
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2D(const char* variable, const char* dirName)
{
// 'pt', 'eta', 'phi' or 'wGammaP' so far... added y(?)
TString var(variable);
TH2 *h = static_cast<TH2*>(this->GetHist(var+="Hist2D", dirName));
if(!h){
this->Info("GetHist2D", "no hist for %s, try %s", var.Data(), variable);
h = static_cast<TH2*>(this->GetHist(variable, dirName));
}
return h;
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DTrack(const char* variablePart)
{
// variable: 'Pt', 'Theta', 'Phi' or 'Length'
// Part: 'K', 'Pi' or 'Pis'
// example: variablePart = "PtPi"
TString var(variablePart);
return static_cast<TH2*>(this->GetHist("trackHist"+var, "tracks"));
}
//________________________________________________
TH3* GFDstarHistsAnalysis::GetHist3D(const char* variablePair)
{
// yeta, wGammaPeta
TString var(variablePair);
return static_cast<TH3*>(this->GetHist(var+="Hist3D"));
}
//________________________________________________
TH1 *GFDstarHistsAnalysis::GetHist2Diff(Bool_t& flip, const char *v1,const char *v2,
const char * nameAdd, const char* dir)
{
flip = kFALSE;
TH1 *h = this->GetHist(Form("%s%s%s", v1, v2, nameAdd), dir);
if(!h) {
this->Info("GetHist2Diff", "try %s%s instead of %s%s", v2, v1, v1, v2);
h = this->GetHist(Form("%s%s%s", v2, v1, nameAdd), dir);
if(h) flip = kTRUE;
}
return h;
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistD0Pt1Cut()
{
return this->GetHist("D0Hist1", "D0");
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistD0Pt2Cut()
{
return this->GetHist("D0Hist2", "D0");
}
//________________________________________________
TH1* GFDstarHistsAnalysis::GetHistD0PtCut()
{
// gives histPt1/Cut or ptHist2Cut corresponding to trigger!
return this->GetTrigger() == 83 ? this->GetHistD0Pt2Cut()
: this->GetHistD0Pt1Cut();
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DD0(const char* variable)
{
// 'pt', 'eta', 'phi' or 'wGammaP' so far...
TString var(variable);
return static_cast<TH2*>(this->GetHist(var+="Hist2DD0", "D0"));
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DTriggers(const char* var)
{
// L4v, L1ac, L1L4; for S83/S84 maybe followed by Ntrig
return static_cast<TH2*>(this->GetHist(Form("trig%s", var), "triggers"));
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DEventClass()
{
return static_cast<TH2*>(this->GetHist("eventClass"));
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DDoubleTrigPt()
{
return static_cast<TH2*>(this->GetHist("doubleTrigPt"));
}
//________________________________________________
TH2* GFDstarHistsAnalysis::GetHist2DPtDsOverSumEt()
{
return static_cast<TH2*>(this->GetHist("PtDstarOverSumETDm"));
}
//________________________________________________
TH1* GFDstarHistsAnalysis::CreateHistYBins(TH2* hist2D, const char* titleAdd)
{
return GFHistManip::CreateHistYBins(hist2D, titleAdd);
}
//________________________________________________
void GFDstarHistsAnalysis::MakeStringName(TString& name) const
{
// replace all 'un-nice' char (and 'GeV'):
// name.ReplaceAll("#{}()[]*/.,;:<>?+-^! ", 0); doesn't work
// and add name of period to the name
// also make name unique via GFHistManip::MakeUniqueName
name += fPeriod;
name.ReplaceAll("#", 0);
name.ReplaceAll("{", 0);
name.ReplaceAll("}", 0);
name.ReplaceAll("(", 0);
name.ReplaceAll(")", 0);
name.ReplaceAll("[", 0);
name.ReplaceAll("]", 0);
name.ReplaceAll("*", 0);
name.ReplaceAll("/", 0);
name.ReplaceAll(".", 0);
name.ReplaceAll(",", 0);
name.ReplaceAll(";", 0);
name.ReplaceAll(":", 0);
name.ReplaceAll("<", 0);
name.ReplaceAll(">", 0);
name.ReplaceAll("?", 0);
name.ReplaceAll("+", 0);
name.ReplaceAll("-", 0);
name.ReplaceAll("^", 0);
name.ReplaceAll("!", 0);
name.ReplaceAll(" ", 0);
name.ReplaceAll("GeV", 0);
GFHistManip::MakeUniqueName(name);
}
// //________________________________________________
// void GFDstarHistsAnalysis::DeleteObject(const char* name) const
// {
// // look for object with name 'name' in gDirectory and deletes it (after warning!)
// if(fgNoDelete) return;
// TObject* obj = gDirectory->GetList()->FindObject(name);
// // TObject* obj = fFileIn->GetList()->FindObject(name);
// if (obj) {
// this->Warning("DeleteObject", "Deleting object %s!", obj->GetName());
// delete obj;
// }
// }
//________________________________________________
void GFDstarHistsAnalysis::SetHistPointerNull()
{
fHistSignalBack = fHist1D = fHistNumDstar = fHist2D = NULL;
}
//________________________________________________
void GFDstarHistsAnalysis::CorrectForBinWidth(TH1* hist) const
{
// divides all bins content and error through the binwidth
GFHistManip::CorrectForBinWidth(hist);
}
//________________________________________________
void GFDstarHistsAnalysis::Normalise(TH1* hist, Bool_t overUnderFlow) const
{
// normalises hist content to 1 (for 1D may take over/underflow into account!)
if(hist){
Double_t integral = 0.;
if(overUnderFlow){
if(!(hist->InheritsFrom(TH2::Class()) || hist->InheritsFrom(TH3::Class()))) {
Int_t overflowBin = hist->GetNbinsX()+1;
integral = hist->Integral(0, overflowBin);
} else {
this->Warning("Normalise", "Cannot care about over/underflow for TH2/TH3 %s!",
hist->GetName());
}
}
if(!integral) integral = hist->Integral();
if(integral) {
hist->Scale(1./ integral);
// hist->SetMaximum(hist->GetMaximum() / integral);
}
}
}
//________________________________________________
void GFDstarHistsAnalysis::FitD0Hists(GFHistArray* arrayD0Hists)
{
Bool_t fillNumD0 = kTRUE;
Int_t numHists = arrayD0Hists->GetEntries();
if(!fHistNumDstar || fHistNumDstar->GetNbinsX() != numHists){
this->Warning("FitD0Hists",
"No fHistNumDstar or wrong number of bins, not filled!");
fillNumD0 = kFALSE;
}
if(!fD0Fitter) fD0Fitter = new D0Fitter;
for(Int_t i = 1; i <= numHists; ++i){
fD0Fitter->Fit(fFitModeD0, (*arrayD0Hists)[i-1]);
if(fillNumD0){
fHistNumDstar->SetBinContent(i, fD0Fitter->GetND0());
fHistNumDstar->SetBinError(i, fD0Fitter->GetND0Err());
}
}
}
//________________________________________________
void GFDstarHistsAnalysis::FitGaus(GFHistArray* arrayGausHists,
TH1* histMean, TH1* histSigma)
{
// fitting a gaus to hists in arrayGausHists
// if histsMean/histsSigma given and correct number of bins:
// will fill mean and sigma into it
if(histMean && histMean->GetNbinsX() != arrayGausHists->GetEntriesFast()){
this->Warning("FitGaus", "Wrong bin number for mean!");
histMean = NULL;
}
if(histSigma && histSigma->GetNbinsX() != arrayGausHists->GetEntriesFast()){
this->Warning("FitGaus", "Wrong bin number for sigma!");
histSigma = NULL;
}
for(Int_t i = 0; i < arrayGausHists->GetEntriesFast(); ++i){
TH1* h = (*arrayGausHists)[i];
cout << "start gaus fit on " << h->GetName() << endl;
if(h->Integral() < 5) {
cout << "only " << h->Integral() << " entries, omit fit!" << endl;
} else h->Fit("gaus","EMI0Q");
TF1* gaus = h->GetFunction("gaus");
if(gaus){
gaus->ResetBit(1<<9);
if(histMean){
histMean->SetBinContent(i+1, gaus->GetParameter(1));
histMean->SetBinError(i+1, gaus->GetParError(1));
}
if(histSigma){
histSigma->SetBinContent(i+1, gaus->GetParameter(2));
histSigma->SetBinError(i+1, gaus->GetParError(2));
}
}
}
}
//________________________________________________
const char* GFDstarHistsAnalysis::GetSampleName() const
{
switch(this->GetTrigger()){
case 83: return "ET33";//"ETAG33";
case 84: return "ET44";//"ETAG44";
case -1: return "DIS_gammaP";
default:
this->Warning("GetSampleName", "trigger %d is no sample!", this->GetTrigger());
return "Unknown";
}
}
//____________________________________
TH1* GFDstarHistsAnalysis::DrawDmVariable(const char* varName, TH2* hWc2D, TH1* hWc, Bool_t perBinW)
{
// assumes that fHist2D and fHist1D are already filled correctly,
// returns N(D*)/binwidth-hist
if(!fHist2D) {
this->Error("DrawDmVariable", "2D hist missing");
return NULL;
}
fHist2D->SetYTitle(varName);// also for title in single hists!
TH1* histNumDstar = fMultiDmFitter->Fit(fHist2D, fHist1D, hWc2D, hWc);
histNumDstar->SetXTitle(varName);
histNumDstar->SetTitle(perBinW ? (TString("N(D*) / #Delta")+= varName).Data() : "N(D*)");
TString name(varName);
this->MakeStringName(name+=fTrigger); // changes 'name'!
histNumDstar->SetName("numDstar"+name);
fHistSignalBack = fMultiDmFitter->GetHistSignalBack();
// fHistSignalBack->SetXTitle(varName);//not necessary in case of MC!
// fHistSignalBack->SetTitle("S/B");//not necessary in case of MC!
fHistSignalBack->SetName("signalBack"+name);//not necessary in case of MC!
fHistNumDstar = histNumDstar;
if(perBinW) GFHistManip::CorrectForBinWidth(fHistNumDstar);
fHistManager->AddHists(fMultiDmFitter->GetHists1D(),0);
fHistManager->AddHist(fHist1D,1);
fHistManager->SetHistsOption("e");
if(fMultiDmFitter->GetHists1DWc()) {
// we need to scale the hists! FIXME ???
if(hWc) fHistManager->AddHistSame(hWc, 1, 0);
fHistManager->AddHistsSame(fMultiDmFitter->GetHists1DWc(), 0);
}
fHistManager->AddHist(fHist2D,2);
fHistManager->SetHistsOption("BOX",2);
fHistManager->SetHistsXTitle("m(K#pi#pi_{s}) - m(K#pi) [GeV]");
fHistManager->AddHist(fHistSignalBack,2);
fHistManager->AddHist(histNumDstar,3);
fHistManager->Draw();
fHistManager->SetHistsLineStyle(1, 0, 1); //solid line, layer 0, all hists 1
this->SetHistPointerNull();
TString title(Form("%s: N(D*)%s in bins of ", this->GetPeriod(), (perBinW?"/binwidth":"")));
this->PrintTable(histNumDstar, title += varName);
return histNumDstar;
}
//________________________________________________
TH1* GFDstarHistsAnalysis::DrawTrackReally(const char* varName)
{
if(!fHist2D || !fHist1D) return NULL;
fHistManager->Clear();
fHist2D->SetYTitle(varName);// also for title in single hists!
TH1* histNumDstar = fMultiDmFitter->Fit(fHist2D, fHist1D);
histNumDstar->SetXTitle(varName);
histNumDstar->SetTitle("N(D*)");
TString name(varName);
this->MakeStringName(name+=fTrigger); // changes 'name'!
histNumDstar->SetName("numDstar"+name);
fHistSignalBack = fMultiDmFitter->GetHistSignalBack();
fHistSignalBack->SetXTitle(varName);
fHistSignalBack->SetTitle("S/B");
fHistSignalBack->SetName("signalBack"+name);
fHistNumDstar = histNumDstar;
// GFHistManip::CorrectForBinWidth(fHistNumDstar); ?????
fHist2D->SetYTitle(varName);
fHistManager->AddHists(fMultiDmFitter->GetHists1D(),0);
fHistManager->AddHist(fHist1D,1);
fHistManager->SetHistsOption("e");
fHistManager->AddHist(fHist2D,2);
fHistManager->SetHistsOption("BOX",2);
fHistManager->SetHistsXTitle("m(K#pi#pi_{s}) - m(K#pi) [GeV]");
fHistManager->AddHist(fHistSignalBack,2);
fHistManager->AddHist(histNumDstar,3);
fHistManager->Draw();
this->SetHistPointerNull();
return histNumDstar;
}
//____________________________________
TArrayD GFDstarHistsAnalysis::CountSignal(TH1* hist, Bool_t subtractBack) const
{
return fMultiDmFitter->CountSignal(hist, NULL, subtractBack);
}
//____________________________________
void GFDstarHistsAnalysis::PrintTable(const TH1* hist, const char* title) const
{
// formated printing the bin content and error, starting with title
GFHistManip::PrintTable(hist, title);
}
//____________________________________
void GFDstarHistsAnalysis::PrintPlusMinus(const TArrayD& values, const char* title) const
{
// priniting <title>: values[0] +/- values[1]
if(values.GetSize() < 2) {
Error("PrintPlusMinus", "array of values too short: %d", values.GetSize());
return;
}
const char line[]
= "------------------------------------------------------------------------------n";
printf("%s| %s: %10.3g +/- %10.3g n%s", line, title, values.At(0), values.At(1), line);
if(values.GetSize() >= 4){
const char txt[] = "asymmetric errors";
for(size_t i = 0; title && i <= strlen(title) - strlen(txt); ++i) printf(" ");
printf(" %s: %+10.3g %+10.3g n%s", txt, values.At(2), values.At(3), line);
}
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.