]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Spelling
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMultiplicityTask.cxx
CommitLineData
6f791cc3 1//====================================================================
2//
3// Base class for classes that calculate the multiplicity in the
4// central region event-by-event
5//
6// Inputs:
7// - AliESDEvent
8//
9// Outputs:
10// - AliAODCentralMult
11//
12// Histograms
13//
14// Corrections used
15#include "AliCentralMultiplicityTask.h"
8449e3e0 16#include "AliCentralCorrectionManager.h"
17#include "AliCentralCorrAcceptance.h"
18#include "AliCentralCorrSecondaryMap.h"
52047b6f 19#include "AliAODForwardMult.h"
6f791cc3 20#include "AliForwardUtil.h"
21#include "AliLog.h"
22#include "AliAODHandler.h"
6f791cc3 23#include "AliAnalysisManager.h"
24#include "AliESDEvent.h"
25#include "AliMultiplicity.h"
26#include <TROOT.h>
2a276c75 27#include <TSystem.h>
6f791cc3 28#include <TFile.h>
3e478dba 29#include <TError.h>
28b4012a 30#include <TSystem.h>
8449e3e0 31#include <TObjArray.h>
6f791cc3 32#include <iostream>
33#include <iomanip>
34
35//====================================================================
36AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
c8b1a7db 37 : AliBaseESDTask(name, "AliCentralMultiplicityTask",
38 &(AliCentralCorrectionManager::Instance())),
39 fInspector("event"),
6f791cc3 40 fAODCentral(kFALSE),
e58000b7 41 fUseSecondary(true),
9453b19e 42 fUseAcceptance(true),
8449e3e0 43 fIvz(0),
44 fNClusterTracklet(0),
c8b1a7db 45 fClusterPerTracklet(0),
46 fNCluster(0),
8449e3e0 47 fNTracklet(0),
48 fVtxList(0),
49 fStore(false),
c8b1a7db 50 fHData(0)
6f791cc3 51{
fb3430ac 52 //
53 // Constructor
54 //
5ca83fee 55 DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
a59cbd24 56 fBranchNames =
57 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
58 "SPDVertex.,PrimaryVertex.";
6f791cc3 59}
60//____________________________________________________________________
9c825779 61AliCentralMultiplicityTask::AliCentralMultiplicityTask()
c8b1a7db 62 : AliBaseESDTask(),
52047b6f 63 fInspector(),
9c825779 64 fAODCentral(),
e58000b7 65 fUseSecondary(true),
9453b19e 66 fUseAcceptance(true),
8449e3e0 67 fIvz(0),
28b4012a 68 fNClusterTracklet(0),
69 fClusterPerTracklet(0),
70 fNCluster(0),
9ecab72f 71 fNTracklet(0),
8449e3e0 72 fVtxList(0),
73 fStore(false),
c8b1a7db 74 fHData(0)
9c825779 75{
fb3430ac 76 //
77 // Constructor
78 //
5ca83fee 79 DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
9c825779 80}
2a276c75 81
9c825779 82//____________________________________________________________________
c8b1a7db 83void
84AliCentralMultiplicityTask::CreateBranches(AliAODHandler* ah)
3e478dba 85{
fb3430ac 86 //
87 // Create output objects
88 //
89 //
6ab100ec 90 DGUARD(fDebug,1,"Create user output in AliCentralMultiplicityTask");
6f791cc3 91
c8b1a7db 92 if (!ah)
8449e3e0 93 //AliFatal("No AOD output handler set in analysis manager");
c8b1a7db 94 return;
52047b6f 95
c8b1a7db 96 TObject* obj = &fAODCentral;
97 ah->AddBranch("AliAODCentralMult", &obj);
52047b6f 98}
52047b6f 99//____________________________________________________________________
c8b1a7db 100Bool_t
101AliCentralMultiplicityTask::Book()
52047b6f 102{
c8b1a7db 103 fNeededCorrections = (AliCentralCorrectionManager::kSecondaryMap|
104 AliCentralCorrectionManager::kAcceptance);
105 return true;
106}
5934a3e3 107
c8b1a7db 108//____________________________________________________________________
109Bool_t
110AliCentralMultiplicityTask::PreData(const TAxis& vertex, const TAxis& eta)
111{
112 // FindEtaLimits();
113 // unsigned short s = 1;
114 TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
115 eta.GetNbins(), eta.GetXmin(), eta.GetXmax(),
116 vertex.GetNbins(),vertex.GetXmin(),
117 vertex.GetXmax());
118 hCoverage->SetDirectory(0);
119 hCoverage->SetXTitle("#eta");
120 hCoverage->SetYTitle("v_{z} [cm]");
121 hCoverage->SetZTitle("n_{bins}");
6f791cc3 122
c8b1a7db 123 fAODCentral.Init(eta);
124
125 UShort_t nVz = vertex.GetNbins();
126 fVtxList = new TObjArray(nVz, 1);
127 fVtxList->SetName("centMultVtxBins");
128 fVtxList->SetOwner();
129
130 // Bool_t store = false;
131 for (Int_t v = 1; v <= nVz; v++) {
132 VtxBin* bin = new VtxBin(v, vertex.GetBinLowEdge(v),
133 vertex.GetBinUpEdge(v));
134 bin->SetupForData(fList, hCoverage, fStore);
135 fVtxList->AddAt(bin, v);
28b4012a 136 }
c8b1a7db 137 fList->Add(hCoverage);
9ecab72f 138
77f97e3f
CHC		 139 // Bins are
140 TArrayD bins;
141 // N-bins, loweset order, higest order, return array
142 AliForwardUtil::MakeLogScale(300, 0, 5, bins);
28b4012a 143 fNClusterTracklet = new TH2D("nClusterVsnTracklet",
9ecab72f 144 "Total number of cluster vs number of tracklets",
77f97e3f
CHC		 145 bins.GetSize()-1, bins.GetArray(),
146 bins.GetSize()-1, bins.GetArray());
28b4012a 147 fNClusterTracklet->SetDirectory(0);
77f97e3f
CHC		 148 fNClusterTracklet->SetXTitle("N_{free cluster}");
149 fNClusterTracklet->SetYTitle("N_{tracklet}");
28b4012a 150 fNClusterTracklet->SetStats(0);
151 fList->Add(fNClusterTracklet);
152
153 Int_t nEta = 80;
154 Double_t lEta = 2;
155 fClusterPerTracklet = new TH2D("clusterPerTracklet",
156 "N_{free cluster}/N_{tracklet} vs. #eta",
157 nEta,-lEta,lEta, 101, -.05, 10.05);
158 fClusterPerTracklet->SetDirectory(0);
159 fClusterPerTracklet->SetXTitle("#eta");
160 fClusterPerTracklet->SetYTitle("N_{free cluster}/N_{tracklet}");
161 fClusterPerTracklet->SetStats(0);
162 fList->Add(fClusterPerTracklet);
163
164 // Cache histograms
165 fNCluster = new TH1D("cacheCluster", "", nEta,-lEta,lEta);
166 fNCluster->SetDirectory(0);
167 fNCluster->Sumw2();
168
169 fNTracklet = new TH1D("cacheTracklet", "", nEta,-lEta,lEta);
170 fNTracklet->SetDirectory(0);
171 fNTracklet->Sumw2();
172
bfab35d9 173 fList->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondary));
174 fList->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
175
176 fHData = static_cast<TH2D*>(fAODCentral.GetHistogram().Clone("d2Ndetadphi"));
177 fHData->SetStats(0);
178 fHData->SetDirectory(0);
179 fList->Add(fHData);
180
28b4012a 181 // Initialize the inspecto
c8b1a7db 182 // fInspector.SetupForData(vertex);
28b4012a 183
bfab35d9 184 fAODCentral.SetBit(AliAODCentralMult::kSecondary, fUseSecondary);
185 fAODCentral.SetBit(AliAODCentralMult::kAcceptance, fUseAcceptance);
52047b6f 186
c8b1a7db 187 return true;
52047b6f 188}
9ecab72f 189//____________________________________________________________________
c8b1a7db 190Bool_t
191AliCentralMultiplicityTask::PreEvent()
9ecab72f 192{
c8b1a7db 193 fAODCentral.Clear("");
194 return true;
9ecab72f 195}
6f791cc3 196//____________________________________________________________________
c8b1a7db 197Bool_t
198AliCentralMultiplicityTask::Event(AliESDEvent& esd)
3e478dba 199{
fb3430ac 200 //
201 // Process each event
202 //
203 // Parameters:
204 // option Not used
205 //
6ab100ec 206 DGUARD(fDebug,1,"Process event in AliCentralMultiplicityTask");
8449e3e0 207 fIvz = 0;
52047b6f 208 Bool_t lowFlux = kFALSE;
209 UInt_t triggers = 0;
210 UShort_t ivz = 0;
5ca83fee 211 TVector3 ip;
52047b6f 212 Double_t cent = -1;
213 UShort_t nClusters = 0;
c8b1a7db 214 UInt_t found = fInspector.Process(&esd, triggers, lowFlux,
5ca83fee 215 ivz, ip, cent, nClusters);
52047b6f 216
217 // No event or no trigger
c8b1a7db 218 if (found & AliFMDEventInspector::kNoEvent) return false;
219 if (found & AliFMDEventInspector::kNoTriggers) return false;
6f791cc3 220
221 // Make sure AOD is filled
52047b6f 222 MarkEventForStore();
223
c8b1a7db 224 if (found & AliFMDEventInspector::kNoSPD) return false;
225 if (found & AliFMDEventInspector::kNoVertex) return false;
226 if (triggers & AliAODForwardMult::kPileUp) return false;
227 if (found & AliFMDEventInspector::kBadVertex) return false;
6f791cc3 228
229 //Doing analysis
c8b1a7db 230 const AliMultiplicity* spdmult = esd.GetMultiplicity();
52047b6f 231
232 TH2D& aodHist = fAODCentral.GetHistogram();
233
8449e3e0 234 VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivz));
c8b1a7db 235 if (!bin) return false;
1e0e4b43 236
237 ProcessESD(aodHist, spdmult);
8449e3e0 238 bin->Correct(aodHist, fUseSecondary, fUseAcceptance);
614f9452 239
bfab35d9 240 if (triggers & AliAODForwardMult::kInel)
241 fHData->Add(&(fAODCentral.GetHistogram()));
242
c8b1a7db 243 return true;
52047b6f 244}
245//____________________________________________________________________
246void
247AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
248 const AliMultiplicity* spdmult) const
249{
6ab100ec 250 DGUARD(fDebug,1,"Process the ESD in AliCentralMultiplicityTask");
28b4012a 251 fNTracklet->Reset();
252 fNCluster->Reset();
253
6f791cc3 254 //Filling clusters in layer 1 used for tracklets...
28b4012a 255 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
256 Double_t eta = spdmult->GetEta(j);
257 fNTracklet->Fill(eta);
258 aodHist.Fill(eta,spdmult->GetPhi(j));
259 }
3e478dba 260
6f791cc3 261 //...and then the unused ones in layer 1
28b4012a 262 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
263 Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
264 fNCluster->Fill(eta);
265 aodHist.Fill(eta, spdmult->GetPhiSingle(j));
266 }
267 fNClusterTracklet->Fill(fNCluster->GetEntries(),
268 fNTracklet->GetEntries());
269
270 fNCluster->Divide(fNTracklet);
271 for (Int_t j = 1; j <= fNCluster->GetNbinsX(); j++)
272 fClusterPerTracklet->Fill(fNCluster->GetXaxis()->GetBinCenter(j),
273 fNCluster->GetBinContent(j));
274
52047b6f 275}
276
6f791cc3 277//____________________________________________________________________
c8b1a7db 278Bool_t
279AliCentralMultiplicityTask::Finalize()
3e478dba 280{
fb3430ac 281 //
282 // End of job
283 //
284 // Parameters:
285 // option Not used
286 //
6ab100ec 287 DGUARD(fDebug,1,"Process merged output in AliCentralMultiplicityTask");
bfab35d9 288
bfab35d9 289 Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
c8b1a7db 290 MakeSimpledNdeta(fList, fResults, nTr, nTrVtx, nAcc);
bfab35d9 291 AliInfoF("\n"
292 "\t# events w/trigger: %f\n"
293 "\t# events w/trigger+vertex: %f\n"
3c69ddbf 294 "\t# events accepted by cuts: %f",
bfab35d9 295 nTr, nTrVtx, nAcc);
c8b1a7db 296 return true;
6f791cc3 297}
bfab35d9 298
299//____________________________________________________________________
300Bool_t
301AliCentralMultiplicityTask::MakeSimpledNdeta(const TList* input,
302 TList* output,
303 Double_t& nTr,
304 Double_t& nTrVtx,
305 Double_t& nAcc)
306{
307 // Get our histograms from the container
308 TH1I* hEventsTr = 0;
309 TH1I* hEventsTrVtx = 0;
310 TH1I* hEventsAcc = 0;
311 TH1I* hTriggers = 0;
312 if (!GetEventInspector().FetchHistograms(input,
313 hEventsTr,
314 hEventsTrVtx,
315 hEventsAcc,
316 hTriggers)) {
317 AliError(Form("Didn't get histograms from event selector "
318 "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
319 hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
320 input->ls();
321 return false;
322 }
323 nTr = hEventsTr->Integral();
324 nTrVtx = hEventsTrVtx->Integral();
325 nAcc = hEventsAcc->Integral();
326 Double_t vtxEff = nTrVtx / nTr;
327 TH2D* hData = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
328 if (!hData) {
329 AliError(Form("Couldn't get our summed histogram from output "
330 "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
331 input->ls();
332 return false;
333 }
334
335 Int_t nY = hData->GetNbinsY();
336 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, nY, "e");
337 TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1, nY, "e");
338 TH1D* norm = hData->ProjectionX("norm", 0, 0, "");
339 TH1D* phi = hData->ProjectionX("phi", nY+1, nY+1, "");
340 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
341 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
342 dNdeta->SetMarkerColor(kRed+1);
343 dNdeta->SetMarkerStyle(20);
344 dNdeta->SetDirectory(0);
345
346 dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
347 dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
348 dNdeta_->SetMarkerColor(kMagenta+1);
349 dNdeta_->SetMarkerStyle(21);
350 dNdeta_->SetDirectory(0);
351
352 norm->SetTitle("Normalization to #eta coverage");
353 norm->SetYTitle("#eta coverage");
354 norm->SetLineColor(kBlue+1);
355 norm->SetMarkerColor(kBlue+1);
356 norm->SetMarkerStyle(21);
357 norm->SetFillColor(kBlue+1);
358 norm->SetFillStyle(3005);
359 norm->SetDirectory(0);
360
361 phi->SetTitle("Normalization to #phi acceptance");
362 phi->SetYTitle("#phi acceptance");
363 phi->SetLineColor(kGreen+1);
364 phi->SetMarkerColor(kGreen+1);
365 phi->SetMarkerStyle(20);
366 phi->SetFillColor(kGreen+1);
367 phi->SetFillStyle(3004);
368 // phi->Scale(1. / nAcc);
369 phi->SetDirectory(0);
370
371 // dNdeta->Divide(norm);
372 dNdeta->Divide(phi);
373 dNdeta->SetStats(0);
374 dNdeta->Scale(vtxEff, "width");
375
376 dNdeta_->Divide(norm);
377 dNdeta_->SetStats(0);
378 dNdeta_->Scale(vtxEff, "width");
379
380 output->Add(dNdeta);
381 output->Add(dNdeta_);
382 output->Add(norm);
383 output->Add(phi);
384
385 return true;
386}
387
c8b1a7db 388#define PF(N,V,...) \
389 AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
390#define PFB(N,FLAG) \
391 do { \
392 AliForwardUtil::PrintName(N); \
393 std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
394 } while(false)
395#define PFV(N,VALUE) \
396 do { \
397 AliForwardUtil::PrintName(N); \
398 std::cout << (VALUE) << std::endl; } while(false)
399
6f791cc3 400//____________________________________________________________________
401void
52047b6f 402AliCentralMultiplicityTask::Print(Option_t* option) const
6f791cc3 403{
fb3430ac 404 //
405 // Print information
406 //
407 // Parameters:
408 // option Not used
409 //
c8b1a7db 410 AliBaseESDTask::Print(option);
411 PFB("Use secondary correction", fUseSecondary);
412 PFB("Use acceptance correction", fUseAcceptance);
8449e3e0 413
414 AliCentralCorrectionManager& ccm =
415 AliCentralCorrectionManager::Instance();
416 if (ccm.IsInit()) {
417 const AliCentralCorrSecondaryMap* secMap = ccm.GetSecondaryMap();
418 if (secMap) {
419 const TAxis& vaxis = secMap->GetVertexAxis();
c8b1a7db 420 // fVtxList->ls();
8449e3e0 421 std::cout << " Eta ranges:\n"
422 << " Vertex | Eta bins\n"
423 << " bin range | \n"
424 << " ----------------+-----------" << std::endl;
425 for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
426 VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(v));
427 if (!bin) continue;
428 bin->Print();
429 }
9ecab72f 430 }
431 }
432
52047b6f 433 gROOT->IncreaseDirLevel();
8449e3e0 434 ccm.Print(option);
c8b1a7db 435 // fInspector.Print(option);
52047b6f 436 gROOT->DecreaseDirLevel();
437
6f791cc3 438}
8449e3e0 439
3e478dba 440//====================================================================
8449e3e0 441AliCentralMultiplicityTask::VtxBin::VtxBin(Int_t iVz,
442 Double_t minIpZ,
443 Double_t maxIpZ)
444 : fId(iVz),
445 fMinIpZ(minIpZ),
446 fMaxIpZ(maxIpZ),
447 fEtaMin(999),
448 fEtaMax(0),
449 fSec(0),
450 fAcc(0),
451 fHits(0)
6f791cc3 452{
6f791cc3 453}
454//____________________________________________________________________
8449e3e0 455AliCentralMultiplicityTask::VtxBin::VtxBin(const VtxBin& o)
456 : TObject(o),
457 fId(o.fId),
458 fMinIpZ(o.fMinIpZ),
459 fMaxIpZ(o.fMaxIpZ),
460 fEtaMin(o.fEtaMin),
461 fEtaMax(o.fEtaMax),
462 fSec(o.fSec),
463 fAcc(o.fAcc),
464 fHits(o.fHits)
fb3430ac 465{
fb3430ac 466}
3e478dba 467//____________________________________________________________________
8449e3e0 468AliCentralMultiplicityTask::VtxBin&
469AliCentralMultiplicityTask::VtxBin::operator=(const VtxBin& o)
3e478dba 470{
8449e3e0 471 if (&o == this) return *this;
472 fId = o.fId;
473 fMinIpZ = o.fMinIpZ;
474 fMaxIpZ = o.fMaxIpZ;
475 fEtaMin = o.fEtaMin;
476 fEtaMax = o.fEtaMax;
477 fSec = o.fSec;
478 fAcc = o.fAcc;
479 fHits = o.fHits;
480
3e478dba 481 return *this;
482}
483
484//____________________________________________________________________
8449e3e0 485const char*
486AliCentralMultiplicityTask::VtxBin::GetName() const
3e478dba 487{
8449e3e0 488 return Form("%c%03d_%c%03d",
489 (fMinIpZ >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fMinIpZ)),
490 (fMaxIpZ >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fMaxIpZ)));
3e478dba 491}
492
493//____________________________________________________________________
8449e3e0 494void
495AliCentralMultiplicityTask::VtxBin::SetupForData(TList* l,
496 TH2* coverage,
497 Bool_t store)
3e478dba 498{
8449e3e0 499 TList* out = 0;
500 if (store) {
501 out = new TList;
502 out->SetName(GetName());
503 out->SetOwner();
504 l->Add(out);
6f791cc3 505 }
3e478dba 506
8449e3e0 507 AliCentralCorrectionManager& ccm =
508 AliCentralCorrectionManager::Instance();
3e478dba 509
8449e3e0 510 // Clean-up
511 if (fSec) {
512 // delete fSec;
513 fSec = 0;
6f791cc3 514 }
8449e3e0 515 if (fAcc) {
516 // delete fAcc;
517 fAcc = 0;
3e478dba 518 }
8449e3e0 519 // Get secondary correction and make a projection onto eta
520 TH2* sec = ccm.GetSecondaryMap()->GetCorrection(UShort_t(fId));
521 TH1* acc = ccm.GetAcceptance()->GetCorrection(UShort_t(fId));
522 fSec = static_cast<TH2*>(sec->Clone());
523 fAcc = static_cast<TH1*>(acc->Clone());
524 fSec->SetDirectory(0);
525 fAcc->SetDirectory(0);
526
527 TH1D* proj = fSec->ProjectionX("secondary");
528 proj->SetDirectory(0);
529 proj->Scale(1. / fSec->GetNbinsY());
530
531 // Find lower bound on eta
532 fEtaMin = proj->GetNbinsX();
533 for (Int_t e = 1; e <= proj->GetNbinsX(); e++) {
534 Double_t c = proj->GetBinContent(e);
535 if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
536 fEtaMin = e;
537 break;
538 }
28b4012a 539 }
8449e3e0 540 // Find upper bound on eta
541 fEtaMax = 1;
542 for (Int_t e = proj->GetNbinsX(); e >= 1; e--) {
543 Double_t c = proj->GetBinContent(e);
544 if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
545 fEtaMax = e;
546 break;
547 }
548 }
549 // Fill our coverage histogram
550 for (Int_t nn = fEtaMin; nn<=fEtaMax; nn++) {
551 coverage->SetBinContent(nn,fId,1);
28b4012a 552 }
553
8449e3e0 554 if (!store) {
555 // If we're not asked to store anything, clean-up, and get out
556 delete proj;
557 return;
28b4012a 558 }
559
8449e3e0 560 // Modify the title of the projection
561 proj->SetTitle(Form("Projection of secondary correction "
562 "for %+5.1f<v_{z}<%+5.1f",fMinIpZ, fMaxIpZ));
563 proj->SetYTitle("#LT 2^{nd} correction#GT");
564 proj->SetMarkerStyle(20);
565 proj->SetMarkerColor(kBlue+1);
566 out->Add(proj);
567
568 // Make some histograms to store diagnostics
569 TH2D* obg = static_cast<TH2D*>(fSec->Clone("secondaryMapFiducial"));
570 obg->SetTitle(Form("%s - fiducial volume", obg->GetTitle()));
571 obg->GetYaxis()->SetTitle("#varphi");
572 obg->SetDirectory(0);
573 out->Add(obg);
574
575 TH1D* after = static_cast<TH1D*>(proj->Clone("secondaryFiducial"));
576 after->SetDirectory(0);
577 after->GetYaxis()->SetTitle("#LT 2^{nd} correction#GT");
578 after->SetTitle(Form("%s - fiducial volume", after->GetTitle()));
579 after->SetMarkerColor(kRed+1);
580 out->Add(after);
581
582 if (fHits) {
583 // delete fHits;
584 fHits = 0;
585 }
586 fHits = static_cast<TH2D*>(fSec->Clone("hitMap"));
587 fHits->SetDirectory(0);
588 fHits->SetTitle(Form("d^{2}N/d#eta d#phi for %+5.1f<v_{z}<%+5.1f",
589 fMinIpZ, fMaxIpZ));
590 fHits->GetYaxis()->SetTitle("#varphi");
591 fHits->GetZaxis()->SetTitle("d^{2}N/d#eta d#varphi");
592 fHits->SetMarkerColor(kBlack);
593 fHits->SetMarkerStyle(1);
594 out->Add(fHits);
595
596 // Get the acceptance, and store that
597 TH1D* accClone = static_cast<TH1D*>(fAcc->Clone("acceptance"));
598 accClone->SetTitle(Form("Acceptance for %+5.1f<v_{z}<%+5.1f",
599 fMinIpZ, fMaxIpZ));
600 accClone->SetDirectory(0);
601 out->Add(accClone);
602
603 // Now zero content outside our eta range
604 for (Int_t e = 1; e < fEtaMin; e++) {
605 after->SetBinContent(e, 0);
606 after->SetBinError(e, 0);
607 for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
608 obg->SetBinContent(e,nn,0);
28b4012a 609 }
28b4012a 610
8449e3e0 611 for (Int_t e = fEtaMax+1; e <= proj->GetNbinsX(); e++) {
612 after->SetBinContent(e, 0);
613 after->SetBinError(e, 0);
614 for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
615 obg->SetBinContent(e,nn,0);
28b4012a 616 }
28b4012a 617}
8449e3e0 618//____________________________________________________________________
619void
620AliCentralMultiplicityTask::VtxBin::Correct(TH2D& aodHist,
621 Bool_t useSecondary,
622 Bool_t useAcceptance,
623 Bool_t sum) const
624{
625 if (useSecondary && fSec) aodHist.Divide(fSec);
626
627 Int_t nY = aodHist.GetNbinsY();
628 for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) {
629 Bool_t fiducial = true;
630 if (ix < fEtaMin || ix > fEtaMax) fiducial = false;
631 // Bool_t etabinSeen = kFALSE;
632
bfab35d9 633 Double_t eta = aodHist.GetXaxis()->GetBinCenter(ix);
634 Int_t iax = fAcc->GetXaxis()->FindBin(eta);
635 Float_t accCor = fAcc->GetBinContent(iax);
8449e3e0 636 // For test
637 // Float_t accErr = fAcc->GetBinError(ix);
638
639 // Loop over phi
640 for(Int_t iy = 1; iy <= nY; iy++) {
641 // If outside our fiducial volume, zero content
642 if (!fiducial) {
643 aodHist.SetBinContent(ix, iy, 0);
644 aodHist.SetBinError(ix, iy, 0);
645 continue;
646 }
647 // Get currrent value
648 Float_t aodValue = aodHist.GetBinContent(ix,iy);
649 Float_t aodErr = aodHist.GetBinError(ix,iy);
6f791cc3 650
8449e3e0 651 // Ignore very small values
652 if (aodValue < 0.000001) {
653 aodHist.SetBinContent(ix,iy, 0);
654 aodHist.SetBinError(ix,iy, 0);
655 continue;
656 }
1e0e4b43 657 if (!useAcceptance) continue;
8449e3e0 658
659 // Acceptance correction
1e0e4b43 660 Float_t accTmp = accCor;
661 if (accTmp < 0.000001) accTmp = 1;
662 Float_t aodNew = aodValue / accTmp ;
8449e3e0 663 aodHist.SetBinContent(ix,iy, aodNew);
664 aodHist.SetBinError(ix,iy,aodErr);
665 // - Test -
666 // Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
667 // TMath::Power(accErr/accCor,2) );
668 // test - aodHist.SetBinError(ix,iy,error);
669 } // for (iy)
670 //Filling underflow bin if we eta bin is in range
671 if (fiducial) {
672 aodHist.SetBinContent(ix,0, 1.);
1e0e4b43 673 aodHist.SetBinContent(ix,nY+1, accCor);
8449e3e0 674 }
675 } // for (ix)
676 if (sum && fHits) fHits->Add(&aodHist);
677}
678
52047b6f 679//____________________________________________________________________
8449e3e0 680void
681AliCentralMultiplicityTask::VtxBin::Print(Option_t* /*option*/) const
52047b6f 682{
8449e3e0 683 std::cout << " "
684 << std::setw(2) << fId << " "
685 << std::setw(5) << fMinIpZ << "-"
686 << std::setw(5) << fMaxIpZ << " | "
687 << std::setw(3) << fEtaMin << "-"
688 << std::setw(3) << fEtaMax << std::endl;
6f791cc3 689}
52047b6f 690
6f791cc3 691//
692// EOF
693//