From 5bb5d1f658f8370e7c1c652174efa7fd10c71943 Mon Sep 17 00:00:00 2001 From: cholm Date: Wed, 6 Apr 2011 13:55:47 +0000 Subject: [PATCH] Mega commit. In general, more diagnostics histograms. Event inspector can now flag #_SPD_Cluster > 0 events New implementation of sharing algorithm that zero's candidates if not used. dN/deta tasks now outputs summary stacks. --- PWG2/FORWARD/analysis2/AliAODForwardMult.cxx | 19 +- PWG2/FORWARD/analysis2/AliAODForwardMult.h | 55 +- PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx | 194 ++-- PWG2/FORWARD/analysis2/AliBasedNdetaTask.h | 12 +- PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx | 63 +- PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h | 45 + PWG2/FORWARD/analysis2/AliFMDCorrector.cxx | 38 +- PWG2/FORWARD/analysis2/AliFMDCorrector.h | 2 +- .../analysis2/AliFMDDensityCalculator.cxx | 240 +++-- .../analysis2/AliFMDDensityCalculator.h | 57 +- .../analysis2/AliFMDEnergyFitterTask.cxx | 43 +- .../analysis2/AliFMDEventInspector.cxx | 112 ++- PWG2/FORWARD/analysis2/AliFMDEventInspector.h | 11 +- .../FORWARD/analysis2/AliFMDHistCollector.cxx | 104 ++- PWG2/FORWARD/analysis2/AliFMDHistCollector.h | 50 +- PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx | 24 +- PWG2/FORWARD/analysis2/AliFMDMCCorrector.h | 25 +- .../analysis2/AliFMDMCEventInspector.cxx | 2 +- .../analysis2/AliFMDMCSharingFilter.cxx | 84 +- .../FORWARD/analysis2/AliFMDMCSharingFilter.h | 18 +- .../FORWARD/analysis2/AliFMDSharingFilter.cxx | 450 +++++++-- PWG2/FORWARD/analysis2/AliFMDSharingFilter.h | 67 +- .../analysis2/AliForwardMCCorrectionsTask.cxx | 185 ++-- .../analysis2/AliForwardMCCorrectionsTask.h | 5 +- .../AliForwardMCMultiplicityTask.cxx | 79 +- .../analysis2/AliForwardMCMultiplicityTask.h | 8 + .../analysis2/AliForwardMultiplicityBase.cxx | 80 ++ .../analysis2/AliForwardMultiplicityBase.h | 13 +- .../analysis2/AliForwardMultiplicityTask.cxx | 43 +- .../analysis2/AliForwardMultiplicityTask.h | 1 + PWG2/FORWARD/analysis2/AliForwardUtil.cxx | 1 + PWG2/FORWARD/analysis2/AliForwardUtil.h | 3 +- .../analysis2/AliForwarddNdetaTask.cxx | 139 ++- PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h | 9 +- PWG2/FORWARD/analysis2/DrawdNdeta.C | 269 +++--- PWG2/FORWARD/analysis2/ForwardAODConfig.C | 38 +- PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C | 851 ++++++++++++++++++ PWG2/FORWARD/analysis2/OtherData.C | 4 +- 38 files changed, 2814 insertions(+), 629 deletions(-) create mode 100644 PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C diff --git a/PWG2/FORWARD/analysis2/AliAODForwardMult.cxx b/PWG2/FORWARD/analysis2/AliAODForwardMult.cxx index acb8908dba5..e6f6029ab6f 100644 --- a/PWG2/FORWARD/analysis2/AliAODForwardMult.cxx +++ b/PWG2/FORWARD/analysis2/AliAODForwardMult.cxx @@ -35,7 +35,8 @@ AliAODForwardMult::AliAODForwardMult() fHist(), fTriggers(0), fIpZ(fgkInvalidIpZ), - fCentrality(-1) + fCentrality(-1), + fNClusters(0) { // // Constructor @@ -49,7 +50,8 @@ AliAODForwardMult::AliAODForwardMult(Bool_t isMC) 200, -4, 6, 20, 0, 2*TMath::Pi()), fTriggers(0), fIpZ(fgkInvalidIpZ), - fCentrality(-1) + fCentrality(-1), + fNClusters(0) { // // Constructor @@ -87,8 +89,9 @@ AliAODForwardMult::Clear(Option_t* option) // option Passed to TH1::Reset // fHist.Reset(option); - fTriggers = 0; - fIpZ = fgkInvalidIpZ; + fTriggers = 0; + fIpZ = fgkInvalidIpZ; + fNClusters = 0; } //____________________________________________________________________ void @@ -160,13 +163,16 @@ AliAODForwardMult::Browse(TBrowser* b) static TObjString ipz; static TObjString trg; static TObjString cnt; + static TObjString ncl; ipz = Form("ip_z=%fcm", fIpZ); trg = GetTriggerString(fTriggers); cnt = Form("%+6.1f%%", fCentrality); + ncl = Form("%d clusters", fNClusters); b->Add(&fHist); b->Add(&ipz); b->Add(&trg); b->Add(&cnt); + b->Add(&ncl); } //____________________________________________________________________ @@ -189,6 +195,7 @@ AliAODForwardMult::GetTriggerString(UInt_t mask) if ((mask & kC) != 0x0) trg.Append("C "); if ((mask & kE) != 0x0) trg.Append("E "); if ((mask & kMCNSD) != 0x0) trg.Append("MCNSD "); + if ((mask & kNClusterGt0) != 0x0) trg.Append("NCluster>0 "); return trg.Data(); } @@ -224,6 +231,7 @@ AliAODForwardMult::MakeTriggerHistogram(const char* name) ret->GetXaxis()->SetBinLabel(kBinMCNSD, "NSD (MC truth)"); ret->GetXaxis()->SetBinLabel(kBinPileUp, "w/Pileup"); ret->GetXaxis()->SetBinLabel(kBinOffline, "w/Offline"); + ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1"); ret->GetXaxis()->SetBinLabel(kWithVertex, "w/Vertex"); ret->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger"); ret->GetXaxis()->SetBinLabel(kAccepted, "Accepted by cut"); @@ -251,6 +259,8 @@ AliAODForwardMult::MakeTriggerMask(const char* what) else if (s.CompareTo("A") == 0) trgMask |= AliAODForwardMult::kA; else if (s.CompareTo("C") == 0) trgMask |= AliAODForwardMult::kC; else if (s.CompareTo("E") == 0) trgMask |= AliAODForwardMult::kE; + else if (s.CompareTo("NCluster>0") == 0) + trgMask |= AliAODForwardMult::kNClusterGt0; else AliWarningGeneral("MakeTriggerMask", Form("Unknown trigger %s", s.Data())); @@ -303,6 +313,7 @@ AliAODForwardMult::CheckEvent(Int_t triggerMask, if (IsTriggerBits(kPileUp)) hist->AddBinContent(kBinPileUp); if (IsTriggerBits(kMCNSD)) hist->AddBinContent(kBinMCNSD); if (IsTriggerBits(kOffline)) hist->AddBinContent(kBinOffline); + if (IsTriggerBits(kNClusterGt0)) hist->AddBinContent(kBinNClusterGt0); } // Check if we have an event of interest. if (!IsTriggerBits(triggerMask|kB)) return false; diff --git a/PWG2/FORWARD/analysis2/AliAODForwardMult.h b/PWG2/FORWARD/analysis2/AliAODForwardMult.h index eb6c96eb0c9..19ab1687cd6 100644 --- a/PWG2/FORWARD/analysis2/AliAODForwardMult.h +++ b/PWG2/FORWARD/analysis2/AliAODForwardMult.h @@ -127,7 +127,9 @@ public: /** true NSD from MC */ kMCNSD = 0x400, /** Offline MB triggered */ - kOffline = 0x800 + kOffline = 0x800, + /** At least one SPD cluster */ + kNClusterGt0 = 0x1000 }; /** * Bin numbers in trigger histograms @@ -144,6 +146,7 @@ public: kBinPileUp, kBinMCNSD, kBinOffline, + kBinNClusterGt0, kWithTrigger, kWithVertex, kAccepted @@ -201,13 +204,27 @@ public: */ void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits /** - * Check if bit(s) are set in the trigger mask + * Check if all bit(s) are set in the trigger mask. Note, this is + * an @e and between the bits. If you need an @e or you should use + * the member function IsTriggerOrBits * * @param bits Bits to test for * - * @return + * @return true if all enabled bits in the argument is also set in + * the trigger word */ Bool_t IsTriggerBits(UInt_t bits) const; + /** + * Check if any of bit(s) are enabled in the trigger word. This is + * an @e or between the selected bits. If you need and @a and you + * should use the member function IsTriggerBits; + * + * @param bits Bits to check for + * + * @return true if any of the enabled bits in the arguments are also + * enabled in the trigger mask + */ + Bool_t IsTriggerOrBits(UInt_t bits) const; /** * Whether we have any trigger bits */ @@ -316,7 +333,18 @@ public: * @return */ Bool_t HasCentrality() const { return !(fCentrality < 0); } - + /** + * Get the number of SPD clusters seen in @f$ |\eta|<1@f$ + * + * @return Number of SPD clusters seen + */ + UShort_t GetNClusters() const { return fNClusters; } + /** + * Set the number of SPD clusters seen in @f$ |\eta|<1@f$ + * + * @param n Number of SPD clusters + */ + void SetNClusters(UShort_t n) { fNClusters = n; } /** * Get the name of the object * @@ -385,14 +413,15 @@ public: */ static UInt_t MakeTriggerMask(const char* what); protected: - Bool_t fIsMC; // Whether this is from MC - TH2D fHist; // Histogram of d^2N_{ch}/(deta dphi) for this event - UInt_t fTriggers; // Trigger bit mask - Float_t fIpZ; // Z coordinate of the interaction point - Float_t fCentrality; // Event centrality + Bool_t fIsMC; // Whether this is from MC + TH2D fHist; // Histogram of d^2N_{ch}/(deta dphi) for this event + UInt_t fTriggers; // Trigger bit mask + Float_t fIpZ; // Z coordinate of the interaction point + Float_t fCentrality; // Event centrality + UShort_t fNClusters; // Number of SPD clusters in |eta|<1 static const Float_t fgkInvalidIpZ; // Invalid IpZ value - ClassDef(AliAODForwardMult,2); // AOD forward multiplicity + ClassDef(AliAODForwardMult,3); // AOD forward multiplicity }; //____________________________________________________________________ @@ -408,6 +437,12 @@ AliAODForwardMult::IsTriggerBits(UInt_t bits) const { return HasTrigger() && ((fTriggers & bits) == bits); } +//____________________________________________________________________ +inline Bool_t +AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const +{ + return HasTrigger() && ((fTriggers & bits) != 0); +} #endif // Local Variables: diff --git a/PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx b/PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx index bb7cab406fb..8595c84f103 100644 --- a/PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx +++ b/PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx @@ -13,6 +13,7 @@ #include "AliForwardUtil.h" #include "AliAODForwardMult.h" #include +#include //____________________________________________________________________ AliBasedNdetaTask::AliBasedNdetaTask() @@ -153,7 +154,6 @@ AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high) // high High cut // CentralityBin* bin = MakeCentralityBin(GetName(), low, high); - AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at)); fListOfCentralities->AddAtAndExpand(bin, at); } @@ -548,10 +548,67 @@ AliBasedNdetaTask::Terminate(Option_t *) // Loop over centrality bins TIter next(fListOfCentralities); CentralityBin* bin = 0; - while ((bin = static_cast(next()))) + gStyle->SetPalette(1); + THStack* dndetaStack = new THStack("dndeta", "dN/d#eta"); + THStack* dndetaStackRebin = new THStack(Form("dndeta_rebin%02d", fRebin), + "dN_{ch}/d#eta"); + THStack* dndetaMCStack = new THStack("dndetaMC", "dN_{ch}/d#eta"); + THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin), + "dN_{ch}/d#eta"); + + while ((bin = static_cast(next()))) { bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr, fTriggerEff, fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges, - fTriggerMask, GetColor(), GetMarker()); + fTriggerMask, GetMarker()); + 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 (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"); + if (dndeta) dndetaStackRebin->Add(dndeta); + if (dndetaSym) dndetaStackRebin->Add(dndetaSym); + if (dndetaMC) dndetaMCStackRebin->Add(dndetaMC); + if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym); + } + } + // Output the stack + fOutput->Add(dndetaStack); + + // If available output rebinned stack + if (!dndetaStackRebin->GetHists() || + dndetaStackRebin->GetHists()->GetEntries() <= 0) { + AliWarning("No rebinned histograms found"); + delete dndetaStackRebin; + dndetaStackRebin = 0; + } + if (dndetaStackRebin) fOutput->Add(dndetaStackRebin); + + // If available, output track-ref stack + if (!dndetaMCStack->GetHists() || + dndetaMCStack->GetHists()->GetEntries() <= 0) { + AliWarning("No MC histograms found"); + delete dndetaMCStack; + dndetaMCStack = 0; + } + if (dndetaMCStack) fOutput->Add(dndetaMCStack); + + // If available, output rebinned track-ref stack + if (!dndetaMCStackRebin->GetHists() || + dndetaMCStackRebin->GetHists()->GetEntries() <= 0) { + AliWarning("No rebinned MC histograms found"); + delete dndetaMCStackRebin; + dndetaMCStackRebin = 0; + } + if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin); // Output collision energy string if (fSNNString) fOutput->Add(fSNNString->Clone()); @@ -901,6 +958,17 @@ AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o) return *this; } //____________________________________________________________________ +Int_t +AliBasedNdetaTask::CentralityBin::GetColor() const +{ + if (IsAllBin()) return kRed+2; + Float_t fc = (fLow+double(fHigh-fLow)/2) / 100; + Int_t nCol = gStyle->GetNumberOfColors(); + Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5)); + Int_t col = gStyle->GetColorPalette(icol); + return col; +} +//____________________________________________________________________ const char* AliBasedNdetaTask::CentralityBin::GetListName() const { @@ -1109,6 +1177,38 @@ AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t, return scaler; } +//________________________________________________________________________ +const char* +AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin, + Bool_t sym, + const char* postfix) const +{ + static TString n; + 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 +{ + if (!fOutput) { + AliWarning(Form("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())); + return 0; + } + return static_cast(o); +} + //________________________________________________________________________ void AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, @@ -1120,7 +1220,6 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, bool symmetrice, Int_t rebin, bool cutEdges, - Int_t color, Int_t marker) { // @@ -1170,9 +1269,13 @@ AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum, copy->Scale(1., "width"); // --- Set some histogram attributes ------------------------------- - SetHistogramAttributes(dndeta, color, marker, Form("ALICE %s", GetName())); - SetHistogramAttributes(accNorm, color, marker, Form("ALICE %s normalisation", - GetName())); + TString post; + if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix); + SetHistogramAttributes(dndeta, GetColor(), marker, + Form("ALICE %s%s", GetName(), post.Data())); + SetHistogramAttributes(accNorm, GetColor(), marker, + Form("ALICE %s normalisation%s", + GetName(), post.Data())); // --- Make symmetric extensions and rebinnings -------------------- if (symmetrice) fOutput->Add(Symmetrice(dndeta)); @@ -1196,7 +1299,6 @@ AliBasedNdetaTask::CentralityBin::End(TList* sums, Bool_t corrEmpty, Bool_t cutEdges, Int_t triggerMask, - Int_t color, Int_t marker) { // @@ -1256,85 +1358,17 @@ AliBasedNdetaTask::CentralityBin::End(TList* sums, // --- Make result and store --------------------------------------- MakeResult(fSum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0, - scaler, symmetrice, rebin, cutEdges, color, marker); + scaler, symmetrice, rebin, cutEdges, marker); // --- Process result from TrackRefs ------------------------------- - if (fSumMC) { + if (fSumMC) MakeResult(fSumMC, "MC", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0, - scaler, symmetrice, rebin, cutEdges, color+2, marker); - } + scaler, symmetrice, rebin, cutEdges, marker+2); // Temporary stuff - if (!IsAllBin()) return; - - TFile* forward = TFile::Open("forward.root", "READ"); - if (!forward) { - AliWarning(Form("No forward.root file found")); - return; - } - - TH1D* shapeCorrProj = 0; - if (shapeCorr) { - shapeCorrProj = static_cast(shapeCorr)->ProjectionX(); - shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY()); - shapeCorrProj->SetDirectory(0); - fOutput->Add(shapeCorrProj); - } - - TList* official = static_cast(forward->Get("official")); - if (official) { - TH1F* histEta = static_cast(official->FindObject("fHistEta")); - if (histEta) { - TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), - histEta->GetXaxis()->GetXmin(), - histEta->GetXaxis()->GetXmax()); - for (Int_t i = 1; i < histEta->GetNbinsX(); i++) { - oEta->SetBinContent(i, histEta->GetBinContent(i)); - oEta->SetBinError(i, histEta->GetBinError(i)); - } - if (shapeCorrProj) oEta->Divide(shapeCorrProj); - oEta->Scale(1./ntotal, "width"); - oEta->SetDirectory(0); - oEta->SetMarkerStyle(marker+4); - oEta->SetMarkerColor(color+5); - fOutput->Add(oEta); - fOutput->Add(Rebin(oEta, rebin, false)); - } - else - AliWarning(Form("Couldn't find histogram fHistEta in list %s", - official->GetName())); - } - else - AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName())); - - TList* tracks = static_cast(forward->Get("tracks")); - if (tracks) { - TH1F* histEta = static_cast(tracks->FindObject("fHistEta")); - if (histEta) { - TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(), - histEta->GetXaxis()->GetXmin(), - histEta->GetXaxis()->GetXmax()); - for (Int_t i = 1; i < histEta->GetNbinsX(); i++) { - oEta->SetBinContent(i, histEta->GetBinContent(i)); - oEta->SetBinError(i, histEta->GetBinError(i)); - } - if (shapeCorrProj) oEta->Divide(shapeCorrProj); - oEta->Scale(1./ntotal, "width"); - oEta->SetDirectory(0); - oEta->SetMarkerStyle(marker); - oEta->SetMarkerColor(color+5); - fOutput->Add(oEta); - fOutput->Add(Rebin(oEta, rebin, false)); - } - else - AliWarning(Form("Couldn't find histogram fHistEta in list %s", - tracks->GetName())); - } - else - AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName())); + // if (!IsAllBin()) return; - forward->Close(); } // diff --git a/PWG2/FORWARD/analysis2/AliBasedNdetaTask.h b/PWG2/FORWARD/analysis2/AliBasedNdetaTask.h index 281937bb28c..3efcee3a2ad 100644 --- a/PWG2/FORWARD/analysis2/AliBasedNdetaTask.h +++ b/PWG2/FORWARD/analysis2/AliBasedNdetaTask.h @@ -463,7 +463,6 @@ protected: * @param symmetrice Whether to make symmetric extensions * @param rebin Whether to rebin * @param cutEdges Whether to cut edges when rebinning - * @param color Color of markers * @param marker Marker style */ virtual void MakeResult(const TH2D* sum, @@ -475,7 +474,6 @@ protected: bool symmetrice, Int_t rebin, bool cutEdges, - Int_t color, Int_t marker); /** * End of processing @@ -491,7 +489,6 @@ protected: * @param corrEmpty Whether to correct for empty bins * @param cutEdges Whether to cut edges when rebinning * @param triggerMask Trigger mask - * @param color Base colour for markers * @param marker Marker style */ virtual void End(TList* sums, @@ -505,7 +502,6 @@ protected: Bool_t corrEmpty, Bool_t cutEdges, Int_t triggerMask, - Int_t color, Int_t marker); /** * @{ @@ -540,6 +536,14 @@ protected: */ TH1I* GetTrigggers() { return fTriggers; } /** @} */ + + Int_t GetColor() const; + TList* GetResults() const { return fOutput; } + const char* GetResultName(Int_t rebin, Bool_t sym, + const char* postfix="") const; + TH1* GetResult(Int_t rebin, Bool_t sym, + const char* postfix="") const; + protected: /** * Create sum histogram diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx b/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx index e3f42187e82..3991420b80e 100644 --- a/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx @@ -512,6 +512,17 @@ AliFMDCorrELossFit::ELossFit::Draw(Option_t* option) //____________________________________________________________________ #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T)) +//____________________________________________________________________ +Double_t +AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f, + Bool_t includeSigma) const +{ + // + // Return + // Delta - f * (xi + sigma) + return fDelta - f * (fXi + (includeSigma ? fSigma : 0)); +} + //____________________________________________________________________ void AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu, @@ -804,24 +815,24 @@ AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Int_t etabin) const return 0; } if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) { - AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", - etabin, 1, fEtaAxis.GetNbins(), d, r)); + // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", + // etabin, 1, fEtaAxis.GetNbins(), d, r)); return 0; } if (etabin > ringArray->GetEntriesFast()) { - AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", - etabin, 1, ringArray->GetEntriesFast(), d, r)); + // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", + // etabin, 1, ringArray->GetEntriesFast(), d, r)); return 0; } else if (etabin >= ringArray->GetEntriesFast()) { - AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " - "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r, - etabin-1)); + // AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " + // "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r, + // etabin-1)); etabin--; } else if (!ringArray->At(etabin)) { - AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", - etabin, d, r, etabin+1)); + // AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", + // etabin, d, r, etabin+1)); etabin++; } return static_cast(ringArray->At(etabin)); @@ -898,6 +909,38 @@ AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r) return static_cast(fRings.At(idx)); } +//____________________________________________________________________ +Double_t +AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Int_t etabin, + Double_t f, Bool_t showErrors, + Bool_t includeSigma) const +{ + ELossFit* fit = GetFit(d, r, etabin); + if (!fit) { + if (showErrors) { + AliWarning(Form("No fit for FMD%d%c @ etabin=%d", d, r, etabin)); + } + return -1024; + } + return fit->GetLowerBound(f, includeSigma); +} + +//____________________________________________________________________ +Double_t +AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Double_t eta, + Double_t f, Bool_t showErrors, + Bool_t includeSigma) const +{ + Int_t bin = FindEtaBin(eta); + if (bin <= 0) { + if (showErrors) + AliError(Form("eta=%f out of bounds for FMD%d%c", eta, d, r)); + return -1024; + } + return GetLowerBound(d, r, bin, f, showErrors, includeSigma); +} + +//____________________________________________________________________ namespace { TH1D* MakeHist(const TAxis& axis, const char* name, const char* title, Int_t color) @@ -916,6 +959,8 @@ namespace { } } + + #define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O') #define IDX2DET(I) (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3)) //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h b/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h index 4ca91e5be8e..185ad3ffc3a 100644 --- a/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h +++ b/PWG2/FORWARD/analysis2/AliFMDCorrELossFit.h @@ -331,6 +331,15 @@ public: * @return */ const Char_t* GetName() const; + /** + * Calculate the lower bound + * + * @param f Width factor + * @param includeSigma Whether to include sigma + * + * @return @f$ \Delta - f (\xi + \sigma)@f$ + */ + Double_t GetLowerBound(Double_t f, Bool_t includeSigma) const; /** * Calculate the quality */ @@ -526,6 +535,42 @@ public: ELossFit* GetFit(UShort_t d, Char_t r, Int_t etabin) const; /* @} */ + /** + * @{ + * @name Lower bounds on fits + */ + /** + * Get the lower validity bound of the fit + * + * @param d Detector + * @param r Ring + * @param etaBin Eta bin (1-based) + * @param f Factor on xi (and sigma) + * @param showErrors Show errors + * @param includeSigma Whether to include sigma + * + * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$ + */ + Double_t GetLowerBound(UShort_t d, Char_t r, Int_t etaBin, + Double_t f, Bool_t showErrors=true, + Bool_t includeSigma=true) const; + /** + * Get the lower validity bound of the fit + * + * @param d Detector + * @param r Ring + * @param eta Eta value + * @param f Factor on xi (and sigma) + * @param showErrors Show errors + * @param includeSigma Whether to include sigma + * + * @return @f$ \Delta_{mp} - f(\xi+\sigma)@f$ + */ + Double_t GetLowerBound(UShort_t d, Char_t r, Double_t eta, + Double_t f, Bool_t showErrors=true, + Bool_t includeSigma=true) const; + + /** * @{ * @name Miscellaneous diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrector.cxx b/PWG2/FORWARD/analysis2/AliFMDCorrector.cxx index d8fb3cc4dba..62a6f29f3c0 100644 --- a/PWG2/FORWARD/analysis2/AliFMDCorrector.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDCorrector.cxx @@ -15,6 +15,7 @@ #include "AliLog.h" #include #include +#include #include #include @@ -178,8 +179,8 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, if (fUseSecondaryMap) { TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb); if (!bg) { - AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d", - d, r, uvb)); + AliWarning(Form("No secondary correction for FMDM%d%c " + "in vertex bin %d", d, r, uvb)); continue; } // Divide by primary/total ratio @@ -198,18 +199,14 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists, if (fUseAcceptance) { TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb); if (!ac) { - AliWarning(Form("No acceptance correction for FMD%d%c in vertex bin %d", - d, r, uvb)); + AliWarning(Form("No acceptance correction for FMD%d%c in " + "vertex bin %d", d, r, uvb)); continue; } // Divide by the acceptance correction h->Divide(ac); } - - - - if (fUseMergingEfficiency) { if (!fcm.GetMergingEfficiency()) { AliWarning("No merging efficiencies"); @@ -307,8 +304,19 @@ AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents) TIter next(&fRingHistos); RingHistos* o = 0; - while ((o = static_cast(next()))) + THStack* sums = new THStack("sums", "Sums of ring results"); + while ((o = static_cast(next()))) { o->ScaleHistograms(d, nEvents); + TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, + o->fDensity->GetNbinsY(),"e"); + sum->Scale(1., "width"); + sum->SetTitle(o->GetName()); + sum->SetDirectory(0); + sum->SetYTitle("#sum N_{ch,primary}"); + sums->Add(sum); + } + d->Add(sums); + } //____________________________________________________________________ void @@ -372,11 +380,13 @@ AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r) // d detector // r ring // - fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r), - Form("in FMD%d%c", d, r), + fDensity = new TH2D("primaryDensity", + "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)", 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi()); fDensity->SetDirectory(0); + fDensity->SetMarkerColor(Color()); + fDensity->Sumw2(); fDensity->SetXTitle("#eta"); fDensity->SetYTitle("#phi [radians]"); fDensity->SetZTitle("Primary N_{ch} density"); @@ -437,7 +447,7 @@ AliFMDCorrector::RingHistos::Output(TList* dir) //____________________________________________________________________ void -AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t /*nEvents*/) +AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents) { // // Scale the histograms to the total number of events @@ -448,8 +458,8 @@ AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t /*nEvents*/) TList* l = GetOutputList(dir); if (!l) return; - //TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data())); - //if (density) density->Scale(1./nEvents); + TH1* density = GetOutputHist(l,"primaryDensity"); + if (density) density->Scale(1./nEvents); } //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliFMDCorrector.h b/PWG2/FORWARD/analysis2/AliFMDCorrector.h index 4471cbd7103..6665a929ad5 100644 --- a/PWG2/FORWARD/analysis2/AliFMDCorrector.h +++ b/PWG2/FORWARD/analysis2/AliFMDCorrector.h @@ -167,7 +167,7 @@ public: * * @param option Not used */ - void Print(Option_t* option="") const; + virtual void Print(Option_t* option="") const; protected: /** * Internal data structure to keep track of the histograms diff --git a/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx b/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx index 8320b9c6b71..d6654c70e7b 100644 --- a/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx @@ -12,6 +12,7 @@ #include "AliLog.h" #include #include +#include #include #include #include @@ -26,11 +27,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator() : TNamed(), fRingHistos(), fMultCut(0), + fNXi(1), + fIncludeSigma(true), fSumOfWeights(0), fWeightedSum(0), fCorrections(0), fMaxParticles(5), fUsePoisson(false), + fUsePhiAcceptance(false), fAccI(0), fAccO(0), fFMD1iMax(0), @@ -38,6 +42,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator() fFMD2oMax(0), fFMD3iMax(0), fFMD3oMax(0), + fMaxWeights(0), + fLowCuts(0), fEtaLumping(1), fPhiLumping(1), fDebug(0) @@ -53,11 +59,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) : TNamed("fmdDensityCalculator", title), fRingHistos(), fMultCut(0), + fNXi(1), + fIncludeSigma(true), fSumOfWeights(0), fWeightedSum(0), fCorrections(0), fMaxParticles(5), fUsePoisson(false), + fUsePhiAcceptance(false), fAccI(0), fAccO(0), fFMD1iMax(0), @@ -65,6 +74,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) fFMD2oMax(0), fFMD3iMax(0), fFMD3oMax(0), + fMaxWeights(0), + fLowCuts(0), fEtaLumping(5), fPhiLumping(1), fDebug(0) @@ -97,6 +108,17 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title) fAccI = GenerateAcceptanceCorrection('I'); fAccO = GenerateAcceptanceCorrection('O'); + + fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use", + 1, 0, 1, 1, 0, 1); + fMaxWeights->SetXTitle("#eta"); + fMaxWeights->SetDirectory(0); + + fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1); + fLowCuts->SetXTitle("#eta"); + fLowCuts->SetDirectory(0); + + for (Int_t i = 0; i < 5; i++) fMultCuts[i] = 0; } //____________________________________________________________________ @@ -105,11 +127,14 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const : TNamed(o), fRingHistos(), fMultCut(o.fMultCut), + fNXi(o.fNXi), + fIncludeSigma(o.fIncludeSigma), fSumOfWeights(o.fSumOfWeights), fWeightedSum(o.fWeightedSum), fCorrections(o.fCorrections), fMaxParticles(o.fMaxParticles), fUsePoisson(o.fUsePoisson), + fUsePhiAcceptance(o.fUsePhiAcceptance), fAccI(o.fAccI), fAccO(o.fAccO), fFMD1iMax(o.fFMD1iMax), @@ -117,6 +142,8 @@ AliFMDDensityCalculator::AliFMDDensityCalculator(const fFMD2oMax(o.fFMD2oMax), fFMD3iMax(o.fFMD3iMax), fFMD3oMax(o.fFMD3oMax), + fMaxWeights(o.fMaxWeights), + fLowCuts(o.fLowCuts), fEtaLumping(o.fEtaLumping), fPhiLumping(o.fPhiLumping), fDebug(o.fDebug) @@ -157,9 +184,12 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o) TNamed::operator=(o); fMultCut = o.fMultCut; + fNXi = o.fNXi; + fIncludeSigma = o.fIncludeSigma; fDebug = o.fDebug; fMaxParticles = o.fMaxParticles; fUsePoisson = o.fUsePoisson; + fUsePhiAcceptance = o.fUsePhiAcceptance; fAccI = o.fAccI; fAccO = o.fAccO; fFMD1iMax = o.fFMD1iMax; @@ -167,6 +197,8 @@ AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o) fFMD2oMax = o.fFMD2oMax; fFMD3iMax = o.fFMD3iMax; fFMD3oMax = o.fFMD3oMax; + fMaxWeights = o.fMaxWeights; + fLowCuts = o.fLowCuts; fEtaLumping = o.fEtaLumping; fPhiLumping = o.fPhiLumping; fRingHistos.Delete(); @@ -220,10 +252,27 @@ AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const return static_cast(fRingHistos.At(idx)); } - + +//____________________________________________________________________ +void +AliFMDDensityCalculator::SetMultCuts(Double_t fmd1i, + Double_t fmd2i, + Double_t fmd2o, + Double_t fmd3i, + Double_t fmd3o) +{ + fMultCuts[0] = fmd1i; + fMultCuts[1] = fmd2i; + fMultCuts[2] = fmd2o; + fMultCuts[3] = fmd3i; + fMultCuts[4] = fmd3o; + fMultCut = (fmd1i+fmd2i+fmd2o+fmd3i+fmd3o) / 5; +} + //____________________________________________________________________ Double_t -AliFMDDensityCalculator::GetMultCut() const +AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta, + Bool_t errors) const { // // Get the multiplicity cut. If the user has set fMultCut (via @@ -233,11 +282,35 @@ AliFMDDensityCalculator::GetMultCut() const // Return: // Lower cut on multiplicity // + Int_t idx = (d == 1 ? 0 : 2*(d - 2) + 1 + ((r=='I' || r=='i') ? 0 : 1)); + if (fMultCuts[idx] > 0) return fMultCuts[idx]; if (fMultCut > 0) return fMultCut; AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); AliFMDCorrELossFit* fits = fcm.GetELossFit(); - return fits->GetLowCut(); + if (fNXi < 0) return fits->GetLowCut(); + + return fits->GetLowerBound(d, r, ieta, fNXi, errors, fIncludeSigma); +} + +//____________________________________________________________________ +Double_t +AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta, + Bool_t /*errors*/) const +{ + // + // Get the multiplicity cut. If the user has set fMultCut (via + // SetMultCut) then that value is used. If not, then the lower + // value of the fit range for the enery loss fits is returned. + // + // Return: + // Lower cut on multiplicity + // + AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); + AliFMDCorrELossFit* fits = fcm.GetELossFit(); + Int_t iEta = fits->FindEtaBin(eta); + + return GetMultCut(d, r, iEta); } //____________________________________________________________________ @@ -276,22 +349,27 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, for (UShort_t s=0; sfTotalStrips->Fill(eta, phi); - if (mult == AliESDFMD::kInvalidMult || mult > 20) { + if (mult == AliESDFMD::kInvalidMult || mult > 20) { rh->fEmptyStrips->Fill(eta,phi); + rh->fEvsM->Fill(mult,0); continue; } - - Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux); + + Double_t n = 0; + if (cut > 0 && mult > cut) + n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux); rh->fEvsN->Fill(mult,n); rh->fEtaVsN->Fill(eta, n); - Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux); + Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux); fCorrections->Fill(c); if (c > 0) n /= c; rh->fEvsM->Fill(mult,n); @@ -302,7 +380,7 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, else rh->fEmptyStrips->Fill(eta,phi); h->Fill(eta,phi,n); - rh->fDensity->Fill(eta,phi,n); + if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n); } // for t } // for s @@ -311,20 +389,32 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) { for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { Double_t eLossV = h->GetBinContent(ieta, iphi); - Float_t eta = h->GetXaxis()->GetBinCenter(ieta); - Float_t phi = h->GetYaxis()->GetBinCenter(iphi); + Double_t eta = h->GetXaxis()->GetBinCenter(ieta); + Double_t phi = h->GetYaxis()->GetBinCenter(iphi); Int_t jeta = rh->fEmptyStrips->GetXaxis()->FindBin(eta); Int_t jphi = rh->fEmptyStrips->GetYaxis()->FindBin(phi); Double_t empty = rh->fEmptyStrips->GetBinContent(jeta, jphi); Double_t total = rh->fTotalStrips->GetBinContent(jeta, jphi); Double_t hits = rh->fBasicHits->GetBinContent(ieta,iphi); + // Mean in region of interest Double_t poissonM = (total <= 0 || empty <= 0 ? 0 : -TMath::Log(empty / total)); + + // Note, that given filled=total-empty, and + // + // m = -log(empty/total) + // = -log(1 - filled/total) + // + // v = m / (1 - exp(-m)) + // = -total/filled * (log(total-filled)-log(total)) + // = -total / (total-empty) * log(empty/total) + // = total (log(total)-log(empty)) / (total-empty) + // Double_t poissonV = hits; if(poissonM > 0) // Correct for counting statistics and weight by counts - poissonV = (hits * poissonM) / (1 - TMath::Exp(-1*poissonM)); + poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM)); Double_t poissonE = TMath::Sqrt(hits); if(poissonV > 0) poissonE = TMath::Sqrt(poissonV); @@ -333,6 +423,7 @@ AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd, if (fUsePoisson) { h->SetBinContent(ieta,iphi,poissonV); h->SetBinError(ieta,iphi,poissonE); + rh->fDensity->Fill(eta, phi, poissonV); } } } @@ -383,12 +474,39 @@ AliFMDDensityCalculator::CacheMaxWeights() fFMD3iMax.Set(nEta); fFMD3oMax.Set(nEta); + fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5); + fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i"); + fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i"); + fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o"); + fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i"); + fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o"); + + AliInfo(Form("Get eta axis with %d bins from %f to %f", + nEta, eta.GetXmin(), eta.GetXmax())); + fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5); + fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i"); + fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i"); + fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o"); + fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i"); + fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o"); + for (Int_t i = 0; i < nEta; i++) { - fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1); - fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1); - fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1); - fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1); - fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1); + Double_t w[5]; + w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1); + w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1); + w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1); + w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1); + w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1); + Double_t l[5]; + l[0] = GetMultCut(1, 'I', i+1, false); + l[1] = GetMultCut(2, 'I', i+1, false); + l[2] = GetMultCut(2, 'O', i+1, false); + l[3] = GetMultCut(3, 'I', i+1, false); + l[4] = GetMultCut(3, 'O', i+1, false); + for (Int_t j = 0; j < 5; j++) { + if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]); + if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]); + } } } @@ -481,7 +599,7 @@ AliFMDDensityCalculator::NParticles(Float_t mult, // Return: // The number of particles // - if (mult <= GetMultCut()) return 0; + // if (mult <= GetMultCut()) return 0; if (lowFlux) return 1; AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); @@ -540,7 +658,8 @@ AliFMDDensityCalculator::Correction(UShort_t d, // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); - Float_t correction = AcceptanceCorrection(r,t); + Float_t correction = 1; + if (fUsePhiAcceptance) correction = AcceptanceCorrection(r,t); if (lowFlux) { TH1D* dblHitCor = 0; if (fcm.GetDoubleHit()) @@ -670,8 +789,18 @@ AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents) TIter next(&fRingHistos); RingHistos* o = 0; - while ((o = static_cast(next()))) + THStack* sums = new THStack("sums", "sums of ring signals"); + while ((o = static_cast(next()))) { o->ScaleHistograms(d, nEvents); + TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, + o->fDensity->GetNbinsY(),"e"); + sum->Scale(1., "width"); + sum->SetTitle(o->GetName()); + sum->SetDirectory(0); + sum->SetYTitle("#sum N_{ch,incl}"); + sums->Add(sum); + } + d->Add(sums); } //____________________________________________________________________ @@ -685,6 +814,7 @@ AliFMDDensityCalculator::DefineOutput(TList* dir) // dir List to write in // TList* d = new TList; + d->SetOwner(); d->SetName(GetName()); dir->Add(d); d->Add(fWeightedSum); @@ -692,6 +822,8 @@ AliFMDDensityCalculator::DefineOutput(TList* dir) d->Add(fCorrections); d->Add(fAccI); d->Add(fAccO); + d->Add(fMaxWeights); + d->Add(fLowCuts); TIter next(&fRingHistos); RingHistos* o = 0; @@ -713,10 +845,19 @@ AliFMDDensityCalculator::Print(Option_t* option) const for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; ind[gROOT->GetDirLevel()] = '\0'; std::cout << ind << ClassName() << ": " << GetName() << '\n' + << std::boolalpha << ind << " Multiplicity cut: " << fMultCut << '\n' + << ind << " # of (xi+sigma) factor: " << fNXi << '\n' + << ind << " Include sigma in cut: " << fIncludeSigma << '\n' + << ind << " Low cut method: " + << (fMultCut > 0 ? "fixed" : + (fNXi >= 0 ? "xi+sigma" : "fit range")) << '\n' << ind << " Max(particles): " << fMaxParticles << '\n' + << ind << " Poisson method: " << fUsePoisson << '\n' + << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n' << ind << " Eta lumping: " << fEtaLumping << '\n' << ind << " Phi lumping: " << fPhiLumping << '\n' + << std::noboolalpha << std::flush; TString opt(option); opt.ToLower(); @@ -791,70 +932,61 @@ AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r) // d detector // r ring // - fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()), - Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive " - "N_{ch} in %s", fName.Data()), - 2500, -.5, 24.5, 2500, -.5, 24.5); - fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()), - Form("#Delta E/#Delta E_{mip} vs corrected inclusive " - "N_{ch} in %s", fName.Data()), - 2500, -.5, 24.5, 2500, -.5, 24.5); + fEvsN = new TH2D("elossVsNnocorr", + "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}", + 250, -.5, 24.5, 250, -.5, 24.5); fEvsN->SetXTitle("#Delta E/#Delta E_{mip}"); fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)"); fEvsN->Sumw2(); fEvsN->SetDirectory(0); - fEvsM->SetXTitle("#Delta E/#Delta E_{mip}"); - fEvsM->SetYTitle("Inclusive N_{ch} (corrected)"); - fEvsM->Sumw2(); + + fEvsM = static_cast(fEvsN->Clone("elossVsNcorr")); + fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}"); fEvsM->SetDirectory(0); - fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()), - Form("Average inclusive N_{ch} vs #eta (uncorrected) " - "in %s", fName.Data()), 200, -4, 6); - fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()), - Form("Average inclusive N_{ch} vs #eta (corrected) " - "in %s", fName.Data()), 200, -4, 6); + fEtaVsN = new TProfile("etaVsNnocorr", + "Average inclusive N_{ch} vs #eta (uncorrected)", + 200, -4, 6); fEtaVsN->SetXTitle("#eta"); fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)"); fEtaVsN->SetDirectory(0); fEtaVsN->SetLineColor(Color()); fEtaVsN->SetFillColor(Color()); - fEtaVsM->SetXTitle("#eta"); + + fEtaVsM = static_cast(fEtaVsN->Clone("etaVsNcorr")); + fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)"); fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)"); fEtaVsM->SetDirectory(0); - fEtaVsM->SetLineColor(Color()); - fEtaVsM->SetFillColor(Color()); - fCorr = new TProfile(Form("%s_corr", fName.Data()), - Form("Average correction in %s", fName.Data()), - 200, -4, 6); + fCorr = new TProfile("corr", "Average correction", 200, -4, 6); fCorr->SetXTitle("#eta"); fCorr->SetYTitle("#LT correction#GT"); fCorr->SetDirectory(0); fCorr->SetLineColor(Color()); fCorr->SetFillColor(Color()); - fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()), - Form("Inclusive N_{ch} density in %s", fName.Data()), + fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density", 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi()); fDensity->SetDirectory(0); + fDensity->Sumw2(); + fDensity->SetMarkerColor(Color()); fDensity->SetXTitle("#eta"); fDensity->SetYTitle("#phi [radians]"); fDensity->SetZTitle("Inclusive N_{ch} density"); - fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()), - Form("N_{ch} from energy loss vs from Poission %s", - fName.Data()), 100, 0, 20, 100, 0, 20); + fELossVsPoisson = new TH2D("elossVsPoisson", + "N_{ch} from energy loss vs from Poission", + 100, 0, 20, 100, 0, 20); fELossVsPoisson->SetDirectory(0); fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE"); fELossVsPoisson->SetYTitle("N_{ch} from Poisson"); fELossVsPoisson->SetZTitle("Correlation"); - fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data()), - Form("# of empty strips vs. total @ # strips in %s", - fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5); + fEmptyVsTotal = new TH2D("emptyVsTotal", + "# of empty strips vs. total # strips", + 21, -.5, 20.5, 21, -0.5, 20.5); fEmptyVsTotal->SetDirectory(0); fEmptyVsTotal->SetXTitle("Total # strips"); fEmptyVsTotal->SetYTitle("Empty # strips"); @@ -1030,7 +1162,7 @@ AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents) TList* l = GetOutputList(dir); if (!l) return; - TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data())); + TH1* density = GetOutputHist(l,"inclDensity"); if (density) density->Scale(1./nEvents); } diff --git a/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h b/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h index fc5391cc151..80309cdb18a 100644 --- a/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h +++ b/PWG2/FORWARD/analysis2/AliFMDDensityCalculator.h @@ -132,13 +132,33 @@ public: * number of particles that has hit within a region. */ void SetUsePoisson(Bool_t u) { fUsePoisson = u; } + /** + * Set whether to use the phi acceptance correction. + * + * @param u If true, use the phi acceptance (default is false) + */ + void SetUsePhiAcceptance(Bool_t u) { fUsePhiAcceptance = u; } /** * Set the lower multiplicity cut. This overrides the setting in * the energy loss fits. * * @param cut Cut to use */ - void SetMultCut(Double_t cut) { fMultCut = cut; } + void SetMultCut(Double_t cut) { SetMultCuts(cut,cut,cut,cut,cut); } + /** + * Set the lower multiplicity cuts + * + * @param fmd1i Lower mulitplicyt cut for FMD1i + * @param fmd2i Lower mulitplicyt cut for FMD2i + * @param fmd2o Lower mulitplicyt cut for FMD2o + * @param fmd3i Lower mulitplicyt cut for FMD3i + * @param fmd3o Lower mulitplicyt cut for FMD3o + */ + void SetMultCuts(Double_t fmd1i, + Double_t fmd2i, + Double_t fmd2o, + Double_t fmd3i, + Double_t fmd3o); /** * Set the luming factors used in the Poisson method * @@ -149,6 +169,29 @@ public: fEtaLumping = (eta < 1 ? 1 : eta); fPhiLumping = (phi < 1 ? 1 : phi); } + /** + * Set the number of landau width to subtract from the most probably + * value to get the low cut. + * + * @param n Number of @f$ \xi@f$ + */ + void SetNXi(Double_t nXi) { fNXi = nXi; } + /** + * Whether to include sigma in the number subtracted from the MPV to + * get the low cut + * + * @param u If true, then low cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$ + */ + void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; } + /** + * Get the multiplicity cut. If the user has set fMultCut (via + * SetMultCut) then that value is used. If not, then the lower + * value of the fit range for the enery loss fits is returned. + * + * @return Lower cut on multiplicity + */ + Double_t GetMultCut(UShort_t d, Char_t r, Double_t eta, + Bool_t errors=true) const; /** * Get the multiplicity cut. If the user has set fMultCut (via * SetMultCut) then that value is used. If not, then the lower @@ -156,7 +199,8 @@ public: * * @return Lower cut on multiplicity */ - Double_t GetMultCut() const; + Double_t GetMultCut(UShort_t d, Char_t r, Int_t ieta, + Bool_t errors=true) const; /** * Print information * @@ -341,11 +385,14 @@ protected: TList fRingHistos; // List of histogram containers Double_t fMultCut; // Low cut on scaled energy loss + Double_t fNXi; // Delta-n(xi+sigma) factor + Bool_t fIncludeSigma; // Whether or not to include sigma in cut TH1D* fSumOfWeights; // Histogram TH1D* fWeightedSum; // Histogram TH1D* fCorrections; // Histogram UShort_t fMaxParticles; // Maximum particle weight to use Bool_t fUsePoisson; // If true, then use poisson statistics + Bool_t fUsePhiAcceptance; // Whether to correct for corners TH1D* fAccI; // Acceptance correction for inner rings TH1D* fAccO; // Acceptance correction for outer rings TArrayI fFMD1iMax; // Array of max weights @@ -353,12 +400,14 @@ protected: TArrayI fFMD2oMax; // Array of max weights TArrayI fFMD3iMax; // Array of max weights TArrayI fFMD3oMax; // Array of max weights + TH2D* fMaxWeights; // Histogram of max weights + TH2D* fLowCuts; // Histogram of low cuts Int_t fEtaLumping; // How to lump eta bins for Poisson Int_t fPhiLumping; // How to lump phi bins for Poisson Int_t fDebug; // Debug level + Double_t fMultCuts[5]; // Per-ring multiplicity cuts - - ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density + ClassDef(AliFMDDensityCalculator,3); // Calculate Nch density }; #endif diff --git a/PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx b/PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx index fb111a97694..02f19642263 100644 --- a/PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx @@ -187,7 +187,8 @@ AliFMDEnergyFitterTask::UserExec(Option_t*) AliMCEvent* mcevent = MCEvent(); if(mcevent) { AliHeader* header = mcevent->Header(); - AliGenHijingEventHeader* hijingHeader = dynamic_cast(header->GenEventHeader()); + AliGenHijingEventHeader* hijingHeader = + dynamic_cast(header->GenEventHeader()); if(hijingHeader) { Float_t b = hijingHeader->ImpactParameter(); if(bfbHigh) return; @@ -229,13 +230,14 @@ AliFMDEnergyFitterTask::UserExec(Option_t*) InitializeSubs(); } - Bool_t lowFlux = kFALSE; - UInt_t triggers = 0; - UShort_t ivz = 0; - Double_t vz = 0; - Double_t cent = 0; - UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent); + Bool_t lowFlux = kFALSE; + UInt_t triggers = 0; + UShort_t ivz = 0; + Double_t vz = 0; + Double_t cent = 0; + UShort_t nClusters = 0; + UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, + ivz, vz, cent, nClusters); if (found & AliFMDEventInspector::kNoEvent) return; if (found & AliFMDEventInspector::kNoTriggers) return; if (found & AliFMDEventInspector::kNoSPD) return; @@ -277,20 +279,6 @@ AliFMDEnergyFitterTask::Terminate(Option_t*) return; } -#if 0 - // Get our histograms from the container - TH1I* hEventsTr = 0;//static_cast(list->FindObject("nEventsTr")); - TH1I* hEventsTrVtx = 0;//static_cast(list->FindObject("nEventsTrVtx")); - TH1I* hTriggers = 0; - if (!fEventInspector.FetchHistograms(list, hEventsTr, - hEventsTrVtx, hTriggers)) { - AliError(Form("Didn't get histograms from event selector " - "(hEventsTr=%p,hEventsTrVtx=%p)", - hEventsTr, hEventsTrVtx)); - list->ls(); - return; - } -#endif AliInfo("Fitting energy loss spectra"); fEnergyFitter.Fit(list); @@ -298,17 +286,6 @@ AliFMDEnergyFitterTask::Terminate(Option_t*) TList* list2 = static_cast(list->Clone(Form("%sResults", list->GetName()))); PostData(2, list2); -#if 0 - // Temporary code to save output to a file - should be disabled once - // the plugin stuff works - list->ls(); - TDirectory* savdir = gDirectory; - TFile* tmp = TFile::Open("elossfits.root", "RECREATE"); - list->Write(list->GetName(), TObject::kSingleKey); - tmp->Write(); - tmp->Close(); - savdir->cd(); -#endif } //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx b/PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx index 3b5153827d3..ece4bb3f1be 100644 --- a/PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx @@ -37,6 +37,7 @@ AliFMDEventInspector::AliFMDEventInspector() : TNamed(), fHEventsTr(0), fHEventsTrVtx(0), + fHEventsAccepted(0), fHTriggers(0), fHType(0), fHWords(0), @@ -49,7 +50,8 @@ AliFMDEventInspector::AliFMDEventInspector() fField(999), fCollisionSystem(kUnknown), fDebug(0), - fCentAxis(0) + fCentAxis(0), + fVtxAxis(10,-10,10) { // // Constructor @@ -61,6 +63,7 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name) : TNamed("fmdEventInspector", name), fHEventsTr(0), fHEventsTrVtx(0), + fHEventsAccepted(0), fHTriggers(0), fHType(0), fHWords(0), @@ -73,7 +76,8 @@ AliFMDEventInspector::AliFMDEventInspector(const char* name) fField(999), fCollisionSystem(kUnknown), fDebug(0), - fCentAxis(0) + fCentAxis(0), + fVtxAxis(10,-10,10) { // // Constructor @@ -88,6 +92,7 @@ AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o) : TNamed(o), fHEventsTr(o.fHEventsTr), fHEventsTrVtx(o.fHEventsTrVtx), + fHEventsAccepted(o.fHEventsAccepted), fHTriggers(o.fHTriggers), fHType(o.fHType), fHWords(o.fHWords), @@ -100,7 +105,8 @@ AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o) fField(o.fField), fCollisionSystem(o.fCollisionSystem), fDebug(0), - fCentAxis(0) + fCentAxis(0), + fVtxAxis(o.fVtxAxis) { // // Copy constructor @@ -134,6 +140,7 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o) TNamed::operator=(o); fHEventsTr = o.fHEventsTr; fHEventsTrVtx = o.fHEventsTrVtx; + fHEventsAccepted = o.fHEventsAccepted; fHTriggers = o.fHTriggers; fHType = o.fHType; fHWords = o.fHWords; @@ -146,6 +153,8 @@ AliFMDEventInspector::operator=(const AliFMDEventInspector& o) fEnergy = o.fEnergy; fField = o.fField; fCollisionSystem = o.fCollisionSystem; + fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), + o.fVtxAxis.GetXmax()); if (fList) { fList->SetName(GetName()); if (fHEventsTr) fList->Add(fHEventsTr); @@ -206,12 +215,14 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis) // ----- 92 number --------- ---- 1 --- TArrayD limits(93); for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i; + + fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax()); fCentAxis = new TAxis(limits.GetSize()-1, limits.GetArray()); fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", - vtxAxis.GetNbins(), - vtxAxis.GetXmin(), - vtxAxis.GetXmax()); + 4*vtxAxis.GetNbins(), + 2*vtxAxis.GetXmin(), + 2*vtxAxis.GetXmax()); fHEventsTr->SetXTitle("v_{z} [cm]"); fHEventsTr->SetYTitle("# of events"); fHEventsTr->SetFillColor(kRed+1); @@ -220,19 +231,26 @@ AliFMDEventInspector::Init(const TAxis& vtxAxis) // fHEventsTr->Sumw2(); fList->Add(fHEventsTr); - fHEventsTrVtx = new TH1I("nEventsTrVtx", - "Number of events w/trigger and vertex", - vtxAxis.GetNbins(), - vtxAxis.GetXmin(), - vtxAxis.GetXmax()); - fHEventsTrVtx->SetXTitle("v_{z} [cm]"); - fHEventsTrVtx->SetYTitle("# of events"); + fHEventsTrVtx = static_cast(fHEventsTr->Clone("nEventsTrVtx")); + fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex"); fHEventsTrVtx->SetFillColor(kBlue+1); - fHEventsTrVtx->SetFillStyle(3001); fHEventsTrVtx->SetDirectory(0); // fHEventsTrVtx->Sumw2(); fList->Add(fHEventsTrVtx); + fHEventsAccepted = new TH1I("nEventsAccepted", + "Number of events w/trigger and vertex in range", + 2*vtxAxis.GetNbins(), + 2*vtxAxis.GetXmin(), + 2*vtxAxis.GetXmax()); + fHEventsAccepted->SetXTitle("v_{z} [cm]"); + fHEventsAccepted->SetYTitle("# of events"); + fHEventsAccepted->SetFillColor(kGreen+1); + fHEventsAccepted->SetFillStyle(3001); + fHEventsAccepted->SetDirectory(0); + // fHEventsAccepted->Sumw2(); + fList->Add(fHEventsAccepted); + fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1); fHTriggers->SetFillColor(kRed+1); @@ -316,7 +334,8 @@ AliFMDEventInspector::Process(const AliESDEvent* event, Bool_t& lowFlux, UShort_t& ivz, Double_t& vz, - Double_t& cent) + Double_t& cent, + UShort_t& nClusters) { // // Process the event @@ -342,7 +361,7 @@ AliFMDEventInspector::Process(const AliESDEvent* event, } // --- Read trigger information from the ESD and store in AOD object - if (!ReadTriggers(event, triggers)) { + if (!ReadTriggers(event, triggers, nClusters)) { if (fDebug > 2) { AliWarning("Failed to read triggers from ESD"); } return kNoTriggers; @@ -386,15 +405,16 @@ AliFMDEventInspector::Process(const AliESDEvent* event, fHEventsTrVtx->Fill(vz); // --- Get the vertex bin ------------------------------------------ - ivz = fHEventsTr->GetXaxis()->FindBin(vz); - if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { + ivz = fVtxAxis.FindBin(vz); + if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) { if (fDebug > 3) { AliWarning(Form("Vertex @ %f outside of range [%f,%f]", - vz, fHEventsTr->GetXaxis()->GetXmin(), - fHEventsTr->GetXaxis()->GetXmax())); } + vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); + } ivz = 0; return kBadVertex; } + fHEventsAccepted->Fill(vz); // --- Check the FMD ESD data -------------------------------------- if (!event->GetFMDData()) { @@ -439,7 +459,8 @@ AliFMDEventInspector::ReadCentrality(const AliESDEvent* esd, //____________________________________________________________________ Bool_t -AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers) +AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers, + UShort_t& nClusters) { // // Read the trigger information from the ESD event @@ -475,6 +496,7 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers) // then the AliPhysicsSelection object would overcount by a // factor of 2! :-( Bool_t offline = ih->IsEventSelected(); + nClusters = 0; if (offline) { triggers |= AliAODForwardMult::kOffline; triggers |= AliAODForwardMult::kInel; @@ -489,14 +511,22 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers) // Check if we have one or more tracklets // in the range -1 < eta < 1 to set the INEL>0 // trigger flag. + // + // Also count tracklets as a single cluster Int_t n = spdmult->GetNumberOfTracklets(); for (Int_t j = 0; j < n; j++) { if(TMath::Abs(spdmult->GetEta(j)) < 1) { triggers |= AliAODForwardMult::kInelGt0; - break; + nClusters++; } } + n = spdmult->GetNumberOfSingleClusters(); + for (Int_t j = 0; j < n; j++) { + Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)); + if (TMath::Abs(eta) < 1) nClusters++; + } } + if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0; } // Analyse some trigger stuff @@ -505,7 +535,8 @@ AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers) triggers |= AliAODForwardMult::kNSD; - //Check pileup + // Check for multiple vertices (pile-up) with at least 3 + // contributors and at least 0.8cm from the primary vertex Bool_t pileup = esd->IsPileupFromSPD(3,0.8); if (pileup) { triggers |= AliAODForwardMult::kPileUp; @@ -628,6 +659,40 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz) // @c true on success, @c false otherwise // vz = 0; +#if 1 + // This is the code used by the 1st physics people + const AliESDVertex* vertex = esd->GetPrimaryVertex(); + if (!vertex || !vertex->GetStatus()) { + if (fDebug > 2) { + AliWarning(Form("No primary vertex (%p) or bad status %d", + vertex, (vertex ? vertex->GetStatus() : -1))); + } + return false; + } + const AliESDVertex* vertexSPD = esd->GetPrimaryVertexSPD(); + if (!vertexSPD || !vertexSPD->GetStatus()) { + if (fDebug > 2) { + AliWarning(Form("No primary SPD vertex (%p) or bad status %d", + vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1))); + } + return false; + } + + // if vertex is from SPD vertexZ, require more stringent cuts + if (vertex->IsFromVertexerZ()) { + if (vertex->GetDispersion() > fMaxVzErr || + vertex->GetZRes() > 1.25 * fMaxVzErr) { + if (fDebug > 2) { + AliWarning(Form("Dispersion %f > %f or resolution %f > %f", + vertex->GetDispersion(), fMaxVzErr, + vertex->GetZRes(), 1.25 * fMaxVzErr)); + } + return false; + } + } + vz = vertex->GetZ(); + return true; +#else // Get the vertex const AliESDVertex* vertex = esd->GetPrimaryVertexSPD(); if (!vertex) { @@ -655,6 +720,7 @@ AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz) // Get the z coordiante vz = vertex->GetZ(); return kTRUE; +#endif } //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliFMDEventInspector.h b/PWG2/FORWARD/analysis2/AliFMDEventInspector.h index ec91eb6a7e1..c5262020300 100644 --- a/PWG2/FORWARD/analysis2/AliFMDEventInspector.h +++ b/PWG2/FORWARD/analysis2/AliFMDEventInspector.h @@ -14,6 +14,7 @@ * @ingroup pwg2_forward_aod */ #include +#include class AliESDEvent; class TH2D; class TH1D; @@ -135,6 +136,7 @@ public: * @param vz On return, the z position of the interaction * @param cent On return, the centrality (in percent) or < 0 * if not found + * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$ * * @return 0 (or kOk) on success, otherwise a bit mask of error codes */ @@ -143,7 +145,8 @@ public: Bool_t& lowFlux, UShort_t& ivz, Double_t& vz, - Double_t& cent); + Double_t& cent, + UShort_t& nClusters); /** * Define the output histograms. These are put in a sub list of the * passed list. The histograms are merged before the parent task calls @@ -224,10 +227,12 @@ protected: * * @param esd ESD event * @param triggers On return, contains the trigger bits + * @param nClusters On return, number of SPD clusters in @f$ |\eta|<1@f$ * * @return @c true on success, @c false otherwise */ - Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers); + Bool_t ReadTriggers(const AliESDEvent* esd, UInt_t& triggers, + UShort_t& nClusters); /** * Read the vertex information from the ESD event * @@ -251,6 +256,7 @@ protected: TH1I* fHEventsTr; //! Histogram of events w/trigger TH1I* fHEventsTrVtx; //! Events w/trigger and vertex + TH1I* fHEventsAccepted; //! Events w/trigger and vertex in range TH1I* fHTriggers; //! Triggers TH1I* fHType; //! Type (low/high flux) of event TH1I* fHWords; //! Trigger words @@ -264,6 +270,7 @@ protected: UShort_t fCollisionSystem; // Collision system Int_t fDebug; // Debug level TAxis* fCentAxis; // Centrality axis used in histograms + TAxis fVtxAxis; ClassDef(AliFMDEventInspector,3); // Inspect the event }; diff --git a/PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx b/PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx index b4c4d50a470..f3b4d4ebc1c 100644 --- a/PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx @@ -33,6 +33,50 @@ ClassImp(AliFMDHistCollector) ; // For Emacs #endif +//____________________________________________________________________ +AliFMDHistCollector::AliFMDHistCollector() + : fNCutBins(0), + fCorrectionCut(0), + fFirstBins(), + fLastBins(), + fDebug(0), + fList(0), + fSumRings(0), + fCoverage(0), + fMergeMethod(kStraightMean), + fFiducialMethod(kByCut) +{} + +//____________________________________________________________________ +AliFMDHistCollector::AliFMDHistCollector(const char* title) + : TNamed("fmdHistCollector", title), + fNCutBins(2), + fCorrectionCut(0.5), + fFirstBins(1), + fLastBins(1), + fDebug(0), + fList(0), + fSumRings(0), + fCoverage(0), + fMergeMethod(kStraightMean), + fFiducialMethod(kByCut) +{ +} +//____________________________________________________________________ +AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o) + : TNamed(o), + fNCutBins(o.fNCutBins), + fCorrectionCut(o.fCorrectionCut), + fFirstBins(o.fFirstBins), + fLastBins(o.fLastBins), + fDebug(o.fDebug), + fList(o.fList), + fSumRings(o.fSumRings), + fCoverage(o.fCoverage), + fMergeMethod(o.fMergeMethod), + fFiducialMethod(o.fFiducialMethod) +{} + //____________________________________________________________________ AliFMDHistCollector::~AliFMDHistCollector() { @@ -63,6 +107,7 @@ AliFMDHistCollector::operator=(const AliFMDHistCollector& o) fCoverage = o.fCoverage; fMergeMethod = o.fMergeMethod; fFiducialMethod = o.fFiducialMethod; + return *this; } @@ -151,7 +196,7 @@ AliFMDHistCollector::Init(const TAxis& vtxAxis, // Store the result for later use fFirstBins[(iVz-1)*5+iIdx] = first; fLastBins[(iVz-1)*5+iIdx] = last; - TH2D* obg = static_cast(bg->Clone()); + TH2D* obg = static_cast(bg->Clone(Form("secMapFMD%d%c", d, r))); obg->SetDirectory(0); obg->Reset(); vtxList->Add(obg); @@ -212,8 +257,10 @@ AliFMDHistCollector::DefineOutput(TList* dir) fList->SetOwner(); fList->SetName(GetName()); dir->Add(fList); + } + //____________________________________________________________________ Int_t AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const @@ -477,7 +524,8 @@ AliFMDHistCollector::MergeBins(Double_t c, Double_t e, //____________________________________________________________________ Bool_t -AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists, +AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists, + AliForwardUtil::Histos& sums, UShort_t vtxbin, TH2D& out) { @@ -499,13 +547,32 @@ AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists, TH2D* h = hists.Get(d,r); TH2D* t = static_cast(h->Clone(Form("FMD%d%c_tmp",d,r))); Int_t i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1)); + TH2D* o = sums.Get(d, r); - // Outer rings have better phi segmentation - rebin to same as inner. - if (q == 1) t->RebinY(2); - + // Get valid range Int_t first = 0; Int_t last = 0; GetFirstAndLast(d, r, vtxbin, first, last); + + // Zero outside valid range + for (Int_t iPhi = 0; iPhi <= t->GetNbinsY()+1; iPhi++) { + // Lower range + for (Int_t iEta = 1; iEta < first; iEta++) { + t->SetBinContent(iEta,iPhi,0); + t->SetBinError(iEta,iPhi,0); + } + for (Int_t iEta = last+1; iEta <= t->GetNbinsX(); iEta++) { + t->SetBinContent(iEta,iPhi,0); + t->SetBinError(iEta,iPhi,0); + } + } + for (Int_t iEta = first; iEta <= last; iEta++) + t->SetBinContent(iEta,0,1); + // Add to our per-ring sum + o->Add(t); + + // Outer rings have better phi segmentation - rebin to same as inner. + if (q == 1) t->RebinY(2); // Now update profile output for (Int_t iEta = first; iEta <= last; iEta++) { @@ -518,7 +585,8 @@ AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists, Float_t noc = overlap >= 0? 0.5 : 1; out.SetBinContent(iEta, 0, ooc + noc); - for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) { + // Should we loop over h or t Y bins - I think it's t + for (Int_t iPhi = 1; iPhi <= t->GetNbinsY(); iPhi++) { Double_t c = t->GetBinContent(iEta,iPhi); Double_t e = t->GetBinError(iEta,iPhi); Double_t ee = t->GetXaxis()->GetBinCenter(iEta); @@ -577,22 +645,26 @@ AliFMDHistCollector::Print(Option_t* /* option */) const case kLeastError: std::cout << "least error\n"; break; } - std::cout << ind << " Bin ranges:\n" << ind << " v_z bin"; + std::cout << ind << " Bin ranges:\n" << ind << " rings "; Int_t nVz = fFirstBins.fN / 5; - for (UShort_t iVz = 1; iVz <= nVz; iVz++) - std::cout << " | " << std::setw(7) << iVz; - std::cout << '\n' << ind << " / ring "; - for (UShort_t iVz = 1; iVz <= nVz; iVz++) - std::cout << "-+--------"; - std::cout << std::endl; - for (Int_t iIdx = 0; iIdx < 5; iIdx++) { UShort_t d = 0; Char_t r = 0; GetDetRing(iIdx, d, r); + std::cout << ind << " | FMD" << d << r << " "; + } + std::cout << '\n' << ind << " /vz_bin "; + for (Int_t iIdx = 0; iIdx < 5; iIdx++) + std::cout << "-+--------"; + std::cout << std::endl; + + for (UShort_t iVz = 1; iVz <= nVz; iVz++) { + std::cout << " " << std::setw(7) << iVz << " "; + for (Int_t iIdx = 0; iIdx < 5; iIdx++) { + UShort_t d = 0; + Char_t r = 0; + GetDetRing(iIdx, d, r); - std::cout << ind << " FMD" << d << r << " "; - for (UShort_t iVz = 1; iVz <= nVz; iVz++) { Int_t first, last; GetFirstAndLast(iIdx, iVz, first, last); std::cout << " | " << std::setw(3) << first << "-" diff --git a/PWG2/FORWARD/analysis2/AliFMDHistCollector.h b/PWG2/FORWARD/analysis2/AliFMDHistCollector.h index 56634c88e22..3b0e0898cd9 100644 --- a/PWG2/FORWARD/analysis2/AliFMDHistCollector.h +++ b/PWG2/FORWARD/analysis2/AliFMDHistCollector.h @@ -97,54 +97,19 @@ public: /** * Constructor */ - AliFMDHistCollector() - : fNCutBins(0), - fCorrectionCut(0), - fFirstBins(), - fLastBins(), - fDebug(0), - fList(0), - fSumRings(0), - fCoverage(0), - fMergeMethod(kStraightMean), - fFiducialMethod(kByCut) - {} + AliFMDHistCollector(); /** * Constructor * * @param title Name of object */ - AliFMDHistCollector(const char* title) - : TNamed("fmdHistCollector", title), - fNCutBins(2), - fCorrectionCut(0.5), - fFirstBins(1), - fLastBins(1), - fDebug(0), - fList(0), - fSumRings(0), - fCoverage(0), - fMergeMethod(kStraightMean), - fFiducialMethod(kByCut) - {} + AliFMDHistCollector(const char* title); /** * Copy constructor * * @param o Object to copy from */ - AliFMDHistCollector(const AliFMDHistCollector& o) - : TNamed(o), - fNCutBins(o.fNCutBins), - fCorrectionCut(o.fCorrectionCut), - fFirstBins(o.fFirstBins), - fLastBins(o.fLastBins), - fDebug(o.fDebug), - fList(o.fList), - fSumRings(o.fSumRings), - fCoverage(o.fCoverage), - fMergeMethod(o.fMergeMethod), - fFiducialMethod(o.fFiducialMethod) - {} + AliFMDHistCollector(const AliFMDHistCollector& o); /** * Destructor @@ -170,13 +135,16 @@ public: * Do the calculations * * @param hists Cache of histograms + * @param sums Cache to sum ring histograms in * @param vtxBin Vertex bin (1 based) * @param out Output histogram * * @return true on successs */ - virtual Bool_t Collect(AliForwardUtil::Histos& hists, UShort_t vtxBin, - TH2D& out); + virtual Bool_t Collect(const AliForwardUtil::Histos& hists, + AliForwardUtil::Histos& sums, + UShort_t vtxBin, + TH2D& out); /** * Output diagnostic histograms to directory * @@ -368,6 +336,7 @@ protected: void MergeBins(Double_t c, Double_t e, Double_t oc, Double_t oe, Double_t& rc, Double_t& re) const; + Int_t fNCutBins; // Number of additional bins to cut away Float_t fCorrectionCut; // Cut-off on secondary corrections @@ -379,6 +348,7 @@ protected: TH2D* fCoverage; // Sum per ring (on y-axis) MergeMethod fMergeMethod; // Merge methiod for overlapping bins FiducialMethod fFiducialMethod; // Fidicual method + ClassDef(AliFMDHistCollector,1); // Calculate Nch density }; diff --git a/PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx b/PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx index d56a10bffb4..e12ec02d197 100644 --- a/PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDMCCorrector.cxx @@ -61,7 +61,7 @@ AliFMDMCCorrector::operator=(const AliFMDMCCorrector& o) // Reference to this object // AliFMDCorrector::operator=(o); - + fSecondaryForMC = o.fSecondaryForMC; return *this; } @@ -80,6 +80,9 @@ AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists, // Return: // true on successs // + if ((!fUseSecondaryMap || !fSecondaryForMC) && fUseVertexBias) + return kTRUE; + AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); UShort_t uvb = vtxbin; @@ -90,7 +93,7 @@ AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists, TH2D* h = hists.Get(d,r); - if (fUseSecondaryMap) { + if (fUseSecondaryMap && fSecondaryForMC) { TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb); if (!bg) { AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d", @@ -247,6 +250,23 @@ AliFMDMCCorrector::DefineOutput(TList* dir) fComps->SetName("esd_mc_comparison"); d->Add(fComps); } +//____________________________________________________________________ +void +AliFMDMCCorrector::Print(Option_t* option) const +{ + // + // Print information + // Parameters: + // option Not used + // + char ind[gROOT->GetDirLevel()+1]; + for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; + ind[gROOT->GetDirLevel()] = '\0'; + AliFMDCorrector::Print(option); + std::cout << std::boolalpha + << ind << " Use sec. map on MC: " << fSecondaryForMC + << std::noboolalpha << std::endl; +} //____________________________________________________________________ // diff --git a/PWG2/FORWARD/analysis2/AliFMDMCCorrector.h b/PWG2/FORWARD/analysis2/AliFMDMCCorrector.h index a0b2d8e9194..7954b452d87 100644 --- a/PWG2/FORWARD/analysis2/AliFMDMCCorrector.h +++ b/PWG2/FORWARD/analysis2/AliFMDMCCorrector.h @@ -58,7 +58,8 @@ public: fFMD2o(0), fFMD3i(0), fFMD3o(0), - fComps(0) + fComps(0), + fSecondaryForMC(true) {} /** * Constructor @@ -72,7 +73,8 @@ public: fFMD2o(0), fFMD3i(0), fFMD3o(0), - fComps(0) + fComps(0), + fSecondaryForMC(true) {} /** * Copy constructor @@ -86,7 +88,8 @@ public: fFMD2o(o.fFMD2o), fFMD3i(o.fFMD3i), fFMD3o(o.fFMD3o), - fComps(0) + fComps(0), + fSecondaryForMC(o.fSecondaryForMC) {} /** * Destructor @@ -100,6 +103,12 @@ public: * @return Reference to this object */ AliFMDMCCorrector& operator=(const AliFMDMCCorrector&); + /** + * If set, then do not do the secondary correction for MC data + * + * @param use + */ + void SetSecondaryForMC(Bool_t use) { fSecondaryForMC = use; } /** * Initialize this object * @@ -133,6 +142,13 @@ public: * @param dir List to write in */ void DefineOutput(TList* dir); + + /** + * Print information + * + * @param option Not used + */ + void Print(Option_t* option="") const; protected: /** * MAke comparison profiles @@ -160,7 +176,8 @@ protected: TProfile2D* fFMD3i; // Comparison TProfile2D* fFMD3o; // Comparison TList* fComps; // List of comparisons - + Bool_t fSecondaryForMC; // Whether to correct MC data + ClassDef(AliFMDMCCorrector,1); // Calculate Nch density }; diff --git a/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx b/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx index edaf4dd1ee2..e933a92d9eb 100644 --- a/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDMCEventInspector.cxx @@ -456,7 +456,7 @@ AliFMDMCEventInspector::IsSingleDiffractive(AliStack* stack, else return false; - // Invariant masses + // Rapidity shift Double_t m02s = 1 - 2 * p1->Energy() / fEnergy; Double_t m12s = 1 - 2 * p2->Energy() / fEnergy; diff --git a/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx b/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx index db3c03b35e5..cc2e86664b7 100644 --- a/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx @@ -10,7 +10,7 @@ // - AliESDFMD object - copy of input, but with signals merged // // Corrections used: -// - None +// - ELoss fits // // Histograms: // - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of @@ -28,6 +28,7 @@ #include #include #include "AliFMDStripIndex.h" +#include "AliFMDFloatMap.h" #include #include #include @@ -47,7 +48,10 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title) fFMD2o(0), fFMD3i(0), fFMD3o(0), - fSumEta(0) + fSumEta(0), + fOperComp(0), + fThetaVsNr(0), + fOnlyPrimary(false) { // // Constructor @@ -84,7 +88,25 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title) fSumEta->SetMarkerStyle(22); fSumEta->SetFillColor(0); fSumEta->SetFillStyle(0); + + fOper = new AliFMDFloatMap(0,0,0,0); + fOperComp = new TH2I("operComp", "Operation vs # track refs", + kMergedInto, kNone-.5, kMergedInto+.5, + 20, -.5, 19.5); + fOperComp->SetXTitle("Operation"); + fOperComp->SetYTitle("# of track refs in sector"); + fOperComp->SetZTitle("Observations"); + fOperComp->GetXaxis()->SetBinLabel(kNone, "None"); + fOperComp->GetXaxis()->SetBinLabel(kCandidate, "Candidate"); + fOperComp->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other"); + fOperComp->GetXaxis()->SetBinLabel(kMergedInto, "Merged into"); + fOperComp->SetDirectory(0); + fThetaVsNr = new TH2D("thetaVsNr", "#theta of track vs # track references", + 360, 0, 360, 20, -.5, 19.5); + fThetaVsNr->SetXTitle("#theta [degrees]"); + fThetaVsNr->SetYTitle("# of track references"); + fThetaVsNr->SetDirectory(0); } //____________________________________________________________________ @@ -95,7 +117,10 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o) fFMD2o(o.fFMD2o), fFMD3i(o.fFMD3i), fFMD3o(o.fFMD3o), - fSumEta(o.fSumEta) + fSumEta(o.fSumEta), + fOperComp(o.fOperComp), + fThetaVsNr(o.fThetaVsNr), + fOnlyPrimary(o.fOnlyPrimary) { // // Copy constructor @@ -133,6 +158,7 @@ AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o) // Reference to this // AliFMDSharingFilter::operator=(o); + fOnlyPrimary = o.fOnlyPrimary; return *this; } @@ -142,6 +168,8 @@ AliFMDMCSharingFilter::StoreParticle(UShort_t d, Char_t r, UShort_t s, UShort_t t, + UShort_t nR, + Double_t theta, AliESDFMD& output) const { // @@ -152,10 +180,15 @@ AliFMDMCSharingFilter::StoreParticle(UShort_t d, // r Ring // s Sector // t Strip + // nR Number of references to this particle in this sector // output Output ESD object // Double_t old = output.Multiplicity(d,r,s,t); if (old == AliESDFMD::kInvalidMult) old = 0; + if (fOper) fOperComp->Fill(fOper->operator()(d,r,s,t), nR); + if (theta < 0) theta += 2*TMath::Pi(); + theta *= 180. / TMath::Pi(); + fThetaVsNr->Fill(theta, nR); output.SetMultiplicity(d,r,s,t,old+1); } @@ -183,6 +216,10 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input, // output.Clear(); + // Increase event count - stored in + // underflow bin + fSumEta->AddBinContent(0); + // Copy eta values to output for (UShort_t ed = 1; ed <= 3; ed++) { UShort_t nq = (ed == 1 ? 1 : 2); @@ -210,17 +247,28 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input, if (!isCharged) continue; Bool_t isPrimary = stack->IsPhysicalPrimary(iTr); + // Pseudo rapidity and azimuthal angle + Double_t eta = particle->Eta(); + Double_t phi = particle->Phi(); + // Fill 'dn/deta' histogram if (isPrimary && iTr < nPrim) { - fSumEta->Fill(particle->Eta()); - primary->Fill(particle->Eta(), particle->Phi()); + // Avoid under count - used to store event count + if (eta >= fSumEta->GetXaxis()->GetXmin()) fSumEta->Fill(eta); + primary->Fill(eta, phi); } + // Bail out if we're only processing primaries - perhaps we should + // track back to the original primary? + if (fOnlyPrimary && !isPrimary) continue; + Int_t nTrRef = particle->GetNumberOfTrackReferences(); Int_t longest = -1; Double_t angle = 0; UShort_t oD = 0, oS = 1024, oT = 1024; Char_t oR = '\0'; + UShort_t nC = 0; + Double_t oTheta = 0; for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { AliTrackReference* ref = particle->GetTrackReference(iTrRef); @@ -231,13 +279,24 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input, if (ref->DetectorId() != AliTrackReference::kFMD) continue; + // Count number of track refs in this sector + nC++; + // Get the detector coordinates UShort_t d, s, t; Char_t r; AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t); - if (oD > 0 && oR != '\0' && (d != oD || r != oR)) { + // If this is a new detector/ring, then reset the other one + if (oD > 0 && oR != '\0' && oS != 1024 && + (d != oD || r != oR || s != oS)) { longest = -1; - StoreParticle(oD, oR, oS, oT, output); + angle = 0; + StoreParticle(oD, oR, oS, oT, nC, oTheta, output); + nC = 0; + oD = 0; + oR = '\0'; + oS = 1024; + oT = 1024; } oD = d; @@ -256,6 +315,7 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input, longest = iTrRef; angle = ang; } + oTheta = theta; } // Loop over track references if (longest < 0) continue; @@ -267,7 +327,7 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input, Char_t r; AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t); - StoreParticle(d,r,s,t,output); + StoreParticle(d,r,s,t,nC,particle->Theta(),output); } // Loop over tracks return kTRUE; } @@ -335,6 +395,8 @@ AliFMDMCSharingFilter::DefineOutput(TList* dir) cd->Add(fFMD3i); cd->Add(fFMD3o); dir->Add(fSumEta); + cd->Add(fOperComp); + cd->Add(fThetaVsNr); } //____________________________________________________________________ @@ -354,7 +416,8 @@ AliFMDMCSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents) AliWarning(Form("No mcSumEta histogram found in output list")); return; } - sumEta->Scale(1. / nEvents, "width"); + Double_t n = nEvents; // sumEta->GetBinContent(0); + sumEta->Scale(1. / n, "width"); } //____________________________________________________________________ @@ -371,6 +434,9 @@ AliFMDMCSharingFilter::Print(Option_t* option) const for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; ind[gROOT->GetDirLevel()] = '\0'; AliFMDSharingFilter::Print(option); + std::cout << std::boolalpha + << ind << " Only primary tracks: " << fOnlyPrimary + << std::noboolalpha << std::endl; } //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h b/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h index 894761c60c0..09aa4e0724c 100644 --- a/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h +++ b/PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h @@ -58,7 +58,10 @@ public: fFMD2o(0), fFMD3i(0), fFMD3o(0), - fSumEta(0) + fSumEta(0), + fOperComp(0), + fThetaVsNr(0), + fOnlyPrimary(false) {} /** * Constructor @@ -80,6 +83,13 @@ public: * @return Reference to this */ AliFMDMCSharingFilter& operator=(const AliFMDMCSharingFilter& o); + + /** + * If set, then only process primary tracks + * + * @param use + */ + void SetOnlyPrimary(Bool_t use) { fOnlyPrimary = use; } /** * Filter the input kinematics and track references, using * some of the ESD information @@ -138,14 +148,16 @@ protected: * @param output Output ESD object */ void StoreParticle(UShort_t d, Char_t r, UShort_t s, UShort_t t, - AliESDFMD& output) const; + UShort_t nr, Double_t theta, AliESDFMD& output) const; TH2D* fFMD1i; // ESD-MC correlation TH2D* fFMD2i; // ESD-MC correlation TH2D* fFMD2o; // ESD-MC correlation TH2D* fFMD3i; // ESD-MC correlation TH2D* fFMD3o; // ESD-MC correlation TH1D* fSumEta; // MC dN/deta - + TH2I* fOperComp; // Operation vs # trackrefs + TH2D* fThetaVsNr; // Theta vs # trackrefs + Bool_t fOnlyPrimary; // Only process primary tracks ClassDef(AliFMDMCSharingFilter,1); // }; diff --git a/PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx b/PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx index 71aba79d3eb..bfa24faeef5 100644 --- a/PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx +++ b/PWG2/FORWARD/analysis2/AliFMDSharingFilter.cxx @@ -31,6 +31,7 @@ #include "AliFMDCorrELossFit.h" #include #include +#include #include #include @@ -39,6 +40,12 @@ ClassImp(AliFMDSharingFilter) ; // This is for Emacs #endif +#define DBG(L,M) \ + do { if (L>fDebug)break; std::cout << (M) << std::flush;} while(false) +#define DBGL(L,M) \ + do { if (L>fDebug)break; std::cout << (M) << std::endl;} while(false) + + //____________________________________________________________________ AliFMDSharingFilter::AliFMDSharingFilter() @@ -47,6 +54,10 @@ AliFMDSharingFilter::AliFMDSharingFilter() fLowCut(0.), fCorrectAngles(kFALSE), fNXi(1), + fIncludeSigma(true), + fSummed(0), + fHighCuts(0), + fOper(0), fDebug(0) { // @@ -61,6 +72,10 @@ AliFMDSharingFilter::AliFMDSharingFilter(const char* title) fLowCut(0.), fCorrectAngles(kFALSE), fNXi(1), + fIncludeSigma(true), + fSummed(0), + fHighCuts(0), + fOper(0), fDebug(0) { // @@ -83,6 +98,10 @@ AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o) fLowCut(o.fLowCut), fCorrectAngles(o.fCorrectAngles), fNXi(o.fNXi), + fIncludeSigma(o.fIncludeSigma), + fSummed(o.fSummed), + fHighCuts(o.fHighCuts), + fOper(o.fOper), fDebug(o.fDebug) { // @@ -124,6 +143,10 @@ AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o) fCorrectAngles = o.fCorrectAngles; fNXi = o.fNXi; fDebug = o.fDebug; + fOper = o.fOper; + fSummed = o.fSummed; + fHighCuts = o.fHighCuts; + fIncludeSigma = o.fIncludeSigma; fRingHistos.Delete(); TIter next(&o.fRingHistos); @@ -158,6 +181,43 @@ AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const return static_cast(fRingHistos.At(idx)); } +//____________________________________________________________________ +void +AliFMDSharingFilter::Init() +{ + // Initialise + AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance(); + + // Get the high cut. The high cut is defined as the + // most-probably-value peak found from the energy distributions, minus + // 2 times the width of the corresponding Landau. + AliFMDCorrELossFit* fits = fcm.GetELossFit(); + const TAxis& eAxis = fits->GetEtaAxis(); + + UShort_t nEta = eAxis.GetNbins(); + fHighCuts->SetBins(nEta, eAxis.GetXmin(), eAxis.GetXmax(), 5, .5, 5.5); + fHighCuts->GetYaxis()->SetBinLabel(1, "FMD1i"); + fHighCuts->GetYaxis()->SetBinLabel(2, "FMD2i"); + fHighCuts->GetYaxis()->SetBinLabel(3, "FMD2o"); + fHighCuts->GetYaxis()->SetBinLabel(4, "FMD3i"); + fHighCuts->GetYaxis()->SetBinLabel(5, "FMD3o"); + + UShort_t ybin = 0; + for (UShort_t d = 1; d <= 3; d++) { + UShort_t nr = (d == 1 ? 1 : 2); + for (UShort_t q = 0; q < nr; q++) { + Char_t r = (q == 0 ? 'I' : 'O'); + ybin++; + for (UShort_t e = 1; e <= nEta; e++) { + Double_t eta = eAxis.GetBinCenter(e); + Double_t cut = GetHighCut(d, r, eta, false); + if (cut <= 0) continue; + fHighCuts->SetBinContent(e, ybin, cut); + } + } + } +} + //____________________________________________________________________ Bool_t AliFMDSharingFilter::Filter(const AliESDFMD& input, @@ -178,45 +238,59 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input, output.Clear(); TIter next(&fRingHistos); RingHistos* o = 0; - Double_t lowCut = GetLowCut(); while ((o = static_cast(next()))) o->Clear(); - // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); + if (fOper) fOper->Reset(0); + Int_t nNone = 0; + Int_t nCandidate = 0; + Int_t nMerged = 0; + Int_t nSummed = 0; + Status status[512]; + for(UShort_t d = 1; d <= 3; d++) { Int_t nRings = (d == 1 ? 1 : 2); for (UShort_t q = 0; q < nRings; q++) { - Char_t r = (q == 0 ? 'I' : 'O'); - UShort_t nsec = (q == 0 ? 20 : 40); - UShort_t nstr = (q == 0 ? 512 : 256); - + Char_t r = (q == 0 ? 'I' : 'O'); + UShort_t nsec = (q == 0 ? 20 : 40); + UShort_t nstr = (q == 0 ? 512 : 256); RingHistos* histos = GetRingHistos(d, r); for(UShort_t s = 0; s < nsec; s++) { - Bool_t usedThis = kFALSE; - Bool_t usedPrev = kFALSE; - + for (UShort_t t = 0; t < nstr; t++) status[t] = kCandidate; +#ifdef USE_OLDER_MERGING + Bool_t usedThis = kFALSE; + Bool_t usedPrev = kFALSE; +#endif for(UShort_t t = 0; t < nstr; t++) { output.SetMultiplicity(d,r,s,t,0.); Float_t mult = SignalInStrip(input,d,r,s,t); + // Get the pseudo-rapidity + Double_t eta = input.Eta(d,r,s,t); + Double_t phi = input.Phi(d,r,s,t) * TMath::Pi() / 180.; + if (s == 0) output.SetEta(d,r,s,t,eta); // Keep dead-channel information. if(mult == AliESDFMD::kInvalidMult) output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult); // If no signal or dead strip, go on. - if (mult == AliESDFMD::kInvalidMult || mult == 0) continue; + if (mult == AliESDFMD::kInvalidMult || mult == 0) { + if (mult == 0) histos->fSum->Fill(eta,phi,mult); + status[t] = kNone; + continue; + } - // Get the pseudo-rapidity - Double_t eta = input.Eta(d,r,s,t); - // Fill the diagnostics histogram histos->fBefore->Fill(mult); // Get next and previous signal - if any Double_t prevE = 0; Double_t nextE = 0; + Status prevStatus = (t == 0 ? kNone : status[t-1]); + Status thisStatus = status[t]; + Status nextStatus = (t == nstr-1 ? kNone : status[t+1]); if (t != 0) { prevE = SignalInStrip(input,d,r,s,t-1); if (prevE == AliESDFMD::kInvalidMult) prevE = 0; @@ -225,20 +299,60 @@ AliFMDSharingFilter::Filter(const AliESDFMD& input, nextE = SignalInStrip(input,d,r,s,t+1); if (nextE == AliESDFMD::kInvalidMult) nextE = 0; } + if (t != 0) histos->fNeighborsBefore->Fill(prevE, mult); +#ifdef USE_OLDER_MERGING Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE, lowFlux,d,r,s,t, usedPrev,usedThis); - if (mergedEnergy > lowCut) histos->Incr(); + status[t] = (usedPrev ? kMergedWithOther : kNone); + if (t != nstr - 1) status[t] = (usedThis ? kMergedWithOther : kNone); +#else + Double_t mergedEnergy = MultiplicityOfStrip(mult, prevE, nextE, + eta, lowFlux, + d, r, s, t, + prevStatus, + thisStatus, + nextStatus); + if (t != 0) status[t-1] = prevStatus; + if (t != nstr-1) status[t+1] = nextStatus; + status[t] = thisStatus; + +#endif + // If we're processing on non-angle corrected data, we + // should do the angle correction here + if (!fCorrectAngles) + mergedEnergy = AngleCorrect(mergedEnergy, eta); + if (mergedEnergy > 0) histos->Incr(); + + if (t != 0) + histos->fNeighborsAfter->Fill(output.Multiplicity(d,r,s,t-1), + mergedEnergy); + histos->fBeforeAfter->Fill(mult, mergedEnergy); histos->fAfter->Fill(mergedEnergy); + histos->fSum->Fill(eta,phi,mergedEnergy); output.SetMultiplicity(d,r,s,t,mergedEnergy); - output.SetEta(d,r,s,t,eta); } // for strip + for (UShort_t t = 0; t < nstr; t++) { + if (fOper) fOper->operator()(d, r, s, t) = status[t]; + switch (status[t]) { + case kNone: nNone++; break; + case kCandidate: nCandidate++; break; + case kMergedWithOther: nMerged++; break; + case kMergedInto: nSummed++; break; + } + } } // for sector } // for ring } // for detector + fSummed->Fill(kNone, nNone); + fSummed->Fill(kCandidate, nCandidate); + fSummed->Fill(kMergedWithOther, nMerged); + fSummed->Fill(kMergedInto, nSummed); + DBGL(1, Form("none=%9d, candidate=%9d, merged=%9d, summed=%9d", + nNone, nCandidate, nMerged, nSummed)); next.Reset(); while ((o = static_cast(next()))) o->Finish(); @@ -268,17 +382,26 @@ AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input, // The energy signal // Double_t mult = input.Multiplicity(d,r,s,t); - if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult; - - if (fCorrectAngles && !input.IsAngleCorrected()) - mult = AngleCorrect(mult, input.Eta(d,r,s,t)); - else if (!fCorrectAngles && input.IsAngleCorrected()) - mult = DeAngleCorrect(mult, input.Eta(d,r,s,t)); + // In case of + // - bad value (invalid or 0) + // - we want angle corrected and data is + // - we don't want angle corrected and data isn't + // just return read value + if (mult == AliESDFMD::kInvalidMult || + mult == 0 || + (fCorrectAngles && input.IsAngleCorrected()) || + (!fCorrectAngles && !input.IsAngleCorrected())) + return mult; + + // If we want angle corrected data, correct it, + // otherwise de-correct it + if (fCorrectAngles) mult = AngleCorrect(mult, input.Eta(d,r,s,t)); + else mult = DeAngleCorrect(mult, input.Eta(d,r,s,t)); return mult; } //_____________________________________________________________________ Double_t -AliFMDSharingFilter::GetLowCut() const +AliFMDSharingFilter::GetLowCut(UShort_t, Char_t, Double_t) const { // // Get the low cut. Normally, the low cut is taken to be the lower @@ -294,7 +417,8 @@ AliFMDSharingFilter::GetLowCut() const //_____________________________________________________________________ Double_t -AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const +AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, + Double_t eta, Bool_t errors) const { // // Get the high cut. The high cut is defined as the @@ -307,22 +431,130 @@ AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const // Get the high cut. The high cut is defined as the // most-probably-value peak found from the energy distributions, minus // 2 times the width of the corresponding Landau. - AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta); - Double_t delta = 1024; - Double_t xi = 1024; - if (!fit) { - AliError(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta)); + AliFMDCorrELossFit* fits = fcm.GetELossFit(); + + return fits->GetLowerBound(d, r, eta, fNXi, errors, fIncludeSigma); +} + +//_____________________________________________________________________ +namespace { + const char* status2String(AliFMDSharingFilter::Status s) + { + switch(s) { + case AliFMDSharingFilter::kNone: return "none "; + case AliFMDSharingFilter::kCandidate: return "candidate"; + case AliFMDSharingFilter::kMergedWithOther: return "merged "; + case AliFMDSharingFilter::kMergedInto: return "summed "; + } + return "unknown "; } - else { - delta = fit->fDelta; - xi = fit->fXi; +} + +//_____________________________________________________________________ +Double_t +AliFMDSharingFilter::MultiplicityOfStrip(Double_t thisE, + Double_t prevE, + Double_t nextE, + Double_t eta, + Bool_t lowFlux, + UShort_t d, + Char_t r, + UShort_t s, + UShort_t t, + Status& prevStatus, + Status& thisStatus, + Status& nextStatus) const +{ + Double_t lowCut = GetLowCut(d, r, eta); + + DBG(3,Form("FMD%d%c[%2d,%3d]: this=%9f(%s), prev=%9f(%s), next=%9f(%s)", + d, r, s, t, + thisE, status2String(thisStatus), + prevE, status2String(prevStatus), + nextE, status2String(nextStatus))); + + // If below cut, then modify zero signal and make sure the next + // strip is considered a candidate. + if (thisE < lowCut || thisE > 20) { + thisStatus = kNone; + DBGL(3,Form(" %9f<%9f || %9f>20, 0'ed", thisE, lowCut, thisE)); + if (prevStatus == kCandidate) prevStatus = kNone; + return 0; + } + // It this strip was merged with the previous strip, then + // make the next strip a candidate and zero the value in this strip. + if (thisStatus == kMergedWithOther) { + DBGL(3,Form(" Merged with other, 0'ed")); + return 0; } - if (delta > 100) { - AliWarning(Form("FMD%d%c, eta=%f, Delta=%f, xi=%f", d, r, eta, delta, xi)); + // Get the upper cut on the sharing + Double_t highCut = GetHighCut(d, r, eta ,false); + if (highCut <= 0) { + thisStatus = kNone; + return 0; + } + + // Only for low-flux events + if (lowFlux) { + // If this is smaller than the next, and the next is larger + // then the cut, mark both as candidates for sharing. + // They should be be merged on the next iteration + if (thisE < nextE && nextE > highCut) { + thisStatus = kCandidate; + nextStatus = kCandidate; + return 0; + } + } + + // Variable to sum signal in + Double_t totalE = thisE; + + // If the previous signal was between the two cuts, and is still + // marked as a candidate , then merge it in here. + if (prevStatus == kCandidate && prevE > lowCut && prevE < highCut) { + totalE += prevE; + prevStatus = kMergedWithOther; + DBG(3, Form(" Prev candidate %9f<%9f<%9f->%9f", + lowCut, prevE, highCut, totalE)); + } + + // If the next signal is between the two cuts, then merge it here + if (nextE > lowCut && nextE < highCut) { + totalE += nextE; + nextStatus = kMergedWithOther; + DBG(3, Form(" Next %9f<%9f<%9f->%9f", lowCut, nextE, highCut,totalE)); + } + + // If the signal is bigger than the high-cut, then return the value + if (totalE > highCut) { + if (totalE > thisE) + thisStatus = kMergedInto; + else + thisStatus = kNone; + DBGL(3, Form(" %9f>%f9 (was %9f) -> %9f %s", + totalE, highCut, thisE, totalE,status2String(thisStatus))); + return totalE; + } + else { + // This is below cut, so flag next as a candidate + DBG(3, Form(" %9f<=%9f, next candidate", totalE, highCut)); + nextStatus = kCandidate; + } + + // If the total signal is smaller than low cut then zero this and kill this + if (totalE < lowCut) { + DBGL(3, Form(" %9f<%9f (was %9f), zeroed", totalE, lowCut, thisE)); + thisStatus = kNone; + return 0; } - Double_t highCut = (delta - fNXi * xi); - return highCut; + + // If total signal not above high cut or lower than low cut, + // mark this as a candidate for merging into the next, and zero signal + DBGL(3, Form(" %9f<%9f<%9f (was %9f), zeroed, candidate", + lowCut, totalE, highCut, thisE)); + thisStatus = kCandidate; + return 0; } //_____________________________________________________________________ @@ -361,7 +593,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult, // IF the multiplicity is very large, or below the cut, reject it, and // flag it as candidate for sharing - Double_t lowCut = GetLowCut(); + Double_t lowCut = GetLowCut(d,r,eta); if(mult > 12 || mult < lowCut) { usedThis = kFALSE; usedPrev = kFALSE; @@ -379,7 +611,12 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult, //analyse and perform sharing on one strip // Get the high cut - Double_t highCut = GetHighCut(d, r, eta); + Double_t highCut = GetHighCut(d, r, eta, false); + if (highCut <= 0) { + usedThis = kFALSE; + usedPrev = kTRUE; + return 0; + } // If this signal is smaller than the next, and the next signal is smaller // than then the high cut, and it's a low-flux event, then mark this and @@ -406,7 +643,7 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult, totalE += prevE; // If the next signal is larger than the low cut, and - // the next signal is smaller than the low cut, then add the next signal + // the next signal is smaller than the high cut, then add the next signal // to this, and marked the next signal as used if(nextE > lowCut && nextE < highCut ) { totalE += nextE; @@ -415,8 +652,8 @@ AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult, // If we're processing on non-angle corrected data, we should do the // angle correction here - if (!fCorrectAngles) - totalE = AngleCorrect(totalE, eta); + // if (!fCorrectAngles) + // totalE = AngleCorrect(totalE, eta); // Fill post processing histogram // if(totalE > fLowCut) @@ -494,8 +731,17 @@ AliFMDSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents) TIter next(&fRingHistos); RingHistos* o = 0; - while ((o = static_cast(next()))) + THStack* sums = new THStack("sums", "Sum of ring signals"); + while ((o = static_cast(next()))) { o->ScaleHistograms(d, nEvents); + TH1D* sum = o->fSum->ProjectionX(o->GetName(), 1, o->fSum->GetNbinsY(),"e"); + sum->Scale(1., "width"); + sum->SetTitle(o->GetName()); + sum->SetDirectory(0); + sum->SetYTitle("#sum #Delta/#Delta_{mip}"); + sums->Add(sum); + } + d->Add(sums); } //____________________________________________________________________ @@ -513,6 +759,26 @@ AliFMDSharingFilter::DefineOutput(TList* dir) TList* d = new TList; d->SetName(GetName()); dir->Add(d); + + fSummed = new TH2I("operations", "Strip operations", + kMergedInto, kNone-.5, kMergedInto+.5, + 51201, -.5, 51200.5); + fSummed->SetXTitle("Operation"); + fSummed->SetYTitle("# of strips"); + fSummed->SetZTitle("Events"); + fSummed->GetXaxis()->SetBinLabel(kNone, "None"); + fSummed->GetXaxis()->SetBinLabel(kCandidate, "Candidate"); + fSummed->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other"); + fSummed->GetXaxis()->SetBinLabel(kMergedInto, "Merged into"); + fSummed->SetDirectory(0); + d->Add(fSummed); + + fHighCuts = new TH2D("highCuts", "High cuts used", 1,0,1, 1,0,1); + fHighCuts->SetXTitle("#eta"); + fHighCuts->SetDirectory(0); + d->Add(fHighCuts); + + TIter next(&fRingHistos); RingHistos* o = 0; while ((o = static_cast(next()))) { @@ -533,10 +799,13 @@ AliFMDSharingFilter::Print(Option_t* /*option*/) const for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; ind[gROOT->GetDirLevel()] = '\0'; std::cout << ind << ClassName() << ": " << GetName() << '\n' + << std::boolalpha + << ind << " Debug: " << fDebug << "\n" << ind << " Low cut: " << fLowCut << '\n' << ind << " N xi factor: " << fNXi << '\n' - << ind << " Use corrected angles: " - << (fCorrectAngles ? "yes" : "no") << std::endl; + << ind << " Include sigma in cut: " << fIncludeSigma << '\n' + << ind << " Use corrected angles: " << fCorrectAngles + << std::noboolalpha << std::endl; } //==================================================================== @@ -544,6 +813,10 @@ AliFMDSharingFilter::RingHistos::RingHistos() : AliForwardUtil::RingHistos(), fBefore(0), fAfter(0), + fBeforeAfter(0), + fNeighborsBefore(0), + fNeighborsAfter(0), + fSum(0), fHits(0), fNHits(0) { @@ -558,6 +831,10 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r) : AliForwardUtil::RingHistos(d,r), fBefore(0), fAfter(0), + fBeforeAfter(0), + fNeighborsBefore(0), + fNeighborsAfter(0), + fSum(0), fHits(0), fNHits(0) { @@ -568,28 +845,51 @@ AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r) // d detector // r ring // - fBefore = new TH1D(Form("%s_ESD_Eloss", fName.Data()), - Form("Energy loss in %s (reconstruction)", fName.Data()), + fBefore = new TH1D("esdEloss", "Energy loss (reconstruction)", 1000, 0, 25); - fAfter = new TH1D(Form("%s_Ana_Eloss", fName.Data()), - Form("Energy loss in %s (sharing corrected)", - fName.Data()), 1000, 0, 25); fBefore->SetXTitle("#Delta E/#Delta E_{mip}"); fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})"); - fBefore->SetFillColor(kRed+1); + fBefore->SetFillColor(Color()); fBefore->SetFillStyle(3001); - // fBefore->Sumw2(); fBefore->SetDirectory(0); - fAfter->SetXTitle("#Delta E/#Delta E_{mip}"); - fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})"); - fAfter->SetFillColor(kBlue+1); - fAfter->SetFillStyle(3001); - // fAfter->Sumw2(); + + fAfter = static_cast(fBefore->Clone("anaEloss")); + fAfter->SetTitle("Energy loss in %s (sharing corrected)"); + fAfter->SetFillColor(Color()+2); fAfter->SetDirectory(0); - fHits = new TH1D(Form("%s_Hits", fName.Data()), - Form("Number of hits in %s", fName.Data()), - 200, 0, 200000); + Double_t max = 15; + Double_t min = -1; + Int_t n = int((max-min) / (max / 300)); + fBeforeAfter = new TH2D("beforeAfter", "Before and after correlation", + n, min, max, n, min, max); + fBeforeAfter->SetXTitle("#Delta E/#Delta E_{mip} before"); + fBeforeAfter->SetYTitle("#Delta E/#Delta E_{mip} after"); + fBeforeAfter->SetZTitle("Correlation"); + fBeforeAfter->SetDirectory(0); + + fNeighborsBefore = static_cast(fBeforeAfter->Clone("neighborsBefore")); + fNeighborsBefore->SetTitle("Correlation of neighbors befire"); + fNeighborsBefore->SetXTitle("#Delta E_{i}/#Delta E_{mip}"); + fNeighborsBefore->SetYTitle("#Delta E_{i+1}/#Delta E_{mip}"); + fNeighborsBefore->SetDirectory(0); + + fNeighborsAfter = + static_cast(fNeighborsBefore->Clone("neighborsAfter")); + fNeighborsAfter->SetTitle("Correlation of neighbors after"); + fNeighborsAfter->SetDirectory(0); + + fSum = new TH2D("summed", "Summed signal", 200, -4, 6, + (fRing == 'I' || fRing == 'i' ? 20 : 40), 0, 2*TMath::Pi()); + fSum->SetDirectory(0); + fSum->Sumw2(); + fSum->SetMarkerColor(Color()); + // fSum->SetFillColor(Color()); + fSum->SetXTitle("#eta"); + fSum->SetYTitle("#varphi [radians]"); + fSum->SetZTitle("#sum #Delta/#Delta_{mip}(#eta,#varphi) "); + + fHits = new TH1D("hits", "Number of hits", 200, 0, 200000); fHits->SetDirectory(0); fHits->SetXTitle("# of hits"); fHits->SetFillColor(kGreen+1); @@ -599,6 +899,10 @@ AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o) : AliForwardUtil::RingHistos(o), fBefore(o.fBefore), fAfter(o.fAfter), + fBeforeAfter(o.fBeforeAfter), + fNeighborsBefore(o.fNeighborsBefore), + fNeighborsAfter(o.fNeighborsAfter), + fSum(o.fSum), fHits(o.fHits), fNHits(o.fNHits) { @@ -631,10 +935,14 @@ AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o) if (fAfter) delete fAfter; if (fHits) delete fHits; - fBefore = static_cast(o.fBefore->Clone()); - fAfter = static_cast(o.fAfter->Clone()); - fHits = static_cast(o.fHits->Clone()); - + fBefore = static_cast(o.fBefore->Clone()); + fAfter = static_cast(o.fAfter->Clone()); + fBeforeAfter = static_cast(o.fBeforeAfter->Clone()); + fNeighborsBefore = static_cast(o.fNeighborsBefore->Clone()); + fNeighborsAfter = static_cast(o.fNeighborsAfter->Clone()); + fHits = static_cast(o.fHits->Clone()); + fSum = static_cast(o.fSum->Clone()); + return *this; } //____________________________________________________________________ @@ -643,9 +951,6 @@ AliFMDSharingFilter::RingHistos::~RingHistos() // // Destructor // - if (fBefore) delete fBefore; - if (fAfter) delete fAfter; - if (fHits) delete fHits; } //____________________________________________________________________ void @@ -673,12 +978,13 @@ AliFMDSharingFilter::RingHistos::ScaleHistograms(const TList* dir, TList* l = GetOutputList(dir); if (!l) return; - TH1D* before = static_cast(l->FindObject(Form("%s_ESD_ELoss", - fName.Data()))); - TH1D* after = static_cast(l->FindObject(Form("%s_Ana_ELoss", - fName.Data()))); + TH1D* before = static_cast(l->FindObject("esdELoss")); + TH1D* after = static_cast(l->FindObject("anaELoss")); + TH2D* summed = static_cast(l->FindObject("summed")); if (before) before->Scale(1./nEvents); if (after) after->Scale(1./nEvents); + if (summed) summed->Scale(1./nEvents); + } //____________________________________________________________________ @@ -692,9 +998,15 @@ AliFMDSharingFilter::RingHistos::Output(TList* dir) // dir where to store // TList* d = DefineOutputList(dir); + d->Add(fBefore); d->Add(fAfter); + d->Add(fBeforeAfter); + d->Add(fNeighborsBefore); + d->Add(fNeighborsAfter); d->Add(fHits); + d->Add(fSum); + dir->Add(d); } diff --git a/PWG2/FORWARD/analysis2/AliFMDSharingFilter.h b/PWG2/FORWARD/analysis2/AliFMDSharingFilter.h index f103bc81035..14bfe90ad91 100644 --- a/PWG2/FORWARD/analysis2/AliFMDSharingFilter.h +++ b/PWG2/FORWARD/analysis2/AliFMDSharingFilter.h @@ -24,6 +24,7 @@ class AliESDFMD; class TAxis; class TList; class TH2; +class AliFMDFloatMap; /** * Class to do the sharing correction. That is, a filter that merges @@ -53,6 +54,12 @@ class TH2; class AliFMDSharingFilter : public TNamed { public: + enum Status { + kNone = 1, + kCandidate = 2, + kMergedWithOther = 3, + kMergedInto = 4 + }; /** * Destructor */ @@ -82,6 +89,11 @@ public: */ AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o); + /** + * Initialize + * + */ + void Init(); /** * Set the low cut used for sharing * @@ -107,14 +119,21 @@ public: * otherwise use de-corrected signals. In the final output, the * signals are always angle corrected. */ - void UseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; } + void SetUseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; } /** * Set the number of landau width to subtract from the most probably * value to get the high cut for the merging algorithm. * * @param n Number of @f$ \xi@f$ */ - void SetNXi(Short_t n) { fNXi = n; } + void SetNXi(Double_t n) { fNXi = n; } + /** + * Whether to include sigma in the number subtracted from the MPV to + * get the high cut + * + * @param u If true, then high cut is @f$ \Delta_{mp} - n(\xi+\sigma)@f$ + */ + void SetIncludeSigma(Bool_t u) { fIncludeSigma = u; } /** * Filter the input AliESDFMD object * @@ -213,6 +232,10 @@ protected: void ScaleHistograms(const TList* dir, Int_t nEvents); TH1D* fBefore; // Distribution of signals before filter TH1D* fAfter; // Distribution of signals after filter + TH2D* fBeforeAfter; // Correlation of before and after + TH2D* fNeighborsBefore; // Correlation of neighbors + TH2D* fNeighborsAfter; // Correlation of neighbors + TH2D* fSum; // Summed signal TH1D* fHits; // Distribution of hit strips. Int_t fNHits; // Number of hit strips per event ClassDef(RingHistos,1); @@ -270,6 +293,18 @@ protected: UShort_t t, Bool_t& usedPrev, Bool_t& usedThis) const; + Double_t MultiplicityOfStrip(Double_t thisE, + Double_t prevE, + Double_t nextE, + Double_t eta, + Bool_t lowFlux, + UShort_t d, + Char_t r, + UShort_t s, + UShort_t t, + Status& prevStatus, + Status& thisStatus, + Status& nextStatus) const; /** * Angle correct the signal * @@ -288,27 +323,45 @@ protected: * @return Angle un-corrected signal */ Double_t DeAngleCorrect(Double_t mult, Double_t eta) const; - /** + /** * Get the high cut. The high cut is defined as the * most-probably-value peak found from the energy distributions, minus * 2 times the width of the corresponding Landau. + * + * @param d Detector + * @param r Ring + * @param eta Eta value + * @param errors If false, do not show errors + * + * @return 0 or less on failure, otherwise @f$\Delta_{mp}-n\xi@f$ */ - virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta) const; + virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta, + Bool_t errors=true) const; /** * Get the low cut. Normally, the low cut is taken to be the lower * value of the fit range used when generating the energy loss fits. * However, if fLowCut is set (using SetLowCit) to a value greater * than 0, then that value is used. + * + * @param d Detector + * @param r Ring + * @param eta Eta value + * + * @return */ - virtual Double_t GetLowCut() const; + virtual Double_t GetLowCut(UShort_t d, Char_t r, Double_t eta) const; TList fRingHistos; // List of histogram containers Double_t fLowCut; // Low cut on sharing Bool_t fCorrectAngles; // Whether to work on angle corrected signals - Short_t fNXi; // Number of xi's from Delta to stop merging + Double_t fNXi; // Number of xi's from Delta to stop merging + Bool_t fIncludeSigma; // Whether to include sigma in cut + TH2* fSummed; // Operations histogram + TH2* fHighCuts; // High cuts used + AliFMDFloatMap* fOper; // Operation done per strip Int_t fDebug; // Debug level - ClassDef(AliFMDSharingFilter,1); // + ClassDef(AliFMDSharingFilter,2); // }; #endif diff --git a/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx b/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx index 1fda2c5c4ba..5376027afc0 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx @@ -23,7 +23,8 @@ #include "AliInputEventHandler.h" #include "AliStack.h" #include "AliMCEvent.h" -//#include "AliFMDStripIndex.h" +#include "AliAODForwardMult.h" +#include "AliFMDStripIndex.h" #include #include #include @@ -56,10 +57,11 @@ namespace { //==================================================================== AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask() : AliAnalysisTaskSE(), + fInspector(), + fFirstEvent(true), fHEvents(0), fHEventsTr(0), fHEventsTrVtx(0), - fHEventsVtx(0), fHTriggers(0), fPrimaryInnerAll(0), fPrimaryOuterAll(0), @@ -89,11 +91,12 @@ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask() //____________________________________________________________________ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name) - : AliAnalysisTaskSE(name), + : AliAnalysisTaskSE(name), + fInspector("eventInspector"), + fFirstEvent(true), fHEvents(0), fHEventsTr(0), fHEventsTrVtx(0), - fHEventsVtx(0), fHTriggers(0), fPrimaryInnerAll(0), fPrimaryOuterAll(0), @@ -120,15 +123,17 @@ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name) // name Name of task // DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); } //____________________________________________________________________ AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o) : AliAnalysisTaskSE(o), + fInspector(o.fInspector), + fFirstEvent(o.fFirstEvent), fHEvents(o.fHEvents), fHEventsTr(o.fHEventsTr), fHEventsTrVtx(o.fHEventsTrVtx), - fHEventsVtx(o.fHEventsVtx), fHTriggers(o.fHTriggers), fPrimaryInnerAll(o.fPrimaryInnerAll), fPrimaryOuterAll(o.fPrimaryOuterAll), @@ -169,6 +174,8 @@ AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o) // Return: // Reference to this object // + fInspector = o.fInspector; + fFirstEvent = o.fFirstEvent; fHEventsTr = o.fHEventsTr; fHEventsTrVtx = o.fHEventsTrVtx; fHTriggers = o.fHTriggers; @@ -333,7 +340,8 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects() // // fList = new TList; - fList->SetName(GetName()); + fList->SetOwner(); + fList->SetName(Form("%sSums", GetName())); fHEvents = new TH1I(GetEventName(false,false), "Number of all events", @@ -371,19 +379,6 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects() fHEventsTrVtx->SetDirectory(0); fList->Add(fHEventsTrVtx); - fHEventsVtx = new TH1I(GetEventName(false,true), - "Number of events w/vertex", - fVtxAxis.GetNbins(), - fVtxAxis.GetXmin(), - fVtxAxis.GetXmax()); - fHEventsVtx->SetXTitle("v_{z} [cm]"); - fHEventsVtx->SetYTitle("# of events"); - fHEventsVtx->SetFillColor(kBlue+1); - fHEventsVtx->SetFillStyle(3001); - fHEventsVtx->SetDirectory(0); - fList->Add(fHEventsVtx); - - fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10); fHTriggers->SetFillColor(kRed+1); fHTriggers->SetFillStyle(3001); @@ -444,6 +439,9 @@ AliForwardMCCorrectionsTask::UserCreateOutputObjects() strips->Add(fStripsFMD3o); fList->Add(strips); + fInspector.Init(fVtxAxis); + fInspector.DefineOutput(fList); + PostData(1, fList); } //____________________________________________________________________ @@ -470,52 +468,57 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) AliWarning("No ESD event found for input event"); return; } - - // Get the particle stack - AliStack* stack = mcEvent->Stack(); - // Get the event generate header - AliHeader* mcHeader = mcEvent->Header(); - AliGenEventHeader* genHeader = mcHeader->GenEventHeader(); - - // Get the generator vertex - TArrayF mcVertex; - genHeader->PrimaryVertex(mcVertex); - Double_t mcVtxZ = mcVertex.At(2); - - // Check the MC vertex is in range - Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ); - if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) { -#ifdef VERBOSE - AliWarning(Form("MC v_z=%f is out of rante [%,%f]", - mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); -#endif - return; + //--- Read run information ----------------------------------------- + if (fFirstEvent && esd->GetESDRun()) { + fInspector.ReadRunDetails(esd); + + AliInfo(Form("Initializing with parameters from the ESD:\n" + " AliESDEvent::GetBeamEnergy() ->%f\n" + " AliESDEvent::GetBeamType() ->%s\n" + " AliESDEvent::GetCurrentL3() ->%f\n" + " AliESDEvent::GetMagneticField()->%f\n" + " AliESDEvent::GetRunNumber() ->%d\n", + esd->GetBeamEnergy(), + esd->GetBeamType(), + esd->GetCurrentL3(), + esd->GetMagneticField(), + esd->GetRunNumber())); + + fFirstEvent = false; } - // UInt_t triggers; - // Bool_t gotTrigggers = false; - Bool_t gotInel = false; - // Double_t vZ; - Bool_t gotVertex = false; -#if 0 - // Use event inspector instead - // Get the triggers - UInt_t triggers = 0; - Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers); - Bool_t gotInel = triggers & AliAODForwardMult::kInel; - - // Get the ESD vertex - Double_t vZ = -1000000; - Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ); -#endif - + // Get the particle stack + AliStack* stack = mcEvent->Stack(); + // Some variables + UInt_t triggers; // Trigger bits + Bool_t lowFlux; // Low flux flag + UShort_t iVz; // Vertex bin from ESD + Double_t vZ; // Z coordinate from ESD + Double_t cent; // Centrality + UShort_t iVzMc; // Vertex bin from MC + Double_t vZMc; // Z coordinate of IP vertex from MC + Double_t b; // Impact parameter + Int_t nPart; // Number of participants + Int_t nBin; // Number of binary collisions + Double_t phiR; // Reaction plane from MC + UShort_t nClusters;// Number of SPD clusters + // Process the data + UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, + cent, nClusters); + /*UInt_t retMC =*/ fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, + b, nPart, nBin, phiR); + + // Bool_t isInelMC = true; // (triggers & AliAODForwardMult::kB); + // Bool_t isNSDMC = (triggers & AliAODForwardMult::kMCNSD); + Bool_t isInel = triggers & AliAODForwardMult::kInel; + // Bool_t isNSD = triggers & AliAODForwardMult::kNSD; + Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk; // Fill the event count histograms - if (gotInel) fHEventsTr->Fill(mcVtxZ); - if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ); - if (gotVertex) fHEventsVtx->Fill(mcVtxZ); - fHEvents->Fill(mcVtxZ); + if (isInel) fHEventsTr->Fill(vZMc); + if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc); + fHEvents->Fill(vZMc); // Cache of the hits AliFMDFloatMap hitsByStrip; @@ -539,10 +542,10 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) Double_t phi = particle->Phi(); if (isCharged && isPrimary) - FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi); + FillPrimary(isInel, hasVtx, vZMc, eta, phi); // For the rest, ignore non-collisions, and events out of vtx range - if (!gotInel || gotVertex) continue; + if (!isInel || !hasVtx) continue; Int_t nTrRef = particle->GetNumberOfTrackReferences(); for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { @@ -558,7 +561,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) // Get the detector coordinates UShort_t d = 0, s = 0, t = 0; Char_t r = '\0'; - // AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t); + AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t); // Check if mother (?) is charged and that we haven't dealt with it // already @@ -572,7 +575,7 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t); // Fill strip information into histograms - FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1); + FillStrip(d, r, vZMc, eta, phi, hitsByStrip(d,r,s,t) == 1); // Set the last processed track number - marking it as done for // this strip @@ -588,20 +591,21 @@ AliForwardMCCorrectionsTask::UserExec(Option_t*) //____________________________________________________________________ void -AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx, - Double_t vz, Double_t eta, Double_t phi) +AliForwardMCCorrectionsTask::FillPrimary(Bool_t isInel, Bool_t hasVtx, + Double_t vz, + Double_t eta, Double_t phi) { // // Fill in primary information // // Parameters: - // gotInel Got INEL trigger from ESD + // isInel Got INEL trigger from ESD // gotVtx Got vertex Z from ESD // vz @f$z@f$ coordinate of interation point // eta Pseudo rapidity // phi Azimuthal angle // - if (gotInel && gotVtx) { + if (isInel && hasVtx) { // This takes the place of hPrimary_FMD__vtx and // Analysed_FMD_vtx fPrimaryInnerTrVtx->Fill(vz,eta,phi); @@ -698,32 +702,32 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*) // option Not used // - TList* list = dynamic_cast(GetOutputData(1)); - if (!list) { + fList = dynamic_cast(GetOutputData(1)); + if (!fList) { AliError("No output list defined"); return; } - TList* primaries = static_cast(list->FindObject("primaries")); - TList* hits = static_cast(list->FindObject("hits")); - TList* strips = static_cast(list->FindObject("strips")); + TList* primaries = static_cast(fList->FindObject("primaries")); + TList* hits = static_cast(fList->FindObject("hits")); + TList* strips = static_cast(fList->FindObject("strips")); if (!primaries || !hits || !strips) { AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)", - list->GetName(), primaries, hits, strips)); + fList->GetName(), primaries, hits, strips)); return; } TH1I* eventsAll = - static_cast(list->FindObject(GetEventName(false,false))); + static_cast(fList->FindObject(GetEventName(false,false))); TH1I* eventsTr = - static_cast(list->FindObject(GetEventName(true,false))); + static_cast(fList->FindObject(GetEventName(true,false))); TH1I* eventsVtx = - static_cast(list->FindObject(GetEventName(false,true))); + static_cast(fList->FindObject(GetEventName(false,true))); TH1I* eventsTrVtx = - static_cast(list->FindObject(GetEventName(true,true))); + static_cast(fList->FindObject(GetEventName(true,true))); if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) { AliError(Form("Could not find all event histograms in %s " - "(a=%p,t=%p,v=%p,tv=%p)", list->GetName(), + "(a=%p,t=%p,v=%p,tv=%p)", fList->GetName(), eventsAll, eventsTr, eventsVtx, eventsTrVtx)); return; } @@ -738,7 +742,7 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*) static_cast(primaries->FindObject(GetPrimaryName('O',true))); if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) { AliError(Form("Could not find all primary particle histograms in %s " - "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(), + "(ai=%p,ao=%p,tvi=%p,tvo=%p)", fList->GetName(), primIAll, primOAll, primITrVtx, primOTrVtx)); return; } @@ -767,6 +771,11 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*) return; } + // Output list + TList* output = new TList; + output->SetOwner(); + output->SetName(Form("%sResults", GetName())); + // Calculate the over-all event selection efficiency TH1D* selEff = new TH1D("selEff", "Event selection efficiency", fVtxAxis.GetNbins(), @@ -778,14 +787,15 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*) selEff->SetFillStyle(3001); selEff->Add(eventsAll); selEff->Divide(eventsTrVtx); - list->Add(selEff); + output->Add(selEff); // Loop over vertex bins and do vertex dependendt stuff for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) { - // Make a sub-list in the output + // Make a sub-fList in the output TList* vl = new TList; + vl->SetOwner(); vl->SetName(Form("vtx%02d", v)); - list->Add(vl); + output->Add(vl); // Get event counts Int_t nEventsAll = Int_t(eventsAll->GetBinContent(v)); @@ -858,10 +868,11 @@ AliForwardMCCorrectionsTask::Terminate(Option_t*) doubleHit->SetFillStyle(3001); doubleHit->Sumw2(); doubleHit->Divide(hStrips); - list->Add(doubleHit); - } - } - } + output->Add(doubleHit); + } // for q + } // for d + } // for v + PostData(2, output); } //____________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.h b/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.h index 80b7f09b4d5..7701145b778 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.h +++ b/PWG2/FORWARD/analysis2/AliForwardMCCorrectionsTask.h @@ -15,8 +15,8 @@ */ #include #include +#include "AliFMDMCEventInspector.h" #include -class AliFMDAnaParameters; class AliESDEvent; class TH2D; class TH3D; @@ -186,10 +186,11 @@ protected: void FillStrip(UShort_t d, Char_t r, Double_t vz, Double_t eta, Double_t phi, Bool_t first); + AliFMDMCEventInspector fInspector; // Event inspector + Bool_t fFirstEvent; // First event flag TH1I* fHEvents; // All Events TH1I* fHEventsTr; // Histogram of events w/trigger TH1I* fHEventsTrVtx; // Events w/trigger and vertex - TH1I* fHEventsVtx; // Events w/vertex TH1I* fHTriggers; // Triggers TH3D* fPrimaryInnerAll; // Distribution of primaries - all events TH3D* fPrimaryOuterAll; // Distribution of primaries - all events diff --git a/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx b/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx index 4e29b526770..5b5d77506b3 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx @@ -38,6 +38,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask() fMCESDFMD(), fMCHistos(), fMCAODFMD(), + fRingSums(), + fMCRingSums(), fPrimary(0), fEventInspector(), fSharingFilter(), @@ -61,6 +63,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name) fMCESDFMD(), fMCHistos(), fMCAODFMD(kTRUE), + fRingSums(), + fMCRingSums(), fPrimary(0), fEventInspector("event"), fSharingFilter("sharing"), @@ -88,6 +92,8 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul fMCESDFMD(o.fMCESDFMD), fMCHistos(o.fMCHistos), fMCAODFMD(o.fMCAODFMD), + fRingSums(o.fRingSums), + fMCRingSums(o.fMCRingSums), fPrimary(o.fPrimary), fEventInspector(o.fEventInspector), fSharingFilter(o.fSharingFilter), @@ -130,6 +136,8 @@ AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o) fAODFMD = o.fAODFMD; fMCHistos = o.fMCHistos; fMCAODFMD = o.fMCAODFMD; + fRingSums = o.fRingSums; + fMCRingSums = o.fMCRingSums; fPrimary = o.fPrimary; fList = o.fList; @@ -152,6 +160,14 @@ AliForwardMCMultiplicityTask::SetDebug(Int_t dbg) fCorrections.SetDebug(dbg); fHistCollector.SetDebug(dbg); } +//____________________________________________________________________ +void +AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use) +{ + fSharingFilter.SetOnlyPrimary(use); + fCorrections.SetSecondaryForMC(!use); +} + //____________________________________________________________________ void AliForwardMCMultiplicityTask::InitializeSubs() @@ -169,6 +185,8 @@ AliForwardMCMultiplicityTask::InitializeSubs() fAODFMD.Init(*pe); fMCHistos.Init(*pe); fMCAODFMD.Init(*pe); + fRingSums.Init(*pe); + fMCRingSums.Init(*pe); fHData = static_cast(fAODFMD.GetHistogram().Clone("d2Ndetadphi")); fHData->SetStats(0); @@ -176,7 +194,42 @@ AliForwardMCMultiplicityTask::InitializeSubs() fList->Add(fHData); + TList* rings = new TList; + rings->SetName("ringSums"); + rings->SetOwner(); + fList->Add(rings); + + rings->Add(fRingSums.Get(1, 'I')); + rings->Add(fRingSums.Get(2, 'I')); + rings->Add(fRingSums.Get(2, 'O')); + rings->Add(fRingSums.Get(3, 'I')); + rings->Add(fRingSums.Get(3, 'O')); + fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I')); + fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I')); + fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O')); + fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I')); + fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O')); + + TList* mcRings = new TList; + mcRings->SetName("mcRingSums"); + mcRings->SetOwner(); + fList->Add(mcRings); + + mcRings->Add(fMCRingSums.Get(1, 'I')); + mcRings->Add(fMCRingSums.Get(2, 'I')); + mcRings->Add(fMCRingSums.Get(2, 'O')); + mcRings->Add(fMCRingSums.Get(3, 'I')); + mcRings->Add(fMCRingSums.Get(3, 'O')); + fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I')); + fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I')); + fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O')); + fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I')); + fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O')); + + + fEventInspector.Init(*pv); + fSharingFilter.Init(); fDensityCalculator.Init(*pe); fCorrections.Init(*pe); fHistCollector.Init(*pv,*pe); @@ -249,13 +302,14 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) fMCAODFMD.Clear(); fPrimary->Reset(); - Bool_t lowFlux = kFALSE; - UInt_t triggers = 0; - UShort_t ivz = 0; - Double_t vz = 0; - Double_t cent = 0; - UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent); + Bool_t lowFlux = kFALSE; + UInt_t triggers = 0; + UShort_t ivz = 0; + Double_t vz = 0; + Double_t cent = 0; + UShort_t nClusters = 0; + UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, + ivz, vz, cent, nClusters); UShort_t ivzMC = 0; Double_t vzMC = 0; Double_t phiR = 0; @@ -281,11 +335,13 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) fAODFMD.SetSNN(fEventInspector.GetEnergy()); fAODFMD.SetSystem(fEventInspector.GetCollisionSystem()); fAODFMD.SetCentrality(cent); + fAODFMD.SetNClusters(nClusters); fMCAODFMD.SetTriggerBits(triggers); fMCAODFMD.SetSNN(fEventInspector.GetEnergy()); fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem()); fMCAODFMD.SetCentrality(cent); + fMCAODFMD.SetNClusters(nClusters); //All events should be stored - HHD //if (isAccepted) MarkEventForStore(); @@ -346,11 +402,13 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) } fCorrections.CompareResults(fHistos, fMCHistos); - if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) { + if (!fHistCollector.Collect(fHistos, fRingSums, + ivz, fAODFMD.GetHistogram())) { AliWarning("Histogram collector failed"); return; } - if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) { + if (!fHistCollector.Collect(fMCHistos, fMCRingSums, + ivz, fMCAODFMD.GetHistogram())) { AliWarning("MC Histogram collector failed"); return; } @@ -411,6 +469,9 @@ AliForwardMCMultiplicityTask::Terminate(Option_t*) list->Add(dNdeta); list->Add(norm); + MakeRingdNdeta(list, "ringSums", list, "ringResults"); + MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24); + fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral())); fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); diff --git a/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h b/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h index 96c9cd910ca..311f2d85394 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h +++ b/PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.h @@ -99,6 +99,12 @@ public: /** * @} */ + /** + * Process only primary MC tracks + * + * @param use + */ + void SetOnlyPrimary(Bool_t use); /** * @{ * @name Access to sub-algorithms @@ -186,6 +192,8 @@ protected: AliESDFMD fMCESDFMD; // MC 'Sharing corrected' ESD object AliForwardUtil::Histos fMCHistos; // MC Cache histograms AliAODForwardMult fMCAODFMD; // MC Output object + AliForwardUtil::Histos fRingSums; // Cache histograms + AliForwardUtil::Histos fMCRingSums; // Cache histograms TH2D* fPrimary; // Per event primary particles AliFMDMCEventInspector fEventInspector; // Algorithm diff --git a/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx b/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx index e238ac83242..7230dc5d588 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.cxx @@ -27,6 +27,7 @@ #include "AliESDEvent.h" #include #include +#include #include #include @@ -207,6 +208,85 @@ AliForwardMultiplicityBase::MarkEventForStore() const ah->SetFillAOD(kTRUE); } +//____________________________________________________________________ +void +AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, + const char* inName, + TList* output, + const char* outName, + Int_t style) const +{ + // Make dN/deta for each ring found in the input list. + // + // A stack of all the dN/deta is also made for easy drawing. + // + // Note, that the distributions are normalised to the number of + // observed events only - they should be corrected for + if (!input) return; + TList* list = static_cast(input->FindObject(inName)); + if (!list) { + AliWarning(Form("No list %s found in %s", inName, input->GetName())); + return; + } + + TList* out = new TList; + out->SetName(outName); + out->SetOwner(); + output->Add(out); + + THStack* dndetaRings = new THStack("all", "dN/d#eta per ring"); + const char* names[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 }; + const char** ptr = names; + + while (*ptr) { + TList* thisList = new TList; + thisList->SetOwner(); + thisList->SetName(*ptr); + out->Add(thisList); + + TH2D* h = static_cast(list->FindObject(Form("%s_cache", *ptr))); + if (!h) { + AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName())); + ptr++; + continue; + } + TH2D* copy = static_cast(h->Clone("sum")); + copy->SetDirectory(0); + thisList->Add(copy); + + TH1D* norm =static_cast(h->ProjectionX("norm", 0, 0, "")); + for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { + for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { + Double_t c = copy->GetBinContent(i, j); + Double_t e = copy->GetBinError(i, j); + Double_t a = norm->GetBinContent(i); + copy->SetBinContent(i, j, a <= 0 ? 0 : c / a); + copy->SetBinError(i, j, a <= 0 ? 0 : e / a); + } + } + + TH1D* res =static_cast(copy->ProjectionX("dndeta",1, + h->GetNbinsY(),"e")); + TH1D* proj =static_cast(h->ProjectionX("proj",1,h->GetNbinsY(),"e")); + res->SetTitle(*ptr); + res->Scale(1., "width"); + copy->Scale(1., "width"); + proj->Scale(1. / norm->GetMaximum(), "width"); + norm->Scale(1. / norm->GetMaximum()); + + res->SetMarkerStyle(style); + norm->SetDirectory(0); + res->SetDirectory(0); + proj->SetDirectory(0); + thisList->Add(norm); + thisList->Add(res); + thisList->Add(proj); + dndetaRings->Add(res); + ptr++; + } + out->Add(dndetaRings); +} + //____________________________________________________________________ void AliForwardMultiplicityBase::Print(Option_t* option) const diff --git a/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.h b/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.h index 270b1369e4e..558c6ca3014 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.h +++ b/PWG2/FORWARD/analysis2/AliForwardMultiplicityBase.h @@ -303,7 +303,18 @@ protected: * */ virtual void MarkEventForStore() const; - + /** + * Make Ring @f$ dN/d\eta @f$ histogram and a stack + * + * @param input List with summed signals + * @param output Output list + * @param stackName Stack name + */ + virtual void MakeRingdNdeta(const TList* input, + const char* inName, + TList* output, + const char* outName, + Int_t style=20) const; Bool_t fEnableLowFlux;// Whether to use low-flux specific code Bool_t fFirstEvent; // Whether the event is the first seen private: diff --git a/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx b/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx index ca348400ace..9b445725f31 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.cxx @@ -34,6 +34,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask() fESDFMD(), fHistos(), fAODFMD(), + fRingSums(), fEventInspector(), fSharingFilter(), fDensityCalculator(), @@ -53,6 +54,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name) fESDFMD(), fHistos(), fAODFMD(false), + fRingSums(), fEventInspector("event"), fSharingFilter("sharing"), fDensityCalculator("density"), @@ -76,6 +78,7 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplic fESDFMD(o.fESDFMD), fHistos(o.fHistos), fAODFMD(o.fAODFMD), + fRingSums(o.fRingSums), fEventInspector(o.fEventInspector), fSharingFilter(o.fSharingFilter), fDensityCalculator(o.fDensityCalculator), @@ -115,6 +118,7 @@ AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o) fHistCollector = o.fHistCollector; fHistos = o.fHistos; fAODFMD = o.fAODFMD; + fRingSums = o.fRingSums; fList = o.fList; return *this; @@ -152,13 +156,31 @@ AliForwardMultiplicityTask::InitializeSubs() fHistos.Init(*pe); fAODFMD.Init(*pe); + fRingSums.Init(*pe); fHData = static_cast(fAODFMD.GetHistogram().Clone("d2Ndetadphi")); fHData->SetStats(0); fHData->SetDirectory(0); fList->Add(fHData); + TList* rings = new TList; + rings->SetName("ringSums"); + rings->SetOwner(); + fList->Add(rings); + + rings->Add(fRingSums.Get(1, 'I')); + rings->Add(fRingSums.Get(2, 'I')); + rings->Add(fRingSums.Get(2, 'O')); + rings->Add(fRingSums.Get(3, 'I')); + rings->Add(fRingSums.Get(3, 'O')); + fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I')); + fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I')); + fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O')); + fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I')); + fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O')); + fEventInspector.Init(*pv); + fSharingFilter.Init(); fDensityCalculator.Init(*pe); fCorrections.Init(*pe); fHistCollector.Init(*pv,*pe); @@ -215,13 +237,14 @@ AliForwardMultiplicityTask::UserExec(Option_t*) fESDFMD.Clear(); fAODFMD.Clear(); - Bool_t lowFlux = kFALSE; - UInt_t triggers = 0; - UShort_t ivz = 0; - Double_t vz = 0; - Double_t cent = -1; - UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, - ivz, vz, cent); + Bool_t lowFlux = kFALSE; + UInt_t triggers = 0; + UShort_t ivz = 0; + Double_t vz = 0; + Double_t cent = -1; + UShort_t nClusters = 0; + UInt_t found = fEventInspector.Process(esd, triggers, lowFlux, + ivz, vz, cent, nClusters); if (found & AliFMDEventInspector::kNoEvent) return; if (found & AliFMDEventInspector::kNoTriggers) return; @@ -231,6 +254,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*) fAODFMD.SetSNN(fEventInspector.GetEnergy()); fAODFMD.SetSystem(fEventInspector.GetCollisionSystem()); fAODFMD.SetCentrality(cent); + fAODFMD.SetNClusters(nClusters); MarkEventForStore(); if (found & AliFMDEventInspector::kNoSPD) return; @@ -267,7 +291,8 @@ AliForwardMultiplicityTask::UserExec(Option_t*) return; } - if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) { + if (!fHistCollector.Collect(fHistos, fRingSums, + ivz, fAODFMD.GetHistogram())) { AliWarning("Histogram collector failed"); return; } @@ -329,6 +354,8 @@ AliForwardMultiplicityTask::Terminate(Option_t*) list->Add(dNdeta); list->Add(norm); + MakeRingdNdeta(list, "ringSums", list, "ringResults"); + fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral())); fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral())); diff --git a/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.h b/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.h index 1d784b82c38..1bbd5e3a265 100644 --- a/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.h +++ b/PWG2/FORWARD/analysis2/AliForwardMultiplicityTask.h @@ -179,6 +179,7 @@ protected: AliESDFMD fESDFMD; // Sharing corrected ESD object AliForwardUtil::Histos fHistos; // Cache histograms AliAODForwardMult fAODFMD; // Output object + AliForwardUtil::Histos fRingSums; // Cache histograms AliFMDEventInspector fEventInspector; // Algorithm AliFMDSharingFilter fSharingFilter; // Algorithm diff --git a/PWG2/FORWARD/analysis2/AliForwardUtil.cxx b/PWG2/FORWARD/analysis2/AliForwardUtil.cxx index e3d29b89ca8..c14619d8efb 100644 --- a/PWG2/FORWARD/analysis2/AliForwardUtil.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardUtil.cxx @@ -842,6 +842,7 @@ AliForwardUtil::RingHistos::DefineOutputList(TList* d) const // if (!d) return 0; TList* list = new TList; + list->SetOwner(); list->SetName(fName.Data()); d->Add(list); return list; diff --git a/PWG2/FORWARD/analysis2/AliForwardUtil.h b/PWG2/FORWARD/analysis2/AliForwardUtil.h index d7f103574f0..59d623d7cce 100644 --- a/PWG2/FORWARD/analysis2/AliForwardUtil.h +++ b/PWG2/FORWARD/analysis2/AliForwardUtil.h @@ -43,7 +43,7 @@ public: static Color_t RingColor(UShort_t d, Char_t r) { return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue)) - + ((r == 'I' || r == 'i') ? 2 : -2)); + + ((r == 'I' || r == 'i') ? 2 : -3)); } //================================================================== /** @@ -606,6 +606,7 @@ public: { return AliForwardUtil::RingColor(fDet, fRing); } + const char* GetName() const { return fName.Data(); } UShort_t fDet; // Detector Char_t fRing; // Ring TString fName; // Name diff --git a/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx b/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx index 0b9484d4541..580b7014ecd 100644 --- a/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx +++ b/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.cxx @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -99,6 +100,57 @@ AliForwarddNdetaTask::UserExec(Option_t* option) bin->ProcessPrimary(forward, fTriggerMask, fVtxMin, fVtxMax, primary); } +//________________________________________________________________________ +void +AliForwarddNdetaTask::Terminate(Option_t *option) +{ + // + // Called at end of event processing.. + // + // This is called once in the master + // + // Parameters: + // option Not used + AliBasedNdetaTask::Terminate(option); + + THStack* truth = new THStack("dndetaTruth", "dN/d#eta MC Truth"); + THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin), + "dN/d#eta MC Truth"); + + TIter next(fListOfCentralities); + CentralityBin* bin = 0; + while ((bin = static_cast(next()))) { + if (fCentAxis && bin->IsAllBin()) continue; + + TList* results = bin->GetResults(); + if (!results) continue; + + TH1* dndeta = static_cast(results->FindObject("dndetaTruth")); + TH1* dndetaRebin = + static_cast(results->FindObject(Form("dndetaTruth_rebin%02d", + fRebin))); + if (dndeta) truth->Add(dndeta); + if (dndetaRebin) truthRebin->Add(dndetaRebin); + } + // If available output rebinned stack + if (!truth->GetHists() || + truth->GetHists()->GetEntries() <= 0) { + AliWarning("No MC truth histograms found"); + delete truth; + truth = 0; + } + if (truth) fOutput->Add(truth); + + // If available output rebinned stack + if (!truthRebin->GetHists() || + truthRebin->GetHists()->GetEntries() <= 0) { + AliWarning("No rebinned MC truth histograms found"); + delete truthRebin; + truthRebin = 0; + } + if (truthRebin) fOutput->Add(truthRebin); + +} //____________________________________________________________________ TH2D* @@ -132,8 +184,8 @@ void AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult* forward, Int_t triggerMask, - Double_t /*vzMin*/, - Double_t /*vzMax*/, + Double_t vzMin, + Double_t vzMax, const TH2D* primary) { // Check the centrality class unless this is the 'all' bin @@ -152,13 +204,18 @@ AliForwarddNdetaTask::CentralityBin::ProcessPrimary(const AliAODForwardMult* // translate real trigger mask to MC trigger mask Int_t mask = AliAODForwardMult::kB; - if (triggerMask == AliAODForwardMult::kNSD) - mask = AliAODForwardMult::kMCNSD; + if (triggerMask == AliAODForwardMult::kNSD) { + mask ^= AliAODForwardMult::kNSD; + mask = AliAODForwardMult::kMCNSD; + } // Now use our normal check, but with the new mask, except - if (!forward->CheckEvent(mask, -1000, -1000, 0, 0, 0)) return; + vzMin = vzMax = -10000; // ignore vertex + if (!forward->CheckEvent(mask, vzMin, vzMax, 0, 0, 0)) return; fSumPrimary->Add(primary); + Int_t n = fSumPrimary->GetBinContent(0,0); + fSumPrimary->SetBinContent(0,0, ++n); } //________________________________________________________________________ @@ -174,7 +231,6 @@ AliForwarddNdetaTask::CentralityBin::End(TList* sums, Bool_t corrEmpty, Bool_t cutEdges, Int_t triggerMask, - Int_t color, Int_t marker) { AliInfo(Form("In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d", @@ -183,28 +239,73 @@ AliForwarddNdetaTask::CentralityBin::End(TList* sums, shapeCorr, trigEff, symmetrice, rebin, rootProj, corrEmpty, cutEdges, - triggerMask, color, marker); + triggerMask, marker); fSumPrimary = static_cast(fSums->FindObject("truth")); if (!fSumPrimary) { AliWarning("No MC truth histogram found"); - fSums->ls(); - return; } - Int_t n = (triggerMask == AliAODForwardMult::kNSD ? - Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) : - Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll))); - AliInfo(Form("Normalising MC truth to %d", n)); - - TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1, + else { +#if 0 + Int_t n = fSumPrimary->GetBinContent(0,0); +#else + Int_t n = (triggerMask == AliAODForwardMult::kNSD ? + Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) : + Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll))); +#endif + AliInfo(Form("Normalising MC truth to %d", n)); + + TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",1, fSumPrimary->GetNbinsY(),"e"); - dndetaTruth->Scale(1./n, "width"); + dndetaTruth->SetDirectory(0); + dndetaTruth->Scale(1./n, "width"); - SetHistogramAttributes(dndetaTruth, color-2, marker+4, "Monte-Carlo truth"); - fOutput->Add(dndetaTruth); - fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges)); + SetHistogramAttributes(dndetaTruth, GetColor(), 30, "Monte-Carlo truth"); + + fOutput->Add(dndetaTruth); + fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges)); + } + + if (!IsAllBin()) return; + TFile* file = TFile::Open("forward.root", "READ"); + if (!file) return; + + TList* forward = static_cast(file->Get("Forward")); + if (!forward) { + AliError("List Forward not found in forward.root"); + return; + } + TList* rings = static_cast(forward->FindObject("ringResults")); + if (!rings) { + AliError("List ringResults not found in forward.root"); + return; + } + THStack* res = static_cast(rings->FindObject("all")); + if (!res) { + AliError(Form("Stack all not found in %s", rings->GetName())); + return; + } + if (!fTriggers) { + AliError("Triggers histogram not set"); + return; + } + Double_t ntotal = 0; + Double_t epsilonT = trigEff; + // TEMPORARY FIX + if (triggerMask == AliAODForwardMult::kNSD) { + // This is a local change + epsilonT = 0.92; + AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT)); + } + AliInfo("Adding per-ring histograms to output"); + Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal); + TIter next(res->GetHists()); + TH1* hist = 0; + while ((hist = static_cast(next()))) hist->Scale(scaler); + res->SetName("dndetaRings"); + fOutput->Add(res); } //________________________________________________________________________ diff --git a/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h b/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h index 0859946f21d..309e4b0abba 100644 --- a/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h +++ b/PWG2/FORWARD/analysis2/AliForwarddNdetaTask.h @@ -51,6 +51,14 @@ public: * @param option Not used */ virtual void UserExec(Option_t* option); + /** + * Called at end of event processing. + * + * This is called once in the master + * + * @param option Not used + */ + virtual void Terminate(Option_t* option); protected: /** * Copy constructor @@ -162,7 +170,6 @@ protected: Bool_t corrEmpty, Bool_t cutEdges, Int_t triggerMask, - Int_t color, Int_t marker); protected: TH2D* fSumPrimary; // Sum of primary histograms diff --git a/PWG2/FORWARD/analysis2/DrawdNdeta.C b/PWG2/FORWARD/analysis2/DrawdNdeta.C index c6afac18bff..33086534ac6 100644 --- a/PWG2/FORWARD/analysis2/DrawdNdeta.C +++ b/PWG2/FORWARD/analysis2/DrawdNdeta.C @@ -3,11 +3,13 @@ * @author Christian Holm Christensen * @date Wed Mar 23 14:07:10 2011 * - * @brief Script to visualise the dN/deta + * @brief Script to visualise the dN/deta for pp and PbPb * * This script is independent of any AliROOT code - and should stay * that way. * + * The script is very long - sigh - the joy of drawing + * things nicely in ROOT * * @ingroup pwg2_forward_dndeta */ @@ -377,10 +379,12 @@ struct dNdetaDrawer sys,snn,trg, centLow,centHigh,onlya)); if (!ret) { +#if 0 Warning("FetchResults", "No other data found for sys=%d, sNN=%d, " "trigger=%d %d%%-%d%% central %s", sys, snn, trg, centLow, centHigh, onlya ? " (ALICE results)" : "all"); +#endif return 0; } thisOther = reinterpret_cast(ret); @@ -410,22 +414,16 @@ struct dNdetaDrawer Error("FetchResults", "Couldn't find list 'all' in %s", list->GetName()); else - FetchResults(all, name, FetchOthers(0,0), -1, 0, max, rmax, amax); + FetchResults(all, name, FetchOthers(0,0), -1000, 0, max, rmax, amax); return; } - Int_t nCol = gStyle->GetNumberOfColors(); for (UShort_t i = 0; i < n; i++) { - UShort_t centLow = fCentAxis->GetXbins()->At(i); - UShort_t centHigh = fCentAxis->GetXbins()->At(i+1); + UShort_t centLow = fCentAxis->GetBinLowEdge(i+1); + UShort_t centHigh = fCentAxis->GetBinUpEdge(i+1); TString lname = Form("cent%03d_%03d", centLow, centHigh); TList* thisCent = static_cast(list->FindObject(lname)); - - Float_t fc = (centLow+double(centHigh-centLow)/2) / 100; - Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5)); - Int_t col = gStyle->GetColorPalette(icol); - Info("FetchResults","Centrality %d-%d color index %d (=%d*%f) -> %d", - centLow, centHigh, icol, nCol, fc, col); + Int_t col = GetCentralityColor(i+1); TString centTxt = Form("%3d%%-%3d%% central", centLow, centHigh); if (!thisCent) @@ -437,6 +435,18 @@ struct dNdetaDrawer } } //__________________________________________________________________ + Int_t GetCentralityColor(Int_t bin) const + { + UShort_t centLow = fCentAxis->GetBinLowEdge(bin); + UShort_t centHigh = fCentAxis->GetBinUpEdge(bin); + Float_t fc = (centLow+double(centHigh-centLow)/2) / 100; + Int_t nCol = gStyle->GetNumberOfColors(); + Int_t icol = TMath::Min(nCol-1,int(fc * nCol + .5)); + Int_t col = gStyle->GetColorPalette(icol); + //Info("GetCentralityColor","%3d: %3d-%3d -> %3d",bin,centLow,centHigh,col); + return col; + } + //__________________________________________________________________ void SetAttributes(TH1* h, Int_t color) { if (!h) return; @@ -457,8 +467,9 @@ struct dNdetaDrawer //__________________________________________________________________ void ModifyTitle(TNamed* h, const char* centTxt) { - if (!centTxt || !h) return; - h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt)); + return; + // if (!centTxt || !h) return; + // h->SetTitle(Form("%s, %s", h->GetTitle(), centTxt)); } //__________________________________________________________________ @@ -491,11 +502,17 @@ struct dNdetaDrawer TH1* tracks = FetchResult(list, "tracks"); if (tracks) tracks->SetTitle("ALICE Tracks"); SetAttributes(dndeta, color); - SetAttributes(dndetaMC, color+2); + SetAttributes(dndetaMC, fCentAxis ? color : color+2); SetAttributes(dndetaTruth,color); SetAttributes(dndetaSym, color); - SetAttributes(dndetaMCSym,color+2); - SetAttributes(tracks, color+3); + SetAttributes(dndetaMCSym,fCentAxis ? color : color+2); + SetAttributes(tracks, fCentAxis ? color : color+3); + if (dndetaMC && fCentAxis) + dndetaMC->SetMarkerStyle(dndetaMC->GetMarkerStyle()+2); + if (dndetaMCSym && fCentAxis) + dndetaMCSym->SetMarkerStyle(dndetaMCSym->GetMarkerStyle()+2); + if (dndetaTruth && fCentAxis) + dndetaTruth->SetMarkerStyle(34); ModifyTitle(dndeta, centTxt); ModifyTitle(dndetaMC, centTxt); ModifyTitle(dndetaTruth,centTxt); @@ -508,7 +525,14 @@ struct dNdetaDrawer max = TMath::Max(max, AddHistogram(fResults, dndetaMC, dndetaMCSym)); max = TMath::Max(max, AddHistogram(fResults, dndeta, dndetaSym)); max = TMath::Max(max, AddHistogram(fResults, tracks)); - + + THStack* rings = static_cast(list->FindObject("dndetaRings")); + if (rings) { + TIter next(rings->GetHists()); + TH1* hist = 0; + while ((hist = static_cast(next()))) + max = TMath::Max(max, AddHistogram(fResults, hist)); + } // Info("FetchResults", "Got %p, %p, %p from %s with name %s, max=%f", // dndeta, dndetaMC, dndetaTruth, list->GetName(), name, max); @@ -540,11 +564,13 @@ struct dNdetaDrawer fRatios->Add(Ratio(dndeta, tracks, rmax)); } + if (dndetaMC) { + fRatios->Add(Ratio(dndeta, dndetaMC, rmax)); + fRatios->Add(Ratio(dndetaSym, dndetaMCSym, rmax)); + } if (dndetaTruth) { fRatios->Add(Ratio(dndeta, dndetaTruth, rmax)); fRatios->Add(Ratio(dndetaSym, dndetaTruth, rmax)); - fRatios->Add(Ratio(dndetaMC, dndetaTruth, rmax)); - fRatios->Add(Ratio(dndetaMCSym, dndetaTruth, rmax)); } } //__________________________________________________________________ @@ -580,12 +606,6 @@ struct dNdetaDrawer c->SetBorderSize(0); c->SetBorderMode(0); -#if 1 - Info("Plot", "y1=%f, y2=%f, y3=%f extra: %s %s", - y1, y2, y2, (fShowRatios ? "ratios" : ""), - (fShowLeftRight ? "right/left" : "")); - Info("Plot", "Maximum is %f", max); -#endif PlotResults(max, y1); c->cd(); @@ -628,27 +648,54 @@ struct dNdetaDrawer Double_t x1, Double_t y1, Double_t x2, Double_t y2) { TLegend* l = new TLegend(x1,y1,x2,y2); - l->SetNColumns(2); + l->SetNColumns(fCentAxis ? 1 : 2); l->SetFillColor(0); l->SetFillStyle(0); l->SetBorderSize(0); l->SetTextFont(132); TIter next(stack->GetHists()); - TObject* hist = 0; - while ((hist = next())) { + TH1* hist = 0; + TObjArray unique; + unique.SetOwner(); + while ((hist = static_cast(next()))) { TString n(hist->GetTitle()); if (n.Contains("mirrored")) continue; - l->AddEntry(hist, hist->GetTitle(), "pl"); + if (unique.FindObject(n.Data())) continue; + TObjString* s = new TObjString(n); + s->SetUniqueID(((hist->GetMarkerStyle() & 0xFFFF) << 16) | + ((hist->GetMarkerColor() & 0xFFFF) << 0)); + unique.Add(s); + // l->AddEntry(hist, hist->GetTitle(), "pl"); } if (mg) { TIter nexto(mg->GetListOfGraphs()); - while ((hist = nexto())) { - TString n(hist->GetTitle()); + TGraph* g = 0; + while ((g = static_cast(nexto()))) { + TString n(g->GetTitle()); if (n.Contains("mirrored")) continue; - l->AddEntry(hist, hist->GetTitle(), "pl"); + if (unique.FindObject(n.Data())) continue; + TObjString* s = new TObjString(n); + s->SetUniqueID(((g->GetMarkerStyle() & 0xFFFF) << 16) | + ((g->GetMarkerColor() & 0xFFFF) << 0)); + unique.Add(s); + // l->AddEntry(hist, hist->GetTitle(), "pl"); } } + unique.ls(); + TIter nextu(&unique); + TObject* s = 0; + Int_t i = 0; + while ((s = nextu())) { + TLegendEntry* dd = l->AddEntry(Form("data%2d", i++), + s->GetName(), "lp"); + Int_t style = (s->GetUniqueID() >> 16) & 0xFFFF; + Int_t color = (s->GetUniqueID() >> 0) & 0xFFFF; + dd->SetLineColor(kBlack); + if (fCentAxis) dd->SetMarkerColor(kBlack); + else dd->SetMarkerColor(color); + dd->SetMarkerStyle(style); + } TLegendEntry* d1 = l->AddEntry("d1", "Data", "lp"); d1->SetLineColor(kBlack); d1->SetMarkerColor(kBlack); @@ -661,6 +708,38 @@ struct dNdetaDrawer l->Draw(); } //__________________________________________________________________ + /** + * Build centrality legend + * + * @param axis Stack to include + * @param x1 Lower X coordinate in the range [0,1] + * @param y1 Lower Y coordinate in the range [0,1] + * @param x2 Upper X coordinate in the range [0,1] + * @param y2 Upper Y coordinate in the range [0,1] + */ + void BuildCentLegend(Double_t x1, Double_t y1, Double_t x2, Double_t y2) + { + if (!fCentAxis) return; + + TLegend* l = new TLegend(x1,y1,x2,y2); + l->SetNColumns(1); + l->SetFillColor(0); + l->SetFillStyle(0); + l->SetBorderSize(0); + l->SetTextFont(132); + + Int_t n = fCentAxis->GetNbins(); + for (Int_t i = 1; i <= n; i++) { + Double_t low = fCentAxis->GetBinLowEdge(i); + Double_t upp = fCentAxis->GetBinUpEdge(i); + TLegendEntry* e = l->AddEntry(Form("dummy%02d", i), + Form("%3d%% - %3d%%", + int(low), int(upp)), "pl"); + e->SetMarkerColor(GetCentralityColor(i)); + } + l->Draw(); + } + //__________________________________________________________________ /** * Plot the results * @@ -703,15 +782,14 @@ struct dNdetaDrawer } // Make a legend in the main result pad - BuildLegend(fResults, fOthers, .15,p1->GetBottomMargin()+.01,.90,.35); -#if 0 - TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35); - l->SetNColumns(2); - l->SetFillColor(0); - l->SetFillStyle(0); - l->SetBorderSize(0); - l->SetTextFont(132); -#endif + BuildCentLegend(.12, 1-p1->GetTopMargin()-.01-.5, + .35, 1-p1->GetTopMargin()-.01-.1); + Double_t x1 = (fCentAxis ? .7 : .15); + Double_t x2 = (fCentAxis ? 1-p1->GetRightMargin()-.01: .90); + Double_t y1 = (fCentAxis ? .5: p1->GetBottomMargin()+.01); + Double_t y2 = (fCentAxis ? 1-p1->GetTopMargin()-.01-.15 : .35); + + BuildLegend(fResults, fOthers, x1, y1, x2, y2); // Put a title on top fTitle.ReplaceAll("@", " "); @@ -725,11 +803,14 @@ struct dNdetaDrawer TString eS; UShort_t snn = fSNNString->GetUniqueID(); const char* sys = fSysString->GetTitle(); + if (snn == 2750) snn = 2760; if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000); else eS = Form("%3dGeV", snn); - TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s", + TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s%s}=%s, %s", sys, + (fCentAxis ? "_{NN}" : ""), eS.Data(), + fCentAxis ? "by centrality" : fTrigString->GetTitle())); tt->SetNDC(); tt->SetTextFont(132); @@ -886,12 +967,17 @@ struct dNdetaDrawer // Make a legend - TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5); - l2->SetNColumns(2); - l2->SetFillColor(0); - l2->SetFillStyle(0); - l2->SetBorderSize(0); - l2->SetTextFont(132); + Double_t xx1 = (fCentAxis ? .7 : .15); + Double_t xx2 = (fCentAxis ? 1-p3->GetRightMargin()-.01 : .90); + Double_t yy1 = p3->GetBottomMargin()+.01; + Double_t yy2 = (fCentAxis ? 1-p3->GetTopMargin()-.01-.15 : .5); + BuildLegend(fLeftRight, 0, xx1, yy1, xx2, yy2); + // TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5); + // l2->SetNColumns(2); + // l2->SetFillColor(0); + // l2->SetFillStyle(0); + // l2->SetBorderSize(0); + // l2->SetTextFont(132); // Make a nice band from 0.9 to 1.1 TGraphErrors* band = new TGraphErrors(2); @@ -936,10 +1022,12 @@ struct dNdetaDrawer if (!list) return 0; TH1* ret = static_cast(list->FindObject(name)); +#if 0 if (!ret) { // all->ls(); Warning("GetResult", "Histogram %s not found", name); } +#endif return ret; } @@ -961,10 +1049,7 @@ struct dNdetaDrawer // Rebin if needed Rebin(hist); - // Info("AddHistogram", "Adding %s to %s", - // hist->GetName(), stack->GetName()); stack->Add(hist, option); - // stack->ls(); return hist->GetMaximum(); } //__________________________________________________________________ @@ -992,8 +1077,6 @@ struct dNdetaDrawer sym = Symmetrice(hist); stack->Add(sym, option); - // Info("AddHistogram", "Adding %s and %s to %s", - // hist->GetName(), sym->GetName(), stack->GetName()); return hist->GetMaximum(); } @@ -1239,32 +1322,34 @@ struct dNdetaDrawer { if (!o1 || !o2) return 0; - TH1* r = 0; + TH1* r = 0; + const TAttMarker* m1 = 0; + const TAttMarker* m2 = 0; const TH1* h1 = dynamic_cast(o1); if (h1) { - // Info("Ratio", "First is a TH1"); + m1 = h1; const TH1* h2 = dynamic_cast(o2); if (h2) { - // Info("Ratio", "Second is a TH1"); - r = RatioHH(h1,h2,max); + m2 = h2; + r = RatioHH(h1,h2); } else { const TGraph* g2 = dynamic_cast(o2); if (g2) { - // Info("Ratio", "Second os a TGraph"); - r = RatioHG(h1,g2,max); + m2 = g2; + r = RatioHG(h1,g2); } } } else { const TGraphAsymmErrors* g1 = dynamic_cast(o1); if (g1) { - // Info("Ratio", "First is a TGraphAsymmErrors"); + m1 = g1; const TGraphAsymmErrors* g2 = dynamic_cast(o2); if (g2) { - // Info("Ratio", "Second is a TGraphAsymmErrors"); - r = RatioGG(g1, g2, max); + m2 = g2; + r = RatioGG(g1, g2); } } } @@ -1278,9 +1363,16 @@ struct dNdetaDrawer delete r; r = 0; } - // for (Int_t bin = 1; bin <= r->GetNbinsX(); bin++) - // if (r->GetBinContent(bin) != 0) return r; - + if (r) { + r->SetMarkerStyle(m1->GetMarkerStyle()); + r->SetMarkerColor(m2->GetMarkerColor()); + r->SetMarkerSize(0.9*m1->GetMarkerSize()); + r->SetName(Form("%s_over_%s", o1->GetName(), o2->GetName())); + r->SetTitle(Form("%s / %s", o1->GetTitle(), o2->GetTitle())); + r->SetDirectory(0); + max = TMath::Max(RatioMax(r), max); + } + return r; } //__________________________________________________________________ @@ -1290,25 +1382,15 @@ struct dNdetaDrawer * * @param h Numerator * @param g Divisor - * @param max Maximum diviation from 1 * * @return h/g */ - TH1* RatioHG(const TH1* h, const TGraph* g, Double_t& max) const + TH1* RatioHG(const TH1* h, const TGraph* g) const { if (!h || !g) return 0; TH1* ret = static_cast(h->Clone("tmp")); - ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName())); - ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle())); ret->Reset(); - ret->SetMarkerStyle(h->GetMarkerStyle()); - ret->SetMarkerColor(g->GetMarkerColor()); - ret->SetMarkerSize(0.9*h->GetMarkerSize()); - // ret->SetMarkerStyle(g->GetMarkerStyle()); - // ret->SetMarkerColor(h->GetMarkerColor()); - // ret->SetMarkerSize(0.9*g->GetMarkerSize()); - ret->SetDirectory(0); Double_t xlow = g->GetX()[0]; Double_t xhigh = g->GetX()[g->GetN()-1]; if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; } @@ -1326,8 +1408,6 @@ struct dNdetaDrawer ret->SetBinContent(i, c / f); ret->SetBinError(i, h->GetBinError(i) / f); } - if (ret->GetEntries() > 0) - max = TMath::Max(RatioMax(ret), max); return ret; } //__________________________________________________________________ @@ -1336,25 +1416,14 @@ struct dNdetaDrawer * * @param h1 First histogram (numerator) * @param h2 Second histogram (denominator) - * @param max Maximum diviation from 1 * * @return h1 / h2 */ - TH1* RatioHH(const TH1* h1, const TH1* h2, Double_t& max) const + TH1* RatioHH(const TH1* h1, const TH1* h2) const { if (!h1 || !h2) return 0; - TH1* t1 = static_cast(h1->Clone(Form("%s_%s", - h1->GetName(), - h2->GetName()))); - t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle())); + TH1* t1 = static_cast(h1->Clone("tmp")); t1->Divide(h2); - // t1->SetMarkerColor(h1->GetMarkerColor()); - // t1->SetMarkerStyle(h2->GetMarkerStyle()); - t1->SetMarkerColor(h2->GetMarkerColor()); - t1->SetMarkerStyle(h1->GetMarkerStyle()); - t1->SetMarkerSize(0.9*h1->GetMarkerSize()); - t1->SetDirectory(0); - max = TMath::Max(RatioMax(t1), max); return t1; } //__________________________________________________________________ @@ -1363,12 +1432,11 @@ struct dNdetaDrawer * * @param g1 Numerator * @param g2 Denominator - * @param max Maximum diviation from 1 * * @return g1 / g2 in a histogram */ TH1* RatioGG(const TGraphAsymmErrors* g1, - const TGraphAsymmErrors* g2, Double_t& max) const + const TGraphAsymmErrors* g2) const { Int_t nBins = g1->GetN(); TArrayF bins(nBins+1); @@ -1383,22 +1451,13 @@ struct dNdetaDrawer if (i == 0) dx = dxi; else if (dxi != dx) dx = 0; } - TString name(Form("%s_%s", g1->GetName(), g2->GetName())); - TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle())); TH1* h = 0; if (dx != 0) { - h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]); + h = new TH1F("tmp", "tmp", nBins, bins[0], bins[nBins]); } else { - h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray); + h = new TH1F("tmp", "tmp", nBins, bins.fArray); } - h->SetMarkerStyle(g1->GetMarkerStyle()); - h->SetMarkerColor(g2->GetMarkerColor()); - h->SetMarkerSize(0.9*g1->GetMarkerSize()); - // h->SetMarkerStyle(g2->GetMarkerStyle()); - // h->SetMarkerColor(g1->GetMarkerColor()); - // h->SetMarkerSize(0.9*g2->GetMarkerSize()); - h->SetDirectory(0); Double_t low = g2->GetX()[0]; Double_t high = g2->GetX()[g2->GetN()-1]; @@ -1413,7 +1472,6 @@ struct dNdetaDrawer h->SetBinContent(i+1, c1 / c2); h->SetBinError(i+1, e1 / c2); } - max = TMath::Max(RatioMax(h), max); return h; } /* @} */ @@ -1543,7 +1601,6 @@ void UpdateRange(dNdetaDrawer::RangeParam* p) Int_t last = p->fMasterAxis->GetLast(); Double_t x1 = p->fMasterAxis->GetBinCenter(first); Double_t x2 = p->fMasterAxis->GetBinCenter(last); - //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2); if (p->fSlave1Axis) { Int_t i1 = p->fSlave1Axis->FindBin(x1); diff --git a/PWG2/FORWARD/analysis2/ForwardAODConfig.C b/PWG2/FORWARD/analysis2/ForwardAODConfig.C index 89528561dc1..50a8e15fe92 100644 --- a/PWG2/FORWARD/analysis2/ForwardAODConfig.C +++ b/PWG2/FORWARD/analysis2/ForwardAODConfig.C @@ -28,6 +28,22 @@ ForwardAODConfig(AliForwardMultiplicityBase* task) // Whether to enable low flux specific code task->SetEnableLowFlux(kFALSE); + // Would like to use dynamic cast but CINT interprets that as a + // static cast - sigh! + Bool_t mc = false; + if (task->IsA() == AliForwardMCMultiplicityTask::Class()) + mc = true; + +#if 0 + if (mc) { + AliForwardMCMultiplicityTask* mcTask = + static_cast(task); + mcTask->SetOnlyPrimary(true); + } +#endif + Double_t nXi = mc ? 1 : .5; + Bool_t includeSigma = true; + // --- Event inspector --------------------------------------------- // Set the number of SPD tracklets for which we consider the event a // low flux event @@ -37,17 +53,29 @@ ForwardAODConfig(AliForwardMultiplicityBase* task) // --- Sharing filter ---------------------------------------------- // Set the low cut used for sharing - overrides settings in eloss fits - task->GetSharingFilter().SetLowCut(0.3); + task->GetSharingFilter().SetLowCut(0.15); // Set the number of xi's (width of landau peak) to stop at - task->GetSharingFilter().SetNXi(1); + task->GetSharingFilter().SetNXi(nXi); + // Set whether or not to include sigma in cut + task->GetSharingFilter().SetIncludeSigma(includeSigma); + // Enable use of angle corrected signals in the algorithm + task->GetSharingFilter().SetUseAngleCorrectedSignals(true); // --- Density calculator // Set the maximum number of particle to try to reconstruct - task->GetDensityCalculator().SetMaxParticles(3); + task->GetDensityCalculator().SetMaxParticles(10); // Wet whether to use poisson statistics to estimate N_ch task->GetDensityCalculator().SetUsePoisson(false); // Set the lower multiplicity cut. Overrides setting in energy loss fits. - task->GetDensityCalculator().SetMultCut(0.3); //was 0.3 + task->GetDensityCalculator().SetMultCut(-1); //was 0.3 + // Set the lower per-ring multiplicity cuts + task->GetDensityCalculator().SetMultCuts(-1,-1,-1,-1,-1); + // USe this many times xi+sigma below MPV + task->GetDensityCalculator().SetNXi(nXi); + // Set whether or not to include sigma in cut + task->GetDensityCalculator().SetIncludeSigma(includeSigma); + // Set whether or not to use the phi acceptance + task->GetDensityCalculator().SetUsePhiAcceptance(true); // --- Corrector --------------------------------------------------- // Whether to use the secondary map correction @@ -82,7 +110,7 @@ ForwardAODConfig(AliForwardMultiplicityBase* task) // Set the overall debug level (1: some output, 3: a lot of output) task->SetDebug(0); // Set the debug level of a single algorithm - // task->GetEventInspector().SetDebug(4); + task->GetSharingFilter().SetDebug(0); // --- Set limits on fits the energy ------------------------------- // Maximum relative error on parameters diff --git a/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C b/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C new file mode 100644 index 00000000000..4e5481c4934 --- /dev/null +++ b/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C @@ -0,0 +1,851 @@ +#ifdef BUILD +#include "AliAnalysisTaskSE.h" +#include "AliAnalysisManager.h" +#include "AliAnalysisDataContainer.h" +#include "AliMCEvent.h" +#include "AliESDEvent.h" +#include "AliStack.h" +#include "AliMultiplicity.h" +#include "AliFMDMCEventInspector.h" +#include "AliAODForwardMult.h" +#include "AliLog.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//==================================================================== +/** + * Task to evaluate trigger bias in pp + * + */ +class EvaluateTrigger : public AliAnalysisTaskSE +{ +public: + enum { + kNone = 0x0, + kESD = 0x1, + kMC = 0x2 + }; + /** + * Constructor + */ + EvaluateTrigger() + : AliAnalysisTaskSE(), + fInel(), + fInelGt0(), + fNSD(), + fInspector(), + fFirstEvent(true), + fData(0), + fTriggers(0), + fTrackletRequirement(kESD), + fVertexRequirement(kESD), + fVertexAxis(0, 0, 0), + fVertexESD(0), + fVertexMC(0), + fM(0) + {} + /** + * Constructor + */ + EvaluateTrigger(const char* /*name*/) + : AliAnalysisTaskSE("evaluateTrigger"), + fInel(AliAODForwardMult::kInel), + fInelGt0(AliAODForwardMult::kInelGt0), + fNSD(AliAODForwardMult::kNSD), + fInspector("eventInspector"), + fFirstEvent(true), + fData(0), + fTriggers(0), + fTrackletRequirement(kESD), + fVertexRequirement(kESD), + fVertexAxis(10, -10, 10), + fVertexESD(0), + fVertexMC(0), + fM(0) + { + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + } + void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; } + void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; } + void SetVertexAxis(Int_t nBins, Double_t low, Double_t high) + { + fVertexAxis.Set(nBins, low, high); + } + /** + * Intialize + */ + void Init() {} + /** + * Create output objects + */ + void UserCreateOutputObjects() + { + fList = new TList; + fList->SetOwner(); + fList->SetName(GetName()); + + Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 1000 }; + Int_t nM = 10; + TAxis mAxis(nM, mb); + TAxis eAxis(200, -4, 6); + TAxis pAxis(40, 0, 2*TMath::Pi()); + + fData = new TH2D("data", "Cache", + eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), + pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax()); + fData->SetDirectory(0); + fData->SetXTitle("#eta"); + fData->SetYTitle("#varphi [radians]"); + fData->SetZTitle("N_{ch}(#eta,#varphi)"); + fData->Sumw2(); + + fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}",nM+1, -0.5, nM+.5); + fM->SetXTitle("N_{ch}|_{|#eta|<1}"); + fM->SetYTitle("Events"); + fM->SetFillColor(kRed+1); + fM->SetFillStyle(3001); + fM->SetDirectory(0); + fList->Add(fM); + + for (Int_t i = 0; i <= nM; i++) { + TString lbl; + if (i == 0) lbl = "all"; + else if (i == nM) lbl = Form("%d+",int(mAxis.GetBinLowEdge(i)+.5)); + else lbl = Form("<%d",int(mAxis.GetBinUpEdge(i)+.5)); + fM->GetXaxis()->SetBinLabel(i+1, lbl); + } + + fTriggers = new TH1I("triggers", "Triggers", 6, -.5, 5.5); + fTriggers->SetDirectory(0); + fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)"); + fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)"); + fTriggers->GetXaxis()->SetBinLabel(3, "INEL>0 (MC)"); + fTriggers->GetXaxis()->SetBinLabel(4, "INEL>0 (ESD)"); + fTriggers->GetXaxis()->SetBinLabel(5, "NSD (MC)"); + fTriggers->GetXaxis()->SetBinLabel(6, "NSD (ESD)"); + fTriggers->SetFillColor(kYellow+1); + fTriggers->SetFillStyle(3001); + fList->Add(fTriggers); + + fVertexESD = new TH1D("vertexESD", "ESD vertex distribution", + fVertexAxis.GetNbins(), + fVertexAxis.GetXmin(), + fVertexAxis.GetXmax()); + fVertexESD->SetDirectory(0); + fVertexESD->SetFillColor(kRed+1); + fVertexESD->SetFillStyle(3001); + fVertexESD->SetXTitle("v_{z} [cm]"); + fVertexESD->SetYTitle("P(v_{z})"); + fList->Add(fVertexESD); + + fVertexMC = new TH1D("vertexMC", "MC vertex distribution", + fVertexAxis.GetNbins(), + fVertexAxis.GetXmin(), + fVertexAxis.GetXmax()); + fVertexMC->SetDirectory(0); + fVertexMC->SetFillColor(kBlue+1); + fVertexMC->SetFillStyle(3001); + fVertexMC->SetXTitle("v_{z} [cm]"); + fVertexMC->SetYTitle("P(v_{z})"); + fList->Add(fVertexMC); + + fInel.CreateObjects(fList, fM, fData); + fInelGt0.CreateObjects(fList, fM, fData); + fNSD.CreateObjects(fList, fM, fData); + + + fInspector.DefineOutput(fList); + fInspector.Init(fVertexAxis); + + PostData(1, fList); + } + /** + * Event processing + */ + void UserExec(Option_t*) + { + // Get the input data - MC event + AliMCEvent* mcEvent = MCEvent(); + if (!mcEvent) { + AliWarning("No MC event found"); + return; + } + + // Get the input data - ESD event + AliESDEvent* esd = dynamic_cast(InputEvent()); + if (!esd) { + AliWarning("No ESD event found for input event"); + return; + } + + if (fFirstEvent && esd->GetESDRun()) { + fInspector.ReadRunDetails(esd); + + AliInfo(Form("Initializing with parameters from the ESD:\n" + " AliESDEvent::GetBeamEnergy() ->%f\n" + " AliESDEvent::GetBeamType() ->%s\n" + " AliESDEvent::GetCurrentL3() ->%f\n" + " AliESDEvent::GetMagneticField()->%f\n" + " AliESDEvent::GetRunNumber() ->%d\n", + esd->GetBeamEnergy(), + esd->GetBeamType(), + esd->GetCurrentL3(), + esd->GetMagneticField(), + esd->GetRunNumber())); + + fFirstEvent = false; + } + + // Get the particle stack + AliStack* stack = mcEvent->Stack(); + + // Some variables + UInt_t triggers; // Trigger bits + Bool_t lowFlux; // Low flux flag + UShort_t iVz; // Vertex bin from ESD + Double_t vZ; // Z coordinate from ESD + Double_t cent; // Centrality + UShort_t iVzMc; // Vertex bin from MC + Double_t vZMc; // Z coordinate of IP vertex from MC + Double_t b; // Impact parameter + Int_t nPart; // Number of participants + Int_t nBin; // Number of binary collisions + Double_t phiR; // Reaction plane from MC + + // Process the data + Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent); + Int_t retMC = fInspector.ProcessMC(mcEvent, triggers, iVzMc, + vZMc, b, nPart, nBin, phiR); + + Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk; + Bool_t hasMCVtx = retMC == AliFMDEventInspector::kOk; + if (hasESDVtx) fVertexESD->Fill(vZ); + if (hasMCVtx) fVertexMC->Fill(vZMc); + + Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB); + Bool_t isMcNSD = (triggers & AliAODForwardMult::kMCNSD); + + Int_t mESD = 0; + const AliMultiplicity* spdmult = esd->GetMultiplicity(); + if (!spdmult) { + AliWarning("No SPD multiplicity"); + } + else { + // Check if we have one or more tracklets + // in the range -1 < eta < 1 to set the INEL>0 + // trigger flag. + Int_t n = spdmult->GetNumberOfTracklets(); + for (Int_t j = 0; j < n; j++) + if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++; + } + + // Reset cache + fData->Reset(); + Int_t mMC = 0; // Number of particles in |eta|<1 + + // Loop over all tracks + Int_t nTracks = mcEvent->GetNumberOfTracks(); + for (Int_t iTr = 0; iTr < nTracks; iTr++) { + AliMCParticle* particle = + static_cast(mcEvent->GetTrack(iTr)); + + // Check the returned particle + if (!particle) continue; + + // Check if this charged and a primary + Bool_t isCharged = particle->Charge() != 0; + Bool_t isPrimary = stack->IsPhysicalPrimary(iTr); + + if (!isCharged || !isPrimary) continue; + + + // Fill (eta,phi) of the particle into histograsm for b + Double_t eta = particle->Eta(); + Double_t phi = particle->Phi(); + + fData->Fill(eta, phi); + if (TMath::Abs(eta) <= 1) mMC++; + } + Int_t m = mESD; + if (fTrackletRequirement == kMC) m = mMC; + fM->Fill(m); + + bool isMcInelGt0 = isMcInel && (mMC > 0); + + bool hasVertex = true; + if (fVertexRequirement & kMC) hasVertex = hasVertex && hasMCVtx; + if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx; + + if (isMcInel) { + fTriggers->Fill(0); + bool triggered = (triggers & AliAODForwardMult::kInel); + if (triggered) fTriggers->Fill(1); + fInel.AddEvent(triggered, hasVertex, m, fData); + } + if (isMcInelGt0) { + fTriggers->Fill(2); + bool triggered = (triggers & AliAODForwardMult::kInelGt0); + if (triggered) fTriggers->Fill(3); + fInelGt0.AddEvent(triggered, hasVertex, m, fData); + } + if (isMcNSD) { + fTriggers->Fill(4); + bool triggered = (triggers & AliAODForwardMult::kNSD); + if (triggered) fTriggers->Fill(5); + fNSD.AddEvent(triggered, hasVertex, m, fData); + } + PostData(1, fList); + } + /** + * End of job processing + */ + void Terminate(Option_t*) + { + fList = dynamic_cast(GetOutputData(1)); + if (!fList) { + AliError(Form("No output list defined (%p)", GetOutputData(1))); + if (GetOutputData(1)) GetOutputData(1)->Print(); + return; + } + + + TList* output = new TList; + output->SetName(GetName()); + output->SetOwner(); + + fVertexMC = static_cast(fList->FindObject("vertexMC")); + fVertexESD = static_cast(fList->FindObject("vertexESD")); + fM = static_cast(fList->FindObject("m")); + if (fVertexMC) { + TH1D* vtxMC = static_cast(fVertexMC->Clone("vertexMC")); + vtxMC->SetDirectory(0); + vtxMC->Scale(1. / vtxMC->GetEntries()); + output->Add(vtxMC); + } + if (fVertexESD) { + TH1D* vtxESD = static_cast(fVertexESD->Clone("vertexESD")); + vtxESD->SetDirectory(0); + vtxESD->Scale(1. / vtxESD->GetEntries()); + output->Add(vtxESD); + } + if (fM) { + TH1D* m = static_cast(fM->Clone("m")); + m->SetDirectory(0); + m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)"); + m->Scale(1. / m->GetBinContent(1)); + output->Add(m); + } + + fInel.Finish(fList, output); + fInelGt0.Finish(fList, output); + fNSD.Finish(fList, output); + + PostData(2, output); + } + +protected: + //__________________________________________________________________ + /** + * Structure to hold per trigger type information + */ + struct TriggerType : public TNamed + { + //________________________________________________________________ + /** + * Structure to hold per multiplicity bin information + */ + struct MBin : public TNamed + { + TH2D* fTruth; + TH2D* fTriggered; + TH2D* fAccepted; + TH1I* fCounts; + UShort_t fLow; + UShort_t fHigh; + /** + * Constructor + */ + MBin() + : fTruth(0), fTriggered(0), fAccepted(0), + fCounts(0), fLow(0), fHigh(1000) {} + /** + * Constructor + * + * @param p Parent list + * @param low Low cut + * @param high High cut + * @param eAxis Eta axis + * @param pAxis Phi axis + */ + MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist) + : fTruth(0), + fTriggered(0), + fAccepted(0), + fCounts(0), + fLow(low), + fHigh(high) + { + if (low >= high) SetName("all"); + else SetName(Form("m%03d_%03d", fLow, fHigh)); + TList* l = new TList; + l->SetOwner(); + l->SetName(GetName()); + p->Add(l); + + fTruth = static_cast(dHist->Clone(("truth"))); + fTruth->SetTitle("MC truth"); + fTruth->SetDirectory(0); + fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)"); + fTruth->Sumw2(); + fTruth->Reset(); + l->Add(fTruth); + + fTriggered = static_cast(fTruth->Clone(("triggered"))); + fTriggered->SetTitle("Triggered"); + fTriggered->SetDirectory(0); + fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)"); + fTriggered->Sumw2(); + fTriggered->Reset(); + l->Add(fTriggered); + + fAccepted = static_cast(fTruth->Clone(("accepted"))); + fAccepted->SetTitle("Accepted"); + fAccepted->SetDirectory(0); + fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)"); + fAccepted->Sumw2(); + fAccepted->Reset(); + l->Add(fAccepted); + + fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5); + fCounts->SetDirectory(0); + fCounts->GetXaxis()->SetBinLabel(1, "Truth"); + fCounts->GetXaxis()->SetBinLabel(2, "Triggered"); + fCounts->GetXaxis()->SetBinLabel(3, "Accepted"); + fCounts->SetYTitle("# events"); + l->Add(fCounts); + } + /** + * Add event observation + * + * @param triggered Whether the event was triggered + * @param event Data for this event + */ + void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event) + { + fCounts->Fill(0); + fTruth->Add(event); + if (triggered) { + fCounts->Fill(1); + fTriggered->Add(event); + if (hasVtx) { + fCounts->Fill(2); + fAccepted->Add(event); + } + } + } + /** + * End of job processing + * + * @param p Parent list + * @param o Output parent list + * @param stack Stack of histograms + * + * @return Trigger efficiency + */ + Double_t Finish(const TList* p, TList* o, THStack* stack) + { + TList* l = dynamic_cast(p->FindObject(GetName())); + if (!l) { + Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName()); + return 0; + } + fTruth = static_cast(l->FindObject("truth")); + fTriggered = static_cast(l->FindObject("triggered")); + fAccepted = static_cast(l->FindObject("accepted")); + fCounts = static_cast(l->FindObject("counts")); + + Int_t nTruth = fCounts->GetBinContent(1); + Int_t nTriggered = fCounts->GetBinContent(2); + Int_t nAccepted = fCounts->GetBinContent(3); + Double_t eff = 0; + if (nTruth > 0) eff = double(nTriggered) / nTruth; + else if (nTriggered == nTruth) eff = 1; + + if (nTruth > 0) fTruth->Scale(1. / nTruth); + if (nTriggered > 0) fTriggered->Scale(1. / nTriggered); + if (nAccepted > 0) fAccepted->Scale(1. / nAccepted); + + if (fLow >= fHigh) + Info("Finish", "%-6s [all] E_X=N_T/N_X=%9d/%-9d=%f " + "E_V=N_A/N_T=%9d/%-9d=%f", + p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered, + (nTriggered > 0 ? double(nAccepted) / nTriggered: 0)); + else + Info("Finish", "%-6s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f " + "E_V=N_A/N_T=%9d/%-9d=%f", + p->GetName(), fLow, fHigh, nTriggered, nTruth, eff, + nAccepted, nTriggered, + (nTriggered > 0 ? double(nAccepted) / nTriggered: 0)); + + TList* out = new TList; + out->SetName(GetName()); + out->SetOwner(); + o->Add(out); + + out->Add(fTruth); + out->Add(fTriggered); + out->Add(new TParameter("eff", eff)); + + TH2D* bias = static_cast(fAccepted->Clone("bias")); + bias->Divide(fTruth); + bias->SetDirectory(0); + bias->SetZTitle("Trigger bias (accepted/truth)"); + out->Add(bias); + + TH1D* truth_px = static_cast(fTruth->ProjectionX("truth_eta")); + truth_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", + fLow, fHigh)); + truth_px->Scale(1. / fTruth->GetNbinsY()); + truth_px->SetDirectory(0); + truth_px->SetLineColor(kRed+1); + truth_px->SetMarkerColor(kRed+1); + truth_px->SetFillColor(kRed+1); + truth_px->SetFillStyle(0); + out->Add(truth_px); + + TH1D* triggered_px = + static_cast(fTriggered->ProjectionX("triggered_eta")); + triggered_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", + fLow, fHigh)); + triggered_px->Scale(1. / fTriggered->GetNbinsY()); + triggered_px->SetDirectory(0); + triggered_px->SetLineColor(kGreen+1); + triggered_px->SetMarkerColor(kGreen+1); + triggered_px->SetFillColor(kGreen+1); + triggered_px->SetFillStyle(0); + out->Add(triggered_px); + + TH1D* accepted_px = + static_cast(fAccepted->ProjectionX("accepted_eta")); + accepted_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", + fLow, fHigh)); + accepted_px->Scale(1. / fAccepted->GetNbinsY()); + accepted_px->SetLineColor(kBlue+1); + accepted_px->SetMarkerColor(kBlue+1); + accepted_px->SetFillColor(kBlue+1); + accepted_px->SetDirectory(0); + out->Add(accepted_px); + + THStack* data = new THStack("data", "Data distributions"); + data->Add(truth_px); + data->Add(triggered_px); + data->Add(accepted_px); + out->Add(data); + +#if 0 + TH1D* bias_px = static_cast(bias->ProjectionX("bias_eta")); + bias_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", + fLow, fHigh)); + bias_px->Scale(1. / bias->GetNbinsY()); +#else + TH1D* bias_px = static_cast(accepted_px->Clone("bias_px")); + bias_px->Divide(truth_px); + bias_px->SetYTitle("Trigger bias (triggered/truth)"); +#endif + bias_px->SetDirectory(0); + bias_px->SetMarkerStyle(20); + bias_px->SetFillStyle(0); + bias_px->SetMinimum(0); + out->Add(bias_px); + + stack->Add(bias_px); + + return eff; + } + ClassDef(MBin,1); + }; + + /** + * Constructor + */ + TriggerType() + : TNamed(), + fMask(0), + fM(0), + fBins(0) + {} + //--- Constructor ------------------------------------------------ + /** + * Constructor + * + * @param mask Trigger mask + */ + TriggerType(UShort_t mask) + : TNamed(AliAODForwardMult::GetTriggerString(mask), ""), + fMask(mask), + fM(0), + fBins(0) + { + } + /** + * Create objects + * + * @param list PArent list + * @param mAxis Multiplicity axis + * @param eAxis Eta axis + * @param pAxis Phi axis + */ + void CreateObjects(TList* list, + const TH1D* mHist, + const TH2D* dHist) + { + TList* ours = new TList; + ours->SetName(GetName()); + ours->SetOwner(); + list->Add(ours); + + fM = static_cast(mHist->Clone("m")); + fM->SetDirectory(0); + ours->Add(fM); + + fBins = new TObjArray; + fBins->AddAt(new MBin(ours, 0, 0, dHist), 0); + for (Int_t i = 1; i <= fM->GetNbinsX(); i++) { + Double_t low = fM->GetXaxis()->GetBinLowEdge(i); + Double_t high = fM->GetXaxis()->GetBinUpEdge(i); + + fBins->AddAt(new MBin(ours, low, high, dHist), i); + } + } + /** + * Find bin corresponding to m + * + * @param m Multiplicity + * + * @return Bin. + */ + MBin* FindBin(UShort_t m) + { + Int_t b = fM->GetXaxis()->FindBin(m); + if (b <= 0) return 0; + if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX(); + return static_cast(fBins->At(b)); + } + /** + * Add event observation + * + * @param triggered IF this is triggered + * @param m Multiplicity + * @param data Observation + */ + void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data) + { + fM->AddBinContent(1); + fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2)); + + MBin* all = static_cast(fBins->At(0)); + all->AddEvent(triggered, hasVtx, data); + + MBin* bin = FindBin(m); + bin->AddEvent(triggered, hasVtx, data); + } + /** + * End of job processing + * + * @param p Parent list + * @param o Parent output list + */ + void Finish(const TList* p, TList* o) + { + TList* l = dynamic_cast(p->FindObject(GetName())); + if (!l) { + Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName()); + return; + } + + TList* ours = new TList; + ours->SetName(GetName()); + ours->SetOwner(); + o->Add(ours); + + fM = static_cast(l->FindObject("m")); + if (!fM) { + Warning("Finish", "Didn't find object 'm' in %s", l->GetName()); + return; + } + TH1D* m = static_cast(fM->Clone("m")); + m->SetDirectory(0); + m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)"); + if (m->GetBinContent(1) > 0) + m->Scale(1. / m->GetBinContent(1)); + ours->Add(m); + + Int_t nBin = fM->GetNbinsX(); + TH1D* effs = static_cast(fM->Clone("effs")); + effs->SetYTitle("#epsilon_{X}"); + effs->SetFillColor(kRed+1); + effs->SetDirectory(0); + effs->SetMinimum(0); + + gStyle->SetPalette(1); + Int_t nCol = gStyle->GetNumberOfColors(); + THStack* stack = new THStack("biases", "Trigger biases"); + for (Int_t i = 0; i <= nBin; i++) { + MBin* bin = static_cast(fBins->At(i)); + effs->SetBinContent(i+1, bin->Finish(l, ours, stack)); + TH1* h = static_cast(stack->GetHists()->At(i)); + Int_t col = kBlack; + if (i != 0) { + Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5)); + col = gStyle->GetColorPalette(icol); + } + h->SetMarkerColor(col); + h->SetFillColor(col); + h->SetLineColor(col); + } + + ours->Add(stack); + ours->Add(effs); + } + UShort_t fMask; + TH1D* fM; + TObjArray* fBins; + ClassDef(TriggerType,1); + }; + TriggerType fInel; + TriggerType fInelGt0; + TriggerType fNSD; + AliFMDMCEventInspector fInspector; + TList* fList; + Bool_t fFirstEvent; + TH2D* fData; + TH1I* fTriggers; + UInt_t fTrackletRequirement; + UInt_t fVertexRequirement; + TAxis fVertexAxis; + TH1D* fVertexESD; + TH1D* fVertexMC; + TH1D* fM; + ClassDef(EvaluateTrigger,1); +}; +#else +//==================================================================== +void MakeEvaluateTriggers(const char* esddir, + Int_t nEvents = -1, + UInt_t vtx = 0x1, + UInt_t trk = 0x1, + UInt_t vz = 10, + Int_t proof = 0) +{ + // --- Libraries to load ------------------------------------------- + gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C"); + + // --- Check for proof mode, and possibly upload pars -------------- + if (proof> 0) { + gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C"); + if (!LoadPars(proof)) { + Error("MakeAOD", "Failed to load PARs"); + return; + } + } + + // --- Our data chain ---------------------------------------------- + gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C"); + TChain* chain = MakeChain("ESD", esddir,true); + // If 0 or less events is select, choose all + if (nEvents <= 0) nEvents = chain->GetEntries(); + + // --- Set the macro path ------------------------------------------ + gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:" + "$ALICE_ROOT/ANALYSIS/macros", + gROOT->GetMacroPath())); + + // --- Creating the manager and handlers --------------------------- + AliAnalysisManager *mgr = new AliAnalysisManager("Triggers", + "Forward multiplicity"); + AliAnalysisManager::SetCommonFileName("triggers.root"); + + // --- ESD input handler ------------------------------------------- + AliESDInputHandler *esdHandler = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdHandler); + + // --- Monte Carlo handler ----------------------------------------- + AliMCEventHandler* mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + mcHandler->SetReadTR(true); + + // --- Add tasks --------------------------------------------------- + // Physics selection + gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)")); + +#if 0 + // --- Fix up physics selection to give proper A,C, and E triggers - + AliInputEventHandler* ih = + static_cast(mgr->GetInputEventHandler()); + AliPhysicsSelection* ps = + static_cast(ih->GetEventSelection()); + // Ignore trigger class when selecting events. This mean that we + // get offline+(A,C,E) events too + ps->SetSkipTriggerClassSelection(true); +#endif + + // --- compile our code -------------------------------------------- + gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 " + "-I${ALICE_ROOT}/ANALYSIS " + "-I${ALICE_ROOT}/include -DBUILD=1"); + gROOT->LoadMacro("${ALICE_ROOT}/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C++g"); + + // --- Make our object --------------------------------------------- + EvaluateTrigger* task = new EvaluateTrigger("triggers"); + mgr->AddTask(task); + task->SetVertexRequirement(vtx); + task->SetTrackletRequirement(trk); + task->SetVertexAxis(10, -vz, vz); + + // --- create containers for input/output -------------------------- + AliAnalysisDataContainer *sums = + mgr->CreateContainer("triggerSums", TList::Class(), + AliAnalysisManager::kOutputContainer, + AliAnalysisManager::GetCommonFileName()); + AliAnalysisDataContainer *output = + mgr->CreateContainer("triggerResults", TList::Class(), + AliAnalysisManager::kParamContainer, + AliAnalysisManager::GetCommonFileName()); + + // --- connect input/output ---------------------------------------- + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, sums); + mgr->ConnectOutput(task, 2, output); + + // --- Run the analysis -------------------------------------------- + TStopwatch t; + if (!mgr->InitAnalysis()) { + Error("MakeAOD", "Failed to initialize analysis train!"); + return; + } + // Skip terminate if we're so requested and not in Proof or full mode + mgr->SetSkipTerminate(false); + // Some informative output + mgr->PrintStatus(); + if (proof) mgr->SetDebugLevel(3); + if (mgr->GetDebugLevel() < 1 && !proof) + mgr->SetUseProgressBar(kTRUE,100); + + // Run the train + t.Start(); + Printf("=== RUNNING ANALYSIS on %9d events ==========================", + nEvents); + mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents); + t.Stop(); + t.Print(); +} +#endif diff --git a/PWG2/FORWARD/analysis2/OtherData.C b/PWG2/FORWARD/analysis2/OtherData.C index 6bffe59d061..47164ffd5db 100644 --- a/PWG2/FORWARD/analysis2/OtherData.C +++ b/PWG2/FORWARD/analysis2/OtherData.C @@ -1150,9 +1150,11 @@ GetData(UShort_t sys, mp->Add(AliceCentralInelGt7000()); } } +#if 0 else Warning("GetData", "No other results for sys=%d, energy=%d", sys, energy); +#endif } else if (sys == 2) { // Nothing for PbPb so far @@ -1162,7 +1164,7 @@ GetData(UShort_t sys, en = Form(", #sqrt{s_{NN}}=%dGeV", energy); else en = Form(", #sqrt{s_{NN}}=%f4.2TeV", float(energy)/1000); - Warning("GetData", "No other data for PbP b yet"); + // Warning("GetData", "No other data for PbPb yet"); } else Warning("GetData", "Unknown system %d", sys); -- 2.39.3