PWGLF/FORWARD/analysis2/AliBaseMCCorrectionsTask.cxx

   1 //
   2 // Calculate the corrections in the base regions
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODBaseMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliBaseMCCorrectionsTask.h"
  15 #include "AliBaseMCTrackDensity.h"
  16 #include "AliCorrectionManagerBase.h"
  17 #include "AliLog.h"
  18 #include "AliESDEvent.h"
  19 #include "AliMCEvent.h"
  20 #include "AliAODForwardMult.h"
  21 #include <TH1.h>
  22 #include <TH2D.h>
  23 #include <TDirectory.h>
  24 #include <TTree.h>
  25 #include <TList.h>
  26 #include <TROOT.h>
  27 #include <iostream>
  28
  29 //====================================================================
  30 namespace {
  31   const char* GetEventName(Bool_t tr, Bool_t vtx)
  32   {
  33     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  34   }
  35 }
  36
  37 //====================================================================
  38 AliBaseMCCorrectionsTask::AliBaseMCCorrectionsTask()
  39   : AliBaseESDTask(),
  40     fInspector(),
  41     fVtxBins(0),
  42     fHEvents(0),
  43     fHEventsTr(0),
  44     fHEventsTrVtx(0),
  45     fVtxAxis(),
  46     fEtaAxis(),
  47     fUseESDVertex(false),
  48     fAfterEventSel(false)
  49 {
  50   //
  51   // Constructor
  52   //
  53   // Parameters:
  54   //    name Name of task
  55   //
  56 }
  57
  58 //____________________________________________________________________
  59 AliBaseMCCorrectionsTask::AliBaseMCCorrectionsTask(const char* name,
  60                                                    AliCorrectionManagerBase* m)
  61   : AliBaseESDTask(name, "", m),
  62     fInspector("eventInspector"),
  63     fVtxBins(0),
  64     fHEvents(0),
  65     fHEventsTr(0),
  66     fHEventsTrVtx(0),
  67     fVtxAxis(10,-10,10),
  68     fEtaAxis(200,-4,6),
  69     fUseESDVertex(false),
  70     fAfterEventSel(false)
  71 {
  72   //
  73   // Constructor
  74   //
  75   // Parameters:
  76   //    name Name of task
  77   fBranchNames =
  78     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  79     "AliESDFMD.,SPDVertex.,PrimaryVertex.";
  80 }
  81
  82
  83 //____________________________________________________________________
  84 void
  85 AliBaseMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
  86                                            Double_t max)
  87 {
  88   //
  89   // Set the vertex axis to use
  90   //
  91   // Parameters:
  92   //    nBins Number of bins
  93   //    vzMin Least @f$z@f$ coordinate of interation point
  94   //    vzMax Largest @f$z@f$ coordinate of interation point
  95   //
  96   if (max < min) {
  97     Double_t tmp = min;
  98     min          = max;
  99     max          = tmp;
 100   }
 101   /*
 102     if (min < -10)
 103     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
 104     if (max > +10)
 105     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
 106   */
 107   fVtxAxis.Set(nBin, min, max);
 108 }
 109 //____________________________________________________________________
 110 void
 111 AliBaseMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
 112 {
 113   //
 114   // Set the vertex axis to use
 115   //
 116   // Parameters:
 117   //    axis Axis
 118   //
 119   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 120 }
 121
 122 //____________________________________________________________________
 123 void
 124 AliBaseMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
 125 {
 126   //
 127   // Set the eta axis to use
 128   //
 129   // Parameters:
 130   //    nBins Number of bins
 131   //    vzMin Least @f$\eta@f$
 132   //    vzMax Largest @f$\eta@f$
 133   //
 134   if (max < min) {
 135     Double_t tmp = min;
 136     min          = max;
 137     max          = tmp;
 138   }
 139   if (min < -4)
 140     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
 141   if (max > +6)
 142     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
 143   fEtaAxis.Set(nBin, min, max);
 144 }
 145 //____________________________________________________________________
 146 void
 147 AliBaseMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
 148 {
 149   //
 150   // Set the eta axis to use
 151   //
 152   // Parameters:
 153   //    axis Axis
 154   //
 155   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
 156 }
 157
 158 //____________________________________________________________________
 159 void
 160 AliBaseMCCorrectionsTask::DefineBins(TList* l)
 161 {
 162   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
 163   if (fVtxBins->GetEntries() > 0) return;
 164
 165   fVtxBins->SetOwner();
 166   for (Int_t  i    = 1; i <= fVtxAxis.GetNbins(); i++) {
 167     Double_t  low  = fVtxAxis.GetBinLowEdge(i);
 168     Double_t  high = fVtxAxis.GetBinUpEdge(i);
 169     VtxBin*   bin  = CreateVtxBin(low, high);
 170     fVtxBins->AddAt(bin, i);
 171     bin->CreateOutputObjects(l);
 172   }
 173 }
 174
 175 //____________________________________________________________________
 176 Bool_t
 177 AliBaseMCCorrectionsTask::Book()
 178 {
 179   //
 180   // Create output objects
 181   //
 182   //
 183   DefineBins(fList);
 184   fNeededCorrections = 0;
 185   fExtraCorrections  = 0;
 186
 187   fHEvents = new TH1I(GetEventName(false,false),
 188                       "Number of all events",
 189                       fVtxAxis.GetNbins(),
 190                       fVtxAxis.GetXmin(),
 191                       fVtxAxis.GetXmax());
 192   fHEvents->SetXTitle("v_{z} [cm]");
 193   fHEvents->SetYTitle("# of events");
 194   fHEvents->SetFillColor(kBlue+1);
 195   fHEvents->SetFillStyle(3001);
 196   fHEvents->SetDirectory(0);
 197   fList->Add(fHEvents);
 198
 199   fHEventsTr = new TH1I(GetEventName(true, false),
 200                         "Number of triggered events",
 201                         fVtxAxis.GetNbins(),
 202                         fVtxAxis.GetXmin(),
 203                         fVtxAxis.GetXmax());
 204   fHEventsTr->SetXTitle("v_{z} [cm]");
 205   fHEventsTr->SetYTitle("# of events");
 206   fHEventsTr->SetFillColor(kRed+1);
 207   fHEventsTr->SetFillStyle(3001);
 208   fHEventsTr->SetDirectory(0);
 209   fList->Add(fHEventsTr);
 210
 211   fHEventsTrVtx = new TH1I(GetEventName(true,true),
 212                            "Number of events w/trigger and vertex",
 213                            fVtxAxis.GetNbins(),
 214                            fVtxAxis.GetXmin(),
 215                            fVtxAxis.GetXmax());
 216   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 217   fHEventsTrVtx->SetYTitle("# of events");
 218   fHEventsTrVtx->SetFillColor(kBlue+1);
 219   fHEventsTrVtx->SetFillStyle(3001);
 220   fHEventsTrVtx->SetDirectory(0);
 221   fList->Add(fHEventsTrVtx);
 222
 223   // Copy axis objects to output
 224   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
 225                           fVtxAxis.GetNbins(),
 226                           fVtxAxis.GetXmin(),
 227                           fVtxAxis.GetXmax());
 228   TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
 229                           fEtaAxis.GetNbins(),
 230                           fEtaAxis.GetXmin(),
 231                           fEtaAxis.GetXmax());
 232   fList->Add(vtxAxis);
 233   fList->Add(etaAxis);
 234
 235   AliInfo(Form("Initialising sub-routines: %p, %p",
 236                &fInspector, &GetTrackDensity()));
 237   GetTrackDensity().CreateOutputObjects(fList);
 238
 239   return true;
 240 }
 241
 242 //____________________________________________________________________
 243 Bool_t
 244 AliBaseMCCorrectionsTask::Event(AliESDEvent& esd)
 245 {
 246   //
 247   // Process each event
 248   //
 249   // Parameters:
 250   //    option Not used
 251   //
 252
 253   // Get the input data - MC event
 254   AliMCEvent*  mcEvent = MCEvent();
 255   if (!mcEvent) {
 256     AliWarning("No MC event found");
 257     return false;
 258   }
 259
 260   // Some variables
 261   UInt_t   triggers  = 0;    // Trigger bits
 262   Bool_t   lowFlux   = true; // Low flux flag
 263   UShort_t iVz       = 0;    // Vertex bin from ESD
 264   TVector3 ip(1024,1024,1000);
 265   Double_t cent      = -1;   // Centrality
 266   UShort_t iVzMc     = 0;    // Vertex bin from MC
 267   Double_t vZMc      = 1000; // Z coordinate of IP vertex from MC
 268   Double_t b         = -1;   // Impact parameter
 269   Double_t cMC       = -1;   // Centrality estimate from b
 270   Int_t    nPart     = -1;   // Number of participants
 271   Int_t    nBin      = -1;   // Number of binary collisions
 272   Double_t phiR      = 100;  // Reaction plane from MC
 273   UShort_t nClusters = 0;    // Number of SPD clusters
 274   // Process the data
 275   UInt_t retESD = fInspector.Process(&esd, triggers, lowFlux, iVz, ip,
 276                                      cent, nClusters);
 277
 278   Bool_t isAccepted = true;
 279   if(fAfterEventSel) {
 280     if (retESD & AliFMDEventInspector::kNoEvent)    isAccepted = false;
 281     if (retESD & AliFMDEventInspector::kNoTriggers) isAccepted = false;
 282     if (retESD & AliFMDEventInspector::kNoVertex)   isAccepted = false;
 283     if (retESD & AliFMDEventInspector::kNoFMD)      isAccepted = false;
 284     if (!isAccepted) return false;
 285     // returns if there is not event , does not pass phys selection ,
 286     // has no veretx and lack of FMD data.
 287     // with good veretx outside z range it will continue
 288   }
 289
 290
 291   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
 292                        b, cMC, nPart, nBin, phiR);
 293
 294   fInspector.CompareResults(ip.Z(),    vZMc,
 295                             cent,      cMC,
 296                             b, nPart,  nBin);
 297
 298   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
 299   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
 300
 301   // Fill the event count histograms
 302   if (isInel)           fHEventsTr->Fill(vZMc);
 303   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
 304   fHEvents->Fill(vZMc);
 305
 306   // Now find our vertex bin object
 307   VtxBin*  bin      = 0;
 308   UShort_t usedZbin = iVzMc;
 309   if (fUseESDVertex) usedZbin = iVz;
 310
 311
 312   if (usedZbin > 0 && usedZbin <= fVtxAxis.GetNbins())
 313     bin = static_cast<VtxBin*>(fVtxBins->At(usedZbin));
 314   if (!bin) {
 315     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
 316     return false;
 317   }
 318
 319   return ProcessESD(esd, *mcEvent, *bin, vZMc);
 320 }
 321
 322 //____________________________________________________________________
 323 Bool_t
 324 AliBaseMCCorrectionsTask::Finalize()
 325 {
 326   //
 327   // End of job
 328   //
 329   // Parameters:
 330   //    option Not used
 331   //
 332   DefineBins(fList);
 333   CreateCorrections(fResults);
 334
 335   TIter     next(fVtxBins);
 336   VtxBin*   bin = 0;
 337   UShort_t  iVz = 1;
 338   while ((bin = static_cast<VtxBin*>(next())))
 339     FinalizeVtxBin(bin, iVz++);
 340
 341   return true;
 342 }
 343
 344 //____________________________________________________________________
 345 void
 346 AliBaseMCCorrectionsTask::Print(Option_t* option) const
 347 {
 348   AliBaseESDTask::Print(option);
 349   std::cout << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
 350             << "  Vertex range:     [" << fVtxAxis.GetXmin()
 351             << "," << fVtxAxis.GetXmax() << "]\n"
 352             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
 353             << "  Eta range:        [" << fEtaAxis.GetXmin()
 354             << "," << fEtaAxis.GetXmax() << "]"
 355             << std::endl;
 356 }
 357
 358 //====================================================================
 359 const char*
 360 AliBaseMCCorrectionsTask::VtxBin::BinName(Double_t low,
 361                                              Double_t high)
 362 {
 363   static TString buf;
 364   buf = Form("vtx%+05.1f_%+05.1f", low, high);
 365   buf.ReplaceAll("+", "p");
 366   buf.ReplaceAll("-", "m");
 367   buf.ReplaceAll(".", "d");
 368   return buf.Data();
 369 }
 370
 371
 372 //____________________________________________________________________
 373 AliBaseMCCorrectionsTask::VtxBin::VtxBin()
 374   : TNamed(),
 375     fPrimary(0),
 376     fCounts(0)
 377 {
 378 }
 379 //____________________________________________________________________
 380 AliBaseMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
 381                                          Double_t     high,
 382                                          const TAxis& axis,
 383                                          UShort_t     nPhi)
 384   : TNamed(BinName(low, high),
 385            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
 386     fPrimary(0),
 387     fCounts(0)
 388 {
 389   fPrimary = new TH2D("primary", "Primaries",
 390                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
 391                       nPhi, 0, 2*TMath::Pi());
 392   fPrimary->SetXTitle("#eta");
 393   fPrimary->SetYTitle("#varphi [radians]");
 394   fPrimary->Sumw2();
 395   fPrimary->SetDirectory(0);
 396
 397   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
 398   fCounts->SetXTitle("Events");
 399   fCounts->SetYTitle("# of Events");
 400   fCounts->SetDirectory(0);
 401 }
 402
 403
 404 //____________________________________________________________________
 405 TList*
 406 AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
 407 {
 408   TList* d = new TList;
 409   d->SetName(GetName());
 410   d->SetOwner();
 411   l->Add(d);
 412
 413   d->Add(fPrimary);
 414   d->Add(fCounts);
 415
 416   return d;
 417 }
 418
 419
 420 //
 421 // EOF
 422 //