PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx

   1 //
   2 // This class inspects the event
   3 //
   4 // Input:
   5 //   - AliESDFMD object possibly corrected for sharing
   6 //
   7 // Output:
   8 //   - A histogram of v_z of events with triggers.
   9 //   - A histogram of v_z of events with vertex and triggers
  10 //   - A histogram of trigger counters
  11 //
  12 // Note, that these are added to the master output list
  13 //
  14 // Corrections used:
  15 //   - None
  16 //
  17 #include "AliFMDEventInspector.h"
  18 #include "AliProdInfo.h"
  19 #include "AliLog.h"
  20 #include "AliESDEvent.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliMCEventHandler.h"
  24 #include "AliInputEventHandler.h"
  25 #include "AliTriggerAnalysis.h"
  26 #include "AliPhysicsSelection.h"
  27 #include "AliOADBPhysicsSelection.h"
  28 #include "AliAODForwardMult.h"
  29 #include "AliForwardUtil.h"
  30 #include "AliCentrality.h"
  31 #include <TH1.h>
  32 #include <TList.h>
  33 #include <TDirectory.h>
  34 #include <TROOT.h>
  35 #include <TParameter.h>
  36 #include <TMatrixDSym.h>
  37 #include <TPRegexp.h>
  38 #include <iostream>
  39 #include <iomanip>
  40 #include "AliMCEvent.h"
  41 #include "AliHeader.h"
  42 #include "AliGenEventHeader.h"
  43 #include "AliCollisionGeometry.h"
  44 #include "AliVVZERO.h"
  45
  46 //====================================================================
  47 const char* AliFMDEventInspector::fgkFolderName = "fmdEventInspector";
  48
  49 //____________________________________________________________________
  50 AliFMDEventInspector::AliFMDEventInspector()
  51   : TNamed(),
  52     fHEventsTr(0),
  53     fHEventsTrVtx(0),
  54     fHEventsAccepted(0),
  55     fHEventsAcceptedXY(0),
  56     fHTriggers(0),
  57     fHTriggerCorr(0),
  58     fHType(0),
  59     fHWords(0),
  60     fHCent(0),
  61     fHCentVsQual(0),
  62     fHStatus(0),
  63     fHVtxStatus(0),
  64     fHTrgStatus(0),
  65     fLowFluxCut(1000),
  66     fMaxVzErr(0.2),
  67     fList(0),
  68     fEnergy(0),
  69     fField(999),
  70     fCollisionSystem(kUnknown),
  71     fDebug(0),
  72     fCentAxis(0),
  73     fVtxAxis(10,-10,10),
  74     fUseFirstPhysicsVertex(false),
  75     fUseV0AND(false),
  76     fMinPileupContrib(3),
  77     fMinPileupDistance(0.8),
  78     fUseDisplacedVertices(false),
  79     fDisplacedVertex(),
  80     fCollWords(),
  81     fBgWords(),
  82     fCentMethod("V0M"),
  83     fMinCent(-1.0),
  84     fMaxCent(-1.0),
  85     fUsepA2012Vertex(false),
  86     fRunNumber(0),
  87     fMC(false),
  88   fProdYear(-1),
  89   fProdLetter('?'),
  90   fProdPass(-1),
  91   fProdSVN(-1),
  92   fProdMC(false)
  93 {
  94   //
  95   // Constructor
  96   //
  97   DGUARD(fDebug,1,"Default CTOR of AliFMDEventInspector");
  98 }
  99
 100 //____________________________________________________________________
 101 AliFMDEventInspector::AliFMDEventInspector(const char* name)
 102   : TNamed(fgkFolderName, name),
 103     fHEventsTr(0),
 104     fHEventsTrVtx(0),
 105     fHEventsAccepted(0),
 106     fHEventsAcceptedXY(0),
 107     fHTriggers(0),
 108     fHTriggerCorr(0),
 109     fHType(0),
 110     fHWords(0),
 111     fHCent(0),
 112     fHCentVsQual(0),
 113     fHStatus(0),
 114     fHVtxStatus(0),
 115     fHTrgStatus(0),
 116     fLowFluxCut(1000),
 117     fMaxVzErr(0.2),
 118     fList(0),
 119     fEnergy(0),
 120     fField(999),
 121     fCollisionSystem(kUnknown),
 122     fDebug(0),
 123     fCentAxis(0),
 124     fVtxAxis(10,-10,10),
 125     fUseFirstPhysicsVertex(false),
 126     fUseV0AND(false),
 127     fMinPileupContrib(3),
 128     fMinPileupDistance(0.8),
 129     fUseDisplacedVertices(false),
 130   fDisplacedVertex(),
 131   fCollWords(),
 132   fBgWords(),
 133   fCentMethod("V0M"),
 134   fMinCent(-1.0),
 135   fMaxCent(-1.0),
 136   fUsepA2012Vertex(false),
 137   fRunNumber(0),
 138   fMC(false),
 139   fProdYear(-1),
 140   fProdLetter('?'),
 141   fProdPass(-1),
 142   fProdSVN(-1),
 143   fProdMC(false)
 144 {
 145   //
 146   // Constructor
 147   //
 148   // Parameters:
 149   //   name Name of object
 150   //
 151   DGUARD(fDebug,1,"Named CTOR of AliFMDEventInspector: %s", name);
 152 }
 153
 154 //____________________________________________________________________
 155 AliFMDEventInspector::~AliFMDEventInspector()
 156 {
 157   //
 158   // Destructor
 159   //
 160   DGUARD(fDebug,1,"DTOR of AliFMDEventInspector");
 161   // if (fList)         delete fList;
 162 }
 163
 164 //____________________________________________________________________
 165 void
 166 AliFMDEventInspector::SetCentralityMethod(ECentMethod m)
 167 {
 168   switch (m) {
 169   case kV0Multiplicity: fCentMethod = "VOM"; break; // VZERO multiplicity
 170   case kV0Amplitude:    fCentMethod = "V0A"; break; // VZERO amplitude
 171   case kV0Charge:       fCentMethod = "V0C"; break; // VZERO charge
 172   case kFMDRough:       fCentMethod = "FMD"; break; // FMD scaled energy l
 173   case kNTracks:        fCentMethod = "TRK"; break; // Number of tracks
 174   case kLTracks:        fCentMethod = "TKL"; break; // Number of tracks
 175   case kCL0:            fCentMethod = "CL0"; break; //
 176   case kCL1:            fCentMethod = "CL1"; break; //
 177   case kCND:            fCentMethod = "CND"; break; //
 178   case kNParticles:     fCentMethod = "NPA"; break; // Neutral particles
 179   case kNeutrons:       fCentMethod = "ZNA"; break; // ZDC neutron amplitu
 180   case kV0vsFMD:        fCentMethod = "V0MvsFMD"; break; // VZERO versus FMD
 181   case kV0vsNTracks:    fCentMethod = "TKLvsVOM"; break; // Tracks versus VZERO
 182   case kZEMvsZDC:       fCentMethod = "ZEMvsZDC"; break; // ZDC
 183   default:              fCentMethod = "V0M"; break;
 184   }
 185 }
 186
 187 //____________________________________________________________________
 188 void
 189 AliFMDEventInspector::SetMinCentrality(Double_t minCent)
 190 {
 191   AliWarning("\n"
 192              "*******************************************************\n"
 193              "* Setting centrality cuts in this stage is deprecated *\n"
 194              "*******************************************************");
 195   fMinCent = minCent;
 196 }
 197 //____________________________________________________________________
 198 void
 199 AliFMDEventInspector::SetMaxCentrality(Double_t maxCent)
 200 {
 201   AliWarning("\n"
 202              "*******************************************************\n"
 203              "* Setting centrality cuts in this stage is deprecated *\n"
 204              "*******************************************************");
 205   fMaxCent = maxCent;
 206 }
 207
 208 //____________________________________________________________________
 209 Bool_t
 210 AliFMDEventInspector::FetchHistograms(const TList* d,
 211                                       TH1I*& hEventsTr,
 212                                       TH1I*& hEventsTrVtx,
 213                                       TH1I*& hEventsAcc,
 214                                       TH1I*& hTriggers) const
 215 {
 216   //
 217   // Fetch our histograms from the passed list
 218   //
 219   // Parameters:
 220   //   d             Input
 221   //   hEventsTr     On return, pointer to histogram, or null
 222   //   hEventsTrVtx  On return, pointer to histogram, or null
 223   //   hTriggers     On return, pointer to histogram, or null
 224   //
 225   // Return:
 226   //    true on success, false otherwise
 227   //
 228   DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
 229   hEventsTr    = 0;
 230   hEventsTrVtx = 0;
 231   hEventsAcc   = 0;
 232   hTriggers    = 0;
 233   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
 234   if (!dd) return kFALSE;
 235
 236   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
 237   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
 238   hEventsAcc   = dynamic_cast<TH1I*>(dd->FindObject("nEventsAccepted"));
 239   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
 240
 241   if (!hEventsTr    ||
 242       !hEventsTrVtx ||
 243       !hEventsAcc   ||
 244       !hTriggers) return kFALSE;
 245   return kTRUE;
 246 }
 247 //____________________________________________________________________
 248 void
 249 AliFMDEventInspector::CacheConfiguredTriggerClasses(TList& cache,
 250                                                     const TList* classes,
 251                                                     AliOADBPhysicsSelection* o)
 252 {
 253   TIter nextClass(classes);
 254   TObjString* trigClass = 0;
 255   // Loop over all trigger classes.  Trigger classes have the format
 256   //
 257   //   class          := positive_words SPACE(s) negative_words
 258   //   positive_words :=
 259   //                  |  '+' words
 260   //   negative_words :=
 261   //                  |  '-' words
 262   //   words          := word
 263   //                  |  word ',' words
 264   //
 265   while ((trigClass = static_cast<TObjString*>(nextClass()))) {
 266     // Tokenize on space to get positive and negative parts
 267     TString     side   = o->GetBeamSide(trigClass->String());
 268     TObjArray*  parts  = trigClass->String().Tokenize(" ");
 269     TObjString* part   = 0;
 270     TIter       nextPart(parts);
 271     while ((part = static_cast<TObjString*>(nextPart()))) {
 272       // We only care about the positive ones
 273       if (part->GetName()[0] != '+') continue;
 274       part->String().Remove(0,1);
 275
 276       // Tokenize on a comma to get the words
 277       TObjArray*  words = part->String().Tokenize(",");
 278       TObjString* word  = 0;
 279       TIter       nextWord(words);
 280       while ((word = static_cast<TObjString*>(nextWord()))) {
 281         TNamed* store = new TNamed(word->String(), side);
 282         cache.Add(store);
 283         DMSG(fDebug,3,"Caching %s trigger word %s",
 284              store->GetTitle(), store->GetName());
 285       } // while (word)
 286       delete words;
 287     }
 288     delete parts;
 289   }
 290 }
 291
 292 //____________________________________________________________________
 293 void
 294 AliFMDEventInspector::SetupForData(const TAxis& vtxAxis)
 295 {
 296   //
 297   // Initialize the object - this is called on the first seen event.
 298   //
 299   // Parameters:
 300   //   vtxAxis Vertex axis in use
 301   //
 302   DGUARD(fDebug,1,"Initialize in AliFMDEventInspector");
 303
 304   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 305
 306   // Get the input handler - should always be there
 307   AliInputEventHandler* ih =
 308     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 309   if (!ih) {
 310     AliWarning("No input handler");
 311     return;
 312   }
 313   // Get the physics selection - should always be there
 314   AliPhysicsSelection* ps =
 315     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
 316   if (!ps) {
 317     AliWarning("No physics selection");
 318     return;
 319   }
 320   // Get the configured triggers
 321   AliOADBPhysicsSelection* oadb =
 322     const_cast<AliOADBPhysicsSelection*>(ps->GetOADBPhysicsSelection());
 323   if (!oadb) {
 324     AliWarning("No OADB physics selection object");
 325     return;
 326   }
 327   // Get the configured trigger words from the physics selection
 328   const TList* collTriggClasses = ps->GetCollisionTriggerClasses();
 329   const TList* bgTriggClasses   = ps->GetBGTriggerClasses();
 330   if (!collTriggClasses) {
 331     AliWarning("No configured collision trigger classes");
 332     return;
 333   }
 334   if (!bgTriggClasses) {
 335     AliWarning("No configured background trigger classes");
 336     return;
 337   }
 338   CacheConfiguredTriggerClasses(fCollWords, collTriggClasses, oadb);
 339   CacheConfiguredTriggerClasses(fBgWords,   bgTriggClasses,   oadb);
 340   // fCollWords.ls();
 341   // fBgWords.ls();
 342
 343
 344   TArrayD limits;
 345   if ((fMinCent < 0 && fMaxCent < 0) || fMaxCent <= fMinCent) {
 346     // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
 347     // ----- 92 number --------- ---- 1 ---
 348     limits.Set(93);
 349     for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
 350     limits[92] = 100.5;
 351   }
 352   else {
 353     Int_t n = fMaxCent-fMinCent+2;
 354     limits.Set(n);
 355     for (Int_t i = 0; i < n; i++) {
 356       limits[i] = fMinCent + i - .5;
 357     }
 358   }
 359
 360   fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
 361
 362   fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
 363   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
 364                         4*vtxAxis.GetNbins(),
 365                         2*vtxAxis.GetXmin(),
 366                         2*vtxAxis.GetXmax());
 367   fHEventsTr->SetXTitle("v_{z} [cm]");
 368   fHEventsTr->SetYTitle("# of events");
 369   fHEventsTr->SetFillColor(kRed+1);
 370   fHEventsTr->SetFillStyle(3001);
 371   fHEventsTr->SetDirectory(0);
 372   // fHEventsTr->Sumw2();
 373   fList->Add(fHEventsTr);
 374
 375   fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
 376   fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex");
 377   fHEventsTrVtx->SetFillColor(kBlue+1);
 378   fHEventsTrVtx->SetDirectory(0);
 379   // fHEventsTrVtx->Sumw2();
 380   fList->Add(fHEventsTrVtx);
 381
 382   fHEventsAccepted = new TH1I("nEventsAccepted",
 383                               "Number of events  w/trigger and vertex in range",
 384                               2*vtxAxis.GetNbins(),
 385                               2*vtxAxis.GetXmin(),
 386                               2*vtxAxis.GetXmax());
 387   fHEventsAccepted->SetXTitle("v_{z} [cm]");
 388   fHEventsAccepted->SetYTitle("# of events");
 389   fHEventsAccepted->SetFillColor(kGreen+1);
 390   fHEventsAccepted->SetFillStyle(3001);
 391   fHEventsAccepted->SetDirectory(0);
 392   // fHEventsAccepted->Sumw2();
 393   fList->Add(fHEventsAccepted);
 394
 395   fHEventsAcceptedXY = new TH2D("nEventsAcceptedXY",
 396                               "XY vertex w/trigger and Z vertex in range",
 397                                 1000,-1,1,1000,-1,1);
 398
 399   fHEventsAcceptedXY->SetXTitle("v_{x} [cm]");
 400   fHEventsAcceptedXY->SetYTitle("v_{y} [cm]");
 401   fHEventsAcceptedXY->SetDirectory(0);
 402   // fHEventsAccepted->Sumw2();
 403   fList->Add(fHEventsAcceptedXY);
 404
 405
 406   fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
 407   fHTriggers->SetFillColor(kRed+1);
 408   fHTriggers->SetFillStyle(3001);
 409   fHTriggers->SetStats(0);
 410   fHTriggers->SetDirectory(0);
 411
 412   fHTriggerCorr = new TH2I("triggerCorr", "Trigger correlation",
 413                            kOffline+1, 0, kOffline+1,
 414                            kOffline+1, 0, kOffline+1);
 415   fHTriggerCorr->SetStats(0);
 416   fHTriggerCorr->SetDirectory(0);
 417   fHTriggerCorr->SetXTitle("Requirement");
 418   fHTriggerCorr->SetYTitle("Companion");
 419
 420   Int_t binNum[] = { kInel   +1,
 421                      kInelGt0+1,
 422                      kNSD    +1,
 423                      kV0AND  +1,
 424                      kEmpty  +1,
 425                      kA      +1,
 426                      kB      +1,
 427                      kC      +1,
 428                      kE      +1,
 429                      kPileUp +1,
 430                      kMCNSD  +1,
 431                      kSatellite+1,
 432                      kOffline+1 };
 433   const char* binLbl[] = { "INEL",
 434                            "INEL>0",
 435                            "NSD",
 436                            "VOAND",
 437                            "Empty",
 438                            "A",
 439                            "B",
 440                            "C",
 441                            "E",
 442                            "Pileup",
 443                            "NSD_{MC}",
 444                            "Satellite",
 445                            "Offline" };
 446   for (Int_t i = 0; i < kOffline+1; i++) {
 447     fHTriggers->GetXaxis()->SetBinLabel(binNum[i], binLbl[i]);
 448     fHTriggerCorr->GetXaxis()->SetBinLabel(binNum[i], binLbl[i]);
 449     fHTriggerCorr->GetYaxis()->SetBinLabel(binNum[i], binLbl[i]);
 450   }
 451   fList->Add(fHTriggers);
 452   fList->Add(fHTriggerCorr);
 453
 454
 455   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)",
 456                                  fLowFluxCut), 2, -.5, 1.5);
 457   fHType->SetFillColor(kRed+1);
 458   fHType->SetFillStyle(3001);
 459   fHType->SetStats(0);
 460   fHType->SetDirectory(0);
 461   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
 462   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
 463   fList->Add(fHType);
 464
 465 #if 0
 466   // This histogram disabled as it causes problems in the merge
 467   fHWords = new TH1I("words", "Trigger words seen", 1, 0, 0);
 468   fHWords->SetFillColor(kBlue+1);
 469   fHWords->SetFillStyle(3001);
 470   fHWords->SetStats(0);
 471   fHWords->SetDirectory(0);
 472   fHWords->SetBit(TH1::kCanRebin);
 473   fList->Add(fHWords);
 474 #endif
 475
 476   fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
 477   fHCent->SetFillColor(kBlue+1);
 478   fHCent->SetFillStyle(3001);
 479   fHCent->SetStats(0);
 480   fHCent->SetDirectory(0);
 481   fHCent->SetXTitle("Centrality [%]");
 482   fHCent->SetYTitle("Events");
 483   fList->Add(fHCent);
 484
 485   fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality",
 486                           5, 0, 5, limits.GetSize()-1, limits.GetArray());
 487   fHCentVsQual->SetXTitle("Quality");
 488   fHCentVsQual->SetYTitle("Centrality [%]");
 489   fHCentVsQual->SetZTitle("Events");
 490   fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
 491   fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
 492   fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
 493   fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
 494   fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
 495   fHCentVsQual->SetDirectory(0);
 496   fList->Add(fHCentVsQual);
 497
 498   fHStatus = new TH1I("status", "Status", 7, 1, 8);
 499   fHStatus->SetFillColor(kBlue+1);
 500   fHStatus->SetFillStyle(3001);
 501   fHStatus->SetStats(0);
 502   fHStatus->SetDirectory(0);
 503   TAxis* xAxis = fHStatus->GetXaxis();
 504   xAxis->SetBinLabel(1, "OK");
 505   xAxis->SetBinLabel(2, "No event");
 506   xAxis->SetBinLabel(3, "No triggers");
 507   xAxis->SetBinLabel(4, "No SPD");
 508   xAxis->SetBinLabel(5, "No FMD");
 509   xAxis->SetBinLabel(6, "No vertex");
 510   xAxis->SetBinLabel(7, "Bad vertex");
 511   fList->Add(fHStatus);
 512
 513   fHVtxStatus = new TH1I("vtxStatus","Vertex Status",
 514                          kNotVtxZ,kVtxOK,kNotVtxZ+1);
 515   fHVtxStatus->SetFillColor(kGreen+1);
 516   fHVtxStatus->SetFillStyle(3001);
 517   fHVtxStatus->SetStats(0);
 518   fHVtxStatus->SetDirectory(0);
 519   xAxis = fHVtxStatus->GetXaxis();
 520   xAxis->SetBinLabel(kVtxOK,      "OK");
 521   xAxis->SetBinLabel(kNoVtx,      "None/bad status");
 522   xAxis->SetBinLabel(kNoSPDVtx,   "No SPD/bad status");
 523   xAxis->SetBinLabel(kFewContrib, "N_{contrib} <= 0");
 524   xAxis->SetBinLabel(kUncertain,  Form("#delta z > %4.2f", fMaxVzErr));
 525   xAxis->SetBinLabel(kNotVtxZ,    "Not Z vertexer");
 526   fList->Add(fHVtxStatus);
 527
 528   fHTrgStatus = new TH1I("trgStatus", "Trigger Status",
 529                          kOther, kNoTrgWords, kOther+1);
 530   fHTrgStatus->SetFillColor(kMagenta+1);
 531   fHTrgStatus->SetFillStyle(3001);
 532   fHTrgStatus->SetStats(0);
 533   fHTrgStatus->SetDirectory(0);
 534   xAxis = fHTrgStatus->GetXaxis();
 535   xAxis->SetBinLabel(kNoTrgWords,       "No words");
 536   xAxis->SetBinLabel(kPP2760Fast,       "FAST in pp@#sqrt{s}=2.76TeV");
 537   xAxis->SetBinLabel(kMUON,             "Muon trigger");
 538   xAxis->SetBinLabel(kTriggered,        "Triggered");
 539   xAxis->SetBinLabel(kMinBias,          "CINT1 (V0A||V0C||FASTOR)");
 540   xAxis->SetBinLabel(kMinBiasNoSPD,     "CINT5 (V0A||V0C)");
 541   xAxis->SetBinLabel(kV0AndTrg,         "CINT7 (V0A&&V0C)");
 542   xAxis->SetBinLabel(kHighMult,         "N>>0");
 543   xAxis->SetBinLabel(kCentral,          "Central");
 544   xAxis->SetBinLabel(kSemiCentral,      "Semi-central");
 545   xAxis->SetBinLabel(kDiffractive,      "Diffractive");
 546   xAxis->SetBinLabel(kUser,             "User");
 547   xAxis->SetBinLabel(kOther,            "Other");
 548   fList->Add(fHTrgStatus);
 549
 550   if (fUseDisplacedVertices) fDisplacedVertex.SetupForData(fList, "", false);
 551 }
 552
 553 //____________________________________________________________________
 554 void
 555 AliFMDEventInspector::StoreInformation()
 556 {
 557   // Write TNamed objects to output list containing information about
 558   // the running conditions
 559   DGUARD(fDebug,2,"Store information from AliFMDEventInspector");
 560   if (!fList) return;
 561
 562
 563   fList->Add(AliForwardUtil::MakeParameter("sys", fCollisionSystem));
 564   fList->Add(AliForwardUtil::MakeParameter("sNN", fEnergy));
 565   fList->Add(AliForwardUtil::MakeParameter("field", fField));
 566   fList->Add(AliForwardUtil::MakeParameter("runNo", fRunNumber));
 567   fList->Add(AliForwardUtil::MakeParameter("lowFlux", fLowFluxCut));
 568   fList->Add(AliForwardUtil::MakeParameter("fpVtx",fUseFirstPhysicsVertex));
 569   fList->Add(AliForwardUtil::MakeParameter("v0and",fUseV0AND));
 570   fList->Add(AliForwardUtil::MakeParameter("nPileUp", fMinPileupContrib));
 571   fList->Add(AliForwardUtil::MakeParameter("dPileup", fMinPileupDistance));
 572   fList->Add(AliForwardUtil::MakeParameter("satellite", fUseDisplacedVertices));
 573   fList->Add(AliForwardUtil::MakeParameter("alirootRev",
 574                                            AliForwardUtil::AliROOTRevision()));
 575   fList->Add(AliForwardUtil::MakeParameter("alirootBranch",
 576                                            AliForwardUtil::AliROOTBranch()));
 577   fList->Add(AliForwardUtil::MakeParameter("mc", fMC));
 578
 579 }
 580
 581 //____________________________________________________________________
 582 void
 583 AliFMDEventInspector::StoreProduction()
 584 {
 585   if (!fList) return;
 586
 587   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 588   AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
 589   if (!inputHandler) {
 590     AliWarning("Got no input handler");
 591     return;
 592   }
 593   TList *uiList = inputHandler->GetUserInfo();
 594   if (!uiList) {
 595     AliWarning("Got no user list from input tree");
 596     return;
 597   }
 598
 599   AliProdInfo p(uiList);
 600   p.List();
 601   if (p.GetAlirootSvnVersion() <= 0) return;
 602
 603   // Make our output list
 604   TList* out = new TList;
 605   out->SetOwner(true);
 606   out->SetName("production");
 607   // out->SetTitle("ESD production details");
 608   fList->Add(out);
 609
 610   TString period            = p.GetLHCPeriod();
 611   // TString aliROOTVersion    = p.GetAlirootVersion();
 612   fProdSVN                  = p.GetAlirootSvnVersion();
 613   // TString rootVersion       = p.GetRootVersion();
 614   Int_t   rootSVN           = p.GetRootSvnVersion();
 615   fProdPass                 = p.GetRecoPass();
 616   fProdMC                   = p.IsMC();
 617
 618   TObjArray* pp = TPRegexp("LHC([0-9]+)([a-z]+)").MatchS(period);
 619   fProdYear     = static_cast<TObjString*>(pp->At(1))->String().Atoi();
 620   fProdLetter   = static_cast<TObjString*>(pp->At(2))->String()[0];
 621   pp->Delete();
 622
 623   out->Add(AliForwardUtil::MakeParameter("year", fProdYear));
 624   out->Add(AliForwardUtil::MakeParameter("letter", Int_t(fProdLetter)));
 625   out->Add(AliForwardUtil::MakeParameter("alirootSVN", fProdSVN));
 626   out->Add(AliForwardUtil::MakeParameter("rootSVN", rootSVN));
 627   out->Add(AliForwardUtil::MakeParameter("pass", fProdPass));
 628   out->Add(AliForwardUtil::MakeParameter("mc", fProdMC));
 629 }
 630
 631
 632
 633
 634 //____________________________________________________________________
 635 void
 636 AliFMDEventInspector::CreateOutputObjects(TList* dir)
 637 {
 638   //
 639   // Define the output histograms.  These are put in a sub list of the
 640   // passed list.   The histograms are merged before the parent task calls
 641   // AliAnalysisTaskSE::Terminate
 642   //
 643   //   dir Directory to add to
 644   //
 645   DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
 646   fList = new TList;
 647   fList->SetName(GetName());
 648   fList->SetOwner();
 649   dir->Add(fList);
 650 }
 651
 652 //____________________________________________________________________
 653 UInt_t
 654 AliFMDEventInspector::Process(const AliESDEvent* event,
 655                               UInt_t&            triggers,
 656                               Bool_t&            lowFlux,
 657                               UShort_t&          ivz,
 658                               TVector3&          ip,
 659                               Double_t&          cent,
 660                               UShort_t&          nClusters)
 661 {
 662   //
 663   // Process the event
 664   //
 665   // Parameters:
 666   //   event     Input event
 667   //   triggers  On return, the triggers fired
 668   //   lowFlux   On return, true if the event is considered a low-flux
 669   //                  event (according to the setting of fLowFluxCut)
 670   //   ivz       On return, the found vertex bin (1-based).  A zero
 671   //                  means outside of the defined vertex range
 672   //   vz        On return, the z position of the interaction
 673   //   cent      On return, the centrality - if not available < 0
 674   //
 675   // Return:
 676   //    0 (or kOk) on success, otherwise a bit mask of error codes
 677   //
 678   DGUARD(fDebug,1,"Process event in AliFMDEventInspector");
 679   // --- Check that we have an event ---------------------------------
 680   if (!event) {
 681     AliWarning("No ESD event found for input event");
 682     fHStatus->Fill(2);
 683     return kNoEvent;
 684   }
 685
 686   // --- Process satellite event information is requested ------------
 687   if (fUseDisplacedVertices) {
 688     if (!fDisplacedVertex.Process(event))
 689       AliWarning("Failed to process satellite event");
 690   }
 691
 692   // --- Read trigger information from the ESD and store in AOD object
 693   if (!ReadTriggers(*event, triggers, nClusters)) {
 694     if (fDebug > 2) {
 695       AliWarning("Failed to read triggers from ESD"); }
 696     fHStatus->Fill(3);
 697     return kNoTriggers;
 698   }
 699
 700   // --- Check if this is a high-flux event --------------------------
 701   const AliMultiplicity* testmult = event->GetMultiplicity();
 702   if (!testmult) {
 703     if (fDebug > 3) {
 704       AliWarning("No central multiplicity object found"); }
 705   }
 706   else
 707     lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 708
 709   fHType->Fill(lowFlux ? 0 : 1);
 710
 711   // --- Get the interaction point -----------------------------------
 712   Bool_t vzOk = ReadVertex(*event, ip);
 713   fHEventsTr->Fill(ip.Z());
 714   if (!vzOk) {
 715     if (fDebug > 3) {
 716       AliWarning("Failed to read vertex from ESD"); }
 717     fHStatus->Fill(6);
 718     return kNoVertex;
 719   }
 720
 721    // --- Read centrality information
 722   cent          = -10;
 723   UShort_t qual = 0;
 724   if (!ReadCentrality(*event, cent, qual)) {
 725     if (fDebug > 3)
 726       AliWarning("Failed to get centrality");
 727   }
 728   // --- check centrality cut
 729
 730   if(fMinCent > -0.0001 && cent < fMinCent) return  kNoEvent;
 731   if(fMaxCent > -0.0001 && cent > fMaxCent) return  kNoEvent;
 732   fHCent->Fill(cent);
 733   if (qual == 0) fHCentVsQual->Fill(0., cent);
 734   else {
 735     for (UShort_t i = 0; i < 4; i++)
 736       if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
 737   }
 738
 739
 740   fHEventsTrVtx->Fill(ip.Z());
 741
 742   // --- Get the vertex bin ------------------------------------------
 743   ivz = fVtxAxis.FindBin(ip.Z());
 744   if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
 745     if (fDebug > 3) {
 746       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 747                       ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
 748     }
 749     ivz = 0;
 750     fHStatus->Fill(7);
 751     return kBadVertex;
 752   }
 753   fHEventsAccepted->Fill(ip.Z());
 754   fHEventsAcceptedXY->Fill(ip.X(),ip.Y());
 755
 756   // --- Check the FMD ESD data --------------------------------------
 757   if (!event->GetFMDData()) {
 758     if (fDebug > 3) {
 759       AliWarning("No FMD data found in ESD"); }
 760     fHStatus->Fill(5);
 761     return kNoFMD;
 762   }
 763
 764   fHStatus->Fill(1);
 765   return kOk;
 766 }
 767
 768 //____________________________________________________________________
 769 Bool_t
 770 AliFMDEventInspector::ReadCentrality(const AliESDEvent& esd,
 771                                      Double_t& cent,
 772                                      UShort_t& qual) const
 773 {
 774   //
 775   // Read centrality from event
 776   //
 777   // Parameters:
 778   //    esd  Event
 779   //    cent On return, the centrality or negative if not found
 780   //
 781   // Return:
 782   //    False on error, true otherwise
 783   //
 784   DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
 785
 786   cent = -1;
 787   qual = 1;
 788   AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
 789   if (centObj) {
 790     cent = centObj->GetCentralityPercentile(fCentMethod);
 791     qual = centObj->GetQuality();
 792   }
 793
 794   // We overwrite with satellite events, so we can be sure to get the
 795   // centrality determination from the satellite vertex selection
 796   if (fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
 797     cent = fDisplacedVertex.GetCentralityPercentile();
 798     qual = 0;
 799   }
 800
 801   return true;
 802 }
 803
 804 //____________________________________________________________________
 805 Bool_t
 806 AliFMDEventInspector::CheckpAExtraV0(const AliESDEvent& esd) const
 807 {
 808   if (fCollisionSystem != AliForwardUtil::kPPb) return true;
 809
 810    AliVVZERO* esdV0 = esd.GetVZEROData();
 811    if ((esdV0->GetV0ADecision()!=1) || (esdV0->GetV0CDecision()!=1))
 812      return false;
 813    return true;
 814 }
 815
 816 //____________________________________________________________________
 817 Bool_t
 818 AliFMDEventInspector::ReadTriggers(const AliESDEvent& esd, UInt_t& triggers,
 819                                    UShort_t& nClusters)
 820 {
 821   //
 822   // Read the trigger information from the ESD event
 823   //
 824   // Parameters:
 825   //   esd        ESD event
 826   //   triggers   On return, contains the trigger bits
 827   //
 828   // Return:
 829   //    @c true on success, @c false otherwise
 830   //
 831   DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
 832   triggers = 0;
 833
 834   // Get the analysis manager - should always be there
 835   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 836   DMSG(fDebug,10,"Got analysis manager %p", am);
 837   if (!am) {
 838     AliWarning("No analysis manager defined!");
 839     return kFALSE;
 840   }
 841
 842   // Get the input handler - should always be there
 843   AliInputEventHandler* ih =
 844     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 845   DMSG(fDebug,10,"Got input handler %p", ih);
 846   if (!ih) {
 847     AliWarning("No input handler");
 848     return kFALSE;
 849   }
 850
 851   // Check if this is a collision candidate (MB)
 852   ///
 853   // Historic remark: Note, that we should use the value cached in the
 854   //   input handler rather than calling IsCollisionCandiate directly
 855   //   on the AliPhysicsSelection obejct.  If we called the latter
 856   //   then the AliPhysicsSelection object would overcount by a factor
 857   //   of 2! :-(
 858   UInt_t  trgMask  = ih->IsEventSelected();
 859   Bool_t  offline  = trgMask;
 860   Bool_t  fastonly = (trgMask & AliVEvent::kFastOnly);
 861   TString trigStr  = esd.GetFiredTriggerClasses();
 862
 863   if (trigStr.IsNull()) fHTrgStatus->Fill(kNoTrgWords);
 864   if (fHWords) fHWords->Fill(trigStr.Data(), 1);
 865
 866   if(fUseDisplacedVertices) {
 867     DMSG(fDebug,3,"Using displaced vertex stuff");
 868     // if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
 869     if (fDisplacedVertex.IsSatellite())
 870       triggers |= AliAODForwardMult::kSatellite;
 871   }
 872
 873   if (CheckFastPartition(fastonly)) {
 874     fHTrgStatus->Fill(kPP2760Fast);
 875     offline = false;
 876   }
 877
 878   if (offline && CheckCosmics(trigStr)) {
 879     fHTrgStatus->Fill(kMUON);
 880     offline = false;
 881   }
 882   if (offline) fHTrgStatus->Fill(kTriggered);
 883   Int_t f = 0;
 884   if (trgMask & AliVEvent::kMB)          f += fHTrgStatus->Fill(kMinBias);
 885   if (trgMask & AliVEvent::kCINT5)       f += fHTrgStatus->Fill(kMinBiasNoSPD);
 886   if (trgMask & AliVEvent::kINT7)        f += fHTrgStatus->Fill(kV0AndTrg);
 887   if (trgMask & AliVEvent::kHighMult)    f += fHTrgStatus->Fill(kHighMult);
 888   if (trgMask & AliVEvent::kCentral)     f += fHTrgStatus->Fill(kCentral);
 889   if (trgMask & AliVEvent::kSemiCentral) f += fHTrgStatus->Fill(kSemiCentral);
 890   if (trgMask & AliVEvent::kDG5)         f += fHTrgStatus->Fill(kDiffractive);
 891   if (trgMask & AliVEvent::kUserDefined) f += fHTrgStatus->Fill(kUser);
 892   if (f <= 0) fHTrgStatus->Fill(kOther);
 893
 894   // if (!CheckpAExtraV0(esd))             offline = false;
 895
 896   DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
 897
 898   if (offline)  triggers |= AliAODForwardMult::kOffline;
 899   // Only flag as in-elastic if a min-bias trigger is here
 900   if (trgMask & AliVEvent::kMB) {
 901     triggers |= AliAODForwardMult::kInel;
 902     CheckINELGT0(esd, nClusters, triggers);
 903   }
 904
 905   CheckNSD(esd,triggers);
 906   CheckPileup(esd, triggers);
 907   CheckEmpty(trigStr, triggers);
 908   // if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
 909   // if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
 910
 911   CheckWords(esd, triggers);
 912
 913 #define TEST_TRIG_BIN(RET,BIN,TRIGGERS) \
 914   do { switch (BIN) { \
 915     case kInel:     RET = triggers & AliAODForwardMult::kInel;      break; \
 916     case kInelGt0:  RET = triggers & AliAODForwardMult::kInelGt0;   break; \
 917     case kNSD:      RET = triggers & AliAODForwardMult::kNSD;       break; \
 918     case kV0AND:    RET = triggers & AliAODForwardMult::kV0AND;     break; \
 919     case kEmpty:    RET = triggers & AliAODForwardMult::kEmpty;     break; \
 920     case kA:        RET = triggers & AliAODForwardMult::kA;         break; \
 921     case kB:        RET = triggers & AliAODForwardMult::kB;         break; \
 922     case kC:        RET = triggers & AliAODForwardMult::kC;         break; \
 923     case kE:        RET = triggers & AliAODForwardMult::kE;         break; \
 924     case kPileUp:   RET = triggers & AliAODForwardMult::kPileUp;    break; \
 925     case kMCNSD:    RET = triggers & AliAODForwardMult::kMCNSD;     break; \
 926     case kSatellite:RET = triggers & AliAODForwardMult::kSatellite; break; \
 927     case kOffline:  RET = triggers & AliAODForwardMult::kOffline;   break; \
 928     default:        RET = false; } } while(false)
 929
 930   if (!fHTriggers) {
 931     AliWarning("Histogram of triggers not defined - has init been called");
 932     return false;
 933   }
 934
 935   for (Int_t i = 0; i < kOffline+1; i++) {
 936     Bool_t hasX = false;
 937     TEST_TRIG_BIN(hasX, i, triggers);
 938     if (!hasX) continue;
 939     fHTriggers->Fill(i+.5);
 940     for (Int_t j = 0; j < kOffline+1; j++) {
 941       Bool_t hasY = false;
 942       TEST_TRIG_BIN(hasY, j, triggers);
 943       if (!hasY) continue;
 944
 945       fHTriggerCorr->Fill(i+.5, j+.5);
 946     }
 947   }
 948   return kTRUE;
 949 }
 950
 951 //____________________________________________________________________
 952 Bool_t
 953 AliFMDEventInspector::CheckFastPartition(bool fastonly) const
 954 {
 955   // For the 2.76 TeV p+p run, the FMD ran in the slow partition
 956   // so it received no triggers from the fast partition. Therefore
 957   // the fast triggers are removed here but not for MC where all
 958   // triggers are fast.
 959   if (TMath::Abs(fEnergy - 2750.) > 20) return false;
 960   if (fCollisionSystem != AliForwardUtil::kPP) return false;
 961   if (fMC) return false; // All MC events for pp @ 2.76TeV are `fast'
 962   if (fastonly)
 963     DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
 964
 965   return fastonly;
 966 }
 967
 968 //____________________________________________________________________
 969 Bool_t
 970 AliFMDEventInspector::CheckCosmics(const TString& trigStr) const
 971 {
 972   // MUON triggers are not strictly minimum bias (MB) so they are
 973   // removed (HHD)
 974   if(trigStr.Contains("CMUS1")) {
 975     DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
 976     return true;
 977   }
 978   return false;
 979 }
 980
 981 //____________________________________________________________________
 982 Bool_t
 983 AliFMDEventInspector::CheckINELGT0(const AliESDEvent& esd,
 984                                    UShort_t& nClusters,
 985                                    UInt_t& triggers) const
 986 {
 987   nClusters = 0;
 988
 989   // If this is inel, see if we have a tracklet
 990   const AliMultiplicity* spdmult = esd.GetMultiplicity();
 991   if (!spdmult) {
 992     AliWarning("No SPD multiplicity");
 993     return false;
 994   }
 995
 996   // Check if we have one or more tracklets
 997   // in the range -1 < eta < 1 to set the INEL>0
 998   // trigger flag.
 999   //
1000   // Also count tracklets as a single cluster
1001   Int_t n = spdmult->GetNumberOfTracklets();
1002   for (Int_t j = 0; j < n; j++) {
1003     if(TMath::Abs(spdmult->GetEta(j)) < 1) nClusters++;
1004   }
1005   // If we at this point had non-zero nClusters, it's INEL>0
1006   if (nClusters > 0) triggers |= AliAODForwardMult::kInelGt0;
1007
1008   // Loop over single clusters
1009   n = spdmult->GetNumberOfSingleClusters();
1010   for (Int_t j = 0; j < n; j++) {
1011     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
1012     if (TMath::Abs(eta) < 1) nClusters++;
1013   }
1014   if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
1015
1016   return triggers & AliAODForwardMult::kNClusterGt0;
1017 }
1018
1019 //____________________________________________________________________
1020 Bool_t
1021 AliFMDEventInspector::CheckNSD(const AliESDEvent& esd, UInt_t& triggers) const
1022 {
1023   // Analyse some trigger stuff
1024   AliTriggerAnalysis ta;
1025   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kV0AND)) {
1026     triggers |= AliAODForwardMult::kV0AND;
1027     if (fUseV0AND)
1028       triggers |= AliAODForwardMult::kNSD;
1029   }
1030   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1))
1031     triggers |= AliAODForwardMult::kNSD;
1032   return triggers & AliAODForwardMult::kNSD;
1033 }
1034 //____________________________________________________________________
1035 Bool_t
1036 AliFMDEventInspector::CheckPileup(const AliESDEvent& esd,
1037                                   UInt_t& triggers) const
1038 {
1039   // Check for multiple vertices (pile-up) with at least 3
1040   // contributors and at least 0.8cm from the primary vertex
1041   // if(fCollisionSystem != AliForwardUtil::kPP) return false;
1042
1043   // Check for standard SPD pile-up
1044   Bool_t spdPileup = esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance);
1045
1046   // Check for multi-vertex pileup
1047   Bool_t mvPileup  = false; // CheckMultiVertex(esd);
1048
1049   // Check for out-of-bunch pileup
1050   Bool_t outPileup = (esd.GetHeader()->GetIRInt2ClosestInteractionMap() != 0 ||
1051                       esd.GetHeader()->GetIRInt1ClosestInteractionMap() != 0);
1052
1053   // Our result
1054   Bool_t pileup    = spdPileup || mvPileup || outPileup;
1055   if (pileup) triggers |= AliAODForwardMult::kPileUp;
1056   return pileup;
1057 }
1058
1059 //____________________________________________________________________
1060 Bool_t
1061 AliFMDEventInspector::CheckMultiVertex(const AliESDEvent& esd,
1062                                        Bool_t             checkOtherBC) const
1063 {
1064   // Adapted from AliAnalysisUtils
1065   //
1066   // Parameters
1067   const Double_t maxChi2nu    = 5;
1068
1069   // Number of vertices
1070   Int_t n = esd.GetNumberOfPileupVerticesTracks();
1071   if (n < 1)
1072     // No additional vertices from tracks -> no pileup
1073     return false;
1074
1075   // Primary vertex
1076   const AliESDVertex* primary = esd.GetPrimaryVertexTracks();
1077   if (primary->GetStatus() != 1)
1078     // No good primary vertex, but many candidates -> pileup
1079     return true;
1080
1081   // Bunch crossing number
1082   Int_t bcPrimary = primary->GetBC();
1083
1084   for (Int_t i = 0; i < n; i++) {
1085     const AliESDVertex* other = esd.GetPileupVertexTracks(i);
1086
1087     if (other->GetNContributors() < fMinPileupContrib)
1088       // Not enough contributors to this vertex
1089       continue;
1090
1091     if (other->GetChi2perNDF() > maxChi2nu)
1092       // Poorly determined vertex
1093       continue;
1094     if (checkOtherBC) {
1095       Int_t bcOther = other->GetBC();
1096       if (bcOther != AliVTrack::kTOFBCNA && TMath::Abs(bcOther-bcPrimary) > 2)
1097         // Pile-up from other BC
1098         return true;
1099     }
1100     // Calculate the distance
1101     Double_t d  = -1;
1102     Double_t dx = primary->GetX() - other->GetX();
1103     Double_t dy = primary->GetY() - other->GetY();
1104     Double_t dz = primary->GetZ() - other->GetZ();
1105     Double_t covPrimary[6], covOther[6];
1106     primary->GetCovarianceMatrix(covPrimary);
1107     other->GetCovarianceMatrix(covOther);
1108     TMatrixDSym v(3);
1109     v(0,0) = covPrimary[0] + covOther[0]; // diagonal
1110     v(1,1) = covPrimary[2] + covOther[2]; // diagonal
1111     v(2,2) = covPrimary[5] + covOther[5]; // diagonal
1112     v(1,0) = v(0,1) = covPrimary[1]+covOther[1]; // Band
1113     v(0,2) = v(1,2) = v(2,0) = v(2,1) = 0; // Off-diagonal+band
1114     // Try to invert
1115     v.InvertFast();
1116     if (!v.IsValid())
1117       // Question if kStatus bit is every set after InvertFast?
1118       continue;
1119     d = (v(0,0) * dx * dx + v(1,1) * dy * dy + v(2,2) * dz * dz +
1120          2 * (v(0,1) * dx * dy + v(0,2) * dx * dz + v(1,2) * dy * dz));
1121
1122     if (d < 0 || TMath::Sqrt(d) < fMinPileupDistance)
1123       // Bad distance, or not fare enough from each other
1124       continue;
1125
1126     // Well separated vertices -> pileup
1127     return true;
1128   }
1129   return false;
1130 }
1131
1132 //____________________________________________________________________
1133 Bool_t
1134 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
1135 {
1136   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
1137     triggers |= AliAODForwardMult::kEmpty;
1138     return true;
1139   }
1140   return false;
1141 }
1142 //____________________________________________________________________
1143 Bool_t
1144 AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const
1145 {
1146   TObject* word = 0;
1147   TIter nextColl(&fCollWords);
1148   while ((word = nextColl())) {
1149     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
1150          word->GetTitle(), word->GetName());
1151     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1152
1153     TString beamSide = word->GetTitle();
1154     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1155
1156     if (!beamSide.EqualTo("B")) continue;
1157     triggers |= AliAODForwardMult::kB;
1158     break; // No more to do here
1159   }
1160   TIter nextBg(&fBgWords);
1161   UInt_t all = (AliAODForwardMult::kA |
1162                 AliAODForwardMult::kC |
1163                 AliAODForwardMult::kE);
1164   while ((word = nextBg())) {
1165     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
1166          word->GetTitle(), word->GetName());
1167     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1168
1169     TString beamSide = word->GetTitle();
1170     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1171
1172     if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
1173     if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
1174     if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
1175
1176     if ((triggers & all) == all) break; // No more to do
1177   }
1178   return true;
1179 }
1180
1181
1182 //____________________________________________________________________
1183 Bool_t
1184 AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, TVector3& ip)
1185 {
1186   //
1187   // Read the vertex information from the ESD event
1188   //
1189   // Parameters:
1190   //   esd  ESD event
1191   //   vz   On return, the vertex Z position
1192   //
1193   // Return:
1194   //    @c true on success, @c false otherwise
1195   //
1196   DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
1197   ip.SetXYZ(1024, 1024, 0);
1198
1199   EVtxStatus s = kNoVtx;
1200   if (fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
1201     s = kVtxOK;
1202     ip.SetZ(fDisplacedVertex.GetVertexZ());
1203   }
1204   else if (fUseFirstPhysicsVertex)
1205     s = CheckPWGUDVertex(esd, ip);
1206   else if (fUsepA2012Vertex)
1207     s = CheckpA2012Vertex(esd,ip);
1208   else
1209     s = CheckVertex(esd, ip);
1210
1211   fHVtxStatus->Fill(s);
1212
1213   return s == kVtxOK;
1214 }
1215
1216 //____________________________________________________________________
1217 AliFMDEventInspector::EVtxStatus
1218 AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd,
1219                                        TVector3& ip)  const
1220 {
1221   // This is the code used by the 1st physics people
1222   const AliESDVertex* vertex    = esd.GetPrimaryVertex();
1223   if (!vertex  || !vertex->GetStatus()) {
1224     DMSG(fDebug,2,"No primary vertex (%p) or bad status %d",
1225          vertex, (vertex ? vertex->GetStatus() : -1));
1226     return kNoVtx;
1227   }
1228   const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
1229   if (!vertexSPD || !vertexSPD->GetStatus()) {
1230     DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d",
1231          vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
1232     return kNoSPDVtx;
1233   }
1234
1235   // if vertex is from SPD vertexZ, require more stringent cuts
1236   if (vertex->IsFromVertexerZ()) {
1237     if (vertex->GetDispersion() > fMaxVzErr ||
1238         vertex->GetZRes() > 1.25 * fMaxVzErr) {
1239       DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
1240            vertex->GetDispersion(), fMaxVzErr,
1241            vertex->GetZRes(), 1.25 * fMaxVzErr);
1242       return kUncertain;
1243     }
1244   }
1245   ip.SetZ(vertex->GetZ());
1246
1247   if(!vertex->IsFromVertexerZ()) {
1248     ip.SetX(vertex->GetX());
1249     ip.SetY(vertex->GetY());
1250   }
1251   return kVtxOK;
1252 }
1253 //____________________________________________________________________
1254 AliFMDEventInspector::EVtxStatus
1255 AliFMDEventInspector::CheckpA2012Vertex(const AliESDEvent& esd,
1256                                         TVector3& ip)  const
1257 {
1258   const AliESDVertex *vertex = esd.GetPrimaryVertexSPD();
1259   if (!vertex) return kNoSPDVtx;
1260   if (vertex->GetNContributors() <= 0) return kFewContrib;
1261
1262   TString vtxTyp = vertex->GetTitle();
1263   if (vtxTyp.Contains("vertexer: Z")) return kNotVtxZ;
1264
1265   if (vertex->GetDispersion() >= 0.04 || vertex->GetZRes()>=0.25)
1266     return kUncertain;
1267
1268   ip.SetX(vertex->GetX());
1269   ip.SetY(vertex->GetY());
1270   ip.SetZ(vertex->GetZ());
1271
1272   return kVtxOK;
1273 }
1274
1275 //____________________________________________________________________
1276 AliFMDEventInspector::EVtxStatus
1277 AliFMDEventInspector::CheckVertex(const AliESDEvent& esd,
1278                                   TVector3& ip) const
1279 {
1280   // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1281   // Get the vertex
1282   const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1283   if (!vertex) {
1284     if (fDebug > 2) {
1285       AliWarning("No SPD vertex found in ESD"); }
1286     return kNoSPDVtx;
1287   }
1288
1289   // #if 0 // Check disabled - seem to kill a lot of PbPb events
1290   // Check that enough tracklets contributed
1291   if(vertex->GetNContributors() <= 0) {
1292     DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1293          vertex->GetNContributors());
1294     ip.SetZ(0);
1295     return kFewContrib;
1296   }
1297   // #endif
1298
1299   // Check that the uncertainty isn't too large
1300   if (vertex->GetZRes() > fMaxVzErr) {
1301     DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f",
1302          vertex->GetZRes(), fMaxVzErr);
1303     return kUncertain;
1304   }
1305
1306   // Get the z coordiante
1307   ip.SetZ(vertex->GetZ());
1308   const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1309
1310
1311   if(!vertexXY->IsFromVertexerZ()) {
1312     ip.SetX(vertexXY->GetX());
1313     ip.SetY(vertexXY->GetY());
1314   }
1315   return kVtxOK;
1316 }
1317
1318 //____________________________________________________________________
1319 Bool_t
1320 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
1321 {
1322   //
1323   // Read the collision system, collision energy, and L3 field setting
1324   // from the ESD
1325   //
1326   // Parameters:
1327   //   esd ESD to get information from
1328   //
1329   // Return:
1330   //    true on success, false
1331   //
1332   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1333   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(),
1334   //           esd->GetMagneticField()));
1335   DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1336   const char* sys  = esd->GetBeamType();
1337   Float_t     cms  = 2 * esd->GetBeamEnergy();
1338   Float_t     fld  = esd->GetMagneticField();
1339   fCollisionSystem = AliForwardUtil::ParseCollisionSystem(sys);
1340   fEnergy          = AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,
1341                                                              cms);
1342   fField           = AliForwardUtil::ParseMagneticField(fld);
1343   fRunNumber       = esd->GetRunNumber();
1344   StoreProduction();
1345   StoreInformation();
1346   if (fCollisionSystem   == AliForwardUtil::kUnknown) {
1347     AliWarningF("Unknown collision system: %s - please check", sys);
1348     return false;
1349   }
1350   if (fEnergy            <= 0) {
1351     AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1352     return false;
1353   }
1354   if (TMath::Abs(fField) >  10) {
1355     AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1356     return false;
1357   }
1358
1359   return true;
1360 }
1361
1362
1363 //____________________________________________________________________
1364 const Char_t*
1365 AliFMDEventInspector::CodeString(UInt_t code)
1366 {
1367   static TString s;
1368   s = "";
1369   if (code & kNoEvent)    s.Append("NOEVENT ");
1370   if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1371   if (code & kNoSPD)      s.Append("NOSPD ");
1372   if (code & kNoFMD)      s.Append("NOFMD ");
1373   if (code & kNoVertex)   s.Append("NOVERTEX ");
1374   if (code & kBadVertex)  s.Append("BADVERTEX ");
1375   return s.Data();
1376 }
1377 #define PF(N,V,...)                                     \
1378   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1379 #define PFB(N,FLAG)                             \
1380   do {                                                                  \
1381     AliForwardUtil::PrintName(N);                                       \
1382     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1383   } while(false)
1384 #define PFV(N,VALUE)                                    \
1385   do {                                                  \
1386     AliForwardUtil::PrintName(N);                       \
1387     std::cout << (VALUE) << std::endl; } while(false)
1388
1389 //____________________________________________________________________
1390 void
1391 AliFMDEventInspector::Print(Option_t*) const
1392 {
1393   //
1394   // Print information
1395   //
1396   //   option Not used
1397   //
1398   AliForwardUtil::PrintTask(*this);
1399   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
1400   sNN.Strip(TString::kBoth, '0');
1401   sNN.ReplaceAll("GeV", " GeV");
1402   TString field(AliForwardUtil::MagneticFieldString(fField));
1403   field.ReplaceAll("p",  "+");
1404   field.ReplaceAll("m",  "-");
1405   field.ReplaceAll("kG", " kG");
1406
1407   gROOT->IncreaseDirLevel();
1408   PFV("Production", "");
1409   gROOT->IncreaseDirLevel();
1410   PF("Period", "LHC%02d%c", fProdYear, fProdLetter);
1411   PFB("MC anchor", fProdMC);
1412   PFV("AliROOT Revision", fProdSVN);
1413   gROOT->DecreaseDirLevel();
1414   PFV("Vertex bins", fVtxAxis.GetNbins());
1415   PF("Vertex Range", "[%+6.1f,%+6.1f]", fVtxAxis.GetXmin(), fVtxAxis.GetXmax());
1416   PFV("Low flux cut",           fLowFluxCut );
1417   PFV("Max(delta v_z)",         fMaxVzErr );
1418   PFV("Min(nContrib_pileup)",   fMinPileupContrib );
1419   PFV("Min(v-pileup)",          fMinPileupDistance );
1420   PFV("System", AliForwardUtil::CollisionSystemString(fCollisionSystem));
1421   PFV("CMS energy per nucleon", sNN);
1422   PFV("Field",                  field);
1423   PFB("Satellite events",       fUseDisplacedVertices);
1424   PFB("Use 2012 pA vertex",     fUsepA2012Vertex );
1425   PFB("Use PWG-UD vertex",      fUseFirstPhysicsVertex);
1426   PFB("Simulation input",       fMC );
1427   PFV("Centrality method",      fCentMethod);
1428   PFV("Centrality axis",        (!fCentAxis ? "none" : ""));
1429   if (!fCentAxis) {return; }
1430   Int_t nBin = fCentAxis->GetNbins();
1431   for (Int_t i = 0; i < nBin; i++) {
1432     if (i != 0 && (i % 10) == 0) std::cout << '\n';
1433     if ((i % 10) == 0)           std::cout << "    ";
1434     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1435   }
1436   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1437   gROOT->DecreaseDirLevel();
1438 }
1439
1440
1441 //
1442 // EOF
1443 //
1444