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