PWG2/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx

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