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   fInspector.Init(fVtxAxis);
 326   fTrackDensity.DefineOutput(fList);
 327
 328   PostData(1, fList);
 329 }
 330 //____________________________________________________________________
 331 void
 332 AliCentralMCCorrectionsTask::UserExec(Option_t*)
 333 {
 334   //
 335   // Process each event
 336   //
 337   // Parameters:
 338   //    option Not used
 339   //
 340
 341   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event");
 342   // Get the input data - MC event
 343   AliMCEvent*  mcEvent = MCEvent();
 344   if (!mcEvent) {
 345     AliWarning("No MC event found");
 346     return;
 347   }
 348
 349   // Get the input data - ESD event
 350   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 351   if (!esd) {
 352     AliWarning("No ESD event found for input event");
 353     return;
 354   }
 355
 356   //--- Read run information -----------------------------------------
 357   if (fFirstEvent && esd->GetESDRun()) {
 358     fInspector.ReadRunDetails(esd);
 359
 360     AliInfo(Form("Initializing with parameters from the ESD:\n"
 361                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 362                  "         AliESDEvent::GetBeamType()     ->%s\n"
 363                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 364                  "         AliESDEvent::GetMagneticField()->%f\n"
 365                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 366                  esd->GetBeamEnergy(),
 367                  esd->GetBeamType(),
 368                  esd->GetCurrentL3(),
 369                  esd->GetMagneticField(),
 370                  esd->GetRunNumber()));
 371
 372     Print();
 373     fFirstEvent = false;
 374   }
 375
 376   // Some variables
 377   UInt_t   triggers; // Trigger bits
 378   Bool_t   lowFlux;  // Low flux flag
 379   UShort_t iVz;      // Vertex bin from ESD
 380   Double_t vZ;       // Z coordinate from ESD
 381   Double_t cent;     // Centrality
 382   UShort_t iVzMc;    // Vertex bin from MC
 383   Double_t vZMc;     // Z coordinate of IP vertex from MC
 384   Double_t b;        // Impact parameter
 385   Double_t cMC;      // Centrality estimate from b
 386   Int_t    nPart;    // Number of participants
 387   Int_t    nBin;     // Number of binary collisions
 388   Double_t phiR;     // Reaction plane from MC
 389   UShort_t nClusters;// Number of SPD clusters
 390   // Process the data
 391   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
 392                                      cent, nClusters);
 393   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
 394                        b, cMC, nPart, nBin, phiR);
 395
 396   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
 397   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
 398
 399   // Fill the event count histograms
 400   if (isInel)           fHEventsTr->Fill(vZMc);
 401   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
 402   fHEvents->Fill(vZMc);
 403
 404   // Now find our vertex bin object
 405   VtxBin* bin = 0;
 406   if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
 407     bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
 408   if (!bin) {
 409     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
 410     return;
 411   }
 412
 413   // Now process our input data and store in own ESD object
 414   fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
 415
 416   // Get the ESD object
 417   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 418
 419   // Count number of tracklets per bin
 420   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 421     bin->fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 422   //...and then the unused clusters in layer 1
 423   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 424      Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 425      bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j));
 426   }
 427
 428   // Count events
 429   bin->fCounts->Fill(0.5);
 430
 431   PostData(1, fList);
 432 }
 433
 434 //____________________________________________________________________
 435 void
 436 AliCentralMCCorrectionsTask::Terminate(Option_t*)
 437 {
 438   //
 439   // End of job
 440   //
 441   // Parameters:
 442   //    option Not used
 443   //
 444   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output");
 445
 446   fList = dynamic_cast<TList*>(GetOutputData(1));
 447   if (!fList) {
 448     AliError("No output list defined");
 449     return;
 450   }
 451
 452   DefineBins(fList);
 453
 454   // Output list
 455   TList* output = new TList;
 456   output->SetOwner();
 457   output->SetName(Form("%sResults", GetName()));
 458
 459   // --- Fill correction object --------------------------------------
 460   AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
 461   corr->SetVertexAxis(fVtxAxis);
 462   // corr->SetEtaAxis(fEtaAxis);
 463
 464  AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance;
 465   acorr->SetVertexAxis(fVtxAxis);
 466
 467   TIter     next(fVtxBins);
 468   VtxBin*   bin = 0;
 469   UShort_t  iVz = 1;
 470   while ((bin = static_cast<VtxBin*>(next())))
 471     bin->Finish(fList, output, iVz++, fEffectiveCorr, corr,acorr);
 472
 473   output->Add(corr);
 474  output->Add(acorr);
 475
 476   PostData(2, output);
 477 }
 478
 479 //____________________________________________________________________
 480 void
 481 AliCentralMCCorrectionsTask::Print(Option_t* option) const
 482 {
 483   std::cout << ClassName() << "\n"
 484             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
 485             << "  Vertex range:     [" << fVtxAxis.GetXmin()
 486             << "," << fVtxAxis.GetXmax() << "]\n"
 487             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
 488             << "  Eta range:        [" << fEtaAxis.GetXmin()
 489             << "," << fEtaAxis.GetXmax() << "]\n"
 490             << "  # of phi bins:    " << fNPhiBins
 491             << std::endl;
 492   gROOT->IncreaseDirLevel();
 493   fInspector.Print(option);
 494   fTrackDensity.Print(option);
 495   gROOT->DecreaseDirLevel();
 496 }
 497
 498 //====================================================================
 499 const char*
 500 AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
 501                                              Double_t high)
 502 {
 503   static TString buf;
 504   buf = Form("vtx%+05.1f_%+05.1f", low, high);
 505   buf.ReplaceAll("+", "p");
 506   buf.ReplaceAll("-", "m");
 507   buf.ReplaceAll(".", "d");
 508   return buf.Data();
 509 }
 510
 511
 512 //____________________________________________________________________
 513 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
 514   : fHits(0),
 515     fClusters(0),
 516     fPrimary(0),
 517     fCounts(0)
 518 {
 519 }
 520 //____________________________________________________________________
 521 AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
 522                                             Double_t     high,
 523                                             const TAxis& axis,
 524                                             UShort_t     nPhi)
 525   : TNamed(BinName(low, high),
 526            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
 527     fHits(0),
 528     fClusters(0),
 529     fPrimary(0),
 530     fCounts(0)
 531 {
 532   fPrimary = new TH2D("primary", "Primaries",
 533                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
 534                       nPhi, 0, 2*TMath::Pi());
 535   fPrimary->SetXTitle("#eta");
 536   fPrimary->SetYTitle("#varphi [radians]");
 537   fPrimary->Sumw2();
 538   fPrimary->SetDirectory(0);
 539
 540   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
 541   fHits->SetTitle("Hits");
 542   fHits->SetDirectory(0);
 543
 544   fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
 545   fClusters->SetTitle("Clusters");
 546   fClusters->SetDirectory(0);
 547
 548   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
 549   fCounts->SetXTitle("Events");
 550   fCounts->SetYTitle("# of Events");
 551   fCounts->SetDirectory(0);
 552 }
 553
 554 //____________________________________________________________________
 555 AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
 556   : TNamed(o),
 557     fHits(0),
 558     fClusters(0),
 559     fPrimary(0),
 560     fCounts(0)
 561 {
 562   if (o.fHits) {
 563     fHits = static_cast<TH2D*>(o.fHits->Clone());
 564     fHits->SetDirectory(0);
 565   }
 566   if (o.fClusters) {
 567     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
 568     fClusters->SetDirectory(0);
 569   }
 570   if (o.fPrimary) {
 571     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 572     fPrimary->SetDirectory(0);
 573   }
 574   if (o.fCounts) {
 575     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 576     fCounts->SetDirectory(0);
 577   }
 578 }
 579
 580 //____________________________________________________________________
 581 AliCentralMCCorrectionsTask::VtxBin&
 582 AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
 583 {
 584   if (&o == this) return *this;
 585   TNamed::operator=(o);
 586   fHits     = 0;
 587   fPrimary  = 0;
 588   fClusters = 0;
 589   fCounts   = 0;
 590   if (o.fHits) {
 591     fHits = static_cast<TH2D*>(o.fHits->Clone());
 592     fHits->SetDirectory(0);
 593   }
 594   if (o.fClusters) {
 595     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
 596     fClusters->SetDirectory(0);
 597   }
 598   if (o.fPrimary) {
 599     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 600     fPrimary->SetDirectory(0);
 601   }
 602   if (o.fCounts) {
 603     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 604     fCounts->SetDirectory(0);
 605   }
 606   return *this;
 607 }
 608
 609 //____________________________________________________________________
 610 void
 611 AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
 612 {
 613   TList* d = new TList;
 614   d->SetName(GetName());
 615   d->SetOwner();
 616   l->Add(d);
 617
 618   d->Add(fHits);
 619   d->Add(fClusters);
 620   d->Add(fPrimary);
 621   d->Add(fCounts);
 622 }
 623 //____________________________________________________________________
 624 void
 625 AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input,
 626                                             TList* output,
 627                                             UShort_t iVz,
 628                                             Bool_t effectiveCorr,
 629                                             AliCentralCorrSecondaryMap* map,
 630                                             AliCentralCorrAcceptance* acorr)
 631 {
 632   TList* out = new TList;
 633   out->SetName(GetName());
 634   out->SetOwner();
 635   output->Add(out);
 636
 637   TList* l = static_cast<TList*>(input->FindObject(GetName()));
 638   if (!l) {
 639     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
 640     return;
 641   }
 642
 643
 644   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
 645   TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
 646   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
 647   if (!hits || !prim) {
 648     AliError(Form("Missing histograms: %p, %p", hits, prim));
 649     return;
 650   }
 651
 652   TH2D* h = 0;
 653   if (effectiveCorr) h = static_cast<TH2D*>(clus->Clone("bgCorr"));
 654   else               h = static_cast<TH2D*>(hits->Clone("bgCorr"));
 655   h->SetDirectory(0);
 656   h->Divide(prim);
 657
 658   TH1D* acc = new TH1D(Form("SPDacc_vrtbin_%d",iVz),
 659                           "Acceptance correction for SPD" ,
 660                           fPrimary->GetXaxis()->GetNbins(),
 661                           fPrimary->GetXaxis()->GetXmin(),
 662                           fPrimary->GetXaxis()->GetXmax());
 663   TH1F* accden = static_cast<TH1F*>(acc->Clone(Form("%s_den",
 664                                                     acc->GetName())));
 665
 666   for(Int_t xx = 1; xx <=h->GetNbinsX(); xx++) {
 667     for(Int_t yy = 1; yy <=h->GetNbinsY(); yy++) {
 668       if(TMath::Abs(h->GetXaxis()->GetBinCenter(xx)) > 1.9) {
 669         h->SetBinContent(xx,yy,0.);
 670         h->SetBinError(xx,yy,0.);
 671       }
 672       if(h->GetBinContent(xx,yy) > 0.9)
 673         acc->Fill(h->GetXaxis()->GetBinCenter(xx));
 674       else {
 675         h->SetBinContent(xx,yy,0.);
 676         h->SetBinError(xx,yy,0.);
 677       }
 678       accden->Fill(h->GetXaxis()->GetBinCenter(xx));
 679
 680     }
 681   }
 682   acc->Divide(accden);
 683
 684   map->SetCorrection(iVz, h);
 685   acorr->SetCorrection(iVz, acc);
 686   out->Add(h);
 687   out->Add(acc);
 688 }
 689
 690 //
 691 // EOF
 692 //