PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx

   1 //
   2 // Calculate the corrections in the central regions
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODCentralMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliCentralMCCorrectionsTask.h"
  15 #include "AliTriggerAnalysis.h"
  16 #include "AliPhysicsSelection.h"
  17 #include "AliLog.h"
  18 #include "AliHeader.h"
  19 #include "AliGenEventHeader.h"
  20 #include "AliESDEvent.h"
  21 #include "AliAODHandler.h"
  22 #include "AliMultiplicity.h"
  23 #include "AliInputEventHandler.h"
  24 #include "AliStack.h"
  25 #include "AliMCEvent.h"
  26 #include "AliAODForwardMult.h"
  27 #include "AliCentralCorrSecondaryMap.h"
  28 #include "AliCentralCorrAcceptance.h"
  29 #include "AliForwardUtil.h"
  30 #include "AliMultiplicity.h"
  31 #include <TH1.h>
  32 #include <TH2D.h>
  33 #include <TDirectory.h>
  34 #include <TList.h>
  35 #include <TROOT.h>
  36 #include <iostream>
  37
  38 //====================================================================
  39 namespace {
  40   const char* GetEventName(Bool_t tr, Bool_t vtx)
  41   {
  42     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  43   }
  44 }
  45
  46 //====================================================================
  47 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
  48   : AliAnalysisTaskSE(),
  49     fInspector(),
  50     fTrackDensity(),
  51     fVtxBins(0),
  52     fFirstEvent(true),
  53     fHEvents(0),
  54     fHEventsTr(0),
  55     fHEventsTrVtx(0),
  56     fVtxAxis(),
  57     fEtaAxis(),
  58     fList(),
  59     fNPhiBins(20),
  60     fEffectiveCorr(true),
  61     fEtaCut(1.9),
  62     fCorrCut(0.6)
  63 {
  64   //
  65   // Constructor
  66   //
  67   // Parameters:
  68   //    name Name of task
  69   //
  70   DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
  71 }
  72
  73 //____________________________________________________________________
  74 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
  75   : AliAnalysisTaskSE(name),
  76     fInspector("eventInspector"),
  77     fTrackDensity("trackDensity"),
  78     fVtxBins(0),
  79     fFirstEvent(true),
  80     fHEvents(0),
  81     fHEventsTr(0),
  82     fHEventsTrVtx(0),
  83     fVtxAxis(10,-10,10),
  84     fEtaAxis(200,-4,6),
  85     fList(),
  86     fNPhiBins(20),
  87     fEffectiveCorr(true),
  88     fEtaCut(1.9),
  89     fCorrCut(0.6)
  90 {
  91   //
  92   // Constructor
  93   //
  94   // Parameters:
  95   //    name Name of task
  96   //
  97   DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
  98   DefineOutput(1, TList::Class());
  99   DefineOutput(2, TList::Class());
 100   fBranchNames =
 101     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
 102     "SPDVertex.,PrimaryVertex.";
 103 }
 104
 105 //____________________________________________________________________
 106 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o)
 107   : AliAnalysisTaskSE(o),
 108     fInspector(o.fInspector),
 109     fTrackDensity(),
 110     fVtxBins(0),
 111     fFirstEvent(o.fFirstEvent),
 112     fHEvents(o.fHEvents),
 113     fHEventsTr(o.fHEventsTr),
 114     fHEventsTrVtx(o.fHEventsTrVtx),
 115     fVtxAxis(10,-10,10),
 116     fEtaAxis(200,-4,6),
 117     fList(o.fList),
 118     fNPhiBins(o.fNPhiBins),
 119     fEffectiveCorr(o.fEffectiveCorr),
 120     fEtaCut(o.fEtaCut),
 121     fCorrCut(o.fCorrCut)
 122 {
 123   //
 124   // Copy constructor
 125   //
 126   // Parameters:
 127   //    o Object to copy from
 128   //
 129   DGUARD(fDebug, 3,"Copy CTOR of AliCentralMCCorrectionsTask");
 130 }
 131
 132 //____________________________________________________________________
 133 AliCentralMCCorrectionsTask&
 134 AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
 135 {
 136   //
 137   // Assignment operator
 138   //
 139   // Parameters:
 140   //    o Object to assign from
 141   //
 142   // Return:
 143   //    Reference to this object
 144   //
 145   DGUARD(fDebug,3,"Assignment of AliCentralMCCorrectionsTask");
 146   if (&o == this) return *this;
 147   fInspector         = o.fInspector;
 148   fTrackDensity      = o.fTrackDensity;
 149   fVtxBins           = o.fVtxBins;
 150   fFirstEvent        = o.fFirstEvent;
 151   fHEvents           = o.fHEvents;
 152   fHEventsTr         = o.fHEventsTr;
 153   fHEventsTrVtx      = o.fHEventsTrVtx;
 154   SetVertexAxis(o.fVtxAxis);
 155   SetEtaAxis(o.fEtaAxis);
 156   fNPhiBins          = o.fNPhiBins;
 157   fEffectiveCorr     = o.fEffectiveCorr;
 158   fEtaCut            = o.fEtaCut;
 159   fCorrCut           = o.fCorrCut;
 160
 161   return *this;
 162 }
 163
 164 //____________________________________________________________________
 165 void
 166 AliCentralMCCorrectionsTask::Init()
 167 {
 168   //
 169   // Initialize the task
 170   //
 171   //
 172   DGUARD(fDebug,1,"Initialize AliCentralMCCorrectionsTask");
 173 }
 174
 175 //____________________________________________________________________
 176 void
 177 AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
 178                                            Double_t max)
 179 {
 180   //
 181   // Set the vertex axis to use
 182   //
 183   // Parameters:
 184   //    nBins Number of bins
 185   //    vzMin Least @f$z@f$ coordinate of interation point
 186   //    vzMax Largest @f$z@f$ coordinate of interation point
 187   //
 188   DGUARD(fDebug,3,"Set vertex axis AliCentralMCCorrectionsTask [%d,%f,%f]",
 189          nBin, min, max);
 190   if (max < min) {
 191     Double_t tmp = min;
 192     min          = max;
 193     max          = tmp;
 194   }
 195   if (min < -10)
 196     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
 197   if (max > +10)
 198     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
 199   fVtxAxis.Set(nBin, min, max);
 200 }
 201 //____________________________________________________________________
 202 void
 203 AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
 204 {
 205   //
 206   // Set the vertex axis to use
 207   //
 208   // Parameters:
 209   //    axis Axis
 210   //
 211   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 212 }
 213
 214 //____________________________________________________________________
 215 void
 216 AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
 217 {
 218   //
 219   // Set the eta axis to use
 220   //
 221   // Parameters:
 222   //    nBins Number of bins
 223   //    vzMin Least @f$\eta@f$
 224   //    vzMax Largest @f$\eta@f$
 225   //
 226   DGUARD(fDebug,3,"Set eta axis AliCentralMCCorrectionsTask [%d,%f,%f]",
 227          nBin, min, max);
 228   if (max < min) {
 229     Double_t tmp = min;
 230     min          = max;
 231     max          = tmp;
 232   }
 233   if (min < -4)
 234     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
 235   if (max > +6)
 236     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
 237   fEtaAxis.Set(nBin, min, max);
 238 }
 239 //____________________________________________________________________
 240 void
 241 AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
 242 {
 243   //
 244   // Set the eta axis to use
 245   //
 246   // Parameters:
 247   //    axis Axis
 248   //
 249   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 250 }
 251
 252 //____________________________________________________________________
 253 void
 254 AliCentralMCCorrectionsTask::DefineBins(TList* l)
 255 {
 256   DGUARD(fDebug,1,"Define bins in AliCentralMCCorrectionsTask");
 257   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
 258   if (fVtxBins->GetEntries() > 0) return;
 259
 260   fVtxBins->SetOwner();
 261   for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
 262     Double_t low  = fVtxAxis.GetBinLowEdge(i);
 263     Double_t high = fVtxAxis.GetBinUpEdge(i);
 264     VtxBin*  bin  = new VtxBin(low, high, fEtaAxis, fNPhiBins);
 265     fVtxBins->AddAt(bin, i);
 266     bin->CreateOutputObjects(l);
 267   }
 268 }
 269
 270 //____________________________________________________________________
 271 void
 272 AliCentralMCCorrectionsTask::UserCreateOutputObjects()
 273 {
 274   //
 275   // Create output objects
 276   //
 277   //
 278   DGUARD(fDebug,1,"Create user output for AliCentralMCCorrectionsTask");
 279   fList = new TList;
 280   fList->SetOwner();
 281   fList->SetName(Form("%sSums", GetName()));
 282
 283   DefineBins(fList);
 284
 285   fHEvents = new TH1I(GetEventName(false,false),
 286                       "Number of all events",
 287                       fVtxAxis.GetNbins(),
 288                       fVtxAxis.GetXmin(),
 289                       fVtxAxis.GetXmax());
 290   fHEvents->SetXTitle("v_{z} [cm]");
 291   fHEvents->SetYTitle("# of events");
 292   fHEvents->SetFillColor(kBlue+1);
 293   fHEvents->SetFillStyle(3001);
 294   fHEvents->SetDirectory(0);
 295   fList->Add(fHEvents);
 296
 297   fHEventsTr = new TH1I(GetEventName(true, false),
 298                         "Number of triggered events",
 299                         fVtxAxis.GetNbins(),
 300                         fVtxAxis.GetXmin(),
 301                         fVtxAxis.GetXmax());
 302   fHEventsTr->SetXTitle("v_{z} [cm]");
 303   fHEventsTr->SetYTitle("# of events");
 304   fHEventsTr->SetFillColor(kRed+1);
 305   fHEventsTr->SetFillStyle(3001);
 306   fHEventsTr->SetDirectory(0);
 307   fList->Add(fHEventsTr);
 308
 309   fHEventsTrVtx = new TH1I(GetEventName(true,true),
 310                            "Number of events w/trigger and vertex",
 311                            fVtxAxis.GetNbins(),
 312                            fVtxAxis.GetXmin(),
 313                            fVtxAxis.GetXmax());
 314   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 315   fHEventsTrVtx->SetYTitle("# of events");
 316   fHEventsTrVtx->SetFillColor(kBlue+1);
 317   fHEventsTrVtx->SetFillStyle(3001);
 318   fHEventsTrVtx->SetDirectory(0);
 319   fList->Add(fHEventsTrVtx);
 320
 321   // Copy axis objects to output
 322   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
 323                           fVtxAxis.GetNbins(),
 324                           fVtxAxis.GetXmin(),
 325                           fVtxAxis.GetXmax());
 326   TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
 327                           fEtaAxis.GetNbins(),
 328                           fEtaAxis.GetXmin(),
 329                           fEtaAxis.GetXmax());
 330   fList->Add(vtxAxis);
 331   fList->Add(etaAxis);
 332
 333   fList->Add(AliForwardUtil::MakeParameter("effective", fEffectiveCorr));
 334   fList->Add(AliForwardUtil::MakeParameter("maxEta", fEtaCut));
 335   fList->Add(AliForwardUtil::MakeParameter("accCut", fCorrCut));
 336
 337   AliInfo(Form("Initialising sub-routines: %p, %p",
 338                &fInspector, &fTrackDensity));
 339   fInspector.CreateOutputObjects(fList);
 340   fTrackDensity.CreateOutputObjects(fList);
 341
 342   PostData(1, fList);
 343 }
 344 //____________________________________________________________________
 345 void
 346 AliCentralMCCorrectionsTask::UserExec(Option_t*)
 347 {
 348   //
 349   // Process each event
 350   //
 351   // Parameters:
 352   //    option Not used
 353   //
 354
 355   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event");
 356   // Get the input data - MC event
 357   AliMCEvent*  mcEvent = MCEvent();
 358   if (!mcEvent) {
 359     AliWarning("No MC event found");
 360     return;
 361   }
 362
 363   // Get the input data - ESD event
 364   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 365   if (!esd) {
 366     AliWarning("No ESD event found for input event");
 367     return;
 368   }
 369
 370   //--- Read run information -----------------------------------------
 371   if (fFirstEvent && esd->GetESDRun()) {
 372     fInspector.ReadRunDetails(esd);
 373
 374     AliInfo(Form("Initializing with parameters from the ESD:\n"
 375                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 376                  "         AliESDEvent::GetBeamType()     ->%s\n"
 377                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 378                  "         AliESDEvent::GetMagneticField()->%f\n"
 379                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 380                  esd->GetBeamEnergy(),
 381                  esd->GetBeamType(),
 382                  esd->GetCurrentL3(),
 383                  esd->GetMagneticField(),
 384                  esd->GetRunNumber()));
 385
 386     fInspector.SetupForData(fVtxAxis);
 387
 388     Print();
 389     fFirstEvent = false;
 390   }
 391
 392   // Some variables
 393   UInt_t   triggers  = 0;    // Trigger bits
 394   Bool_t   lowFlux   = true; // Low flux flag
 395   UShort_t iVz       = 0;    // Vertex bin from ESD
 396   TVector3 ip;               // Z coordinate from ESD
 397   Double_t cent      = -1;   // Centrality
 398   UShort_t iVzMc     = 0;    // Vertex bin from MC
 399   Double_t vZMc      = 1000; // Z coordinate of IP vertex from MC
 400   Double_t b         = -1;   // Impact parameter
 401   Double_t cMC       = -1;   // Centrality estimate from b
 402   Int_t    nPart     = -1;   // Number of participants
 403   Int_t    nBin      = -1;   // Number of binary collisions
 404   Double_t phiR      = 1000; // Reaction plane from MC
 405   UShort_t nClusters = 0;    // Number of SPD clusters
 406   // Process the data
 407   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip,
 408                                      cent, nClusters);
 409   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
 410                        b, cMC, nPart, nBin, phiR);
 411
 412   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
 413   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
 414
 415   // Fill the event count histograms
 416   if (isInel)           fHEventsTr->Fill(vZMc);
 417   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
 418   fHEvents->Fill(vZMc);
 419
 420   // Now find our vertex bin object
 421   VtxBin* bin = 0;
 422   if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
 423     bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
 424   if (!bin) {
 425     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
 426     return;
 427   }
 428
 429   // Now process our input data and store in own ESD object
 430   fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
 431
 432   // Get the ESD object
 433   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 434
 435   // Count number of tracklets per bin
 436   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 437     bin->fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 438   //...and then the unused clusters in layer 1
 439   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 440      Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 441      bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j));
 442   }
 443
 444   // Count events
 445   bin->fCounts->Fill(0.5);
 446
 447   PostData(1, fList);
 448 }
 449
 450 //____________________________________________________________________
 451 void
 452 AliCentralMCCorrectionsTask::Terminate(Option_t*)
 453 {
 454   //
 455   // End of job
 456   //
 457   // Parameters:
 458   //    option Not used
 459   //
 460   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output");
 461
 462   fList = dynamic_cast<TList*>(GetOutputData(1));
 463   if (!fList) {
 464     AliError("No output list defined");
 465     return;
 466   }
 467
 468   DefineBins(fList);
 469
 470   // Output list
 471   TList* output = new TList;
 472   output->SetOwner();
 473   output->SetName(Form("%sResults", GetName()));
 474
 475   // --- Fill correction object --------------------------------------
 476   AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
 477   corr->SetVertexAxis(fVtxAxis);
 478   // corr->SetEtaAxis(fEtaAxis);
 479
 480   AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance;
 481   acorr->SetVertexAxis(fVtxAxis);
 482
 483   TIter     next(fVtxBins);
 484   VtxBin*   bin = 0;
 485   UShort_t  iVz = 1;
 486   while ((bin = static_cast<VtxBin*>(next())))
 487     bin->Terminate(fList, output, iVz++, fEffectiveCorr,
 488                    fEtaCut, fCorrCut, corr,acorr);
 489
 490   output->Add(corr);
 491   output->Add(acorr);
 492
 493   PostData(2, output);
 494 }
 495
 496 //____________________________________________________________________
 497 void
 498 AliCentralMCCorrectionsTask::Print(Option_t* option) const
 499 {
 500   std::cout << ClassName() << "\n"
 501             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
 502             << "  Vertex range:     [" << fVtxAxis.GetXmin()
 503             << "," << fVtxAxis.GetXmax() << "]\n"
 504             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
 505             << "  Eta range:        [" << fEtaAxis.GetXmin()
 506             << "," << fEtaAxis.GetXmax() << "]\n"
 507             << "  # of phi bins:    " << fNPhiBins << "\n"
 508             << "  Effective corr.:  " << fEffectiveCorr << "\n"
 509             << "  Eta cut-off:      " << fEtaCut << "\n"
 510             << "  Acceptance cut:   " << fCorrCut
 511             << std::endl;
 512   gROOT->IncreaseDirLevel();
 513   fInspector.Print(option);
 514   fTrackDensity.Print(option);
 515   gROOT->DecreaseDirLevel();
 516 }
 517
 518 //====================================================================
 519 const char*
 520 AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
 521                                              Double_t high)
 522 {
 523   static TString buf;
 524   buf = Form("vtx%+05.1f_%+05.1f", low, high);
 525   buf.ReplaceAll("+", "p");
 526   buf.ReplaceAll("-", "m");
 527   buf.ReplaceAll(".", "d");
 528   return buf.Data();
 529 }
 530
 531
 532 //____________________________________________________________________
 533 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
 534   : fHits(0),
 535     fClusters(0),
 536     fPrimary(0),
 537     fCounts(0)
 538 {
 539 }
 540 //____________________________________________________________________
 541 AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
 542                                             Double_t     high,
 543                                             const TAxis& axis,
 544                                             UShort_t     nPhi)
 545   : TNamed(BinName(low, high),
 546            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
 547     fHits(0),
 548     fClusters(0),
 549     fPrimary(0),
 550     fCounts(0)
 551 {
 552   fPrimary = new TH2D("primary", "Primaries",
 553                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
 554                       nPhi, 0, 2*TMath::Pi());
 555   fPrimary->SetXTitle("#eta");
 556   fPrimary->SetYTitle("#varphi [radians]");
 557   fPrimary->Sumw2();
 558   fPrimary->SetDirectory(0);
 559
 560   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
 561   fHits->SetTitle("Hits");
 562   fHits->SetDirectory(0);
 563
 564   fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
 565   fClusters->SetTitle("Clusters");
 566   fClusters->SetDirectory(0);
 567
 568   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
 569   fCounts->SetXTitle("Events");
 570   fCounts->SetYTitle("# of Events");
 571   fCounts->SetDirectory(0);
 572 }
 573
 574 //____________________________________________________________________
 575 AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
 576   : TNamed(o),
 577     fHits(0),
 578     fClusters(0),
 579     fPrimary(0),
 580     fCounts(0)
 581 {
 582   if (o.fHits) {
 583     fHits = static_cast<TH2D*>(o.fHits->Clone());
 584     fHits->SetDirectory(0);
 585   }
 586   if (o.fClusters) {
 587     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
 588     fClusters->SetDirectory(0);
 589   }
 590   if (o.fPrimary) {
 591     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 592     fPrimary->SetDirectory(0);
 593   }
 594   if (o.fCounts) {
 595     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 596     fCounts->SetDirectory(0);
 597   }
 598 }
 599
 600 //____________________________________________________________________
 601 AliCentralMCCorrectionsTask::VtxBin&
 602 AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
 603 {
 604   if (&o == this) return *this;
 605   TNamed::operator=(o);
 606   fHits     = 0;
 607   fPrimary  = 0;
 608   fClusters = 0;
 609   fCounts   = 0;
 610   if (o.fHits) {
 611     fHits = static_cast<TH2D*>(o.fHits->Clone());
 612     fHits->SetDirectory(0);
 613   }
 614   if (o.fClusters) {
 615     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
 616     fClusters->SetDirectory(0);
 617   }
 618   if (o.fPrimary) {
 619     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 620     fPrimary->SetDirectory(0);
 621   }
 622   if (o.fCounts) {
 623     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 624     fCounts->SetDirectory(0);
 625   }
 626   return *this;
 627 }
 628
 629 //____________________________________________________________________
 630 void
 631 AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
 632 {
 633   TList* d = new TList;
 634   d->SetName(GetName());
 635   d->SetOwner();
 636   l->Add(d);
 637
 638   d->Add(fHits);
 639   d->Add(fClusters);
 640   d->Add(fPrimary);
 641   d->Add(fCounts);
 642 }
 643 //____________________________________________________________________
 644 void
 645 AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input,
 646                                                TList* output,
 647                                                UShort_t iVz,
 648                                                Bool_t effectiveCorr,
 649                                                Double_t etaCut,
 650                                                Double_t accCut,
 651                                                AliCentralCorrSecondaryMap* map,
 652                                                AliCentralCorrAcceptance* acorr)
 653 {
 654   TList* out = new TList;
 655   out->SetName(GetName());
 656   out->SetOwner();
 657   output->Add(out);
 658
 659   TList* l = static_cast<TList*>(input->FindObject(GetName()));
 660   if (!l) {
 661     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
 662     return;
 663   }
 664
 665   // Get the sums
 666   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
 667   TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
 668   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
 669   if (!hits || !prim) {
 670     AliError(Form("Missing histograms: %p, %p", hits, prim));
 671     return;
 672   }
 673
 674   // Clone cluster and hit map
 675   TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
 676   TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
 677   secMapEff->SetTitle("2^{nd} map from clusters");
 678   secMapEff->SetDirectory(0);
 679   secMapHit->SetTitle("2^{nd} map from MC hits");
 680   secMapHit->SetDirectory(0);
 681
 682   // Divide cluster and hit map with number of primaries
 683   secMapEff->Divide(prim);
 684   secMapHit->Divide(prim);
 685
 686   // Create acceptance histograms
 687   TH1D* accEff = new TH1D("accEff",
 688                           "Acceptance correction for SPD (from 2^{nd} map)" ,
 689                           fPrimary->GetXaxis()->GetNbins(),
 690                           fPrimary->GetXaxis()->GetXmin(),
 691                           fPrimary->GetXaxis()->GetXmax());
 692   TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
 693   accHit->SetTitle("Acceptance correction for SPD (from clusters)");
 694
 695   // Diagnostics histogra,
 696   TH2*  dia    = static_cast<TH2D*>(clus->Clone("diagnostics"));
 697   dia->SetTitle("Scaled cluster density");
 698
 699   // Total number of channels along phi and # of eta bins
 700   Int_t nTotal = secMapHit->GetNbinsY();
 701   Int_t nEta   = secMapHit->GetNbinsX();
 702
 703   for(Int_t xx = 1; xx <= nEta; xx++) {
 704     Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
 705     Bool_t   ins = TMath::Abs(eta) <= etaCut;
 706     Double_t mm  = 0;
 707     if (ins) {
 708       // Find the maximum cluster signal in this phi range
 709       for (Int_t yy = 1; yy <= nTotal; yy++) {
 710         Double_t c = clus->GetBinContent(xx,yy);
 711         mm         = TMath::Max(mm, c);
 712       }
 713     }
 714     // Count number of phi bins with enough clusters or high enough
 715     // correction.
 716     Int_t nOKEff    = 0;
 717     Int_t nOKHit    = 0;
 718     for(Int_t yy = 1; yy <=nTotal; yy++) {
 719       if (!ins) { // Not inside Eta cut
 720         secMapEff->SetBinContent(xx,yy,0.);
 721         secMapEff->SetBinError(xx,yy,0.);
 722         secMapHit->SetBinContent(xx,yy,0.);
 723         secMapHit->SetBinError(xx,yy,0.);
 724         dia->SetBinContent(xx,yy,0);
 725         continue;
 726       }
 727
 728       // Check if the background correction is large enough, or zero map
 729       if(secMapEff->GetBinContent(xx,yy) > 0.9) {
 730         // acc->Fill(h->GetXaxis()->GetBinCenter(xx));
 731         nOKEff++;
 732       }
 733       else {
 734         secMapEff->SetBinContent(xx,yy,0.);
 735         secMapEff->SetBinError(xx,yy,0.);
 736       }
 737
 738       // Check if the number of cluster is large enough, or zero map
 739       Double_t c = clus->GetBinContent(xx,yy);
 740       Double_t s = (mm < 1e-6) ? 0 : c / mm;
 741       dia->SetBinContent(xx,yy,s);
 742       if (s >= accCut) {
 743         nOKHit++;
 744       }
 745       else {
 746         secMapHit->SetBinContent(xx,yy,0);
 747         secMapHit->SetBinError(xx,yy,0);
 748       }
 749     }
 750
 751     // Calculate acceptance as ratio of bins with enough clusters and
 752     // total number of phi bins.
 753     Double_t accXX = float(nOKHit) / nTotal;
 754     if (accXX < 0.2) accXX = 0;
 755     accHit->SetBinContent(xx, accXX);
 756
 757     // Calculate acceptance as ratio of bins with large enough
 758     // correction and total number of phi bins.
 759     accXX = float(nOKEff) / nTotal;
 760     if (accXX < 0.2) accXX = 0;
 761     accEff->SetBinContent(xx, accXX);
 762   }
 763
 764   TH2D* secMap    = (effectiveCorr ? secMapEff : secMapHit);
 765   TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
 766   TH1D* acc       = (effectiveCorr ? accEff    : accHit);
 767   TH1D* accAlt    = (effectiveCorr ? accHit    : accEff);
 768   out->Add(secMap->Clone("secMap"));
 769   out->Add(secMapAlt->Clone());
 770   out->Add(acc->Clone("acc"));
 771   out->Add(accAlt->Clone());
 772   out->Add(dia->Clone());
 773
 774   map->SetCorrection(iVz, secMap);
 775   acorr->SetCorrection(iVz, acc);
 776 }
 777
 778 //
 779 // EOF
 780 //