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