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