#include <AliLog.h>
#include "AliCorrectionMatrix3D.h"
+#include "AliPWG0Helper.h"
//____________________________________________________________________
ClassImp(AliCorrectionMatrix3D)
AliCorrectionMatrix::SaveHistograms();
- WriteProjections(GetMeasuredHistogram());
- WriteProjections(GetGeneratedHistogram());
+ AliPWG0Helper::CreateProjections(GetMeasuredHistogram());
+ AliPWG0Helper::CreateProjections(GetGeneratedHistogram());
if (GetCorrectionHistogram())
- WriteProjections(GetCorrectionHistogram());
-}
-
-//____________________________________________________________________
-void AliCorrectionMatrix3D::WriteProjections(TH3F* hist)
-{
- // write some projections to disk
-
- TH1* proj = hist->Project3D("yx");
- proj->SetXTitle(hist->GetXaxis()->GetTitle());
- proj->SetYTitle(hist->GetYaxis()->GetTitle());
-
- proj = hist->Project3D("zx");
- proj->SetXTitle(hist->GetXaxis()->GetTitle());
- proj->SetYTitle(hist->GetZaxis()->GetTitle());
-
- proj = hist->Project3D("zy");
- proj->SetXTitle(hist->GetYaxis()->GetTitle());
- proj->SetYTitle(hist->GetZaxis()->GetTitle());
+ AliPWG0Helper::CreateProjections(GetCorrectionHistogram());
}
void RemoveEdges(Float_t cut=2, Int_t nBinsXedge = 0, Int_t nBinsYedge = 0, Int_t nBinsZedge = 0);
virtual void SaveHistograms();
- void WriteProjections(TH3F* hist);
protected:
ClassDef(AliCorrectionMatrix3D,1)
#include <TParticle.h>
#include <TParticlePDG.h>
+#include <TH1.h>
+#include <TH3F.h>
#include <AliLog.h>
#include <AliESD.h>
return kTRUE;
}
+
+//____________________________________________________________________
+void AliPWG0Helper::CreateProjections(TH3F* hist)
+{
+ // create projections of 3d hists to all 2d combinations
+ // the histograms are not returned, just use them from memory or use this to create them in a file
+
+ TH1* proj = hist->Project3D("yx");
+ proj->SetXTitle(hist->GetXaxis()->GetTitle());
+ proj->SetYTitle(hist->GetYaxis()->GetTitle());
+
+ proj = hist->Project3D("zx");
+ proj->SetXTitle(hist->GetXaxis()->GetTitle());
+ proj->SetYTitle(hist->GetZaxis()->GetTitle());
+
+ proj = hist->Project3D("zy");
+ proj->SetXTitle(hist->GetYaxis()->GetTitle());
+ proj->SetYTitle(hist->GetZaxis()->GetTitle());
+}
class AliESD;
class TParticle;
+class TH3F;
class AliPWG0Helper : public TObject
{
static Bool_t IsVertexReconstructed(AliESD* aEsd);
static Bool_t IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries);
+ static void CreateProjections(TH3F* hist);
+
protected:
ClassDef(AliPWG0Helper, 0)
};
// Constructor. Initialization of pointers
//
- AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
+ //AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
}
AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
return kFALSE;
}
+ if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
+ return kTRUE;
+
if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
return kTRUE;
Int_t nGoodTracks = list->GetEntries();
Float_t vertexRecoCorr = fdNdEtaCorrection->GetVertexRecoCorrection(vtx[2], nGoodTracks);
+ if (vertexRecoCorr <= 0)
+ {
+ AliDebug(AliLog::kError, Form("INFO: Skipping event because vertexRecoCorr is <= 0 (%f)", vertexRecoCorr));
+ delete list;
+ return kTRUE;
+ }
// loop over esd tracks
for (Int_t t=0; t<nGoodTracks; t++)
Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
Float_t pt = vector.Pt();
- // TODO pt cut
-
Float_t track2particleCorr = fdNdEtaCorrection->GetTrack2ParticleCorrection(vtx[2], eta, pt);
- fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, vertexRecoCorr * track2particleCorr);
+ Float_t weight = vertexRecoCorr * track2particleCorr;
+ if (weight <= 0)
+ {
+ AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f)", t, weight));
+ continue;
+ }
+ fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, weight);
} // end of track loop
delete list;
#include "AlidNdEtaCorrection.h"
#include <TCanvas.h>
-#include <TH2F.h>
+#include <TH3F.h>
+#include <TH1D.h>
//____________________________________________________________________
ClassImp(AlidNdEtaCorrection)
fTriggerBiasCorrection ->DrawHistograms();
}
+
+//____________________________________________________________________
+Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
+{
+ // calculates the fraction of particles measured (some are missed due to the pt cut off)
+ // uses the generated particle histogram from fTrack2ParticleCorrection
+
+ TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
+
+ // find eta borders, if eta is negative assume -0.8 ... 0.8
+ Int_t etaBegin = 0;
+ Int_t etaEnd = 0;
+ if (eta < 0)
+ {
+ etaBegin = generated->GetYaxis()->FindBin(-0.8);
+ etaEnd = generated->GetYaxis()->FindBin(0.8);
+ }
+ else
+ {
+ etaBegin = generated->GetYaxis()->FindBin(eta);
+ etaEnd = etaBegin;
+ }
+
+ Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
+ Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
+
+ TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
+
+ Int_t ptBin = ptProj->FindBin(ptCutOff);
+ Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
+ Float_t all = ptProj->Integral();
+
+ if (all == 0)
+ return -1;
+
+ Float_t fraction = abovePtCut / all;
+
+ if (debug)
+ {
+ new TCanvas;
+ ptProj->Draw();
+ }
+
+ return fraction;
+}
+
Float_t GetTrack2ParticleCorrection(Float_t vtx, Float_t eta, Float_t pt)
{return fTrack2ParticleCorrection->GetCorrection(vtx, eta, pt);}
-
+
Float_t GetVertexRecoCorrection(Float_t vtx, Float_t n) {return fVertexRecoCorrection->GetCorrection(vtx, n);}
Float_t GetTriggerBiasCorrection(Float_t eta, Float_t pt=0) {return fTriggerBiasCorrection->GetCorrection(eta, pt);}
+ Float_t GetMeasuredFraction(Float_t ptCutOff, Float_t eta = -1, Bool_t debug = kFALSE);
+
protected:
AliCorrectionMatrix3D* fTrack2ParticleCorrection; // handles the track-to-particle correction, function of vtx_z, eta, pt
AliCorrectionMatrix2D* fVertexRecoCorrection; // handles the vertex reconstruction efficiency, function of n_clustersITS and vtx_z
Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
- // check if the event was triggered
Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
fdNdEtaCorrection->IncreaseEventCount();
fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
} // end of track loop
- fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
- if (vertexReconstructed)
- fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
+
+ if (eventTriggered)
+ {
+ fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
+ if (vertexReconstructed)
+ fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
+ }
return kTRUE;
}
{
gSystem->Load("libPWG0base");
- dNdEtaCorrection* dNdEtaCorrection = 0;
+ AlidNdEtaCorrection* dNdEtaCorrection = 0;
if (correct)
{
- dNdEtaCorrection = new dNdEtaCorrection();
+ dNdEtaCorrection = new AlidNdEtaCorrection();
dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
- dNdEtaCorrection->RemoveEdges(2, 0, 2);
+ //dNdEtaCorrection->RemoveEdges(2, 0, 2);
}
fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
}
fdNdEtaAnalysis->LoadHistograms();
- fdNdEtaAnalysis->Finish(dNdEtaCorrection);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3);
fdNdEtaAnalysis->DrawHistograms();
}
#include <TFile.h>
#include <TH3F.h>
+#include <TH2D.h>
#include <TH1D.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TLegend.h>
#include "AlidNdEtaCorrection.h"
+#include "AliPWG0Helper.h"
//____________________________________________________________________
ClassImp(dNdEtaAnalysis)
{
// constructor
- fData = new TH3F(Form("%s_analysis", name),"",80,-20,20,120,-6,6,100, 0, 10);
+ fData = new TH3F(Form("%s_analysis", name),"dNdEtaAnalysis",80,-20,20,120,-6,6,100, 0, 10);
fData->SetXTitle("vtx z [cm]");
fData->SetYTitle("#eta");
fData->SetZTitle("p_{T}");
{
// fills an event into the histograms
- fVtx->Fill(vtx, weight); // TODO vtx distribution with or without weight?
+ fVtx->Fill(vtx, weight);
fNEvents += weight;
}
//____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
{
// correct with correction values if available
if (!correction)
printf("INFO: No correction applied\n");
- // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
- /* for (Int_t iVtx=0; iVtx<=fDataUncorrected->GetNbinsX(); iVtx++)
- {
- for (Int_t iEta=0; iEta<=fDataUncorrected->GetNbinsY(); iEta++)
- {
- Float_t correctionValue = 1;
- if (correction)
- correctionValue = correction->GetTrack2ParticleCorrection(fDataUncorrected->GetXaxis()->GetBinCenter(iVtx), fDataUncorrected->GetYaxis()->GetBinCenter(iEta), 1.0);
-
- Float_t value = fDataUncorrected->GetBinContent(iVtx, iEta);
- Float_t error = fDataUncorrected->GetBinError(iVtx, iEta);
-
- Float_t correctedValue = value * correctionValue;
- Float_t correctedError = error * correctionValue;
+ // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
- if (correctedValue != 0)
- {
- fData->SetBinContent(iVtx, iEta, correctedValue);
- fData->SetBinError(iVtx, iEta, correctedError);
- }
- }
+ // integrate over pt (with pt cut)
+ fData->GetZaxis()->SetRange(fData->GetZaxis()->FindBin(ptCut), fData->GetZaxis()->GetNbins());
+ TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2"));
+ if (vtxVsEta == 0)
+ {
+ printf("ERROR: pt integration failed\n");
+ return;
}
- for (Int_t iEta=0; iEta<=fData->GetNbinsY(); iEta++)
+ new TCanvas;
+ vtxVsEta->Draw("COLZ");
+
+ for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
{
// do we have several histograms for different vertex positions?
Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
Float_t sumError2 = 0;
for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
{
- if (fData->GetBinContent(iVtx, iEta) != 0)
+ if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
{
- sum = sum + fData->GetBinContent(iVtx, iEta);
- sumError2 = sumError2 + TMath::Power(fData->GetBinError(iVtx, iEta),2);
+ sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
+ sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
}
}
- Float_t dndeta = sum / totalEvents;
- Float_t error = TMath::Sqrt(sumError2) / totalEvents;
+ Float_t ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
+ //ptCutOffCorrection = 1;
+ if (ptCutOffCorrection <= 0)
+ {
+ printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
+ continue;
+ }
+
+ Float_t dndeta = sum / totalEvents / ptCutOffCorrection;
+ Float_t error = TMath::Sqrt(sumError2) / totalEvents / ptCutOffCorrection;
dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
error = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
fdNdEta[vertexPos]->SetBinError(iEta, error);
}
- }*/
+ }
}
//____________________________________________________________________
gDirectory->cd(GetName());
fData ->Write();
+ AliPWG0Helper::CreateProjections(fData);
fDataUncorrected->Write();
+ AliPWG0Helper::CreateProjections(fDataUncorrected);
+
fVtx ->Write();
for (Int_t i=0; i<kVertexBinning; ++i)
fdNdEta[i] ->Write();
void FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight);
void FillEvent(Float_t vtx, Float_t weight);
- void Finish(AlidNdEtaCorrection* correction);
+ void Finish(AlidNdEtaCorrection* correction, Float_t ptCut);
void DrawHistograms();
void LoadHistograms();
}
}
- chain->SetProof(proof);
+ //chain->SetProof(proof);
}
// ########################################################
if (proof)
proof->AddInput(esdTrackCuts);
+ if (aMC == kFALSE)
+ {
+ AlidNdEtaCorrection* dNdEtaCorrection = new AlidNdEtaCorrection();
+ dNdEtaCorrection->LoadHistograms("correction_map.root","dndeta_correction");
+ //dNdEtaCorrection->RemoveEdges(2, 0, 2);
+
+ chain->GetUserInfo()->Add(dNdEtaCorrection);
+ if (proof)
+ proof->AddInput(dNdEtaCorrection);
+ }
+
TString selectorName = ((aMC == kFALSE) ? "AlidNdEtaAnalysisESDSelector" : "AlidNdEtaAnalysisMCSelector");
AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
TStopwatch timer;
timer.Start();
- Long64_t result = chain->Process(selectorName);
+ Long64_t result = -1;
+
+ if (proof != kFALSE)
+ result = chain->MakeTDSet()->Process(selectorName);
+ else
+ result = chain->Process(selectorName);
+
if (result != 0)
{
printf("ERROR: Executing process failed with %d.\n", result);