// Author: Gero Flucke <mailto:flucke@mail.desy.de>
//____________________________________
// GFDstarAnalyse
// Author: Gero Flucke
// Date: May 31st, 2002
// last update: $Date: 2006/01/11 11:33:27 $
// by: $Author: flucke $
//
// C++
#include <iostream>
//#include <cstdlib> //?
// RooT
#include "TString.h"
#include "TMath.h"
#include "TArrayD.h"
#include "TArrayI.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile2D.h"
#include "TFile.h"
#include "TRandom.h"
#include "H1JetProfile2D.h"
// H1OO
#include "H1Steering/H1SteerManager.h"
#include "H1Skeleton/H1Tree.h"
#include "H1Skeleton/H1EventList.h"
#include "H1UserDstar/H1UserDstarEvent.h"
#include "H1Pointers/H1BytePtr.h"
#include "H1Pointers/H1ShortPtr.h"
#include "H1Pointers/H1IntPtr.h"
#include "H1Pointers/H1FloatPtr.h"
#include "H1Mods/H1PartSelTrack.h"
#include "H1Mods/H1PartSelTrackArrayPtr.h"
#include "H1Mods/H1PartDstar.h"
#include "H1Mods/H1PartDstarArrayPtr.h"
#include "H1Mods/H1PartCand.h"
#include "H1Mods/H1PartCandArrayPtr.h"
#include "H1Mods/H1PartJet.h"
#include "H1Mods/H1PartJetArrayPtr.h"
#include "H1Mods/H1PartEm.h"
#include "H1Mods/H1PartEmArrayPtr.h"
#include "H1Mods/H1InclHfsIterator.h"
#include "H1OOBanks/ODS/H1Dmis.h"
#include "H1OOBanks/ODS/H1DmisBankPtr.h"
#include "H1OOBanks/ODS/H1Yecl.h"
#include "H1OOBanks/ODS/H1YeclBankPtr.h"
#include "H1UserDstar/H1UserDstarAccess.h"
#include "H1PhysUtils/H1DedxNLH.h"
#include "H1Geom/H1Dedx.h"
#include "H1Geom/H1DBManager.h"
#include "H1Tracks/H1CombinedFittedTrack.h"
// #include "H1JetCalibration/H1JetCalibration.h"
// private
#include "GFData/GFDstarAnalyse.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFMath.h"
#include "GFUtils/GFDstarRunPeriods.h"
#include "GFData/GFSteerDstarAnalyse.h"
using namespace std;
ClassImp(GFDstarAnalyse)
const Int_t GFDstarAnalyse::kDmHistNumBins = 30;
const Double_t GFDstarAnalyse::kDmHistLowEdge = 0.1375;
const Double_t GFDstarAnalyse::kDmHistHighEdge = 0.1675;
const Int_t GFDstarAnalyse::kD0HistNumBins = 50;
const Float_t GFDstarAnalyse::kD0HistLowEdge = TDatabasePDG::Instance()
->GetParticle(421)->Mass() -.25;
const Float_t GFDstarAnalyse::kD0HistHighEdge = TDatabasePDG::Instance()
->GetParticle(421)->Mass() +.25;
const Double_t GFDstarAnalyse::kCmEnergySqr = 101200.; // s = 4*E_e*E_p, 98-2000
const Int_t GFDstarAnalyse::kNumEventClass = 32;
const Int_t GFDstarAnalyse::kL4FinderBit = 11; // D*-> D0pi ->Kpi pi
const Int_t GFDstarAnalyse::kNumMaxWarn = 10;
const Int_t GFDstarAnalyse::kNumDefaultHists = 16;
//_____________________________________________________
GFDstarAnalyse::GFDstarAnalyse()
{
// get cuts etc. from steering class GFSteerDstarAnalyse
// nad initialise other pointers
if(!gH1SteerManager) new H1SteerManager;
fSteer = static_cast<GFSteerDstarAnalyse*>
(gH1SteerManager->GetSteer(GFSteerDstarAnalyse::Class()));
fY33Bins = fSteer->GetYBins33();
fEventCutyMinS83 = fY33Bins[0];
fEventCutyMaxS83 = fY33Bins[fY33Bins.GetSize()-1];
fY44Bins = fSteer->GetYBins44();
fEventCutyMinS84 = fY44Bins[0];
fEventCutyMaxS84 = fY44Bins[fY44Bins.GetSize()-1];
fWgammaP33Bins = *fSteer->GetWgammaPBins33();
fWgammaP44Bins = *fSteer->GetWgammaPBins44();
if(fWgammaP44Bins.GetSize() < 2){
fWgammaP44Bins.Set(2);
fWgammaP44Bins[0] = TMath::Sqrt(fEventCutyMinS84*this->GetCmEnergySqr());
fWgammaP44Bins[1] = TMath::Sqrt(fEventCutyMaxS84*this->GetCmEnergySqr());
}
fDsCutPt1 = TMath::Min(fSteer->GetPtBins44()->At(0), fSteer->GetPtBins33()->At(0));
fDsCutPt2 = TMath::Max(fSteer->GetPtBins44()->At(0), fSteer->GetPtBins33()->At(0));
// fJetEtaBins = fSteer->GetJetEtaBins();
// fJet1PtBins = fSteer->GetJet1PtBins();
// fJet2PtBins = fSteer->GetJet2PtBins();
fXgamProfile = 0.6;
fIgnoreDedx = kFALSE;
fIgnoreSumEt = kFALSE;
fIgnoreSumPtKPi = kFALSE;
fDmHists = fDmHistsS83 = fDmHistsS84 = NULL;
fD0CheckHists = fD0CheckHistsS83 = NULL;
fD0Hists = fD0HistsS83 = fD0HistsS84 = NULL;
fTrackHists = fTrackHistsS83 = fTrackHistsS84= NULL;
fSumEtHistsS83 = fSumEtHistsS84 = NULL;
fL4NoHistsS83 = fL4NoHistsS84 = fL4YesHistsS83 = fL4YesHistsS84 = NULL;
fL4YesNotDsHistsS83 = fL4YesNotDsHistsS84 = NULL;
fL4NoDstarHistsS83 = fL4NoDstarHistsS84 = fL4YesDstarHistsS83 = fL4YesDstarHistsS84 = NULL;
fEventClassHists = fEventClassHistsS83 = fEventClassHistsS84 = NULL;
fBgFinderHists = NULL;
fTriggersHists = fOpenCharmFindHists = NULL;
fTrackKindHists = fCentralNotCombinedTrackHists = NULL;
fJet1ProfileHists = fJet1ProfileHistsS83 = fJet1ProfileHistsS84 = NULL;
fJet2ProfileHists = fJet2ProfileHistsS83 = fJet2ProfileHistsS84 = NULL;
fJetCompHists = fJetCompHistsS83 = fJetCompHistsS84 = NULL;
// // fNonDstarJetsHists = fNonDstarJetsHistsS83 = fNonDstarJetsHistsS84 = NULL;
// fJetProfiler.SetOwner();
// fJetProfilerS83.SetOwner();
// fJetProfilerS84.SetOwner();
fTrackClusHists = fTrackClusHistsBar = fTrackClusHistsForw = fTrackClusHistsForwBar
= fTrackClusHistsSpac = NULL;
fTrackErrHists = NULL;
fDstar1JetHists = fDstar1JetHistsS83 = fDstar1JetHistsS84
= fDstarBack1JetHists = fDstarBack1JetHistsS83 = fDstarBack1JetHistsS84
= fDstarForw1JetHists = fDstarForw1JetHistsS83 = fDstarForw1JetHistsS84 = NULL;
fDstarDiJetHists = fDstarDiJetHistsS83 = fDstarDiJetHistsS84 = NULL;
fJetMultHistsS83 = NULL;
fTestHists = NULL;
fHistFile = NULL;
fPeriod = new TString(gH1SteerManager->GetCurrentNameSpace() ?
gH1SteerManager->GetCurrentNameSpace() : "Unknown");
fOpenHat = fOpenMods = kTRUE;
fOpenOds = kFALSE;
fJetFlag = 0;
fList = NULL;
// fList = new H1EventList("badDijetRes");
// fList->AddChainInfo("DSTAR", new H1UserDstarEvent);
}
//_____________________________________________________
GFDstarAnalyse::~GFDstarAnalyse()
{
fAllHistArrays.Delete();
delete fList;
fJetProfiler.Delete();
fJetProfilerS83.Delete();
fJetProfilerS84.Delete();
delete fPeriod;
delete fHistFile;
}
//_____________________________________________________
void GFDstarAnalyse::StartJob(){}
void GFDstarAnalyse::StartRun(Int_t run){}
void GFDstarAnalyse::EndRun(Int_t run){}
//_____________________________________________________
void GFDstarAnalyse::EndJob()
{
//
// no normalisation: have to normalise to N(D*) in the end!
// for(Int_t i = 0; i < fJetProfiler.GetEntriesFast(); ++i){
// static_cast<H1JetProfile2D*>(fJetProfiler[i])->NormaliseHistograms();
// }
// for(Int_t i = 0; i < fJetProfilerS83.GetEntriesFast(); ++i){
// static_cast<H1JetProfile2D*>(fJetProfilerS83[i])->NormaliseHistograms();
// }
// for(Int_t i = 0; i < fJetProfilerS84.GetEntriesFast(); ++i){
// static_cast<H1JetProfile2D*>(fJetProfilerS84[i])->NormaliseHistograms();
// }
if(fHistFile && fHistFile->IsOpen()) {
if(fList) {
fHistFile->cd();
fList->Write();
cout << "nWritten event list " << fList->GetName() << " to file "
<< fHistFile->GetName()
<< endl;
}
fHistFile->Write();
// fHistFile->Close(); // not before all histograms are deleted
cout << "nWritten histograms to file " << fHistFile->GetName()
<< endl;
} else {
this->Warning("EndJob","No file open!");
}
}
//_____________________________________________________
void GFDstarAnalyse::PrepareEvent()
{
// do things that should manipulate the event and should stay for the whole event
// e.g. calibrate jets!
// nothing done currently...
return;
/*
static H1PartJetArrayPtr inclJets("InclusiveKtJets");
TObjArray unCalibPartHFS;
H1InclHfsIterator hfs;
for(Int_t i = 0; i < hfs.GetEntries(); ++i) {
unCalibPartHFS.Add(hfs[i]);
}
TObjArray outputJetsCalib;
H1JetCalibration::Instance()->GetHadroo2KtJetCalibration
(inclJets.GetArray(), &outputJetsCalib, &unCalibPartHFS);
for(Int_t i = 0; i < inclJets.GetEntries(); ++i){
TLorentzVector uncalib = inclJets[i]->GetFourVector();
TLorentzVector calib = *static_cast<TLorentzVector*>(outputJetsCalib.At(i));
const_cast<H1PartJet*>(inclJets[i])->SetFourVector(calib.Px(), calib.Py(), calib.Pz());
}
outputJetsCalib.Delete();
static H1PartDstarArrayPtr dstars;
for(Int_t iDs = 0; iDs < dstars.GetEntries(); ++iDs){
const TObjArray* dsJets = H1UserDstarAccess::DstarJets(dstars[iDs]);
if(!dsJets || !dsJets->GetEntriesFast()) continue;
TObjArray unCalibHfsDstar = unCalibPartHFS;
Int_t inHfs = 0; // following casts due to constness:
if(unCalibHfsDstar.Remove((TObject*)dstars[iDs]->GetKaon()->GetParticle())) ++inHfs;
if(unCalibHfsDstar.Remove((TObject*)dstars[iDs]->GetPion()->GetParticle())) ++inHfs;
if(unCalibHfsDstar.Remove((TObject*)dstars[iDs]->GetSlowPion()->GetParticle())) ++inHfs;
if(inHfs != 3){
this->Warning("GFDstarAnalyse", "missed %d D* daughter%s", 3-inHfs, inHfs==2?"":"s");
}
unCalibHfsDstar.Add(dstars[iDs]);
unCalibHfsDstar.Compress();
H1JetCalibration::Instance()->GetHadroo2KtJetCalibration
(dsJets, &outputJetsCalib, &unCalibHfsDstar);
for(Int_t i = 0; i < dsJets->GetEntriesFast(); ++i){
H1PartJet* jet = static_cast<H1PartJet*>(dsJets->At(i));
TLorentzVector uncalib = jet->GetFourVector();
TLorentzVector calib = *static_cast<TLorentzVector*>(outputJetsCalib.At(i));
jet->SetFourVector(calib.Px(), calib.Py(), calib.Pz());
}
outputJetsCalib.Delete();
}
*/
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsDm(GFHistArray*& dmHists, GFHistArray*& dmHistsS83,
GFHistArray*& dmHistsS84, const TString& nameAdd)
{
// foreseen is nameAdd = "Wc" or ""
dmHists = this->CreateHists("2D"+nameAdd, "#Deltam vs. <y-axis>", -1, 2);
dmHistsS83 = this->CreateHists("2D"+nameAdd, "#Deltam vs. <y-axis>", 83, 2);
dmHistsS84 = this->CreateHists("2D"+nameAdd, "#Deltam vs. <y-axis>", 84, 2);
TDirectory* oldDir = this->OpenHistsFile();
TString title1D("p_{t}(D*) > ");
dmHists->Add(new TH1F("dmHist1"+nameAdd, (TString(title1D)+=fDsCutPt1)+=" GeV", //NDefHist+0
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
dmHists->Add(new TH1F("dmHist2"+nameAdd, (TString(title1D)+=fDsCutPt2)+=" GeV", // +1
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
// dmHists->Add(new TH1F("mDstarHist"+nameAdd, "m(K#pi#pi_{s})", // +2
// 2*kDmHistNumBins, 1.75, 2.25));
dmHistsS83->Add(new TH1F("dmHist1"+nameAdd+"S83", (TString(title1D)+=fDsCutPt1)+=" GeV",//+0
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
dmHistsS83->Add(new TH1F("dmHist2"+nameAdd+"S83", (TString(title1D)+=fDsCutPt2)+=" GeV",// +1
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
// dmHistsS83->Add(new TH1F("mDstarHist"+nameAdd+"S83", "m(K#pi#pi_{s})", // +2
// 2*kDmHistNumBins, 1.75, 2.25));
dmHistsS84->Add(new TH1F("dmHist1"+nameAdd+"S84", (TString(title1D)+=fDsCutPt1)+=" GeV",// +0
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
dmHistsS84->Add(new TH1F("dmHist2"+nameAdd+"S84", (TString(title1D)+=fDsCutPt2)+=" GeV",// +1
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
// dmHistsS84->Add(new TH1F("mDstarHist"+nameAdd+"S84", "m(K#pi#pi_{s})", // +2
// 2*kDmHistNumBins, 1.75, 2.25));
Double_t help[kDmHistNumBins+1];
for(Int_t i = 0; i <= kDmHistNumBins; ++i){
help[i] = i*(kDmHistHighEdge-kDmHistLowEdge)/kDmHistNumBins + kDmHistLowEdge;
}
dmHistsS83->Add(new TH3F("yetaHist3D"+nameAdd+"S83","#Deltam vs. y vs #eta (S83)",//nDefHi+2
kDmHistNumBins, help,
fY33Bins.GetSize()-1, fY33Bins.GetArray(),
fSteer->GetEtaBins33()->GetSize()-1,
fSteer->GetEtaBins33()->GetArray()));
dmHistsS84->Add(new TH3F("yetaHist3D"+nameAdd+"S84","#Deltam vs. y vs #eta (S84)",// +2
kDmHistNumBins, help,
fY33Bins.GetSize()-1, fY33Bins.GetArray(),
fSteer->GetEtaBins44()->GetSize()-1,
fSteer->GetEtaBins44()->GetArray()));
dmHistsS83->Add(new TH3F("wGammaPetaHist3D"+nameAdd+"S83",
"#Deltam vs. W_{#gammap} vs #eta (S83)",//+3
kDmHistNumBins, help,
fWgammaP33Bins.GetSize()-1,fWgammaP33Bins.GetArray(),
fSteer->GetEtaBins33()->GetSize()-1,
fSteer->GetEtaBins33()->GetArray()));
dmHistsS84->Add(new TH3F("wGammaPetaHist3D"+nameAdd+"S84",
"#Deltam vs. W_{#gammap} vs #eta (S84)",//+3
kDmHistNumBins, help,
fWgammaP44Bins.GetSize()-1,fWgammaP44Bins.GetArray(),
fSteer->GetEtaBins44()->GetSize()-1,
fSteer->GetEtaBins44()->GetArray()));
oldDir->cd();
fAllHistArrays.Add(dmHists);
fAllHistArrays.Add(dmHistsS83);
fAllHistArrays.Add(dmHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsD0()
{
TDirectory* oldDir = this->MakeHistDir("D0", "D0 histograms");
fD0Hists = GFHistManip::ClearCreateArray(fD0Hists);
fD0HistsS83 = GFHistManip::ClearCreateArray(fD0HistsS83);
fD0HistsS84 = GFHistManip::ClearCreateArray(fD0HistsS84);
TString title1D("p_{t}(D*) > ");
fD0Hists->Add(new TH1F("D0Hist1", (TString(title1D)+=fDsCutPt1)+=" GeV", // 0
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge));
fD0Hists->Add(new TH1F("D0Hist2", (title1D+=fDsCutPt2)+=" GeV", // 1
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge));
fD0Hists->Add(new TH2F("ptHist2DD0", "m(K#pi) vs. p_{t}(D*)", // 2
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetPtBins44()->GetSize()-1,
fSteer->GetPtBins44()->GetArray()));
fD0Hists->Add(new TH2F("etaHist2DD0", "m(K#pi) vs. #eta(D*)", // 3
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetEtaBins33()->GetSize()-1,
fSteer->GetEtaBins33()->GetArray()));
fD0Hists->Add(new TH2F("phiHist2DD0", "m(K#pi) vs. #phi(D*)", // 5
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetPhiBins()->GetSize()-1, fSteer->GetPhiBins()->GetArray()));
GFHistManip::CreateAnalogHists(fD0Hists, fD0HistsS83, "S83");
GFHistManip::CreateAnalogHists(fD0Hists, fD0HistsS84, "S84");
// ugly:
delete (*fD0HistsS83)[2];
fD0HistsS83->AddAt(new TH2F("ptHist2DD0S83", "m(K#pi) vs. p_{t}(D*) S83",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetPtBins33()->GetSize()-1,
fSteer->GetPtBins33()->GetArray()), 2);
delete (*fD0HistsS84)[3];
fD0HistsS84->AddAt(new TH2F("etaHist2DD0S84", "m(K#pi) vs. #eta(D*) S84",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetEtaBins44()->GetSize()-1,
fSteer->GetEtaBins44()->GetArray()), 3);
fD0HistsS83->Add(new TH2F("wGammaPHist2DD0S83","m(K#pi) vs. W_{#gammap} (S83)",// 5
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fWgammaP33Bins.GetSize()-1,
fWgammaP33Bins.GetArray()));
fD0HistsS84->Add(new TH2F("wGammaPHist2DD0S84","m(K#pi) vs. W_{#gammap} (S84)",// 5
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fWgammaP44Bins.GetSize()-1,
fWgammaP44Bins.GetArray()));
fD0HistsS83->Add(new TH2F("yHist2DD0S83","m(K#pi) vs. y (S83)",// 6
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fY33Bins.GetSize()-1, fY33Bins.GetArray()));
fD0HistsS84->Add(new TH2F("yHist2DD0S84","m(K#pi) vs. y (S84)",// 6
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fY44Bins.GetSize()-1, fY44Bins.GetArray()));
oldDir->cd();
fAllHistArrays.Add(fD0Hists);
fAllHistArrays.Add(fD0HistsS83);
fAllHistArrays.Add(fD0HistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTracks()
{
TDirectory* oldDir = this->MakeHistDir("tracks", "K, pi, pi_s quantities");
fTrackHists = GFHistManip::ClearCreateArray(fTrackHists);
const Float_t minTheta = 20. * TMath::DegToRad();
const Float_t maxTheta = 160. * TMath::DegToRad();
const Int_t nBins = 24; //can be downsized to 12,8,6,4,3,2,1!
fTrackHists->Add(new TH2F("trackHistPtK", "p_{t} (K)", //0
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, fSteer->GetPtCutK(), 4.+fSteer->GetPtCutK()));
fTrackHists->Add(new TH2F("trackHistThetaK", "#theta (K)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, minTheta, maxTheta));
fTrackHists->Add(new TH2F("trackHistPhiK", "#phi (K)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, -TMath::Pi(), TMath::Pi()));
fTrackHists->Add(new TH2F("trackHistLengthK", "track length (K)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, fSteer->GetLengthCutKPi(), 66));
fTrackHists->Add(new TH2F("trackHistPtPi", "p_{t} (#pi)", //4
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, fSteer->GetPtCutPi(), 4.+fSteer->GetPtCutPi()));
fTrackHists->Add(new TH2F("trackHistThetaPi", "#theta (#pi)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, minTheta, maxTheta));
fTrackHists->Add(new TH2F("trackHistPhiPi", "#phi (#pi)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, -TMath::Pi(), TMath::Pi()));
fTrackHists->Add(new TH2F("trackHistLengthPi", "track length (#pi)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, fSteer->GetLengthCutKPi(), 66));
// const Double_t ptPisBins[nBins+1] = {0.12, 0.13, 0.14,
// 0.15, 0.16, 0.17,
// 0.18, 0.19, 0.20,
// 0.21, 0.22, 0.23,
// 0.24, 0.26, 0.28,
// 0.3, 0.33, 0.365,
// 0.4, 0.43, 0.465,
// 0.5, 0.565, 0.63, 0.7};
const Double_t ptPisBins[nBins+1] = {0.12, 0.13, 0.14,
0.15, 0.161, 0.173,
0.186, 0.20, 0.215,
0.231, 0.248, 0.266,// here
0.285, 0.305, 0.327,
0.351, 0.377, 0.405,
0.435, 0.467, 0.501,
0.537, 0.577, 0.627, 0.7};
fTrackHists->Add(new TH2F("trackHistPtPis", "p_{t} (#pi_{s})", //8
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, ptPisBins));//0.12, 1.08));
fTrackHists->Add(new TH2F("trackHistThetaPis", "#theta (#pi_{s})",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, minTheta, maxTheta));
fTrackHists->Add(new TH2F("trackHistPhiPis", "#phi (#pi_{s})",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, -TMath::Pi(), TMath::Pi()));
fTrackHists->Add(new TH2F("trackHistLengthPis", "track length (#pi_{s})",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, fSteer->GetLengthCutPis(), 66));
const Double_t rStartLow = 20.;
const Double_t rStartHigh = 45.;
const Float_t dedxLimit = 4.5;
const Int_t dedxN = 200;
const Int_t dedxNbinsP = 50;
Double_t dedxBinsMom[dedxNbinsP+1];
Double_t dedxBinsMomPis[dedxNbinsP+1];
if(!GFMath::EquidistLogBins(dedxBinsMom, dedxNbinsP,
TMath::Min(fSteer->GetPtCutK(), fSteer->GetPtCutPi()), 5.)
|| !GFMath::EquidistLogBins(dedxBinsMomPis, dedxNbinsP, 0.12, 5.)){
this->Error("CreateHistsTracks", "Problem with dedx bins for momentum");
}
fTrackHists->Add(new TH2F("trackHistStartK", "r_{start}(K)", //12
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, rStartLow, rStartHigh));
fTrackHists->Add(new TH2F("trackHistStartPi", "r_{start}(#pi)", //13
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, rStartLow, rStartHigh));
fTrackHists->Add(new TH2F("trackHistStartPis", "r_{start}(#pi_{s})", //14
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, rStartLow, rStartHigh));
fTrackHists->Add(new TH2F("trackHistDcaK", "dca' (K)", //15
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -1., 1.));
fTrackHists->Add(new TH2F("trackHistDcaPi", "dca' (#pi)", //16
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -1., 1.));
fTrackHists->Add(new TH2F("trackHistDcaPis", "dca' (#pi_{s})", //17
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -1.5, 1.5));
// #ifdef H1OO24
fTrackHists->Add(new TH2F("trackHistNhitK", "N_{hit}^{CJC} (K)", //18
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
fTrackHists->Add(new TH2F("trackHistNhitPi", "N_{hit}^{CJC} (#pi)", //19
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
fTrackHists->Add(new TH2F("trackHistNhitPis", "N_{hit}^{CJC} (#pi_{s})", //20
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
// #else
// fTrackHists->Add(new TH2F("trackHistNhitK", "N_{hit}^{dEdx} (K)", //18
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
// fTrackHists->Add(new TH2F("trackHistNhitPi", "N_{hit}^{dEdx} (#pi)", //19
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
// fTrackHists->Add(new TH2F("trackHistNhitPis", "N_{hit}^{dEdx} (#pi_{s})", //20
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, 10., 58.));
// #endif
fTrackHists->Add(new TH2F("trackHistDz0K", "dz_{0}' (K)", //21
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -10., 10.));
fTrackHists->Add(new TH2F("trackHistDz0Pi", "dz_{0}' (#pi)", //22
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -10., 10.));
fTrackHists->Add(new TH2F("trackHistDz0Pis", "dz_{0}' (#pi_{s})", //23
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBins, -10., 10.));
const Int_t nBinsLH = 110; // pi_s often 2*110...
fTrackHists->Add(new TH2F("trackHistLHK", "LH to be K (K)", //24
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBinsLH, .0, 1.1));
fTrackHists->Add(new TH2F("trackHistLHPi", "LH to be #pi (#pi)", //25
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBinsLH, .0, 1.1));
fTrackHists->Add(new TH2F("trackHistLHPis", "LH to be #pi #pi_{s})", //26
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 2*nBinsLH,.0,1.1));
fTrackHists->Add(new TH2F("trackHistDedxK", "dE/dx vs p (K)", //27
dedxNbinsP, dedxBinsMom, dedxN, 0.7, dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxPi", "dE/dx vs p (#pi)", //28
dedxNbinsP, dedxBinsMom, dedxN, 0.7,dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxPis", "dE/dx vs p (#pi_{s})", //29
dedxNbinsP, dedxBinsMomPis, dedxN, 0.7, dedxLimit+3.));
fTrackHists->Add(new TH2F("trackHistDedxLHKK", "dE/dx vs p (K) LH(K) OK", //30
dedxNbinsP, dedxBinsMom, dedxN, 0.7, dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxLHPiPi", "dE/dx vs p (#pi) LH(#pi) OK", //31
dedxNbinsP, dedxBinsMom, dedxN, 0.7,dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxLHPiPis", "dE/dx vs p (#pi_{s}) LH(#pi) OK", //32
dedxNbinsP, dedxBinsMomPis, dedxN, 0.7, dedxLimit+3.));
fTrackHists->Add(new TH2F("trackHistDedxNoLHKK", "dE/dx vs p (K) LH(K) not OK", //33
dedxNbinsP, dedxBinsMom, dedxN, 0.7, dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxNoLHPiPi", "dE/dx vs p (#pi) LH(#pi) not OK", //34
dedxNbinsP, dedxBinsMom, dedxN, 0.7, dedxLimit));
fTrackHists->Add(new TH2F("trackHistDedxNoLHPiPis", "dE/dx vs p (#pi_{s}) LH(#pi) not OK",
dedxNbinsP, dedxBinsMomPis, dedxN, 0.7, dedxLimit+3.)); // 35
const Int_t nBinsDedxMin = 40;
fTrackHists->Add(new TH2F("trackHistDedxOKmExpK", "(dE/dx-Param) vs p (K) LH(K) OK",
dedxNbinsP, dedxBinsMom, nBinsDedxMin, -.8, .8)); //36
fTrackHists->Add(new TH2F("trackHistDedxOKmExpPi", "(dE/dx-Param) vs p (#pi) LH(#pi) OK",
dedxNbinsP, dedxBinsMom, nBinsDedxMin, -.8, .8)); // 37
fTrackHists->Add(new TH2F("trackHistDedxOKmExpPis","(dE/dx-Param) vs p (#pi_{s}) LH(K) OK",
dedxNbinsP, dedxBinsMomPis, nBinsDedxMin, -.8, .8)); //38
fTrackHists->Add(new TH2F("trackHistDedxNotOKmExpK", "(dE/dx-Param) vs p (K) LH(K) not OK",
dedxNbinsP, dedxBinsMom, nBinsDedxMin, -.8, .8)); // 39
fTrackHists->Add(new TH2F("trackHistDedxNotOKmExpPi", // 40
"(dE/dx-Param) vs p (#pi) LH(#pi) not OK",
dedxNbinsP, dedxBinsMom, nBinsDedxMin, -.8, .8));
fTrackHists->Add(new TH2F("trackHistDedxNotOKmExpPis", // 41
"(dE/dx-Param) vs p (#pi_{s}) LH(#pi) not OK",
dedxNbinsP, dedxBinsMomPis, nBinsDedxMin, -.8, .8));
fTrackHists->Add(new TH2F("trackHistNLHK", "NLH to be K (K)", //42
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBinsLH, 0., 1.1));
fTrackHists->Add(new TH2F("trackHistNLHPi", "NLH to be #pi (#pi)", //43
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBinsLH, 0., 1.1));
fTrackHists->Add(new TH2F("trackHistNLHPis", "NLH to be #pi #pi_{s})", //44
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, nBinsLH, 0., 1.1));
Double_t pBins[9] = {fSteer->GetPtCutK(), 0.5, .8, 1., 1.3, 1.6, 2., 4., 8.};
Double_t dmBins[kDmHistNumBins+1];
GFMath::EquidistBins(dmBins, kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge);
Double_t lhBins[nBinsLH+1];
GFMath::EquidistBins(lhBins, nBinsLH, 0., 1.1);
fTrackHists->Add(new TH3F("trackHistLHKvsP", "LH to be K (K) vs p vs #Deltam", //45
kDmHistNumBins, dmBins, 8, pBins, nBinsLH, lhBins));
fTrackHists->Add(new TH3F("trackHistLHPivsP", "LH to be #pi (#pi) vs p vs #Deltam", //46
kDmHistNumBins, dmBins, 8, pBins, nBinsLH, lhBins));
const Double_t pBinsPis[6] = {0.12, 0.2, 0.3, 0.5, .8, 1.2};
fTrackHists->Add(new TH3F("trackHistLHPisvsP", "LH to be #pi (#pi_{s}) vs p vs #Deltam",//47
kDmHistNumBins,dmBins,5,pBinsPis,nBinsLH,lhBins));//220, -.0, 1.1));
// FIXME: bins for ET44??
fTrackHists->Add(new TH3F("trackHistLHKvsPtDs", "LH to be K (K) vs p_{t}(D*) vs #Deltam",//48
kDmHistNumBins, dmBins, fSteer->GetPtBins33()->GetSize()-2,
fSteer->GetPtBins33()->GetArray(), nBinsLH, lhBins));
fTrackHists->Add(new TH3F("trackHistLHPivsPtDs", "LH to be #pi (#pi) vs p_{t}(D*) vs #Deltam",
kDmHistNumBins, dmBins, fSteer->GetPtBins33()->GetSize()-2, //49
fSteer->GetPtBins33()->GetArray(), nBinsLH, lhBins));
// nothing for pi_s here!
// all tracks hists:
fTrackHists->Add(new TH2F("trackHistDedxAll", "dE/dx vs p, D* in SR ", //50
dedxNbinsP, dedxBinsMomPis, dedxN, 0.7, dedxLimit+3.));
fTrackHists->Add(new TH2F("trackHistDedxPis2", "dE/dx vs p, #pi_{s}, D* in SR ", //51
dedxNbinsP, dedxBinsMomPis, dedxN, 0.7, dedxLimit+3.));
GFHistManip::CreateAnalogHists(fTrackHists, fTrackHistsS83, "S83");
GFHistManip::CreateAnalogHists(fTrackHists, fTrackHistsS84, "S84");
oldDir->cd();
fAllHistArrays.Add(fTrackHists);
fAllHistArrays.Add(fTrackHistsS83);
fAllHistArrays.Add(fTrackHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsSumEt()
{
TDirectory* oldDir = this->OpenHistsFile();
fSumEtHistsS83 = GFHistManip::ClearCreateArray(fSumEtHistsS83);
fSumEtHistsS84 = GFHistManip::ClearCreateArray(fSumEtHistsS84);
fSumEtHistsS83->Add(new TH1F("sumET", "#sum(E_{t}^{#theta>10 degree})", 40, 0, 40));
fSumEtHistsS83->Add(new TH2F("sumETDm", "#sum(E_{t}^{#theta>10 degree})",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
10, 0, 40));
fSumEtHistsS83->Add(new TH1F("PtDstarOverSumET",
"p_{t}(D*) / #sum (E_{t}^{#theta > 10})", 50, 0, .5));
fSumEtHistsS83->Add(new TH2F("PtDstarOverSumETDm",
"p_{t}(D*) / #sum (E_{t}^{#theta > 10})",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
10, 0, .5));
GFHistManip::CreateAnalogHists(fSumEtHistsS83, fSumEtHistsS84, "S84");
GFHistManip::AddToHistsName("S83", fSumEtHistsS83);
oldDir->cd();
fAllHistArrays.Add(fSumEtHistsS83);
fAllHistArrays.Add(fSumEtHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTriggers()
{
TDirectory* oldDir = this->MakeHistDir("triggers", "fired subtriggers");
fTriggersHists = GFHistManip::ClearCreateArray(fTriggersHists);
// all D*
fTriggersHists->Add(new TH2F("trigL4v", "L4 verified ST, all D*",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1ac", "L1 actual ST, all D*",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1L4", "L1 actual AND L4 verified, all D*",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
// S83 visible range (via y_jb and anti DIS)
fTriggersHists->Add(new TH2F("trigL4vNtrigS83", "L4 verified ST, D* of S83 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1acNtrigS83", "L1 actual ST, D* of S83 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1L4NtrigS83",
"L1 actual AND L4 verified, D* of S83 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
// IsS83
fTriggersHists->Add(new TH2F("trigL4vS83", "L4 verified ST, D* of S83",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1acS83", "L1 actual ST, D* of S83",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1L4S83", "L1 actual AND L4 verified, D* of S83",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
// S84 visible range (via y_jb and anti DIS)
fTriggersHists->Add(new TH2F("trigL4vNtrigS84", "L4 verified ST, D* of S84 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1acNtrigS84", "L1 actual ST, D* of S84 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1L4NtrigS84",
"L1 actual AND L4 verified, D* of S84 visible",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
// IsS84
fTriggersHists->Add(new TH2F("trigL4vS84", "L4 verified ST, D* of S84",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1acS84", "L1 actual ST, D* of S84",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
fTriggersHists->Add(new TH2F("trigL1L4S84", "L1 actual AND L4 verified, D* of S84",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
128, 0, 128));
oldDir->cd();
fAllHistArrays.Add(fTriggersHists);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsEventClass()
{
TDirectory* oldDir = this->OpenHistsFile();
fEventClassHists = GFHistManip::ClearCreateArray(fEventClassHists);
fEventClassHistsS83 = GFHistManip::ClearCreateArray(fEventClassHistsS83);
fEventClassHistsS84 = GFHistManip::ClearCreateArray(fEventClassHistsS84);
// Take care:
// Same event (weight) will be filled in several times if it belongs to several classes!
fEventClassHists->Add(new TH3F("eventClass", "event classes", // 0
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
kNumEventClass, 0, kNumEventClass,
50, 0, 50)); // weights
fEventClassHists->Add(new TH2F("eventClassNo15", "event classes if not 15", // 1
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
kNumEventClass, 0, kNumEventClass));
fEventClassHists->Add(new TH2F("eventClassNotL4", "event classes, not L4 found", // 2
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
kNumEventClass, 0, kNumEventClass));
fEventClassHists->Add(new TH2F("eventWeight", "event weight", //3
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
50, 0, 50)); // weights
GFHistManip::CreateAnalogHists(fEventClassHists, fEventClassHistsS83, "S83");
GFHistManip::CreateAnalogHists(fEventClassHists, fEventClassHistsS84, "S84");
// fEventClassHists->Add(new TH2F("eventClassS83", "event classes if(S83)",
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// kEventClassBits, 0, kEventClassBits));
// fEventClassHists->Add(new TH2F("eventClassS84", "event classes if(S84)",
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// kEventClassBits, 0, kEventClassBits));
fAllHistArrays.Add(fEventClassHists);
fAllHistArrays.Add(fEventClassHistsS83);
fAllHistArrays.Add(fEventClassHistsS84);
fOpenCharmFindHists = GFHistManip::ClearCreateArray(fOpenCharmFindHists);
fOpenCharmFindHists->Add(new TH2F("openCharmL4tgP", "L4 finder: tagged #gammap",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
32, 0, 32));
fOpenCharmFindHists->Add(new TH2F("openCharmL4tgPS83", "L4 finder: tagged #gammap (S83)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
32, 0, 32));
fOpenCharmFindHists->Add(new TH2F("openCharmL4tgPS84", "L4 finder: tagged #gammap (S84)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
32, 0, 32));
oldDir->cd();
fAllHistArrays.Add(fOpenCharmFindHists);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTrackKind()
{
TDirectory* oldDir = this->OpenHistsFile();
fTrackKindHists = GFHistManip::ClearCreateArray(fTrackKindHists);
fTrackKindHists->Add(new TH1F("3centTracks", "all 3 tracks central",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("2cent1combTrack", "2 central, 1 combined track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("1cent2combTrack", "1 central, 2 combined track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("2cent1forwTrack", "2 central, 1 forward track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("1cent2forwTrack", "1 central, 2 forward track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("3combTracks", "all 3 tracks combined",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("2comb1forwTrack", "2 combined, 1 forward track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("1comb2forwTrack", "1 combined, 2 forward track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("3forwTracks", "all 3 tracks forward",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fTrackKindHists->Add(new TH1F("1cent1comb1forwTrack",
"1 central, 1 forward, 1 combined track",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
oldDir->cd();
fAllHistArrays.Add(fTrackKindHists);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsCentralNotCombinedTrack()
{
if(!gH1Tree || NULL == gH1Tree->GetChain(H1Tree::kOdsTree)){
this->Warning("CreateHistsCentralNotCombinedTrack",
"H1Tree not yet existing or no ODS tree, ignore!");
return;
}
TDirectory* oldDir = this->OpenHistsFile();
fCentralNotCombinedTrackHists
= GFHistManip::ClearCreateArray(fCentralNotCombinedTrackHists);
fCentralNotCombinedTrackHists->
Add(new TH1F("usualTracksNoCut","tracks as usual, no cuts", //0
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fCentralNotCombinedTrackHists->
Add(new TH1F("centralNotCombinedNoCut", //1
"combined track replaced by central, no cuts",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fCentralNotCombinedTrackHists->
Add(new TH1F("usualTracksCutted","tracks as usual, default cuts",//2
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fCentralNotCombinedTrackHists->
Add(new TH1F("centralNotCombinedCutted",
"combined track replaced by central, default cuts",//3
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
oldDir->cd();
fAllHistArrays.Add(fCentralNotCombinedTrackHists);
}
//_____________________________________________________
void GFDstarAnalyse::CreateRunDepAnalog(TH1* input, GFHistArray* output, Int_t trigger)
{
if(trigger != 84){
this->Warning("CreateRunDepAnalog", "trigger %d not valid, only 84 implemented!", trigger);
return;
}
TString newName(input->GetName());
Size_t samplePosition = newName.Index("S84");
if(samplePosition != kNPOS) newName.ReplaceAll("S84", "_");
TString newTitle(input->GetTitle()); newTitle+= ", run ";
GFDstarRunPeriods periods(trigger, fPeriod->Data());
for(Int_t i = 0; i < periods.GetNumPeriods(); ++i){
TH1* newHist = static_cast<TH1*>(input->Clone(newName+((Long_t)i)));//runPeriods84[i-1]));
newHist->SetTitle((newTitle + ((Long_t)periods.GetFirstRun(i)))
+ " - " + ((Long_t)periods.GetLastRun(i)));
if(samplePosition != kNPOS) {
newHist->SetName(TString(newHist->GetName()) += "S84");
}
output->Add(newHist);
}
}
//_____________________________________________________
GFHistArray* GFDstarAnalyse::CreateHists(const char* nameAdd,
const char* title,
Int_t subTrig, UInt_t dimension, TDirectory* dir)
{
// create a GFHistArray of hists with identifier 'nameAdd', 'title' and
// for 'subTrig' (if subtrig == -1: subTrig is not put to name,
// binning is mixed!)
// 'dimension' is one or two so far:
// 1: TH1 with x-axis for pt, eta, phi, ptSumEt, wGammaP, y, nJet, ptJet1/2, etaJet1/2, xGamma
// 2: x-axis is for delta m, y axes are for pt, eta, phi, ptSumEt, wGammaP, etc.
// if(dir != NULL) dir->cd() before creating hists
// which hists are created has to be in sync with 'FillHists' and ...Mc::FillHistsGen!
const TArrayD* wGammaPBins = NULL;
const TArrayD* ptBins = NULL;
const TArrayD* etaBins = NULL;
const TArrayD* yBins = NULL;
if(subTrig == 83){
ptBins = fSteer->GetPtBins33();
etaBins = fSteer->GetEtaBins33();
// phi is same!
wGammaPBins = &fWgammaP33Bins;
yBins = &fY33Bins;
} else if(subTrig == 84){
ptBins = fSteer->GetPtBins44();
etaBins = fSteer->GetEtaBins44();
// phi is same!
wGammaPBins = &fWgammaP44Bins;
yBins = &fY44Bins;
} else if(subTrig == -1){
ptBins = fSteer->GetPtBins44(); // pt from 44: wider range
etaBins = fSteer->GetEtaBins33(); // eta from 33: ----"----
// phi is same!
TArrayD w(9);
TArrayD y(w.GetSize());
for(Int_t i = 0; i < w.GetSize(); ++i){
y[i] = fY44Bins[0] + i*(fY33Bins[fY33Bins.GetSize()-1]-fY44Bins[0])/ (w.GetSize()-1);
w[i] = fWgammaP44Bins[0] + i*(fWgammaP33Bins[fWgammaP33Bins.GetSize()-1]
-fWgammaP44Bins[0])/(w.GetSize()-1);
}
wGammaPBins = new TArrayD(w);
yBins = new TArrayD(y);
} else {
this->Warning("CreateHists",
"No histograms to build for subTrig = %d, ignore!", subTrig);
return NULL;
}
TString name("Hist");
name += nameAdd;
if(subTrig >= 0) (name += "S") += subTrig;
GFHistArray* result = NULL;
GFHistManip::ClearCreateArray(result);
TDirectory* oldDir = NULL;
if(dir) {
oldDir = gDirectory;
dir->cd();
} else {
oldDir = this->OpenHistsFile();
}
const Int_t nNumBinsPtSumEt = 10;
const Float_t maxPtSumEt = 0.5;
const Int_t maxNumJets = 5;
const Float_t maxMinZvtx = 55.;
// for double diff:
Double_t eta[] = {etaBins->At(0), etaBins->At(0)+0.6, 0.2, etaBins->At(etaBins->GetSize()-1)};
// Double_t pt[] = {ptBins->At(0), 3., 4.25, 7.};
Double_t pt[] = {ptBins->At(0), 3., 4.5, 8.5};
// Double_t eta[] = {-1.5, -0.9, -0.2, 0.6, etaBins->At(etaBins->GetSize()-1)};
// Double_t pt[] = {ptBins->At(0), 4., 7.};
const Double_t w0 = wGammaPBins->At(0), w1 = wGammaPBins->At(wGammaPBins->GetSize()-1);
Double_t w[] = {w0, (w1-w0)/3.+w0, 2.*(w1-w0)/3.+w0, w1};
Double_t zDs[] = {0., 0.2, 0.45, 0.75};
//0, 0.1, 0.2, 0.3, 0.45, 0.6, 0.8
if(dimension == 1){
result->Add(new TH1F("pt"+name, title, //0
ptBins->GetSize()-1, ptBins->GetArray()));
result->Add(new TH1F("eta"+name, title, //1
etaBins->GetSize()-1, etaBins->GetArray()));
result->Add(new TH1F("phi"+name, title, //2
fSteer->GetPhiBins()->GetSize()-1,
fSteer->GetPhiBins()->GetArray()));
result->Add(new TH1F("ptSumEt"+name, title, nNumBinsPtSumEt, 0., maxPtSumEt));//3
result->Add(new TH1F("wGammaP"+name, title, //4
wGammaPBins->GetSize()-1, wGammaPBins->GetArray()));
result->Add(new TH1F("y"+name, title, //5
yBins->GetSize()-1, yBins->GetArray()));
result->Add(new TH1F("zVtx"+name, title, 110, -maxMinZvtx, maxMinZvtx)); //6
result->Add(new TH1F("ptKPi"+name, title,
nNumBinsPtSumEt, fSteer->GetPtCutK() + fSteer->GetPtCutPi(),5.35));//7
result->Add(new TH1F("nTrack"+name, title, 8, 3., 27.)); //8
result->Add(new TH1F("nJet"+name, title, maxNumJets, 0, maxNumJets)); //9
result->Add(new TH1F("zDs"+name, title, //10
fSteer->GetZDsBins()->GetSize()-1, fSteer->GetZDsBins()->GetArray()));
result->Add(new TH2F("etapt"+name, title, sizeof(eta)/sizeof(eta[0])-1, eta,
sizeof(pt)/sizeof(pt[0])-1, pt)); //11
result->Add(new TH2F("wGammaPpt"+name, title, sizeof(w)/sizeof(w[0])-1, w,
sizeof(pt)/sizeof(pt[0])-1, pt)); //12
result->Add(new TH2F("wGammaPeta"+name, title, sizeof(w)/sizeof(w[0])-1, w,
sizeof(eta)/sizeof(eta[0])-1, eta)); //13
result->Add(new TH2F("zDseta"+name, title, sizeof(zDs)/sizeof(zDs[0])-1, zDs,
sizeof(eta)/sizeof(eta[0])-1, eta)); //14
result->Add(new TH2F("zDspt"+name, title, sizeof(zDs)/sizeof(zDs[0])-1, zDs,
sizeof(pt)/sizeof(pt[0])-1, pt)); //14
//// result->Add(new TH1F("EmPz"+name, title, 9, 0., 45.)); //10
// result->Add(new TH1F("ptJet1"+name, title,
// fSteer->GetJet1PtBins()->GetSize()-1,
// fSteer->GetJet1PtBins()->GetArray())); //10
// result->Add(new TH1F("ptJet2"+name, title,
// fSteer->GetJet2PtBins()->GetSize()-1,
// fSteer->GetJet2PtBins()->GetArray())); //11
// result->Add(new TH1F("etaJet1"+name, title,
// fSteer->GetJetEtaBins()->GetSize()-1,
// fSteer->GetJetEtaBins()->GetArray())); //12
// result->Add(new TH1F("etaJet2"+name, title,
// fSteer->GetJetEtaBins()->GetSize()-1,
// fSteer->GetJetEtaBins()->GetArray())); //13
// result->Add(new TH1F("xGamma"+name, title, xGamBins->GetSize()-1, xGamBins->GetArray()));//14
// not more than 'kNumDefaultHists'!
} else if(dimension == 2){
result->Add(new TH2F("pt"+name, title, //0
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBins->GetSize()-1, ptBins->GetArray()));
result->Add(new TH2F("eta"+name, title, //1
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBins->GetSize()-1, etaBins->GetArray()));
result->Add(new TH2F("phi"+name, title, //2
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
fSteer->GetPhiBins()->GetSize()-1,
fSteer->GetPhiBins()->GetArray()));
result->Add(new TH2F("ptSumEt"+name, title, //3
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nNumBinsPtSumEt, 0., maxPtSumEt));
result->Add(new TH2F("wGammaP"+name, title, //4
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
wGammaPBins->GetSize()-1, wGammaPBins->GetArray()));
result->Add(new TH2F("y"+name, title, //5
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
yBins->GetSize()-1, yBins->GetArray()));
result->Add(new TH2F("zVtx"+name, title, //6
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
20, -maxMinZvtx, maxMinZvtx));
result->Add(new TH2F("ptKPi"+name, title, // 7
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nNumBinsPtSumEt, fSteer->GetPtCutK() + fSteer->GetPtCutPi(), 5.35));
result->Add(new TH2F("nTrack"+name, title,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 8, 3., 27.)); //8
result->Add(new TH2F("nJet"+name, title, //9
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
maxNumJets, 0., maxNumJets));
result->Add(new TH2F("zDs"+name, title, //10
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
fSteer->GetZDsBins()->GetSize()-1, fSteer->GetZDsBins()->GetArray()));
Double_t dmBins[kDmHistNumBins+1];
GFMath::EquidistBins(dmBins, kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge);
result->Add(new TH3F("etapt"+name, title, kDmHistNumBins, dmBins, //11
sizeof(eta)/sizeof(eta[0])-1, eta, sizeof(pt)/sizeof(pt[0])-1, pt));
result->Add(new TH3F("wGammaPpt"+name, title, kDmHistNumBins, dmBins,
sizeof(w)/sizeof(w[0])-1, w, sizeof(pt)/sizeof(pt[0])-1, pt)); //12
result->Add(new TH3F("wGammaPeta"+name, title, kDmHistNumBins, dmBins,
sizeof(w)/sizeof(w[0])-1, w, sizeof(eta)/sizeof(eta[0])-1, eta)); //13
result->Add(new TH3F("zDseta"+name, title, kDmHistNumBins, dmBins,
sizeof(zDs)/sizeof(zDs[0])-1, zDs,
sizeof(eta)/sizeof(eta[0])-1, eta)); //14
result->Add(new TH3F("zDspt"+name, title, kDmHistNumBins, dmBins,
sizeof(zDs)/sizeof(zDs[0])-1, zDs,
sizeof(pt)/sizeof(pt[0])-1, pt)); //15
// // result->Add(new TH2F("EmPz"+name, title,
// // kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// // 9, 0., 45.)); //10
// result->Add(new TH2F("ptJet1"+name, title, //10
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// fSteer->GetJet1PtBins()->GetSize()-1,
// fSteer->GetJet1PtBins()->GetArray()));
// result->Add(new TH2F("ptJet2"+name, title, //11
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// fSteer->GetJet2PtBins()->GetSize()-1,
// fSteer->GetJet2PtBins()->GetArray()));
// result->Add(new TH2F("etaJet1"+name, title, //12
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// fSteer->GetJetEtaBins()->GetSize()-1,
// fSteer->GetJetEtaBins()->GetArray()));
// result->Add(new TH2F("etaJet2"+name, title, //13
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// fSteer->GetJetEtaBins()->GetSize()-1,
// fSteer->GetJetEtaBins()->GetArray()));
// result->Add(new TH2F("xGamma"+name, title, //14
// kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
// xGamBins->GetSize()-1, xGamBins->GetArray()));
// not more than 'kNumDefaultHists'!
} else {
this->Warning("CreateHists",
"No histograms in %d dimensions, ignore!", dimension);
delete result;
result = NULL;
}
if(result && kNumDefaultHists != result->GetEntriesFast()){
this->Fatal("CreateHists", "number of hists %d, should be %d", result->GetEntriesFast(),
kNumDefaultHists);
}
oldDir->cd();
if(subTrig == -1){
delete wGammaPBins;
delete yBins;
}
return result;
}
//_____________________________________________________
void GFDstarAnalyse::FillHists(GFHistArray* hists, const H1PartDstar* ds,
Int_t subTrig, Double_t weightExtern)
{
// filling histograms in hists with weight by GetWeight(subTrig):
// 1st with pt, 2nd with eta, 3rd with phi of 'ds'
// 4th with pt(D*)/ sum(E_t)
// 5th with wgammap for corresponding 'subTrig'
// (83 or 84, else 5th hist ignored!).
// If a hist is 2D: x-axis will be filled with ds->GetDm().
// if weightExtern == -1. (default), use this->GetWeight(subTrig) as weight, else weightExtern
// which hists are expected has to be in sync with 'CreateHists'!
const Double_t weight = (weightExtern == -1. ? this->GetWeight(subTrig, ds) : weightExtern);
const Double_t dm = ds->GetDm();
// stuff for all triggers: pt,eta, phi
TH1* hist = (*hists)[0];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, ds->GetPt(), weight);
else hist->Fill(ds->GetPt(), weight);
hist = (*hists)[1];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, ds->GetEta(), weight);
else hist->Fill(ds->GetEta(), weight);
hist = (*hists)[2];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, ds->GetPhi(), weight);
else hist->Fill(ds->GetPhi(), weight);
const Double_t sumEt = this->GetSumEt();
if(sumEt){
hist = (*hists)[3];
if(hist->InheritsFrom(TH2::Class()))
static_cast<TH2*>(hist)->Fill(dm, ds->GetPt()/sumEt, weight);
else hist->Fill(ds->GetPt()/sumEt, weight);
} else { // if no user tree, no wonder !
this->Warning("FillHists", "sumEt = 0, run %d evt %d",
gH1Tree->GetRunNumber(), gH1Tree->GetEventNumber());
}
// now looking for wGammaP and y:
const Double_t wGammaP = this->GetWgammaP(subTrig);
hist = (*hists)[4];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, wGammaP, weight);
else hist->Fill(wGammaP, weight);
const Double_t y = this->GetY(subTrig);
hist = (*hists)[5];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, y, weight);
else hist->Fill(y, weight);
static H1FloatPtr zVtx("VtxZ");
hist = (*hists)[6];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, *zVtx, weight);
else hist->Fill(*zVtx, weight);
const Double_t sumPtKPi = ds->GetKaon()->GetPt() + ds->GetPion()->GetPt();
hist = (*hists)[7];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, sumPtKPi, weight);
else hist->Fill(sumPtKPi, weight);
static H1ShortPtr numTr("NumCentralTracks");
hist = (*hists)[8];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, *numTr, weight);
else hist->Fill(*numTr, weight);
// now with jets:
TObjArray jets;
this->GetNonDstarJets(ds, jets);
const Int_t nJets = jets.GetEntriesFast();
// this->RemoveLowPt(&jets, 2);
hist = (*hists)[9];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, nJets, weight);
else hist->Fill(nJets, weight);
const Double_t zDs = this->GetZ(ds);
hist = (*hists)[10];
if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, zDs,weight);
else hist->Fill(this->GetZ(ds), weight);
hist = (*hists)[11];
if(hist->InheritsFrom(TH3::Class())) {
static_cast<TH3*>(hist)->Fill(dm, ds->GetEta(), ds->GetPt(), weight);
} else {
static_cast<TH2*>(hist)->Fill(ds->GetEta(), ds->GetPt(), weight);
}
hist = (*hists)[12];
if(hist->InheritsFrom(TH3::Class())) {
static_cast<TH3*>(hist)->Fill(dm, wGammaP, ds->GetPt(), weight);
} else {
static_cast<TH2*>(hist)->Fill(wGammaP, ds->GetPt(), weight);
}
hist = (*hists)[13];
if(hist->InheritsFrom(TH3::Class())) {
static_cast<TH3*>(hist)->Fill(dm, wGammaP, ds->GetEta(), weight);
} else {
static_cast<TH2*>(hist)->Fill(wGammaP, ds->GetEta(), weight);
}
hist = (*hists)[14];
if(hist->InheritsFrom(TH3::Class())) {
static_cast<TH3*>(hist)->Fill(dm, zDs, ds->GetEta(), weight);
} else {
static_cast<TH2*>(hist)->Fill(zDs, ds->GetEta(), weight);
}
hist = (*hists)[15];
if(hist->InheritsFrom(TH3::Class())) {
static_cast<TH3*>(hist)->Fill(dm, zDs, ds->GetPt(), weight);
} else {
static_cast<TH2*>(hist)->Fill(zDs, ds->GetPt(), weight);
}
// // static H1FloatPtr q2e("Q2e");
// // if(*q2e){
// // static H1FloatPtr empz("Epz");
// // hist = (*hists)[10];
// // if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, *empz, weight);
// // else hist->Fill(*empz, weight);
// // }
// if(nJets > 0){ // same in ...Mc::FillHistsGen!
// H1PartJet* jet1 = static_cast<H1PartJet*>(jets[0]);
// hist = (*hists)[10];
// if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, jet1->GetPt(),
// weight);
// else hist->Fill(jet1->GetPt(), weight);
// hist = (*hists)[11;
// if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, jet1->GetEta(),
// weight);
// else hist->Fill(jet1->GetEta(), weight);
// const H1PartJet* jetHighPt = static_cast<H1PartJet*>(jets[this->GetHighestPt(jets)]);
// hist = (*hists)[14];
// const Double_t xGamma = this->GetXgammaMass(ds, jetHighPt, this->GetYjb());
// if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, xGamma, weight);
// else hist->Fill(xGamma, weight);
// }
// if(nJets > 1){
// H1PartJet* jet2 = static_cast<H1PartJet*>(jets[1]);
// hist = (*hists)[11;
// if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, jet2->GetPt(),
// weight);
// else hist->Fill(jet2->GetPt(), weight);
// hist = (*hists)[13];
// if(hist->InheritsFrom(TH2::Class())) static_cast<TH2*>(hist)->Fill(dm, jet2->GetEta(),
// weight);
// else hist->Fill(jet2->GetEta(), weight);
// }
// not more than 'kNumDefaultHists'!
}
//_____________________________________________________
Bool_t GFDstarAnalyse::FillHistsDstar(const H1PartDstar* dstar)
{
// returns true if dstar survives all cuts (only looser pt and eta cuts neccesary)
// (subtrigger don't care for return value!)
if(!fDmHists) this->CreateHistsDm(fDmHists, fDmHistsS83, fDmHistsS84, "");
return this->FillHistsDstar(dstar, fDmHists, fDmHistsS83, fDmHistsS84, kTRUE); // L4 required!
}
//_____________________________________________________
Bool_t GFDstarAnalyse::FillHistsDstar(const H1PartDstar* dstar,
GFHistArray* dmHists,
GFHistArray* dmHistsS83,
GFHistArray* dmHistsS84, Bool_t requireL4)
{
// returns true if dstar survives all cuts (only looser pt and eta cuts neccesary)
// (subtrigger don't care for return value!)
// first argument may be NULL, if the others also ==> crash!
if(!this->DstarCutLoose(dstar)) return kFALSE;
const Double_t dm = dstar->GetDm();
const Double_t etaDs = dstar->GetEta();
const Bool_t s83 = (requireL4 ? this->IsS83(dstar) : this->IsS83NoL4Check(dstar));
const Bool_t s84 = (requireL4 ? this->IsS84(dstar) : this->IsS84NoL4Check(dstar));
const Double_t weight83 = this->GetWeight(83, dstar);
const Double_t weight84 = this->GetWeight(84, dstar);
if(dmHists){
this->FillHists(dmHists, dstar, -1);
(*dmHists)[kNumDefaultHists + 0]->Fill(dm);
}
const Double_t origPtCut = fSteer->GetPtBins33()->At(0);
const_cast<TArrayD*>(fSteer->GetPtBins33())->AddAt(0., 0);
if((requireL4 ? this->IsS83(dstar) : this->IsS83NoL4Check(dstar))){
(*dmHistsS83)[kNumDefaultHists + 0]->Fill(dm, weight83);
}
const_cast<TArrayD*>(fSteer->GetPtBins33())->AddAt(origPtCut, 0);
if(dstar->GetPt() > fDsCutPt2) {
if(dmHists) (*dmHists)[kNumDefaultHists + 1]->Fill(dm);
if(s84) (*dmHistsS84)[kNumDefaultHists + 1]->Fill(dm, weight84);
}
if(s84) {
this->FillHists(dmHistsS84, dstar, 84);
(*dmHistsS84)[kNumDefaultHists + 0]->Fill(dm, weight84);
static_cast<TH3*>((*dmHistsS84)[kNumDefaultHists + 2])->Fill(dm, this->GetY44(),
etaDs, weight84);
static_cast<TH3*>((*dmHistsS84)[kNumDefaultHists + 3])->Fill(dm, this->GetWgammaP44(),
etaDs, weight84);
}
if(s83) {
this->FillHists(dmHistsS83, dstar, 83);
(*dmHistsS83)[kNumDefaultHists + 1]->Fill(dm, weight83);
static_cast<TH3*>((*dmHistsS83)[kNumDefaultHists + 2])->Fill(dm, this->GetY33(),
etaDs, weight83);
static_cast<TH3*>((*dmHistsS83)[kNumDefaultHists + 3])->Fill(dm, this->GetWgammaP33(),
etaDs, weight83);
}
return kTRUE;
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTracks(const H1PartDstar* dstar)
{
// all track quantities are filled!
if(!fTrackHists) this->CreateHistsTracks();
this->FillHistsTracks(dstar, fTrackHists, this->GetWeight(-1, dstar));
if(this->IsS83(dstar)) this->FillHistsTracks(dstar, fTrackHistsS83,this->GetWeight(83,dstar));
if(this->IsS84(dstar)) this->FillHistsTracks(dstar, fTrackHistsS84,this->GetWeight(84,dstar));
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTracks(const H1PartDstar* dstar,
GFHistArray* trackHists, Double_t w)// const
{
// fill with weight 'w' (except dedx scatters... )
const Bool_t oldIgnoreDedx = this->SetIgnoreDedx(kFALSE);
const Bool_t kDedxOK = this->IsDedxOK(dstar, 0);
const Bool_t piDedxOK = this->IsDedxOK(dstar, 1);
const Bool_t pisDedxOK = this->IsDedxOK(dstar, 2);
this->SetIgnoreDedx(oldIgnoreDedx);
const Double_t dm = dstar->GetDm();
const Int_t numDstar = H1UserDstarAccess::ArrayIndex(dstar);
if(oldIgnoreDedx || (kDedxOK && piDedxOK && pisDedxOK)){
static_cast<TH2*>((*trackHists)[0])->Fill(dm, dstar->GetKaon()->GetPt(), w);
static_cast<TH2*>((*trackHists)[1])->Fill(dm, dstar->GetKaon()->GetTheta(), w);
static_cast<TH2*>((*trackHists)[2])->Fill(dm, dstar->GetKaon()->GetPhi(), w);
static_cast<TH2*>((*trackHists)[3])->Fill(dm, dstar->GetKaon()->GetRadLength(), w);
static_cast<TH2*>((*trackHists)[4])->Fill(dm, dstar->GetPion()->GetPt(), w);
static_cast<TH2*>((*trackHists)[5])->Fill(dm, dstar->GetPion()->GetTheta(), w);
static_cast<TH2*>((*trackHists)[6])->Fill(dm, dstar->GetPion()->GetPhi(), w);
static_cast<TH2*>((*trackHists)[7])->Fill(dm, dstar->GetPion()->GetRadLength(), w);
static_cast<TH2*>((*trackHists)[8])->Fill(dm, dstar->GetSlowPion()->GetPt(), w);
static_cast<TH2*>((*trackHists)[9])->Fill(dm, dstar->GetSlowPion()->GetTheta(), w);
static_cast<TH2*>((*trackHists)[10])->Fill(dm, dstar->GetSlowPion()->GetPhi(), w);
static_cast<TH2*>((*trackHists)[11])->Fill(dm, dstar->GetSlowPion()->GetRadLength(), w);
static_cast<TH2*>((*trackHists)[12])->Fill(dm, H1UserDstarAccess::TrackStartK(numDstar), w);
static_cast<TH2*>((*trackHists)[13])->Fill(dm, H1UserDstarAccess::TrackStartPi(numDstar),w);
static_cast<TH2*>((*trackHists)[14])->Fill(dm,H1UserDstarAccess::TrackStartPis(numDstar),w);
static_cast<TH2*>((*trackHists)[15])->Fill(dm, dstar->GetKaon()->GetDcaPrime(), w);
static_cast<TH2*>((*trackHists)[16])->Fill(dm, dstar->GetPion()->GetDcaPrime(), w);
static_cast<TH2*>((*trackHists)[17])->Fill(dm, dstar->GetSlowPion()->GetDcaPrime(), w);
#ifdef H1OO24
(*trackHists)[18]->Fill(dm, dstar->GetKaon()->GetNHitCjc());
(*trackHists)[19]->Fill(dm, dstar->GetPion()->GetNHitCjc());
(*trackHists)[20]->Fill(dm, dstar->GetSlowPion()->GetNHitCjc());
#else
// static_cast<TH2*>((*trackHists)[18])->Fill(dm, dstar->GetKaon()->GetNHitDedx(), w);
// static_cast<TH2*>((*trackHists)[19])->Fill(dm, dstar->GetPion()->GetNHitDedx(), w);
// static_cast<TH2*>((*trackHists)[20])->Fill(dm, dstar->GetSlowPion()->GetNHitDedx(), w);
static_cast<TH2*>((*trackHists)[18])->Fill(dm,H1UserDstarAccess::TrackNhitCjcK(numDstar),w);
static_cast<TH2*>((*trackHists)[19])->
Fill(dm, H1UserDstarAccess::TrackNhitCjcPi(numDstar), w);
static_cast<TH2*>((*trackHists)[20])->
Fill(dm, H1UserDstarAccess::TrackNhitCjcPis(numDstar), w);
#endif
static_cast<TH2*>((*trackHists)[21])->Fill(dm, H1UserDstarAccess::TrackDz0K(numDstar), w);
static_cast<TH2*>((*trackHists)[22])->Fill(dm, H1UserDstarAccess::TrackDz0Pi(numDstar), w);
static_cast<TH2*>((*trackHists)[23])->Fill(dm, H1UserDstarAccess::TrackDz0Pis(numDstar), w);
}
const Double_t momK = dstar->GetKaon()->GetP();
const Double_t momPi = dstar->GetPion()->GetP();
const Double_t momPis = dstar->GetSlowPion()->GetP();
// likelihoods filled if all others are OK, but user tree values outdated!
// const Float_t lhKOld = H1UserDstarAccess::TrackLHKofK(numDstar);
// const Float_t lhPiOld = H1UserDstarAccess::TrackLHPiofPi(numDstar);
// const Float_t lhPisOld = H1UserDstarAccess::TrackLHPiofPis(numDstar);
const Float_t lhK = dstar->GetKaon()->DedxLikelihood(H1DedxNLH::kKaon);
const Float_t lhPi = dstar->GetPion()->DedxLikelihood(H1DedxNLH::kPion);
const Float_t lhPis = dstar->GetSlowPion()->DedxLikelihood(H1DedxNLH::kPion);
// static Bool_t first = kTRUE;
// if(first){
// this->Info("FillHistsTracks", "manipulated LH plots");
// first = kFALSE;
// }
if(oldIgnoreDedx || (piDedxOK && pisDedxOK)){
Float_t intLhK = (lhK < 0. ? 1.09 : lhK);
// if(intLhK >= fSteer->GetLikeliHoodCut()){
static_cast<TH2*>((*trackHists)[24])->Fill(dm, intLhK, w);
static_cast<TH2*>((*trackHists)[42])->Fill(dm, this->GetNLH(dstar, 0), w);
static_cast<TH3*>((*trackHists)[45])->Fill(dm, momK, intLhK, w);
static_cast<TH3*>((*trackHists)[48])->Fill(dm, dstar->GetPt(), intLhK, w);
// }
}
if(oldIgnoreDedx || (kDedxOK && pisDedxOK)){
Float_t intLhPi = (lhPi < 0. ? 1.09 : lhPi);
// if(intLhPi >= fSteer->GetLikeliHoodCut()){
static_cast<TH2*>((*trackHists)[25])->Fill(dm, intLhPi, w);
static_cast<TH2*>((*trackHists)[43])->Fill(dm, this->GetNLH(dstar, 1), w);
static_cast<TH3*>((*trackHists)[46])->Fill(dm, momPi, intLhPi, w);
static_cast<TH3*>((*trackHists)[49])->Fill(dm, dstar->GetPt(), intLhPi, w);
// }
}
if(oldIgnoreDedx || (kDedxOK && piDedxOK)){
Float_t intLhPis = (lhPis < 0. ? 1.09 : lhPis);
static_cast<TH2*>((*trackHists)[26])->Fill(dm, intLhPis, w);
static_cast<TH2*>((*trackHists)[44])->Fill(dm, this->GetNLH(dstar, 2), w);
static_cast<TH3*>((*trackHists)[47])->Fill(dm, momPis, intLhPis, w);
// no pt(D*) plot here for pi_s
}
const H1DedxNLH dedxNLH;
if(oldIgnoreDedx || (piDedxOK && pisDedxOK)){ // kaon
const Float_t dedxK = dstar->GetKaon()->GetDedx();
(*trackHists)[27]->Fill(momK, dedxK);
if(kDedxOK){
(*trackHists)[33]->Fill(momK, dedxK);
(*trackHists)[36]->Fill(momK, dedxK - dedxNLH.RefDedxKaon(momK));
} else {
(*trackHists)[30]->Fill(momK, dedxK);
(*trackHists)[39]->Fill(momK, dedxK - dedxNLH.RefDedxKaon(momK));
}
}
if(oldIgnoreDedx || (kDedxOK && pisDedxOK)){ // pion
const Float_t dedxPi = dstar->GetPion()->GetDedx();
(*trackHists)[28]->Fill(momPi, dedxPi);
if(piDedxOK){
(*trackHists)[34]->Fill(momPi, dedxPi);
(*trackHists)[37]->Fill(momPi, dedxPi - dedxNLH.RefDedxPion(momPi));
} else {
(*trackHists)[31]->Fill(momPi, dedxPi);
(*trackHists)[40]->Fill(momPi, dedxPi - dedxNLH.RefDedxPion(momPi));
}
}
if(oldIgnoreDedx || (kDedxOK && piDedxOK)){ // slow pion
const Float_t dedxPis = dstar->GetSlowPion()->GetDedx();
(*trackHists)[29]->Fill(momPis, dedxPis);
if(pisDedxOK){
(*trackHists)[35]->Fill(momPis, dedxPis);
(*trackHists)[38]->Fill(momPis, dedxPis - dedxNLH.RefDedxPion(momPis));
} else {
(*trackHists)[32]->Fill(momPis, dedxPis);
(*trackHists)[41]->Fill(momPis, dedxPis - dedxNLH.RefDedxPion(momPis));
}
}
if(!this->IsSignalRegionDm(dstar)) return;
(*trackHists)[51]->Fill(momPis, dstar->GetSlowPion()->GetDedx());
static H1PartSelTrackArrayPtr tracks;
for(Int_t i = 0; i < tracks.GetEntries(); ++i){
if(tracks[i]->IsCentralTrk() && tracks[i]->GetDedx() > 0.){
// fill only those far away from pion hypothesis, others randomly downscaled:
const Float_t dedx = tracks[i]->GetDedx();
const Double_t p = tracks[i]->GetP();
Float_t diff = TMath::Abs(dedx - dedxNLH.RefDedxPion(p));
if(diff < 0.1) diff = 0.1;
if(diff>=1. || gRandom->Rndm() < diff) (*trackHists)[50]->Fill(p,dedx);
else if(p < 0.17){ // take also additional slow particles:
const Double_t lowEdge = (*trackHists)[50]->GetBinLowEdge(1);
if((p - lowEdge) < (0.17-lowEdge)*gRandom->Rndm()*gRandom->Rndm()){
(*trackHists)[50]->Fill(p,dedx);
}
}
}
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsD0(const H1PartDstar* dstar)
{
if(!fD0Hists) this->CreateHistsD0();
if(this->D0CutLoose(dstar)){
const Bool_t s83 = this->IsS83(dstar);
const Bool_t s84 = this->IsS84(dstar);
const Double_t w83 = this->GetWeight(83, dstar);
const Double_t w84 = this->GetWeight(84, dstar);
const Double_t dsPt = dstar->GetPt();
const Double_t etaDs = dstar->GetEta();
const Double_t phiDs = dstar->GetPhi();
const Double_t mD0 = dstar->GetD0Mass();
(*fD0Hists)[0]->Fill(mD0);
(*fD0Hists)[2]->Fill(mD0, dsPt);
(*fD0Hists)[3]->Fill(mD0, etaDs);
(*fD0Hists)[4]->Fill(mD0, phiDs);
// fixme: same gymnastics for fDsCutPt2 and fDsCutPt1 as for dstar hists?
if(s83) (*fD0HistsS83)[0]->Fill(mD0, w83);
if(dsPt > fDsCutPt2){
(*fD0Hists)[1]->Fill(mD0);
if(s84) (*fD0HistsS84)[1]->Fill(mD0, w84);
}
if(s84) {
(*fD0HistsS84)[0]->Fill(mD0, w84);
static_cast<TH2*>((*fD0HistsS84)[2])->Fill(mD0, dsPt, w84);
static_cast<TH2*>((*fD0HistsS84)[3])->Fill(mD0, etaDs, w84);
static_cast<TH2*>((*fD0HistsS84)[4])->Fill(mD0, phiDs, w84);
static_cast<TH2*>((*fD0HistsS84)[5])->Fill(mD0, this->GetWgammaP44(), w84);
static_cast<TH2*>((*fD0HistsS84)[6])->Fill(mD0, this->GetY44(), w84);
}
if(s83) {
(*fD0HistsS83)[1]->Fill(mD0, w83);
static_cast<TH2*>((*fD0HistsS83)[2])->Fill(mD0, dsPt, w83);
static_cast<TH2*>((*fD0HistsS83)[3])->Fill(mD0, etaDs, w83);
static_cast<TH2*>((*fD0HistsS83)[4])->Fill(mD0, phiDs, w83);
static_cast<TH2*>((*fD0HistsS83)[5])->Fill(mD0, this->GetWgammaP33(), w83);
static_cast<TH2*>((*fD0HistsS83)[6])->Fill(mD0, this->GetY33(), w83);
}
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTriggers(const H1PartDstar* dstar)
{
// assuming dstar fulfills cuts that you want
// only pt2 cut for S83 and eta cut for s84!
if(!fTriggersHists) this->CreateHistsTriggers();
const Bool_t s83 = this->IsS83(dstar);
const Bool_t s84 = this->IsS84(dstar);
const Float_t yjb = this->GetYjb();
static H1FloatPtr elecE("ElecE");
const Bool_t s83visible = (*elecE < 10. // remove 'nice' DIS
&& yjb >= fEventCutyMinS83 && yjb <= fEventCutyMaxS83
&& dstar->GetPt() >= fSteer->GetPtBins33()->At(0));
const Bool_t s84visible = (*elecE < 10. // remove 'nice' DIS
&& yjb >= fEventCutyMinS84 && yjb <= fEventCutyMaxS84
&& dstar->GetEta() >= fSteer->GetEtaBins44()->At(0));
const Double_t dm = dstar->GetDm();
static H1BytePtr l4vst("Il4vst");
static H1BytePtr l1ac("Il1ac");
for(Int_t i = 0; i < 128; ++i){
if(l4vst[i] > 0){
(*fTriggersHists)[0]->Fill(dm, i);
if(s83visible) (*fTriggersHists)[3]->Fill(dm, i);
if(s83) (*fTriggersHists)[6]->Fill(dm, i);
if(s84visible) (*fTriggersHists)[9]->Fill(dm, i);
if(s84) (*fTriggersHists)[12]->Fill(dm, i);
}
if(l1ac[i] > 0){
(*fTriggersHists)[1]->Fill(dm, i);
if(s83visible) (*fTriggersHists)[4]->Fill(dm, i);
if(s83) (*fTriggersHists)[7]->Fill(dm, i);
if(s84visible) (*fTriggersHists)[10]->Fill(dm, i);
if(s84) (*fTriggersHists)[13]->Fill(dm, i);
if(l4vst[i] > 0){
(*fTriggersHists)[2]->Fill(dm, i);
if(s83visible) (*fTriggersHists)[5]->Fill(dm, i);
if(s83) (*fTriggersHists)[8]->Fill(dm, i);
if(s84visible) (*fTriggersHists)[11]->Fill(dm, i);
if(s84) (*fTriggersHists)[14]->Fill(dm, i);
}
}
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsBgFinder()
{
TDirectory* oldDir = this->MakeHistDir("bg", "BG_finder");
fBgFinderHists = GFHistManip::ClearCreateArray(fBgFinderHists);
fBgFinderHists->Add(new TH2F("bgBits", "BG finder bits (0-9 Ibg, 20-26 Ibgam)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 27, 0., 27.));
fBgFinderHists->Add(new TH2F("bgBitsS83", "BG finder bits (0-9 Ibg, 20-26 Ibgam) S83",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 27, 0., 27.));
fBgFinderHists->Add(new TH2F("bgBitsS84", "BG finder bits (0-9 Ibg, 20-26 Ibgam) S84",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 27, 0., 27.));
fAllHistArrays.Add(fBgFinderHists);
oldDir->cd();
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsBgFinder(const H1PartDstar* dstar)
{
if(!fBgFinderHists) this->CreateHistsBgFinder();
static H1IntPtr bg("Ibg");
for(Int_t i = 0; i < 10; ++i){
if(*bg & (1 << i)){
fBgFinderHists->At(0)->Fill(dstar->GetDm(), i);
if(this->IsS83(dstar)) static_cast<TH2*>(fBgFinderHists->At(1))
->Fill(dstar->GetDm(), i, this->GetWeight(83, dstar));
if(this->IsS84(dstar)) static_cast<TH2*>(fBgFinderHists->At(2))
->Fill(dstar->GetDm(), i, this->GetWeight(84, dstar));
}
}
static H1IntPtr bgam("Ibgam");
const Int_t o = 20; // offset for Ibgam
for(Int_t i = 0; i < 7; ++i){
if(*bgam & (1 << i)){
fBgFinderHists->At(0)->Fill(dstar->GetDm(), o + i);
if(this->IsS83(dstar)) static_cast<TH2*>(fBgFinderHists->At(1))
->Fill(dstar->GetDm(), o+i,this->GetWeight(83, dstar));
if(this->IsS84(dstar)) static_cast<TH2*>(fBgFinderHists->At(2))
->Fill(dstar->GetDm(), o+i,this->GetWeight(84, dstar));
}
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsEventClass(const H1PartDstar* dstar)
{
// assuming dstar fulfills cuts that you want
// only pt2 cut for S83 / eta for S84!
if(!fEventClassHists) this->CreateHistsEventClass();
const Double_t dm = dstar->GetDm();
const Bool_t s83 = this->IsS83NoL4Check(dstar);//this->IsS83(dstar);
const Bool_t s84 = this->IsS84NoL4Check(dstar);//this->IsS84(dstar);
// filling L4 open charm finder output:
Bool_t isL4Dstar = kFALSE;
static H1YeclBankPtr yecl; // output of original open charm finder on L4
if(yecl.GetEntries() <= 0){
// yecl.SetPointer("myYECL");
// if(yecl.GetEntries() <= 0){ // use for begin 99 e+ (if not excluded):
static Int_t numWarn = 0;
if(numWarn < kNumMaxWarn) {
if(numWarn == kNumMaxWarn-1)
this->Warning("FillHistsEventClass", "Following warning will not be shown further!");
this->Warning("FillHistsEventClass", "no YECL bank in run/event %d / %d, filling bit 0!",
gH1Tree->GetRunNumber(), gH1Tree->GetEventNumber());
}
++numWarn;
(*fOpenCharmFindHists)[0]->Fill(dm, 0);
if(s83) (*fOpenCharmFindHists)[1]->Fill(dm, 0);
if(s84) (*fOpenCharmFindHists)[2]->Fill(dm, 0);
isL4Dstar = kTRUE; // not to put in events without D* finding on L4... ??? FIXME
// }
// yecl.SetPointer("YECL");
} else {
for(Int_t bit = 1; bit < 32; ++bit){ // leave out '0' (empty anyway!)
if(this->IsL4Found(bit, NULL)){// NULL: no way out for high pt D*!
if(bit == kL4FinderBit) isL4Dstar = kTRUE;
(*fOpenCharmFindHists)[0]->Fill(dm, bit);
if(s83) (*fOpenCharmFindHists)[1]->Fill(dm, bit);
if(s84) (*fOpenCharmFindHists)[2]->Fill(dm, bit);
}
}
}
// filling event classes:
static H1FloatPtr weightL4("Weight1");
static H1BytePtr eventClass("Iclass");
for(Int_t i = 0; i < kNumEventClass; ++i){
if(eventClass[i] > 0){
static_cast<TH3*>((*fEventClassHists)[0])->Fill(dm, i, *weightL4);
if(eventClass[15] == 0) (*fEventClassHists)[1]->Fill(dm, i);
if(!isL4Dstar) (*fEventClassHists)[2]->Fill(dm, i);
if(s83) {
static_cast<TH3*>((*fEventClassHistsS83)[0])->Fill(dm, i, *weightL4);
if(eventClass[15] == 0) (*fEventClassHistsS83)[1]->Fill(dm, i); // if not L4 class 15
if(!isL4Dstar) (*fEventClassHistsS83)[2]->Fill(dm, i);
}
if(s84) {
static_cast<TH3*>((*fEventClassHistsS84)[0])->Fill(dm, i, *weightL4);
if(eventClass[15] == 0) (*fEventClassHistsS84)[1]->Fill(dm, i);
if(!isL4Dstar) (*fEventClassHistsS84)[2]->Fill(dm, i);
}
}
}
// filling weights:
(*fEventClassHists)[3]->Fill(dm, *weightL4);
if(s83) (*fEventClassHistsS83)[3]->Fill(dm, *weightL4);
if(s84) (*fEventClassHistsS84)[3]->Fill(dm, *weightL4);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsSumEt(const H1PartDstar* dstar, Double_t sumEt)
{
// if sumEt == 0. (default) take it from H1FloatPtr newSumEt("SumEt"), entry [1]
if(!fSumEtHistsS83) this->CreateHistsSumEt();
if(sumEt == 0) sumEt = this->GetSumEt();
const Double_t ptDstar = dstar->GetPt();
if(this->IsS83(dstar)){
const Double_t weight = this->GetWeight(83, dstar);
const Double_t dm = dstar->GetDm();
(*fSumEtHistsS83)[0]->Fill(sumEt, weight);
static_cast<TH2*>((*fSumEtHistsS83)[1])->Fill(dm, sumEt, weight);
if(sumEt > 0.){
(*fSumEtHistsS83)[2]->Fill(ptDstar/sumEt, weight);
static_cast<TH2*>((*fSumEtHistsS83)[3])->Fill(dm, ptDstar/sumEt, weight);
} else {
this->Warning("FillHistsSumEt", "transverse Energy = %f", sumEt);
}
}
if(this->IsS84(dstar)){
const Double_t weight = this->GetWeight(84, dstar);
const Double_t dm = dstar->GetDm();
(*fSumEtHistsS84)[0]->Fill(sumEt, weight);
static_cast<TH2*>((*fSumEtHistsS84)[1])->Fill(dm, sumEt, weight);
if(sumEt > 0.){
(*fSumEtHistsS84)[2]->Fill(ptDstar/sumEt, weight);
static_cast<TH2*>((*fSumEtHistsS84)[3])->Fill(dm, ptDstar/sumEt, weight);
} else {
this->Warning("FillHistsSumEt", "transverse Energy = %f", sumEt);
}
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTrackKind(const H1PartDstar* dstar)
{
if(!fTrackKindHists) this->CreateHistsTrackKind();
Int_t numCent = 0, numComb = 0, numForw = 0;
const H1PartSelTrack* kaon = dstar->GetKaon();
const H1PartSelTrack* pion = dstar->GetPion();
const H1PartSelTrack* slowPion = dstar->GetSlowPion();
if(kaon->IsCentralTrk()) ++numCent;
else if(kaon->IsCombinedTrk()) ++numComb;
else if(kaon->IsForwardTrk()) ++numForw;
if(pion->IsCentralTrk()) ++numCent;
else if(pion->IsCombinedTrk()) ++numComb;
else if(pion->IsForwardTrk()) ++numForw;
if(slowPion->IsCentralTrk()) ++numCent;
else if(slowPion->IsCombinedTrk()) ++numComb;
else if(slowPion->IsForwardTrk()) ++numForw;
if(numCent == 3) (*fTrackKindHists)[0]->Fill(dstar->GetDm());
else if(numCent == 2 && numComb == 1)(*fTrackKindHists)[1]->Fill(dstar->GetDm());
else if(numCent == 1 && numComb == 2)(*fTrackKindHists)[2]->Fill(dstar->GetDm());
else if(numCent == 2 && numForw == 1)(*fTrackKindHists)[3]->Fill(dstar->GetDm());
else if(numCent == 1 && numForw == 2)(*fTrackKindHists)[4]->Fill(dstar->GetDm());
else if(numComb == 3)(*fTrackKindHists)[5]->Fill(dstar->GetDm());
else if(numComb == 2 && numForw == 1)(*fTrackKindHists)[6]->Fill(dstar->GetDm());
else if(numComb == 1 && numForw == 2)(*fTrackKindHists)[7]->Fill(dstar->GetDm());
else if(numForw == 3)(*fTrackKindHists)[8]->Fill(dstar->GetDm());
else if(numCent == 1 && numComb == 1 && numForw == 1)
(*fTrackKindHists)[2]->Fill(dstar->GetDm());
else {
this->Warning("FillHistsTrackKind", "numCent = %d, numComb = %d, numForw = %d",
numCent, numComb, numForw);
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsCentralNotCombinedTrack(const H1PartDstar* dstar)
{
if(NULL == gH1Tree->GetChain(H1Tree::kOdsTree)){
this->Warning("FillHistsCentralNotCombinedTrack",
"No ODS tree, ignore!");
return;
}
const H1PartCand* kaon = dstar->GetKaon()->GetParticle();
const H1PartCand* pion = dstar->GetPion()->GetParticle();
const H1PartCand* slowPion = dstar->GetSlowPion()->GetParticle();
const Double_t mKaon = TDatabasePDG::Instance()->GetParticle(321)->Mass();
const Double_t mPion = TDatabasePDG::Instance()->GetParticle(211)->Mass();
TLorentzVector pKaon = dstar->GetKaon()->GetFourVector(mKaon);
TLorentzVector pPion = dstar->GetPion()->GetFourVector(mPion);
TLorentzVector pSlowPion = dstar->GetSlowPion()->GetFourVector(mPion);
Bool_t replaced = kFALSE;
if(kaon->IsCombinedTrk()){
TVector3 p3 = static_cast<const H1CombinedFittedTrack*>(kaon->GetTrack())
->GetCentralFittedTrack()->GetMomentum();
pKaon.SetVectM(p3, mKaon);
replaced = kTRUE;
}
if(pion->IsCombinedTrk()){
TVector3 p3 = static_cast<const H1CombinedFittedTrack*>(pion->GetTrack())
->GetCentralFittedTrack()->GetMomentum();
pPion.SetVectM(p3, mPion);
replaced = kTRUE;
}
if(slowPion->IsCombinedTrk()){
TVector3 p3 = static_cast<const H1CombinedFittedTrack*>(slowPion->GetTrack())
->GetCentralFittedTrack()->GetMomentum();
pPion.SetVectM(p3, mPion);
replaced = kTRUE;
}
if(replaced){
TLorentzVector pD0 = pKaon + pPion;
TLorentzVector pDstar = pD0 + pSlowPion;
(*fCentralNotCombinedTrackHists)[0]->Fill(dstar->GetDm());
(*fCentralNotCombinedTrackHists)[1]->Fill(pDstar.M() - pD0.M());
if(this->DstarCutLoose(dstar) && dstar->GetPt() > fDsCutPt2){
(*fCentralNotCombinedTrackHists)[2]->Fill(dstar->GetDm());
(*fCentralNotCombinedTrackHists)[3]->Fill(pDstar.M() - pD0.M());
}
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsL4Eff2()
{
fL4NoHistsS83 = this->CreateHists("NoL4Trig", "not L4 triggered", 83, 2); //2-D
fL4NoHistsS84 = this->CreateHists("NoL4Trig", "not L4 triggered", 84, 2); //2-D
fL4YesHistsS83 = this->CreateHists("L4Trig", "L4 triggered", 83, 2); //2-D
fL4YesHistsS84 = this->CreateHists("L4Trig", "L4 triggered", 84, 2); //2-D
fL4YesNotDsHistsS83 = this->CreateHists("L4TrigNotDs", "L4 triggered", 83, 2); //2-D
fL4YesNotDsHistsS84 = this->CreateHists("L4TrigNotDs", "L4 triggered", 84, 2); //2-D
fL4NoDstarHistsS83 = this->CreateHists("NoL4DsTrig", "not L4 (D*) triggered", 83, 2); //2-D
fL4NoDstarHistsS84 = this->CreateHists("NoL4DsTrig", "not L4 (D*) triggered", 84, 2); //2-D
fL4YesDstarHistsS83 = this->CreateHists("L4DsTrig", "L4 (D*) triggered", 83, 2); //2-D
fL4YesDstarHistsS84 = this->CreateHists("L4DsTrig", "L4 (D*) triggered", 84, 2); //2-D
TDirectory* oldDir = this->OpenHistsFile();
const Int_t nBinsTr = 7;
const Int_t firstBinTr = 5;
const Int_t lastBinTr = 40;
// no 2D behind, cf. GFDstarAnalyseData::CreateHistsL4Eff2()
fL4NoHistsS83->Add(new TH1F("nTrackEvHistNoL4TrigS83", "no L4! 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4NoHistsS84->Add(new TH1F("nTrackEvHistNoL4TrigS84", "no L4! 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4YesHistsS83->Add(new TH1F("nTrackEvHistL4TrigS83", "L4 && 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4YesHistsS84->Add(new TH1F("nTrackEvHistL4TrigS84", "L4 && 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4NoDstarHistsS83->Add(new TH1F("nTrackEvHistNoL4DsTrigS83",
"no L4 (D*)! 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4NoDstarHistsS84->Add(new TH1F("nTrackEvHistNoL4DsTrigS84",
"no L4 (D*)! 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4YesDstarHistsS83->Add(new TH1F("nTrackEvHistL4DsTrigS83",
"L4 (D*) && 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
fL4YesDstarHistsS84->Add(new TH1F("nTrackEvDHistL4DsTrigS84",
"L4 (D*) && 1 D* signal region",
nBinsTr, firstBinTr, lastBinTr));
oldDir->cd();
fAllHistArrays.Add(fL4NoHistsS83);
fAllHistArrays.Add(fL4NoHistsS84);
fAllHistArrays.Add(fL4YesHistsS83);
fAllHistArrays.Add(fL4YesHistsS84);
fAllHistArrays.Add(fL4YesNotDsHistsS83);
fAllHistArrays.Add(fL4YesNotDsHistsS84);
fAllHistArrays.Add(fL4NoDstarHistsS83);
fAllHistArrays.Add(fL4NoDstarHistsS84);
fAllHistArrays.Add(fL4YesDstarHistsS83);
fAllHistArrays.Add(fL4YesDstarHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsL4Eff2(const TObjArray* dstars)
{
if(!fL4NoHistsS83) this->CreateHistsL4Eff2(); // build hists
const Bool_t l4Ref83 = this->IsL4Ref(83);
const Bool_t l4Ref84 = this->IsL4Ref(84);
if(!l4Ref83 && !l4Ref84) return; // nothing to do...
static H1BytePtr l4class("Iclass");
Bool_t anyS83 = kFALSE;// event contains D* fullfilling cuts && L4!
Bool_t anyS84 = kFALSE;// dito
Bool_t anyS83Not = kFALSE; // ----------------"-------------- && not L4!
Bool_t anyS84Not = kFALSE; // dito
Bool_t anyS83Dstar = kFALSE;// now same for L4 D* finder
Bool_t anyS84Dstar = kFALSE;
Bool_t anyS83DstarNot = kFALSE;
Bool_t anyS84DstarNot = kFALSE;
TIter next(dstars);
while(H1PartDstar* ds = static_cast<H1PartDstar*>(next())){
if(this->DstarCutLoose(ds)){
if(this->IsS83NoL4Check(ds) && l4Ref83){ //
const Bool_t isDmSR = this->IsSignalRegionDm(ds);
if(l4class[15]){ // just check L4 class
this->FillHists(fL4YesHistsS83, ds, 83);
if(isDmSR) anyS83 = kTRUE;
if(!this->IsL4Found(kL4FinderBit, ds)) this->FillHists(fL4YesNotDsHistsS83, ds, 83);
} else {
this->FillHists(fL4NoHistsS83, ds, 83);
if(isDmSR) anyS83Not = kTRUE;
}
if(this->IsL4Found(kL4FinderBit, ds)){ // check L4 D* finder (in gamma p)
// if(!l4class[15])
// this->Error("FillHistsL4Eff2", "L4 D* bit, but not L4-class (83, run %d, evt %d)!",
// gH1Tree->GetRunNumber(), gH1Tree->GetEventNumber());
this->FillHists(fL4YesDstarHistsS83, ds, 83);
if(isDmSR) anyS83Dstar = kTRUE;
} else {
this->FillHists(fL4NoDstarHistsS83, ds, 83);
if(isDmSR) anyS83DstarNot = kTRUE;
}
}
if(this->IsS84NoL4Check(ds) && l4Ref84){
const Bool_t isDmSR = this->IsSignalRegionDm(ds);
if(l4class[15]){ // just check L4 class
this->FillHists(fL4YesHistsS84, ds, 84);
if(isDmSR) anyS84 = kTRUE;
if(!this->IsL4Found(kL4FinderBit, ds)) this->FillHists(fL4YesNotDsHistsS84, ds, 84);
} else {
this->FillHists(fL4NoHistsS84, ds, 84);
if(isDmSR) anyS84Not = kTRUE;
}
if(this->IsL4Found(kL4FinderBit, ds)){ // check L4 D* finder (in gamma p)
// if(!l4class[15])
// this->Error("FillHistsL4Eff2", "L4 D* bit, but not L4-class (84, run %d, evt %d)!",
// gH1Tree->GetRunNumber(), gH1Tree->GetEventNumber());
this->FillHists(fL4YesDstarHistsS84, ds, 84);
if(isDmSR) anyS84Dstar = kTRUE;
} else {
this->FillHists(fL4NoDstarHistsS84, ds, 84);
if(isDmSR) anyS84DstarNot = kTRUE;
}
}
} // D* cut
} // loop over D*
// event wise hist-filling for nTracks
// first L4class
// FIXME: proper treatment of rec-weigt in MC?
static H1ShortPtr numTracks("NumCentralTracks");
if(anyS83Not){
(*fL4NoHistsS83)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(83, NULL));
}
if(anyS84Not){
(*fL4NoHistsS84)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(84, NULL));
}
if(anyS83){
(*fL4YesHistsS83)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(83, NULL));
}
if(anyS84){
(*fL4YesHistsS84)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(84, NULL));
}
// now check for L4 D*:
if(anyS83DstarNot){
(*fL4NoDstarHistsS83)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(83, NULL));
}
if(anyS84DstarNot){
(*fL4NoDstarHistsS84)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(84, NULL));
}
if(anyS83Dstar){
(*fL4YesDstarHistsS83)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(83, NULL));
}
if(anyS84Dstar){
(*fL4YesDstarHistsS84)[kNumDefaultHists]->Fill(*numTracks, this->GetWeight(84, NULL));
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsJetProf()
{
TDirectory* oldDir = this->MakeHistDir("jetProf", "jet profile hitsograms");
GFHistManip::ClearCreateArray(fJet1ProfileHists);
GFHistManip::ClearCreateArray(fJet2ProfileHists);
const Int_t nBinsEta = 13;//16;
const Float_t minEta = -2.5;//-3.;
const Float_t maxEta = 4.5;//5.;
const Int_t nBinsPhi = 8;//11;
// overall profiles 1st jet
fJet1ProfileHists->Add(new TH2D("jet1EtaProfHist2D", "Jet profile (#eta): 1st jet",
nBinsEta, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D("jet1PhiProfHist2D", "Jet profile (#phi): 1st jet",
nBinsPhi, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
// pt bins for 1st jet
for(Int_t n = 1; n < fSteer->GetDstar1JetPtBins()->GetSize(); ++n){
fJet1ProfileHists->Add(new TH2D(Form("jet1EtaPt%dProfHist2D", n),
Form("#eta-profile: 1st jet, %.1f < p_{t}(jet) < %.1f",
fSteer->GetDstar1JetPtBins()->At(n-1),
fSteer->GetDstar1JetPtBins()->At(n)),
nBinsEta-((n+1)/4)*2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D(Form("jet1PhiPt%dProfHist2D", n),
Form("#phi-profile: 1st jet, %.1f < p_{t}(jet) < %.1f",
fSteer->GetDstar1JetPtBins()->At(n-1),
fSteer->GetDstar1JetPtBins()->At(n)),
// nBinsPhi-((n+1)/4)*2, -TMath::Pi(), TMath::Pi(),
nBinsPhi-((n+1)/4),0,TMath::Pi(),//-TMath::Pi(),TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
}
// overall profiles 2nd jet
fJet2ProfileHists->Add(new TH2D("jet2EtaProfHist2D", "Jet profile (#eta): 2nd jet",
nBinsEta, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D("jet2PhiProfHist2D", "Jet profile (#phi): 2nd jet",
nBinsPhi, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
// pt bins for 2nd jet
for(Int_t n = 1; n < fSteer->GetDstar1JetPtBins()->GetSize(); ++n){
fJet2ProfileHists->Add(new TH2D(Form("jet2EtaPt%dProfHist2D", n),
Form("#eta-profile: 2nd jet, %.1f < p_{t}(jet) < %.1f",
fSteer->GetDstar1JetPtBins()->At(n-1),
fSteer->GetDstar1JetPtBins()->At(n)),
nBinsEta-((n+1)/4)*2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D(Form("jet2PhiPt%dProfHist2D", n),
Form("#phi-profile: 2nd jet, %.1f < p_{t}(jet) < %.1f",
fSteer->GetDstar1JetPtBins()->At(n-1),
fSteer->GetDstar1JetPtBins()->At(n)),
// nBinsPhi-((n+1)/4)*2, -TMath::Pi(), TMath::Pi(),
nBinsPhi-((n+1)/4),0,TMath::Pi(),//-TMath::Pi(),TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
}
// eta bins 1st and 2nd jet
for(Int_t n = 1; n < fSteer->GetDstar1JetEtaBins()->GetSize(); ++n){
fJet1ProfileHists->Add(new TH2D(Form("jet1EtaEta%dProfHist2D", n),
Form("#eta-profile: 1st jet, %.1f < #eta(jet) < %.1f",
fSteer->GetDstar1JetEtaBins()->At(n-1),
fSteer->GetDstar1JetEtaBins()->At(n)),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D(Form("jet1PhiEta%dProfHist2D", n),
Form("#phi-profile: 1st jet, %.1f < #eta(jet) < %.1f",
fSteer->GetDstar1JetEtaBins()->At(n-1),
fSteer->GetDstar1JetEtaBins()->At(n)),
// nBinsPhi-2, -TMath::Pi(), TMath::Pi(),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D(Form("jet2EtaEta%dProfHist2D", n),
Form("#eta-profile: 2nd jet, %.1f < #eta(jet) < %.1f",
fSteer->GetDstar1JetEtaBins()->At(n-1),
fSteer->GetDstar1JetEtaBins()->At(n)),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D(Form("jet2PhiEta%dProfHist2D", n),
Form("#phi-profile: 2nd jet, %.1f < #eta(jet) < %.1f",
fSteer->GetDstar1JetEtaBins()->At(n-1),
fSteer->GetDstar1JetEtaBins()->At(n)),
// nBinsPhi-2, -TMath::Pi(), TMath::Pi(),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
}
fJet1ProfileHists->Add(new TH2D("jet1EtaxGam1ProfHist2D",
Form("#eta-profile: 1st jet, x_{#gamma}^{D*-jet} < %.2f",
fXgamProfile),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D("jet1PhixGam1ProfHist2D",
Form("#phi-profile: 1st jet, x_{#gamma}^{D*-jet} < %.2f",
fXgamProfile),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D("jet1EtaxGam2ProfHist2D",
Form("#eta-profile: 1st jet, x_{#gamma}^{D*-jet} > %.2f",
fXgamProfile),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet1ProfileHists->Add(new TH2D("jet1PhixGam2ProfHist2D",
Form("#phi-profile: 1st jet, x_{#gamma}^{D*-jet} > %.2f",
fXgamProfile),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D("jet2EtaxGam1ProfHist2D",
Form("#eta-profile: 2nd jet, x_{#gamma}^{D*-jet} < %.2f",
fXgamProfile),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D("jet2PhixGam1ProfHist2D",
Form("#phi-profile: 2nd jet, x_{#gamma}^{D*-jet} < %.2f",
fXgamProfile),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D("jet2EtaxGam2ProfHist2D",
Form("#eta-profile: 2nd jet, x_{#gamma}^{D*-jet} > %.2f",
fXgamProfile),
nBinsEta-2, minEta, maxEta,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
fJet2ProfileHists->Add(new TH2D("jet2PhixGam2ProfHist2D",
Form("#phi-profile: 2nd jet, x_{#gamma}^{D*-jet} > %.2f",
fXgamProfile),
nBinsPhi-1, 0., TMath::Pi(),//-TMath::Pi(), TMath::Pi(),
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
GFHistManip::CreateAnalogHists(fJet1ProfileHists, fJet1ProfileHistsS83, "S83");
GFHistManip::CreateAnalogHists(fJet2ProfileHists, fJet2ProfileHistsS83, "S83");
GFHistManip::CreateAnalogHists(fJet1ProfileHists, fJet1ProfileHistsS84, "S84");
GFHistManip::CreateAnalogHists(fJet2ProfileHists, fJet2ProfileHistsS84, "S84");
GFHistManip::CallSumw2(fJet1ProfileHists); GFHistManip::CallSumw2(fJet2ProfileHists);
GFHistManip::CallSumw2(fJet1ProfileHistsS83); GFHistManip::CallSumw2(fJet2ProfileHistsS83);
GFHistManip::CallSumw2(fJet1ProfileHistsS84); GFHistManip::CallSumw2(fJet2ProfileHistsS84);
for(Int_t i = 0; i < fJet1ProfileHists->GetEntriesFast(); i+=2){
// all
H1JetProfile2D* prof = new H1JetProfile2D; // first jet
prof->SetHistsEtaPhi((TH2*)(*fJet1ProfileHists)[i], (TH2*)(*fJet1ProfileHists)[i+1]);
fJetProfiler.Add(prof);
prof = new H1JetProfile2D; // 2nd jet
prof->SetHistsEtaPhi((TH2*)(*fJet2ProfileHists)[i], (TH2*)(*fJet2ProfileHists)[i+1]);
fJetProfiler.Add(prof);
// S83
prof = new H1JetProfile2D; // first jet
prof->SetHistsEtaPhi((TH2*)(*fJet1ProfileHistsS83)[i], (TH2*)(*fJet1ProfileHistsS83)[i+1]);
fJetProfilerS83.Add(prof);
prof = new H1JetProfile2D; // 2nd jet
prof->SetHistsEtaPhi((TH2*)(*fJet2ProfileHistsS83)[i], (TH2*)(*fJet2ProfileHistsS83)[i+1]);
fJetProfilerS83.Add(prof);
// S84
prof = new H1JetProfile2D; // first jet
prof->SetHistsEtaPhi((TH2*)(*fJet1ProfileHistsS84)[i], (TH2*)(*fJet1ProfileHistsS84)[i+1]);
fJetProfilerS84.Add(prof);
prof = new H1JetProfile2D; // 2nd jet
prof->SetHistsEtaPhi((TH2*)(*fJet2ProfileHistsS84)[i], (TH2*)(*fJet2ProfileHistsS84)[i+1]);
fJetProfilerS84.Add(prof);
}
oldDir->cd();
fAllHistArrays.Add(fJet1ProfileHists);
fAllHistArrays.Add(fJet1ProfileHistsS83);
fAllHistArrays.Add(fJet1ProfileHistsS84);
fAllHistArrays.Add(fJet2ProfileHists);
fAllHistArrays.Add(fJet2ProfileHistsS83);
fAllHistArrays.Add(fJet2ProfileHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsJetProf(const H1PartDstar* ds)
{
if(!fJet1ProfileHists) this->CreateHistsJetProf();
TObjArray parts;
H1InclHfsIterator hfsIter;
while(H1PartCand* pCand = hfsIter.Next()){
parts.Add(pCand);
}
// I do not believe in any electrons:
static H1FloatPtr q2("Q2e");
if(*q2 > 0.){ // true if any electron found...
static H1PartEmArrayPtr ems;
for(Int_t i = 0; i < ems.GetEntries(); ++i){
if(ems[i]->IsScatElec()){
parts.Add(const_cast<H1PartCand*>(ems[i]->GetParticle()));
break;
}
}
}
// TObjArray jets = this->GetJets(ds);
TObjArray jets;
this->GetNonDstarJets(ds, jets);
const Int_t nJets = jets.GetEntriesFast();
if(nJets == 0) return;
this->RemoveLowPt(&jets, 2);
H1PartJet *jet1 = static_cast<H1PartJet*>(jets[0]);
H1PartJet *jet2 = (nJets > 1 ? static_cast<H1PartJet*>(jets[1]) : NULL);
const Bool_t s83 = this->IsS83(ds);
const Bool_t s84 = this->IsS84(ds);
const Double_t dm = ds->GetDm();
const Double_t ptJet1 = jet1->GetPt();
const Double_t ptJet2 = (jet2 ? jet2->GetPt() : 0.);
const Double_t etaJet1 = jet1->GetEta();
const Double_t etaJet2 = (jet2 ? jet2->GetEta() : 0.);
const Double_t xGamJet1 = this->GetXgammaMassHelp(ds, jet1);
const Double_t xGamJet2 = (jet2 ? this->GetXgammaMassHelp(ds, jet2) : 0.);
Int_t pt1Bin = -1;
for(Int_t i = 1; i < fSteer->GetDstar1JetPtBins()->GetSize(); ++i){
if(fSteer->GetDstar1JetPtBins()->At(i-1) < ptJet1
&& ptJet1 <= fSteer->GetDstar1JetPtBins()->At(i)){
pt1Bin = i;
break;
}
}
Int_t pt2Bin = -1;
for(Int_t i = 1; jet2 && i < fSteer->GetDstar1JetPtBins()->GetSize(); ++i){
if(fSteer->GetDstar1JetPtBins()->At(i-1) < ptJet2
&& ptJet2 <= fSteer->GetDstar1JetPtBins()->At(i)){
pt2Bin = i;
break;
}
}
Int_t eta1Bin = -1;
for(Int_t i = 1; i < fSteer->GetDstar1JetEtaBins()->GetSize(); ++i){
if(fSteer->GetDstar1JetEtaBins()->At(i-1) < etaJet1
&& etaJet1 <= fSteer->GetDstar1JetEtaBins()->At(i)){
eta1Bin = i;
break;
}
}
Int_t eta2Bin = -1;
for(Int_t i = 1; jet2 && i < fSteer->GetDstar1JetEtaBins()->GetSize(); ++i){
if(fSteer->GetDstar1JetEtaBins()->At(i-1) < etaJet2
&& etaJet2 <= fSteer->GetDstar1JetEtaBins()->At(i)){
eta2Bin = i;
break;
}
}
const Int_t offsetEta = fSteer->GetDstar1JetPtBins()->GetSize() * 2; // pt bins and overall
const Int_t offsetXgam = offsetEta - 2 + fSteer->GetDstar1JetEtaBins()->GetSize() * 2;
static_cast<H1JetProfile2D*>(fJetProfiler[0])->FillProfile(jet1, &parts, dm);
if(jet2) static_cast<H1JetProfile2D*>(fJetProfiler[1])->FillProfile(jet2, &parts, dm);
if(pt1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfiler[pt1Bin*2])->FillProfile(jet1, &parts, dm);
}
if(pt2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfiler[pt2Bin*2+1])->FillProfile(jet2, &parts, dm);
}
if(eta1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfiler[offsetEta +eta1Bin*2 -2]) // bins start with 1
->FillProfile(jet1, &parts, dm);
}
if(eta2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfiler[offsetEta +eta2Bin*2 -1]) // dito
->FillProfile(jet2, &parts, dm);
}
if(xGamJet1 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfiler[offsetXgam])->FillProfile(jet1, &parts, dm);
} else {
static_cast<H1JetProfile2D*>(fJetProfiler[offsetXgam+1])->FillProfile(jet1, &parts, dm);
}
if(jet2 && xGamJet2 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfiler[offsetXgam+2])->FillProfile(jet2, &parts, dm);
} else if(jet2){
static_cast<H1JetProfile2D*>(fJetProfiler[offsetXgam+3])->FillProfile(jet2, &parts, dm);
}
if(s83){ // s83
const Double_t w = this->GetWeight(83, ds);
static_cast<H1JetProfile2D*>(fJetProfilerS83[0])->FillProfile(jet1, &parts, dm, w);
if(jet2) static_cast<H1JetProfile2D*>(fJetProfilerS83[1])->FillProfile(jet2, &parts, dm, w);
if(pt1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS83[pt1Bin*2])->FillProfile(jet1, &parts, dm, w);
}
if(pt2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS83[pt2Bin*2+1])->FillProfile(jet2, &parts,dm,w);
}
if(eta1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetEta +eta1Bin*2 -2])//bins start with 1
->FillProfile(jet1, &parts, dm, w);
}
if(eta2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetEta +eta2Bin*2 -1]) // dito
->FillProfile(jet2, &parts, dm, w);
}
if(xGamJet1 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetXgam])->FillProfile(jet1, &parts, dm);
} else {
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetXgam+1])->FillProfile(jet1,&parts, dm);
}
if(jet2 && xGamJet2 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetXgam+2])->FillProfile(jet2,&parts, dm);
} else if(jet2){
static_cast<H1JetProfile2D*>(fJetProfilerS83[offsetXgam+3])->FillProfile(jet2,&parts, dm);
}
}
if(s84){ // S84
const Double_t w = this->GetWeight(84, ds);
static_cast<H1JetProfile2D*>(fJetProfilerS84[0])->FillProfile(jet1, &parts, dm, w);
if(jet2) static_cast<H1JetProfile2D*>(fJetProfilerS84[1])->FillProfile(jet2, &parts, dm, w);
if(pt1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS84[pt1Bin*2])->FillProfile(jet1, &parts, dm, w);
}
if(pt2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS84[pt2Bin*2+1])->FillProfile(jet2, &parts,dm,w);
}
if(eta1Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetEta +eta1Bin*2 -2])//bins start with 1
->FillProfile(jet1, &parts, dm, w);
}
if(eta2Bin != -1) {
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetEta +eta2Bin*2 -1]) // dito
->FillProfile(jet2, &parts, dm, w);
}
if(xGamJet1 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetXgam])->FillProfile(jet1, &parts, dm);
} else {
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetXgam+1])->FillProfile(jet1,&parts, dm);
}
if(jet2 && xGamJet2 < fXgamProfile){
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetXgam+2])->FillProfile(jet2,&parts, dm);
} else if(jet2){
static_cast<H1JetProfile2D*>(fJetProfilerS84[offsetXgam+3])->FillProfile(jet2,&parts, dm);
}
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTrackClus()
{
TDirectory* oldDir = this->OpenHistsFile();
GFHistManip::ClearCreateArray(fTrackClusHists);
const Int_t nBins = 30;
const Float_t maxRatio = 2.5;
const Float_t maxE = 5;
fTrackClusHists->Add(new TH1F("ratioTrCls", "E_{cl} / E_{tr}", nBins, 0., maxRatio));
fTrackClusHists->Add(new TH2F("TrCls2D", "E_{cl} vs E_{tr}",
nBins, 0., maxE, nBins, 0., maxE));
fTrackClusHists->Add(new TH1F("ratioTrClsEm", "E_{cl} / E_{tr} (em clus)",
nBins, 0., maxRatio));
fTrackClusHists->Add(new TH2F("TrCls2DEm", "E_{cl} vs E_{tr} (em clus)",
nBins, 0., maxE, nBins, 0., maxE));
fTrackClusHists->Add(new TH1F("ratioTrClsHad", "E_{cl} / E_{tr} (had clus)",
nBins, 0., maxRatio));
fTrackClusHists->Add(new TH2F("TrCls2DHad", "E_{cl} vs E_{tr} (had clus)",
nBins, 0., maxE, nBins, 0., maxE));
GFHistManip::CreateAnalogHists(fTrackClusHists, fTrackClusHistsForw, "Forw");
GFHistManip::CreateAnalogHists(fTrackClusHists, fTrackClusHistsForwBar, "ForwBar");
GFHistManip::CreateAnalogHists(fTrackClusHists, fTrackClusHistsBar, "Bar");
GFHistManip::CreateAnalogHists(fTrackClusHists, fTrackClusHistsSpac, "Spac");
oldDir->cd();
fAllHistArrays.Add(fTrackClusHists);
fAllHistArrays.Add(fTrackClusHistsForw);
fAllHistArrays.Add(fTrackClusHistsForwBar);
fAllHistArrays.Add(fTrackClusHistsBar);
fAllHistArrays.Add(fTrackClusHistsSpac);
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTrackErr()
{
TDirectory *oldDir = this->MakeHistDir("trackErr", "track error histograms");
GFHistManip::ClearCreateArray(fTrackErrHists);
fTrackErrHists->Add(new TH1F("sigmaP", "#sigma(p);#sigma(p)", 30, 0., 0.3));
fTrackErrHists->Add(new TH1F("sigmaPrel", "#sigma(p)/p;#sigma(p)/p", 30, 0., .15));
fTrackErrHists->Add(new TH2F("sigmaPvsPt",
"#sigma_{p}(p_{t});p_{t}(track);#sigma(p)",
30, 0.12, 3.12, 30, 0, 0.3));
fTrackErrHists->Add(new TH2F("sigmaPrelvsPt",
"#sigma_{p}/p(p_{t});p_{t}(track);#sigma(p)/p",
30, 0.12, 3.12, 30, 0, .15));
fTrackErrHists->Add(new TH2F("sigmaPvsTheta",
"#sigma_{p}(#theta);#theta(track);#sigma(p)",
30, 0., TMath::Pi(), 30, 0, 0.3));
fTrackErrHists->Add(new TH2F("sigmaPrelvsTheta",
"#sigma_{p}/p(#theta);#theta(track);#sigma(p)/p",
30, 0., TMath::Pi(), 30, 0, .15));
fTrackErrHists->Add(new TProfile2D("sigmaPvsPtTheta",
"#LT#sigma_{p}#GT(p_{t},#theta);p_{t}(track};#theta(track);#sigma(p)/p",
30, 0.12, 3.12, 30, 0., TMath::Pi()));
fTrackErrHists->Add(new TProfile2D("sigmaPrelvsPtTheta",
"#LT#sigma_{p}/p#GT(p_{t},#theta);p_{t}(track};#theta(track);#sigma(p)/p",
30, 0.12, 3.12, 30, 0., TMath::Pi()));
oldDir->cd();
fAllHistArrays.Add(fTrackErrHists);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTrackErr(const TObjArray *dstars)
{
Bool_t dstarSR = kFALSE;
for (Int_t i = 0; i < dstars->GetEntriesFast(); ++i) {
const H1PartDstar *ds = static_cast<H1PartDstar*>(dstars->At(i));
if (this->IsSignalRegionDm(ds) && this->DstarCutLoose(ds) && this->IsS83(ds)) {
dstarSR = kTRUE;
break;
}
}
if (!dstarSR) return;
if (!fTrackErrHists) this->CreateHistsTrackErr();
static H1PartSelTrackArrayPtr tracks;
for (Int_t i = 0; i < tracks.GetEntries(); ++i) {
const H1PartSelTrack *tr = tracks[i];
if (!tr->IsCentralTrk()) continue;
const Double_t deltaP = tr->GetDp(); // it is a float...
const Double_t deltaPrel = deltaP/tr->GetP();
fTrackErrHists->At(0)->Fill(deltaP);
fTrackErrHists->At(1)->Fill(deltaPrel);
fTrackErrHists->At(2)->Fill(tr->GetPt(), deltaP);
fTrackErrHists->At(3)->Fill(tr->GetPt(), deltaPrel);
fTrackErrHists->At(4)->Fill(tr->GetTheta(), deltaP);
fTrackErrHists->At(5)->Fill(tr->GetTheta(), deltaPrel);
static_cast<TProfile2D*>(fTrackErrHists->At(6))
->Fill(tr->GetPt(), tr->GetTheta(), deltaP);
static_cast<TProfile2D*>(fTrackErrHists->At(7))
->Fill(tr->GetPt(), tr->GetTheta(), deltaPrel);
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTrackClus()
{
if(!fTrackClusHists) this->CreateHistsTrackClus();
const Double_t mPion = TDatabasePDG::Instance()->GetParticle(211)->Mass();
static H1PartSelTrackArrayPtr tracks;
Int_t nTrCls = 0;
Int_t nCents = 0;
Int_t nHfsNoMu = 0;
Int_t nCls = 0;
for(Int_t i = 0; i < tracks.GetEntries(); ++i){
if(tracks[i]->IsFromSecondary() || !tracks[i]->IsCentralTrk()) continue;
nCents++;
const H1PartCand* part = tracks[i]->GetParticle();
if(!part->IsHFS() || part->IsMuon() || part->IsEm()) continue;
++nHfsNoMu;
const H1Track* track = part->GetTrack();
if(!track->InheritsFrom(H1VertexFittedTrack::Class())){
tracks[i]->Print();
part->Print();
track->Print();
this->Error("FillHistsTrackClus", "Non vertex fitted track %d, skip", i);
continue;
}
const TObjArray* clusters = part->GetClusters();
if(!clusters) continue;
++nCls;
// part->Print();
// if(clusters->GetEntriesFast() != 1){
// part->Print();
// this->Warning("FillHistsTrackClus", "%d clusters, add", clusters->GetEntriesFast());
// }
// H1CaloGeometry::kCaloType calo = H1CaloGeometry::kUnknown;
Int_t calo = H1CaloGeometry::kUnknown;
Double_t eClus = 0;
for(Int_t i = 0; i < clusters->GetEntriesFast(); ++i){
H1Cluster* cl = static_cast<H1Cluster*>(clusters->At(i));
Int_t currentCalo = cl->GetCalo();
if(calo != H1CaloGeometry::kUnknown && currentCalo != calo){
this->Warning("FillHistsTrackClus", "Different calos! skip cluster");
continue;
}
calo = currentCalo;
eClus += cl->GetEnergy(part->IsHFSEm() ? H1Cell::kEm : H1Cell::kFinal);
}
Double_t eTrack = track->GetFourVector(mPion).E();
Double_t ratio = (eTrack ? eClus/eTrack : -1.);
(*fTrackClusHists)[0]->Fill(ratio);
(*fTrackClusHists)[1]->Fill(eClus, eTrack);
if(calo == H1CaloGeometry::kSpaCal){
(*fTrackClusHistsSpac)[0]->Fill(ratio);
(*fTrackClusHistsSpac)[1]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 55.) {
(*fTrackClusHistsBar)[0]->Fill(ratio);
(*fTrackClusHistsBar)[1]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 30.){
(*fTrackClusHistsForwBar)[0]->Fill(ratio);
(*fTrackClusHistsForwBar)[1]->Fill(eClus, eTrack);
} else {
(*fTrackClusHistsForw)[0]->Fill(ratio);
(*fTrackClusHistsForw)[1]->Fill(eClus, eTrack);
}
if(part->IsHFSEm()) {
(*fTrackClusHists)[2]->Fill(ratio);
(*fTrackClusHists)[3]->Fill(eClus, eTrack);
if(calo == H1CaloGeometry::kSpaCal){
(*fTrackClusHistsSpac)[2]->Fill(ratio);
(*fTrackClusHistsSpac)[3]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 55.) {
(*fTrackClusHistsBar)[2]->Fill(ratio);
(*fTrackClusHistsBar)[3]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 30.) {
(*fTrackClusHistsForwBar)[2]->Fill(ratio);
(*fTrackClusHistsForwBar)[3]->Fill(eClus, eTrack);
} else {
(*fTrackClusHistsForw)[2]->Fill(ratio);
(*fTrackClusHistsForw)[3]->Fill(eClus, eTrack);
}
} else {
(*fTrackClusHists)[4]->Fill(ratio);
(*fTrackClusHists)[5]->Fill(eClus, eTrack);
if(calo == H1CaloGeometry::kSpaCal){
(*fTrackClusHistsSpac)[4]->Fill(ratio);
(*fTrackClusHistsSpac)[5]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 55.) {
(*fTrackClusHistsBar)[4]->Fill(ratio);
(*fTrackClusHistsBar)[5]->Fill(eClus, eTrack);
} else if(track->GetTheta() * TMath::RadToDeg() > 30.) {
(*fTrackClusHistsForwBar)[4]->Fill(ratio);
(*fTrackClusHistsForwBar)[5]->Fill(eClus, eTrack);
} else {
(*fTrackClusHistsForw)[4]->Fill(ratio);
(*fTrackClusHistsForw)[5]->Fill(eClus, eTrack);
}
}
++nTrCls;
}
cout << tracks.GetEntries() << " "
<< nCents << " "
<< nHfsNoMu << " "
<< nCls << " "
<< nTrCls << endl;
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsJetComp()
{
TDirectory* oldDir = this->MakeHistDir("jetComp","jet E_t compositions");
GFHistManip::ClearCreateArray(fJetCompHists);
const Int_t nBins = 10;
const Float_t lower = 0.;
const Float_t upper = 1.000001;// to make 1.000 stay in upper bin
fJetCompHists->Add(new TH2F("jetCompEtTrack", "E_{t}(tr)/E_{t}",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F("jetCompEtEmCl", "E_{t}(em-cl)/E_{t}",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F("jetCompEtHadCl", "E_{t}(had-cl)/E_{t}",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH1F("jetCompEtTrackSR1D", "E_{t}(tr)/E_{t} (SR)",
nBins, lower, upper));
fJetCompHists->Add(new TH1F("jetCompEtEmClSR1D", "E_{t}(em-cl)/E_{t} (SR)",
nBins, lower, upper));
fJetCompHists->Add(new TH1F("jetCompEtHadClSR1D", "E_{t}(had-cl)/E_{t} (SR)",
nBins, lower, upper));
// fJetCompHists->Add(new TH1F("jetCompEtHadClSR1D",
for(Int_t i = 1; i < fSteer->GetDstar1JetEtaBins()->GetSize(); ++i){
TString range(Form("%.2f #leq #eta_{jet} < %.2f", fSteer->GetDstar1JetEtaBins()->At(i-1),
fSteer->GetDstar1JetEtaBins()->At(i)));
fJetCompHists->Add(new TH2F(Form("jetCompEtTrack%d", i), "E_{t}(tr)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F(Form("jetCompEtEmCl%d", i), "E_{t}(em-cl)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F(Form("jetCompEtHadCl%d", i), "E_{t}(had-cl)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtTrackSR1D%d", i), "E_{t}(tr)/E_{t}, " + range,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtEmClSR1D%d", i), "E_{t}(em-cl)/E_{t}, " + range,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtHadClSR1D%d", i),"E_{t}(had-cl)/E_{t}, " + range,
nBins, lower, upper));
}
for(Int_t i = 1; i < fSteer->GetDstar1JetPtBins()->GetSize(); ++i){
TString range(Form("%.2f #leq p_{t}^{jet} < %.2f", fSteer->GetDstar1JetPtBins()->At(i-1),
fSteer->GetDstar1JetPtBins()->At(i)));
fJetCompHists->Add(new TH2F(Form("jetCompEtTrackPt%d", i), "E_{t}(tr)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F(Form("jetCompEtEmClPt%d", i), "E_{t}(em-cl)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH2F(Form("jetCompEtHadClPt%d", i), "E_{t}(had-cl)/E_{t}, " + range,
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtTrackSR1DPt%d", i), "E_{t}(tr)/E_{t}, " + range,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtEmClSR1DPt%d",i), "E_{t}(em-cl)/E_{t}, " + range,
nBins, lower, upper));
fJetCompHists->Add(new TH1F(Form("jetCompEtHadClSR1DPt%d",i), "E_{t}(had-cl)/E_{t}, "+range,
nBins, lower, upper));
}
GFHistManip::CreateAnalogHists(fJetCompHists, fJetCompHistsS83, "S83", "");// no add to title!
GFHistManip::CreateAnalogHists(fJetCompHists, fJetCompHistsS84, "S84", "");// dito
GFHistManip::CallSumw2(fJetCompHistsS83);
GFHistManip::CallSumw2(fJetCompHistsS84);
oldDir->cd();
fAllHistArrays.Add(fJetCompHists);
fAllHistArrays.Add(fJetCompHistsS83);
fAllHistArrays.Add(fJetCompHistsS84);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsJetComp(const H1PartDstar* ds)
{
if(!fJetCompHists) this->CreateHistsJetComp();
// TObjArray jets = this->GetJets(ds);
TObjArray jets;
this->GetNonDstarJets(ds, jets);
if(this->RemoveLowPt(&jets, 2) == 0) return;
H1PartJet *jet1 = static_cast<H1PartJet*>(jets[0]);
// H1PartJet *jet2 = jets.GetEntries() > 1 ? static_cast<H1PartJet*>(jets[1]) : NULL;
const Double_t dm = ds->GetDm();
this->FillHistsJetComp(fJetCompHists, jet1, dm, 1.);
// if(jet2) this->FillHistsJetComp(fJetCompHists, jet2, dm);
if(this->IsS83(ds)) {
this->FillHistsJetComp(fJetCompHistsS83, jet1, dm, this->GetWeight(83, ds));
// if(jet2) this->FillHistsJetComp(fJetCompHistsS83, jet2, dm);
}
if(this->IsS84(ds)) {
this->FillHistsJetComp(fJetCompHistsS84, jet1, dm, this->GetWeight(84, ds));
// if(jet2) this->FillHistsJetComp(fJetCompHistsS84, jet2, dm);
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsJetComp(GFHistArray* hists, const H1PartJet* jet,
Double_t dm, Double_t weight)
{
Double_t ratioTrack = 0., ratioEmClus = 0., ratioHadClus = 0.;
this->GetEtRatios(jet, ratioTrack, ratioEmClus, ratioHadClus);
const Bool_t sigRegion = (dm < fSteer->GetD0DeltaMHigh() && dm >= fSteer->GetD0DeltaMLow());
static_cast<TH2*>((*hists)[0])->Fill(dm, ratioTrack, weight);
static_cast<TH2*>((*hists)[1])->Fill(dm, ratioEmClus, weight);
static_cast<TH2*>((*hists)[2])->Fill(dm, ratioHadClus, weight);
if(sigRegion){
(*hists)[3]->Fill(ratioTrack, weight);
(*hists)[4]->Fill(ratioEmClus, weight);
(*hists)[5]->Fill(ratioHadClus, weight);
}
const Double_t eta = jet->GetEta();
for(Int_t i = 1; i < fSteer->GetDstar1JetEtaBins()->GetSize(); ++i){
if(eta>=fSteer->GetDstar1JetEtaBins()->At(i-1) && eta<fSteer->GetDstar1JetEtaBins()->At(i)){
static_cast<TH2*>((*hists)[i*6 + 0])->Fill(dm, ratioTrack, weight);
static_cast<TH2*>((*hists)[i*6 + 1])->Fill(dm, ratioEmClus, weight);
static_cast<TH2*>((*hists)[i*6 + 2])->Fill(dm, ratioHadClus, weight);
if(sigRegion){
(*hists)[i*6 + 3]->Fill(ratioTrack, weight);
(*hists)[i*6 + 4]->Fill(ratioEmClus, weight);
(*hists)[i*6 + 5]->Fill(ratioHadClus, weight);
}
}
}
const Double_t pt = jet->GetPt();
for(Int_t i = 1; i < fSteer->GetDstar1JetPtBins()->GetSize(); ++i){
if(pt >= fSteer->GetDstar1JetPtBins()->At(i-1) && pt < fSteer->GetDstar1JetPtBins()->At(i)){
const Int_t offset = (i + fSteer->GetDstar1JetEtaBins()->GetSize() - 1) * 6;
static_cast<TH2*>((*hists)[offset + 0])->Fill(dm, ratioTrack, weight);
static_cast<TH2*>((*hists)[offset + 1])->Fill(dm, ratioEmClus, weight);
static_cast<TH2*>((*hists)[offset + 2])->Fill(dm, ratioHadClus, weight);
if(sigRegion){
(*hists)[offset + 3]->Fill(ratioTrack, weight);
(*hists)[offset + 4]->Fill(ratioEmClus, weight);
(*hists)[offset + 5]->Fill(ratioHadClus, weight);
}
}
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsDstarJets()
{
TDirectory* oldDir = this->MakeHistDir("dstarJets","D* and (other) jet hists");
fDstar1JetHists = this->CreateHistsDstar1Jet(0, 2); // no order, 2D
fDstarBack1JetHists = this->CreateHistsDstar1Jet(-1, 2); // D* is back, 2D
fDstarForw1JetHists = this->CreateHistsDstar1Jet(+1, 2); // D* is forw, 2D
// Do not forget to call Sumw2 in MC for new arrays!
GFHistManip::CreateAnalogHists(fDstar1JetHists, fDstar1JetHistsS83, "S83");
GFHistManip::CreateAnalogHists(fDstar1JetHists, fDstar1JetHistsS84, "S84");
GFHistManip::CreateAnalogHists(fDstarBack1JetHists, fDstarBack1JetHistsS83, "S83");
GFHistManip::CreateAnalogHists(fDstarBack1JetHists, fDstarBack1JetHistsS84, "S84");
GFHistManip::CreateAnalogHists(fDstarForw1JetHists, fDstarForw1JetHistsS83, "S83");
GFHistManip::CreateAnalogHists(fDstarForw1JetHists, fDstarForw1JetHistsS84, "S84");
oldDir->cd();
fAllHistArrays.Add(fDstar1JetHists);
fAllHistArrays.Add(fDstar1JetHistsS83);
fAllHistArrays.Add(fDstar1JetHistsS84);
fAllHistArrays.Add(fDstarBack1JetHists);
fAllHistArrays.Add(fDstarBack1JetHistsS83);
fAllHistArrays.Add(fDstarBack1JetHistsS84);
fAllHistArrays.Add(fDstarForw1JetHists);
fAllHistArrays.Add(fDstarForw1JetHistsS83);
fAllHistArrays.Add(fDstarForw1JetHistsS84);
}
//_____________________________________________________
GFHistArray* GFDstarAnalyse::CreateHistsDstar1Jet(Int_t forwBack, Int_t dim) const
{
// creating and fill array for D* + 1 jet hists, hist names
// (parallel to FillHistsDstar1Jet)
// - start with 'dstar1Jet',
// - then the variable follows ('', 'Pt/Eta'+'Ds/Jet', 'Dphi', 'Deta', 'Dpt')
// - then the identifier of eta-ordered set follows: 'Forw', 'Back', ''
//
// which of the last and the chosen binning depends on forwBack flag:
// -1: 'Back', i.e. eta(D*) < eta(jet)
// 0: '', i.e. no matter
// +1: 'Forw, i.e. eta(D*) > eta(jet)
//
// if dim = 2: x-axis is delta m(D*), else 1D hists
const TArrayD *ptBinsJet = NULL, *etaBinsJet = NULL, *ptBinsDs = NULL, *etaBinsDs = NULL;
const TArrayD *deltaEtaBins = fSteer->GetDstar1JetDetaBins();
const TArrayD *deltaPhiBins = fSteer->GetDstar1JetDphiBins();
const TArrayD* xGammaBins = fSteer->GetXgammaBins();
const TArrayD* binsDsJetPt = fSteer->GetDstar1JetDsJetPtBins();
const TArrayD* binsDsJetPt2 = fSteer->GetDstar1JetDsJetPt2Bins();
const TArrayD* binsDsJetMass = fSteer->GetDstar1JetDsJetMassBins();
const TArrayD* binsCosTh = fSteer->GetDstar1JetDsJetCosThBins();
const TArrayD* binsDsJetPtOrth = fSteer->GetDstar1JetDsJetPtOrthBins();
const TArrayD* logXpBins = fSteer->GetLogXpBins();
TArrayD deltaPhiNloBins(deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray());
deltaPhiNloBins[deltaPhiNloBins.GetSize()-1] = deltaPhiBins->At(deltaPhiBins->GetSize()-1);
TString id("");
switch(forwBack){
case -1:
id = "Back";
ptBinsDs = fSteer->GetDstarPtBins1JetB();
etaBinsDs = fSteer->GetDstarEtaBins1JetB();
ptBinsJet = fSteer->GetDstar1JetPtBinsB();
etaBinsJet = fSteer->GetDstar1JetEtaBinsB();
break;
case 0:
ptBinsDs = fSteer->GetDstarPtBins1Jet();
etaBinsDs = fSteer->GetDstarEtaBins1Jet();
ptBinsJet = fSteer->GetDstar1JetPtBins();
etaBinsJet = fSteer->GetDstar1JetEtaBins();
break;
case +1:
id = "Forw";
ptBinsDs = fSteer->GetDstarPtBins1JetF();
etaBinsDs = fSteer->GetDstarEtaBins1JetF();
ptBinsJet = fSteer->GetDstar1JetPtBinsF();
etaBinsJet = fSteer->GetDstar1JetEtaBinsF();
break;
default:
this->Error("CreateHistsDstar1Jet", "Don't know 'forwBack' = %d", forwBack);
return NULL;
}
GFHistArray* result = NULL;
GFHistManip::ClearCreateArray(result);
TString name("dstar1Jet");
const Double_t xGamBins[] = {0., 0.6, 1.}; // for double diff
if(dim == 2){ // dimensions of hists (2 means: one is delta M)
result->Add(new TH1F(name + id, "#Delta m, D* + 1 jet",
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
result->Add(new TH2F(name + "PtDs" + id, "p_{t}(D*), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBinsDs->GetSize()-1, ptBinsDs->GetArray()));
result->Add(new TH2F(name + "EtaDs" + id, "#eta(D*), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBinsDs->GetSize()-1, etaBinsDs->GetArray()));
result->Add(new TH2F(name + "PtJet" + id, "p_{t}(jet), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH2F(name + "EtaJet" + id, "#eta(jet), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH2F(name + "Dphi" + id, "#Delta #phi, D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray()));
result->Add(new TH2F(name + "Deta" + id, "#eta(D*) - #eta(jet), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaEtaBins->GetSize()-1, deltaEtaBins->GetArray()));
result->Add(new TH2F(name + "Dpt" + id, "#p_{t}(D*) - p_{t}(jet), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 7, -7., 3.5));
result->Add(new TH2F(name + "xGam" + id, "x_{#gamma}^{OBS}, D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
xGammaBins->GetSize()-1, xGammaBins->GetArray()));
result->Add(new TH2F(name + "Pt2DsJet" + id, "p_{t}^{2}(D*+j), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsDsJetPt2->GetSize()-1, binsDsJetPt2->GetArray()));
result->Add(new TH2F(name + "PtDsJet" + id, "p_{t}(D*+j), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsDsJetPt->GetSize()-1, binsDsJetPt->GetArray()));
result->Add(new TH2F(name + "MDsJet" + id, "M(D*+j), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsDsJetMass->GetSize()-1, binsDsJetMass->GetArray()));
result->Add(new TH2F(name + "CosTh" + id, "cos(thet*), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsCosTh->GetSize()-1, binsCosTh->GetArray()));
result->Add(new TH2F(name + "DeltaR" + id,"#Delta R = #sqrt{#eta^{2}+#phi^{2}}, D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 10, 0., 4.5));
result->Add(new TH2F(name + "PtDsJetOrth" + id, "p_{xy}(D*+j)#perp p_{xy}(j), D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsDsJetPtOrth->GetSize()-1, binsDsJetPtOrth->GetArray()));
result->Add(new TH2F(name + "xP" + id, "x_{p}(D*+j)^{OBS}, D* + 1 jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
logXpBins->GetSize()-1, logXpBins->GetArray()));
Double_t dmBins[kDmHistNumBins+1];
GFMath::EquidistBins(dmBins, kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge);
result->Add(new TH3F(name + "DphixGam" + id, "#Delta #phi vs x_{#gamma}, D* + 1 jet",
kDmHistNumBins, dmBins,
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray(),
sizeof(xGamBins)/sizeof(xGamBins[0])-1, xGamBins));
result->Add(new TH2F(name + "DphiNlo" + id, "#Delta #phi, D* + 1 jet (NLO-bins)",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaPhiNloBins.GetSize()-1, deltaPhiNloBins.GetArray()));
} else {
result->Add(new TH1F(name + id, "#Delta m, D* + 1 jet",
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
result->Add(new TH1F(name + "PtDs" + id, "p_{t}(D*), D* + 1 jet",
ptBinsDs->GetSize()-1, ptBinsDs->GetArray()));
result->Add(new TH1F(name + "EtaDs" + id, "#eta(D*), D* + 1 jet",
etaBinsDs->GetSize()-1, etaBinsDs->GetArray()));
result->Add(new TH1F(name + "PtJet" + id, "p_{t}(jet), D* + 1 jet",
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH1F(name + "EtaJet" + id, "#eta(jet), D* + 1 jet",
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH1F(name + "Dphi" + id, "#Delta #phi, D* + 1 jet",
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray()));
result->Add(new TH1F(name + "Deta" + id, "#eta(D*) - #eta(jet), D* + 1 jet",
deltaEtaBins->GetSize()-1, deltaEtaBins->GetArray()));
result->Add(new TH1F(name + "Dpt" + id, "#p_{t}(D*) - p_{t}(jet), D* + 1 jet", 7, -7.,3.5));
result->Add(new TH1F(name + "xGam" + id, "x_{#gamma}^{OBS}, D* + 1 jet",
xGammaBins->GetSize()-1, xGammaBins->GetArray()));
result->Add(new TH1F(name + "PtDsJet" + id, "p_{t}(D*+j), D* + 1 jet",
binsDsJetPt->GetSize()-1, binsDsJetPt->GetArray()));
result->Add(new TH1F(name + "Pt2DsJet" + id, "p_{t}^{2}(D*+j), D* + 1 jet",
binsDsJetPt2->GetSize()-1, binsDsJetPt2->GetArray()));
result->Add(new TH1F(name + "MDsJet" + id, "M(D*+j), D* + 1 jet",
binsDsJetMass->GetSize()-1, binsDsJetMass->GetArray()));
result->Add(new TH1F(name + "CosTh" + id, "cos(thet*), D* + 1 jet",
binsCosTh->GetSize()-1, binsCosTh->GetArray()));
result->Add(new TH1F(name + "DeltaR" + id,"#Delta R = #sqrt{#eta^{2}+#phi^{2}}, D* + 1 jet",
10, 0., 4.5));
result->Add(new TH1F(name + "PtDsJetOrth" + id, "p_{xy}(D*+j)#perp p_{xy}(j), D* + 1 jet",
binsDsJetPtOrth->GetSize()-1, binsDsJetPtOrth->GetArray()));
result->Add(new TH1F(name + "xP" + id, "x_{p}(D*+j)^{OBS}, D* + 1 jet",
logXpBins->GetSize()-1, logXpBins->GetArray()));
result->Add(new TH2F(name + "DphixGam" + id, "#Delta #phi vs x_{#gamma}, D* + 1 jet",
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray(),
sizeof(xGamBins)/sizeof(xGamBins[0])-1, xGamBins));
result->Add(new TH1F(name + "DphiNlo" + id, "#Delta #phi, D* + 1 jet (NLO-bins)",
deltaPhiNloBins.GetSize()-1, deltaPhiNloBins.GetArray()));
}
if(forwBack == +1){
GFHistManip::AddToHistsName(NULL, result, " #eta(D*) > #eta(jet)");
} else if(forwBack == -1){
GFHistManip::AddToHistsName(NULL, result, " #eta(D*) < #eta(jet)");
}
return result;
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsDstarJets(const H1PartDstar* ds)
{
// hists for D* and (other) jet(s)
if(!fDstar1JetHists) this->CreateHistsDstarJets();
TObjArray jets;
const Int_t posDstar = this->GetNonDstarJets(ds, jets);
if(jets.GetEntriesFast() == 0) return; //no jets
const Bool_t s83 = this->IsS83(ds);
const Bool_t s84 = this->IsS84(ds);
const Double_t dm = ds->GetDm();
const Double_t weight = this->GetWeight(-1, ds);
const Double_t weight83 = this->GetWeight(83, ds);
const Double_t weight84 = this->GetWeight(84, ds);
const Int_t iJetHighPt = this->GetHighestPt(jets);
const H1PartJet* jetHighPt = static_cast<H1PartJet*>(jets[iJetHighPt]);
this->FillHistsDstar1Jet(ds, jetHighPt, fDstar1JetHists, weight, dm);
if(s83) this->FillHistsDstar1Jet(ds, jetHighPt, fDstar1JetHistsS83, weight83, dm);
if(s84) this->FillHistsDstar1Jet(ds, jetHighPt, fDstar1JetHistsS84, weight84, dm);
if(iJetHighPt >= posDstar){
this->FillHistsDstar1Jet(ds, jetHighPt, fDstarBack1JetHists, weight, dm);
if(s83) this->FillHistsDstar1Jet(ds, jetHighPt, fDstarBack1JetHistsS83, weight83, dm);
if(s84) this->FillHistsDstar1Jet(ds, jetHighPt, fDstarBack1JetHistsS84, weight84, dm);
} else {
this->FillHistsDstar1Jet(ds, jetHighPt, fDstarForw1JetHists, weight, dm);
if(s83) this->FillHistsDstar1Jet(ds, jetHighPt, fDstarForw1JetHistsS83, weight83, dm);
if(s84) this->FillHistsDstar1Jet(ds, jetHighPt, fDstarForw1JetHistsS84, weight84, dm);
}
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsDstar1Jet(const H1Part* ds, const H1PartJet* jet,
GFHistArray* hists, Double_t weight,
Double_t dm)
{
// fill hists with D* + 1 jet quantities with the weight
// if(dm>0) assuming 2D hist where x is filled with 'dm'
// (parallel to CreateHistsDstar1Jet)
const Double_t xGamma = this->GetXgammaMassHelp(ds, jet);
const Double_t logXp = TMath::Log10(this->GetXpMassHelp(ds, jet));
const TLorentzVector ds4vec = ds->GetFourVector();
const TLorentzVector jet4vec = jet->GetFourVector();
const Double_t cosTheta = TMath::TanH((ds4vec.Eta() - jet4vec.Eta())/2.);
const TVector2 jetXYvect = jet4vec.Vect().XYvector();
const TVector2 dsJetXYvect = jetXYvect + ds4vec.Vect().XYvector();
const Double_t dPhi = TMath::Abs(jet4vec.DeltaPhi(ds4vec))*TMath::RadToDeg();
if(dm > 0.){
hists->At(0)->Fill(dm, weight);
static_cast<TH2*>(hists->At(1))->Fill(dm, ds4vec.Pt(), weight);
static_cast<TH2*>(hists->At(2))->Fill(dm, ds4vec.Eta(), weight);
static_cast<TH2*>(hists->At(3))->Fill(dm, jet4vec.Pt(), weight);
static_cast<TH2*>(hists->At(4))->Fill(dm, jet4vec.Eta(), weight);
static_cast<TH2*>(hists->At(5))->Fill(dm, dPhi, weight);
static_cast<TH2*>(hists->At(6))->Fill(dm, ds4vec.Eta() - jet4vec.Eta(), weight);
static_cast<TH2*>(hists->At(7))->Fill(dm, ds4vec.Pt() - jet4vec.Pt(), weight);
static_cast<TH2*>(hists->At(8))->Fill(dm, xGamma, weight);
static_cast<TH2*>(hists->At(9))->Fill(dm, (jet4vec + ds4vec).Pt(), weight);
static_cast<TH2*>(hists->At(10))->Fill(dm, (jet4vec + ds4vec).Perp2(), weight);
static_cast<TH2*>(hists->At(11))->Fill(dm, (jet4vec + ds4vec).M(), weight);
static_cast<TH2*>(hists->At(12))->Fill(dm, cosTheta, weight);
static_cast<TH2*>(hists->At(13))->Fill(dm, jet4vec.DeltaR(ds4vec), weight);
static_cast<TH2*>(hists->At(14))->Fill(dm, dsJetXYvect.Norm(jetXYvect).Mod2(), weight);
static_cast<TH2*>(hists->At(15))->Fill(dm, logXp, weight);
static_cast<TH3*>(hists->At(16))->Fill(dm, dPhi, xGamma, weight);
static_cast<TH2*>(hists->At(17))->Fill(dm, dPhi, weight);
} else {
hists->At(0)->Fill(0.145, weight);
hists->At(1)->Fill(ds4vec.Pt(), weight);
hists->At(2)->Fill(ds4vec.Eta(), weight);
hists->At(3)->Fill(jet4vec.Pt(), weight);
hists->At(4)->Fill(jet4vec.Eta(), weight);
hists->At(5)->Fill(dPhi, weight);
hists->At(6)->Fill(ds4vec.Eta() - jet4vec.Eta(), weight);
hists->At(7)->Fill(ds4vec.Pt() - jet4vec.Pt(), weight);
hists->At(8)->Fill(xGamma, weight);
hists->At(9)->Fill((jet4vec + ds4vec).Pt(), weight);
hists->At(10)->Fill((jet4vec + ds4vec).Perp2(), weight);
hists->At(11)->Fill((jet4vec + ds4vec).M(), weight);
hists->At(12)->Fill(cosTheta, weight);
hists->At(13)->Fill(jet4vec.DeltaR(ds4vec), weight);
hists->At(14)->Fill(dsJetXYvect.Norm(jetXYvect).Mod2(), weight);
hists->At(15)->Fill(logXp, weight);
static_cast<TH2*>(hists->At(16))->Fill(dPhi, xGamma, weight);
hists->At(17)->Fill(dPhi, weight);
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsDstarDiJet()
{
// hists for FillHistsDstarDiJets(const H1PartDstar* ds)
TDirectory* oldDir = this->MakeHistDir("dsDiJets","D*-dijet hists");
// Do not forget to call Sumw2 in MC for new arrays!
fDstarDiJetHists = this->CreateHistsDstarDiJetInt(2);// 2D
GFHistManip::CreateAnalogHists(fDstarDiJetHists, fDstarDiJetHistsS83, "S83");
GFHistManip::CreateAnalogHists(fDstarDiJetHists, fDstarDiJetHistsS84, "S84");
fAllHistArrays.Add(fDstarDiJetHists);
fAllHistArrays.Add(fDstarDiJetHistsS83);
fAllHistArrays.Add(fDstarDiJetHistsS84);
oldDir->cd();
}
//_____________________________________________________
GFHistArray* GFDstarAnalyse::CreateHistsDstarDiJetInt(Int_t dim) const
{
// creating and fill array for D*-dijet hists, hist names
// (parallel to FillHistsDstarDiJet)
// - start with 'dstarDiJet',
// - then the variable follows ('', 'Pt/Eta'+'Ds/DJet/OthJet', 'Dphi', 'Deta', 'Dpt')
//
// if dim = 2: x-axis is delta m(D*), else 1D hists
const TArrayD *ptBinsDs = fSteer->GetDstarPtBins1Jet();
const TArrayD *ptBinsJet = fSteer->GetDstar1JetPtBins();
const TArrayD *etaBinsDs = fSteer->GetDstarEtaBins1Jet();
const TArrayD *etaBinsJet = fSteer->GetDstar1JetEtaBins();
const TArrayD *binsPtDiJet = fSteer->GetDstar1JetDsJetPtBins();
const TArrayD *binsPt2DiJet = fSteer->GetDstar1JetDsJetPt2Bins();
const TArrayD *deltaPhiBins = fSteer->GetDstar1JetDphiBins();
const TArrayD *deltaEtaBins = fSteer->GetDstar1JetDetaBins();
const TArrayD *binsCosTh = fSteer->GetDstar1JetDsJetCosThBins();
const TArrayD *xGammaBins = fSteer->GetXgammaDijetBins();
TArrayD deltaPhiNloBins(deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray());
deltaPhiNloBins[deltaPhiNloBins.GetSize()-1] = deltaPhiBins->At(deltaPhiBins->GetSize()-1);
GFHistArray* result = NULL;
GFHistManip::ClearCreateArray(result);
TString name("dstarDiJet");
const Double_t deltaRbins[] = {0., 0.2, 0.4, 0.6, 0.9, 4.};
if(dim == 2){ // dimensions of hists (2 means: one is delta M)
result->Add(new TH1F(name, "#Delta m, D*-dijet",
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
result->Add(new TH2F(name + "PtDs", "p_{t}(D*), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBinsDs->GetSize()-1, ptBinsDs->GetArray()));
result->Add(new TH2F(name + "PtDJet", "p_{t}(D*-jet), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH2F(name + "PtOthJet", "p_{t}(other jet), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH2F(name + "EtaDs", "#eta(D*), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBinsDs->GetSize()-1, etaBinsDs->GetArray()));
result->Add(new TH2F(name + "EtaDJet", "#eta(D*-jet), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH2F(name + "EtaOthJet", "#eta(other jet), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH2F(name + "DphiDijet", "#Delta#phi(D*j,oj), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray()));
result->Add(new TH2F(name + "PtDijet", "p_{t}(D*j+oj), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsPtDiJet->GetSize()-1, binsPtDiJet->GetArray()));
result->Add(new TH2F(name + "Pt2Dijet", "p_{t}^{2}(D*j+oj), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsPt2DiJet->GetSize()-1, binsPt2DiJet->GetArray()));
result->Add(new TH2F(name + "DetaDijet", "#eta(D*j)-#eta(oj), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaEtaBins->GetSize()-1, deltaEtaBins->GetArray()));
result->Add(new TH2F(name + "CosThDijet", "cos(#theta*), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
binsCosTh->GetSize()-1, binsCosTh->GetArray()));
result->Add(new TH2F(name + "xGamDijet", "x_{#gamma}, D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
xGammaBins->GetSize()-1, xGammaBins->GetArray()));
result->Add(new TH2F(name + "DeltaRDijet","#Delta R(D*,D*jet), D*-dijet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
sizeof(deltaRbins)/sizeof(deltaRbins[0])-1, deltaRbins));//10,0.,4.5
result->Add(new TH2F(name + "DeltaRNoDjet","#Delta R(D*,D*jet), D*+2 other jet",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 10, 0., 4.5));
result->Add(new TH2F(name + "DphiDijetNlo", "#Delta#phi(D*j,oj), D*-dijet, NLO-bins",
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
deltaPhiNloBins.GetSize()-1, deltaPhiNloBins.GetArray()));
// log(xp)
} else {
result->Add(new TH1F(name, "#Delta m, D*-dijet",
2*kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
result->Add(new TH1F(name + "PtDs", "p_{t}(D*), D*-dijet",
ptBinsDs->GetSize()-1, ptBinsDs->GetArray()));
result->Add(new TH1F(name + "PtDJet", "p_{t}(D*-jet), D*-dijet",
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH1F(name + "PtOthJet", "p_{t}(other jet), D*-dijet",
ptBinsJet->GetSize()-1, ptBinsJet->GetArray()));
result->Add(new TH1F(name + "EtaDs", "#eta(D*), D*-dijet",
etaBinsDs->GetSize()-1, etaBinsDs->GetArray()));
result->Add(new TH1F(name + "EtaDJet", "#eta(D*-jet), D*-dijet",
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH1F(name + "EtaOthJet", "#eta(other jet), D*-dijet",
etaBinsJet->GetSize()-1, etaBinsJet->GetArray()));
result->Add(new TH1F(name + "DphiDijet", "#Delta#phi(D*j,oj), D*-dijet",
deltaPhiBins->GetSize()-1, deltaPhiBins->GetArray()));
result->Add(new TH1F(name + "PtDijet", "p_{t}(D*j+oj), D*-dijet",
binsPtDiJet->GetSize()-1, binsPtDiJet->GetArray()));
result->Add(new TH1F(name + "Pt2Dijet", "p_{t}^{2}(D*j+oj), D*-dijet",
binsPt2DiJet->GetSize()-1, binsPt2DiJet->GetArray()));
result->Add(new TH1F(name + "DetaDijet", "#eta(D*j)-#eta(oj), D*-dijet",
deltaEtaBins->GetSize()-1, deltaEtaBins->GetArray()));
result->Add(new TH1F(name + "CosThDijet", "cos(#theta*), D*-dijet",
binsCosTh->GetSize()-1, binsCosTh->GetArray()));
result->Add(new TH1F(name + "xGamDijet", "x_{#gamma}, D*-dijet",
xGammaBins->GetSize()-1, xGammaBins->GetArray()));
result->Add(new TH1F(name + "DeltaRDijet","#Delta R(D*,D*jet), D*-dijet",
sizeof(deltaRbins)/sizeof(deltaRbins[0])-1, deltaRbins));//10,0.,4.5
result->Add(new TH1F(name + "DeltaRNoDjet","#Delta R(D*,D*jet), D*+2 other jet",
10, 0., 4.5));
result->Add(new TH1F(name + "DphiDijetNlo", "#Delta#phi(D*j,oj), D*-dijet, NLO-bins",
deltaPhiNloBins.GetSize()-1, deltaPhiNloBins.GetArray()));
// log(xp)
}
return result;
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsDstarDiJet(const H1PartDstar* ds)
{
// dijet histos
if(!fDstarDiJetHists) this->CreateHistsDstarDiJet();
const H1PartJet *dstarJet = NULL;
const H1PartJet *oJet = NULL;
const Bool_t realDsJet = this->GetDiJets(ds, dstarJet, oJet);
if(!dstarJet || !oJet) {
return; // need both jets...
}
const Double_t dm = ds->GetDm();
this->FillHistsDstarDiJet(ds, dstarJet, oJet, realDsJet, fDstarDiJetHists, this->GetWeight(-1, ds), dm);
if(this->IsS83(ds)) this->FillHistsDstarDiJet(ds, dstarJet, oJet, realDsJet, fDstarDiJetHistsS83, this->GetWeight(83, ds),dm);
if(this->IsS84(ds)) this->FillHistsDstarDiJet(ds, dstarJet, oJet, realDsJet, fDstarDiJetHistsS84, this->GetWeight(84, ds),dm);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsDstarDiJet(const H1Part* ds, const H1PartJet* dsJet,
const H1PartJet* oJet, Bool_t realDsJet, GFHistArray* hists, Double_t weight,
Double_t dm)
{
// fill hists with D* -dijet quantities with the weight
// if(dm>0) assuming 2D hist where x is filled with 'dm'
// (parallel to CreateHistsDstarDiJet)
const TLorentzVector ds4vec = ds->GetFourVector();
const TLorentzVector dsJet4vec = dsJet->GetFourVector();
const TLorentzVector oJet4vec = oJet->GetFourVector();
const Double_t cosTheta = TMath::TanH((dsJet4vec.Eta() - oJet4vec.Eta())/2.);
const Double_t xGamma = this->GetXgammaMassDiJetHelp(ds, dsJet, oJet);
const Double_t dPhi = TMath::Abs(dsJet4vec.DeltaPhi(oJet4vec))*TMath::RadToDeg();
if(dm > 0.){
hists->At(0)->Fill(dm, weight);
static_cast<TH2*>(hists->At(1))->Fill(dm, ds4vec.Pt(), weight);
static_cast<TH2*>(hists->At(2))->Fill(dm, dsJet4vec.Pt(), weight);
static_cast<TH2*>(hists->At(3))->Fill(dm, oJet4vec.Pt(), weight);
static_cast<TH2*>(hists->At(4))->Fill(dm, ds4vec.Eta(), weight);
static_cast<TH2*>(hists->At(5))->Fill(dm, dsJet4vec.Eta(), weight);
static_cast<TH2*>(hists->At(6))->Fill(dm, oJet4vec.Eta(), weight);
static_cast<TH2*>(hists->At(7))->Fill(dm, dPhi, weight);
static_cast<TH2*>(hists->At(8))->Fill(dm, (dsJet4vec+oJet4vec).Pt(), weight);
static_cast<TH2*>(hists->At(9))->Fill(dm, (dsJet4vec+oJet4vec).Perp2(), weight);
static_cast<TH2*>(hists->At(10))->Fill(dm, dsJet4vec.Eta() - oJet4vec.Eta(), weight);
static_cast<TH2*>(hists->At(11))->Fill(dm, cosTheta, weight);
static_cast<TH2*>(hists->At(12))->Fill(dm, xGamma, weight);
static_cast<TH2*>(hists->At(13))->Fill(dm, dsJet4vec.DeltaR(ds4vec), weight);
if(!realDsJet) static_cast<TH2*>(hists->At(14))->Fill(dm, dsJet4vec.DeltaR(ds4vec), weight);
static_cast<TH2*>(hists->At(15))->Fill(dm, dPhi, weight);
// static_cast<TH2*>(hists->At(15))->Fill(dm, logXp, weight);
} else {
hists->At(0)->Fill(0.145, weight);
hists->At(1)->Fill(ds4vec.Pt(), weight);
hists->At(2)->Fill(dsJet4vec.Pt(), weight);
hists->At(3)->Fill(oJet4vec.Pt(), weight);
hists->At(4)->Fill(ds4vec.Eta(), weight);
hists->At(5)->Fill(dsJet4vec.Eta(), weight);
hists->At(6)->Fill(oJet4vec.Eta(), weight);
hists->At(7)->Fill(dPhi, weight);
hists->At(8)->Fill((dsJet4vec+oJet4vec).Pt(), weight);
hists->At(9)->Fill((dsJet4vec+oJet4vec).Perp2(), weight);
hists->At(10)->Fill(dsJet4vec.Eta() - oJet4vec.Eta(), weight);
hists->At(11)->Fill(cosTheta, weight);
hists->At(12)->Fill(xGamma, weight);
hists->At(13)->Fill(dsJet4vec.DeltaR(ds4vec), weight);
if(!realDsJet) hists->At(14)->Fill(dsJet4vec.DeltaR(ds4vec), weight);
hists->At(15)->Fill(dPhi, weight);
// hists->At(15)->Fill(logXp, weight);
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsD0Checks()
{
TDirectory* oldDir = this->MakeHistDir("D0Checks", "check plots for D0 witdh");
fAllHistArrays.Add(GFHistManip::ClearCreateArray(fD0CheckHistsS83));
fD0CheckHistsS83->Add(new TH1F("D0noL4cutS83", "m(K#pi), no L4 cut",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge));
TH1* hWeight = new TH1F("D0L4weightS83", "m(K#pi), no L4 cut, L4 weight",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge);
hWeight->Sumw2();
fD0CheckHistsS83->Add(hWeight);
fD0CheckHistsS83->Add(new TH1F("D0not15L4S83",
"m(K#pi), L4 class 4,5,6,7,8,9,10,11,12,13,14,16,17 or 18",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge));
//
fAllHistArrays.Add(GFHistManip::ClearCreateArray(fD0CheckHists));
fD0CheckHists->Add(new TH2F("ptHist2DnoL4biasD0",
"m(K#pi), L4 class 4,5,6,7,8,9,10,11,12,13,14,16,17 or 18)",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetPtBins33()->GetSize()-1,
fSteer->GetPtBins33()->GetArray()));
fD0CheckHists->Add(new TH2F("etaHist2DnoL4biasD0",
"m(K#pi), L4 class 4,5,6,7,8,9,10,11,12,13,14,16,17 or 18)",
kD0HistNumBins, kD0HistLowEdge, kD0HistHighEdge,
fSteer->GetEtaBins33()->GetSize()-1,
fSteer->GetEtaBins33()->GetArray()));
oldDir->cd();
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsD0Checks(const H1PartDstar* ds)
{
if(!fD0CheckHistsS83) this->CreateHistsD0Checks();
if(!this->D0CutLoose(ds)) return;
const Int_t help[] = {4,5,6,7,8,9,10,11,12,13,14,16,17,18};// change also hist title
const TArrayI otherClasses(sizeof(help)/sizeof(help[0]), help);
Bool_t otherL4 = this->OneL4ClassGiven(otherClasses);
if(!otherL4){
const H1PartDstar *dsTmp = NULL; // NULL No way out for high pt D*
if(this->IsL4Found(kL4FinderBit, dsTmp, 0)) otherL4 = kTRUE;// 0: DIS mode = wide D0 window!
}
const Double_t mkpi = ds->GetD0Mass();
if(this->IsS83NoL4Check(ds)){
static H1FloatPtr wL4("Weight1");
fD0CheckHistsS83->At(0)->Fill(mkpi);
fD0CheckHistsS83->At(1)->Fill(mkpi, *wL4);
if(otherL4) fD0CheckHistsS83->At(2)->Fill(mkpi);
}
if(otherL4){
fD0CheckHists->At(0)->Fill(mkpi, ds->GetPt());
fD0CheckHists->At(1)->Fill(mkpi, ds->GetEta());
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsJetMult()
{
TDirectory* oldDir = this->MakeHistDir("jetMul", "jet multiplicity");
GFHistManip::ClearCreateArray(fJetMultHistsS83);
fJetMultHistsS83->Add(new TH2F("nJet", "N(jet) ", //0
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 7, 0, 7));
fJetMultHistsS83->Add(new TH2F("nJetNonDstar", "N(jet w/o D*) ", //1
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 7, 0, 7));
fJetMultHistsS83->Add(new TH2F("dsPlusOtherJet","0: not 2 jets, 1: 2 jets",//2
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 2, 0, 2));
fJetMultHistsS83->Add(new TH1F("dsNotIn2Jets","D*+jet with 2 jets, but D* not in!",//3
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge));
const Double_t bins[] = {1.,2.,3.,4.,5.,6.,7.,8.,10.,12.,15.,20.};
fJetMultHistsS83->Add(new TH2F("nDaughtDsJet","N(daughters) in D*-jet",//4
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge,
sizeof(bins)/sizeof(bins[0])-1,bins));
fJetMultHistsS83->Add(new TH2F("ptDsByPtJet","p_{t}(D*)/p_{t}(jet) in D*-jet",//5
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 10, 0., 1.));
fJetMultHistsS83->Add(new TH2F("pDsByPJet","p(D*)/p(jet) in D*-jet",//6
kDmHistNumBins, kDmHistLowEdge, kDmHistHighEdge, 10, 0., 1.));
GFHistManip::AddToHistsName("S83", fJetMultHistsS83);
fAllHistArrays.Add(fJetMultHistsS83);
// if you have more GFHistArray's, don't forget to add them in
// GFDstarAnalyseMc::CreateHistsJetMult()!
oldDir->cd();
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsJetMult(const H1PartDstar* ds, Bool_t accessPartCand)
{
if(!this->IsS83(ds)) return;
if(!fJetMultHistsS83) this->CreateHistsJetMult();
TObjArray jets = this->GetJets(ds);
const Int_t numJets = jets.GetEntriesFast();
TObjArray jetsNonDstar;
this->GetNonDstarJets(ds, jetsNonDstar);
const Int_t numJetsNonDstar = jetsNonDstar.GetEntriesFast();
const Double_t w83 = this->GetWeight(83, ds);
const Double_t dm = ds->GetDm();
static_cast<TH2*>(fJetMultHistsS83->At(0))->Fill(dm, numJets, w83);
static_cast<TH2*>(fJetMultHistsS83->At(1))->Fill(dm, numJetsNonDstar, w83);
if(numJetsNonDstar > 0){
TH2 *h2 = static_cast<TH2*>(fJetMultHistsS83->At(2));
if(numJets < 2) h2->Fill(dm, 0);
else {
h2->Fill(dm, 1, w83);
this->RemoveLowPt(&jets, 2);
if(!this->JetContains(static_cast<H1PartJet*>(jets[0]), ds)
&& !this->JetContains(static_cast<H1PartJet*>(jets[1]), ds)){
fJetMultHistsS83->At(3)->Fill(dm, w83);
}
if(accessPartCand){
H1PartJet *jet = static_cast<H1PartJet*>(jets[0]);
if(!this->JetContains(jet, ds)){
jet = static_cast<H1PartJet*>(jets[1]);
if(!this->JetContains(jet, ds)) jet = NULL;
}
if(jet){
static_cast<TH2*>(fJetMultHistsS83->At(4))->Fill(dm, jet->GetNumOfParticles(), w83);
static_cast<TH2*>(fJetMultHistsS83->At(5))->Fill(dm, ds->GetPt()/jet->GetPt(), w83);
static_cast<TH2*>(fJetMultHistsS83->At(6))->Fill(dm, ds->GetP()/jet->GetP(), w83);
}
}
}
}
}
//_____________________________________________________
void GFDstarAnalyse::CreateHistsTest()
{
return;
// TDirectory* oldDir = this->MakeHistDir("testHists", "test");
// GFHistManip::ClearCreateArray(fTestHists);
// fTestHists->Add(new TH1F("mDstar", "m(K#pi#pi_{s})", 20, 1.8, 2.2));
// fTestHists->Add(new TH1F("mDstarS83", "m(K#pi#pi_{s})", 20, 1.8, 2.2));
// fTestHists->Add(new TH1F("mDstarS84", "m(K#pi#pi_{s})", 20, 1.8, 2.2));
// oldDir->cd();
// fAllHistArrays.Add(fTestHists);
}
//_____________________________________________________
void GFDstarAnalyse::FillHistsTest(const H1PartDstar* ds)
{
// if(!fTestHists) this->CreateHistsTest();
// if(this->DstarCutLoose(ds)){
// const Double_t dm = ds->GetDm();
// if(dm >= fSteer->GetD0DeltaMLow() && dm < fSteer->GetD0DeltaMHigh()){
// const Double_t m = ds->GetMass();
// fTestHists->At(0)->Fill(m);
// if(this->IsS83(ds)) fTestHists->At(1)->Fill(m);
// if(this->IsS84(ds)) fTestHists->At(2)->Fill(m);
// }
// }
}
//_____________________________________________________
//_____________________________________________________
//_____________________________________________________
//_____________________________________________________
//_____________________________________________________
TDirectory* GFDstarAnalyse::OpenHistsFile(const char* fileName)
{
// opens file, does file->cd() and return previous (!) gDirectory
// please do '->cd()' on return value not to change gDirectory!
TDirectory* oldDir = gDirectory;
if(!fHistFile){
TString name(fileName ? fileName :
(("DstarAnalysisHists"+ *fPeriod)+=".root").Data());
fHistFile = new TFile(name,"RECREATE");
} else if(fileName){
this->Warning("OpenHistsFile",
"Already opened (name %s), new name (%s) ignored",
fHistFile->GetName(), fileName);
}
fHistFile->cd();
return oldDir;
}
//_____________________________________________________
TDirectory* GFDstarAnalyse::MakeHistDir(const char* name, const char* title)
{
// make new directory in hist output file and cd into it
// the old directory is returned
TDirectory* oldDir = this->OpenHistsFile();
TDirectory* dir = new TDirectory(name, title);
dir->cd();
return oldDir;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsSignalRegionDm(const H1PartDstar* ds) const
{
// true if D* lies in SR of delta m distribution
const Double_t dm = (ds ? ds->GetDm() : -1.);
if(dm >= fSteer->GetD0DeltaMHigh() || dm < fSteer->GetD0DeltaMLow()){
return kFALSE;
}
return kTRUE;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::DstarCutLoose(const H1PartDstar* dstar) const
{
// all cuts common to both taggers (well, not really all...)
return (this->DstarCutLooseNoTrackCut(dstar)
&& this->AreTracksOk(dstar)
);
}
//_____________________________________________________
Bool_t GFDstarAnalyse::DstarCutLooseNoTrackCut(const H1PartDstar* dstar) const
{
// if(dstar->GetCharge() > 0.) return kFALSE; // only negative D*
const Double_t pt = dstar->GetPt();
if(pt <= fDsCutPt1) return kFALSE;
const Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
const Double_t md0 = dstar->GetD0Mass();
if(pt < 6.){
if(TMath::Abs(md0 - mD0PDG) > fSteer->GetDstarD0Window()) return kFALSE;
} else if(pt < 8.5){
if(TMath::Abs(md0 - mD0PDG) > 0.100) return kFALSE;
} else if(TMath::Abs(md0 - mD0PDG) > 0.140){
return kFALSE;
}
if(dstar->GetDm() >= kDmHistHighEdge) return kFALSE;
if(!this->IsSumPtKPiOK(dstar, fSteer->GetPtKPiCut())) return kFALSE;
const TArrayD* eta33bins = fSteer->GetEtaBins33();
const TArrayD* eta44bins = fSteer->GetEtaBins44();
const Double_t eta = dstar->GetEta();
if(!(eta >= eta33bins->At(0) && eta < eta33bins->At(eta33bins->GetSize()-1)) &&
!(eta >= eta44bins->At(0) && eta < eta44bins->At(eta44bins->GetSize()-1))) return kFALSE;
if(!this->IsSumEt(dstar)) return kFALSE;
if(!this->IsZvtxCut()) return kFALSE;
return kTRUE;
}
// #include "H1Mods/H1PartMCArrayPtr.h"
// #include "H1Mods/H1PartMC.h"
// #include "H1UserDstar/H1UserDstarUtil.h"
//_____________________________________________________
Bool_t GFDstarAnalyse::D0CutLoose(const H1PartDstar* dstar) const
{
// static Bool_t first = kTRUE;
// if(first){
// this->Warning("D0CutLoose", "manipulated to require generated D* -> K pi pi_s");
// first = kFALSE;
// }
// DO NOT FORGET to change also ...Mc::DstarCutLoose!!
// Bool_t genDs = kFALSE;
// static H1PartMCArrayPtr mcs;
// for(Int_t i = 0; i < mcs.GetEntries(); ++i){
// if(H1UserDstarUtil::IsDstarKPiPis(mcs[i], mcs.GetArray())){
// // with p_t(D*) > 1.0 GeV and 15 < theta(D*) < 165 degrees
// if(mcs[i]->GetPt() > 1.
// && mcs[i]->GetTheta() > 15.*TMath::DegToRad()
// && mcs[i]->GetTheta() < 165.*TMath::DegToRad()){
// genDs = kTRUE;
// break;
// }
// }
// }
// if(!genDs) return kFALSE;
const TArrayD* eta44bins = fSteer->GetEtaBins44();
const TArrayD* eta33bins = fSteer->GetEtaBins33();
const Double_t eta = dstar->GetEta();
if(!this->IsSumPtKPiOK(dstar, fSteer->GetPtKPiCut())) return kFALSE;
return (this->IsSignalRegionDm(dstar)
&& dstar->GetPt() > fDsCutPt1
&& ((eta > eta33bins->At(0) && eta < eta33bins->At(eta33bins->GetSize()-1))
|| (eta > eta44bins->At(0) && eta < eta44bins->At(eta44bins->GetSize()-1)))
&& this->AreTracksOk(dstar)
&& this->IsSumEt(dstar)
&& this->IsZvtxCut()
);
}
//_____________________________________________________
Bool_t GFDstarAnalyse::AreTracksOk(const H1PartDstar* dstar) const
{
// all cuts on single track quantities
// cut on radial length (K, pi) and likelihoods (K, pi, pi_s)
if(!dstar) return kFALSE;
const H1PartSelTrack* kaon = dstar->GetKaon();
const H1PartSelTrack* pion = dstar->GetPion();
const H1PartSelTrack* slowPion = dstar->GetSlowPion();
if(kaon->GetPt()<=fSteer->GetPtCutK() || pion->GetPt()<=fSteer->GetPtCutPi()) return kFALSE;
if(!(kaon->IsCentralTrk() && pion->IsCentralTrk() && slowPion->IsCentralTrk())) return kFALSE;
if(kaon->GetRadLength() <= fSteer->GetLengthCutKPi() ||
pion->GetRadLength() <= fSteer->GetLengthCutKPi() ||
slowPion->GetRadLength() <= fSteer->GetLengthCutPis()) return kFALSE;
// static Bool_t first = kTRUE;
// if (first) {
// first = kFALSE;
// this->Info("AreTracksOk", "Request one daughter with chi^2/ndf > 50");
// }
// const H1VertexFittedTrack *vTrkK =
// static_cast<const H1VertexFittedTrack*>(kaon->GetParticle()->GetTrack());
// const H1VertexFittedTrack *vTrkPi =
// static_cast<const H1VertexFittedTrack*>(pion->GetParticle()->GetTrack());
// const H1VertexFittedTrack *vTrkPis =
// static_cast<const H1VertexFittedTrack*>(slowPion->GetParticle()->GetTrack());
// if (vTrkK ->GetFitChi2()/vTrkK ->GetFitNdf() < 50. &&
// vTrkPi ->GetFitChi2()/vTrkPi ->GetFitNdf() < 50. &&
// vTrkPis->GetFitChi2()/vTrkPis->GetFitNdf() < 50.) {
// return kFALSE;
// }
return this->IsDedxOK(dstar);
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsDedxOK(const H1PartDstar* dstar, Int_t nDaughter) const
{
// check whether dedx is OK:
// nDaughter < 0: check all K, pi, pis (default)
// = 0: check K only
// = 1: check pi only
// = 2: check pi_s only
if(fIgnoreDedx) return kTRUE;
Float_t likeliHoodK = -1., likeliHoodPi = -1.;//, likeliHoodPis = -1.;
// if(dstar->IsWrongCharge()){
// const Int_t index = H1UserDstarAccess::ArrayIndexWc(dstar);
// // if(dstar->GetKaon()->GetP() < fSteer->GetPMaxForDedxCut())
// likeliHoodK = H1UserDstarAccess::TrackLHKofKWc(index);
// // if(dstar->GetPion()->GetP() < fSteer->GetPMaxForDedxCut())
// likeliHoodPi = H1UserDstarAccess::TrackLHPiofPiWc(index);
// // if(dstar->GetSlowPion()->GetP() < fSteer->GetPMaxForDedxCut())
// // likeliHoodPis = H1UserDstarAccess::TrackLHPiofPisWc(index);
// } else {
// const Int_t index = H1UserDstarAccess::ArrayIndex(dstar);
// // if(dstar->GetKaon()->GetP() < fSteer->GetPMaxForDedxCut())
// likeliHoodK = H1UserDstarAccess::TrackLHKofK(index);
// // if(dstar->GetPion()->GetP() < fSteer->GetPMaxForDedxCut())
// likeliHoodPi = H1UserDstarAccess::TrackLHPiofPi(index);
// // if(dstar->GetSlowPion()->GetP() < fSteer->GetPMaxForDedxCut())
// // likeliHoodPis = H1UserDstarAccess::TrackLHPiofPis(index);
// }
likeliHoodK = dstar->GetKaon()->DedxLikelihood(H1DedxNLH::kKaon);
likeliHoodPi = dstar->GetPion()->DedxLikelihood(H1DedxNLH::kPion);
if(nDaughter <= 0){
if(likeliHoodK >= 0. && likeliHoodK < fSteer->GetLikeliHoodCut()
// && this->GetNLH(dstar, 0) < fSteer->GetNLHCut()
) {
if(dstar->GetKaon()->GetP() < fSteer->GetPMinSpecialLHK() ||
dstar->GetKaon()->GetP() > fSteer->GetPMaxSpecialLHK() ||
likeliHoodK < fSteer->GetSpecialLHCutK()) {
return kFALSE;
}
}
}
if(nDaughter < 0 || nDaughter == 1){
if(likeliHoodPi >= 0. && likeliHoodPi < fSteer->GetLikeliHoodCut()
// && this->GetNLH(dstar, 1) < fSteer->GetNLHCut()
) return kFALSE;
}
if(nDaughter < 0 || nDaughter == 2){
H1DedxNLH nlh;
if(dstar->GetSlowPion()->GetDedx() > 0. && // no kFALSE for no measurement
dstar->GetSlowPion()->GetDedx() > nlh.RefDedxPion(dstar->GetSlowPion()->GetP())
+ fSteer->GetDedxOffsetCutPis()
// if(likeliHoodPis >= 0. && likeliHoodPis < fSteer->GetLikeliHoodCut()
// // && this->GetNLH(dstar, 2) < fSteer->GetNLHCut()
) return kFALSE;
}
return kTRUE;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsSumPtKPiOK(const H1PartDstar* ds, Double_t ptCut) const
{
if(fIgnoreSumPtKPi) return kTRUE;
if(!ds) return kFALSE;
return (ds->GetKaon()->GetPt() + ds->GetPion()->GetPt() > ptCut);
}
//_____________________________________________________
Float_t GFDstarAnalyse::GetNLH(const H1PartDstar* dstar, Int_t part) const
{
// get normalised likelihood taking into account K, pi and proton for a D* daughter:
// part: 0 = K, 1 = pi, 2 = pis (particle fluxes of pi, K and p assumed to be equal)
// 1.09 if no measurement available, FLT_MAX if problems...
const H1PartSelTrack* track = NULL;
switch(part){
case 0: // kaon
track = dstar->GetKaon();
break;
case 1: // pion
track = dstar->GetPion();
break;
case 2: // slow pion
track = dstar->GetSlowPion();
break;
default:
this->Error("GetNLH", "D* has no particle %d", part);
return FLT_MAX;
}
const Float_t dedx = track->GetDedx();
const Short_t nHitDedx = track->GetNHitDedx();
if(nHitDedx < 0 || dedx < 0.) return 1.09; // no measurement
const Float_t deltaP = track->GetDp();
const Double_t mom = track->GetP();
H1Dedx* dedxDB = static_cast<H1Dedx*>(H1DBManager::Instance()->GetDBEntry(H1Dedx::Class()));
const Float_t relDedxErr =
dedxDB->GetDedxResolution(nHitDedx, track->GetPhi(), track->GetCharge());
const Float_t fluxes[] = {0., 1., 1., 1., 0.};//mu, K, pi, proton, deuteron
const H1DedxNLH dedxLikeli(TArrayF(sizeof(fluxes)/sizeof(fluxes[0]), fluxes));
switch(part){
case 0: // kaon
return dedxLikeli.GetNLHKaon(dedx, dedx * relDedxErr, mom, deltaP);
case 1: // pion
case 2: // or slow pion
return dedxLikeli.GetNLHPion(dedx, dedx * relDedxErr, mom, deltaP);
}
return FLT_MAX; // unreached
}
//________________________________________________
void GFDstarAnalyse::Print(const H1PartDstar* dstar) const
{
dstar->Print();
cout << "delta M = " << dstar->GetDm() << endl;
const H1PartSelTrack* k = dstar->GetKaon();
const H1PartSelTrack* pi = dstar->GetPion();
const H1PartSelTrack* pis = dstar->GetSlowPion();
cout << "KPiPis:t\tpttphi\ttheta"
<< "nKaon:t\t" << k->GetPt() << "t" << k->GetPhi() << "t" <<k->GetTheta()
<< "nPion:t\t" << pi->GetPt() << "t" << pi->GetPhi() << "t" <<pi->GetTheta()
<< "nslowPion:t" << pis->GetPt() << "t" << pis->GetPhi() << "t" <<pis->GetTheta()
<< endl;
}
//________________________________________________
Double_t GFDstarAnalyse::CalcSumEt(Double_t thetaCutDegree) const
{
static H1PartCandArrayPtr cands;
Double_t sumEt = 0.;
Int_t numCands = cands.GetEntries();
for(Int_t i = 0; i < numCands; ++i){
if(cands[i]->GetTheta() > thetaCutDegree/180.*TMath::Pi())
sumEt += cands[i]->GetE() * TMath::Sin(cands[i]->GetTheta());//->GetPt();
}
return sumEt;
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetCmEnergySqr() const
{
// square of center of mass energy (e+p) [i.e. s ~= 4 * E_e * E_p]
// if H1Tree exists first try dstarUser tree for beam energies
// if nothing helps: take hardcoded default
if(gH1Tree){
static H1FloatPtr beamProtE("BeamEnergyProt");
static H1FloatPtr beamElecE("BeamEnergyElec");
return 4.* *beamProtE * *beamElecE;
}
static Bool_t first = kTRUE;
if(first){
this->Warning("GetCmEnergySqr",
"no H1Tree, using default result %f",
kCmEnergySqr);
first = kFALSE;
}
return kCmEnergySqr;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsL4Found(Int_t hqsel45Bit, Int_t hqselMode,
const H1BankPtr& yrclYecl) const
{
// true if hqsel45Bit is on in
// word hqselMode (mode 0: DIS, 1: very low Q2, 2: tagged, 3: untagged photoproduction,
// -1: any of them, -2: just check L4-class 15, ignore hqsel45Bit and yrclYecl)
// of given bank yrclYecl (which is supposed to be a YECL [online] or YRCL [offline] bank)
if(hqselMode == -2){
static H1BytePtr l4class("Iclass");
return (l4class[15] != 0);
}
if(yrclYecl.GetEntries() > 0){
if(hqselMode < -1 || hqselMode > 3){
this->Error("IsL4Found", "hqselMode %d not known", hqselMode);
return kFALSE;
}
if(hqselMode == -1){
for(Int_t i = 0; i < 3; ++i){
const Int_t wordInt = yrclYecl[0]->GetInt(i);
if((wordInt & (1 << hqsel45Bit))) return kTRUE;
}
return kFALSE;
}
const Int_t word = yrclYecl[0]->GetInt(hqselMode);
bool result = (word & (1 << hqsel45Bit)); // before 3.03_09 not Bool_t: (Bool_t) 2048 = 0!
return result;
} else { // use for begin 99 e+:
static Int_t numWarn = 0;
if(numWarn < kNumMaxWarn) {
if(numWarn == kNumMaxWarn-1) {
this->Warning("IsL4Found", "Following warning will not be shown further!");
}
this->Warning("IsL4Found", "no %s bank in run/event %d / %d, return false",
yrclYecl.GetBranch() ? yrclYecl.GetBranch()->GetName() : "branchless" ,
gH1Tree->GetRunNumber(), gH1Tree->GetEventNumber());
}
++numWarn;
return kFALSE;
}
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsZvtxCut(Float_t cut) const // 0.: default!
{
// cutting on abs(zVtx) < cut
// cut == 0 (default) cut is set to fZvtxCut
if(cut == 0.) cut = fSteer->GetZvtxCut();
static H1FloatPtr zVtx("VtxZ");
return (TMath::Abs(*zVtx) < cut);
}
//_____________________________________________________
Bool_t GFDstarAnalyse::IsSumEt(const H1PartDstar* ds) const
{
if(fIgnoreSumEt) return kTRUE;
if(!ds) return kFALSE;
return ((ds->GetPt() / this->GetSumEt()) > fSteer->GetSumEtCut());
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetSumEt(Int_t threshold) const
{
// if(!fOpenUserTree) return 0.;
static H1FloatPtr sumEt("SumEt");
return sumEt[threshold];
// static H1FloatPtr sumTrackPt("SumTrackPt");
// return *sumTrackPt;
}
//_____________________________________________________
Float_t GFDstarAnalyse::GetYjb() const
{
static H1FloatPtr yjb("Yjb");
return *yjb;
}
//_____________________________________________________
Float_t GFDstarAnalyse::GetY(Int_t subTrig) const
{
switch(subTrig){
case 83:
return this->GetY33();
case 84:
return this->GetY44();
default:
return this->GetYjb();
}
}
//_____________________________________________________
Float_t GFDstarAnalyse::GetWgammaPjb() const
{
static H1FloatPtr wjb("WgammaPjb");
return *wjb;
}
//_____________________________________________________
Float_t GFDstarAnalyse::GetWgammaP(Int_t subTrig) const
{
switch(subTrig){
case 83:
return this->GetWgammaP33();
case 84:
return this->GetWgammaP44();
default:
return this->GetWgammaPjb();
}
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetZ(const H1PartDstar *ds) const
{
static H1FloatPtr beamElecE("BeamEnergyElec");
// const Double_t empzDstar = H1UserDstarAccess::EminusPz(ds);
// const Double_t empzDstar2 = ds->GetE() - ds->GetPz();
// TLorentzVector dsHfs = ds->GetKaon()->GetParticle()->GetFourVector();
// dsHfs += ds->GetPion()->GetParticle()->GetFourVector();
// dsHfs += ds->GetSlowPion()->GetParticle()->GetFourVector();
// const Double_t empzDstar4 = dsHfs.E() - dsHfs.Pz();
// TLorentzVector dsPionTrack = ds->GetKaon()->GetFourVector();
// dsPionTrack += ds->GetPion()->GetFourVector();
// dsPionTrack += ds->GetSlowPion()->GetFourVector();
// const Double_t empzDstar3 = dsPionTrack.E() - dsPionTrack.Pz();
// if(!GFMath::IsEqualFloat(empzDstar, empzDstar3)){
// cout << "n**************************************n"
// << "E-pz(D*) HFS: " << empzDstar
// << "n tracks: " << empzDstar3
// << "n D*: " << empzDstar2
// << "n HFS new: " << empzDstar4
// << endl;
// ds->GetKaon()->GetParticle()->Print();
// ds->GetPion()->GetParticle()->Print();
// ds->GetSlowPion()->GetParticle()->Print();
// }
// we should use the D* instead of its HFS:
// - There is a GetMass()-bug in H1PartCand from tracks in 2.6.X
// - Even without a bug we neglect the kaon mass compared to D* mass.
const Double_t empzDstar = ds->GetE() - ds->GetPz();
// FIXME: S84???
const Double_t zDs = empzDstar / (2. * this->GetY33() * (*beamElecE));
return zDs;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::OneL4ClassGiven(const TArrayI& classes) const
{
static H1BytePtr eventClass("Iclass");
for(Int_t i = 0; i < classes.GetSize(); ++i){
if(eventClass[classes.At(i)] != 0) return kTRUE;
}
return kFALSE;
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetXgamma(const H1Part* part1, const H1Part* part2,
Double_t y) const
{
// x_gamma from pt, eta, y
Double_t xga = part1->GetPt() * TMath::Exp(-part1->GetEta())
+ part2->GetPt() * TMath::Exp(-part2->GetEta());
static H1FloatPtr elecBeamE("BeamEnergyElec");
return (y ? xga / (2. * *elecBeamE * y) : -111.);
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetXgamma(const H1Part* part1, const H1Part* part2,
Int_t trig) const
{
Double_t y = 0.;
switch(trig){
case 83:
y = this->GetY33();
break;
case 84:
y = this->GetY44();
break;
default:
y = this->GetYjb();
}
return this->GetXgamma(part1, part2, y);
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetXgammaMass(const H1PartDstar* dstar,const H1PartJet* jet,
Double_t yHfs) const
{
// xgamma taking into account masses:
// [(E-pz)_{D*_HFS} + sum_{i=1}^{N(jet daught)}(E-pz)_{i}]/2*yHfs*E_{e}
if(!dstar || !jet) return -1.;
Float_t eMinPzJet = 0.;
switch(fJetFlag){
case 0:
eMinPzJet = H1UserDstarAccess::EminusPz(jet, dstar);
break;
case 1:
eMinPzJet = jet->GetE() - jet->GetPz();
break;
default:
this->Error("GetXgammaMass", "Do not know jet flag %d", fJetFlag);
}
#ifdef H1OO25
const Float_t eMinPzDs = dstar->GetE() - dstar->GetMomentum().Pz();
#else
// dstar E-pz from its particle candidates to remove biases:
const Float_t eMinPzDs = H1UserDstarAccess::EminusPz(dstar);
// const TLorentzVector pDstarHFS = dstar->GetKaon()->GetParticle()->GetFourVector()
// + dstar->GetPion()->GetParticle()->GetFourVector()
// + dstar->GetSlowPion()->GetParticle()->GetFourVector();
// if(TMath::Abs((dstar->H1PartComp::GetFourVector().E() - dstar->H1PartComp::GetFourVector().Pz()) - eMinPzDs) > 0.1){
// dstar->Print("all");
// dstar->GetKaon()->GetParticle()->Print();
// dstar->GetPion()->GetParticle()->Print();
// dstar->GetSlowPion()->GetParticle()->Print();
// cout << "E-pz D*, tracks " << dstar->GetE() - dstar->GetMomentum().Pz()
// << ", HFS " << eMinPzDs
// << ", ratio " << (dstar->GetE() - dstar->GetMomentum().Pz()) / eMinPzDs
// << endl;
// cout << "E-pz D*, tracks as pi " << (dstar->H1PartComp::GetFourVector().E()
// - dstar->H1PartComp::GetFourVector().Pz())
// << ", ratio " << (dstar->H1PartComp::GetFourVector().E()
// - dstar->H1PartComp::GetFourVector().Pz()) / eMinPzDs
// << "nrun/evt: " << gH1Tree->GetRunNumber() << " / " << gH1Tree->GetEventNumber()
// << endl << endl;
// }
#endif
static H1FloatPtr beamElecE("BeamEnergyElec");
return (eMinPzDs + eMinPzJet)/(2. * yHfs * (*beamElecE));
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetXgammaMassDiJet(const H1PartDstar *dstar,
const H1PartJet *dsJet, const H1PartJet *othJet,
Double_t yHfs) const
{
// xgamma taking into account masses:
// (sum_{i=1}^{N(D*-jet daught)}(E-pz)_{i} + sum_{i=1}^{N(other jet daught)}(E-pz)_{i}]) /
// 2*yHfs*E_{e}
if(!dstar || !dsJet || !othJet) return -1.;
Double_t eMinPzJets = 0.;
switch(fJetFlag){
case 0:
eMinPzJets = H1UserDstarAccess::EminusPz(dsJet, dstar);// D* tracks already replaced by HFS
eMinPzJets += H1UserDstarAccess::EminusPz(othJet, dstar);
break;
case 1:
eMinPzJets = dsJet->GetE() - dsJet->GetPz();
eMinPzJets += othJet->GetE() - othJet->GetPz();
break;
default:
this->Error("GetXgammaMassDiJet", "Do not know jet flag %d", fJetFlag);
}
static H1FloatPtr beamElecE("BeamEnergyElec");
return eMinPzJets/(2. * yHfs * (*beamElecE));
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetXpMass(const H1PartDstar* dstar,const H1PartJet* jet) const
{
// SLOW: looping jet daughters!
// x_p taking into account masses:
// [(E+pz)_{D*_HFS} + sum_{i=1}^{N(jet daught)}(E+pz)_{i}]/2*E_{p}
if(!dstar || !jet) return -1.;
// switch(fJetFlag){
// case 0:
// eMinPzJet = H1UserDstarAccess::EminusPz(jet, dstar);
// break;
// case 1:
// eMinPzJet = jet->GetE() - jet->GetPz();
// break;
// default:
// this->Error("GetXgammaMass", "Do not know jet flag %d", fJetFlag);
// }
Int_t nDaugh = jet->GetNumOfParticles();
Float_t ePlusPz = 0.;
for(Int_t i = 0; i < nDaugh; ++i){
if(!jet->GetParticle(i)->InheritsFrom(H1Part::Class())){
this->Error("GetXpMass", "part %d is %s (r/e=%d%d)", i,
jet->GetParticle(i)->IsA()->GetName(), gH1Tree->GetRunNumber(),
gH1Tree->GetEventNumber());
continue;
}
ePlusPz += jet->GetParticle(i)->GetE();
ePlusPz += jet->GetParticle(i)->GetPz();
}
// dstar E-pz from its particle candidates to remove biases:
// const Float_t eMinPzDs = H1UserDstarAccess::EminusPz(dstar);
for(Int_t i = 0; i < 3; ++i){
ePlusPz += dstar->GetPartCand(i)->GetE();
ePlusPz += dstar->GetPartCand(i)->GetPz();
}
static H1FloatPtr beamProtE("BeamEnergyProt");
return ePlusPz/(2. * *beamProtE);
}
//_____________________________________________________
Int_t GFDstarAnalyse::GetHighestPt(const TObjArray& parts) const
{
// parts must contain H1Part, return index of that with highest pt,
// -1 if 'parts' empty
Int_t index = -1;
Double_t maxPt = -DBL_MAX;
for(Int_t i = 0; i < parts.GetEntriesFast(); ++i){
const H1Part *part = static_cast<H1Part*>(parts[i]);
if(part && maxPt < part->GetPt()){
maxPt = part->GetPt();
index = i;
}
}
return index;
}
//_____________________________________________________
TObjArray GFDstarAnalyse::GetJets(const H1PartDstar* ds, Int_t* dsJet) const
{
// get array of jets, input defined by fJetFlag
// output is ordered in eta: starting with negative eta
TObjArray* jetsArray = NULL;
static H1PartJetArrayPtr inclJets("InclusiveKtJets"); // needed only for no D*
if(ds == NULL){
jetsArray = inclJets.GetArray();
} else {
switch(fJetFlag){
case 0:
jetsArray = H1UserDstarAccess::DstarJets(ds);
break;
case 1:
jetsArray = inclJets.GetArray();
break;
default:
this->Warning("GetJets", "Don't know fJetFlag = %d", fJetFlag);
}
}
if(dsJet) *dsJet = -1;
if(!jetsArray) return TObjArray(0);
TObjArray jets(jetsArray->GetEntriesFast());
TArrayD etas(jetsArray->GetEntriesFast());
// fill array with jets passing cuts
for(Int_t i = 0; i < jetsArray->GetEntriesFast(); ++i){
H1PartJet* jet = static_cast<H1PartJet*>(jetsArray->At(i));
if(this->JetCut(jet)){
etas[jets.GetEntriesFast()] = jet->GetEta(); // before Add(jet)!!
jets.Add(jet);
}
}
// now sort in eta:
etas.Set(jets.GetEntriesFast()); // shorten
TArrayI indices(etas.GetSize());
TMath::Sort(etas.GetSize(), etas.GetArray(), indices.GetArray(), kFALSE); // upwart!
//fill output
TObjArray result(etas.GetSize());
for(Int_t i = 0; i < etas.GetSize(); ++i){
result.Add(jets[indices[i]]);
if(dsJet && this->JetContains(static_cast<H1PartJet*>(result.Last()), ds)){
if(*dsJet != -1) this->Error("GetJets", "2 jets contain D*!");
*dsJet = i;
}
}
return result;
}
//_____________________________________________________
Int_t GFDstarAnalyse::GetNonDstarJets(const H1PartDstar* dstar,
TObjArray& jets) const
{
// array of jets passing the minimal criteria and NOT containing 'dstar'
// will be put into 'jets'
// ordered in eta: most backward first!
//
// returns the index of the jet, that comes after the D*
// if we put the D* into the list, too:
//
//-1: no jets left
// 0: D* more backwards than most backward jet
// 1: D* between most backward and next more forward jet
// 2: etc.
Int_t posDstar = -1;
jets.Clear();
// if(!dstar) return posDstar;
// if you change something here, do it analog in GFDstarAnalyseMc::GetNonDstarGenJets(...) !
Int_t dsJet = 0;
TObjArray allJets = this->GetJets(dstar, &dsJet);
for(Int_t i = 0; i < allJets.GetEntriesFast(); ++i){
if(i != dsJet) {
H1PartJet* jet = static_cast<H1PartJet*>(allJets[i]);
// if(!dstar|| TMath::Abs(jet->GetMomentum().DeltaPhi(dstar->GetMomentum()))>TMath::Pi()/2.){
jets.Add(jet);
// } else {
// static Int_t first = 0;
// static H1PartJet *oldJet = NULL;
// if(jet->GetPt() > 3. && oldJet != jet && first++ < 5){
// this->Info("GetNonDstarJets",
// "exclude jets: delta R(D*+j) %.3f , dphi %.3f, deta %.3f (%d)",
// jet->GetMomentum().DeltaR(dstar->GetMomentum()),
// jet->GetMomentum().DeltaPhi(dstar->GetMomentum()),
// jet->GetEta() - dstar->GetEta(), first);
// cout << "rec D*: ";//theta = " << dstar->GetTheta()*TMath::RadToDeg() << endl;
// dstar->Print();
// cout << "rec jet: theta = " << jet->GetTheta()*TMath::RadToDeg() << endl;
// jet->Print();
// cout << endl << endl;
// oldJet = jet;
// }
// }
}
}
// static UInt_t count = 0;
// if(jets.GetEntriesFast()){
// const Int_t iJetHighPt = this->GetHighestPt(jets);
// const H1PartJet* jetHighPt = static_cast<H1PartJet*>(jets[iJetHighPt]);
// if(jetHighPt->GetMomentum().DeltaR(dstar->GetMomentum()) < 1.5){
// if(count < 10){
// this->Warning("GetNonDstarJets", "empty %d jets due to highest pt is dR=%.3f",
// jets.GetEntriesFast(),
// jetHighPt->GetMomentum().DeltaR(dstar->GetMomentum()));
// ++count;
// }
// jets.Clear();
// return posDstar;
// }
// }
if(jets.GetEntriesFast()) posDstar = 0;
const Double_t etaDs = dstar ? dstar->GetEta() : -1.e7;
for(Int_t i = 0; i < jets.GetEntriesFast(); ++i){
H1PartJet* jet = static_cast<H1PartJet*>(jets[i]);
if(etaDs > jet->GetEta()) ++posDstar;
}
return posDstar;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::GetDiJets(const H1PartDstar* dstar, const H1PartJet *&dsJet,
const H1PartJet *&oJet) const
{
// dsJet and oJet are set (possibly only one of them!)
// return kFALSE if dstar not in dsJet (i.e. dsJet is closer than oJet to D*)
dsJet = oJet = NULL;
//Int_t dsJet = 0;
TObjArray jets = this->GetJets(dstar, NULL);//&dsJet);
const Int_t nJet = this->RemoveLowPt(&jets, 2);
if(nJet == 0) return kFALSE;
H1PartJet *jet1 = static_cast<H1PartJet*>(jets[0]);
H1PartJet *jet2 = (nJet > 1 ? static_cast<H1PartJet*>(jets[1]) : NULL);
Bool_t realDsJet = kTRUE;
if(this->JetContains(jet1, dstar)){
dsJet = jet1;
oJet = jet2;
} else if(jet2 && this->JetContains(jet2, dstar)){
dsJet = jet2;
oJet = jet1;
} else {
realDsJet = kFALSE;
if(!jet2){
dsJet = jet2;// jet2 == NULL anyway
oJet = jet1;
} else {
// this->Warning("GetDiJets", "D* not in first two jets");
if(jet1->GetMomentum().DeltaR(dstar->GetMomentum())
< jet2->GetMomentum().DeltaR(dstar->GetMomentum())){
dsJet = jet1;
oJet = jet2;
} else {
dsJet = jet2;
oJet = jet1;
}
}
}
if ((dsJet && dsJet->GetPt() > fSteer->GetDstar1JetPtBins()->At(0) + 1.)
|| (oJet && oJet->GetPt() > fSteer->GetDstar1JetPtBins()->At(0) + 1.)) {
return realDsJet;
} else {
dsJet = oJet = NULL;
return kFALSE;
}
}
//_____________________________________________________
Bool_t GFDstarAnalyse::JetCut(const H1PartJet* jet) const//, Bool_t firstJet) const
{
// true if jet survives cuts, if 'firstJet': require higher pt-cut
// only pt and eta cuts
if(!jet) return kFALSE;
const Double_t eta = jet->GetEta();
if(eta < fSteer->GetDstar1JetEtaBins()->At(0)
|| eta > fSteer->GetDstar1JetEtaBins()->At(fSteer->GetDstar1JetEtaBins()->GetSize()-1)){
return kFALSE;
}
const Double_t pt = jet->GetPt();
// if(firstJet && pt < fSteer->GetJetHigherPtCut()) return kFALSE;
if(pt < fSteer->GetDstar1JetPtBins()->At(0)) return kFALSE;
if(jet->GetNumOfParticles()){
Bool_t recoJet = kFALSE;
for(Int_t i = 0; i < jet->GetNumOfParticles(); ++i){
const H1Part *part = jet->GetParticle(i);
if(part->InheritsFrom(H1PartCand::Class()) || part->InheritsFrom(H1PartDstar::Class())){
recoJet = kTRUE;
break;
}
}
if(recoJet && jet->GetNumOfParticles() == 1
&& !jet->GetParticle(0)->InheritsFrom(H1PartDstar::Class())){
// static Int_t nWarn = 0;
// if(nWarn < 10){
// nWarn++;
// this->Info("JetCut", "removal of single particle jets (%d/10)", nWarn);
// }
return kFALSE;
}
} else {
this->Error("JetCut", "jet without particles!");
}
return kTRUE;
}
//_____________________________________________________
Bool_t GFDstarAnalyse::JetContains(const H1PartJet* jet, const H1PartDstar* ds) const
{
// true if jet contains ds
#ifdef H1OO25
return (jet && jet->Contains(ds));
#else
if(!jet || !ds){
this->Warning("JetContains", "Jet (%p) or D* (%p) not given", jet, ds);
return kFALSE;
}
const Int_t dstarJet = H1UserDstarAccess::JetIndex(ds);
if(dstarJet < 0){
return kFALSE;
} else {
const TObjArray* jets = H1UserDstarAccess::DstarJets(ds);
if(!jets){
this->Error("JetContains", "no jets found for D*");
return kFALSE;
}
return (jet == const_cast<const TObject*>(jets->At(dstarJet)));
}
#endif
}
//_____________________________________________________
Double_t GFDstarAnalyse::GetRecEnergies(Double_t& forward, Double_t& forwardBar,
Double_t& barrel, Double_t& spacal) const
{
// fills energie sums in detector region defined by some angle (cf.code)
// returns sum of all regions
// angles should be the same as in GFDstarAnalyseMc::GetGenEnergies(...)!!!
const Double_t thetaSpac = 177. * TMath::DegToRad();
const Double_t thetaSpacBar = 153.5 * TMath::DegToRad();
const Double_t thetaBarForwBar = 55. * TMath::DegToRad();
const Double_t thetaForwBarForw = 30. * TMath::DegToRad();
// static Bool_t first = kTRUE;
// if(first){
// first = kFALSE;
// this->Info("GetRecEnergies", "manipulated angles of detector regions");
// }
// const Double_t thetaSpacBar = 150.0 * TMath::DegToRad();
// const Double_t thetaBarForwBar = 30. * TMath::DegToRad();
// const Double_t thetaForwBarForw = 15. * TMath::DegToRad();
const Double_t thetaForw = 4.5 * TMath::DegToRad();
forward = forwardBar = barrel = spacal = 0.;
static H1PartCandArrayPtr part;
for(Int_t i = 0; i < part.GetEntries(); ++i){
if(part[i]->GetP() < 1.e-10) continue; // some tracks discarded by HFS, some muons
// if(!TMath::Abs(part[i]->GetCharge())) continue;
const Double_t theta = part[i]->GetTheta();
const Double_t en = part[i]->GetE();
if(theta < thetaForw) continue;
else if(theta < thetaForwBarForw){
forward += en;
} else if(theta < thetaBarForwBar){
forwardBar += en;
} else if(theta < thetaSpacBar){
barrel += en;
} else if(theta < thetaSpac){
// if(part[i]->IsHFSChargedCls() || part[i]->IsHFSUnidentifiedCls()){
// spacal += (1.65 * en);
// } else {
spacal += en;
// }
}
}
return forward + forwardBar + barrel + spacal;
}
//_____________________________________________________
void GFDstarAnalyse::GetEtRatios(const H1PartJet* jet, Double_t& ratioTrack,
Double_t& ratioEmClus, Double_t& ratioHadClus)
{
ratioTrack = ratioEmClus = ratioHadClus = 0.;
for(Int_t n = 0; n < jet->GetNumOfParticles(); ++n){
const H1Part* p = jet->GetParticle(n);
if(p->InheritsFrom(H1PartCand::Class())){
const H1PartCand* pCand = static_cast<const H1PartCand*>(p);
if(pCand->IsHFSChargedGoodTrk() || pCand->IsHFSChargedTrk() || pCand->IsMuon()){
ratioTrack += pCand->GetPt();
} else if(pCand->IsHFSEm() || pCand->IsEm()){
ratioEmClus += pCand->GetPt();
} else {
ratioHadClus += pCand->GetPt();
}
} else if(p->InheritsFrom(H1PartDstar::Class()) ||
p->InheritsFrom(H1PartSelTrack::Class())){// in case we build jets from sel tracks
ratioTrack += p->GetPt();
} else {
this->Warning("GetEtRatios", "particle neither H1PartCand nor H1PartDstar!");
}
}
const Double_t jetPt = jet->GetPt();
if(TMath::Abs(ratioTrack + ratioEmClus + ratioHadClus - jetPt) > 1.e-5){
this->Info("GetEtRatios", "ptTrack + ptEmClus + ptHadClus = %f, jetPt = %f (r/e=%d/%d)",
ratioTrack + ratioEmClus + ratioHadClus, jetPt, gH1Tree->GetRunNumber(),
gH1Tree->GetEventNumber());
}
ratioTrack /= jetPt;
ratioEmClus /= jetPt;
ratioHadClus /= jetPt;
}
//_____________________________________________________
const H1Part* GFDstarAnalyse::GetCloser(const H1Part* reference, const H1Part* first,
const H1Part* second, Float_t* dist)
{
// static method
// returns that H1Part from first/second that is closer in eta phi to reference,
// if dist points to any memory it is filled with the distance
if(!reference) return NULL;
const TVector3 ref = reference->GetMomentum();
const Double_t dist1 = (first ? ref.DeltaR(first->GetMomentum()) : DBL_MAX);
const Double_t dist2 = (second ? ref.DeltaR(second->GetMomentum()) : DBL_MAX);
if(dist1 > dist2){
if(dist) *dist = dist2;
return second;
} else {
if(dist) *dist = dist1;
return first;
}
}
//_____________________________________________________
Int_t GFDstarAnalyse::RemoveLowPt(TObjArray* parts, Int_t nStay) const
{
// remove the lowest pt particles in 'parts' such that 'nStay' are left
// does not change the order
// return number of actually particles left
const Int_t nInput = (parts ? parts->GetEntriesFast() : 0);
if(nInput <= nStay) return nInput; // nothing to do
if(nStay < 0) { // check for nonsense input
if(parts) parts->Clear();
return 0;
}
TArrayD pts(nInput);
for(Int_t i = 0; i < nInput; ++i){
pts[i] = static_cast<H1Part*>(parts->At(i))->GetPt();
}
TArrayI indices(pts.GetSize());
TMath::Sort(pts.GetSize(), pts.GetArray(), indices.GetArray());
for(Int_t i = nStay; i < pts.GetSize(); ++i){
parts->RemoveAt(indices[i]);
}
parts->Compress();
return parts->GetEntriesFast();
}
//_____________________________________________________
Int_t GFDstarAnalyse::RemoveLowPt(TObjArray* parts, Double_t minPt) const
{
// remove particles in 'parts' that have pt lower than 'minPt'
// does not change the order
// return number of actually particles left
const Int_t nInput = (parts ? parts->GetEntriesFast() : 0);
for(Int_t i = 0; i < nInput; ++i){
H1Part* part = static_cast<H1Part*>(parts->At(i));
if(part->GetPt() < minPt) parts->RemoveAt(i);
}
parts->Compress();
return parts->GetEntriesFast();
}
//_____________________________________________________
Double_t GFDstarAnalyse::Associate(const H1Part* ref1, const H1Part* ref2,
H1PartJet*& outCloseTo1, H1PartJet*& outCloseTo2)
{
// taking outCloseTo1/2 as input, which combination
// (ref1 <-> outCloseTo1 + ref2 <-> outCloseTo2)
// or (ref1 <-> outCloseTo2 + ref2 <-> outCloseTo1)
// gives the minimal sum of eta-phi distances.
// outCloseTo1/2 are changed that 1 is the one for ref1, analog for2,
// return sum of distances.
const TVector3 ref1Mom = ref1->GetMomentum();
const TVector3 ref2Mom = ref2->GetMomentum();
const TVector3 out1Mom = outCloseTo1->GetMomentum();
const TVector3 out2Mom = outCloseTo2->GetMomentum();
const Double_t dist11_22 = ref1Mom.DeltaR(out1Mom) + ref2Mom.DeltaR(out2Mom);
const Double_t dist12_21 = ref1Mom.DeltaR(out2Mom) + ref2Mom.DeltaR(out1Mom);
if(dist11_22 <= dist12_21){
return dist11_22;
} else {
H1PartJet* temp = outCloseTo1;
outCloseTo1 = outCloseTo2;
outCloseTo2 = temp;
return dist12_21;
}
}
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.