TH2* measuredEvents = GetEventCorrection()->GetMeasuredHistogram();
TH2* generatedEvents = GetEventCorrection()->GetGeneratedHistogram();
+ Printf("AliCorrection::PrintInfo: Whole phasespace:");
+
Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %1.f, events generated %1.f", measured->Integral(), generated->Integral(), measuredEvents->Integral(), generatedEvents->Integral());
- // normalize to number of events;
- measured->Scale(1.0 / measuredEvents->Integral());
- generated->Scale(1.0 / generatedEvents->Integral());
+ Printf("AliCorrection::PrintInfo: Values in |eta| < 0.8, |vtx-z| < 10 and pt > %.2f:", ptCut);
+
+ Float_t tracksM = measured->Integral(measured->GetXaxis()->FindBin(-9.9), measured->GetXaxis()->FindBin(9.9), measured->GetYaxis()->FindBin(-0.79), measured->GetYaxis()->FindBin(0.79), measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
+ Float_t tracksG = generated->Integral(generated->GetXaxis()->FindBin(-9.9), generated->GetXaxis()->FindBin(9.9), generated->GetYaxis()->FindBin(-0.79), generated->GetYaxis()->FindBin(0.79), generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
+
+ Float_t eventsM = measuredEvents->Integral(measuredEvents->GetXaxis()->FindBin(-9.9), measuredEvents->GetXaxis()->FindBin(9.9), 1, measuredEvents->GetNbinsY());
+ Float_t eventsG = generatedEvents->Integral(generatedEvents->GetXaxis()->FindBin(-9.9), generatedEvents->GetXaxis()->FindBin(9.9), 1, generatedEvents->GetNbinsY());
+
+ Printf("tracks measured: %.1f tracks generated: %.1f, events measured: %1.f, events generated %1.f", tracksM, tracksG, eventsM, eventsG);
- Float_t nMeasured = measured->Integral(-1, -1, -1, -1, measured->GetZaxis()->FindBin(ptCut), measured->GetZaxis()->GetNbins());
- Float_t nGenerated = generated->Integral(-1, -1, -1, -1, generated->GetZaxis()->FindBin(ptCut), generated->GetZaxis()->GetNbins());
+ if (tracksM > 0)
+ Printf("Effective average correction factor for TRACKS: %.3f", tracksG / tracksM);
+ if (eventsM > 0)
+ Printf("Effective average correction factor for EVENTS: %.3f", eventsG / eventsM);
- Printf("%.2f tracks/event measured, %.2f tracks after correction --> effective average correction factor is %.2f (pt cut %.2f GeV/c)", nMeasured, nGenerated, nGenerated / nMeasured, ptCut);
+ if (eventsM > 0 && eventsG > 0)
+ {
+ // normalize to number of events;
+ tracksM /= eventsM;
+ tracksG /= eventsG;
+
+ Printf("%.2f tracks/event measured, %.2f tracks/event after correction --> effective average correction factor is %.3f (pt cut %.2f GeV/c)", tracksM, tracksG, tracksG / tracksM, ptCut);
+ }
}
for (Int_t x=1; x<=fhCorr->GetNbinsX(); ++x)
for (Int_t y=1; y<=fhCorr->GetNbinsY(); ++y)
for (Int_t z=1; z<=fhCorr->GetNbinsZ(); ++z)
+ {
fhCorr->SetBinContent(x, y, z, 1);
+ fhCorr->SetBinError(x, y, z, 0);
+ }
}
//____________________________________________________________________
etaEnd = etaBegin;
}
+ //Printf("AlidNdEtaCorrection::GetMeasuredFraction: Using vtx range of +- 10 cm");
Int_t vertexBegin = generated->GetXaxis()->FindBin(-9.99);
Int_t vertexEnd = generated->GetXaxis()->FindBin(9.99);
return fraction;
}
+//____________________________________________________________________
+TH1* AlidNdEtaCorrection::GetMeasuredEventFraction(CorrectionType correctionType, Int_t multCut)
+{
+ // calculates the fraction of events above multCut (but including it)
+ //
+ // uses the generated event histogram from the correction passed, e.g. pass GetTrack2ParticleCorrection()
+
+ if (!GetCorrection(correctionType))
+ return 0;
+
+ const TH2F* generated = GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
+
+ TH1* allEvents = generated->ProjectionX(Form("%s_all", generated->GetName()), 1, generated->GetNbinsY());
+ TH1* aboveEvents = generated->ProjectionX(Form("%s_above", generated->GetName()), generated->GetYaxis()->FindBin(multCut), generated->GetNbinsY());
+
+ aboveEvents->Divide(aboveEvents, allEvents, 1, 1, "B");
+
+ return aboveEvents;
+}
+
//____________________________________________________________________
void AlidNdEtaCorrection::ReduceInformation()
{
#include "AliPWG0Helper.h"
class AliCorrection;
+class TH1;
class AlidNdEtaCorrection : public TNamed
{
void DrawOverview(const char* canvasName = 0);
Float_t GetMeasuredFraction(CorrectionType correctionType, Float_t ptCutOff, Float_t eta = -100, Bool_t debug = kFALSE);
+ TH1* GetMeasuredEventFraction(CorrectionType correctionType, Int_t multCut);
void ReduceInformation();
for (Int_t i=0; i<kVertexBinning; ++i)
{
+ fdNdEtaNotEventCorrected[i] = 0;
fdNdEta[i] = 0;
fdNdEtaPtCutOffCorrected[i] = 0;
}
fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 40, -2, 2);
+ fdNdEtaNotEventCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_noteventcorrected", fdNdEta[0]->GetName())));
fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
for (Int_t i=1; i<kVertexBinning; ++i)
{
fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
+ fdNdEtaNotEventCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaNotEventCorrected[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
}
for (Int_t i=0; i<kVertexBinning; ++i)
{
+ if (fdNdEtaNotEventCorrected[i])
+ {
+ delete fdNdEtaNotEventCorrected[i];
+ fdNdEtaNotEventCorrected[i] = 0;
+ }
if (fdNdEta[i])
{
delete fdNdEta[i];
for (Int_t i=0; i<kVertexBinning; ++i)
{
+ target.fdNdEtaNotEventCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaNotEventCorrected[i]->Clone());
target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
}
}
//____________________________________________________________________
-void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType)
+void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, Int_t multCut)
{
//
// correct with the given correction values and calculate dNdEta and pT distribution
// the corrections that are applied can be steered by the flag correctionType
+ // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
//
// TODO put tag somewhere which corrections have been applied
+ if (multCut > 1)
+ {
+ Printf("ERROR: A bigger multiplicity cut than 1 is not possible in the current implementation");
+ return;
+ }
+
// set corrections to 1
fData->SetCorrectionToUnity();
printf("INFO: No correction applied\n");
fData->Multiply();
+ fData->PrintInfo(ptCut);
TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
// integrate multiplicity axis out (include under/overflow bins!!!)
TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
- TH1D* vertexHist = tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
+
+ // multiplicity cut correction
+ // correct for not-measured events with too small multiplicity
+ // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
+ TH1D* vertexHistUncorrected = tmp->ProjectionX("_px", tmp->GetYaxis()->FindBin(multCut), tmp->GetNbinsY() + 1, "e");
+ TH1D* vertexHist = (TH1D*) vertexHistUncorrected->Clone("vertexHist");
+
+ Printf("Using %d events (above a cut of %d)", (Int_t) vertexHistUncorrected->Integral(), multCut);
+ if (correction && correctionType >= AlidNdEtaCorrection::kVertexReco)
+ {
+ TH1* eventFraction = correction->GetMeasuredEventFraction(correctionType, multCut);
+ vertexHist->Divide(eventFraction);
+ Printf("Multiplicity cut correction: Corrected from %d events to %d events", (Int_t) vertexHistUncorrected->Integral(), (Int_t) vertexHist->Integral());
+ }
// create pt hist
{
//printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
+ Float_t totalEventsUncorrected = vertexHistUncorrected->Integral(vertexBinBegin, vertexBinEnd - 1);
+ if (totalEventsUncorrected == 0)
+ {
+ printf("WARNING: No events (uncorrected) for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
+ continue;
+ }
+
Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
if (totalEvents == 0)
{
Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
if (bin > 0 && bin < fdNdEta[vertexPos]->GetNbinsX())
{
- Float_t dndeta = sum / totalEvents;
- Float_t error = TMath::Sqrt(sumError2) / totalEvents;
+ Float_t dndeta = sum / totalEventsUncorrected;
+ Float_t error = TMath::Sqrt(sumError2) / totalEventsUncorrected;
+
+ dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
+ error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
+
+ fdNdEtaNotEventCorrected[vertexPos]->SetBinContent(bin, dndeta);
+ fdNdEtaNotEventCorrected[vertexPos]->SetBinError(bin, error);
+
+ dndeta = sum / totalEvents;
+ error = TMath::Sqrt(sumError2) / totalEvents;
dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
for (Int_t i=0; i<kVertexBinning; ++i)
{
+ if (fdNdEtaNotEventCorrected[i])
+ fdNdEtaNotEventCorrected[i]->Write();
+ else
+ printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaNotEventCorrected[%d] is 0\n", i);
+
if (fdNdEta[i])
fdNdEta[i]->Write();
else
for (Int_t i=0; i<kVertexBinning; ++i)
{
+ if (dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaNotEventCorrected[i]->GetName())))
+ fdNdEtaNotEventCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaNotEventCorrected[i]->GetName()));
+
fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
}
TObject* obj;
// sub collections
- const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
+ const Int_t nCollections = 3 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
TList* collections[nCollections];
for (Int_t i=0; i<nCollections; ++i)
collections[i] = new TList;
{
collections[2+i]->Add(entry->fdNdEta[i]);
collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
+ collections[2+2*kVertexBinning+i]->Add(entry->fdNdEtaNotEventCorrected[i]);
}
++count;
for (Int_t i=0; i<kVertexBinning; ++i)
{
fdNdEta[i]->Merge(collections[2+i]);
- fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
+ fdNdEtaPtCutOffCorrected[i]->Merge(collections[2+kVertexBinning+i]);
+ fdNdEtaNotEventCorrected[i]->Merge(collections[2+kVertexBinning+i]);
}
for (Int_t i=0; i<nCollections; ++i)
void FillTrack(Float_t vtx, Float_t eta, Float_t pt);
void FillEvent(Float_t vtx, Float_t n);
- void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType);
+ void Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, Int_t multCut = 0);
void DrawHistograms(Bool_t simple = kFALSE);
void LoadHistograms(const Char_t* dir = 0);
TH1F* fPtDist; // pt distribution
- TH1F* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
- TH1F* fdNdEtaPtCutOffCorrected[kVertexBinning]; // dndeta results for different vertex bins (0 = full range), pt cut off corrected
+ TH1F* fdNdEtaNotEventCorrected[kVertexBinning]; // dndeta results for different vertex bins (0 = full range)
+ TH1F* fdNdEta[kVertexBinning]; // dndeta results for different vertex bins (0 = full range), mult cut off corrected
+ TH1F* fdNdEtaPtCutOffCorrected[kVertexBinning]; // dndeta results for different vertex bins (0 = full range), mult + pt cut off corrected
ClassDef(dNdEtaAnalysis, 1)
};
// --- end of helpers --- begin functions ---
-void DrawOverview(const char* fileName, const char* dirName)
+void DrawOverview(const char* fileName = "correction_map.root", const char* dirName = "dndeta_correction")
{
loadlibs();
if (!TFile::Open(fileName))
return diffFullRange;
}
-void dNdEta(Bool_t onlyESD = kFALSE)
+void dNdEta(Bool_t onlyESD = kFALSE, Bool_t save = kTRUE)
{
TFile* file = TFile::Open("analysis_esd.root");
TH1* histESD = (TH1*) file->Get("dndeta/dNdEta_corrected");
TH1* histESDMBVtxNoPt = (TH1*) file->Get("dndetaTrVtx/dNdEta");
TH1* histESDMBTracksNoPt = (TH1*) file->Get("dndetaTracks/dNdEta");
+ TH1* histESDNoEvCorr = (TH1*) file->Get("dndeta/dNdEta_noteventcorrected");
+
Prepare1DPlot(histESD);
Prepare1DPlot(histESDMB);
Prepare1DPlot(histESDMBVtx);
Prepare1DPlot(histESDMBVtxNoPt);
Prepare1DPlot(histESDMBTracksNoPt);
+ Prepare1DPlot(histESDNoEvCorr);
+
histESD->SetLineWidth(0);
histESDMB->SetLineWidth(0);
histESDMBVtx->SetLineWidth(0);
histESDMB->SetMarkerColor(2);
histESDMBVtx->SetMarkerColor(3);
+ histESD->SetLineColor(1);
+ histESDMB->SetLineColor(2);
+ histESDMBVtx->SetLineColor(3);
+
histESDNoPt->SetMarkerColor(1);
histESDMBNoPt->SetMarkerColor(2);
histESDMBVtxNoPt->SetMarkerColor(3);
histESDMBTracksNoPt->SetMarkerColor(4);
+ histESDNoEvCorr->SetMarkerColor(6);
+
histESD->SetMarkerStyle(20);
histESDMB->SetMarkerStyle(21);
histESDMBVtx->SetMarkerStyle(22);
histESDMBNoPt->SetMarkerStyle(21);
histESDMBVtxNoPt->SetMarkerStyle(22);
histESDMBTracksNoPt->SetMarkerStyle(23);
+
+ histESDNoEvCorr->SetMarkerStyle(29);
TH2F* dummy = new TH2F("dummy", "", 100, -1.5, 1.5, 1000, 0, histESDMBVtx->GetMaximum() * 1.1);
Prepare1DPlot(dummy);
dummy->SetYTitle("dN_{ch}/d#eta");
dummy->GetYaxis()->SetTitleOffset(1);
- Float_t etaLimit = 0.7999;
+ Float_t etaLimit = 1.2999;
histESDMBVtx->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESDMB->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
histESDMB->Draw("SAME");
histESD->Draw("SAME");
- canvas->SaveAs("dNdEta1.gif");
- canvas->SaveAs("dNdEta1.eps");
+ if (save)
+ {
+ canvas->SaveAs("dNdEta1.gif");
+ canvas->SaveAs("dNdEta1.eps");
+ }
if (onlyESD)
return;
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
fdNdEtaAnalysis->LoadHistograms();
+ dNdEtaAnalysis* fdNdEtaAnalysis2 = (dNdEtaAnalysis*) fdNdEtaAnalysis->Clone();
+
fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
- TH1* histMCPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
+ TH1* histMCPtCut = (TH1*) fdNdEtaAnalysis->GetdNdEtaHistogram(0)->Clone("histMCPtCut");
+
+ fdNdEtaAnalysis2->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
+ TH1* histMCPtCutNoEvCorr = (TH1*) fdNdEtaAnalysis2->GetdNdEtaHistogram(0)->Clone("histMCPtCutNoEvCorr");
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
fdNdEtaAnalysis->LoadHistograms();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
fdNdEtaAnalysis->LoadHistograms();
- fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
TH1* histMCTrVtxPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
fdNdEtaAnalysis->LoadHistograms();
- fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
TH1* histMCTracksPtCut = fdNdEtaAnalysis->GetdNdEtaHistogram(0);
TCanvas* canvas2 = new TCanvas("dNdEta2", "dNdEta2", 500, 500);
Prepare1DPlot(histMCTrVtx);
Prepare1DPlot(histMCPtCut);
+ Prepare1DPlot(histMCPtCutNoEvCorr);
Prepare1DPlot(histMCTrPtCut);
Prepare1DPlot(histMCTrVtxPtCut);
if (histMCTracksPtCut)
histMCTrVtxPtCut->SetLineColor(3);
if (histMCTracksPtCut)
histMCTracksPtCut->SetLineColor(4);
+ histMCPtCutNoEvCorr->SetLineColor(6);
TH2* dummy2 = (TH2F*) dummy->Clone("dummy2");
Prepare1DPlot(dummy2);
histESDMBNoPt->Draw("SAME");
histESDMBVtxNoPt->Draw("SAME");
histESDMBTracksNoPt->Draw("SAME");
+ histESDNoEvCorr->Draw("SAME");
histMCPtCut->Draw("SAME");
+ histMCPtCutNoEvCorr->Draw("SAME");
histMCTrPtCut->Draw("SAME");
histMCTrVtxPtCut->Draw("SAME");
if (histMCTracksPtCut)
histMCTracksPtCut->Draw("SAME");
- canvas2->SaveAs("dNdEta2.gif");
- canvas2->SaveAs("dNdEta2.eps");
+ if (save)
+ {
+ canvas2->SaveAs("dNdEta2.gif");
+ canvas2->SaveAs("dNdEta2.eps");
+ }
TH1* ratio = (TH1*) histMC->Clone("ratio");
TH1* ratioNoPt = (TH1*) histMCPtCut->Clone("ratioNoPt");
ratio->Divide(histESD);
ratioNoPt->Divide(histESDNoPt);
- ratio->GetXaxis()->SetRangeUser(-0.7999, 0.7999);
+ ratio->GetXaxis()->SetRangeUser(-etaLimit, etaLimit);
ratio->SetLineColor(1);
ratioNoPt->SetLineColor(2);
canvas3->Modified();
- canvas3->SaveAs("dNdEta.gif");
- canvas3->SaveAs("dNdEta.eps");
+ if (save)
+ {
+ canvas3->SaveAs("dNdEta.gif");
+ canvas3->SaveAs("dNdEta.eps");
+ }
TCanvas* canvas4 = new TCanvas("ratio", "ratio", 700, 500);
TH3* hist1 = (TH3*) dNdEtaCorrection->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("mc");
hist1->SetTitle("mc");
+ Printf("mc contains %f entries", hist1->Integral());
+ Printf("mc contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
TH3* hist2 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd");
hist2->SetTitle("esd");
+ Printf("esd contains %f entries", hist2->Integral());
+ Printf("esd contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
AliPWG0Helper::CreateDividedProjections(hist1, hist2);
new TCanvas; gROOT->FindObject("mc_yx_div_esd_yx")->Draw("COLZ");
new TCanvas; gROOT->FindObject("mc_zx_div_esd_zx")->Draw("COLZ");
new TCanvas; gROOT->FindObject("mc_zy_div_esd_zy")->Draw("COLZ");
+}
+
+void CompareMeasured2Measured(const char* dataInput = "analysis_esd_raw.root", const char* dataInput2 = "analysis_esd_raw.root")
+{
+ loadlibs();
+
+ TFile* file = TFile::Open(dataInput);
+
+ if (!file)
+ {
+ cout << "Error. File not found" << endl;
+ return;
+ }
+
+ dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
+ fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
+
+ TFile* file = TFile::Open(dataInput2);
+
+ if (!file)
+ {
+ cout << "Error. File not found" << endl;
+ return;
+ }
+
+ dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis("dndeta2", "dndeta2");
+ fdNdEtaAnalysis2->LoadHistograms("fdNdEtaAnalysisESD");
+
+ gROOT->cd();
+
+ TH3* hist1 = (TH3*) fdNdEtaAnalysis->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd1");
+ hist1->SetTitle("esd1");
+ Printf("esd1 contains %f entries", hist1->GetEntries());
+ Printf("esd1 contains %f entries in |vtx-z| < 10, pt > 0.3", hist1->Integral(hist1->GetXaxis()->FindBin(-9.9), hist1->GetXaxis()->FindBin(9.9), 1, hist1->GetNbinsY(), hist1->GetZaxis()->FindBin(0.301), hist1->GetNbinsZ()));
+
+ TH3* hist2 = (TH3*) fdNdEtaAnalysis2->GetData()->GetTrackCorrection()->GetMeasuredHistogram()->Clone("esd2");
+ hist2->SetTitle("esd2");
+ Printf("esd2 contains %f entries", hist2->GetEntries());
+ Printf("esd2 contains %f entries in |vtx-z| < 10, pt > 0.3", hist2->Integral(hist2->GetXaxis()->FindBin(-9.9), hist2->GetXaxis()->FindBin(9.9), 1, hist2->GetNbinsY(), hist2->GetZaxis()->FindBin(0.301), hist2->GetNbinsZ()));
+
+ AliPWG0Helper::CreateDividedProjections(hist1, hist2);
+
+ new TCanvas; gROOT->FindObject("esd1_yx_div_esd2_yx")->Draw("COLZ");
+ new TCanvas; gROOT->FindObject("esd1_zx_div_esd2_zx")->Draw("COLZ");
+ new TCanvas; gROOT->FindObject("esd1_zy_div_esd2_zy")->Draw("COLZ");
+
+ TH2* event1 = (TH2*) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event1");
+ TH2* event2 = (TH2*) fdNdEtaAnalysis2->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Clone("event2");
+
+ Printf("event1 contains %f entries", event1->GetEntries());
+ Printf("event2 contains %f entries", event2->GetEntries());
+ Printf("event1 integral is %f", event1->Integral());
+ Printf("event2 integral is %f", event2->Integral());
+ Printf("event1 contains %f entries in |vtx-z| < 10", event1->Integral(event1->GetXaxis()->FindBin(-9.9), event1->GetXaxis()->FindBin(9.9), 1, event1->GetNbinsY()));
+ Printf("event2 contains %f entries in |vtx-z| < 10", event2->Integral(event2->GetXaxis()->FindBin(-9.9), event2->GetXaxis()->FindBin(9.9), 1, event2->GetNbinsY()));
+
+ projx1 = event1->ProjectionX();
+ projx2 = event2->ProjectionX();
+
+ new TCanvas; projx1->DrawCopy(); projx2->SetLineColor(2); projx2->DrawCopy("SAME");
+
+ projx1->Divide(projx2);
+ new TCanvas; projx1->Draw();
+
+ event1->Divide(event2);
+ new TCanvas; event1->Draw("COLZ");
}
-
+
Track2Particle1DComposition(fileNames, nCompositions, folderNames);
}
-Double_t ConvSigma1To2D(Double_t sigma)
-{
- return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
-}
-
-Double_t ConvDistance1DTo2D(Double_t distance)
-{
- return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
-}
-
-Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
-{
- Double_t count = 0;
-
- //nSigma = ConvSigma1To2D(nSigma);
-
- for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
- for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
- {
- Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
- Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
-
- Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
-
- d = ConvDistance1DTo2D(d);
-
- if (d < nSigma)
- count += tracks->GetBinContent(x, y);
- }
-
- return count;
-}
-
-TH2F* Sigma2VertexGaussianTracksHist()
-{
- TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
-
- TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
- gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
-
- for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
- for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
- tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
-
- //normalize
- tracks->Scale(1.0 / tracks->Integral());
-
- return tracks;
-}
-
-TH1F* Sigma2VertexGaussian()
-{
- TH2F* tracks = Sigma2VertexGaussianTracksHist();
-
- TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
- canvas->Divide(2, 2);
-
- canvas->cd(1);
- tracks->Draw("COLZ");
-
- TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
- for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
- ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
- ratio->SetMarkerStyle(21);
-
- canvas->cd(2);
- ratio->DrawCopy("P");
-
- TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
- Double_t sigma3 = Sigma2VertexCount(tracks, 4);
- for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
- ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
- ratio2->SetMarkerStyle(21);
-
- canvas->cd(3);
- ratio2->DrawCopy("P");
-
- canvas->SaveAs("Sigma2Vertex.eps");
-
- return ratio2;
-}
-
-TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
-{
- TFile* file = TFile::Open(fileName);
-
- TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
- TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
- if (!sigmavertex || !sigmavertexPrim)
- {
- printf("Could not read histogram(s)\n");
- return;
- }
-
- // calculate ratio
- TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
-
- // calculate contamination
- TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
- contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
-
- for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
- {
- ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
- contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
- }
-
- // print stats
- for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
- {
- Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
- Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
- Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
- Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
- }
-
- TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
- canvas->Divide(2, 1);
-
- canvas->cd(1);
- sigmavertex->SetMarkerStyle(21);
- sigmavertex->Draw("P");
-
- canvas->cd(2);
- ratio->SetMarkerStyle(21);
- ratio->DrawCopy("P");
-
- contamination->DrawCopy("SAME");
-
- TH1F** returnContainer = new TH1F*[2];
- returnContainer[0] = ratio;
- returnContainer[1] = contamination;
-
- return returnContainer;
-}
-
-void Sigma2VertexCompare(const char* fileName = "systematics.root")
-{
- TH1F* ratio1 = Sigma2VertexGaussian();
-
- TH1F** hists = Sigma2VertexSimulation(fileName);
- TH1F* ratio2 = hists[0];
- TH1F* contamination = hists[1];
-
- ratio1->SetStats(kFALSE);
- ratio2->SetStats(kFALSE);
-
- ratio1->SetMarkerStyle(0);
- ratio2->SetMarkerStyle(0);
-
- ratio1->SetLineWidth(2);
- ratio2->SetLineWidth(2);
-
- TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
- legend->SetFillColor(0);
- legend->AddEntry(ratio1, "Gaussian");
- legend->AddEntry(ratio2, "Simulation");
- legend->AddEntry(contamination, "1 + Contamination");
-
- ratio2->SetTitle("");
- ratio2->GetYaxis()->SetTitleOffset(1.5);
- ratio2->GetXaxis()->SetRangeUser(2, 5);
-
- TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
- InitPad();
-
- ratio2->SetMarkerStyle(21);
- ratio1->SetMarkerStyle(22);
-
- ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
- ratio2->SetLineColor(kRed);
- ratio2->SetMarkerColor(kRed);
- ratio2->Draw("PL");
- ratio1->Draw("SAMEPL");
-
- contamination->Draw("SAME");
-
- legend->Draw();
-
- canvas->SaveAs("Sigma2VertexCompare.eps");
-}
void drawSystematics()
{
gSystem->Load("libPWG0base");
}
+void InitPad()
+{
+ if (!gPad)
+ return;
+
+ gPad->Range(0, 0, 1, 1);
+ gPad->SetLeftMargin(0.15);
+ //gPad->SetRightMargin(0.05);
+ //gPad->SetTopMargin(0.13);
+ //gPad->SetBottomMargin(0.1);
+
+ gPad->SetGridx();
+ gPad->SetGridy();
+}
+
void DrawpiKpAndCombinedZOnly(Float_t upperPtLimit=0.99)
{
//gROOT->ProcessLine(".L drawPlots.C");
const char* legendNames[] = { "#pi", "K", "p", "standard" };
Int_t folderCount = 3;
-
-
TH2F* h2DCorrections[4];
TH1F* h1DCorrections[4];
for (Int_t i=0; i<4; i++) {
void drawdNdEtaRatios(const char* canvasName, Int_t n, const char** files, const char** dirs, const char** names, Int_t* histID)
{
- gSystem->Load("libPWG0base");
+ loadlibs();
+ gROOT->ProcessLine(".L $ALICE_ROOT/PWG0/dNdEta/drawPlots.C");
TCanvas* canvas = new TCanvas(canvasName, canvasName, 1000, 500);
canvas->Divide(2, 1);
+ canvas->cd(2)->SetGridx();
+ canvas->cd(2)->SetGridy();
TLegend* legend = new TLegend(0.63, 0.73, 0.98, 0.98);
legend->SetFillColor(0);
if (i != 0)
{
+ PrintIntegratedDeviation(hist, base, names[i]);
+
canvas->cd(2);
- hist->Divide(hist, base, 1, 1, "B");
- hist->GetYaxis()->SetRangeUser(0.98, 1.02);
+ hist->Divide(hist, base, 1, 1);
+ hist->GetYaxis()->SetRangeUser(0.9, 1.1);
hist->DrawCopy((i == 1) ? "" : "SAME");
}
}
canvas->SaveAs(Form("%s.eps", canvas->GetName()));
}
+void MaterialBudgetChange()
+{
+ const char* files[] =
+ { "Material-normal-mcvtx/analysis_esd.root",
+ "Material-increased-mcvtx/analysis_esd.root",
+ "Material-decreased-mcvtx/analysis_esd.root" };
+
+ const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+ const char* names[] = { "no change", "+ 10 % material", "- 10 % material" };
+ Int_t hist[] = { 0, 0, 0 };
+
+ drawdNdEtaRatios("MaterialBudgetChange", 3, files, dirs, names, hist);
+}
+
+void MisalignmentChange()
+{
+ const char* files[] =
+ { "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-sim-mcvtx/analysis_esd.root",
+ "maps/v4-09-Release/tpc-only/fullC-simrec--fullA-simrec-mcvtx/analysis_esd.root" };
+
+ const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+ const char* names[] = { "no change", "increased in sim/rec", "increased in sim" };
+ Int_t hist[] = { 0, 0, 0 };
+
+ drawdNdEtaRatios("MisalignmentChange", 2, files, dirs, names, hist);
+}
+
+void dNdEtaVertexRanges()
+{
+ const char* files[] =
+ { "analysis_esd.root",
+ "analysis_esd.root",
+ "analysis_esd.root" };
+
+ const char* dirs[] = { "dndeta", "dndeta", "dndeta"};
+ const char* names[] = { "|vtx-z| < 10 cm", "-10 cm < vtx-z < 0 cm", "0 cm < vtx-z < 10 cm" };
+ Int_t hist[] = { 0, 1, 2 };
+
+ drawdNdEtaRatios("dNdEtaVertexRanges", 3, files, dirs, names, hist);
+}
+
void vertexShiftStudy(Int_t histID)
{
const char* files[] = { "maps/idealA/mc-vertex/analysis_esd.root", "results/idealC-idealA/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.05/analysis_esd.root", "maps/idealA/mc-vertex-shift-0.1/analysis_esd.root", "maps/idealA/mc-vertex-shift-dep/analysis_esd.root" };
canvas->SaveAs(Form("%s.eps", canvas->GetName()));
}
-void CompareRawTrackPlots(Float_t ptCut = 0.0)
+void CompareRawTrackPlots(const char* fileName1, const char* fileName2, Float_t ptCut = 0.0)
{
loadlibs();
const char* dirName = "fdNdEtaAnalysisESD";
- TFile* file = TFile::Open("fullA-simrec/MB2/analysis_esd_raw.root");
+ TFile* file = TFile::Open(fileName1);
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis(dirName, dirName);
fdNdEtaAnalysis->LoadHistograms();
- //TFile* file2 = TFile::Open("fullA-sim/analysis_esd_raw.root");
- TFile* file2 = TFile::Open("fullA-simrecexcepttpc/analysis_esd_raw.root");
+ TFile* file2 = TFile::Open(fileName2);
dNdEtaAnalysis* fdNdEtaAnalysis2 = new dNdEtaAnalysis(dirName, dirName);
fdNdEtaAnalysis2->LoadHistograms();
track1->Scale(1.0 / event1Count);
track2->Scale(1.0 / event2Count);
- Float_t nTrack1 = track1->Integral(-1, -1, -1, -1, track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
- Float_t nTrack2 = track2->Integral(-1, -1, -1, -1, track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
+ Float_t nTrack1 = track1->Integral(track1->GetXaxis()->FindBin(-9.9), track1->GetXaxis()->FindBin(9.9), track1->GetYaxis()->FindBin(-0.79), track1->GetYaxis()->FindBin(0.79), track1->GetZaxis()->FindBin(ptCut), track1->GetZaxis()->GetNbins());
+ Float_t nTrack2 = track2->Integral(track2->GetXaxis()->FindBin(-9.9), track2->GetXaxis()->FindBin(9.9), track2->GetYaxis()->FindBin(-0.79), track2->GetYaxis()->FindBin(0.79), track2->GetZaxis()->FindBin(ptCut), track2->GetZaxis()->GetNbins());
Printf("There are %.2f tracks/event in the first sample and %.2f tracks/event in the second sample. %.2f %% difference (with pt cut at %.2f GeV/c)", nTrack1, nTrack2, 100.0 * (nTrack1 - nTrack2) / nTrack1, ptCut);
new TCanvas; gROOT->FindObject("track1_zx_div_track2_zx")->Draw("COLZ");
new TCanvas; gROOT->FindObject("track1_zy_div_track2_zy")->Draw("COLZ");
- new TCanvas;
- proj1 = track1->Project3D("ze2");
- proj2 = track2->Project3D("ze2");
- AliPWG0Helper::NormalizeToBinWidth(proj1);
- AliPWG0Helper::NormalizeToBinWidth(proj2);
-
- proj1->DrawCopy();
- proj2->SetLineColor(2);
- proj2->SetMarkerColor(2);
- proj2->DrawCopy("SAME");
-
- AliPWG0Helper::CreateDividedProjections(track1, track2, "ze2");
- TH1* pt = gROOT->FindObject("track1_ze2_div_track2_ze2");
- new TCanvas; pt->Draw();
+ for (Int_t i=0; i<3; i++)
+ {
+ char c = 'x' + (char) i;
+
+ /*proj1 = track1->Project3D(Form("%ce2", c));
+ proj2 = track2->Project3D(Form("%ce2", c));
+ AliPWG0Helper::NormalizeToBinWidth(proj1);
+ AliPWG0Helper::NormalizeToBinWidth(proj2);
+
+ new TCanvas;
+ proj1->DrawCopy();
+ proj2->SetLineColor(2);
+ proj2->SetMarkerColor(2);
+ proj2->DrawCopy("SAME");*/
+
+ AliPWG0Helper::CreateDividedProjections(track1, track2, Form("%ce2", c));
+ TH1* pt = gROOT->FindObject(Form("track1_%ce2_div_track2_%ce2", c, c));
+ new TCanvas; pt->DrawCopy();
+ gPad->SetGridx(); gPad->SetGridy();
+ }
}
void MagnitudeOfCorrection(const char* fileName, const char* dirName = "dndeta", Float_t ptCut = 0.3)
fdNdEtaAnalysis->GetData()->PrintInfo(ptCut);
}
+Double_t ConvSigma1To2D(Double_t sigma)
+{
+ return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
+}
+
+Double_t ConvDistance1DTo2D(Double_t distance)
+{
+ return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
+}
+
+Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
+{
+ Double_t count = 0;
+
+ //nSigma = ConvSigma1To2D(nSigma);
+
+ for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+ for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+ {
+ Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
+ Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
+
+ Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
+
+ d = ConvDistance1DTo2D(d);
+
+ if (d < nSigma)
+ count += tracks->GetBinContent(x, y);
+ }
+
+ return count;
+}
+
+TH2F* Sigma2VertexGaussianTracksHist()
+{
+ TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
+
+ TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
+ gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
+
+ for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
+ for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
+ tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
+
+ //normalize
+ tracks->Scale(1.0 / tracks->Integral());
+
+ return tracks;
+}
+
+TH1F* Sigma2VertexGaussian()
+{
+ TH2F* tracks = Sigma2VertexGaussianTracksHist();
+
+ TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
+ canvas->Divide(2, 2);
+
+ canvas->cd(1);
+ tracks->Draw("COLZ");
+
+ TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 50, 0.05, 5.05);
+ for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
+ ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
+ ratio->SetMarkerStyle(21);
+
+ canvas->cd(2);
+ ratio->DrawCopy("P");
+
+ TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 4 sigma / % included n sigma", 50, 0.05, 5.05);
+ Double_t sigma3 = Sigma2VertexCount(tracks, 4);
+ for (Double_t nSigma = 0.1; nSigma < 5.05; nSigma += 0.1)
+ ratio2->Fill(nSigma, sigma3 / ratio->GetBinContent(ratio->FindBin(nSigma)));
+ ratio2->SetMarkerStyle(21);
+
+ canvas->cd(3);
+ ratio2->DrawCopy("P");
+
+ canvas->SaveAs("Sigma2Vertex.eps");
+
+ return ratio2;
+}
+
+TH1F** Sigma2VertexSimulation(const char* fileName = "systematics.root")
+{
+ TFile* file = TFile::Open(fileName);
+
+ TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertexTracks"));
+ TH1F* sigmavertexPrim = dynamic_cast<TH1F*> (file->Get("fSigmaVertexPrim"));
+ if (!sigmavertex || !sigmavertexPrim)
+ {
+ printf("Could not read histogram(s)\n");
+ return;
+ }
+
+ // calculate ratio
+ TH1F* ratio = new TH1F("sigmavertexsimulation_ratio", "sigmavertexsimulation_ratio;N#sigma;% included in 4 #sigma / % included in N#sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
+
+ // calculate contamination
+ TH1F* contamination = ratio->Clone("sigmavertexsimulation_contamination");
+ contamination->SetTitle("sigmavertexsimulation_contamination;N#sigma;1 + N_{secondaries} / N_{all}");
+
+ for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
+ {
+ ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(4)) / sigmavertex->GetBinContent(i));
+ contamination->SetBinContent(i, 1 + (sigmavertex->GetBinContent(i) - sigmavertexPrim->GetBinContent(i)) / sigmavertex->GetBinContent(i));
+ }
+
+ // print stats
+ for (Float_t sigma = 2.0; sigma < 5.25; sigma += 0.5)
+ {
+ Float_t error1 = 1 - ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma - 0.5));
+ Float_t error2 = -1 + ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma)) / ratio->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma + 0.5));
+ Float_t cont = -1 + contamination->GetBinContent(sigmavertex->GetXaxis()->FindBin(sigma));
+ Printf("%.2f sigma --> syst. error = - %.2f %% + %.2f %%, cont. = %.2f %%", sigma, error1 * 100, error2 * 100, cont * 100);
+ }
+
+ TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
+ canvas->Divide(2, 1);
+
+ canvas->cd(1);
+ sigmavertex->SetMarkerStyle(21);
+ sigmavertex->Draw("P");
+
+ canvas->cd(2);
+ ratio->SetMarkerStyle(21);
+ ratio->DrawCopy("P");
+
+ contamination->DrawCopy("SAME");
+
+ TH1F** returnContainer = new TH1F*[2];
+ returnContainer[0] = ratio;
+ returnContainer[1] = contamination;
+
+ return returnContainer;
+}
+
+void Sigma2VertexCompare(const char* fileName = "systematics.root")
+{
+ TH1F* ratio1 = Sigma2VertexGaussian();
+
+ TH1F** hists = Sigma2VertexSimulation(fileName);
+ TH1F* ratio2 = hists[0];
+ TH1F* contamination = hists[1];
+
+ ratio1->SetStats(kFALSE);
+ ratio2->SetStats(kFALSE);
+
+ ratio1->SetMarkerStyle(0);
+ ratio2->SetMarkerStyle(0);
+
+ ratio1->SetLineWidth(2);
+ ratio2->SetLineWidth(2);
+
+ TLegend* legend = new TLegend(0.6, 0.8, 0.95, 0.95);
+ legend->SetFillColor(0);
+ legend->AddEntry(ratio1, "Gaussian");
+ legend->AddEntry(ratio2, "Simulation");
+ legend->AddEntry(contamination, "1 + Contamination");
+
+ ratio2->SetTitle("");
+ ratio2->GetYaxis()->SetTitleOffset(1.5);
+ ratio2->GetXaxis()->SetRangeUser(2, 5);
+
+ TCanvas* canvas = new TCanvas("Sigma2VertexCompare", "Sigma2VertexCompare", 500, 500);
+ InitPad();
+
+ ratio2->SetMarkerStyle(21);
+ ratio1->SetMarkerStyle(22);
+
+ ratio2->GetYaxis()->SetRangeUser(0.8, 1.2);
+ ratio2->SetLineColor(kRed);
+ ratio2->SetMarkerColor(kRed);
+ ratio2->Draw("PL");
+ ratio1->Draw("SAMEPL");
+
+ contamination->Draw("SAME");
+
+ legend->Draw();
+
+ canvas->SaveAs("Sigma2VertexCompare.eps");
+}
dNdEtaAnalysis* fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kINEL, 1);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
TFile* file2 = TFile::Open(dataOutput, "RECREATE");
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kVertexReco, 1);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle);
+ fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0.3, AlidNdEtaCorrection::kTrack2Particle, 1);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
file->cd();
fdNdEtaAnalysis = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks");
fdNdEtaAnalysis->LoadHistograms("fdNdEtaAnalysisESD");
- fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone, 1);
//fdNdEtaAnalysis->DrawHistograms(kTRUE);
file2->cd();
fdNdEtaAnalysis->SaveHistograms();
//fdNdEtaAnalysis->Finish(dNdEtaCorrection, 0, AlidNdEtaCorrection::kTrack2Particle);
}
else
- fdNdEtaAnalysis->Finish(0, 0, AlidNdEtaCorrection::kNone);
+ fdNdEtaAnalysis->Finish(0, 0.3, AlidNdEtaCorrection::kNone);
fdNdEtaAnalysis->DrawHistograms(simple);
Int_t label = TMath::Abs(track->GetLabel());
if (stack->IsPhysicalPrimary(label) == kTRUE)
{
- if (label >= max)
- {
- Printf("Warning label %d is higher than number of primaries %d", label, max);
- continue;
- }
-
- if (fTrackCutsPrimaries->AcceptTrack(track))
- {
- primAcc[label]++;
- }
- else
- primRej[label]++;
+ if (label >= max)
+ {
+ Printf("Warning label %d is higher than number of primaries %d", label, max);
+ continue;
+ }
+
+ if (fTrackCutsPrimaries->AcceptTrack(track))
+ {
+ primAcc[label]++;
+ }
+ else
+ primRej[label]++;
}
else
for (Int_t i=0; i<max; i++) {
if (primAcc[i] == 1) {
- fPrimStats->Fill(1);
+ fPrimStats->Fill(1);
} else if (primAcc[i] > 1) {
- fPrimStats->Fill(2);
- fPrimStats->Fill(3, primAcc[i]);
+ fPrimStats->Fill(2);
+ fPrimStats->Fill(3, primAcc[i]);
}
if (primRej[i] > 0) {
- if (primAcc[i] == 0) {
- fPrimStats->Fill(4);
- fPrimStats->Fill(5, primRej[i]);
- } else if (primAcc[i] > 0) {
- fPrimStats->Fill(6);
- fPrimStats->Fill(7, primRej[i]);
- }
+ if (primAcc[i] == 0) {
+ fPrimStats->Fill(4);
+ fPrimStats->Fill(5, primRej[i]);
+ } else if (primAcc[i] > 0) {
+ fPrimStats->Fill(6);
+ fPrimStats->Fill(7, primRej[i]);
+ }
}
}
void draw()
{
- const char* files[] = { "trackCuts_normal.root", "trackCuts_increased.root", "trackCuts_decreased.root" };
+ //const char* files[] = { "trackCuts_normal.root", "trackCuts_increased.root", "trackCuts_decreased.root" };
+ const char* files[] =
+ { "Material-normal/trackCuts.root",
+ "Material-increased-mcvtx/trackCuts.root",
+ "Material-decreased-mcvtx/trackCuts.root" };
const char* titles[] = { "default geometry", "+ 10% material", "- 10% material" };
Int_t colors[] = { 1, 2, 4 };
+ Int_t markers[] = {24, 25, 26, 27 };
TCanvas* c = new TCanvas;
+ TCanvas* c2 = new TCanvas;
- TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
+ TLegend* legend = new TLegend(0.6, 0.5, 0.9, 0.9);
+ TLegend* legend2 = new TLegend(0.7, 0.7, 0.9, 0.9);
for (Int_t i=0; i<3; i++) {
- TFile::Open(files[i]);
-
+ if (!TFile::Open(files[i]))
+ return;
+
TH1* ptPrim = gFile->Get("fTrackCutsPrimaries/before_cuts/pt");
TH1* ptPrimCut = gFile->Get("fTrackCutsPrimaries/after_cuts/pt_cut");
TH1* ptSec = gFile->Get("fTrackCutsSecondaries/before_cuts/pt");
TH1* ptSecCut = gFile->Get("fTrackCutsSecondaries/after_cuts/pt_cut");
+ TH1* vertex = gFile->Get("fVertex");
+ Int_t nEvents = vertex->GetEntries();
+
ptPrim->Add(ptSec);
ptPrimCut->Add(ptSecCut);
ptPrimCut->Sumw2();
ptSec->Sumw2();
ptSecCut->Sumw2();
-
+
Printf("%s", titles[i]);
+ Printf("Total particles: %d", (Int_t) (ptPrim->GetEntries() + ptSec->GetEntries()));
+ Printf("Total particles/event: %.2f", (ptPrim->GetEntries() + ptSec->GetEntries()) / nEvents);
+ Printf("Primaries/event: %.2f", ptPrim->GetEntries() / nEvents);
+ Printf("Secondaries/event: %.2f", ptSec->GetEntries() / nEvents);
+ Printf("Primaries > 0.2 GeV/c/event: %.2f", ptPrim->Integral(ptPrim->GetXaxis()->FindBin(0.21), ptPrim->GetNbinsX()) / nEvents);
+ Printf("Secondaries > 0.2 GeV/c/event: %.2f", ptSec->Integral(ptSec->GetXaxis()->FindBin(0.21), ptSec->GetNbinsX()) / nEvents);
+ Printf("Primaries after cuts > 0.2 GeV/c/event: %.2f +- %.2f", ptPrimCut->Integral(ptPrimCut->GetXaxis()->FindBin(0.21), ptPrimCut->GetNbinsX()) / nEvents, TMath::Sqrt(ptPrimCut->Integral(ptPrimCut->GetXaxis()->FindBin(0.21), ptPrimCut->GetNbinsX())) / nEvents);
+ Printf("Secondaries after cuts > 0.2 GeV/c/event: %.2f +- %.2f", ptSecCut->Integral(ptSecCut->GetXaxis()->FindBin(0.21), ptSecCut->GetNbinsX()) / nEvents, TMath::Sqrt(ptSecCut->Integral(ptSecCut->GetXaxis()->FindBin(0.21), ptSecCut->GetNbinsX())) / nEvents);
Printf("%.2f %% secondaries before cuts", 100.0 * ptSec->GetEntries() / ptPrim->GetEntries());
Printf("%.2f %% secondaries after cuts", 100.0 * ptSecCut->GetEntries() / ptPrimCut->GetEntries());
Printf("");
- ptSec->Divide(ptSec, ptPrim, 1, 1, "B");
- ptSecCut->Divide(ptSecCut, ptPrimCut, 1, 1, "B");
-
+ ptPrim->SetLineColor(colors[i]);
+ ptPrimCut->SetLineColor(colors[i]);
ptSec->SetLineColor(colors[i]);
ptSecCut->SetLineColor(colors[i]);
+
+ ptPrim->SetMarkerColor(colors[i]);
+ ptPrimCut->SetMarkerColor(colors[i]);
+ ptSec->SetMarkerColor(colors[i]);
+ ptSecCut->SetMarkerColor(colors[i]);
+
+ ptPrim->SetStats(kFALSE);
+ ptPrim->SetTitle("");
+ ptPrim->GetYaxis()->SetTitle("N");
ptSec->SetStats(kFALSE);
+ ptSec->SetTitle("");
+
+ ptPrim->SetMarkerStyle(markers[0]);
+ ptPrimCut->SetMarkerStyle(markers[1]);
+ ptSec->SetMarkerStyle(markers[2]);
+ ptSecCut->SetMarkerStyle(markers[3]);
+ if (i == 0) {
+ legend->AddEntry(ptPrim->Clone(), "Primaries");
+ legend->AddEntry(ptPrimCut->Clone(), "Primaries after cuts");
+ legend->AddEntry(ptSec->Clone(), "Secondaries");
+ legend->AddEntry(ptSecCut->Clone(), "Secondaries after cuts");
+ }
+
+ ptPrim->GetXaxis()->SetRangeUser(0, 2);
ptSec->GetXaxis()->SetRangeUser(0, 2);
+ //ptPrim->GetYaxis()->SetRangeUser(1e-5, ptSec->GetMaximum() * 1.1);
+
+ c->cd();
+ ptPrim->DrawCopy((i > 0) ? "SAME" : "");
+ ptPrimCut->DrawCopy("SAME");
+ ptSec->DrawCopy("SAME");
+ ptSecCut->DrawCopy("SAME");
+
+ ptSec->Divide(ptSec, ptPrim, 1, 1, "B");
+ ptSecCut->Divide(ptSecCut, ptPrimCut, 1, 1, "B");
+
+ ptSec->SetMarkerStyle(1);
+ ptSecCut->SetMarkerStyle(1);
+
ptSec->GetYaxis()->SetRangeUser(0, 1);
- ptSec->SetTitle("");
ptSec->GetYaxis()->SetTitle("N_{Secondaries} / N_{All}");
-
+
+ c2->cd();
ptSec->DrawCopy((i > 0) ? "SAME" : "");
ptSecCut->DrawCopy("SAME");
legend->AddEntry(ptSec, titles[i]);
+ legend2->AddEntry(ptSec, titles[i]);
}
+ c->cd();
legend->Draw();
+ c->SaveAs("changedmaterial_absolute.gif");
- c->SaveAs("secondaries_changedmaterial.gif");
+ c2->cd();
+ legend2->Draw();
+ c2->SaveAs("changedmaterial_relative.gif");
}
void drawStats(const char* fileName = "trackCuts.root")
Printf("UNEXPECTED: %f != %f", fPrimStats->GetBinContent(5) + fPrimStats->GetBinContent(7), fStatsPrim->GetBinContent(2));
Float_t notLostPrimaries = 100.0 * fPrimStats->GetBinContent(7) / (fPrimStats->GetBinContent(5) + fPrimStats->GetBinContent(7));
-
+
Printf("Accepted %d out of %d primary tracks (%.2f %%)", tracksPrimAc, tracksPrim, 100.0 * tracksPrimAc / tracksPrim);
+ Printf("Per Event: %.2f out of %.2f", (Float_t) tracksPrimAc / eventsVertex, (Float_t) tracksPrim / eventsVertex);
Printf("Among the non accepted ones %.2f %% are from primaries that have been found with other tracks", notLostPrimaries);
Printf("Accepted %d out of %d secondary tracks (%.2f %%)", tracksSecAc, tracksSec, 100.0 * tracksSecAc / tracksSec);
+ Printf("Per Event: %.2f out of %.2f", (Float_t) tracksSecAc / eventsVertex, (Float_t) tracksSec / eventsVertex);
Printf("Before cuts: %.2f %% of the tracks are secondaries", 100.0 * tracksSec / (tracksPrim + tracksSec));
Printf("After cuts: %.2f %% of the tracks are secondaries", 100.0 * tracksSecAc / (tracksPrimAc + tracksSecAc));
}
-void ComparePlots(const char* fileName1, const char* fileName2, const char* dir = "AliESDtrackCuts/before_cuts")
+void ComparePlots(const char* fileName1, const char* fileName2, const char* dirName1 = "AliESDtrackCuts/before_cuts", const char* dirName2 = 0)
{
file1 = TFile::Open(fileName1);
file2 = TFile::Open(fileName2);
- dir1 = file1->GetDirectory(dir);
- dir2 = file2->GetDirectory(dir);
+ if (!dirName2)
+ dirName1 = dirName2;
+
+ dir1 = file1->GetDirectory(dirName1);
+ dir2 = file2->GetDirectory(dirName2);
keyList = dir1->GetListOfKeys();
TIter iter(keyList);
TString name(key->GetName());
- c = new TCanvas(key->GetName(), key->GetName(), 600, 400);
+ c = new TCanvas(Form("%s_canvas", key->GetName()), key->GetName(), 600, 400);
+
hist1->Draw();
hist2->SetLineColor(2);
hist2->Draw("SAME");
+ /* this does not work, although given here: http://pcroot.cern.ch/root/html/THistPainter.html ;-)
+ TPaveStats *st = (TPaveStats*) hist1->FindObject("stats");
+ st->SetTextColor(2);
+ st->SetY1NDC(st->GetY1NDC() - 0.2);
+ st->SetY2NDC(st->GetY2NDC() - 0.2);
+ */
+
Float_t min = 0.9 * TMath::Min(hist1->GetMinimum(), hist2->GetMinimum());
if (name.Contains("cov") || name.Contains("dXY") || name.Contains("dZ")) {
hist1->GetYaxis()->SetRangeUser(min, 1.1 * TMath::Max(hist1->GetMaximum(), hist2->GetMaximum()));
c->SaveAs(Form("%s.gif", key->GetName()));
-
+
//break;
}
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff