From: rbertens Date: Thu, 27 Mar 2014 16:48:37 +0000 (+0100) Subject: implemented weighted merging of centrality bins X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=c03f7598a197378002de2298039346fc67f053db;p=u%2Fmrichter%2FAliRoot.git implemented weighted merging of centrality bins --- diff --git a/PWG/FLOW/Tasks/AliJetFlowTools.cxx b/PWG/FLOW/Tasks/AliJetFlowTools.cxx index c96e3804463..7b5a0546818 100644 --- a/PWG/FLOW/Tasks/AliJetFlowTools.cxx +++ b/PWG/FLOW/Tasks/AliJetFlowTools.cxx @@ -81,8 +81,8 @@ AliJetFlowTools::AliJetFlowTools() : fRefreshInput (kTRUE), fOutputFileName ("UnfoldedSpectra.root"), fOutputFile (0x0), - fCentralityBin (0), fCentralityArray (0x0), + fCentralityWeights (0x0), fDetectorResponse (0x0), fJetFindingEff (0x0), fBetaIn (.1), @@ -366,7 +366,7 @@ TH1D* AliJetFlowTools::UnfoldWrapper( else if(fUnfoldingAlgorithm == kNone) { TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone")); clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data())); - return clone; + return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE); } else return 0x0; // do the actual unfolding with the selected function @@ -865,23 +865,27 @@ Bool_t AliJetFlowTools::PrepareForUnfolding() // clear minuit state to avoid constraining the fit with the results of the previous iteration for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.); } - // extract the spectra - TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin)); - fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data())); - if(!fJetPtDeltaPhi) { + // extract the spectra + TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0))); + if(!fInputList->FindObject(spectrumName.Data())) { printf(" Couldn't find spectrum %s ! \n", spectrumName.Data()); return kFALSE; } - if(fCentralityArray) { - printf(" merging centralities \n"); - for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { - spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i)); - printf( " searching for %s \n", spectrumName.Data()); - fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data()))); - } - } + // get the first scaled spectrum + fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data()); + // clone the spectrum on the heap. this is necessary since scale or add change the + // contents of the original histogram fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE); + fJetPtDeltaPhi->Scale(fCentralityWeights->At(0)); + printf("Extracted %i wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0)); + // merge subsequent bins (if any) + for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { + spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i)); + printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i)); + fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i)); + } + // in plane spectrum if(!fDphiUnfolding) { fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e"); @@ -896,12 +900,10 @@ Bool_t AliJetFlowTools::PrepareForUnfolding() } // normalize spectra to event count if requested if(fNormalizeSpectra) { - TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin))); - if(fCentralityArray) { - for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { - printf(" merging centralities \n"); - rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i)))); - } + TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0)))); + rho->Scale(fCentralityWeights->At(0)); + for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { + rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i)); } if(!rho) return 0x0; Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE; @@ -932,7 +934,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding() } // extract the delta pt matrices TString deltaptName(""); - deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin); + deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0)); fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data())); if(!fDeltaPtDeltaPhi) { printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data()); @@ -940,16 +942,17 @@ Bool_t AliJetFlowTools::PrepareForUnfolding() fRefreshInput = kTRUE; return kTRUE; } - if(fCentralityArray) { - printf(" merging centralities \n"); - for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { - deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i)); - printf(" searching for %s \n ", deltaptName.Data()); - fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data()))); - } - } + // clone the distribution on the heap and if requested merge with other centralities fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE); + fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0)); + printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0)); + for(Int_t i(1); i < fCentralityArray->GetSize(); i++) { + deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i)); + printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i)); + fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i)); + } + // in plane delta pt distribution if(!fDphiUnfolding) { fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e"); @@ -982,12 +985,12 @@ Bool_t AliJetFlowTools::PrepareForUnfolding() } } fDptIn = new TH2D(*rfIn); - fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin)); + fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0))); fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]"); fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]"); fDptIn = ProtectHeap(fDptIn); fDptOut = new TH2D(*rfOut); - fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin)); + fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0))); fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]"); fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]"); fDptOut = ProtectHeap(fDptOut); @@ -1014,7 +1017,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) { for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.); } // extract the spectra - TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin)); + TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0))); fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data())); if(!fJetPtDeltaPhi) { printf(" Couldn't find spectrum %s ! \n", spectrumName.Data()); @@ -1031,7 +1034,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) { fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e"); // extract the delta pt matrices TString deltaptName(""); - deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin); + deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0)); fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data())); if(!fDeltaPtDeltaPhi) { printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data()); @@ -1059,7 +1062,7 @@ Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) { } } fDptIn = new TH2D(*rfIn); - fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin)); + fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0))); fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]"); fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]"); fDptIn = ProtectHeap(fDptIn); @@ -1400,7 +1403,7 @@ void AliJetFlowTools::Style(TH1* h, EColor col, histoType type) h->SetLineColor(col); h->SetMarkerColor(col); h->SetLineWidth(2.); - h->SetMarkerSize(2.); + h->SetMarkerSize(1.); h->SetTitle(""); h->GetYaxis()->SetLabelSize(0.05); h->GetXaxis()->SetLabelSize(0.05); @@ -1408,31 +1411,36 @@ void AliJetFlowTools::Style(TH1* h, EColor col, histoType type) h->GetXaxis()->SetTitleOffset(1.5); h->GetYaxis()->SetTitleSize(.05); h->GetXaxis()->SetTitleSize(.05); - h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); switch (type) { case kInPlaneSpectrum : { h->SetTitle("IN PLANE"); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}"); } break; case kOutPlaneSpectrum : { h->SetTitle("OUT OF PLANE"); h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}"); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); } break; case kUnfoldedSpectrum : { h->SetTitle("UNFOLDED"); h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}"); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); } break; case kFoldedSpectrum : { h->SetTitle("FOLDED"); h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}"); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); } break; case kMeasuredSpectrum : { h->SetTitle("MEASURED"); h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}"); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); } break; case kBar : { h->SetFillColor(col); h->SetBarWidth(.6); + h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} [GeV/#it{c}]"); h->SetBarOffset(0.2); } default : break; @@ -1445,7 +1453,7 @@ void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type) h->SetLineColor(col); h->SetMarkerColor(col); h->SetLineWidth(2.); - h->SetMarkerSize(2.); + h->SetMarkerSize(1.); h->SetTitle(""); h->SetFillColor(kYellow); h->GetYaxis()->SetLabelSize(0.05); @@ -1505,7 +1513,6 @@ void AliJetFlowTools::GetNominalValues( printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <"); readMe->ls(); TFile* output(new TFile(outFile.Data(), "RECREATE")); // create new output file - // get some placeholders, have to be initialized but will be deleted ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()); TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray())); @@ -1532,7 +1539,6 @@ void AliJetFlowTools::GetNominalValues( // close the open files, reclaim ownership of histograms which are necessary outside of the file scope ratio->SetDirectory(0); // disassociate from current gDirectory readMe->Close(); -// output->Close(); } //_____________________________________________________________________________ void AliJetFlowTools::GetCorrelatedUncertainty( @@ -1706,8 +1712,8 @@ void AliJetFlowTools::GetCorrelatedUncertainty( upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1); } else { // the in and out of plane correlated errors will be fully correlated, so take the correlation coefficient equal to 1 - lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) - 2.*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1); - upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 2.*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1); + lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) - 1.*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1); + upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 1.*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1); } // set the errors ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1); @@ -1748,7 +1754,7 @@ void AliJetFlowTools::GetCorrelatedUncertainty( relativeErrorInLow, relativeErrorOutUp, relativeErrorOutLow, - 1.)); + .5)); // pass the nominal values to the pointer references corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone(); TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}")); @@ -3164,11 +3170,11 @@ Bool_t AliJetFlowTools::SetRawInput ( fSpectrumOut->Scale(1./((double)fEventCount)); } fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError); - fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin)); + fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0))); fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]"); fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]"); fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError); - fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin)); + fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0))); fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]"); fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]"); @@ -3251,20 +3257,15 @@ TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name) Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.); Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.); - cout << " GetV2 " << endl; - cout << " prefactor " << pre << endl; for(Int_t i(1); i <= h1->GetNbinsX(); i++) { binCent = h1->GetXaxis()->GetBinCenter(i); binWidth = h1->GetXaxis()->GetBinWidth(i); if(h2->GetBinContent(i) > 0.) { in = h1->GetBinContent(i); - cout << " yield in " << in << endl; ein = h1->GetBinError(i); out = h2->GetBinContent(i); - cout << " yield out " << out << endl; eout = h2->GetBinError(i); ratio = pre*((in-out)/(in+out)); - cout << " v2 " << ratio << endl; error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout; error2 = error2*pre*pre; if(error2 > 0) error2 = TMath::Sqrt(error2); @@ -3394,8 +3395,8 @@ void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) summary->GetXaxis()->SetBinLabel(1, "fBetaIn"); summary->SetBinContent(2, fBetaOut); summary->GetXaxis()->SetBinLabel(2, "fBetaOut"); - summary->SetBinContent(3, fCentralityBin); - summary->GetXaxis()->SetBinLabel(3, "fCentralityBin"); + summary->SetBinContent(3, fCentralityArray->At(0)); + summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]"); summary->SetBinContent(4, (int)convergedIn); summary->GetXaxis()->SetBinLabel(4, "convergedIn"); summary->SetBinContent(5, (int)convergedOut); diff --git a/PWG/FLOW/Tasks/AliJetFlowTools.h b/PWG/FLOW/Tasks/AliJetFlowTools.h index a85f4cc4c9a..30596dd52a4 100644 --- a/PWG/FLOW/Tasks/AliJetFlowTools.h +++ b/PWG/FLOW/Tasks/AliJetFlowTools.h @@ -78,10 +78,20 @@ class AliJetFlowTools { fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data()); fActiveDir->cd(); } - void SetCentralityBin(Int_t bin) {fCentralityBin = bin;} + void SetCentralityBin(Int_t bin) { + // in case of one centraltiy + fCentralityArray = new TArrayI(1); + fCentralityArray->AddAt(bin, 0); + // for one centrality there's no need for weights + fCentralityWeights = new TArrayD(1); + fCentralityWeights->AddAt(1., 0); + } void SetCentralityBin(TArrayI* bins) { fCentralityArray = bins; - fCentralityBin = fCentralityArray->At(0); + } + void SetCentralityWeight(TArrayD* weights) { + fCentralityWeights = weights; + if(!fCentralityArray) printf(" > Warning: centrality weights set, but bins are not defined! \n"); } void SetDetectorResponse(TH2D* dr) {fDetectorResponse = dr;} void SetJetFindingEfficiency(TH1D* e) {fJetFindingEff = e;} @@ -312,8 +322,8 @@ class AliJetFlowTools { Bool_t fRefreshInput; // re-read the input (called automatically if input list changes) TString fOutputFileName; // output file name TFile* fOutputFile; // output file - Int_t fCentralityBin; // centrality bin TArrayI* fCentralityArray; // array of bins that are merged + TArrayD* fCentralityWeights; // array of centrality weights TH2D* fDetectorResponse; // detector response TH1D* fJetFindingEff; // jet finding efficiency Double_t fBetaIn; // regularization strength, in plane unfolding