PWGLF/FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx

   1 //
   2 // Calculate the corrections in the forward regions
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODForwardMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliForwardMCCorrectionsTask.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 "AliFMDStripIndex.h"
  28 #include "AliFMDCorrSecondaryMap.h"
  29 #include <TH1.h>
  30 #include <TH2D.h>
  31 #include <TDirectory.h>
  32 #include <TTree.h>
  33 #include <TList.h>
  34 #include <TROOT.h>
  35 #include <iostream>
  36
  37 //====================================================================
  38 namespace {
  39   const char* GetEventName(Bool_t tr, Bool_t vtx)
  40   {
  41     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  42   }
  43 #if 0
  44   const char* GetHitsName(UShort_t d, Char_t r)
  45   {
  46     return Form("hitsFMD%d%c", d, r);
  47   }
  48   const char* GetStripsName(UShort_t d, Char_t r)
  49   {
  50     return Form("stripsFMD%d%c", d, r);
  51   }
  52   const char* GetPrimaryName(Char_t r, Bool_t trVtx)
  53   {
  54     return Form("primaries%s%s",
  55                 ((r == 'I' || r == 'i') ? "Inner" : "Outer"),
  56                 (trVtx ? "TrVtx" : "All"));
  57   }
  58 #endif
  59 }
  60
  61 //====================================================================
  62 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask()
  63   : AliAnalysisTaskSE(),
  64     fInspector(),
  65     fTrackDensity(),
  66     fESDFMD(),
  67     fVtxBins(0),
  68     fFirstEvent(true),
  69     fHEvents(0),
  70     fHEventsTr(0),
  71     fHEventsTrVtx(0),
  72     fVtxAxis(),
  73     fEtaAxis(),
  74     fList(),
  75     fUseESDVertexCoordinate(false),
  76     fCalculateafterESDeventcuts(false)
  77 {
  78   //
  79   // Constructor
  80   //
  81   // Parameters:
  82   //    name Name of task
  83   //
  84 }
  85
  86 //____________________________________________________________________
  87 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name)
  88   : AliAnalysisTaskSE(name),
  89     fInspector("eventInspector"),
  90     fTrackDensity("trackDensity"),
  91     fESDFMD(),
  92     fVtxBins(0),
  93     fFirstEvent(true),
  94     fHEvents(0),
  95     fHEventsTr(0),
  96     fHEventsTrVtx(0),
  97     fVtxAxis(10,-10,10),
  98     fEtaAxis(200,-4,6),
  99     fList(),
 100     fUseESDVertexCoordinate(false),
 101     fCalculateafterESDeventcuts(false)
 102
 103 {
 104   //
 105   // Constructor
 106   //
 107   // Parameters:
 108   //    name Name of task
 109   //
 110   DefineOutput(1, TList::Class());
 111   DefineOutput(2, TList::Class());
 112   fBranchNames =
 113     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
 114     "AliESDFMD.,SPDVertex.,PrimaryVertex.";
 115 }
 116
 117 //____________________________________________________________________
 118 AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o)
 119   : AliAnalysisTaskSE(o),
 120     fInspector(o.fInspector),
 121     fTrackDensity(),
 122     fESDFMD(o.fESDFMD),
 123     fVtxBins(0),
 124     fFirstEvent(o.fFirstEvent),
 125     fHEvents(o.fHEvents),
 126     fHEventsTr(o.fHEventsTr),
 127     fHEventsTrVtx(o.fHEventsTrVtx),
 128     fVtxAxis(10,-10,10),
 129     fEtaAxis(200,-4,6),
 130     fList(o.fList),
 131     fUseESDVertexCoordinate(o.fUseESDVertexCoordinate),
 132     fCalculateafterESDeventcuts(o.fCalculateafterESDeventcuts)
 133
 134 {
 135   //
 136   // Copy constructor
 137   //
 138   // Parameters:
 139   //    o Object to copy from
 140   //
 141 }
 142
 143 //____________________________________________________________________
 144 AliForwardMCCorrectionsTask&
 145 AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o)
 146 {
 147   //
 148   // Assignment operator
 149   //
 150   // Parameters:
 151   //    o Object to assign from
 152   //
 153   // Return:
 154   //    Reference to this object
 155   //
 156   if (&o == this) return *this;
 157   fInspector         = o.fInspector;
 158   fTrackDensity      = o.fTrackDensity;
 159   fESDFMD            = o.fESDFMD;
 160   fVtxBins           = o.fVtxBins;
 161   fFirstEvent        = o.fFirstEvent;
 162   fHEvents           = o.fHEvents;
 163   fHEventsTr         = o.fHEventsTr;
 164   fHEventsTrVtx      = o.fHEventsTrVtx;
 165   SetVertexAxis(o.fVtxAxis);
 166   SetEtaAxis(o.fEtaAxis);
 167   fUseESDVertexCoordinate=o.fUseESDVertexCoordinate;
 168   fCalculateafterESDeventcuts=o.fCalculateafterESDeventcuts;
 169
 170   return *this;
 171 }
 172
 173 //____________________________________________________________________
 174 void
 175 AliForwardMCCorrectionsTask::Init()
 176 {
 177   //
 178   // Initialize the task
 179   //
 180   //
 181 }
 182
 183 //____________________________________________________________________
 184 void
 185 AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
 186                                            Double_t max)
 187 {
 188   //
 189   // Set the vertex axis to use
 190   //
 191   // Parameters:
 192   //    nBins Number of bins
 193   //    vzMin Least @f$z@f$ coordinate of interation point
 194   //    vzMax Largest @f$z@f$ coordinate of interation point
 195   //
 196   if (max < min) {
 197     Double_t tmp = min;
 198     min          = max;
 199     max          = tmp;
 200   }
 201 /*
 202   if (min < -10)
 203     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
 204   if (max > +10)
 205     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
 206 */
 207   fVtxAxis.Set(nBin, min, max);
 208 }
 209 //____________________________________________________________________
 210 void
 211 AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
 212 {
 213   //
 214   // Set the vertex axis to use
 215   //
 216   // Parameters:
 217   //    axis Axis
 218   //
 219   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 220 }
 221
 222 //____________________________________________________________________
 223 void
 224 AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
 225 {
 226   //
 227   // Set the eta axis to use
 228   //
 229   // Parameters:
 230   //    nBins Number of bins
 231   //    vzMin Least @f$\eta@f$
 232   //    vzMax Largest @f$\eta@f$
 233   //
 234   if (max < min) {
 235     Double_t tmp = min;
 236     min          = max;
 237     max          = tmp;
 238   }
 239   if (min < -4)
 240     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
 241   if (max > +6)
 242     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
 243   fEtaAxis.Set(nBin, min, max);
 244 }
 245 //____________________________________________________________________
 246 void
 247 AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
 248 {
 249   //
 250   // Set the eta axis to use
 251   //
 252   // Parameters:
 253   //    axis Axis
 254   //
 255   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 256 }
 257
 258 //____________________________________________________________________
 259 void
 260 AliForwardMCCorrectionsTask::DefineBins(TList* l)
 261 {
 262   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
 263   if (fVtxBins->GetEntries() > 0) return;
 264
 265   fVtxBins->SetOwner();
 266   for (Int_t  i    = 1; i <= fVtxAxis.GetNbins(); i++) {
 267     Double_t  low  = fVtxAxis.GetBinLowEdge(i);
 268     Double_t  high = fVtxAxis.GetBinUpEdge(i);
 269     VtxBin*   bin  = new VtxBin(low, high, fEtaAxis);
 270     fVtxBins->AddAt(bin, i);
 271     bin->CreateOutputObjects(l);
 272   }
 273 }
 274
 275 //____________________________________________________________________
 276 void
 277 AliForwardMCCorrectionsTask::UserCreateOutputObjects()
 278 {
 279   //
 280   // Create output objects
 281   //
 282   //
 283   fList = new TList;
 284   fList->SetOwner();
 285   fList->SetName(Form("%sSums", GetName()));
 286
 287   DefineBins(fList);
 288
 289   fHEvents = new TH1I(GetEventName(false,false),
 290                       "Number of all events",
 291                       fVtxAxis.GetNbins(),
 292                       fVtxAxis.GetXmin(),
 293                       fVtxAxis.GetXmax());
 294   fHEvents->SetXTitle("v_{z} [cm]");
 295   fHEvents->SetYTitle("# of events");
 296   fHEvents->SetFillColor(kBlue+1);
 297   fHEvents->SetFillStyle(3001);
 298   fHEvents->SetDirectory(0);
 299   fList->Add(fHEvents);
 300
 301   fHEventsTr = new TH1I(GetEventName(true, false),
 302                         "Number of triggered events",
 303                         fVtxAxis.GetNbins(),
 304                         fVtxAxis.GetXmin(),
 305                         fVtxAxis.GetXmax());
 306   fHEventsTr->SetXTitle("v_{z} [cm]");
 307   fHEventsTr->SetYTitle("# of events");
 308   fHEventsTr->SetFillColor(kRed+1);
 309   fHEventsTr->SetFillStyle(3001);
 310   fHEventsTr->SetDirectory(0);
 311   fList->Add(fHEventsTr);
 312
 313   fHEventsTrVtx = new TH1I(GetEventName(true,true),
 314                            "Number of events w/trigger and vertex",
 315                            fVtxAxis.GetNbins(),
 316                            fVtxAxis.GetXmin(),
 317                            fVtxAxis.GetXmax());
 318   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 319   fHEventsTrVtx->SetYTitle("# of events");
 320   fHEventsTrVtx->SetFillColor(kBlue+1);
 321   fHEventsTrVtx->SetFillStyle(3001);
 322   fHEventsTrVtx->SetDirectory(0);
 323   fList->Add(fHEventsTrVtx);
 324
 325   // Copy axis objects to output
 326   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
 327                           fVtxAxis.GetNbins(),
 328                           fVtxAxis.GetXmin(),
 329                           fVtxAxis.GetXmax());
 330   TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
 331                           fEtaAxis.GetNbins(),
 332                           fEtaAxis.GetXmin(),
 333                           fEtaAxis.GetXmax());
 334   fList->Add(vtxAxis);
 335   fList->Add(etaAxis);
 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 AliForwardMCCorrectionsTask::UserExec(Option_t*)
 347 {
 348   //
 349   // Process each event
 350   //
 351   // Parameters:
 352   //    option Not used
 353   //
 354
 355   // Get the input data - MC event
 356   AliMCEvent*  mcEvent = MCEvent();
 357   if (!mcEvent) {
 358     AliWarning("No MC event found");
 359     return;
 360   }
 361
 362   // Get the input data - ESD event
 363   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 364   if (!esd) {
 365     AliWarning("No ESD event found for input event");
 366     return;
 367   }
 368
 369   // --- Read in the data --------------------------------------------
 370   LoadBranches();
 371
 372   //--- Read run information -----------------------------------------
 373   if (fFirstEvent && esd->GetESDRun()) {
 374     fInspector.ReadRunDetails(esd);
 375
 376     AliInfo(Form("Initializing with parameters from the ESD:\n"
 377                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 378                  "         AliESDEvent::GetBeamType()     ->%s\n"
 379                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 380                  "         AliESDEvent::GetMagneticField()->%f\n"
 381                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 382                  esd->GetBeamEnergy(),
 383                  esd->GetBeamType(),
 384                  esd->GetCurrentL3(),
 385                  esd->GetMagneticField(),
 386                  esd->GetRunNumber()));
 387
 388     fInspector.SetupForData(fVtxAxis);
 389
 390     Print();
 391     fFirstEvent = false;
 392   }
 393
 394   // Some variables
 395   UInt_t   triggers  = 0;    // Trigger bits
 396   Bool_t   lowFlux   = true; // Low flux flag
 397   UShort_t iVz       = 0;    // Vertex bin from ESD
 398   TVector3 ip(1024,1024,1000);
 399   Double_t cent      = -1;   // Centrality
 400   UShort_t iVzMc     = 0;    // Vertex bin from MC
 401   Double_t vZMc      = 1000; // Z coordinate of IP vertex from MC
 402   Double_t b         = -1;   // Impact parameter
 403   Double_t cMC       = -1;   // Centrality estimate from b
 404   Int_t    nPart     = -1;   // Number of participants
 405   Int_t    nBin      = -1;   // Number of binary collisions
 406   Double_t phiR      = 100;  // Reaction plane from MC
 407   UShort_t nClusters = 0;    // Number of SPD clusters
 408   // Process the data
 409   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip,
 410                                      cent, nClusters);
 411
 412   Bool_t isAccepted = true;
 413   if(fCalculateafterESDeventcuts) {
 414     if (retESD & AliFMDEventInspector::kNoEvent)    isAccepted = false;
 415     if (retESD & AliFMDEventInspector::kNoTriggers) isAccepted = false;
 416     if (retESD & AliFMDEventInspector::kNoVertex)   isAccepted = false;
 417     if (retESD & AliFMDEventInspector::kNoFMD)      isAccepted = false;
 418     if (!isAccepted) return;
 419     // returns if there is not event , does not pass phys selection ,
 420     // has no veretx and lack of FMD data.
 421     // with good veretx outside z range it will continue
 422   }
 423
 424
 425   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
 426                        b, cMC, nPart, nBin, phiR);
 427
 428   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
 429   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
 430
 431   // Fill the event count histograms
 432   if (isInel)           fHEventsTr->Fill(vZMc);
 433   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
 434   fHEvents->Fill(vZMc);
 435
 436   // Now find our vertex bin object
 437   VtxBin* bin = 0;
 438   UShort_t usedZbin=iVzMc;
 439   if(fUseESDVertexCoordinate)
 440         usedZbin=iVz;
 441
 442
 443   if (usedZbin > 0 && usedZbin <= fVtxAxis.GetNbins())
 444     bin = static_cast<VtxBin*>(fVtxBins->At(usedZbin));
 445   if (!bin) {
 446     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
 447     return;
 448   }
 449
 450   // Clear our ESD object
 451   fESDFMD.Clear();
 452
 453   // Get FMD data
 454   AliESDFMD* esdFMD = esd->GetFMDData();
 455
 456   // Now process our input data and store in own ESD object
 457   fTrackDensity.Calculate(*esdFMD, *mcEvent, vZMc, fESDFMD, bin->fPrimary);
 458   bin->fCounts->Fill(0.5);
 459
 460   // And then bin the data in our vtxbin
 461   for (UShort_t d=1; d<=3; d++) {
 462     UShort_t nr = (d == 1 ? 1 : 2);
 463     for (UShort_t q=0; q<nr; q++) {
 464       Char_t      r = (q == 0 ? 'I' : 'O');
 465       UShort_t    ns= (q == 0 ?  20 :  40);
 466       UShort_t    nt= (q == 0 ? 512 : 256);
 467       TH2D*       h = bin->fHists.Get(d,r);
 468
 469       for (UShort_t s=0; s<ns; s++) {
 470         for (UShort_t t=0; t<nt; t++) {
 471           Float_t mult = fESDFMD.Multiplicity(d,r,s,t);
 472
 473           if (mult == 0 || mult > 20) continue;
 474
 475           Float_t phi = fESDFMD.Phi(d,r,s,t) / 180 * TMath::Pi();
 476           Float_t eta = fESDFMD.Eta(d,r,s,t);
 477           h->Fill(eta,phi,mult);
 478         } // for t
 479       } // for s
 480     } // for q
 481   } // for d
 482
 483   PostData(1, fList);
 484 }
 485 //____________________________________________________________________
 486 void
 487 AliForwardMCCorrectionsTask::Terminate(Option_t*)
 488 {
 489   //
 490   // End of job
 491   //
 492   // Parameters:
 493   //    option Not used
 494   //
 495
 496   fList = dynamic_cast<TList*>(GetOutputData(1));
 497   if (!fList) {
 498     AliError("No output list defined");
 499     return;
 500   }
 501
 502   DefineBins(fList);
 503
 504   // Output list
 505   TList* output = new TList;
 506   output->SetOwner();
 507   output->SetName(Form("%sResults", GetName()));
 508
 509   // --- Fill correction object --------------------------------------
 510   AliFMDCorrSecondaryMap* corr = new AliFMDCorrSecondaryMap;
 511   corr->SetVertexAxis(fVtxAxis);
 512   corr->SetEtaAxis(fEtaAxis);
 513
 514   TIter     next(fVtxBins);
 515   VtxBin*   bin = 0;
 516   UShort_t  iVz = 1;
 517   while ((bin = static_cast<VtxBin*>(next())))
 518     bin->Terminate(fList, output, iVz++, corr);
 519
 520   output->Add(corr);
 521
 522   PostData(2, output);
 523 }
 524
 525 //____________________________________________________________________
 526 void
 527 AliForwardMCCorrectionsTask::Print(Option_t* option) const
 528 {
 529   std::cout << ClassName() << "\n"
 530             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
 531             << "  Vertex range:     [" << fVtxAxis.GetXmin()
 532             << "," << fVtxAxis.GetXmax() << "]\n"
 533             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
 534             << "  Eta range:        [" << fEtaAxis.GetXmin()
 535             << "," << fEtaAxis.GetXmax() << "]"
 536             << std::endl;
 537   gROOT->IncreaseDirLevel();
 538   fInspector.Print(option);
 539   fTrackDensity.Print(option);
 540   gROOT->DecreaseDirLevel();
 541 }
 542
 543 //====================================================================
 544 const char*
 545 AliForwardMCCorrectionsTask::VtxBin::BinName(Double_t low,
 546                                              Double_t high)
 547 {
 548   static TString buf;
 549   buf = Form("vtx%+05.1f_%+05.1f", low, high);
 550   buf.ReplaceAll("+", "p");
 551   buf.ReplaceAll("-", "m");
 552   buf.ReplaceAll(".", "d");
 553   return buf.Data();
 554 }
 555
 556
 557 //____________________________________________________________________
 558 AliForwardMCCorrectionsTask::VtxBin::VtxBin()
 559   : fHists(),
 560     fPrimary(0),
 561     fCounts(0)
 562 {
 563 }
 564 //____________________________________________________________________
 565 AliForwardMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
 566                                             Double_t high,
 567                                             const TAxis& axis)
 568   : TNamed(BinName(low, high),
 569            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
 570     fHists(),
 571     fPrimary(0),
 572     fCounts(0)
 573 {
 574   fHists.Init(axis);
 575
 576   fPrimary = new TH2D("primary", "Primaries",
 577                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
 578                       40, 0, 2*TMath::Pi());
 579   fPrimary->SetXTitle("#eta");
 580   fPrimary->SetYTitle("#varphi [radians]");
 581   fPrimary->Sumw2();
 582   fPrimary->SetDirectory(0);
 583
 584   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
 585   fCounts->SetXTitle("Events");
 586   fCounts->SetYTitle("# of Events");
 587   fCounts->SetDirectory(0);
 588 }
 589
 590 //____________________________________________________________________
 591 AliForwardMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
 592   : TNamed(o),
 593     fHists(o.fHists),
 594     fPrimary(0),
 595     fCounts(0)
 596 {
 597   if (o.fPrimary) {
 598     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 599     fPrimary->SetDirectory(0);
 600   }
 601   if (o.fCounts) {
 602     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 603     fCounts->SetDirectory(0);
 604   }
 605 }
 606
 607 //____________________________________________________________________
 608 AliForwardMCCorrectionsTask::VtxBin&
 609 AliForwardMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
 610 {
 611   if (&o == this) return *this;
 612   TNamed::operator=(o);
 613   fHists   = o.fHists;
 614   fPrimary = 0;
 615   fCounts  = 0;
 616   if (o.fPrimary) {
 617     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 618     fPrimary->SetDirectory(0);
 619   }
 620   if (o.fCounts) {
 621     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 622     fCounts->SetDirectory(0);
 623   }
 624   return *this;
 625 }
 626
 627 //____________________________________________________________________
 628 void
 629 AliForwardMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
 630 {
 631   TList* d = new TList;
 632   d->SetName(GetName());
 633   d->SetOwner();
 634   l->Add(d);
 635
 636   d->Add(fHists.fFMD1i);
 637   d->Add(fHists.fFMD2i);
 638   d->Add(fHists.fFMD2o);
 639   d->Add(fHists.fFMD3i);
 640   d->Add(fHists.fFMD3o);
 641
 642   d->Add(fPrimary);
 643   d->Add(fCounts);
 644 }
 645
 646 //____________________________________________________________________
 647 TH2D*
 648 AliForwardMCCorrectionsTask::VtxBin::MakeBg(const TH2D* hits,
 649                                             const TH2D* primary) const
 650 {
 651   TH2D* h = static_cast<TH2D*>(hits->Clone());
 652   h->SetDirectory(0);
 653   TString n(h->GetName());
 654   n.ReplaceAll("_cache", "");
 655   h->SetName(n);
 656   h->Divide(primary);
 657
 658   return h;
 659 }
 660
 661 //____________________________________________________________________
 662 void
 663 AliForwardMCCorrectionsTask::VtxBin::Terminate(const TList* input,
 664                                             TList* output,
 665                                             UShort_t iVz,
 666                                             AliFMDCorrSecondaryMap* map)
 667 {
 668   TList* out = new TList;
 669   out->SetName(GetName());
 670   out->SetOwner();
 671   output->Add(out);
 672
 673   TList* l = static_cast<TList*>(input->FindObject(GetName()));
 674   if (!l) {
 675     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
 676     return;
 677   }
 678
 679   TH2D*   fmd1i = static_cast<TH2D*>(l->FindObject("FMD1I_cache"));
 680   TH2D*   fmd2i = static_cast<TH2D*>(l->FindObject("FMD2I_cache"));
 681   TH2D*   fmd2o = static_cast<TH2D*>(l->FindObject("FMD2O_cache"));
 682   TH2D*   fmd3i = static_cast<TH2D*>(l->FindObject("FMD3I_cache"));
 683   TH2D*   fmd3o = static_cast<TH2D*>(l->FindObject("FMD3O_cache"));
 684   TH2D*   primO = static_cast<TH2D*>(l->FindObject("primary"));
 685   if (!fmd1i || !fmd2i || !fmd2o || !fmd3i || !fmd3o || !primO) {
 686     AliError(Form("Missing histogram(s): %p,%p,%p,%p,%p,%p",
 687                   fmd1i, fmd2i, fmd2o, fmd3i, fmd3o, primO));
 688     return;
 689   }
 690
 691   // Half coverage in phi for inners
 692   TH2D*   primI = static_cast<TH2D*>(primO->Clone());
 693   primI->SetDirectory(0);
 694   primI->RebinY(2);
 695
 696   TH2D* bg1i = MakeBg(fmd1i, primI);
 697   TH2D* bg2i = MakeBg(fmd2i, primI);
 698   TH2D* bg2o = MakeBg(fmd2o, primO);
 699   TH2D* bg3i = MakeBg(fmd3i, primI);
 700   TH2D* bg3o = MakeBg(fmd3o, primO);
 701   map->SetCorrection(1, 'I', iVz, bg1i);
 702   map->SetCorrection(2, 'I', iVz, bg2i);
 703   map->SetCorrection(2, 'O', iVz, bg2o);
 704   map->SetCorrection(3, 'I', iVz, bg3i);
 705   map->SetCorrection(3, 'O', iVz, bg3o);
 706   out->Add(bg1i);
 707   out->Add(bg2i);
 708   out->Add(bg2o);
 709   out->Add(bg3i);
 710   out->Add(bg3o);
 711
 712 }
 713
 714 //
 715 // EOF
 716 //