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