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) {
 181     Double_t tmp = min;
 182     min          = max;
 183     max          = tmp;
 184   }
 185   if (min < -10)
 186     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
 187   if (max > +10)
 188     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
 189   fVtxAxis.Set(nBin, min, max);
 190 }
 191 //____________________________________________________________________
 192 void
 193 AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
 194 {
 195   //
 196   // Set the vertex axis to use
 197   //
 198   // Parameters:
 199   //    axis Axis
 200   //
 201   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 202 }
 203
 204 //____________________________________________________________________
 205 void
 206 AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
 207 {
 208   //
 209   // Set the eta axis to use
 210   //
 211   // Parameters:
 212   //    nBins Number of bins
 213   //    vzMin Least @f$\eta@f$
 214   //    vzMax Largest @f$\eta@f$
 215   //
 216   if (max < min) {
 217     Double_t tmp = min;
 218     min          = max;
 219     max          = tmp;
 220   }
 221   if (min < -4)
 222     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
 223   if (max > +6)
 224     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
 225   fEtaAxis.Set(nBin, min, max);
 226 }
 227 //____________________________________________________________________
 228 void
 229 AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
 230 {
 231   //
 232   // Set the eta axis to use
 233   //
 234   // Parameters:
 235   //    axis Axis
 236   //
 237   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 238 }
 239
 240 //____________________________________________________________________
 241 void
 242 AliCentralMCCorrectionsTask::DefineBins(TList* l)
 243 {
 244   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
 245   if (fVtxBins->GetEntries() > 0) return;
 246
 247   fVtxBins->SetOwner();
 248   for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
 249     Double_t low  = fVtxAxis.GetBinLowEdge(i);
 250     Double_t high = fVtxAxis.GetBinUpEdge(i);
 251     VtxBin*  bin  = new VtxBin(low, high, fEtaAxis, fNPhiBins);
 252     fVtxBins->AddAt(bin, i);
 253     bin->DefineOutput(l);
 254   }
 255 }
 256
 257 //____________________________________________________________________
 258 void
 259 AliCentralMCCorrectionsTask::UserCreateOutputObjects()
 260 {
 261   //
 262   // Create output objects
 263   //
 264   //
 265   fList = new TList;
 266   fList->SetOwner();
 267   fList->SetName(Form("%sSums", GetName()));
 268
 269   DefineBins(fList);
 270
 271   fHEvents = new TH1I(GetEventName(false,false),
 272                       "Number of all events",
 273                       fVtxAxis.GetNbins(),
 274                       fVtxAxis.GetXmin(),
 275                       fVtxAxis.GetXmax());
 276   fHEvents->SetXTitle("v_{z} [cm]");
 277   fHEvents->SetYTitle("# of events");
 278   fHEvents->SetFillColor(kBlue+1);
 279   fHEvents->SetFillStyle(3001);
 280   fHEvents->SetDirectory(0);
 281   fList->Add(fHEvents);
 282
 283   fHEventsTr = new TH1I(GetEventName(true, false),
 284                         "Number of triggered events",
 285                         fVtxAxis.GetNbins(),
 286                         fVtxAxis.GetXmin(),
 287                         fVtxAxis.GetXmax());
 288   fHEventsTr->SetXTitle("v_{z} [cm]");
 289   fHEventsTr->SetYTitle("# of events");
 290   fHEventsTr->SetFillColor(kRed+1);
 291   fHEventsTr->SetFillStyle(3001);
 292   fHEventsTr->SetDirectory(0);
 293   fList->Add(fHEventsTr);
 294
 295   fHEventsTrVtx = new TH1I(GetEventName(true,true),
 296                            "Number of events w/trigger and vertex",
 297                            fVtxAxis.GetNbins(),
 298                            fVtxAxis.GetXmin(),
 299                            fVtxAxis.GetXmax());
 300   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 301   fHEventsTrVtx->SetYTitle("# of events");
 302   fHEventsTrVtx->SetFillColor(kBlue+1);
 303   fHEventsTrVtx->SetFillStyle(3001);
 304   fHEventsTrVtx->SetDirectory(0);
 305   fList->Add(fHEventsTrVtx);
 306
 307   // Copy axis objects to output
 308   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
 309                           fVtxAxis.GetNbins(),
 310                           fVtxAxis.GetXmin(),
 311                           fVtxAxis.GetXmax());
 312   TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
 313                           fEtaAxis.GetNbins(),
 314                           fEtaAxis.GetXmin(),
 315                           fEtaAxis.GetXmax());
 316   fList->Add(vtxAxis);
 317   fList->Add(etaAxis);
 318
 319   AliInfo(Form("Initialising sub-routines: %p, %p",
 320                &fInspector, &fTrackDensity));
 321   fInspector.DefineOutput(fList);
 322   fInspector.Init(fVtxAxis);
 323   fTrackDensity.DefineOutput(fList);
 324
 325   PostData(1, fList);
 326 }
 327 //____________________________________________________________________
 328 void
 329 AliCentralMCCorrectionsTask::UserExec(Option_t*)
 330 {
 331   //
 332   // Process each event
 333   //
 334   // Parameters:
 335   //    option Not used
 336   //
 337
 338   // Get the input data - MC event
 339   AliMCEvent*  mcEvent = MCEvent();
 340   if (!mcEvent) {
 341     AliWarning("No MC event found");
 342     return;
 343   }
 344
 345   // Get the input data - ESD event
 346   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 347   if (!esd) {
 348     AliWarning("No ESD event found for input event");
 349     return;
 350   }
 351
 352   //--- Read run information -----------------------------------------
 353   if (fFirstEvent && esd->GetESDRun()) {
 354     fInspector.ReadRunDetails(esd);
 355
 356     AliInfo(Form("Initializing with parameters from the ESD:\n"
 357                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 358                  "         AliESDEvent::GetBeamType()     ->%s\n"
 359                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 360                  "         AliESDEvent::GetMagneticField()->%f\n"
 361                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 362                  esd->GetBeamEnergy(),
 363                  esd->GetBeamType(),
 364                  esd->GetCurrentL3(),
 365                  esd->GetMagneticField(),
 366                  esd->GetRunNumber()));
 367
 368     Print();
 369     fFirstEvent = false;
 370   }
 371
 372   // Some variables
 373   UInt_t   triggers; // Trigger bits
 374   Bool_t   lowFlux;  // Low flux flag
 375   UShort_t iVz;      // Vertex bin from ESD
 376   Double_t vZ;       // Z coordinate from ESD
 377   Double_t cent;     // Centrality
 378   UShort_t iVzMc;    // Vertex bin from MC
 379   Double_t vZMc;     // Z coordinate of IP vertex from MC
 380   Double_t b;        // Impact parameter
 381   Int_t    nPart;    // Number of participants
 382   Int_t    nBin;     // Number of binary collisions
 383   Double_t phiR;     // Reaction plane from MC
 384   UShort_t nClusters;// Number of SPD clusters
 385   // Process the data
 386   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
 387                                      cent, nClusters);
 388   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
 389                        b, nPart, nBin, phiR);
 390
 391   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
 392   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
 393
 394   // Fill the event count histograms
 395   if (isInel)           fHEventsTr->Fill(vZMc);
 396   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
 397   fHEvents->Fill(vZMc);
 398
 399   // Now find our vertex bin object
 400   VtxBin* bin = 0;
 401   if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
 402     bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
 403   if (!bin) {
 404     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
 405     return;
 406   }
 407
 408   // Now process our input data and store in own ESD object
 409   fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
 410   bin->fCounts->Fill(0.5);
 411
 412   PostData(1, fList);
 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   DefineBins(fList);
 433
 434   // Output list
 435   TList* output = new TList;
 436   output->SetOwner();
 437   output->SetName(Form("%sResults", GetName()));
 438
 439   // --- Fill correction object --------------------------------------
 440   AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
 441   corr->SetVertexAxis(fVtxAxis);
 442   // corr->SetEtaAxis(fEtaAxis);
 443
 444   TIter     next(fVtxBins);
 445   VtxBin*   bin = 0;
 446   UShort_t  iVz = 1;
 447   while ((bin = static_cast<VtxBin*>(next())))
 448     bin->Finish(fList, output, iVz++, corr);
 449
 450   output->Add(corr);
 451
 452   PostData(2, output);
 453 }
 454
 455 //____________________________________________________________________
 456 void
 457 AliCentralMCCorrectionsTask::Print(Option_t* option) const
 458 {
 459   std::cout << ClassName() << "\n"
 460             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
 461             << "  Vertex range:     [" << fVtxAxis.GetXmin()
 462             << "," << fVtxAxis.GetXmax() << "]\n"
 463             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
 464             << "  Eta range:        [" << fEtaAxis.GetXmin()
 465             << "," << fEtaAxis.GetXmax() << "]\n"
 466             << "  # of phi bins:    " << fNPhiBins
 467             << std::endl;
 468   gROOT->IncreaseDirLevel();
 469   fInspector.Print(option);
 470   fTrackDensity.Print(option);
 471   gROOT->DecreaseDirLevel();
 472 }
 473
 474 //====================================================================
 475 const char*
 476 AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
 477                                              Double_t high)
 478 {
 479   static TString buf;
 480   buf = Form("vtx%+05.1f_%+05.1f", low, high);
 481   buf.ReplaceAll("+", "p");
 482   buf.ReplaceAll("-", "m");
 483   buf.ReplaceAll(".", "d");
 484   return buf.Data();
 485 }
 486
 487
 488 //____________________________________________________________________
 489 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
 490   : fHits(0),
 491     fPrimary(0),
 492     fCounts(0)
 493 {
 494 }
 495 //____________________________________________________________________
 496 AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
 497                                             Double_t     high,
 498                                             const TAxis& axis,
 499                                             UShort_t     nPhi)
 500   : TNamed(BinName(low, high),
 501            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
 502     fHits(0),
 503     fPrimary(0),
 504     fCounts(0)
 505 {
 506   fPrimary = new TH2D("primary", "Primaries",
 507                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
 508                       nPhi, 0, 2*TMath::Pi());
 509   fPrimary->SetXTitle("#eta");
 510   fPrimary->SetYTitle("#varphi [radians]");
 511   fPrimary->Sumw2();
 512   fPrimary->SetDirectory(0);
 513
 514   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
 515   fHits->SetTitle("Hits");
 516   fHits->SetDirectory(0);
 517
 518   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
 519   fCounts->SetXTitle("Events");
 520   fCounts->SetYTitle("# of Events");
 521   fCounts->SetDirectory(0);
 522 }
 523
 524 //____________________________________________________________________
 525 AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
 526   : TNamed(o),
 527     fHits(0),
 528     fPrimary(0),
 529     fCounts(0)
 530 {
 531   if (o.fHits) {
 532     fHits = static_cast<TH2D*>(o.fHits->Clone());
 533     fHits->SetDirectory(0);
 534   }
 535   if (o.fPrimary) {
 536     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 537     fPrimary->SetDirectory(0);
 538   }
 539   if (o.fCounts) {
 540     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 541     fCounts->SetDirectory(0);
 542   }
 543 }
 544
 545 //____________________________________________________________________
 546 AliCentralMCCorrectionsTask::VtxBin&
 547 AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
 548 {
 549   if (&o == this) return *this;
 550   TNamed::operator=(o);
 551   fHits    = 0;
 552   fPrimary = 0;
 553   fCounts  = 0;
 554   if (o.fHits) {
 555     fHits = static_cast<TH2D*>(o.fHits->Clone());
 556     fHits->SetDirectory(0);
 557   }
 558   if (o.fPrimary) {
 559     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
 560     fPrimary->SetDirectory(0);
 561   }
 562   if (o.fCounts) {
 563     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
 564     fCounts->SetDirectory(0);
 565   }
 566   return *this;
 567 }
 568
 569 //____________________________________________________________________
 570 void
 571 AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
 572 {
 573   TList* d = new TList;
 574   d->SetName(GetName());
 575   d->SetOwner();
 576   l->Add(d);
 577
 578   d->Add(fHits);
 579   d->Add(fPrimary);
 580   d->Add(fCounts);
 581 }
 582 //____________________________________________________________________
 583 void
 584 AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input,
 585                                             TList* output,
 586                                             UShort_t iVz,
 587                                             AliCentralCorrSecondaryMap* map)
 588 {
 589   TList* out = new TList;
 590   out->SetName(GetName());
 591   out->SetOwner();
 592   output->Add(out);
 593
 594   TList* l = static_cast<TList*>(input->FindObject(GetName()));
 595   if (!l) {
 596     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
 597     return;
 598   }
 599
 600
 601   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
 602   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
 603   if (!hits || !prim) {
 604     AliError(Form("Missing histograms: %p, %p", hits, prim));
 605     return;
 606   }
 607
 608   TH2D* h = static_cast<TH2D*>(hits->Clone("bgCorr"));
 609   h->SetDirectory(0);
 610   h->Divide(prim);
 611
 612   map->SetCorrection(iVz, h);
 613
 614   out->Add(h);
 615 }
 616
 617 //
 618 // EOF
 619 //