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