1 //====================================================================
2 #include "AliForwarddNdetaTask.h"
9 #include <AliAnalysisManager.h>
10 #include <AliAODEvent.h>
11 #include <AliAODHandler.h>
12 #include <AliAODInputHandler.h>
13 #include "AliForwardUtil.h"
14 #include "AliAODForwardMult.h"
16 //____________________________________________________________________
17 AliForwarddNdetaTask::AliForwarddNdetaTask()
23 DGUARD(fDebug, 3, "Default CTOR of AliForwarddNdetaTask");
26 //____________________________________________________________________
27 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
28 : AliBasedNdetaTask("Forward")
36 DGUARD(fDebug, 3, "Named CTOR of AliForwarddNdetaTask");
39 //____________________________________________________________________
40 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
41 : AliBasedNdetaTask(o)
46 DGUARD(fDebug, 3, "Copy CTOR of AliForwarddNdetaTask");
49 //____________________________________________________________________
50 AliBasedNdetaTask::CentralityBin*
51 AliForwarddNdetaTask::MakeCentralityBin(const char* name, Short_t l,Short_t h)
55 // Make a new centrality bin
58 // name Histogram names
63 // Newly allocated object (of our type)
66 "Make a centrality bin for AliForwarddNdetaTask: %s [%d,%d]",
68 return new AliForwarddNdetaTask::CentralityBin(name, l, h);
72 //____________________________________________________________________
74 AliForwarddNdetaTask::GetHistogram(const AliAODEvent& aod, Bool_t mc)
77 // Retrieve the histogram
81 // mc Whether to get the MC histogram or not
84 // Retrieved histogram or null
86 // We should have a forward object at least
87 AliAODForwardMult* forward = GetForward(aod, mc, !mc);
88 if (!forward) return 0;
89 return &(forward->GetHistogram());
91 //____________________________________________________________________
93 AliForwarddNdetaTask::CheckEventData(Double_t vtx,
97 // Check if this is satellite
98 // if (!fSatelliteVertices) return;
99 Double_t aVtx = TMath::Abs(vtx);
100 if (aVtx < 37.5 || aVtx > 400) return;
102 TH2* hists[] = { data, dataMC };
104 // In satellite vertices FMD2i is cut away manually at this point
105 // for certain vertices. It could be done in the ESDs, but as of
106 // this writing not for specific vertices.
108 // cholm comment: It would be difficult to setup the filter in the
109 // reconstruction pass, but it could perhaps be done in the AOD
112 // This is what was done for
113 // the Pb-Pb paper (arXiv:1304.0347).
114 for (Int_t iX = 0; iX<=data->GetNbinsX(); iX++) {
115 // Do all checks up front - as soon as we can - branching is
117 Double_t x = data->GetXaxis()->GetBinCenter(iX);
119 if (((vtx > 60 && vtx < 90) && x < 3) ||
120 ((vtx > 330 && vtx < 350) && x > -2.5) ||
121 ((vtx < 100 || vtx > 305) && TMath::Abs(x) < 4.5) ||
122 (vtx < 50 && TMath::Abs(x) < 4.75))
126 for (Int_t iH = 0; iH < 2; iH++) {
127 if (!hists[iH]) continue;
128 // if (iX > hists[iH]->GetNbinsX()+1) continue;
129 // Also zero coverage and phi acceptance for this
130 for (Int_t iY = 0; iY<=hists[iH]->GetNbinsY()+1; iY++) {
131 hists[iH]->SetBinContent(iX, iY, 0);
132 hists[iH]->SetBinError(iX, iY, 0);
138 // Now, since we have some dead areas in FMD2i (sectors 16 and
139 // 17), we need to remove the corresponding bins from the
140 // histogram. However, it is not obvious which bins (in eta) to
141 // remove, so remove everything starting from the most negative to
142 // the middle of the histogram.
144 // This hack was first introduced by HHD, but was done at the end of
145 // the event processing (CentralityBin::MakeResults). That is,
146 // however, not very practical, as we'd like to normalize to the phi
147 // acceptance rather than the eta coverage and then correct for
148 // empty bins. Since the only way to really update the phi
149 // acceptance stored in the overflow bin is on the event level, we
150 // should really do it here.
151 const Int_t phiBin1 = 17; // Sector 16
152 const Int_t phiBin2 = 18; // Sector 17
153 for (Int_t iH = 0; iH < 2; iH++) {
154 if (!hists[iH]) continue;
156 Int_t midX = hists[iH]->GetNbinsX() / 2;
157 // Int_t nY = hists[iH]->GetNbinsY();
158 for (Int_t i = 1; i <= midX; i++) {
159 hists[iH]->SetBinContent(i, phiBin1, 0);
160 hists[iH]->SetBinContent(i, phiBin2, 0);
161 hists[iH]->SetBinError(i, phiBin1, 0);
162 hists[iH]->SetBinError(i, phiBin2, 0);
164 // Here, we should also modify the overflow bin to reflect the
165 // new phi acceptance. First get the old phi acceptance -
166 // then multiply this on the number of bins. This gives us -
167 // roughly - the number of sectors we had. Then take out two
168 // from that number, and then calculate the new phi
169 // Acceptance. Note, if the sectors where already taken out in
170 // the AOD production, we _will_ end up with a wrong number,
171 // so we should _not_ do that in the AOD production. This is
172 // tricky and may not work at all. For now, we should rely on
173 // the old way of correcting to the eta coverage and
174 // correcting for empty bins.
179 //========================================================================
181 AliForwarddNdetaTask::CentralityBin::End(TList* sums,
184 const TH2F* shapeCorr,
198 DGUARD(fDebug, 1,"In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d",
199 GetName(), corrEmpty, cutEdges, rootProj);
200 AliBasedNdetaTask::CentralityBin::End(sums, results, scheme,
201 shapeCorr, trigEff, trigEff0,
203 rootProj, corrEmpty, cutEdges,
204 triggerMask, marker, color, mclist,
207 if (!IsAllBin()) return;
208 TFile* file = TFile::Open("forward.root", "READ");
211 TList* forward = static_cast<TList*>(file->Get("ForwardSums"));
213 AliError("List Forward not found in forward.root");
216 TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
218 AliError("List ringResults not found in forward.root");
221 THStack* res = static_cast<THStack*>(rings->FindObject("all"));
223 AliError(Form("Stack all not found in %s", rings->GetName()));
227 AliError("Triggers histogram not set");
232 Double_t epsilonT = trigEff;
235 if (triggerMask == AliAODForwardMult::kNSD) {
236 // This is a local change
238 AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
241 AliInfo("Adding per-ring histograms to output");
243 Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
244 TIter next(res->GetHists());
246 while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
247 res->SetName("dndetaRings");
249 fOutput->Add(new TNamed("normCalc", text.Data()));
252 //________________________________________________________________________