X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGLF%2FFORWARD%2Fanalysis2%2FAliBasedNdetaTask.cxx;h=fa2337e378d69bcaf3c6601a743c6189f5a83664;hb=959254d6c1a60993de84fc881a86fbb8a7e6715e;hp=2407298b8132018b95c53d9fd96031783c92afc7;hpb=09d5920fcfc6bd42b5941c9dbeb76e00be842cd5;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx index 2407298b813..fa2337e378d 100644 --- a/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliBasedNdetaTask.cxx @@ -13,32 +13,25 @@ #include "AliAODForwardMult.h" #include #include +#include //____________________________________________________________________ AliBasedNdetaTask::AliBasedNdetaTask() - : AliAnalysisTaskSE(), - fSums(0), // Container of sums - fOutput(0), // Container of output - fVtxMin(0), // Minimum v_z - fVtxMax(0), // Maximum v_z - fTriggerMask(0), // Trigger mask + : AliBaseAODTask(), fRebin(0), // Rebinning factor fCutEdges(false), fSymmetrice(true), fCorrEmpty(true), fUseROOTProj(false), fTriggerEff(1), - fTriggerEff0(1), + fTriggerEff0(1), fShapeCorr(0), fListOfCentralities(0), - fSNNString(0), - fSysString(0), - fCent(0), - fCentAxis(0), fNormalizationScheme(kFull), - fSchemeString(0), - fTriggerString(0), - fFinalMCCorrFile("") + fFinalMCCorrFile(""), + fSatelliteVertices(0), + fglobalempiricalcorrection(0), + fmeabsignalvscentr(0) { // // Constructor @@ -48,12 +41,7 @@ AliBasedNdetaTask::AliBasedNdetaTask() //____________________________________________________________________ AliBasedNdetaTask::AliBasedNdetaTask(const char* name) - : AliAnalysisTaskSE(name), - fSums(0), // Container of sums - fOutput(0), // Container of output - fVtxMin(-10), // Minimum v_z - fVtxMax(10), // Maximum v_z - fTriggerMask(AliAODForwardMult::kInel), + : AliBaseAODTask(Form("%sdNdeta", name)), fRebin(5), // Rebinning factor fCutEdges(false), fSymmetrice(true), @@ -63,60 +51,27 @@ AliBasedNdetaTask::AliBasedNdetaTask(const char* name) fTriggerEff0(1), fShapeCorr(0), fListOfCentralities(0), - fSNNString(0), - fSysString(0), - fCent(0), - fCentAxis(0), fNormalizationScheme(kFull), - fSchemeString(0), - fTriggerString(0), - fFinalMCCorrFile("") + fFinalMCCorrFile(""), + fSatelliteVertices(0), + fglobalempiricalcorrection(0), + fmeabsignalvscentr(0) { // // Constructor // DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name); + + fTriggerMask = AliAODForwardMult::kInel; + fMinIpZ = -10; + fMaxIpZ = +10; fListOfCentralities = new TObjArray(1); // Set the normalisation scheme SetNormalizationScheme(kFull); - // Set the trigger mask - SetTriggerMask(AliAODForwardMult::kInel); - - // Output slot #1 writes into a TH1 container - DefineOutput(1, TList::Class()); - DefineOutput(2, TList::Class()); } -//____________________________________________________________________ -AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o) - : AliAnalysisTaskSE(o), - fSums(o.fSums), // TList* - Container of sums - fOutput(o.fOutput), // Container of output - fVtxMin(o.fVtxMin), // Double_t - Minimum v_z - fVtxMax(o.fVtxMax), // Double_t - Maximum v_z - fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask - fRebin(o.fRebin), // Int_t - Rebinning factor - fCutEdges(o.fCutEdges), // Bool_t - Whether to cut edges when rebinning - fSymmetrice(o.fSymmetrice), - fCorrEmpty(o.fCorrEmpty), - fUseROOTProj(o.fUseROOTProj), - fTriggerEff(o.fTriggerEff), - fTriggerEff0(o.fTriggerEff0), - fShapeCorr(o.fShapeCorr), - fListOfCentralities(o.fListOfCentralities), - fSNNString(o.fSNNString), - fSysString(o.fSysString), - fCent(o.fCent), - fCentAxis(o.fCentAxis), - fNormalizationScheme(o.fNormalizationScheme), - fSchemeString(o.fSchemeString), - fTriggerString(o.fTriggerString), - fFinalMCCorrFile(o.fFinalMCCorrFile) -{ - DGUARD(fDebug, 3,"Copy CTOR of AliBasedNdetaTask"); -} //____________________________________________________________________ AliBasedNdetaTask::~AliBasedNdetaTask() @@ -125,16 +80,6 @@ AliBasedNdetaTask::~AliBasedNdetaTask() // Destructor // DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask"); - if (fSums) { - fSums->Delete(); - delete fSums; - fSums = 0; - } - if (fOutput) { - fOutput->Delete(); - delete fOutput; - fOutput = 0; - } } //________________________________________________________________________ @@ -149,22 +94,6 @@ AliBasedNdetaTask::SetDebugLevel(Int_t lvl) } } -//________________________________________________________________________ -void -AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins) -{ - DGUARD(fDebug,3,"Set centrality axis, %d bins", n); - if (!fCentAxis) { - fCentAxis = new TAxis(); - fCentAxis->SetName("centAxis"); - fCentAxis->SetTitle("Centrality [%]"); - } - TArrayD dbins(n+1); - for (UShort_t i = 0; i <= n; i++) - dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]); - fCentAxis->Set(n, dbins.GetArray()); -} - //________________________________________________________________________ void AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high) @@ -184,6 +113,7 @@ AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high) GetName(), low, high, at); return; } + bin->SetSatelliteVertices(fSatelliteVertices); bin->SetDebugLevel(fDebug); fListOfCentralities->AddAtAndExpand(bin, at); } @@ -284,32 +214,7 @@ AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme) { DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme); fNormalizationScheme = scheme; - if (fSchemeString) delete fSchemeString; - fSchemeString = AliForwardUtil::MakeParameter("scheme", scheme); -} -//________________________________________________________________________ -void -AliBasedNdetaTask::SetTriggerMask(const char* mask) -{ - // - // Set the trigger maskl - // - // Parameters: - // mask Trigger mask - // - DGUARD(fDebug,3,"Set the trigger mask: %s", mask); - SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask)); } -//________________________________________________________________________ -void -AliBasedNdetaTask::SetTriggerMask(UShort_t mask) -{ - DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask); - fTriggerMask = mask; - if (fTriggerString) delete fTriggerString; - fTriggerString = AliForwardUtil::MakeParameter("trigger", fTriggerMask); -} - //________________________________________________________________________ void AliBasedNdetaTask::SetShapeCorrection(const TH2F* c) @@ -326,7 +231,6 @@ AliBasedNdetaTask::SetShapeCorrection(const TH2F* c) fShapeCorr = static_cast(c->Clone()); fShapeCorr->SetDirectory(0); } - //________________________________________________________________________ void AliBasedNdetaTask::InitializeCentBins() @@ -335,17 +239,17 @@ AliBasedNdetaTask::InitializeCentBins() // Automatically add 'all' centrality bin if nothing has been defined. AddCentralityBin(0, 0, 0); - if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) { - const TArrayD* bins = fCentAxis->GetXbins(); - Int_t nbin = fCentAxis->GetNbins(); + if (HasCentrality() && fCentAxis.GetXbins()) { + const TArrayD* bins = fCentAxis.GetXbins(); + Int_t nbin = fCentAxis.GetNbins(); for (Int_t i = 0; i < nbin; i++) AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1])); } } //________________________________________________________________________ -void -AliBasedNdetaTask::UserCreateOutputObjects() +Bool_t +AliBasedNdetaTask::Book() { // // Create output objects. @@ -353,23 +257,13 @@ AliBasedNdetaTask::UserCreateOutputObjects() // This is called once per slave process // DGUARD(fDebug,1,"Create user ouput object"); - fSums = new TList; - fSums->SetName(Form("%s_sums", GetName())); - fSums->SetOwner(); - InitializeCentBins(); - if (fCentAxis) fSums->Add(fCentAxis); - - fSums->Add(AliForwardUtil::MakeParameter("alirootRev", - AliForwardUtil::AliROOTRevision())); - fSums->Add(AliForwardUtil::MakeParameter("alirootBranch", - AliForwardUtil::AliROOTBranch())); + fSums->Add(AliForwardUtil::MakeParameter("empirical", + fglobalempiricalcorrection != 0)); + fSums->Add(AliForwardUtil::MakeParameter("scheme", fNormalizationScheme)); - // Centrality histogram - fCent = new TH1D("cent", "Centrality", 100, 0, 100); - fCent->SetDirectory(0); - fCent->SetXTitle(0); - fSums->Add(fCent); + // Make our centrality bins + InitializeCentBins(); // Loop over centrality bins TIter next(fListOfCentralities); @@ -377,15 +271,28 @@ AliBasedNdetaTask::UserCreateOutputObjects() while ((bin = static_cast(next()))) bin->CreateOutputObjects(fSums, fTriggerMask); + fmeabsignalvscentr=new TH2D("meanAbsSignalVsCentr", + "Mean absolute signal versus centrality", + 400, 0, 20, 100, 0, 100); + fSums->Add(fmeabsignalvscentr); - // Post data for ALL output slots >0 here, to get at least an empty - // histogram - PostData(1, fSums); + return true; } - + //____________________________________________________________________ -void -AliBasedNdetaTask::UserExec(Option_t *) +Bool_t +AliBasedNdetaTask::CheckEvent(const AliAODForwardMult& fwd) +{ + AliBaseAODTask::CheckEvent(fwd); + // Here, we always return true, as the centrality bins will do + // their own checks on the events - this is needed for event + // normalization etc. + return true; +} + +//____________________________________________________________________ +Bool_t +AliBasedNdetaTask::Event(AliAODEvent& aod) { // // Process a single event @@ -395,59 +302,53 @@ AliBasedNdetaTask::UserExec(Option_t *) // // Main loop DGUARD(fDebug,1,"Analyse the AOD event"); - AliAODEvent* aod = AliForwardUtil::GetAODEvent(this); - if (!aod) return; - - TObject* obj = aod->FindListObject("Forward"); - if (!obj) { - AliWarning("No forward object found"); - return; - } - AliAODForwardMult* forward = static_cast(obj); + AliAODForwardMult* forward = GetForward(aod); + if (!forward) return false; // Fill centrality histogram - Float_t cent = forward->GetCentrality(); - fCent->Fill(cent); + + Double_t vtx = forward->GetIpZ(); + TH2D* data = GetHistogram(aod, false); + TH2D* dataMC = GetHistogram(aod, true); + if (!data) return false; + + CheckEventData(vtx, data, dataMC); + + if (!ApplyEmpiricalCorrection(forward,data)) + return false; - // Get the histogram(s) - TH2D* data = GetHistogram(aod, false); - TH2D* dataMC = GetHistogram(aod, true); Bool_t isZero = ((fNormalizationScheme & kZeroBin) && !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0)); - - + Bool_t taken = false; + // Loop over centrality bins CentralityBin* allBin = static_cast(fListOfCentralities->At(0)); - allBin->ProcessEvent(forward, fTriggerMask, isZero, - fVtxMin, fVtxMax, data, dataMC); + if (allBin->ProcessEvent(forward, fTriggerMask, isZero, + fMinIpZ, fMaxIpZ, data, dataMC)) taken = true; // Find this centrality bin - if (fCentAxis && fCentAxis->GetNbins() > 0) { - Int_t icent = fCentAxis->FindBin(cent); + if (HasCentrality()) { + Double_t cent = forward->GetCentrality(); + Int_t icent = fCentAxis.FindBin(cent); CentralityBin* thisBin = 0; - if (icent >= 1 && icent <= fCentAxis->GetNbins()) + if (icent >= 1 && icent <= fCentAxis.GetNbins()) thisBin = static_cast(fListOfCentralities->At(icent)); if (thisBin) - thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin, fVtxMax, - data, dataMC); + if (thisBin->ProcessEvent(forward, fTriggerMask, isZero, fMinIpZ, + fMaxIpZ, data, dataMC)) taken = true; } + + return taken; +} - // Here, we get the update - if (!fSNNString) { - fSNNString = AliForwardUtil::MakeParameter("sNN", forward->GetSNN()); - fSysString = AliForwardUtil::MakeParameter("sys", forward->GetSystem()); - - fSums->Add(fSNNString); - fSums->Add(fSysString); - fSums->Add(fSchemeString); - fSums->Add(fTriggerString); - - // Print(); - } - PostData(1, fSums); +//________________________________________________________________________ +void AliBasedNdetaTask::CheckEventData(Double_t, + TH2*, + TH2*) +{ } //________________________________________________________________________ @@ -470,7 +371,10 @@ AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker, h->SetMarkerStyle(marker); h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1); h->SetFillStyle(0); - h->SetYTitle(ytitle); + TString ytit; + if (ytitle && ytitle[0] != '\0') ytit = ytitle; + ytit = "#frac{1}{#it{N}}#frac{d#it{N}_{ch}}{d#it{#eta}}"; + h->SetYTitle(ytit); h->GetXaxis()->SetTitleFont(132); h->GetXaxis()->SetLabelFont(132); h->GetXaxis()->SetNdivisions(10); @@ -572,7 +476,7 @@ AliBasedNdetaTask::ProjectX(const TH2D* h, } // Loop over X bins - // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName())); + //DMSG(fDebug,3,"Projecting bins [%d,%d] of %s", first, last, h->GetName())); Int_t ybins = (last-first+1); for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) { Double_t content = 0; @@ -581,8 +485,8 @@ AliBasedNdetaTask::ProjectX(const TH2D* h, for (Int_t ybin = first; ybin <= last; ybin++) { - Double_t c1 = h->GetCellContent(xbin, ybin); - Double_t e1 = h->GetCellError(xbin, ybin); + Double_t c1 = h->GetBinContent(h->GetBin(xbin, ybin)); + Double_t e1 = h->GetBinError(h->GetBin(xbin, ybin)); // Ignore empty bins if (c1 < 1e-12) continue; @@ -614,10 +518,10 @@ AliBasedNdetaTask::ProjectX(const TH2D* h, return ret; } - + //________________________________________________________________________ -void -AliBasedNdetaTask::Terminate(Option_t *) +Bool_t +AliBasedNdetaTask::Finalize() { // // Called at end of event processing.. @@ -631,26 +535,20 @@ AliBasedNdetaTask::Terminate(Option_t *) // once at the end of the query DGUARD(fDebug,1,"Process final merged results"); - fSums = dynamic_cast (GetOutputData(1)); - if(!fSums) { - AliError("Could not retrieve TList fSums"); - return; - } - - fOutput = new TList; - fOutput->SetName(Form("%s_result", GetName())); - fOutput->SetOwner(); - - fSNNString = fSums->FindObject("sNN"); - fSysString = fSums->FindObject("sys"); - fCentAxis = static_cast(fSums->FindObject("centAxis")); - fSchemeString = fSums->FindObject("scheme"); - fTriggerString = fSums->FindObject("trigger"); + UShort_t sNN; + UShort_t sys; + ULong_t trig; + UShort_t scheme; + AliForwardUtil::GetParameter(fSums->FindObject("sNN"), sNN); + AliForwardUtil::GetParameter(fSums->FindObject("sys"), sys); + AliForwardUtil::GetParameter(fSums->FindObject("trigger"), trig); + AliForwardUtil::GetParameter(fSums->FindObject("scheme"), scheme); + + TAxis* cA = static_cast(fSums->FindObject("centAxis")); + if (cA) cA->Copy(fCentAxis); - if(fSysString && fSNNString && - fSysString->GetUniqueID() == AliForwardUtil::kPP) - LoadNormalizationData(fSysString->GetUniqueID(), - fSNNString->GetUniqueID()); + if(sNN > 0 && sys == AliForwardUtil::kPP) + LoadNormalizationData(sys, sNN); InitializeCentBins(); @@ -683,27 +581,31 @@ AliBasedNdetaTask::Terminate(Option_t *) Int_t style = GetMarker(); Int_t color = GetColor(); - AliInfo(Form("Marker style=%d, color=%d", style, color)); + DMSG(fDebug,3,"Marker style=%d, color=%d", style, color); while ((bin = static_cast(next()))) { - - bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, + bin->End(fSums, fResults, fNormalizationScheme, fShapeCorr, fTriggerEff, fTriggerEff0, fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, fTriggerMask, style, color, mclist, truthlist); - if (fCentAxis && bin->IsAllBin()) continue; - TH1* dndeta = bin->GetResult(0, false, ""); - TH1* dndetaSym = bin->GetResult(0, true, ""); - TH1* dndetaMC = bin->GetResult(0, false, "MC"); - TH1* dndetaMCSym = bin->GetResult(0, true, "MC"); + if (HasCentrality() && bin->IsAllBin()) + // If we have centrality bins, do not add the min-bias + // distribution to the output stack. + continue; + TH1* dndeta = bin->GetResult(0, false, ""); + TH1* dndetaSym = fSymmetrice ? bin->GetResult(0, true, "") : 0; + TH1* dndetaMC = bin->GetResult(0, false, "MC", false); + TH1* dndetaMCSym = fSymmetrice ? bin->GetResult(0, true, "MC", false) : 0; + DMSG(fDebug,2,"Results: bare=%p sym=%p mcbare=%p mcsym=%p", + dndeta, dndetaSym, dndetaMC, dndetaMCSym); if (dndeta) dndetaStack->Add(dndeta); if (dndetaSym) dndetaStack->Add(dndetaSym); if (dndetaMC) dndetaMCStack->Add(dndetaMC); if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym); if (fRebin > 1) { - dndeta = bin->GetResult(fRebin, false, ""); - dndetaSym = bin->GetResult(fRebin, true, ""); - dndetaMC = bin->GetResult(fRebin, false, "MC"); - dndetaMCSym = bin->GetResult(fRebin, true, "MC"); + dndeta = bin->GetResult(fRebin, false, ""); + dndetaSym = fSymmetrice ? bin->GetResult(fRebin, true, "") : 0; + dndetaMC = bin->GetResult(fRebin, false, "MC", false); + dndetaMCSym = fSymmetrice ? bin->GetResult(fRebin, true, "MC", false): 0; if (dndeta) dndetaStackRebin->Add(dndeta); if (dndetaSym) dndetaStackRebin->Add(dndetaSym); if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC); @@ -711,84 +613,81 @@ AliBasedNdetaTask::Terminate(Option_t *) } } // Output the stack - fOutput->Add(dndetaStack); + fResults->Add(dndetaStack); // If available output rebinned stack - if (!dndetaStackRebin->GetHists() || - dndetaStackRebin->GetHists()->GetEntries() <= 0) { + if (fRebin > 0 && (!dndetaStackRebin->GetHists() || + dndetaStackRebin->GetHists()->GetEntries() <= 0)) { AliWarning("No rebinned histograms found"); delete dndetaStackRebin; dndetaStackRebin = 0; } - if (dndetaStackRebin) fOutput->Add(dndetaStackRebin); + if (dndetaStackRebin) fResults->Add(dndetaStackRebin); // If available, output track-ref stack if (!dndetaMCStack->GetHists() || dndetaMCStack->GetHists()->GetEntries() <= 0) { - AliWarning("No MC histograms found"); + // AliWarning("No MC histograms found"); delete dndetaMCStack; dndetaMCStack = 0; } - if (dndetaMCStack) fOutput->Add(dndetaMCStack); + if (dndetaMCStack) fResults->Add(dndetaMCStack); // If available, output rebinned track-ref stack - if (!dndetaMCStackRebin->GetHists() || - dndetaMCStackRebin->GetHists()->GetEntries() <= 0) { + if ((fRebin > 0 && dndetaMCStack) + && (!dndetaMCStackRebin->GetHists() || + dndetaMCStackRebin->GetHists()->GetEntries() <= 0)) { AliWarning("No rebinned MC histograms found"); delete dndetaMCStackRebin; dndetaMCStackRebin = 0; } - if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin); + if (dndetaMCStackRebin) fResults->Add(dndetaMCStackRebin); // Output collision energy string - if (fSNNString) { - UShort_t sNN = fSNNString->GetUniqueID(); - TNamed* sNNObj = new TNamed(fSNNString->GetName(), + if (sNN > 0) { + TNamed* sNNObj = new TNamed("sNN", AliForwardUtil::CenterOfMassEnergyString(sNN)); sNNObj->SetUniqueID(sNN); - fOutput->Add(sNNObj); // fSNNString->Clone()); + fResults->Add(sNNObj); // fSNNString->Clone()); } // Output collision system string - if (fSysString) { - UShort_t sys = fSysString->GetUniqueID(); - TNamed* sysObj = new TNamed(fSysString->GetName(), + if (sys > 0) { + TNamed* sysObj = new TNamed("sys", AliForwardUtil::CollisionSystemString(sys)); sysObj->SetUniqueID(sys); - fOutput->Add(sysObj); // fSysString->Clone()); + fResults->Add(sysObj); // fSysString->Clone()); } // Output centrality axis - if (fCentAxis) fOutput->Add(fCentAxis); + fResults->Add(&fCentAxis); // Output trigger string - if (fTriggerString) { - UShort_t mask = fTriggerString->GetUniqueID(); - TNamed* maskObj = new TNamed(fTriggerString->GetName(), - AliAODForwardMult::GetTriggerString(mask)); - maskObj->SetUniqueID(mask); - fOutput->Add(maskObj); // fTriggerString->Clone()); + if (trig) { + TNamed* maskObj = new TNamed("trigger", + AliAODForwardMult::GetTriggerString(trig)); + maskObj->SetUniqueID(trig); + fResults->Add(maskObj); // fTriggerString->Clone()); } // Normalization string - if (fSchemeString) { - UShort_t scheme = fSchemeString->GetUniqueID(); - TNamed* schemeObj = new TNamed(fSchemeString->GetName(), + if (scheme > 0) { + TNamed* schemeObj = new TNamed("scheme", NormalizationSchemeString(scheme)); schemeObj->SetUniqueID(scheme); - fOutput->Add(schemeObj); // fSchemeString->Clone()); + fResults->Add(schemeObj); // fSchemeString->Clone()); } // Output vertex axis - TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax); + TAxis* vtxAxis = new TAxis(1,fMinIpZ,fMaxIpZ); vtxAxis->SetName("vtxAxis"); - vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax)); - fOutput->Add(vtxAxis); + vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fMinIpZ,fMaxIpZ)); + fResults->Add(vtxAxis); - // Output trigger efficiency and shape correction - fOutput->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff)); - fOutput->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0)); - if (fShapeCorr) fOutput->Add(fShapeCorr); + // Output trigger efficiency and shape correction + fResults->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff)); + fResults->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0)); + if (fShapeCorr) fResults->Add(fShapeCorr); TNamed* options = new TNamed("options",""); TString str; @@ -796,9 +695,9 @@ AliBasedNdetaTask::Terminate(Option_t *) str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not ")); str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not ")); options->SetTitle(str); - fOutput->Add(options); + fResults->Add(options); - PostData(2, fOutput); + return true; } //________________________________________________________________________ void @@ -812,21 +711,32 @@ AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy) if(energy == 7000) snn.Form("7000"); if(energy == 2750) snn.Form("2750"); - if(fShapeCorr && (fTriggerEff != 1)) { - AliInfo("Objects already set for normalization - no action taken"); + // Check if shape correction/trigger efficiency was requsted and not + // already set + Bool_t needShape = ((fNormalizationScheme & kShape) && !fShapeCorr); + Bool_t needEff = ((fNormalizationScheme & kTriggerEfficiency) && + ((1 - fTriggerEff) < 1e-6) && fTriggerEff > 0); + if (needShape) DMSG(fDebug, 0, "Will load shape correction"); + if (needEff) DMSG(fDebug, 0, "Will load trigger efficiency, was=%f, %f", + fTriggerEff, fTriggerEff0); + if(!needShape) { // && !needEff) { + DMSG(fDebug,1,"Objects already set for normalization - no action taken"); return; } - - TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/" - "Normalization/normalizationHists_%s_%s.root", - type.Data(),snn.Data())); + + TString fname(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/" + "Normalization/normalizationHists_%s_%s.root", + type.Data(),snn.Data())); + AliWarningF("Using old-style corrections from %s", fname.Data()); + TFile* fin = TFile::Open(fname, "READ"); if(!fin) { - AliWarning(Form("no file for normalization of %d/%d", sys, energy)); + AliWarningF("no file for normalization of %d/%d (%s)", + sys, energy, fname.Data()); return; } // Shape correction - if ((fNormalizationScheme & kShape) && !fShapeCorr) { + if (needShape) { TString trigName("All"); if (fTriggerMask == AliAODForwardMult::kInel || fTriggerMask == AliAODForwardMult::kNClusterGt0) @@ -849,90 +759,84 @@ AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy) } // Trigger efficiency - TString effName(Form("%sTriggerEff", - fTriggerMask == AliAODForwardMult::kInel ? "inel" : - fTriggerMask == AliAODForwardMult::kNSD ? "nsd" : - fTriggerMask == AliAODForwardMult::kInelGt0 ? - "inelgt0" : "all")); - - Double_t trigEff = 1; - if (fNormalizationScheme & kTriggerEfficiency) { + if (needEff) { + TString effName(Form("%sTriggerEff", + fTriggerMask == AliAODForwardMult::kInel ? "inel" : + fTriggerMask == AliAODForwardMult::kNSD ? "nsd" : + fTriggerMask == AliAODForwardMult::kInelGt0 ? + "inelgt0" : "all")); + Double_t trigEff = 1; TObject* eff = fin->Get(effName); if (eff) AliForwardUtil::GetParameter(eff, trigEff); - } - if (fTriggerEff != 1) SetTriggerEff(trigEff); - if (fTriggerEff < 0) fTriggerEff = 1; - // Trigger efficiency - TString eff0Name(Form("%sTriggerEff0", - fTriggerMask == AliAODForwardMult::kInel ? "inel" : - fTriggerMask == AliAODForwardMult::kNSD ? "nsd" : - fTriggerMask == AliAODForwardMult::kInelGt0 ? - "inelgt0" : "all")); - - Double_t trigEff0 = 1; - if (fNormalizationScheme & kTriggerEfficiency) { - TObject* eff = fin->Get(eff0Name); - if (eff) AliForwardUtil::GetParameter(eff, trigEff0); + if (trigEff <= 0) + AliWarningF("Retrieved trigger efficiency %s is %f<=0, ignoring", + effName.Data(), trigEff); + else + SetTriggerEff(trigEff); + + // Trigger efficiency + TString eff0Name(effName); + eff0Name.Append("0"); + + Double_t trigEff0 = 1; + TObject* eff0 = fin->Get(eff0Name); + if (eff0) AliForwardUtil::GetParameter(eff, trigEff0); + if (trigEff0 < 0) + AliWarningF("Retrieved trigger efficiency %s is %f<0, ignoring", + eff0Name.Data(), trigEff0); + else + SetTriggerEff0(trigEff0); } - if (fTriggerEff0 != 1) SetTriggerEff0(trigEff0); - if (fTriggerEff0 < 0) fTriggerEff0 = 1; - + // TEMPORARY FIX // Rescale the shape correction by the trigger efficiency if (fShapeCorr) { AliWarning(Form("Rescaling shape correction by trigger efficiency: " - "1/E_X=1/%f", trigEff)); - fShapeCorr->Scale(1. / trigEff); + "1/E_X=1/%f", fTriggerEff)); + fShapeCorr->Scale(1. / fTriggerEff); } + if (fin) fin->Close(); // Print - out - if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization."); + if (fDebug > 1 && fShapeCorr && fTriggerEff) + DMSG(fDebug, 1, "Loaded objects for normalization."); } +#define PF(N,V,...) \ + AliForwardUtil::PrintField(N,V, ## __VA_ARGS__) +#define PFB(N,FLAG) \ + do { \ + AliForwardUtil::PrintName(N); \ + std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \ + } while(false) +#define PFV(N,VALUE) \ + do { \ + AliForwardUtil::PrintName(N); \ + std::cout << (VALUE) << std::endl; } while(false) + //________________________________________________________________________ void -AliBasedNdetaTask::Print(Option_t*) const +AliBasedNdetaTask::Print(Option_t* option) const { // // Print information // - std::cout << this->ClassName() << ": " << this->GetName() << "\n" - << std::boolalpha - << " Trigger: " << (fTriggerString ? - fTriggerString->GetTitle() : - "none") << "\n" - << " Vertex range: [" << fVtxMin << ":" - << fVtxMax << "]\n" - << " Rebin factor: " << fRebin << "\n" - << " Cut edges: " << fCutEdges << "\n" - << " Symmertrice: " << fSymmetrice << "\n" - << " Use TH2::ProjectionX: " << fUseROOTProj << "\n" - << " Correct for empty: " << fCorrEmpty << "\n" - << " Normalization scheme: " << (fSchemeString ? - fSchemeString->GetTitle() : - "none") <<"\n" - << " Trigger efficiency: " << fTriggerEff << "\n" - << " Bin-0 Trigger efficiency: " << fTriggerEff0 << "\n" - << " Shape correction: " << (fShapeCorr ? - fShapeCorr->GetName() : - "none") << "\n" - << " sqrt(s_NN): " << (fSNNString ? - fSNNString->GetTitle() : - "unknown") << "\n" - << " Collision system: " << (fSysString ? - fSysString->GetTitle() : - "unknown") << "\n" - << " Centrality bins: " << (fCentAxis ? "" : "none"); - if (fCentAxis) { - Int_t nBins = fCentAxis->GetNbins(); - const Double_t* bins = fCentAxis->GetXbins()->GetArray(); - for (Int_t i = 0; i <= nBins; i++) - std::cout << (i==0 ? "" : "-") << bins[i]; - } - std::cout << std::noboolalpha << std::endl; - + AliBaseAODTask::Print(option); + TString schemeString(NormalizationSchemeString(fNormalizationScheme)); + + gROOT->IncreaseDirLevel(); + PFV("Rebin factor", fRebin); + PFB("Cut edges", fCutEdges); + PFB("Symmertrice", fSymmetrice); + PFB("Use TH2::ProjectionX", fUseROOTProj); + PFB("Correct for empty", fCorrEmpty); + PFV("Normalization scheme", schemeString ); + PFV("Trigger efficiency", fTriggerEff); + PFV("Bin-0 Trigger efficiency", fTriggerEff0); + PFV("Shape correction", (fShapeCorr?fShapeCorr->GetName():"none"));; + gROOT->DecreaseDirLevel(); } //________________________________________________________________________ @@ -963,6 +867,8 @@ AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges) TH1D* tmp = static_cast(h->Clone(Form("%s_rebin%02d", h->GetName(), rebin))); tmp->Rebin(rebin); + // Hist should be reset, as it otherwise messes with the cutEdges option + tmp->Reset(); tmp->SetDirectory(0); // The new number of bins @@ -975,7 +881,11 @@ AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges) for(Int_t j = 1; j<=rebin;j++) { Int_t bin = (i-1)*rebin + j; Double_t c = h->GetBinContent(bin); - if (c <= 0) continue; + if (c <= 0) { + continue; // old TODO: check + //content = -1; // new + //break; // also new + } if (cutEdges) { if (h->GetBinContent(bin+1)<=0 || @@ -1137,11 +1047,15 @@ AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col) //____________________________________________________________________ TString -AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post) +AliBasedNdetaTask::Sum::GetHistName(const char* /*name*/, + Int_t what, const char* post) { - TString n(name); - if (what == 1) n.Append("0"); - else if (what == 2) n.Append("Events"); + TString n; + switch (what) { + case 0: n = "sum"; break; + case 1: n = "sum0"; break; + case 2: n = "events"; break; + } if (post && post[0] != '\0') n.Append(post); return n; } @@ -1175,19 +1089,20 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, { DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName()); + // The return value `ret' is not scaled in anyway TH2D* ret = static_cast(fSum->Clone(fSum->GetName())); ret->SetDirectory(0); - ret->Reset(); Int_t n = Int_t(fEvents->GetBinContent(1)); Int_t n0 = Int_t(fEvents->GetBinContent(2)); - - AliInfoF("Adding histograms %s(%d) and %s(%d) with weights %f and %f resp.", - fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon, 1./epsilon0); - DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.", - fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0); - // Generate merged histogram - ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon); - ntotal = n / epsilon + n0 / epsilon0; + ntotal = n; + if (n0 > 0) { + ret->Reset(); + DMSG(fDebug,1, + "Adding histograms %s(%d) and %s(%d) w/weights %f and %f resp.", + fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon,1./epsilon0); + ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon); + ntotal = n / epsilon + n0 / epsilon0; + } TList* out = new TList; out->SetOwner(); @@ -1197,6 +1112,8 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, output->Add(out); // Now make copies, normalize them, and store in output list + // Note, these are the only ones normalized here + // These are mainly for diagnostics TH2D* sumCopy = static_cast(fSum->Clone("sum")); TH2D* sum0Copy = static_cast(fSum0->Clone("sum0")); TH2D* retCopy = static_cast(ret->Clone("sumAll")); @@ -1212,41 +1129,60 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, TH1D* norm = ProjectX(fSum, "norm", o, o, rootProj, corrEmpty, false); TH1D* norm0 = ProjectX(fSum0, "norm0", o, o, rootProj, corrEmpty, false); TH1D* normAll = ProjectX(ret, "normAll", o, o, rootProj, corrEmpty, false); + norm->SetTitle("#eta coverage - >0-bin"); + norm0->SetTitle("#eta coverage - 0-bin"); + normAll->SetTitle("#eta coverage"); norm->SetDirectory(0); norm0->SetDirectory(0); normAll->SetDirectory(0); - ScaleToCoverage(sumCopy, norm); - ScaleToCoverage(sum0Copy, norm0); - ScaleToCoverage(retCopy, normAll); - - TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1, nY,rootProj,corrEmpty); - TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1, nY,rootProj,corrEmpty); - TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1, nY,rootProj,corrEmpty); - sumCopyPx->SetDirectory(0); + TH1D* sumCopyPx = ProjectX(sumCopy, "average", 1,nY,rootProj,corrEmpty); + TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0", 1,nY,rootProj,corrEmpty); + TH1D* retCopyPx = ProjectX(retCopy, "averageAll", 1,nY,rootProj,corrEmpty); + sumCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sumCopy->GetTitle())); + sum0CopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sum0Copy->GetTitle())); + retCopyPx-> SetTitle(Form("#sum_{i}^{N_{#phi}}%s", retCopy->GetTitle())); + sumCopyPx-> SetDirectory(0); sum0CopyPx->SetDirectory(0); - retCopyPx->SetDirectory(0); - - TH1D* phi = ProjectX(fSum, "phi", nY+1, nY+1,rootProj,corrEmpty); - TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1, nY+1,rootProj,corrEmpty); - TH1D* phiAll = ProjectX(ret, "phiAll", nY+1, nY+1,rootProj,corrEmpty); - phi->SetDirectory(0); - phi0->SetDirectory(0); + retCopyPx-> SetDirectory(0); + + TH1D* phi = ProjectX(fSum, "phi", nY+1,nY+1,rootProj,corrEmpty,false); + TH1D* phi0 = ProjectX(fSum0, "phi0", nY+1,nY+1,rootProj,corrEmpty,false); + TH1D* phiAll = ProjectX(ret, "phiAll", nY+1,nY+1,rootProj,corrEmpty,false); + phi ->SetTitle("#phi acceptance from dead strips - >0-bin"); + phi0 ->SetTitle("#phi acceptance from dead strips - 0-bin"); + phiAll->SetTitle("#phi acceptance from dead strips"); + phi ->SetDirectory(0); + phi0 ->SetDirectory(0); phiAll->SetDirectory(0); + const TH1D* cov = (corrEmpty ? norm : phi); + const TH1D* cov0 = (corrEmpty ? norm0 : phi0); + const TH1D* covAll = (corrEmpty ? normAll : phiAll); + + // Here, we scale to the coverage (or phi acceptance) + ScaleToCoverage(sumCopy, cov); + ScaleToCoverage(sum0Copy, cov0); + ScaleToCoverage(retCopy, covAll); + // Scale our 1D histograms - sumCopyPx->Scale(1., "width"); + sumCopyPx ->Scale(1., "width"); sum0CopyPx->Scale(1., "width"); - retCopyPx->Scale(1., "width"); + retCopyPx ->Scale(1., "width"); - AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(), - sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum())); + DMSG(fDebug,2,"Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(), + sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()); // Scale the normalization - they should be 1 at the maximum - norm->Scale(n > 0 ? 1. / n : 1); - norm0->Scale(n0 > 0 ? 1. / n0 : 1); + norm ->Scale(n > 0 ? 1. / n : 1); + norm0 ->Scale(n0 > 0 ? 1. / n0 : 1); normAll->Scale(ntotal > 0 ? 1. / ntotal : 1); + // Scale the normalization - they should be 1 at the maximum + phi ->Scale(n > 0 ? 1. / n : 1); + phi0 ->Scale(n0 > 0 ? 1. / n0 : 1); + phiAll->Scale(ntotal > 0 ? 1. / ntotal : 1); + out->Add(sumCopy); out->Add(sum0Copy); out->Add(retCopy); @@ -1260,10 +1196,16 @@ AliBasedNdetaTask::Sum::CalcSum(TList* output, out->Add(phi0); out->Add(phiAll); - AliInfoF("Returning (1/%f * %s + 1/%f * %s), " + if (fDebug >= 1) { + if (n0 > 0) + DMSG(fDebug,1,"Returning (1/%f * %s + 1/%f * %s), " "1/%f * %d + 1/%f * %d = %d", epsilon0, fSum0->GetName(), epsilon, fSum->GetName(), epsilon0, n0, epsilon, n, int(ntotal)); + else + DMSG(fDebug,1, "Returning %s, %d", fSum->GetName(), int(ntotal)); + } + #if 0 for (Int_t i = 1; i <= ret->GetNbinsX(); i++) { Double_t nc = sum->GetBinContent(i, 0); @@ -1283,9 +1225,11 @@ AliBasedNdetaTask::CentralityBin::CentralityBin() fSum(0), fSumMC(0), fTriggers(0), + fStatus(0), fLow(0), fHigh(0), - fDoFinalMCCorrection(false), + fDoFinalMCCorrection(false), + fSatelliteVertices(false), fDebug(0) { // @@ -1302,9 +1246,11 @@ AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name, fSum(0), fSumMC(0), fTriggers(0), + fStatus(0), fLow(low), fHigh(high), fDoFinalMCCorrection(false), + fSatelliteVertices(false), fDebug(0) { // @@ -1336,9 +1282,11 @@ AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o) fSum(o.fSum), fSumMC(o.fSumMC), fTriggers(o.fTriggers), + fStatus(o.fStatus), fLow(o.fLow), fHigh(o.fHigh), - fDoFinalMCCorrection(o.fDoFinalMCCorrection), + fDoFinalMCCorrection(o.fDoFinalMCCorrection), + fSatelliteVertices(o.fSatelliteVertices), fDebug(o.fDebug) { // @@ -1383,9 +1331,11 @@ AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o) fSum = o.fSum; fSumMC = o.fSumMC; fTriggers = o.fTriggers; + fStatus = o.fStatus; fLow = o.fLow; fHigh = o.fHigh; fDoFinalMCCorrection = o.fDoFinalMCCorrection; + fSatelliteVertices = o.fSatelliteVertices; return *this; } @@ -1431,7 +1381,12 @@ AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir, Int_t mask) fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers", mask); fTriggers->SetDirectory(0); + + fStatus = AliAODForwardMult::MakeStatusHistogram("status"); + fStatus->SetDirectory(0); + fSums->Add(fTriggers); + fSums->Add(fStatus); } //____________________________________________________________________ void @@ -1450,15 +1405,16 @@ AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc) TString sn = Sum::GetHistName(GetName(),0,post); TString sn0 = Sum::GetHistName(GetName(),1,post); TString ev = Sum::GetHistName(GetName(),2,post); - TH2D* sum = static_cast(list->FindObject(sn)); - TH2D* sum0 = static_cast(list->FindObject(sn0)); - TH1I* events = static_cast(list->FindObject(ev)); + TH2D* sum = static_cast(list->FindObject(sn)); + TH2D* sum0 = static_cast(list->FindObject(sn0)); + TH1I* events = static_cast(list->FindObject(ev)); if (!sum || !sum0 || !events) { - AliWarningF("Failed to find one or more histograms: " - "%s (%p) %s (%p) %s (%p)", - sn.Data(), sum, - sn0.Data(), sum0, - ev.Data(), events); + if (!mc) + AliWarningF("Failed to find one or more histograms: " + "%s (%p) %s (%p) %s (%p)", + sn.Data(), sum, + sn0.Data(), sum0, + ev.Data(), events); return false; } Sum* ret = new Sum(GetName(), post); @@ -1518,12 +1474,13 @@ AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward, DGUARD(fDebug,2,"Check the event"); // We do not check for centrality here - it's already done - return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, fTriggers); + return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0, + fTriggers, fStatus); } //____________________________________________________________________ -void +Bool_t AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward, Int_t triggerMask, Bool_t isZero, Double_t vzMin, Double_t vzMax, @@ -1542,12 +1499,14 @@ AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward, // DGUARD(fDebug,1,"Process one event for %s a given centrality bin", data ? data->GetName() : "(null)"); - if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return; - if (!data) return; + if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return false; + if (!data) return false; if (!fSum) CreateSums(data, mc); fSum->Add(data, isZero); if (mc) fSumMC->Add(mc, isZero); + + return true; } //________________________________________________________________________ @@ -1567,8 +1526,8 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, // trigEff From MC // ntotal On return, contains the number of events. // - DGUARD(fDebug,1,"Normalize centrality bin %s with %s", - GetName(), t.GetName()); + DGUARD(fDebug,1,"Normalize centrality bin %s [%3d-%3d%%] with %s", + GetName(), fLow, fHigh, t.GetName()); Double_t nAll = t.GetBinContent(AliAODForwardMult::kBinAll); Double_t nB = t.GetBinContent(AliAODForwardMult::kBinB); Double_t nA = t.GetBinContent(AliAODForwardMult::kBinA); @@ -1599,10 +1558,10 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, if (scheme & kEventLevel) { ntotal = nAccepted / vtxEff; scaler = vtxEff; - AliInfoF("Calculating event normalisation as\n" - " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)", - Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex), - ntotal, scaler); + DMSG(fDebug,0,"Calculating event normalisation as\n" + " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)", + Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex), + ntotal, scaler); if (scheme & kBackground) { // 1 E_V E_V // s = --------- = ------------- = ------------ @@ -1621,23 +1580,24 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, // scaler = 1. / (1. / vtxEff - beta / nWithVertex); // A simpler expresion scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689 - AliInfo(Form("Calculating event normalisation as\n" - " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" - " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)", - Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), - nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta), - Int_t(nWithVertex), ntotal, scaler)); + DMSG(fDebug,0,"Calculating event normalisation as\n" + " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" + " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)", + Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), + nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta), + Int_t(nWithVertex), ntotal, scaler); rhs.Append("(1/eps_V - beta/N_vtx)"); } // Background else rhs.Append("/eps_V"); } // Event-level if (scheme & kTriggerEfficiency) { - ntotal /= trigEff; - scaler *= trigEff; - AliInfo(Form("Correcting for trigger efficiency:\n" - " N = 1 / E_X * N = 1 / %f * %d = %f (%f)", - trigEff, Int_t(ntotal), ntotal / trigEff, scaler)); + Int_t old = Int_t(ntotal); + ntotal /= trigEff; + scaler *= trigEff; + DMSG(fDebug,0,"Correcting for trigger efficiency:\n" + " N = 1 / E_X * N = 1 / %f * %d = %f (%f)", + trigEff, old, ntotal, scaler); rhs.Append("/eps_T"); } // Trigger efficiency } @@ -1656,51 +1616,54 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, - beta / nWithVertex)); scaler = nWithVertex / (nWithVertex + 1/trigEff * (nTriggered-nWithVertex-beta)); - AliInfo(Form("Calculating event normalisation as\n" - " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" - " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = " - "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)", - Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), - Int_t(nAccepted), trigEff, Int_t(nTriggered), - Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex), - ntotal, scaler)); + DMSG(fDebug,0,"Calculating event normalisation as\n" + " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n" + " N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = " + "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)", + Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta), + Int_t(nAccepted), trigEff, Int_t(nTriggered), + Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex), + ntotal, scaler); rhs.Append("(1+1/eps_T(1/eps_V-1-beta/N_vtx))"); } if (text) { - text->Append(Form("%-40s = %d\n", "N_all", UInt_t(nAll))); - text->Append(Form("%-40s = %d\n", "N_acc", UInt_t(nAccepted))); - text->Append(Form("%-40s = %d\n", "N_trg", UInt_t(nTriggered))); - text->Append(Form("%-40s = %d\n", "N_vtx", UInt_t(nWithVertex))); - text->Append(Form("%-40s = %d\n", "N_B", UInt_t(nB))); - text->Append(Form("%-40s = %d\n", "N_A", UInt_t(nA))); - text->Append(Form("%-40s = %d\n", "N_C", UInt_t(nC))); - text->Append(Form("%-40s = %d\n", "N_E", UInt_t(nE))); + text->Append(Form("%-40s = %d\n", "N_all", UInt_t(nAll))); + text->Append(Form("%-40s = %d\n", "N_acc", UInt_t(nAccepted))); + text->Append(Form("%-40s = %d\n", "N_trg", UInt_t(nTriggered))); + text->Append(Form("%-40s = %d\n", "N_vtx", UInt_t(nWithVertex))); + text->Append(Form("%-40s = %d\n", "N_B", UInt_t(nB))); + text->Append(Form("%-40s = %d\n", "N_A", UInt_t(nA))); + text->Append(Form("%-40s = %d\n", "N_C", UInt_t(nC))); + text->Append(Form("%-40s = %d\n", "N_E", UInt_t(nE))); text->Append(Form("%-40s = %d\n", "beta = N_A + N_C - 2N_E",UInt_t(beta))); text->Append(Form("%-40s = %f\n", "eps_V = N_vtx/N_trg",vtxEff)); - text->Append(Form("%-40s = %f\n", "eps_T", trigEff)); - text->Append(Form("%-40s = %f\n", rhs.Data(), ntotal)); + text->Append(Form("%-40s = %f\n", "eps_T", trigEff)); + text->Append(Form("%-40s = %f\n", rhs.Data(), ntotal)); } - - AliInfo(Form("\n" - " Total of %9d events for %s\n" - " of these %9d have an offline trigger\n" - " of these N_T = %9d has the selected trigger\n" - " of these N_V = %9d has a vertex\n" - " of these N_A = %9d were in the selected range\n" - " Triggers by hardware type:\n" - " N_b = %9d\n" - " N_ac = %9d (%d+%d)\n" + TString tN = t.GetXaxis()->GetBinLabel(AliAODForwardMult::kWithTrigger); + tN.ReplaceAll("w/Selected trigger (",""); + tN.ReplaceAll(")", ""); + DMSG(fDebug,0,"\n" + " Total of %9d events for %s\n" + " of these %9d have an offline trigger\n" + " of these N_T = %9d has the selected trigger (%s)\n" + " of these N_V = %9d has a vertex\n" + " of these N_A = %9d were in the selected range\n" + " Triggers by hardware type:\n" + " N_b = %9d\n" + " N_ac = %9d (%d+%d)\n" " N_e = %9d\n" - " Vertex efficiency: %f\n" - " Trigger efficiency: %f\n" - " Total number of events: N = %f\n" - " Scaler (N_A/N): %f\n" - " %25s = %f", - Int_t(nAll), GetTitle(), Int_t(nOffline), - Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted), - Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE), - vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal)); + " Vertex efficiency: %f\n" + " Trigger efficiency: %f\n" + " Total number of events: N = %f\n" + " Scaler (N_A/N): %f\n" + " %25s = %f", + Int_t(nAll), GetTitle(), Int_t(nOffline), + Int_t(nTriggered), tN.Data(), + Int_t(nWithVertex), Int_t(nAccepted), + Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE), + vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal); return scaler; } @@ -1711,26 +1674,33 @@ AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin, const char* postfix) const { static TString n; - n = Form("dndeta%s%s",GetName(), postfix); + n = GetName(); + n.ReplaceAll("dNdeta", ""); + n.Prepend("dndeta"); + n.Append(postfix); + // n = Form("dndeta%s%s",GetName(), postfix); if (rebin > 1) n.Append(Form("_rebin%02d", rebin)); if (sym) n.Append("_mirror"); return n.Data(); } //________________________________________________________________________ TH1* -AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin, - Bool_t sym, - const char* postfix) const +AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin, + Bool_t sym, + const char* postfix, + Bool_t verbose) const { if (!fOutput) { - AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(), - fLow, fHigh)); + AliWarningF("No output list defined in %s [%3d,%3d]", GetName(), + fLow, fHigh); return 0; } TString n = GetResultName(rebin, sym, postfix); TObject* o = fOutput->FindObject(n.Data()); if (!o) { - // AliWarning(Form("Object %s not found in output list", n.Data())); + if (verbose) + AliWarningF("Object %s not found in output list of %s", + n.Data(), GetName()); return 0; } return static_cast(o); @@ -1767,17 +1737,32 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, // cutEdges Whether to cut edges when rebinning // DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName()); - TH2D* copy = static_cast(sum->Clone(Form("d2Ndetadphi%s%s", - GetName(), postfix))); + TString base(GetName()); + base.ReplaceAll("dNdeta", ""); + base.Append(postfix); + TH2D* copy = static_cast(sum->Clone(Form("d2Ndetadphi%s", + base.Data()))); + + TH1D* accNorm = 0; Int_t nY = sum->GetNbinsY(); + // Hack HHD Hans test code to manually remove FMD2 dead channel (but + // it is on outer?) + // + // cholm comment: The original hack has been moved to + // AliForwarddNdetaTask::CheckEventData - this simplifies things a + // great deal, and we could in principle use the new phi acceptance. + // + // However, since we may have filtered out the dead sectors in the + // AOD production already, we can't be sure we can recalculate the + // phi acceptance correctly, so for now, we rely on fCorrEmpty being set. Int_t o = (corrEmpty ? 0 : nY+1); - TH1D* accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), o, o, - rootProj, corrEmpty, false); + accNorm = ProjectX(sum, Form("norm%s",base.Data()), o, o, + rootProj, corrEmpty, false); accNorm->SetDirectory(0); // ---- Scale by shape correction ---------------------------------- if (shapeCorr) copy->Divide(shapeCorr); - else AliInfo("No shape correction specified, or disabled"); + else DMSG(fDebug,1,"No shape correction specified, or disabled"); // --- Normalize to the coverage ----------------------------------- if (corrEmpty) { @@ -1787,7 +1772,7 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, } // --- Project on X axis ------------------------------------------- - TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix), + TH1D* dndeta = ProjectX(copy, Form("dndeta%s",base.Data()), 1, nY, rootProj, corrEmpty); dndeta->SetDirectory(0); // Event-level normalization @@ -1798,32 +1783,76 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, dndeta->Scale(1., "width"); copy->Scale(1., "width"); - TH1D* dndetaMCCorrection = 0; - TList* centlist = 0; - TH1D* dndetaMCtruth = 0; - TList* truthcentlist = 0; + TH1D* dndetaMCCorrection = 0; + TH1D* dndetaMCtruth = 0; + TList* centlist = 0; + TList* truthcentlist = 0; - // Possible final correction to / - if(mclist) + // --- Possible final correction to / ----- + // we get the rebinned distribution for satellite to make the correction + TString rebinSuf(fSatelliteVertices ? "_rebin05" : ""); + if(mclist) { centlist = static_cast (mclist->FindObject(GetListName())); - if(centlist) - dndetaMCCorrection = - static_cast(centlist->FindObject(Form("dndeta%s%s", - GetName(), postfix))); - if(truthlist) - truthcentlist =static_cast(truthlist->FindObject(GetListName())); - if(truthcentlist) - dndetaMCtruth =static_cast(truthcentlist->FindObject("dndetaTruth")); - - if(dndetaMCCorrection && dndetaMCtruth) { + if(centlist) + dndetaMCCorrection = + static_cast(centlist->FindObject(Form("dndeta%s%s", + base.Data(), + rebinSuf.Data()))); + } + if (truthlist) { + truthcentlist = static_cast(truthlist->FindObject(GetListName())); + if (truthcentlist) + // TODO here new is "dndetaTruth" + dndetaMCtruth = + static_cast(truthcentlist->FindObject(Form("dndetaMCTruth%s", + rebinSuf.Data()))); + } + + if (dndetaMCCorrection && dndetaMCtruth) { AliInfo("Correcting with final MC correction"); + TString testString(dndetaMCCorrection->GetName()); + + // We take the measured MC dN/deta and divide with truth dndetaMCCorrection->Divide(dndetaMCtruth); dndetaMCCorrection->SetTitle("Final MC correction"); dndetaMCCorrection->SetName("finalMCCorr"); - dndeta->Divide(dndetaMCCorrection); + for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) { + if(dndetaMCCorrection->GetBinContent(m) < 0.5 || + dndetaMCCorrection->GetBinContent(m) > 1.75) { + dndetaMCCorrection->SetBinContent(m,1.); + dndetaMCCorrection->SetBinError(m,0.1); + } + } + // Applying the correction + if (!fSatelliteVertices) + // For non-satellites we took the same binning, so we do a straight + // division + dndeta->Divide(dndetaMCCorrection); + else { + // For satelitte events, we took coarser binned histograms, so + // we need to do a bit more + for(Int_t m = 1; m <= dndeta->GetNbinsX(); m++) { + if(dndeta->GetBinContent(m) <= 0.01 ) continue; + + Double_t eta = dndeta->GetXaxis()->GetBinCenter(m); + Int_t bin = dndetaMCCorrection->GetXaxis()->FindBin(eta); + Double_t mccorr = dndetaMCCorrection->GetBinContent(bin); + Double_t mcerror = dndetaMCCorrection->GetBinError(bin); + if (mccorr < 1e-6) { + dndeta->SetBinContent(m, 0); + dndeta->SetBinError(m, 0); + } + Double_t value = dndeta->GetBinContent(m); + Double_t error = dndeta->GetBinError(m); + Double_t sumw2 = (error * error * mccorr * mccorr + + mcerror * mcerror * value * value); + dndeta->SetBinContent(m,value/mccorr) ; + dndeta->SetBinError(m,TMath::Sqrt(sumw2)/mccorr/mccorr); + } + } } else - AliInfo("No final MC correction applied"); + DMSG(fDebug,1,"No final MC correction applied"); // --- Set some histogram attributes ------------------------------- TString post; @@ -1836,13 +1865,55 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, GetName(), post.Data())); // --- Make symmetric extensions and rebinnings -------------------- - if (symmetrice) fOutput->Add(Symmetrice(dndeta)); + if (symmetrice) fOutput->Add(Symmetrice(dndeta)); fOutput->Add(dndeta); fOutput->Add(accNorm); fOutput->Add(copy); fOutput->Add(Rebin(dndeta, rebin, cutEdges)); - if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges))); + if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges))); if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection); + + // HHD Test of dN/deta in phi bins add flag later? + // + // cholm comment: We disable this for now +#if 0 + for (Int_t nn=1; nn <= sum->GetNbinsY(); nn++) { + TH1D* dndeta_phi = ProjectX(copy, Form("dndeta%s_phibin%d", + base.Data(), nn), + nn, nn, rootProj, corrEmpty); + dndeta_phi->SetDirectory(0); + // Event-level normalization + dndeta_phi->Scale(TMath::Pi()/10., "width"); + + if(centlist) + dndetaMCCorrection = + static_cast(centlist->FindObject(Form("dndeta%s_phibin%d", + base.Data(),nn))); + if(truthcentlist) + dndetaMCtruth + = static_cast(truthcentlist->FindObject("dndetaMCTruth")); + + if (dndetaMCCorrection && dndetaMCtruth) { + AliInfo("Correcting with final MC correction"); + TString testString(dndetaMCCorrection->GetName()); + dndetaMCCorrection->Divide(dndetaMCtruth); + dndetaMCCorrection->SetTitle(Form("Final_MC_correction_phibin%d",nn)); + dndetaMCCorrection->SetName(Form("Final_MC_correction_phibin%d",nn)); + for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) { + if(dndetaMCCorrection->GetBinContent(m) < 0.25 || + dndetaMCCorrection->GetBinContent(m) > 1.75) { + dndetaMCCorrection->SetBinContent(m,1.); + dndetaMCCorrection->SetBinError(m,0.1); + } + } + //Applying the correction + dndeta_phi->Divide(dndetaMCCorrection); + } + fOutput->Add(dndeta_phi); + fOutput->Add(Rebin(dndeta_phi, rebin, cutEdges)); + if(dndetaMCCorrection) fOutput->Add(dndetaMCCorrection); + } // End of phi +#endif } //________________________________________________________________________ @@ -1909,37 +1980,13 @@ AliBasedNdetaTask::CentralityBin::End(TList* sums, // --- Get normalization scaler ------------------------------------ Double_t epsilonT = trigEff; Double_t epsilonT0 = trigEff0; - AliInfoF("Using epsilonT=%f, epsilonT0=%f for %d", - epsilonT, epsilonT0, triggerMask); -#if 0 - // TEMPORARY FIX - if (triggerMask == AliAODForwardMult::kNSD) { - // This is a local change - epsilonT = 0.96; - AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT)); - } - else if (triggerMask == AliAODForwardMult::kInel) { - // This is a local change - epsilonT = 0.934; - AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT)); - } - if (scheme & kZeroBin) { - if (triggerMask==AliAODForwardMult::kInel) - epsilonT0 = 0.785021; // 0.100240; - else if (triggerMask==AliAODForwardMult::kInelGt0) - epsilonT0 = 0; - else if (triggerMask==AliAODForwardMult::kNSD) - epsilonT0 = .706587; - epsilonT = 1; - AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f", - epsilonT0)); - } -#endif + DMSG(fDebug,2,"Using epsilonT=%f, epsilonT0=%f for 0x%x", + epsilonT, epsilonT0, triggerMask); // Get our histograms Double_t nSum = 0; TH2D* sum = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT, - marker, rootProj, corrEmpty); + marker, rootProj, corrEmpty); Double_t nSumMC = 0; TH2D* sumMC = 0; if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC, @@ -1977,6 +2024,37 @@ AliBasedNdetaTask::CentralityBin::End(TList* sums, // if (!IsAllBin()) return; } +//____________________________________________________________________ +Bool_t +AliBasedNdetaTask::ApplyEmpiricalCorrection(const AliAODForwardMult* aod, + TH2D* data) +{ + if (!fglobalempiricalcorrection || !data) + return true; + Float_t zvertex=aod->GetIpZ(); + Int_t binzvertex=fglobalempiricalcorrection->GetXaxis()->FindBin(zvertex); + if(binzvertex<1||binzvertex>fglobalempiricalcorrection->GetNbinsX()) + return false; + for (int i=1;i<=data->GetNbinsX();i++) { + Int_t bincorrection=fglobalempiricalcorrection->GetYaxis() + ->FindBin(data->GetXaxis()->GetBinCenter(i)); + if(bincorrection<1||bincorrection>fglobalempiricalcorrection->GetNbinsY()) + return false; + Float_t correction=fglobalempiricalcorrection + ->GetBinContent(binzvertex,bincorrection); + if(correction<0.001) { + data->SetBinContent(i,0,0); + data->SetBinContent(i,data->GetNbinsY()+1,0); + } + for(int j=1;j<=data->GetNbinsY();j++) { + if (data->GetBinContent(i,j)>0.0) { + data->SetBinContent(i,j,data->GetBinContent(i,j)*correction); + data->SetBinError(i,j,data->GetBinError(i,j)*correction); + } + } + } + return true; +} // // EOF