PWGLF/FORWARD/analysis2/AliAODForwardMult.cxx

   1 //
   2 // Class that contains the forward multiplicity data per event
   3 //
   4 // This class contains a histogram of
   5 // @f[
   6 //   \frac{d^2N_{ch}}{d\eta d\phi}\quad,
   7 // @f]
   8 // as well as a trigger mask for each analysed event.
   9 //
  10 // The eta acceptance of the event is stored in the underflow bins of
  11 // the histogram.  So to build the final histogram, one needs to
  12 // correct for this acceptance (properly weighted by the events), and
  13 // the vertex efficiency.  This simply boils down to defining a 2D
  14 // histogram and summing the event histograms in that histogram.  One
  15 // should of course also do proper book-keeping of the accepted event.
  16 //
  17 #include "AliAODForwardMult.h"
  18 #include <TBrowser.h>
  19 #include <iostream>
  20 #include <TMath.h>
  21 #include <TObjString.h>
  22 #include <TObjArray.h>
  23 #include "AliLog.h"
  24 ClassImp(AliAODForwardMult)
  25 #ifdef DOXY_INPUT
  26 ; // For Emacs
  27 #endif
  28
  29 //____________________________________________________________________
  30 const Float_t AliAODForwardMult::fgkInvalidIpZ = 1e6;
  31
  32 //____________________________________________________________________
  33 AliAODForwardMult::AliAODForwardMult()
  34   : fIsMC(false),
  35     fHist(),
  36     fTriggers(0),
  37     fIpZ(fgkInvalidIpZ),
  38     fCentrality(-1),
  39     fNClusters(0)
  40 {
  41   //
  42   // Constructor
  43   //
  44 }
  45
  46 //____________________________________________________________________
  47 AliAODForwardMult::AliAODForwardMult(Bool_t isMC)
  48   : fIsMC(isMC),
  49     fHist("forwardMult", "d^{2}N_{ch}/d#etad#varphi in the forward regions",
  50           200, -4, 6, 20, 0, 2*TMath::Pi()),
  51     fTriggers(0),
  52     fIpZ(fgkInvalidIpZ),
  53     fCentrality(-1),
  54     fNClusters(0)
  55 {
  56   //
  57   // Constructor
  58   //
  59   // Parameters:
  60   //  isMC   If set to true this is for MC data (effects branch name)
  61   //
  62   fHist.SetXTitle("#eta");
  63   fHist.SetYTitle("#varphi [radians]");
  64   fHist.SetZTitle("#frac{d^{2}N_{ch}}{d#etad#varphi}");
  65   fHist.SetDirectory(0);
  66   fHist.Sumw2();
  67 }
  68
  69 //____________________________________________________________________
  70 void
  71 AliAODForwardMult::Init(const TAxis& etaAxis)
  72 {
  73   // Initialize the histogram with an eta axis
  74   //
  75   // Parameters:
  76   //   etaAxis       Eta axis to use
  77   //
  78   fHist.SetBins(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
  79                 20, 0, 2*TMath::Pi());
  80 }
  81
  82 //____________________________________________________________________
  83 void
  84 AliAODForwardMult::Clear(Option_t* option)
  85 {
  86   // Clear (or reset) internal values
  87   //
  88   // Parameters:
  89   //  option   Passed to TH1::Reset
  90   //
  91   fHist.Reset(option);
  92   fTriggers  = 0;
  93   fIpZ       = fgkInvalidIpZ;
  94   fNClusters = 0;
  95 }
  96 //____________________________________________________________________
  97 void
  98 AliAODForwardMult::SetSNN(UShort_t snn)
  99 {
 100   // set the center of mass energy per nucleon pair (GeV).
 101   // This is stored in bin (0,0) of the histogram
 102   //
 103   // Parameters:
 104   //   sNN   Center of mass energy per nuclean
 105   fHist.SetBinContent(0,0,snn);
 106 }
 107 //____________________________________________________________________
 108 void
 109 AliAODForwardMult::SetSystem(UShort_t sys)
 110 {
 111   // set the center of mass energy per nucleon pair (GeV).
 112   // This is stored in bin (N+1,0) of the histogram
 113   //
 114   // Parameters:
 115   //   sys   Collision system number
 116   fHist.SetBinContent(fHist.GetNbinsX()+1,0,sys);
 117 }
 118
 119 //____________________________________________________________________
 120 Bool_t
 121 AliAODForwardMult::HasIpZ() const
 122 {
 123   // Check if we have valid z coordinate of the interaction point
 124   //
 125   // Return:
 126   //   true if the z coordinate of the interaction point is valid
 127   //
 128   return TMath::Abs(fIpZ - fgkInvalidIpZ) > 1;
 129 }
 130
 131 //____________________________________________________________________
 132 UShort_t
 133 AliAODForwardMult::GetSNN() const
 134 {
 135   // set the center of mass energy per nucleon pair (GeV).
 136   // This is stored in bin (0,0) of the histogram
 137   //
 138   // Parameters:
 139   //   sNN   Center of mass energy per nuclean
 140   return UShort_t(fHist.GetBinContent(0,0));
 141 }
 142
 143 //____________________________________________________________________
 144 UShort_t
 145 AliAODForwardMult::GetSystem() const
 146 {
 147   // set the center of mass energy per nucleon pair (GeV).
 148   // This is stored in bin (N+1,0) of the histogram
 149   //
 150   // Parameters:
 151   //   sNN   Center of mass energy per nuclean
 152   return UShort_t(fHist.GetBinContent(fHist.GetNbinsX()+1,0));
 153 }
 154
 155 //____________________________________________________________________
 156 void
 157 AliAODForwardMult::Browse(TBrowser* b)
 158 {
 159   // Browse this object
 160   //
 161   // Parameters:
 162   //   b   Browser to use
 163   static TObjString ipz;
 164   static TObjString trg;
 165   static TObjString cnt;
 166   static TObjString ncl;
 167   ipz = Form("ip_z=%fcm", fIpZ);
 168   trg = GetTriggerString(fTriggers);
 169   cnt = Form("%+6.1f%%", fCentrality);
 170   ncl = Form("%d clusters", fNClusters);
 171   b->Add(&fHist);
 172   b->Add(&ipz);
 173   b->Add(&trg);
 174   b->Add(&cnt);
 175   b->Add(&ncl);
 176 }
 177
 178 namespace {
 179   void AppendAnd(TString& trg, const TString& what)
 180   {
 181     if (!trg.IsNull()) trg.Append(" & ");
 182     trg.Append(what);
 183   }
 184 }
 185 //____________________________________________________________________
 186 const Char_t*
 187 AliAODForwardMult::GetTriggerString(UInt_t mask)
 188 {
 189   // Get a string that describes the triggers
 190   //
 191   // Parameters:
 192   //   mask  Bit pattern of triggers
 193   // Return:
 194   //   Character string representation of mask
 195   static TString trg;
 196   trg = "";
 197   if ((mask & kInel)        != 0x0) AppendAnd(trg, "MBOR");
 198   if ((mask & kInelGt0)     != 0x0) AppendAnd(trg, "INEL>0");
 199   if ((mask & kNSD)         != 0x0) AppendAnd(trg, "MBAND");
 200   if ((mask & kV0AND)       != 0x0) AppendAnd(trg, "V0AND");
 201   if ((mask & kA)           != 0x0) AppendAnd(trg, "A");
 202   if ((mask & kB)           != 0x0) AppendAnd(trg, "B");
 203   if ((mask & kC)           != 0x0) AppendAnd(trg, "C");
 204   if ((mask & kE)           != 0x0) AppendAnd(trg, "E");
 205   if ((mask & kMCNSD)       != 0x0) AppendAnd(trg, "MCNSD");
 206   if ((mask & kNClusterGt0) != 0x0) AppendAnd(trg, "NCluster>0");
 207   if ((mask & kSatellite)   != 0x0) AppendAnd(trg, "Satellite");
 208   if ((mask & kOffline)     != 0x0) AppendAnd(trg, "Offline");
 209   return trg.Data();
 210 }
 211
 212 //____________________________________________________________________
 213 TH1I*
 214 AliAODForwardMult::MakeTriggerHistogram(const char* name, Int_t mask)
 215 {
 216   //
 217   // Make a histogram to record triggers in.
 218   //
 219   // The bins defined by the trigger enumeration in this class.  One
 220   // can use this enumeration to retrieve the number of triggers for
 221   // each class.
 222   //
 223   // Parameters:
 224   //    name Name of the histogram
 225   //
 226   // Return:
 227   //    Newly allocated histogram
 228   //
 229   TString sel("");
 230   TString andSel("");
 231   if (mask > 0) {
 232     sel    = GetTriggerString(mask);
 233     andSel = GetTriggerString(mask & ~kB);
 234     andSel.Prepend(" & ");
 235   }
 236   TH1I* ret = new TH1I(name, "Triggers", kAccepted+1, -.5, kAccepted+.5);
 237   ret->SetYTitle("Events");
 238   ret->SetFillColor(kRed+1);
 239   ret->SetFillStyle(3001);
 240   ret->GetXaxis()->SetBinLabel(kBinAll,         "All events");
 241   ret->GetXaxis()->SetBinLabel(kBinB,           Form("B (Coll.)%s",
 242                                                      andSel.Data()));
 243   ret->GetXaxis()->SetBinLabel(kBinA,           Form("A%s", andSel.Data()));
 244   ret->GetXaxis()->SetBinLabel(kBinC,           Form("C%s", andSel.Data()));
 245   ret->GetXaxis()->SetBinLabel(kBinE,           Form("E%s", andSel.Data()));
 246   ret->GetXaxis()->SetBinLabel(kBinInel,        "Coll. & MBOR");
 247   ret->GetXaxis()->SetBinLabel(kBinInelGt0,     "Coll. & MBOR&&nTracklet>0");
 248   ret->GetXaxis()->SetBinLabel(kBinNSD,         "Coll. & V0AND||FASTOR>5");
 249   ret->GetXaxis()->SetBinLabel(kBinV0AND,       "Coll. & V0AND");
 250   ret->GetXaxis()->SetBinLabel(kBinMCNSD,       "NSD (MC truth)");
 251   ret->GetXaxis()->SetBinLabel(kBinSatellite,   "Satellite");
 252   ret->GetXaxis()->SetBinLabel(kBinPileUp,      "w/Pileup");
 253   ret->GetXaxis()->SetBinLabel(kBinOffline,     "w/Offline");
 254   ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1");
 255   ret->GetXaxis()->SetBinLabel(kWithVertex,     "w/Vertex");
 256   ret->GetXaxis()->SetBinLabel(kWithTrigger,    Form("w/Selected trigger (%s)",
 257                                                      sel.Data()));
 258   ret->GetXaxis()->SetBinLabel(kAccepted,       "Accepted by cut");
 259   ret->GetXaxis()->SetNdivisions(kAccepted, false);
 260   ret->SetStats(0);
 261
 262   return ret;
 263 }
 264
 265 //____________________________________________________________________
 266 TH1I*
 267 AliAODForwardMult::MakeStatusHistogram(const char* name)
 268 {
 269   //
 270   // Make a histogram to record status in.
 271   //
 272   // The bins defined by the status enumeration in this class.
 273   //
 274   // Parameters:
 275   //    name Name of the histogram
 276   //
 277   // Return:
 278   //    Newly allocated histogram
 279   //
 280   Int_t nBins = kOutlierEvent;
 281   TH1I* ret = new TH1I(name, "Event selection status", nBins+1, -.5, nBins+.5);
 282   ret->SetYTitle("Events");
 283   ret->SetFillColor(kBlue+1);
 284   ret->SetFillStyle(3001);
 285   ret->GetXaxis()->SetBinLabel(kGoodEvent+1,       "Good");
 286   ret->GetXaxis()->SetBinLabel(kWrongCentrality+1, "Out-of-range centrality");
 287   ret->GetXaxis()->SetBinLabel(kWrongTrigger+1,    "Wrong trigger");
 288   ret->GetXaxis()->SetBinLabel(kIsPileup+1,        "Pile-up");
 289   ret->GetXaxis()->SetBinLabel(kNoVertex+1,        "No IP_{z}");
 290   ret->GetXaxis()->SetBinLabel(kWrongVertex+1,     "Out-or-range IP_{z}");
 291   ret->GetXaxis()->SetBinLabel(kOutlierEvent+1,    "SPD Outlier");
 292   ret->GetXaxis()->SetNdivisions(nBins, false);
 293   ret->SetStats(0);
 294   return ret;
 295 }
 296
 297 //____________________________________________________________________
 298 UInt_t
 299 AliAODForwardMult::MakeTriggerMask(const char* what)
 300 {
 301   UShort_t    trgMask = 0;
 302   TString     trgs(what);
 303   trgs.ToUpper();
 304   TObjArray*  parts = trgs.Tokenize("&");
 305   TObjString* trg;
 306   TIter       next(parts);
 307   while ((trg = static_cast<TObjString*>(next()))) {
 308     TString s(trg->GetString());
 309     s.Strip(TString::kBoth, ' ');
 310     s.ToUpper();
 311     if      (s.IsNull()) continue;
 312     if      (s.CompareTo("INEL")       == 0) trgMask |= kInel;
 313     else if (s.CompareTo("MBOR")       == 0) trgMask |= kInel;
 314     else if (s.CompareTo("INEL>0")     == 0) trgMask |= kInelGt0;
 315     else if (s.CompareTo("INELGT0")    == 0) trgMask |= kInelGt0;
 316     else if (s.CompareTo("MBAND")      == 0) trgMask |= kNSD;
 317     else if (s.CompareTo("NSD")        == 0) trgMask |= kV0AND;
 318     else if (s.CompareTo("V0AND")      == 0) trgMask |= kV0AND;
 319     else if (s.CompareTo("MCNSD")      == 0) trgMask |= kMCNSD;
 320     else if (s.CompareTo("B")          == 0) trgMask |= kB;
 321     else if (s.CompareTo("A")          == 0) trgMask |= kA;
 322     else if (s.CompareTo("C")          == 0) trgMask |= kC;
 323     else if (s.CompareTo("SAT")        == 0) trgMask |= kSatellite;
 324     else if (s.CompareTo("E")          == 0) trgMask |= kE;
 325     else if (s.CompareTo("NCLUSTER>0") == 0) trgMask |= kNClusterGt0;
 326     else if (s.CompareTo("CENT")       == 0) trgMask |= kInel;
 327     else if (s.CompareTo("OFFLINE")    == 0) trgMask |= kOffline;
 328     // trgMask &= ~(kInel|kInelGt0|kNSD|kV0AND|kMCNSD);
 329     else
 330       AliWarningGeneral("MakeTriggerMask",
 331                         Form("Unknown trigger %s", s.Data()));
 332   }
 333   delete parts;
 334   return trgMask;
 335 }
 336
 337 //____________________________________________________________________
 338 Bool_t
 339 AliAODForwardMult::CheckEvent(Int_t    triggerMask,
 340                               Double_t vzMin, Double_t vzMax,
 341                               UShort_t cMin,  UShort_t cMax,
 342                               TH1*     hist,  TH1*     status,
 343                               Bool_t   removePileup) const
 344 {
 345   //
 346   // Check if event meets the passses requirements.
 347   //
 348   // It returns true if @e all of the following is true
 349   //
 350   // - The trigger is within the bit mask passed.
 351   // - The vertex is within the specified limits.
 352   // - The centrality is within the specified limits, or if lower
 353   //   limit is equal to or larger than the upper limit.
 354   //
 355   // If a histogram is passed in the last parameter, then that
 356   // histogram is filled with the trigger bits.
 357   //
 358   // Parameters:
 359   //    triggerMask  Trigger mask
 360   //    vzMin        Minimum @f$ v_z@f$ (in centimeters)
 361   //    vzMax        Maximum @f$ v_z@f$ (in centimeters)
 362   //    cMin         Minimum centrality (in percent)
 363   //    cMax         Maximum centrality (in percent)
 364   //    hist         Histogram to fill
 365   //
 366   // Return:
 367   //    @c true if the event meets the requirements
 368   //
 369   if (cMin < cMax && (cMin > fCentrality || cMax <= fCentrality)) {
 370     if (status) status->Fill(kWrongCentrality);
 371     return false;
 372   }
 373
 374   if (hist) {
 375     Int_t tmp = triggerMask & ~kB;
 376     hist->AddBinContent(kBinAll);
 377     if (IsTriggerBits(kB|tmp))          hist->AddBinContent(kBinB);
 378     if (IsTriggerBits(kA|tmp))          hist->AddBinContent(kBinA);
 379     if (IsTriggerBits(kC|tmp))          hist->AddBinContent(kBinC);
 380     if (IsTriggerBits(kE|tmp))          hist->AddBinContent(kBinE);
 381     if (IsTriggerBits(kB|kInel))        hist->AddBinContent(kBinInel);
 382     if (IsTriggerBits(kB|kInelGt0))     hist->AddBinContent(kBinInelGt0);
 383     if (IsTriggerBits(kB|kNSD))         hist->AddBinContent(kBinNSD);
 384     if (IsTriggerBits(kB|kV0AND))       hist->AddBinContent(kBinV0AND);
 385     if (IsTriggerBits(kPileUp))         hist->AddBinContent(kBinPileUp);
 386     if (IsTriggerBits(kMCNSD))          hist->AddBinContent(kBinMCNSD);
 387     if (IsTriggerBits(kOffline))        hist->AddBinContent(kBinOffline);
 388     if (IsTriggerBits(kNClusterGt0))    hist->AddBinContent(kBinNClusterGt0);
 389     if (IsTriggerBits(kSatellite))      hist->AddBinContent(kBinSatellite);
 390     if (IsTriggerBits(triggerMask) && !IsTriggerBits(kB|tmp))
 391       Warning("CheckEvent", "event: 0x%x, mask: 0x%x, tmp: 0x%x, tmp|b: 0x%x",
 392              fTriggers, triggerMask, tmp, tmp|kB);
 393   }
 394   // Check if we have an event of interest.
 395   Int_t mask = triggerMask; //|kB
 396   if (!IsTriggerBits(mask)) {
 397     if (status) status->Fill(kWrongTrigger);
 398     return false;
 399   }
 400   // Check if event is SPD outlier
 401   if (IsTriggerBits(kSPDOutlier)) {
 402     if (status) status->Fill(kOutlierEvent);
 403     return false;
 404   }
 405   // Check for pileup
 406   if (IsTriggerBits(kPileUp)) {
 407     if (status) status->Fill(kIsPileup);
 408     if (removePileup) return false;
 409   }
 410   if (hist) hist->AddBinContent(kWithTrigger);
 411
 412   // Check that we have a valid vertex
 413   if (vzMin < vzMax && !HasIpZ()) {
 414     if (status) status->Fill(kNoVertex);
 415     return false;
 416   }
 417   if (hist) hist->AddBinContent(kWithVertex);
 418
 419   // Check that vertex is within cuts
 420   if (vzMin < vzMax && !InRange(vzMin, vzMax)) {
 421     if (status) status->Fill(kWrongVertex);
 422     return false;
 423   }
 424   if (hist) hist->AddBinContent(kAccepted);
 425
 426   if (status) status->Fill(kGoodEvent);
 427   return true;
 428 }
 429
 430 //____________________________________________________________________
 431 void
 432 AliAODForwardMult::Print(Option_t* option) const
 433 {
 434   // Print this object
 435   //
 436   // Parameters:
 437   //  option   Passed to TH1::Print
 438   fHist.Print(option);
 439   UShort_t sys = GetSystem();
 440   TString  str = "unknown";
 441   switch (sys) {
 442   case 1:  str = "pp"; break;
 443   case 2:  str = "PbPb"; break;
 444   case 3:  str = "pPb" ; break;
 445   }
 446   std::cout << "Ipz:         " << fIpZ << "cm " << (HasIpZ() ? "" : "in")
 447             << "valid\n"
 448             << "Triggers:    " << GetTriggerString(fTriggers)  << "\n"
 449             << "sNN:         " << GetSNN() << "GeV\n"
 450             << "System:      " << str << "\n"
 451             << "Centrality:  " << fCentrality << "%"
 452             << std::endl;
 453 }
 454
 455 //____________________________________________________________________
 456 //
 457 // EOF
 458 //