From 0a17397874c92ec0002dcbe1cc5352eeee6c0b73 Mon Sep 17 00:00:00 2001 From: jgrosseo Date: Wed, 24 Jan 2007 16:55:06 +0000 Subject: [PATCH] introducing multiplicity measurement with the ITS or TPC adding class that stores raw input, correction maps and function to correct the spectrum --- PWG0/PWG0baseLinkDef.h | 1 + PWG0/dNdEta/AliMultiplicityCorrection.cxx | 668 ++++++++++++++++++++++ PWG0/dNdEta/AliMultiplicityCorrection.h | 76 +++ PWG0/dNdEta/AliMultiplicityMCSelector.cxx | 182 +++++- PWG0/dNdEta/AliMultiplicityMCSelector.h | 12 +- PWG0/dNdEta/runMultiplicitySelector.C | 30 +- PWG0/libPWG0base.pkg | 3 +- 7 files changed, 928 insertions(+), 44 deletions(-) create mode 100644 PWG0/dNdEta/AliMultiplicityCorrection.cxx create mode 100644 PWG0/dNdEta/AliMultiplicityCorrection.h diff --git a/PWG0/PWG0baseLinkDef.h b/PWG0/PWG0baseLinkDef.h index 315ea262285..3f4137c1b54 100644 --- a/PWG0/PWG0baseLinkDef.h +++ b/PWG0/PWG0baseLinkDef.h @@ -21,6 +21,7 @@ #pragma link C++ class AliCorrection+; +#pragma link C++ class AliMultiplicityCorrection+; #endif diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.cxx b/PWG0/dNdEta/AliMultiplicityCorrection.cxx new file mode 100644 index 00000000000..c3f79167578 --- /dev/null +++ b/PWG0/dNdEta/AliMultiplicityCorrection.cxx @@ -0,0 +1,668 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +// This class is used to store correction maps, raw input and results of the multiplicity +// measurement with the ITS or TPC +// It also contains functions to correct the spectrum using different methods. +// +// Author: Jan.Fiete.Grosse-Oetringhaus@cern.ch + +#include "AliMultiplicityCorrection.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +ClassImp(AliMultiplicityCorrection) + +const Int_t AliMultiplicityCorrection::fgMaxParams = 110; +TH1* AliMultiplicityCorrection::fCurrentMinuitESD = 0; +TH1* AliMultiplicityCorrection::fCurrentMinuitCorrelation = 0; + +//____________________________________________________________________ +AliMultiplicityCorrection::AliMultiplicityCorrection() : + TNamed() +{ + // + // default constructor + // + + for (Int_t i = 0; i < kESDHists; ++i) + fMultiplicityESD[i] = 0; + + for (Int_t i = 0; i < kMCHists; ++i) + fMultiplicityMC[i] = 0; + + for (Int_t i = 0; i < kCorrHists; ++i) + { + fCorrelation[i] = 0; + fMultiplicityESDCorrected[i] = 0; + } +} + +//____________________________________________________________________ +AliMultiplicityCorrection::AliMultiplicityCorrection(const Char_t* name, const Char_t* title) : + TNamed(name, title) +{ + // + // named constructor + // + + // do not add this hists to the directory + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + Float_t binLimitsVtx[] = {-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10}; + Float_t binLimitsN[] = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, + 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, + 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, + 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, + 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, + 50.5, 55.5, 60.5, 65.5, 70.5, 75.5, 80.5, 85.5, 90.5, 95.5, + 100.5, 105.5, 110.5, 115.5, 120.5, 125.5, 130.5, 135.5, 140.5, 145.5, + 150.5, 160.5, 170.5, 180.5, 190.5, 200.5, 210.5, 220.5, 230.5, 240.5, + 250.5, 275.5, 300.5, 325.5, 350.5, 375.5, 400.5, 425.5, 450.5, 475.5, + 500.5, 525.5, 550.5, 575.5, 600.5, 625.5, 650.5, 675.5, 700.5, 725.5, + 750.5, 775.5, 800.5, 825.5, 850.5, 875.5, 900.5, 925.5, 950.5, 975.5, + 1000.5 }; + + #define BINNING 110, binLimitsN + + for (Int_t i = 0; i < kESDHists; ++i) + fMultiplicityESD[i] = new TH2F(Form("fMultiplicityESD%d", i), "fMultiplicityESD;vtx-z;Ntracks;Count", 10, binLimitsVtx, BINNING); + + for (Int_t i = 0; i < kMCHists; ++i) + { + fMultiplicityMC[i] = dynamic_cast (fMultiplicityESD[0]->Clone(Form("fMultiplicityMC%d", i))); + fMultiplicityMC[i]->SetTitle("fMultiplicityMC"); + } + + for (Int_t i = 0; i < kCorrHists; ++i) + { + fCorrelation[i] = new TH3F(Form("fCorrelation%d", i), "fCorrelation;vtx-z;Npart;Ntracks", 10, binLimitsVtx, BINNING, BINNING); + fMultiplicityESDCorrected[i] = new TH1F(Form("fMultiplicityESDCorrected%d", i), "fMultiplicityESDCorrected;Npart;dN/dN", BINNING); + } + + TH1::AddDirectory(oldStatus); +} + +//____________________________________________________________________ +AliMultiplicityCorrection::~AliMultiplicityCorrection() +{ + // + // Destructor + // +} + +//____________________________________________________________________ +Long64_t AliMultiplicityCorrection::Merge(TCollection* list) +{ + // Merge a list of AliMultiplicityCorrection objects with this (needed for + // PROOF). + // Returns the number of merged objects (including this). + + if (!list) + return 0; + + if (list->IsEmpty()) + return 1; + + TIterator* iter = list->MakeIterator(); + TObject* obj; + + // collections of all histograms + TList collections[kESDHists+kMCHists+kCorrHists*2]; + + Int_t count = 0; + while ((obj = iter->Next())) { + + AliMultiplicityCorrection* entry = dynamic_cast (obj); + if (entry == 0) + continue; + + for (Int_t i = 0; i < kESDHists; ++i) + collections[i].Add(entry->fMultiplicityESD[i]); + + for (Int_t i = 0; i < kMCHists; ++i) + collections[kESDHists+i].Add(entry->fMultiplicityMC[i]); + + for (Int_t i = 0; i < kCorrHists; ++i) + collections[kESDHists+kMCHists+i].Add(entry->fCorrelation[i]); + + for (Int_t i = 0; i < kCorrHists; ++i) + collections[kESDHists+kMCHists+kCorrHists+i].Add(entry->fMultiplicityESDCorrected[i]); + + count++; + } + + for (Int_t i = 0; i < kESDHists; ++i) + fMultiplicityESD[i]->Merge(&collections[i]); + + for (Int_t i = 0; i < kMCHists; ++i) + fMultiplicityMC[i]->Merge(&collections[kESDHists+i]); + + for (Int_t i = 0; i < kCorrHists; ++i) + fCorrelation[i]->Merge(&collections[kESDHists+kMCHists+i]); + + for (Int_t i = 0; i < kCorrHists; ++i) + fMultiplicityESDCorrected[i]->Merge(&collections[kESDHists+kMCHists+kCorrHists+i]); + + delete iter; + + return count+1; +} + +//____________________________________________________________________ +Bool_t AliMultiplicityCorrection::LoadHistograms(const Char_t* dir) +{ + // + // loads the histograms from a file + // if dir is empty a directory with the name of this object is taken (like in SaveHistogram) + // + + if (!dir) + dir = GetName(); + + if (!gDirectory->cd(dir)) + return kFALSE; + + // TODO memory leak. old histograms needs to be deleted. + + Bool_t success = kTRUE; + + for (Int_t i = 0; i < kESDHists; ++i) + { + fMultiplicityESD[i] = dynamic_cast (gDirectory->Get(fMultiplicityESD[i]->GetName())); + if (!fMultiplicityESD[i]) + success = kFALSE; + } + + for (Int_t i = 0; i < kMCHists; ++i) + { + fMultiplicityMC[i] = dynamic_cast (gDirectory->Get(fMultiplicityMC[i]->GetName())); + if (!fMultiplicityMC[i]) + success = kFALSE; + } + + for (Int_t i = 0; i < kCorrHists; ++i) + { + fCorrelation[i] = dynamic_cast (gDirectory->Get(fCorrelation[i]->GetName())); + if (!fCorrelation[i]) + success = kFALSE; + fMultiplicityESDCorrected[i] = dynamic_cast (gDirectory->Get(fMultiplicityESDCorrected[i]->GetName())); + if (!fMultiplicityESDCorrected[i]) + success = kFALSE; + } + + gDirectory->cd(".."); + + return success; +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::SaveHistograms() +{ + // + // saves the histograms + // + + gDirectory->mkdir(GetName()); + gDirectory->cd(GetName()); + + for (Int_t i = 0; i < kESDHists; ++i) + if (fMultiplicityESD[i]) + fMultiplicityESD[i]->Write(); + + for (Int_t i = 0; i < kMCHists; ++i) + if (fMultiplicityMC[i]) + fMultiplicityMC[i]->Write(); + + for (Int_t i = 0; i < kCorrHists; ++i) + { + if (fCorrelation[i]) + fCorrelation[i]->Write(); + if (fMultiplicityESDCorrected[i]) + fMultiplicityESDCorrected[i]->Write(); + } + + gDirectory->cd(".."); +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll) +{ + // + // Fills an event from MC + // + + fMultiplicityMC[0]->Fill(vtx, generated05); + fMultiplicityMC[1]->Fill(vtx, generated10); + fMultiplicityMC[2]->Fill(vtx, generated15); + fMultiplicityMC[3]->Fill(vtx, generated20); + fMultiplicityMC[4]->Fill(vtx, generatedAll); +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) +{ + // + // Fills an event from ESD + // + + fMultiplicityESD[0]->Fill(vtx, measured05); + fMultiplicityESD[1]->Fill(vtx, measured10); + fMultiplicityESD[2]->Fill(vtx, measured15); + fMultiplicityESD[3]->Fill(vtx, measured20); +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20) +{ + // + // Fills an event into the correlation map with the information from MC and ESD + // + + fCorrelation[0]->Fill(vtx, generated05, measured05); + fCorrelation[1]->Fill(vtx, generated10, measured10); + fCorrelation[2]->Fill(vtx, generated15, measured15); + fCorrelation[3]->Fill(vtx, generated20, measured20); + + fCorrelation[4]->Fill(vtx, generatedAll, measured05); + fCorrelation[5]->Fill(vtx, generatedAll, measured10); + fCorrelation[6]->Fill(vtx, generatedAll, measured15); + fCorrelation[7]->Fill(vtx, generatedAll, measured20); +} + +//____________________________________________________________________ +Double_t AliMultiplicityCorrection::MinuitHomogenityPol0(Double_t *params) +{ + // homogenity term for minuit fitting + // pure function of the parameters + // prefers constant function (pol0) + + Double_t chi2 = 0; + + for (Int_t i=1; iGetBinWidth(i+1); + Double_t left = params[i-1] / fCurrentMinuitESD->GetBinWidth(i); + + Double_t diff = (right - left) / right; + chi2 += diff * diff; + } + + return chi2; +} + +//____________________________________________________________________ +Double_t AliMultiplicityCorrection::MinuitHomogenityPol1(Double_t *params) +{ + // homogenity term for minuit fitting + // pure function of the parameters + // prefers linear function (pol1) + + Double_t chi2 = 0; + + for (Int_t i=2; iGetBinWidth(i+1); + Double_t middle = params[i-1] / fCurrentMinuitESD->GetBinWidth(i); + Double_t left = params[i-2] / fCurrentMinuitESD->GetBinWidth(i-1); + + Double_t der1 = (right - middle) / fCurrentMinuitESD->GetBinWidth(i); + Double_t der2 = (middle - left) / fCurrentMinuitESD->GetBinWidth(i-1); + + Double_t diff = der1 - der2; + + // TODO give additonal weight to big bins + chi2 += diff * diff * fCurrentMinuitESD->GetBinWidth(i) * fCurrentMinuitESD->GetBinWidth(i-1); + + //printf("%d --> %f\n", i, diff); + } + + return chi2 / 1e5 / 2; +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t) +{ + // + // fit function for minuit + // + + // TODO take errors into account + + static Int_t callCount = 0; + + Double_t chi2FromFit = 0; + + // loop over multiplicity (x axis of fMultiplicityESD) + for (Int_t i=1; i<=fCurrentMinuitESD->GetNbinsX(); ++i) + { + if (i > fCurrentMinuitCorrelation->GetNbinsY()) + break; + + Double_t sum = 0; + //Double_t error = 0; + // loop over generated particles (x axis of fCorrelation) that resulted in reconstruction of i tracks + for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsX(); ++j) + { + if (j > fgMaxParams) + break; + + sum += fCurrentMinuitCorrelation->GetBinContent(j, i) * params[j-1]; + + //if (params[j-1] > 0) + // error += fCurrentMinuitCorrelation->GetBinError(j, i) * fCurrentMinuitCorrelation->GetBinError(j, i) * params[j-1]; + //printf("%f ", sum); + } + + Double_t diff = fCurrentMinuitESD->GetBinContent(i) - sum; + if (fCurrentMinuitESD->GetBinContent(i) > 0) + diff /= fCurrentMinuitESD->GetBinContent(i); + else + diff /= fCurrentMinuitESD->Integral(); + //error = TMath::Sqrt(error) + fCurrentMinuitESD->GetBinError(i); + //if (error <= 0) + // error = 1; + + //Double_t tmp = diff / error; + //chi2 += tmp * tmp; + chi2FromFit += diff * diff; + + //printf("\nExpected sum = %f; Diff for bin %d is %f\n**********************************\n", fCurrentMinuitESD->GetBinContent(i), i, diff); + //printf("Diff for bin %d is %f\n", i, diff); + } + + Double_t penaltyVal = MinuitHomogenityPol1(params); + + chi2 = chi2FromFit * chi2FromFit + penaltyVal * penaltyVal; + + if ((callCount++ % 100) == 0) + printf("%f %f --> %f\n", chi2FromFit, penaltyVal, chi2); +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace) +{ + // + // correct spectrum using minuit chi2 method + // + + Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); + Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); + + fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY(); + fCurrentMinuitESD->Sumw2(); + + // empty under/overflow bins in x, otherwise Project3D takes them into account + TH3* hist = fCorrelation[correlationID]; + for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) + { + for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) + { + hist->SetBinContent(0, y, z, 0); + hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); + } + } + + fCurrentMinuitCorrelation = hist->Project3D("zy"); + fCurrentMinuitCorrelation->Sumw2(); + + // normalize correction for given nPart + for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i) + { + Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY()); + if (sum <= 0) + continue; + for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j) + { + // npart sum to 1 + fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum); + fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum); + } + } + + // small hack to get around charge conservation for full phase space ;-) + /*if (fullPhaseSpace) + { + for (Int_t i=2; i<=50; i+=2) + for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j) + fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i-1, j)); + }*/ + + TCanvas* canvas1 = new TCanvas("ApplyMinuitFit", "ApplyMinuitFit", 800, 400); + canvas1->Divide(2, 1); + canvas1->cd(1); + fCurrentMinuitESD->DrawCopy(); + canvas1->cd(2); + fCurrentMinuitCorrelation->DrawCopy("COLZ"); + + // Initialize TMinuit via generic fitter interface + TVirtualFitter *minuit = TVirtualFitter::Fitter(0, fgMaxParams); + + minuit->SetFCN(MinuitFitFunction); + + Double_t results[fgMaxParams]; + + //TH1* proj = fMultiplicityMC[mcTarget]->ProjectionY(); + for (Int_t i=0; iGetBinContent(i+1); + //results[i] = proj->GetBinContent(i+1); + minuit->SetParameter(i, Form("param%d", i), results[i], ((results[i] > 1) ? (results[i] / 10) : 0), 0, fCurrentMinuitESD->GetMaximum() * 100); + } + + Int_t dummy; + Double_t chi2; + MinuitFitFunction(dummy, 0, chi2, results, 0); + printf("Chi2 of right solution is = %f\n", chi2); + + //return; + + Double_t arglist[100]; + arglist[0] = 100000; + //minuit->ExecuteCommand("SET PRINT", arglist, 1); + minuit->ExecuteCommand("MIGRAD", arglist, 1); + //minuit->ExecuteCommand("MINOS", arglist, 0); + + for (Int_t i=0; iGetParameter(i); + fMultiplicityESDCorrected[correlationID]->SetBinContent(i+1, results[i]); + fMultiplicityESDCorrected[correlationID]->SetBinError(i+1, 0); + } + + printf("Penalty is %f\n", MinuitHomogenityPol1(results)); + + fMultiplicityESDCorrected[correlationID]->GetXaxis()->SetRangeUser(0, fgMaxParams); + + DrawComparison(mcTarget, correlationID); +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::NormalizeToBinWidth(TH1* hist) +{ + // + // normalizes a 1-d histogram to its bin width + // + + for (Int_t i=1; i<=hist->GetNbinsX(); ++i) + { + hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); + hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); + } +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::NormalizeToBinWidth(TH2* hist) +{ + // + // normalizes a 2-d histogram to its bin width (x width * y width) + // + + for (Int_t i=1; i<=hist->GetNbinsX(); ++i) + for (Int_t j=1; j<=hist->GetNbinsY(); ++j) + { + Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); + hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); + hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); + } +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::ApplyMinuitFitAll() +{ + // + // fit all eta ranges + // + + for (Int_t i=0; iDivide(3, 2); + for (Int_t i = 0; i < kESDHists; ++i) + { + canvas1->cd(i+1); + fMultiplicityESD[i]->DrawCopy("COLZ"); + printf("%d --> %f\n", i, (Float_t) fMultiplicityESD[i]->ProjectionY()->GetMean()); + } + + printf("MC:\n"); + + TCanvas* canvas2 = new TCanvas("fMultiplicityMC", "fMultiplicityMC", 900, 600); + canvas2->Divide(3, 2); + for (Int_t i = 0; i < kMCHists; ++i) + { + canvas2->cd(i+1); + fMultiplicityMC[i]->DrawCopy("COLZ"); + printf("%d --> %f\n", i, fMultiplicityMC[i]->ProjectionY()->GetMean()); + } + + TCanvas* canvas3 = new TCanvas("fCorrelation", "fCorrelation", 900, 900); + canvas3->Divide(3, 3); + for (Int_t i = 0; i < kCorrHists; ++i) + { + canvas3->cd(i+1); + TH1* proj = fCorrelation[i]->Project3D("zy"); + proj->DrawCopy("COLZ"); + } +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::DrawComparison(Int_t mcID, Int_t esdCorrId) +{ + TString name; + name.Form("DrawComparison_%d_%d", mcID, esdCorrId); + + TCanvas* canvas1 = new TCanvas(name, name, 900, 300); + canvas1->Divide(3, 1); + + canvas1->cd(1); + TH1* proj = fMultiplicityMC[mcID]->ProjectionY(); + NormalizeToBinWidth(proj); + NormalizeToBinWidth(fMultiplicityESDCorrected[esdCorrId]); + proj->GetXaxis()->SetRangeUser(0, 200); + proj->DrawCopy("HIST"); + gPad->SetLogy(); + + fMultiplicityESDCorrected[esdCorrId]->SetMarkerStyle(3); + fMultiplicityESDCorrected[esdCorrId]->DrawCopy("SAMEP"); + + canvas1->cd(2); + fMultiplicityESDCorrected[esdCorrId]->DrawCopy(""); + + canvas1->cd(3); + TH1* clone = dynamic_cast (proj->Clone("clone")); + clone->Divide(fMultiplicityESDCorrected[esdCorrId]); + clone->SetTitle("Ratio;Npart;MC/ESD"); + clone->GetYaxis()->SetRangeUser(0.8, 1.2); + clone->GetXaxis()->SetRangeUser(0, 200); + clone->DrawCopy("HIST"); + + /*TLegend* legend = new TLegend(0.6, 0.7, 0.85, 0.85); + legend->AddEntry(fMultiplicityESDCorrected, "ESD corrected"); + legend->AddEntry(fMultiplicityMC, "MC"); + legend->AddEntry(fMultiplicityESD, "ESD"); + legend->Draw();*/ +} + +//____________________________________________________________________ +void AliMultiplicityCorrection::ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace) +{ + // + // correct spectrum using bayesian method + // + + Int_t correlationID = inputRange + ((fullPhaseSpace == kFALSE) ? 0 : 4); + Int_t mcTarget = ((fullPhaseSpace == kFALSE) ? inputRange : 4); + + fCurrentMinuitESD = fMultiplicityESD[inputRange]->ProjectionY(); + fCurrentMinuitESD->Sumw2(); + + // empty under/overflow bins in x, otherwise Project3D takes them into account + TH3* hist = fCorrelation[correlationID]; + for (Int_t y=1; y<=hist->GetYaxis()->GetNbins(); ++y) + { + for (Int_t z=1; z<=hist->GetZaxis()->GetNbins(); ++z) + { + hist->SetBinContent(0, y, z, 0); + hist->SetBinContent(hist->GetXaxis()->GetNbins()+1, y, z, 0); + } + } + + fCurrentMinuitCorrelation = hist->Project3D("zy"); + fCurrentMinuitCorrelation->Sumw2(); + + // normalize correction for given nPart + for (Int_t i=1; i<=fCurrentMinuitCorrelation->GetNbinsX(); ++i) + { + Double_t sum = fCurrentMinuitCorrelation->Integral(i, i, 1, fCurrentMinuitCorrelation->GetNbinsY()); + if (sum <= 0) + continue; + for (Int_t j=1; j<=fCurrentMinuitCorrelation->GetNbinsY(); ++j) + { + // npart sum to 1 + fCurrentMinuitCorrelation->SetBinContent(i, j, fCurrentMinuitCorrelation->GetBinContent(i, j) / sum); + fCurrentMinuitCorrelation->SetBinError(i, j, fCurrentMinuitCorrelation->GetBinError(i, j) / sum); + } + } + + // TODO to be implemented +} diff --git a/PWG0/dNdEta/AliMultiplicityCorrection.h b/PWG0/dNdEta/AliMultiplicityCorrection.h new file mode 100644 index 00000000000..350304223c8 --- /dev/null +++ b/PWG0/dNdEta/AliMultiplicityCorrection.h @@ -0,0 +1,76 @@ +/* $Id$ */ + +#ifndef ALIMULTIPLICITYCORRECTION_H +#define ALIMULTIPLICITYCORRECTION_H + +#include "TNamed.h" + +// +// class that contains the correction matrix and the functions for +// correction the multiplicity spectrum +// + +class TH1; +class TH2; +class TH1F; +class TH2F; +class TH3F; + +class AliMultiplicityCorrection : public TNamed { + public: + AliMultiplicityCorrection(); + AliMultiplicityCorrection(const Char_t* name, const Char_t* title); + virtual ~AliMultiplicityCorrection(); + + virtual Long64_t Merge(TCollection* list); + + void FillMeasured(Float_t vtx, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20); + void FillGenerated(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll); + + void FillCorrection(Float_t vtx, Int_t generated05, Int_t generated10, Int_t generated15, Int_t generated20, Int_t generatedAll, Int_t measured05, Int_t measured10, Int_t measured15, Int_t measured20); + + Bool_t LoadHistograms(const Char_t* dir); + void SaveHistograms(); + void DrawHistograms(); + void DrawComparison(Int_t mcID, Int_t esdCorrId); + + void ApplyMinuitFit(Int_t inputRange, Bool_t fullPhaseSpace); + void ApplyMinuitFitAll(); + + void ApplyBayesianMethod(Int_t inputRange, Bool_t fullPhaseSpace); + + TH2F* GetMultiplicityESD(Int_t i) { return fMultiplicityESD[i]; } + TH2F* GetMultiplicityMC(Int_t i) { return fMultiplicityMC[i]; } + TH3F* GetCorrelation(Int_t i) { return fCorrelation[i]; } + + TH1F* GetMultiplicityESDCorrected(Int_t i) { return fMultiplicityESDCorrected[i]; } + + static void NormalizeToBinWidth(TH1* hist); + static void NormalizeToBinWidth(TH2* hist); + + protected: + enum { kESDHists = 4, kMCHists = 5, kCorrHists = 8 }; + + static const Int_t fgMaxParams; // number of fit params + + static Double_t MinuitHomogenityPol0(Double_t *params); + static Double_t MinuitHomogenityPol1(Double_t *params); + static void MinuitFitFunction(Int_t&, Double_t*, Double_t& chi2, Double_t *params, Int_t); + + static TH1* fCurrentMinuitESD; // static variable to be accessed by MINUIT + static TH1* fCurrentMinuitCorrelation; // static variable to be accessed by MINUIT + + TH2F* fMultiplicityESD[kESDHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2 (0..3) + TH2F* fMultiplicityMC[kMCHists]; // multiplicity histogram: vtx vs multiplicity; array: |eta| < 0.5, 1, 1.5, 2, inf (0..4) + + TH3F* fCorrelation[kCorrHists]; // vtx vs. (gene multiplicity) vs. (meas multiplicity); array: |eta| < 0.5, 1, 1.5, 2 (0..3 and 4..7), the first corrects to the eta range itself, the second to full phase space + TH1F* fMultiplicityESDCorrected[kCorrHists]; // corrected histograms + + private: + AliMultiplicityCorrection(const AliMultiplicityCorrection&); + AliMultiplicityCorrection& operator=(const AliMultiplicityCorrection&); + + ClassDef(AliMultiplicityCorrection, 1); +}; + +#endif diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.cxx b/PWG0/dNdEta/AliMultiplicityMCSelector.cxx index 64352309fa4..effafbbbf4e 100644 --- a/PWG0/dNdEta/AliMultiplicityMCSelector.cxx +++ b/PWG0/dNdEta/AliMultiplicityMCSelector.cxx @@ -8,24 +8,26 @@ #include #include #include -#include #include +#include #include #include #include #include +#include +#include +#include #include "esdTrackCuts/AliESDtrackCuts.h" #include "AliPWG0Helper.h" +#include "dNdEta/AliMultiplicityCorrection.h" ClassImp(AliMultiplicityMCSelector) AliMultiplicityMCSelector::AliMultiplicityMCSelector() : AliSelectorRL(), - fMultiplicityESD(0), - fMultiplicityMC(0), - fCorrelation(0), + fMultiplicity(0), fEsdTrackCuts(0) { // @@ -76,10 +78,25 @@ void AliMultiplicityMCSelector::SlaveBegin(TTree* tree) ReadUserObjects(tree); - fMultiplicityESD = new TH1F("fMultiplicityESD", "fMultiplicityESD;Ntracks;Count", 201, -0.5, 200.5); - fMultiplicityMC = new TH1F("fMultiplicityMC", "fMultiplicityMC;Npart;Count", 201, -0.5, 200.5); + fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); +} + +void AliMultiplicityMCSelector::Init(TTree* tree) +{ + // read the user objects + + AliSelectorRL::Init(tree); + + // TODO Enable only the needed branches + /*if (tree) + { + tree->SetBranchStatus("*", 0); + tree->SetBranchStatus("fTriggerMask", 1); + tree->SetBranchStatus("fSPDVertex*", 1); + tree->SetBranchStatus("fTracks.fLabel", 1); - fCorrelation = new TH2F("fCorrelation", "fCorrelation;Npart;Ntracks", 201, -0.5, 200.5, 201, -0.5, 200.5); + AliESDtrackCuts::EnableNeededBranches(tree); + }*/ } Bool_t AliMultiplicityMCSelector::Process(Long64_t entry) @@ -118,6 +135,13 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry) return kFALSE; } + AliHeader* header = GetHeader(); + if (!header) + { + AliDebug(AliLog::kError, "Header not available"); + return kFALSE; + } + AliStack* stack = GetStack(); if (!stack) { @@ -131,16 +155,37 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry) if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE) return kTRUE; - // TODO put pt cut + // get the ESD vertex + const AliESDVertex* vtxESD = fESD->GetVertex(); + Double_t vtx[3]; + vtxESD->GetXYZ(vtx); + + // get the MC vertex + AliGenEventHeader* genHeader = header->GenEventHeader(); + TArrayF vtxMC(3); + genHeader->PrimaryVertex(vtxMC); + + // get tracklets + const AliMultiplicity* mult = fESD->GetMultiplicity(); + if (!mult) + { + AliDebug(AliLog::kError, "AliMultiplicity not available"); + return kFALSE; + } - // get number of "good" tracks from ESD - Int_t nESDTracks = fEsdTrackCuts->CountAcceptedTracks(fESD); - fMultiplicityESD->Fill(nESDTracks); + //const Float_t kPtCut = 0.3; // get number of tracks from MC // loop over mc particles Int_t nPrim = stack->GetNprimary(); - Int_t nMCTracks = 0; + + // tracks in different eta ranges + Int_t nMCTracks05 = 0; + Int_t nMCTracks10 = 0; + Int_t nMCTracks15 = 0; + Int_t nMCTracks20 = 0; + Int_t nMCTracksAll = 0; + for (Int_t iMc = 0; iMc < nPrim; ++iMc) { AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc)); @@ -156,17 +201,98 @@ Bool_t AliMultiplicityMCSelector::Process(Long64_t entry) if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) continue; - Float_t eta = particle->Eta(); + //if (particle->Pt() < kPtCut) + // continue; - if (TMath::Abs(eta) > 0.9) - continue; + if (TMath::Abs(particle->Eta()) < 0.5) + nMCTracks05++; - nMCTracks++; + if (TMath::Abs(particle->Eta()) < 1.0) + nMCTracks10++; + + if (TMath::Abs(particle->Eta()) < 1.5) + nMCTracks15++; + + if (TMath::Abs(particle->Eta()) < 2.0) + nMCTracks20++; + + nMCTracksAll++; }// end of mc particle - fMultiplicityMC->Fill(nMCTracks); + // FAKE: put back vtxMC[2] + fMultiplicity->FillGenerated(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll); + + Int_t nESDTracks05 = 0; + Int_t nESDTracks10 = 0; + Int_t nESDTracks15 = 0; + Int_t nESDTracks20 = 0; + + // get multiplicity from ITS tracklets + for (Int_t i=0; iGetNumberOfTracklets(); ++i) + { + //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); + + // this removes non-tracklets. Very bad solution. SPD guys are working on better solution... + if (mult->GetDeltaPhi(i) < -1000) + continue; + + Float_t theta = mult->GetTheta(i); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + + if (TMath::Abs(eta) < 0.5) + nESDTracks05++; + + if (TMath::Abs(eta) < 1.0) + nESDTracks10++; + + if (TMath::Abs(eta) < 1.5) + nESDTracks15++; + + if (TMath::Abs(eta) < 2.0) + nESDTracks20++; + } + + // get multiplicity from ESD tracks + /*TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD); + Int_t nGoodTracks = list->GetEntries(); + // loop over esd tracks + for (Int_t i=0; i (list->At(i)); + if (!esdTrack) + { + AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i)); + continue; + } + + Double_t p[3]; + esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist + TVector3 vector(p); + + Float_t theta = vector.Theta(); + Float_t eta = -TMath::Log(TMath::Tan(theta/2.)); + Float_t pt = vector.Pt(); + + if (pt < kPtCut) + continue; + + if (TMath::Abs(eta) < 0.5) + nESDTracks05++; - fCorrelation->Fill(nMCTracks, nESDTracks); + if (TMath::Abs(eta) < 1.0) + nESDTracks10++; + + if (TMath::Abs(eta) < 1.5) + nESDTracks15++; + + if (TMath::Abs(eta) < 2.0) + nESDTracks20++; + }*/ + + fMultiplicity->FillMeasured(vtxMC[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); + + // TODO should this be vtx[2] or vtxMC[2] ? + fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20); return kTRUE; } @@ -186,9 +312,7 @@ void AliMultiplicityMCSelector::SlaveTerminate() return; } - fOutput->Add(fMultiplicityESD); - fOutput->Add(fMultiplicityMC); - fOutput->Add(fCorrelation); + fOutput->Add(fMultiplicity); } void AliMultiplicityMCSelector::Terminate() @@ -199,23 +323,19 @@ void AliMultiplicityMCSelector::Terminate() AliSelector::Terminate(); - fOutput->Print(); - - fMultiplicityESD = dynamic_cast (fOutput->FindObject("fMultiplicityESD")); - fMultiplicityMC = dynamic_cast (fOutput->FindObject("fMultiplicityMC")); - fCorrelation = dynamic_cast (fOutput->FindObject("fCorrelation")); + fMultiplicity = dynamic_cast (fOutput->FindObject("Multiplicity")); - if (!fMultiplicityESD || !fMultiplicityMC || !fCorrelation) + if (!fMultiplicity) { - AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p %p", (void*) fMultiplicityESD, (void*) fMultiplicityMC, (void*) fCorrelation)); + AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity)); return; } TFile* file = TFile::Open("multiplicityMC.root", "RECREATE"); - fMultiplicityMC->Write(); - fMultiplicityESD->Write(); - fCorrelation->Write(); + fMultiplicity->SaveHistograms(); file->Close(); + + fMultiplicity->DrawHistograms(); } diff --git a/PWG0/dNdEta/AliMultiplicityMCSelector.h b/PWG0/dNdEta/AliMultiplicityMCSelector.h index 1e139150041..675643b8cae 100644 --- a/PWG0/dNdEta/AliMultiplicityMCSelector.h +++ b/PWG0/dNdEta/AliMultiplicityMCSelector.h @@ -6,8 +6,9 @@ #include "AliSelectorRL.h" class AliESDtrackCuts; -class TH1F; class TH2F; +class TH3F; +class AliMultiplicityCorrection; class AliMultiplicityMCSelector : public AliSelectorRL { public: @@ -16,6 +17,7 @@ class AliMultiplicityMCSelector : public AliSelectorRL { virtual void Begin(TTree* tree); virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); virtual Bool_t Process(Long64_t entry); virtual void SlaveTerminate(); virtual void Terminate(); @@ -23,12 +25,8 @@ class AliMultiplicityMCSelector : public AliSelectorRL { protected: void ReadUserObjects(TTree* tree); - TH1F* fMultiplicityESD; // multiplicity histogram - TH1F* fMultiplicityMC; // multiplicity histogram - - TH2F* fCorrelation; // (gene multiplicity) vs (meas multiplicity) - - AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts + AliMultiplicityCorrection* fMultiplicity; // object containing the extracted data + AliESDtrackCuts* fEsdTrackCuts; // Object containing the parameters of the esd track cuts private: AliMultiplicityMCSelector(const AliMultiplicityMCSelector&); diff --git a/PWG0/dNdEta/runMultiplicitySelector.C b/PWG0/dNdEta/runMultiplicitySelector.C index 8f376d0da62..dc8eaf7370c 100644 --- a/PWG0/dNdEta/runMultiplicitySelector.C +++ b/PWG0/dNdEta/runMultiplicitySelector.C @@ -55,11 +55,31 @@ void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_ printf("ERROR: Executing process failed with %d.\n", result); return; } +} - // and draw it - if (aMC != kFALSE) - MultiplicityMC(); - else - MultiplicityESD(); +void draw(const char* fileName = "multiplicityMC.root") +{ + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileName); + mult->LoadHistograms("Multiplicity"); + + mult->DrawHistograms(); } + +void* fit(const char* fileName = "multiplicityMC.root") +{ + gSystem->Load("libPWG0base"); + + AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity"); + + TFile::Open(fileName); + mult->LoadHistograms("Multiplicity"); + + mult->ApplyMinuitFit(3, kFALSE); + + return mult; +} diff --git a/PWG0/libPWG0base.pkg b/PWG0/libPWG0base.pkg index bec2d192470..d78431bd889 100644 --- a/PWG0/libPWG0base.pkg +++ b/PWG0/libPWG0base.pkg @@ -7,7 +7,8 @@ HDRS = dNdEta/dNdEtaAnalysis.h \ AliCorrectionMatrix2D.h \ AliCorrectionMatrix3D.h \ AliCorrection.h \ - AliPWG0Helper.h + AliPWG0Helper.h \ + dNdEta/AliMultiplicityCorrection.h SRCS = $(HDRS:.h=.cxx) -- 2.39.3