PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx

   1 //====================================================================
   2 //
   3 // Base class for classes that calculate the multiplicity in the
   4 // central region event-by-event
   5 //
   6 // Inputs:
   7 //   - AliESDEvent
   8 //
   9 // Outputs:
  10 //   - AliAODCentralMult
  11 //
  12 // Histograms
  13 //
  14 // Corrections used
  15 #include "AliCentralMultiplicityTask.h"
  16 #include "AliAODForwardMult.h"
  17 #include "AliForwardUtil.h"
  18 #include "AliLog.h"
  19 #include "AliAODHandler.h"
  20 #include "AliAnalysisManager.h"
  21 #include "AliESDEvent.h"
  22 #include "AliMultiplicity.h"
  23 #include <TROOT.h>
  24 #include <TSystem.h>
  25 #include <TFile.h>
  26 #include <TError.h>
  27 #include <TSystem.h>
  28 #include <iostream>
  29 #include <iomanip>
  30
  31 //====================================================================
  32 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
  33   : AliAnalysisTaskSE(name),
  34     fInspector("centralEventInspector"),
  35     fData(0),
  36     fList(0),
  37     fHits(0),
  38     fAODCentral(kFALSE),
  39     fManager(),
  40     fUseSecondary(true),
  41     fUseAcceptance(true),
  42     fFirstEventSeen(false),
  43     fIvz(0),
  44     fNClusterTracklet(0),
  45     fClusterPerTracklet(0),
  46     fNCluster(0),
  47     fNTracklet(0),
  48     fEtaMin(0),
  49     fEtaMax(0)
  50 {
  51   //
  52   // Constructor
  53   //
  54   DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
  55   DefineOutput(1, TList::Class());
  56   fBranchNames =
  57     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  58     "SPDVertex.,PrimaryVertex.";
  59 }
  60 //____________________________________________________________________
  61 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
  62   : AliAnalysisTaskSE(),
  63     fInspector(),
  64     fData(0),
  65     fList(0),
  66     fHits(0),
  67     fAODCentral(),
  68     fManager(),
  69     fUseSecondary(true),
  70     fUseAcceptance(true),
  71     fFirstEventSeen(false),
  72     fIvz(0),
  73     fNClusterTracklet(0),
  74     fClusterPerTracklet(0),
  75     fNCluster(0),
  76     fNTracklet(0),
  77     fEtaMin(0),
  78     fEtaMax(0)
  79 {
  80   //
  81   // Constructor
  82   //
  83   DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
  84 }
  85 //____________________________________________________________________
  86 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
  87   : AliAnalysisTaskSE(o),
  88     fInspector(o.fInspector),
  89     fData(o.fData),
  90     fList(o.fList),
  91     fHits(o.fHits),
  92     fAODCentral(o.fAODCentral),
  93     fManager(o.fManager),
  94     fUseSecondary(o.fUseSecondary),
  95     fUseAcceptance(o.fUseAcceptance),
  96     fFirstEventSeen(o.fFirstEventSeen),
  97     fIvz(0),
  98     fNClusterTracklet(o.fNClusterTracklet),
  99     fClusterPerTracklet(o.fClusterPerTracklet),
 100     fNCluster(o.fNCluster),
 101     fNTracklet(o.fNTracklet),
 102     fEtaMin(o.fEtaMin),
 103     fEtaMax(o.fEtaMax)
 104 {
 105   //
 106   // Copy constructor
 107   //
 108   DGUARD(fDebug, 3,"COPY CTOR of AliCentralMultiplicityTask");
 109
 110 }
 111 //____________________________________________________________________
 112 AliCentralMultiplicityTask&
 113 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
 114 {
 115   //
 116   // Assignment operator
 117   //
 118   DGUARD(fDebug,3,"Assignment of AliCentralMultiplicityTask");
 119   if (&o == this) return *this;
 120   fInspector         = o.fInspector;
 121   fData              = o.fData;
 122   fList              = o.fList;
 123   fHits              = o.fHits;
 124   fAODCentral        = o.fAODCentral;
 125   fManager           = o.fManager;
 126   fUseSecondary      = o.fUseSecondary;
 127   fUseAcceptance     = o.fUseAcceptance;
 128   fFirstEventSeen    = o.fFirstEventSeen;
 129   fIvz               = 0;
 130   fNClusterTracklet  = o.fNClusterTracklet;
 131   fClusterPerTracklet= o.fClusterPerTracklet;
 132   fNCluster          = o.fNCluster;
 133   fNTracklet         = o.fNTracklet;
 134   fEtaMin            = o.fEtaMin;
 135   fEtaMax            = o.fEtaMax;
 136   return *this;
 137 }
 138 //____________________________________________________________________
 139 Bool_t
 140 AliCentralMultiplicityTask::Configure(const char* macro)
 141 {
 142   // --- Configure the task ------------------------------------------
 143   TString macroPath(gROOT->GetMacroPath());
 144   if (!macroPath.Contains("$(ALICE_ROOT)/PWGLF/FORWARD/analysis2")) {
 145     macroPath.Append(":$(ALICE_ROOT)/PWGLF/FORWARD/analysis2");
 146     gROOT->SetMacroPath(macroPath);
 147   }
 148   const char* config = gSystem->Which(gROOT->GetMacroPath(),macro);
 149   if (!config) {
 150     AliWarningF("%s not found in %s", macro, gROOT->GetMacroPath());
 151     return false;
 152   }
 153
 154   AliInfoF("Loading configuration of '%s' from %s", ClassName(), config);
 155   gROOT->Macro(Form("%s((AliCentralMultiplicityTask*)%p)", config, this));
 156   delete config;
 157
 158   return true;
 159 }
 160
 161 //____________________________________________________________________
 162 void AliCentralMultiplicityTask::UserCreateOutputObjects()
 163 {
 164   //
 165   // Create output objects
 166   //
 167   //
 168   DGUARD(fDebug,1,"Create user output in AliCentralMultiplicityTask");
 169
 170   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 171   AliAODHandler*      ah =
 172     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 173   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 174
 175
 176   TObject* obj = &fAODCentral;
 177   ah->AddBranch("AliAODCentralMult", &obj);
 178
 179
 180
 181   fList = new TList();
 182   fList->SetOwner();
 183
 184   fInspector.CreateOutputObjects(fList);
 185
 186   PostData(1,fList);
 187 }
 188
 189 //____________________________________________________________________
 190 AliESDEvent*
 191 AliCentralMultiplicityTask::GetESDEvent()
 192 {
 193   //
 194   // Get the ESD event. IF this is the first event, initialise
 195   //
 196   DGUARD(fDebug,1,"Get ESD event in AliCentralMultiplicityTask");
 197   if (IsZombie()) return 0;
 198   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 199   if (!esd) {
 200     AliWarning("No ESD event found for input event");
 201     return 0;
 202   }
 203
 204   // Load in the data needed
 205   LoadBranches();
 206
 207   // IF we've read the first event already, just return the event
 208   if (fFirstEventSeen) return esd;
 209
 210   // Read the details of the rung
 211   fInspector.ReadRunDetails(esd);
 212
 213   // If we weren't initialised before (i.e., in the setup), do so now.
 214   if (!GetManager().IsInit()) {
 215     GetManager().Init(fInspector.GetCollisionSystem(),
 216                       fInspector.GetEnergy(),
 217                       fInspector.GetField());
 218     //AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
 219   }
 220   Bool_t ok = true;
 221   if (!GetManager().HasSecondaryCorrection()) {
 222     ok = false;
 223     AliError("No secondary correction defined!");
 224   }
 225   if (!GetManager().HasAcceptanceCorrection()) {
 226     ok = false;
 227     AliError("No acceptance correction defined!");
 228   }
 229   // If the corrections are not seen, make this a zombie, and prevent
 230   // further execution of this task.
 231   if (!ok) {
 232     AliError("Missing corrections, make this a zombie");
 233     SetZombie(true);
 234     esd = 0;
 235     fFirstEventSeen = true;
 236     return esd;
 237   }
 238
 239   // Check for existence and get secondary map
 240   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
 241   const TAxis& vaxis = secMap->GetVertexAxis();
 242
 243   FindEtaLimits();
 244
 245   fNClusterTracklet = new TH2D("nClusterVsnTracklet",
 246                                "Total number of cluster vs number of tracklets",
 247                                100, 0, 10000, 100, 0, 10000);
 248   fNClusterTracklet->SetDirectory(0);
 249   fNClusterTracklet->SetXTitle("# of free clusters");
 250   fNClusterTracklet->SetYTitle("# of tracklets");
 251   fNClusterTracklet->SetStats(0);
 252   fList->Add(fNClusterTracklet);
 253
 254   Int_t    nEta = 80;
 255   Double_t lEta = 2;
 256   fClusterPerTracklet = new TH2D("clusterPerTracklet",
 257                                  "N_{free cluster}/N_{tracklet} vs. #eta",
 258                                  nEta,-lEta,lEta, 101, -.05, 10.05);
 259   fClusterPerTracklet->SetDirectory(0);
 260   fClusterPerTracklet->SetXTitle("#eta");
 261   fClusterPerTracklet->SetYTitle("N_{free cluster}/N_{tracklet}");
 262   fClusterPerTracklet->SetStats(0);
 263   fList->Add(fClusterPerTracklet);
 264
 265   // Cache histograms
 266   fNCluster = new TH1D("cacheCluster", "", nEta,-lEta,lEta);
 267   fNCluster->SetDirectory(0);
 268   fNCluster->Sumw2();
 269
 270   fNTracklet = new TH1D("cacheTracklet", "", nEta,-lEta,lEta);
 271   fNTracklet->SetDirectory(0);
 272   fNTracklet->Sumw2();
 273
 274   // Initialize the inspecto
 275   fInspector.SetupForData(vaxis);
 276   fFirstEventSeen = kTRUE;
 277
 278   // Print some information
 279   Print();
 280
 281   return esd;
 282 }
 283 //____________________________________________________________________
 284 void
 285 AliCentralMultiplicityTask::MarkEventForStore() const
 286 {
 287   // Make sure the AOD tree is filled
 288   DGUARD(fDebug,1,"Mark AOD event for store in AliCentralMultiplicityTask");
 289   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 290   AliAODHandler*      ah =
 291     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 292   if (!ah)
 293     AliFatal("No AOD output handler set in analysis manager");
 294
 295   ah->SetFillAOD(kTRUE);
 296 }
 297
 298 //____________________________________________________________________
 299 void AliCentralMultiplicityTask::FindEtaLimits()
 300 {
 301   // Find our pseudo-rapidity limits
 302   //
 303   // Uses the secondary map to do so.
 304   DGUARD(fDebug,1,"Find eta limits in AliCentralMultiplicityTask");
 305   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
 306
 307   const TAxis& vaxis = secMap->GetVertexAxis();
 308
 309   fEtaMin.Set(vaxis.GetNbins());
 310   fEtaMax.Set(vaxis.GetNbins());
 311
 312   fHits = new TList;
 313   fHits->SetOwner();
 314   fHits->SetName("hitMaps");
 315   fList->Add(fHits);
 316
 317   TList* secs = new TList;
 318   secs->SetOwner();
 319   secs->SetName("secondaryMaps");
 320   fList->Add(secs);
 321   unsigned short s = 1;
 322   TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
 323                              secMap->GetCorrection(s)->GetXaxis()->GetNbins(),
 324                              secMap->GetCorrection(s)->GetXaxis()->GetXmin(),
 325                              secMap->GetCorrection(s)->GetXaxis()->GetXmax(),
 326                              vaxis.GetNbins(),vaxis.GetXmin(),vaxis.GetXmax());
 327   hCoverage->SetDirectory(0);
 328   hCoverage->SetXTitle("#eta");
 329   hCoverage->SetYTitle("v_{z} [cm]");
 330   hCoverage->SetZTitle("n_{bins}");
 331   fList->Add(hCoverage);
 332
 333   fAODCentral.Init(*(secMap->GetCorrection(s)->GetXaxis()));
 334
 335   for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
 336     TH2D* corr = secMap->GetCorrection(UShort_t(v));
 337     TH1D* proj = corr->ProjectionX(Form("secCor%02d", v));
 338     proj->Scale(1. / corr->GetNbinsY());
 339     proj->SetTitle(Form("Projection of secondary correction "
 340                         "for %+5.1f<v_{z}<%+5.1f",
 341                         vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
 342     proj->SetYTitle("#LT 2^{nd} correction#GT");
 343     proj->SetDirectory(0);
 344     proj->SetMarkerStyle(20);
 345     proj->SetMarkerColor(kBlue+1);
 346     secs->Add(proj);
 347
 348     TH2D* obg = static_cast<TH2D*>(corr->Clone(Form("secCor2DFiducial%02d",v)));
 349     obg->SetDirectory(0);
 350     secs->Add(obg);
 351
 352     TH1D* after = static_cast<TH1D*>(proj->Clone(Form("secCorFiducial%02d",v)));
 353     after->SetDirectory(0);
 354     after->SetMarkerColor(kRed+1);
 355     secs->Add(after);
 356
 357     TH2D* data = static_cast<TH2D*>(corr->Clone(Form("hitMap%02d",v)));
 358     //d->SetTitle(Form("hitMap%02d",v));
 359     data->SetTitle(Form("d^{2}N/d#eta d#phi "
 360                         "for %+5.1f<v_{z}<%+5.1f",
 361                         vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
 362     data->GetZaxis()->SetTitle("");
 363     data->SetMarkerColor(kBlack);
 364     data->SetMarkerStyle(1);
 365     fHits->Add(data);
 366
 367     TH1D* hAcceptance = fManager.GetAcceptanceCorrection(v);
 368     TH1D* accClone   = static_cast<TH1D*>(hAcceptance->Clone(Form("acceptance%02d",v)));
 369     secs->Add(accClone);
 370
 371     // Double_t prev = 0;
 372     for (Int_t e = 1; e <= proj->GetNbinsX(); e++) {
 373       Double_t c = proj->GetBinContent(e);
 374       if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
 375         fEtaMin[v-1] = e;
 376         break;
 377       }
 378       // prev = c;
 379       after->SetBinContent(e, 0);
 380       after->SetBinError(e, 0);
 381       for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
 382         obg->SetBinContent(e,nn,0);
 383
 384     }
 385     for (Int_t e = proj->GetNbinsX(); e >= 1; e--) {
 386       Double_t c = proj->GetBinContent(e);
 387       if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
 388         fEtaMax[v-1] = e;
 389         break;
 390       }
 391       // prev = c;
 392       after->SetBinContent(e, 0);
 393       after->SetBinError(e, 0);
 394       for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
 395         obg->SetBinContent(e,nn,0);
 396
 397     }
 398
 399     for (Int_t nn = fEtaMin[v-1]; nn<=fEtaMax[v-1]; nn++) {
 400       hCoverage->SetBinContent(nn,v,1);
 401     }
 402
 403   }
 404 }
 405
 406 //____________________________________________________________________
 407 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
 408 {
 409   //
 410   // Process each event
 411   //
 412   // Parameters:
 413   //    option Not used
 414   //
 415   DGUARD(fDebug,1,"Process event in AliCentralMultiplicityTask");
 416   fAODCentral.Clear("");
 417   fIvz = 0;
 418
 419   AliESDEvent* esd = GetESDEvent();
 420   if (!esd) return;
 421
 422   Bool_t   lowFlux   = kFALSE;
 423   UInt_t   triggers  = 0;
 424   UShort_t ivz       = 0;
 425   TVector3 ip;
 426   Double_t cent      = -1;
 427   UShort_t nClusters = 0;
 428   UInt_t   found     = fInspector.Process(esd, triggers, lowFlux,
 429                                           ivz, ip, cent, nClusters);
 430
 431   // No event or no trigger
 432   if (found & AliFMDEventInspector::kNoEvent)    return;
 433   if (found & AliFMDEventInspector::kNoTriggers) return;
 434
 435   // Make sure AOD is filled
 436   MarkEventForStore();
 437
 438   if (found == AliFMDEventInspector::kNoSPD)      return;
 439   if (found == AliFMDEventInspector::kNoVertex)   return;
 440   if (triggers & AliAODForwardMult::kPileUp)      return;
 441   if (found == AliFMDEventInspector::kBadVertex)  return; // Out of range
 442
 443   //Doing analysis
 444   fIvz = ivz;
 445   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 446
 447   TH2D& aodHist = fAODCentral.GetHistogram();
 448
 449   ProcessESD(aodHist, spdmult);
 450   CorrectData(aodHist, ivz);
 451   //Producing hit maps
 452   // TList* hitList = static_cast<TList*>(fList->FindObject("hitMaps"));
 453   TH2D* data  = static_cast<TH2D*>(fHits->At(ivz-1));
 454   if(data) data->Add(&aodHist);
 455
 456   PostData(1,fList);
 457 }
 458 //____________________________________________________________________
 459 void
 460 AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
 461                                        const AliMultiplicity* spdmult) const
 462 {
 463   DGUARD(fDebug,1,"Process the ESD in AliCentralMultiplicityTask");
 464   fNTracklet->Reset();
 465   fNCluster->Reset();
 466
 467   //Filling clusters in layer 1 used for tracklets...
 468   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
 469     Double_t eta = spdmult->GetEta(j);
 470     fNTracklet->Fill(eta);
 471     aodHist.Fill(eta,spdmult->GetPhi(j));
 472   }
 473
 474   //...and then the unused ones in layer 1
 475   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 476     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 477     fNCluster->Fill(eta);
 478     aodHist.Fill(eta, spdmult->GetPhiSingle(j));
 479   }
 480   fNClusterTracklet->Fill(fNCluster->GetEntries(),
 481                           fNTracklet->GetEntries());
 482
 483   fNCluster->Divide(fNTracklet);
 484   for (Int_t j = 1; j <= fNCluster->GetNbinsX(); j++)
 485     fClusterPerTracklet->Fill(fNCluster->GetXaxis()->GetBinCenter(j),
 486                               fNCluster->GetBinContent(j));
 487
 488 }
 489
 490 //____________________________________________________________________
 491 void
 492 AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const
 493 {
 494   // Corrections
 495   DGUARD(fDebug,1,"Correct data in AliCentralMultiplicityTask");
 496   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
 497   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
 498
 499   if (!hSecMap)     AliFatal("No secondary map!");
 500   if (!hAcceptance) AliFatal("No acceptance!");
 501
 502   if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
 503
 504   Int_t nY = aodHist.GetNbinsY();
 505   for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) {
 506     Float_t accCor = hAcceptance->GetBinContent(ix);
 507     Float_t accErr = hAcceptance->GetBinError(ix);
 508
 509     Bool_t fiducial = true;
 510     if (ix < fEtaMin[vtxbin-1] || ix > fEtaMax[vtxbin-1])
 511       fiducial = false;
 512     //  Bool_t etabinSeen = kFALSE;
 513     for(Int_t iy = 1; iy <= nY; iy++) {
 514 #if 1
 515       if (!fiducial) {
 516         aodHist.SetBinContent(ix, iy, 0);
 517         aodHist.SetBinError(ix, iy, 0);
 518         continue;
 519       }
 520 #endif
 521       // Get currrent value
 522       Float_t aodValue = aodHist.GetBinContent(ix,iy);
 523       Float_t aodErr   = aodHist.GetBinError(ix,iy);
 524
 525 #if 0 // This is done once in the FindEtaBins function
 526       // Set underflow bin
 527       Float_t secCor   = 0;
 528       if(hSecMap)       secCor     = hSecMap->GetBinContent(ix,iy);
 529       if (secCor > 0.5) etabinSeen = kTRUE;
 530 #endif
 531       if (aodValue < 0.000001) {
 532         aodHist.SetBinContent(ix,iy, 0);
 533         aodHist.SetBinError(ix,iy, 0);
 534         continue;
 535       }
 536       if (!fUseAcceptance) continue;
 537
 538       // Acceptance correction
 539       if (accCor   < 0.000001) accCor = 1;
 540       Float_t aodNew   = aodValue / accCor ;
 541       Float_t error    = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
 542                                             TMath::Power(accErr/accCor,2) );
 543       aodHist.SetBinContent(ix,iy, aodNew);
 544       //test
 545       aodHist.SetBinError(ix,iy,error);
 546       aodHist.SetBinError(ix,iy,aodErr);
 547     }
 548     //Filling underflow bin if we eta bin is in range
 549     if (fiducial) {
 550       aodHist.SetBinContent(ix,0, 1.);
 551       aodHist.SetBinContent(ix,nY+1, 1.);
 552     }
 553     // if (etabinSeen) aodHist.SetBinContent(ix,0, 1.);
 554   }
 555 }
 556
 557 //____________________________________________________________________
 558 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
 559 {
 560   //
 561   // End of job
 562   //
 563   // Parameters:
 564   //    option Not used
 565   //
 566   DGUARD(fDebug,1,"Process merged output in AliCentralMultiplicityTask");
 567 }
 568 //____________________________________________________________________
 569 void
 570 AliCentralMultiplicityTask::Print(Option_t* option) const
 571 {
 572   //
 573   // Print information
 574   //
 575   // Parameters:
 576   //    option Not used
 577   //
 578   std::cout << ClassName() << ": " << GetName() << "\n"
 579             << std::boolalpha
 580             << "  Use secondary correction:  " << fUseSecondary << '\n'
 581             << "  Use acceptance correction: " << fUseAcceptance << '\n'
 582             << "  Off-line trigger mask:  0x"
 583             << std::hex     << std::setfill('0')
 584             << std::setw (8) << fOfflineTriggerMask
 585             << std::dec     << std::setfill (' ')
 586             << std::noboolalpha << std::endl;
 587   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
 588   if (secMap) {
 589     const TAxis& vaxis = secMap->GetVertexAxis();
 590     std::cout << "  Eta ranges:\n"
 591             << "     Vertex        | Eta bins\n"
 592               << "   bin     range   | \n"
 593               << "   ----------------+-----------" << std::endl;
 594     for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
 595       std::cout << "   " << std::setw(2) << v << "  "
 596                 << std::setw(5) << vaxis.GetBinLowEdge(v) << "-"
 597                 << std::setw(5) << vaxis.GetBinUpEdge(v) << " | ";
 598       if (fEtaMin.GetSize() <= 0)
 599         std::cout << " ? -   ?";
 600       else
 601         std::cout << std::setw(3) << fEtaMin[v-1] << "-"
 602                   << std::setw(3) << fEtaMax[v-1];
 603       std::cout << std::endl;
 604     }
 605   }
 606
 607   gROOT->IncreaseDirLevel();
 608   fManager.Print(option);
 609   fInspector.Print(option);
 610   gROOT->DecreaseDirLevel();
 611
 612 }
 613 //====================================================================
 614 AliCentralMultiplicityTask::Manager::Manager() :
 615   fAcceptancePath("$ALICE_ROOT/PWGLF/FORWARD/corrections/CentralAcceptance"),
 616   fSecMapPath("$ALICE_ROOT/PWGLF/FORWARD/corrections/CentralSecMap"),
 617   fAcceptance(),
 618   fSecmap(),
 619   fAcceptanceName("centralacceptance"),
 620   fSecMapName("centralsecmap"),
 621   fIsInit(kFALSE)
 622 {
 623   //
 624   // Constructor
 625   //
 626 }
 627 //____________________________________________________________________
 628 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
 629   :fAcceptancePath(o.fAcceptancePath),
 630    fSecMapPath(o.fSecMapPath),
 631    fAcceptance(o.fAcceptance),
 632    fSecmap(o.fSecmap),
 633    fAcceptanceName(o.fAcceptanceName),
 634    fSecMapName(o.fSecMapName),
 635    fIsInit(o.fIsInit)
 636 {
 637   //
 638   // Copy Constructor
 639   //
 640 }
 641 //____________________________________________________________________
 642 AliCentralMultiplicityTask::Manager&
 643 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
 644 {
 645   //
 646   // Assignment operator
 647   //
 648   if (&o == this) return *this;
 649   fAcceptancePath = o.fAcceptancePath;
 650   fSecMapPath     = o.fSecMapPath;
 651   fAcceptance     = o.fAcceptance;
 652   fSecmap         = o.fSecmap;
 653   fAcceptanceName = o.fAcceptanceName;
 654   fSecMapName     = o.fSecMapName;
 655   fIsInit         = o.fIsInit;
 656   return *this;
 657 }
 658
 659 //____________________________________________________________________
 660 const char*
 661 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
 662                                                      UShort_t sys,
 663                                                      UShort_t sNN,
 664                                                      Short_t  field) const
 665 {
 666   //
 667   // Get full path name to object file
 668   //
 669   // Parameters:
 670   //    what   What to get
 671   //    sys    Collision system
 672   //    sNN    Center of mass energy
 673   //    field  Magnetic field
 674   //
 675   // Return:
 676   //
 677   //
 678   return Form("%s/%s",
 679               what == 0 ? GetSecMapPath() : GetAcceptancePath(),
 680               GetFileName(what, sys, sNN, field));
 681 }
 682
 683 //____________________________________________________________________
 684 const char*
 685 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
 686                                                  UShort_t  sys,
 687                                                  UShort_t  sNN,
 688                                                  Short_t   field) const
 689 {
 690   //
 691   // Get the full path name
 692   //
 693   // Parameters:
 694   //    what   What to get
 695   //    sys    Collision system
 696   //    sNN    Center of mass energy
 697   //    field  Magnetic field
 698   //
 699   // Return:
 700   //
 701   //
 702   // Must be static - otherwise the data may disappear on return from
 703   // this member function
 704   static TString fname = "";
 705
 706   switch(what) {
 707   case 0:  fname = fSecMapName;     break;
 708   case 1:  fname = fAcceptanceName; break;
 709   default:
 710     ::Error("GetFileName",
 711             "Invalid indentifier %d for central object, must be 0 or 1!", what);
 712     break;
 713   }
 714   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
 715                     AliForwardUtil::CollisionSystemString(sys),
 716                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
 717
 718   return fname.Data();
 719 }
 720
 721 //____________________________________________________________________
 722 TH2D*
 723 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
 724 {
 725   //
 726   // Get the secondary map
 727   //
 728   // Parameters:
 729   //    vtxbin
 730   //
 731   // Return:
 732   //
 733   //
 734   if (!fSecmap) {
 735     ::Warning("GetSecMapCorrection","No secondary map defined");
 736     return 0;
 737   }
 738   return fSecmap->GetCorrection(vtxbin);
 739 }
 740 //____________________________________________________________________
 741 TH1D*
 742 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
 743   const
 744 {
 745   //
 746   // Get the acceptance correction
 747   //
 748   // Parameters:
 749   //    vtxbin
 750   //
 751   // Return:
 752   //
 753   //
 754   if (!fAcceptance) {
 755     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
 756     return 0;
 757   }
 758   return fAcceptance->GetCorrection(vtxbin);
 759 }
 760
 761 //____________________________________________________________________
 762 void
 763 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
 764                                           UShort_t  sNN,
 765                                           Short_t   field)
 766 {
 767   //
 768   // Initialize
 769   //
 770   // Parameters:
 771   //    sys    Collision system (1: pp, 2: PbPb, 3: pPb)
 772   //    sNN    Center of mass energy per nucleon pair [GeV]
 773   //    field  Magnetic field [kG]
 774   //
 775   if(fIsInit) ::Warning("Init","Already initialised - overriding...");
 776
 777   TFile fsec(GetFullFileName(0,sys,sNN,field));
 778   fSecmap =
 779     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
 780   if(!fSecmap) {
 781     ::Error("Init", "no central Secondary Map found!") ;
 782     return;
 783   }
 784   TFile facc(GetFullFileName(1,sys,sNN,field));
 785   fAcceptance =
 786     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
 787   if(!fAcceptance) {
 788     ::Error("Init", "no central Acceptance found!") ;
 789     return;
 790   }
 791
 792   if(fSecmap && fAcceptance) {
 793     fIsInit = kTRUE;
 794     ::Info("Init",
 795            "Central Manager initialised for %s, energy %dGeV, field %dkG",
 796            sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "unknown", sNN,field);
 797   }
 798 }
 799 //____________________________________________________________________
 800 Bool_t
 801 AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what,
 802                                                UShort_t sys,
 803                                                UShort_t sNN,
 804                                                Short_t  fld,
 805                                                TObject* obj,
 806                                                Bool_t   full) const
 807 {
 808   //
 809   // Write correction output to (a temporary) file
 810   //
 811   // Parameters:
 812   //   What     What to write
 813   //   sys      Collision system (1: pp, 2: PbPb, 3: pPb)
 814   //   sNN      Center of mass energy per nucleon (GeV)
 815   //   fld      Field (kG)
 816   //   obj      Object to write
 817   //   full     if true, write to full path, otherwise locally
 818   //
 819   // Return:
 820   //   true on success.
 821   TString ofName;
 822   if (!full)
 823     ofName = GetFileName(what, sys, sNN, fld);
 824   else
 825     ofName = GetFullFileName(what, sys, sNN, fld);
 826   if (ofName.IsNull()) {
 827     AliErrorGeneral("Manager",Form("Unknown object type %d", what));
 828     return false;
 829   }
 830   TFile* output = TFile::Open(ofName, "RECREATE");
 831   if (!output) {
 832     AliErrorGeneral("Manager",Form("Failed to open file %s", ofName.Data()));
 833     return false;
 834   }
 835
 836   TString oName(GetObjectName(what));
 837   Int_t ret = obj->Write(oName);
 838   if (ret <= 0) {
 839     AliErrorGeneral("Manager",Form("Failed to write %p to %s/%s (%d)",
 840                                    obj, ofName.Data(), oName.Data(), ret));
 841     return false;
 842   }
 843
 844   ret = output->Write();
 845   if (ret < 0) {
 846     AliErrorGeneral("Manager",
 847                     Form("Failed to write %s to disk (%d)", ofName.Data(),ret));
 848     return false;
 849   }
 850   // output->ls();
 851   output->Close();
 852
 853 #if 0
 854   TString cName(obj->IsA()->GetName());
 855   AliInfoGeneral("Manager",
 856                  Form("Wrote %s object %s to %s\n",
 857                       cName.Data(),oName.Data(), ofName.Data()));
 858   if (!full) {
 859     TString dName(GetFileDir(what));
 860     AliInfoGeneral("Manager",
 861                    Form("\n  %s should be copied to %s\n"
 862                         "Do for example\n\t"
 863                         "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/"
 864                         "MoveCorrections.C\\(%d\\)\nor\n\t"
 865                         "cp %s %s/",
 866                         ofName.Data(),dName.Data(),
 867                         what, ofName.Data(),
 868                         gSystem->ExpandPathName(dName.Data())));
 869
 870
 871   }
 872 #endif
 873   return true;
 874 }
 875
 876 //____________________________________________________________________
 877 void
 878 AliCentralMultiplicityTask::Manager::Print(Option_t* option) const
 879 {
 880   //
 881   // Print information to standard output
 882   //
 883   std::cout << " AliCentralMultiplicityTask::Manager\n"
 884             << std::boolalpha
 885             << "  Initialized:     " << fIsInit << '\n'
 886             << "  Acceptance path: " << fAcceptancePath << '\n'
 887             << "  Acceptance name: " << fAcceptanceName << '\n'
 888             << "  Acceptance:      " << fAcceptance << '\n'
 889             << "  Secondary path:  " << fSecMapPath << '\n'
 890             << "  Secondary name:  " << fSecMapName << '\n'
 891             << "  Secondary map:   " << fSecmap
 892             << std::noboolalpha << std::endl;
 893   if (fAcceptance) fAcceptance->Print(option);
 894   if (fSecmap)     fSecmap->Print(option);
 895 }
 896
 897 //
 898 // EOF
 899 //