PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx

   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"
  16 #include "AliCentralCorrectionManager.h"
  17 #include "AliCentralCorrAcceptance.h"
  18 #include "AliCentralCorrSecondaryMap.h"
  19 #include "AliAODForwardMult.h"
  20 #include "AliForwardUtil.h"
  21 #include "AliLog.h"
  22 #include "AliAODHandler.h"
  23 #include "AliAnalysisManager.h"
  24 #include "AliESDEvent.h"
  25 #include "AliMultiplicity.h"
  26 #include <TROOT.h>
  27 #include <TSystem.h>
  28 #include <TFile.h>
  29 #include <TError.h>
  30 #include <TSystem.h>
  31 #include <TObjArray.h>
  32 #include <iostream>
  33 #include <iomanip>
  34
  35 //====================================================================
  36 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
  37   : AliBaseESDTask(name, "AliCentralMultiplicityTask",
  38                    &(AliCentralCorrectionManager::Instance())),
  39     fInspector("event"),
  40     fAODCentral(kFALSE),
  41     fUseSecondary(true),
  42     fUseAcceptance(true),
  43   fIvz(0),
  44   fNClusterTracklet(0),
  45   fClusterPerTracklet(0),
  46   fNCluster(0),
  47   fNTracklet(0),
  48     fVtxList(0),
  49     fStore(false),
  50     fHData(0)
  51 {
  52   //
  53   // Constructor
  54   //
  55   DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
  56   fBranchNames =
  57     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  58     "SPDVertex.,PrimaryVertex.";
  59 }
  60 //____________________________________________________________________
  61 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
  62   : AliBaseESDTask(),
  63     fInspector(),
  64     fAODCentral(),
  65     fUseSecondary(true),
  66     fUseAcceptance(true),
  67   fIvz(0),
  68     fNClusterTracklet(0),
  69     fClusterPerTracklet(0),
  70     fNCluster(0),
  71     fNTracklet(0),
  72     fVtxList(0),
  73     fStore(false),
  74     fHData(0)
  75 {
  76   //
  77   // Constructor
  78   //
  79   DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
  80 }
  81
  82 //____________________________________________________________________
  83 void
  84 AliCentralMultiplicityTask::CreateBranches(AliAODHandler* ah)
  85 {
  86   //
  87   // Create output objects
  88   //
  89   //
  90   DGUARD(fDebug,1,"Create user output in AliCentralMultiplicityTask");
  91
  92   if (!ah)
  93     //AliFatal("No AOD output handler set in analysis manager");
  94     return;
  95
  96   TObject* obj = &fAODCentral;
  97   ah->AddBranch("AliAODCentralMult", &obj);
  98 }
  99 //____________________________________________________________________
 100 Bool_t
 101 AliCentralMultiplicityTask::Book()
 102 {
 103   fNeededCorrections = (AliCentralCorrectionManager::kSecondaryMap|
 104                         AliCentralCorrectionManager::kAcceptance);
 105   return true;
 106 }
 107
 108 //____________________________________________________________________
 109 Bool_t
 110 AliCentralMultiplicityTask::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}");
 122
 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);
 136   }
 137   fList->Add(hCoverage);
 138
 139   // Bins are
 140   TArrayD bins;
 141   // N-bins, loweset order, higest order, return array
 142   AliForwardUtil::MakeLogScale(300, 0, 5, bins);
 143   fNClusterTracklet = new TH2D("nClusterVsnTracklet",
 144                                "Total number of cluster vs number of tracklets",
 145                                bins.GetSize()-1, bins.GetArray(),
 146                                bins.GetSize()-1, bins.GetArray());
 147   fNClusterTracklet->SetDirectory(0);
 148   fNClusterTracklet->SetXTitle("N_{free cluster}");
 149   fNClusterTracklet->SetYTitle("N_{tracklet}");
 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
 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
 181   // Initialize the inspecto
 182   // fInspector.SetupForData(vertex);
 183
 184   fAODCentral.SetBit(AliAODCentralMult::kSecondary,  fUseSecondary);
 185   fAODCentral.SetBit(AliAODCentralMult::kAcceptance, fUseAcceptance);
 186
 187   return true;
 188 }
 189 //____________________________________________________________________
 190 Bool_t
 191 AliCentralMultiplicityTask::PreEvent()
 192 {
 193   fAODCentral.Clear("");
 194   return true;
 195 }
 196 //____________________________________________________________________
 197 Bool_t
 198 AliCentralMultiplicityTask::Event(AliESDEvent& esd)
 199 {
 200   //
 201   // Process each event
 202   //
 203   // Parameters:
 204   //    option Not used
 205   //
 206   DGUARD(fDebug,1,"Process event in AliCentralMultiplicityTask");
 207   fIvz               = 0;
 208   Bool_t   lowFlux   = kFALSE;
 209   UInt_t   triggers  = 0;
 210   UShort_t ivz       = 0;
 211   TVector3 ip;
 212   Double_t cent      = -1;
 213   UShort_t nClusters = 0;
 214   UInt_t   found     = fInspector.Process(&esd, triggers, lowFlux,
 215                                           ivz, ip, cent, nClusters);
 216
 217   // No event or no trigger
 218   if (found & AliFMDEventInspector::kNoEvent)    return false;
 219   if (found & AliFMDEventInspector::kNoTriggers) return false;
 220
 221   // Make sure AOD is filled
 222   MarkEventForStore();
 223
 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;
 228
 229   //Doing analysis
 230   const AliMultiplicity* spdmult = esd.GetMultiplicity();
 231
 232   TH2D& aodHist = fAODCentral.GetHistogram();
 233
 234   VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivz));
 235   if (!bin) return false;
 236
 237   ProcessESD(aodHist, spdmult);
 238   bin->Correct(aodHist, fUseSecondary, fUseAcceptance);
 239
 240   if (triggers & AliAODForwardMult::kInel)
 241     fHData->Add(&(fAODCentral.GetHistogram()));
 242
 243   return true;
 244 }
 245 //____________________________________________________________________
 246 void
 247 AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
 248                                        const AliMultiplicity* spdmult) const
 249 {
 250   DGUARD(fDebug,1,"Process the ESD in AliCentralMultiplicityTask");
 251   fNTracklet->Reset();
 252   fNCluster->Reset();
 253
 254   //Filling clusters in layer 1 used for tracklets...
 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   }
 260
 261   //...and then the unused ones in layer 1
 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
 275 }
 276
 277 //____________________________________________________________________
 278 Bool_t
 279 AliCentralMultiplicityTask::Finalize()
 280 {
 281   //
 282   // End of job
 283   //
 284   // Parameters:
 285   //    option Not used
 286   //
 287   DGUARD(fDebug,1,"Process merged output in AliCentralMultiplicityTask");
 288
 289   Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
 290   MakeSimpledNdeta(fList, fResults, nTr, nTrVtx, nAcc);
 291   AliInfoF("\n"
 292            "\t# events w/trigger:                 %f\n"
 293            "\t# events w/trigger+vertex:          %f\n"
 294            "\t# events accepted y cuts:           %f",
 295            nTr, nTrVtx, nAcc);
 296   return true;
 297 }
 298
 299 //____________________________________________________________________
 300 Bool_t
 301 AliCentralMultiplicityTask::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
 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
 400 //____________________________________________________________________
 401 void
 402 AliCentralMultiplicityTask::Print(Option_t* option) const
 403 {
 404   //
 405   // Print information
 406   //
 407   // Parameters:
 408   //    option Not used
 409   //
 410   AliBaseESDTask::Print(option);
 411   PFB("Use secondary correction", fUseSecondary);
 412   PFB("Use acceptance correction", fUseAcceptance);
 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();
 420       // fVtxList->ls();
 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       }
 430     }
 431   }
 432
 433   gROOT->IncreaseDirLevel();
 434   ccm.Print(option);
 435   // fInspector.Print(option);
 436   gROOT->DecreaseDirLevel();
 437
 438 }
 439
 440 //====================================================================
 441 AliCentralMultiplicityTask::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)
 452 {
 453 }
 454 //____________________________________________________________________
 455 AliCentralMultiplicityTask::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)
 465 {
 466 }
 467 //____________________________________________________________________
 468 AliCentralMultiplicityTask::VtxBin&
 469 AliCentralMultiplicityTask::VtxBin::operator=(const VtxBin& o)
 470 {
 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
 481   return *this;
 482 }
 483
 484 //____________________________________________________________________
 485 const char*
 486 AliCentralMultiplicityTask::VtxBin::GetName() const
 487 {
 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)));
 491 }
 492
 493 //____________________________________________________________________
 494 void
 495 AliCentralMultiplicityTask::VtxBin::SetupForData(TList* l,
 496                                                  TH2* coverage,
 497                                                  Bool_t store)
 498 {
 499   TList* out = 0;
 500   if (store) {
 501     out = new TList;
 502     out->SetName(GetName());
 503     out->SetOwner();
 504     l->Add(out);
 505   }
 506
 507   AliCentralCorrectionManager& ccm =
 508     AliCentralCorrectionManager::Instance();
 509
 510   // Clean-up
 511   if (fSec) {
 512     // delete fSec;
 513     fSec = 0;
 514   }
 515   if (fAcc) {
 516     // delete fAcc;
 517     fAcc = 0;
 518   }
 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     }
 539   }
 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);
 552   }
 553
 554   if (!store) {
 555     // If we're not asked to store anything, clean-up, and get out
 556     delete proj;
 557     return;
 558   }
 559
 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);
 609   }
 610
 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);
 616   }
 617 }
 618 //____________________________________________________________________
 619 void
 620 AliCentralMultiplicityTask::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
 633     Double_t eta    = aodHist.GetXaxis()->GetBinCenter(ix);
 634     Int_t    iax    = fAcc->GetXaxis()->FindBin(eta);
 635     Float_t  accCor = fAcc->GetBinContent(iax);
 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);
 650
 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       }
 657       if (!useAcceptance) continue;
 658
 659       // Acceptance correction
 660       Float_t accTmp = accCor;
 661       if (accTmp   < 0.000001) accTmp = 1;
 662       Float_t aodNew   = aodValue / accTmp ;
 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.);
 673       aodHist.SetBinContent(ix,nY+1, accCor);
 674     }
 675   } // for (ix)
 676   if (sum && fHits) fHits->Add(&aodHist);
 677 }
 678
 679 //____________________________________________________________________
 680 void
 681 AliCentralMultiplicityTask::VtxBin::Print(Option_t* /*option*/) const
 682 {
 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;
 689 }
 690
 691 //
 692 // EOF
 693 //