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   if (period.Length() > 0) {
 646     TObjArray* pp = TPRegexp("LHC([0-9]+)([a-z]+)").MatchS(period);
 647     fProdYear     = static_cast<TObjString*>(pp->At(1))->String().Atoi();
 648     fProdLetter   = static_cast<TObjString*>(pp->At(2))->String()[0];
 649     pp->Delete();
 650   }
 651
 652   out->Add(AliForwardUtil::MakeParameter("year", fProdYear));
 653   out->Add(AliForwardUtil::MakeParameter("letter", Int_t(fProdLetter)));
 654   out->Add(AliForwardUtil::MakeParameter("alirootSVN", fProdSVN));
 655   out->Add(AliForwardUtil::MakeParameter("rootSVN", rootSVN));
 656   out->Add(AliForwardUtil::MakeParameter("pass", fProdPass));
 657   out->Add(AliForwardUtil::MakeParameter("mc", fProdMC));
 658 }
 659
 660
 661
 662
 663 //____________________________________________________________________
 664 void
 665 AliFMDEventInspector::CreateOutputObjects(TList* dir)
 666 {
 667   //
 668   // Define the output histograms.  These are put in a sub list of the
 669   // passed list.   The histograms are merged before the parent task calls
 670   // AliAnalysisTaskSE::Terminate
 671   //
 672   //   dir Directory to add to
 673   //
 674   DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
 675   fList = new TList;
 676   fList->SetName(GetName());
 677   fList->SetOwner();
 678   dir->Add(fList);
 679 }
 680
 681 //____________________________________________________________________
 682 UInt_t
 683 AliFMDEventInspector::Process(const AliESDEvent* event,
 684                               UInt_t&            triggers,
 685                               Bool_t&            lowFlux,
 686                               UShort_t&          ivz,
 687                               TVector3&          ip,
 688                               Double_t&          cent,
 689                               UShort_t&          nClusters)
 690 {
 691   //
 692   // Process the event
 693   //
 694   // Parameters:
 695   //   event     Input event
 696   //   triggers  On return, the triggers fired
 697   //   lowFlux   On return, true if the event is considered a low-flux
 698   //                  event (according to the setting of fLowFluxCut)
 699   //   ivz       On return, the found vertex bin (1-based).  A zero
 700   //                  means outside of the defined vertex range
 701   //   vz        On return, the z position of the interaction
 702   //   cent      On return, the centrality - if not available < 0
 703   //
 704   // Return:
 705   //    0 (or kOk) on success, otherwise a bit mask of error codes
 706   //
 707   DGUARD(fDebug,1,"Process event in AliFMDEventInspector");
 708   // --- Check that we have an event ---------------------------------
 709   if (!event) {
 710     AliWarning("No ESD event found for input event");
 711     fHStatus->Fill(2);
 712     return kNoEvent;
 713   }
 714
 715   // --- Process satellite event information is requested ------------
 716   if (fUseDisplacedVertices) {
 717     if (!fDisplacedVertex.Process(event))
 718       AliWarning("Failed to process satellite event");
 719   }
 720
 721   // --- Read trigger information from the ESD and store in AOD object
 722   if (!ReadTriggers(*event, triggers, nClusters)) {
 723     if (fDebug > 2) {
 724       AliWarning("Failed to read triggers from ESD"); }
 725     fHStatus->Fill(3);
 726     return kNoTriggers;
 727   }
 728
 729   // --- Check if this is a high-flux event --------------------------
 730   const AliMultiplicity* testmult = event->GetMultiplicity();
 731   if (!testmult) {
 732     if (fDebug > 3) {
 733       AliWarning("No central multiplicity object found"); }
 734   }
 735   else
 736     lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 737
 738   fHType->Fill(lowFlux ? 0 : 1);
 739
 740   // --- Get the interaction point -----------------------------------
 741   Bool_t vzOk = ReadVertex(*event, ip);
 742   fHEventsTr->Fill(ip.Z());
 743   if (!vzOk) {
 744     if (fDebug > 3) {
 745       AliWarning("Failed to read vertex from ESD"); }
 746     fHStatus->Fill(6);
 747     return kNoVertex;
 748   }
 749
 750    // --- Read centrality information
 751   cent          = -10;
 752   UShort_t qual = 0;
 753   if (!ReadCentrality(*event, cent, qual)) {
 754     if (fDebug > 3)
 755       AliWarning("Failed to get centrality");
 756   }
 757   // --- check centrality cut
 758
 759   if(fMinCent > -0.0001 && cent < fMinCent) return  kNoEvent;
 760   if(fMaxCent > -0.0001 && cent > fMaxCent) return  kNoEvent;
 761   if (cent >= 0) fHCent->Fill(cent);
 762   if (qual == 0) fHCentVsQual->Fill(0., cent);
 763   else {
 764     fHStatus->Fill(8);
 765     for (UShort_t i = 0; i < 4; i++)
 766       if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
 767   }
 768
 769
 770   fHEventsTrVtx->Fill(ip.Z());
 771
 772   // --- Get the vertex bin ------------------------------------------
 773   ivz = fVtxAxis.FindBin(ip.Z());
 774   if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
 775     if (fDebug > 3) {
 776       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 777                       ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
 778     }
 779     ivz = 0;
 780     fHStatus->Fill(7);
 781     return kBadVertex;
 782   }
 783   fHEventsAccepted->Fill(ip.Z());
 784   fHEventsAcceptedXY->Fill(ip.X(),ip.Y());
 785
 786   // --- Check the FMD ESD data --------------------------------------
 787   if (!event->GetFMDData()) {
 788     if (fDebug > 3) {
 789       AliWarning("No FMD data found in ESD"); }
 790     fHStatus->Fill(5);
 791     return kNoFMD;
 792   }
 793
 794   fHStatus->Fill(1);
 795   return kOk;
 796 }
 797
 798 //____________________________________________________________________
 799 Bool_t
 800 AliFMDEventInspector::ReadCentrality(const AliESDEvent& esd,
 801                                      Double_t& cent,
 802                                      UShort_t& qual) const
 803 {
 804   //
 805   // Read centrality from event
 806   //
 807   // Parameters:
 808   //    esd  Event
 809   //    cent On return, the centrality or negative if not found
 810   //
 811   // Return:
 812   //    False on error, true otherwise
 813   //
 814   DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
 815
 816   cent = -1;
 817   qual = 1;
 818   AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
 819   if (centObj) {
 820     cent = centObj->GetCentralityPercentile(fCentMethod);
 821     qual = centObj->GetQuality();
 822   }
 823
 824   // We overwrite with satellite events, so we can be sure to get the
 825   // centrality determination from the satellite vertex selection
 826   if (fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
 827     cent = fDisplacedVertex.GetCentralityPercentile();
 828     qual = 0;
 829   }
 830
 831   return true;
 832 }
 833
 834 //____________________________________________________________________
 835 Bool_t
 836 AliFMDEventInspector::CheckpAExtraV0(const AliESDEvent& esd) const
 837 {
 838   if (fCollisionSystem != AliForwardUtil::kPPb) return true;
 839
 840    AliVVZERO* esdV0 = esd.GetVZEROData();
 841    if ((esdV0->GetV0ADecision()!=1) || (esdV0->GetV0CDecision()!=1))
 842      return false;
 843    return true;
 844 }
 845
 846 //____________________________________________________________________
 847 Bool_t
 848 AliFMDEventInspector::ReadTriggers(const AliESDEvent& esd, UInt_t& triggers,
 849                                    UShort_t& nClusters)
 850 {
 851   //
 852   // Read the trigger information from the ESD event
 853   //
 854   // Parameters:
 855   //   esd        ESD event
 856   //   triggers   On return, contains the trigger bits
 857   //
 858   // Return:
 859   //    @c true on success, @c false otherwise
 860   //
 861   DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
 862   triggers = 0;
 863
 864   // Get the analysis manager - should always be there
 865   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 866   DMSG(fDebug,10,"Got analysis manager %p", am);
 867   if (!am) {
 868     AliWarning("No analysis manager defined!");
 869     return kFALSE;
 870   }
 871
 872   // Get the input handler - should always be there
 873   AliInputEventHandler* ih =
 874     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 875   DMSG(fDebug,10,"Got input handler %p", ih);
 876   if (!ih) {
 877     AliWarning("No input handler");
 878     return kFALSE;
 879   }
 880
 881   // Check if this is a collision candidate (MB)
 882   ///
 883   // Historic remark: Note, that we should use the value cached in the
 884   //   input handler rather than calling IsCollisionCandiate directly
 885   //   on the AliPhysicsSelection obejct.  If we called the latter
 886   //   then the AliPhysicsSelection object would overcount by a factor
 887   //   of 2! :-(
 888   UInt_t  trgMask  = ih->IsEventSelected();
 889   Bool_t  offline  = trgMask;
 890   Bool_t  fastonly = (trgMask & AliVEvent::kFastOnly);
 891   TString trigStr  = esd.GetFiredTriggerClasses();
 892
 893   if (trigStr.IsNull()) fHTrgStatus->Fill(kNoTrgWords);
 894   if (fHWords) fHWords->Fill(trigStr.Data(), 1);
 895
 896   if(fUseDisplacedVertices) {
 897     DMSG(fDebug,3,"Using displaced vertex stuff");
 898     // if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
 899     if (fDisplacedVertex.IsSatellite())
 900       triggers |= AliAODForwardMult::kSatellite;
 901   }
 902
 903   if (CheckFastPartition(fastonly)) {
 904     fHTrgStatus->Fill(kPP2760Fast);
 905     offline = false;
 906   }
 907
 908   if (offline && CheckCosmics(trigStr)) {
 909     fHTrgStatus->Fill(kMUON);
 910     offline = false;
 911   }
 912   if (offline) fHTrgStatus->Fill(kTriggered);
 913   Int_t f = 0;
 914   if (trgMask & AliVEvent::kMB)          f += fHTrgStatus->Fill(kMinBias);
 915   if (trgMask & AliVEvent::kCINT5)       f += fHTrgStatus->Fill(kMinBiasNoSPD);
 916   if (trgMask & AliVEvent::kINT7)        f += fHTrgStatus->Fill(kV0AndTrg);
 917   if (trgMask & AliVEvent::kHighMult)    f += fHTrgStatus->Fill(kHighMult);
 918   if (trgMask & AliVEvent::kCentral)     f += fHTrgStatus->Fill(kCentral);
 919   if (trgMask & AliVEvent::kSemiCentral) f += fHTrgStatus->Fill(kSemiCentral);
 920   if (trgMask & AliVEvent::kDG5)         f += fHTrgStatus->Fill(kDiffractive);
 921   if (trgMask & AliVEvent::kUserDefined) f += fHTrgStatus->Fill(kUser);
 922   if (f <= 0) fHTrgStatus->Fill(kOther);
 923
 924   // if (!CheckpAExtraV0(esd))             offline = false;
 925
 926   DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
 927
 928   if (offline)  triggers |= AliAODForwardMult::kOffline;
 929   // Only flag as in-elastic if a min-bias trigger is here
 930   if (trgMask & AliVEvent::kMB) {
 931     triggers |= AliAODForwardMult::kInel;
 932     CheckINELGT0(esd, nClusters, triggers);
 933   }
 934
 935   CheckNSD(esd,triggers);
 936   CheckPileup(esd, triggers);
 937   CheckEmpty(trigStr, triggers);
 938   // if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
 939   // if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
 940
 941   CheckWords(esd, triggers);
 942
 943 #define TEST_TRIG_BIN(RET,BIN,TRIGGERS) \
 944   do { switch (BIN) { \
 945     case kInel:     RET = triggers & AliAODForwardMult::kInel;      break; \
 946     case kInelGt0:  RET = triggers & AliAODForwardMult::kInelGt0;   break; \
 947     case kNSD:      RET = triggers & AliAODForwardMult::kNSD;       break; \
 948     case kV0AND:    RET = triggers & AliAODForwardMult::kV0AND;     break; \
 949     case kEmpty:    RET = triggers & AliAODForwardMult::kEmpty;     break; \
 950     case kA:        RET = triggers & AliAODForwardMult::kA;         break; \
 951     case kB:        RET = triggers & AliAODForwardMult::kB;         break; \
 952     case kC:        RET = triggers & AliAODForwardMult::kC;         break; \
 953     case kE:        RET = triggers & AliAODForwardMult::kE;         break; \
 954     case kPileUp:   RET = triggers & AliAODForwardMult::kPileUp;    break; \
 955     case kMCNSD:    RET = triggers & AliAODForwardMult::kMCNSD;     break; \
 956     case kSatellite:RET = triggers & AliAODForwardMult::kSatellite; break; \
 957     case kOffline:  RET = triggers & AliAODForwardMult::kOffline;   break; \
 958     default:        RET = false; } } while(false)
 959
 960   if (!fHTriggers) {
 961     AliWarning("Histogram of triggers not defined - has init been called");
 962     return false;
 963   }
 964
 965   for (Int_t i = 0; i < kOffline+1; i++) {
 966     Bool_t hasX = false;
 967     TEST_TRIG_BIN(hasX, i, triggers);
 968     if (!hasX) continue;
 969     fHTriggers->Fill(i+.5);
 970     for (Int_t j = 0; j < kOffline+1; j++) {
 971       Bool_t hasY = false;
 972       TEST_TRIG_BIN(hasY, j, triggers);
 973       if (!hasY) continue;
 974
 975       fHTriggerCorr->Fill(i+.5, j+.5);
 976     }
 977   }
 978   return kTRUE;
 979 }
 980
 981 //____________________________________________________________________
 982 Bool_t
 983 AliFMDEventInspector::CheckFastPartition(bool fastonly) const
 984 {
 985   // For the 2.76 TeV p+p run, the FMD ran in the slow partition
 986   // so it received no triggers from the fast partition. Therefore
 987   // the fast triggers are removed here but not for MC where all
 988   // triggers are fast.
 989   if (TMath::Abs(fEnergy - 2750.) > 20) return false;
 990   if (fCollisionSystem != AliForwardUtil::kPP) return false;
 991   if (fMC) return false; // All MC events for pp @ 2.76TeV are `fast'
 992   if (fastonly)
 993     DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
 994
 995   return fastonly;
 996 }
 997
 998 //____________________________________________________________________
 999 Bool_t
1000 AliFMDEventInspector::CheckCosmics(const TString& trigStr) const
1001 {
1002   // MUON triggers are not strictly minimum bias (MB) so they are
1003   // removed (HHD)
1004   if(trigStr.Contains("CMUS1")) {
1005     DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
1006     return true;
1007   }
1008   return false;
1009 }
1010
1011 //____________________________________________________________________
1012 Bool_t
1013 AliFMDEventInspector::CheckINELGT0(const AliESDEvent& esd,
1014                                    UShort_t& nClusters,
1015                                    UInt_t& triggers) const
1016 {
1017   nClusters = 0;
1018
1019   // If this is inel, see if we have a tracklet
1020   const AliMultiplicity* spdmult = esd.GetMultiplicity();
1021   if (!spdmult) {
1022     AliWarning("No SPD multiplicity");
1023     return false;
1024   }
1025
1026   // Check if we have one or more tracklets
1027   // in the range -1 < eta < 1 to set the INEL>0
1028   // trigger flag.
1029   //
1030   // Also count tracklets as a single cluster
1031   Int_t n = spdmult->GetNumberOfTracklets();
1032   for (Int_t j = 0; j < n; j++) {
1033     if(TMath::Abs(spdmult->GetEta(j)) < 1) nClusters++;
1034   }
1035   // If we at this point had non-zero nClusters, it's INEL>0
1036   if (nClusters > 0) triggers |= AliAODForwardMult::kInelGt0;
1037
1038   // Loop over single clusters
1039   n = spdmult->GetNumberOfSingleClusters();
1040   for (Int_t j = 0; j < n; j++) {
1041     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
1042     if (TMath::Abs(eta) < 1) nClusters++;
1043   }
1044   if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
1045
1046   return triggers & AliAODForwardMult::kNClusterGt0;
1047 }
1048
1049 //____________________________________________________________________
1050 Bool_t
1051 AliFMDEventInspector::CheckNSD(const AliESDEvent& esd, UInt_t& triggers) const
1052 {
1053   // Analyse some trigger stuff
1054   AliTriggerAnalysis ta;
1055   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kV0AND)) {
1056     triggers |= AliAODForwardMult::kV0AND;
1057     if (fUseV0AND)
1058       triggers |= AliAODForwardMult::kNSD;
1059   }
1060   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1))
1061     triggers |= AliAODForwardMult::kNSD;
1062   return triggers & AliAODForwardMult::kNSD;
1063 }
1064 //____________________________________________________________________
1065 Bool_t
1066 AliFMDEventInspector::CheckPileup(const AliESDEvent& esd,
1067                                   UInt_t& triggers) const
1068 {
1069   // Check for multiple vertices (pile-up) with at least 3
1070   // contributors and at least 0.8cm from the primary vertex
1071   // if(fCollisionSystem != AliForwardUtil::kPP) return false;
1072
1073   // Check for standard SPD pile-up
1074   Bool_t spdPileup = esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance);
1075
1076   // Check for multi-vertex pileup
1077   Bool_t mvPileup  = false; // CheckMultiVertex(esd);
1078
1079   // Check for out-of-bunch pileup
1080   Bool_t outPileup = (esd.GetHeader()->GetIRInt2ClosestInteractionMap() != 0 ||
1081                       esd.GetHeader()->GetIRInt1ClosestInteractionMap() != 0);
1082
1083   // Our result
1084   Bool_t pileup    = spdPileup || mvPileup || outPileup;
1085   if (pileup) triggers |= AliAODForwardMult::kPileUp;
1086   return pileup;
1087 }
1088
1089 //____________________________________________________________________
1090 Bool_t
1091 AliFMDEventInspector::CheckMultiVertex(const AliESDEvent& esd,
1092                                        Bool_t             checkOtherBC) const
1093 {
1094   // Adapted from AliAnalysisUtils
1095   //
1096   // Parameters
1097   const Double_t maxChi2nu    = 5;
1098
1099   // Number of vertices
1100   Int_t n = esd.GetNumberOfPileupVerticesTracks();
1101   if (n < 1)
1102     // No additional vertices from tracks -> no pileup
1103     return false;
1104
1105   // Primary vertex
1106   const AliESDVertex* primary = esd.GetPrimaryVertexTracks();
1107   if (primary->GetStatus() != 1)
1108     // No good primary vertex, but many candidates -> pileup
1109     return true;
1110
1111   // Bunch crossing number
1112   Int_t bcPrimary = primary->GetBC();
1113
1114   for (Int_t i = 0; i < n; i++) {
1115     const AliESDVertex* other = esd.GetPileupVertexTracks(i);
1116
1117     if (other->GetNContributors() < fMinPileupContrib)
1118       // Not enough contributors to this vertex
1119       continue;
1120
1121     if (other->GetChi2perNDF() > maxChi2nu)
1122       // Poorly determined vertex
1123       continue;
1124     if (checkOtherBC) {
1125       Int_t bcOther = other->GetBC();
1126       if (bcOther != AliVTrack::kTOFBCNA && TMath::Abs(bcOther-bcPrimary) > 2)
1127         // Pile-up from other BC
1128         return true;
1129     }
1130     // Calculate the distance
1131     Double_t d  = -1;
1132     Double_t dx = primary->GetX() - other->GetX();
1133     Double_t dy = primary->GetY() - other->GetY();
1134     Double_t dz = primary->GetZ() - other->GetZ();
1135     Double_t covPrimary[6], covOther[6];
1136     primary->GetCovarianceMatrix(covPrimary);
1137     other->GetCovarianceMatrix(covOther);
1138     TMatrixDSym v(3);
1139     v(0,0) = covPrimary[0] + covOther[0]; // diagonal
1140     v(1,1) = covPrimary[2] + covOther[2]; // diagonal
1141     v(2,2) = covPrimary[5] + covOther[5]; // diagonal
1142     v(1,0) = v(0,1) = covPrimary[1]+covOther[1]; // Band
1143     v(0,2) = v(1,2) = v(2,0) = v(2,1) = 0; // Off-diagonal+band
1144     // Try to invert
1145     v.InvertFast();
1146     if (!v.IsValid())
1147       // Question if kStatus bit is every set after InvertFast?
1148       continue;
1149     d = (v(0,0) * dx * dx + v(1,1) * dy * dy + v(2,2) * dz * dz +
1150          2 * (v(0,1) * dx * dy + v(0,2) * dx * dz + v(1,2) * dy * dz));
1151
1152     if (d < 0 || TMath::Sqrt(d) < fMinPileupDistance)
1153       // Bad distance, or not fare enough from each other
1154       continue;
1155
1156     // Well separated vertices -> pileup
1157     return true;
1158   }
1159   return false;
1160 }
1161
1162 //____________________________________________________________________
1163 Bool_t
1164 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
1165 {
1166   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
1167     triggers |= AliAODForwardMult::kEmpty;
1168     return true;
1169   }
1170   return false;
1171 }
1172 //____________________________________________________________________
1173 Bool_t
1174 AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const
1175 {
1176   TObject* word = 0;
1177   TIter nextColl(&fCollWords);
1178   while ((word = nextColl())) {
1179     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
1180          word->GetTitle(), word->GetName());
1181     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1182
1183     TString beamSide = word->GetTitle();
1184     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1185
1186     if (!beamSide.EqualTo("B")) continue;
1187     triggers |= AliAODForwardMult::kB;
1188     break; // No more to do here
1189   }
1190   TIter nextBg(&fBgWords);
1191   UInt_t all = (AliAODForwardMult::kA |
1192                 AliAODForwardMult::kC |
1193                 AliAODForwardMult::kE);
1194   while ((word = nextBg())) {
1195     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
1196          word->GetTitle(), word->GetName());
1197     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1198
1199     TString beamSide = word->GetTitle();
1200     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1201
1202     if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
1203     if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
1204     if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
1205
1206     if ((triggers & all) == all) break; // No more to do
1207   }
1208   return true;
1209 }
1210
1211
1212 //____________________________________________________________________
1213 Bool_t
1214 AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, TVector3& ip)
1215 {
1216   //
1217   // Read the vertex information from the ESD event
1218   //
1219   // Parameters:
1220   //   esd  ESD event
1221   //   vz   On return, the vertex Z position
1222   //
1223   // Return:
1224   //    @c true on success, @c false otherwise
1225   //
1226   DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
1227   ip.SetXYZ(1024, 1024, 0);
1228
1229   EVtxStatus s = kNoVtx;
1230   if (fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
1231     s = kVtxOK;
1232     ip.SetZ(fDisplacedVertex.GetVertexZ());
1233   }
1234   else if (fUseFirstPhysicsVertex)
1235     s = CheckPWGUDVertex(esd, ip);
1236   else if (fUsepA2012Vertex)
1237     s = CheckpA2012Vertex(esd,ip);
1238   else
1239     s = CheckVertex(esd, ip);
1240
1241   fHVtxStatus->Fill(s);
1242
1243   return s == kVtxOK;
1244 }
1245
1246 //____________________________________________________________________
1247 AliFMDEventInspector::EVtxStatus
1248 AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd,
1249                                        TVector3& ip)  const
1250 {
1251   // This is the code used by the 1st physics people
1252   const AliESDVertex* vertex    = esd.GetPrimaryVertex();
1253   if (!vertex  || !vertex->GetStatus()) {
1254     DMSG(fDebug,2,"No primary vertex (%p) or bad status %d",
1255          vertex, (vertex ? vertex->GetStatus() : -1));
1256     return kNoVtx;
1257   }
1258   const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
1259   if (!vertexSPD || !vertexSPD->GetStatus()) {
1260     DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d",
1261          vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
1262     return kNoSPDVtx;
1263   }
1264
1265   // if vertex is from SPD vertexZ, require more stringent cuts
1266   if (vertex->IsFromVertexerZ()) {
1267     if (vertex->GetDispersion() > fMaxVzErr ||
1268         vertex->GetZRes() > 1.25 * fMaxVzErr) {
1269       DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
1270            vertex->GetDispersion(), fMaxVzErr,
1271            vertex->GetZRes(), 1.25 * fMaxVzErr);
1272       return kUncertain;
1273     }
1274   }
1275   ip.SetZ(vertex->GetZ());
1276
1277   if(!vertex->IsFromVertexerZ()) {
1278     ip.SetX(vertex->GetX());
1279     ip.SetY(vertex->GetY());
1280   }
1281   return kVtxOK;
1282 }
1283 //____________________________________________________________________
1284 AliFMDEventInspector::EVtxStatus
1285 AliFMDEventInspector::CheckpA2012Vertex(const AliESDEvent& esd,
1286                                         TVector3& ip)  const
1287 {
1288   const AliESDVertex *vertex = esd.GetPrimaryVertexSPD();
1289   if (!vertex) return kNoSPDVtx;
1290   if (vertex->GetNContributors() <= 0) return kFewContrib;
1291
1292   TString vtxTyp = vertex->GetTitle();
1293   if (vtxTyp.Contains("vertexer: Z")) return kNotVtxZ;
1294
1295   if (vertex->GetDispersion() >= 0.04 || vertex->GetZRes()>=0.25)
1296     return kUncertain;
1297
1298   ip.SetX(vertex->GetX());
1299   ip.SetY(vertex->GetY());
1300   ip.SetZ(vertex->GetZ());
1301
1302   return kVtxOK;
1303 }
1304
1305 //____________________________________________________________________
1306 AliFMDEventInspector::EVtxStatus
1307 AliFMDEventInspector::CheckVertex(const AliESDEvent& esd,
1308                                   TVector3& ip) const
1309 {
1310   // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1311   // Get the vertex
1312   const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1313   if (!vertex) {
1314     if (fDebug > 2) {
1315       AliWarning("No SPD vertex found in ESD"); }
1316     return kNoSPDVtx;
1317   }
1318
1319   // #if 0 // Check disabled - seem to kill a lot of PbPb events
1320   // Check that enough tracklets contributed
1321   if(vertex->GetNContributors() <= 0) {
1322     DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1323          vertex->GetNContributors());
1324     ip.SetZ(0);
1325     return kFewContrib;
1326   }
1327   // #endif
1328
1329   // Check that the uncertainty isn't too large
1330   if (vertex->GetZRes() > fMaxVzErr) {
1331     DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f",
1332          vertex->GetZRes(), fMaxVzErr);
1333     return kUncertain;
1334   }
1335
1336   // Get the z coordiante
1337   ip.SetZ(vertex->GetZ());
1338   const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1339
1340
1341   if(!vertexXY->IsFromVertexerZ()) {
1342     ip.SetX(vertexXY->GetX());
1343     ip.SetY(vertexXY->GetY());
1344   }
1345   return kVtxOK;
1346 }
1347
1348 //____________________________________________________________________
1349 Bool_t
1350 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
1351 {
1352   //
1353   // Read the collision system, collision energy, and L3 field setting
1354   // from the ESD
1355   //
1356   // Parameters:
1357   //   esd ESD to get information from
1358   //
1359   // Return:
1360   //    true on success, false
1361   //
1362   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1363   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(),
1364   //           esd->GetMagneticField()));
1365   DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1366   const char* sys  = esd->GetBeamType();
1367   Float_t     cms  = 2 * esd->GetBeamEnergy();
1368   Float_t     fld  = esd->GetMagneticField();
1369   fCollisionSystem = AliForwardUtil::ParseCollisionSystem(sys);
1370   fEnergy          = AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,
1371                                                              cms);
1372   fField           = AliForwardUtil::ParseMagneticField(fld);
1373   fRunNumber       = esd->GetRunNumber();
1374   StoreProduction();
1375   StoreInformation();
1376   if (fCollisionSystem   == AliForwardUtil::kUnknown) {
1377     AliWarningF("Unknown collision system: %s - please check", sys);
1378     return false;
1379   }
1380   if (fEnergy            <= 0) {
1381     AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1382     return false;
1383   }
1384   if (TMath::Abs(fField) >  10) {
1385     AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1386     return false;
1387   }
1388
1389   return true;
1390 }
1391
1392
1393 //____________________________________________________________________
1394 const Char_t*
1395 AliFMDEventInspector::CodeString(UInt_t code)
1396 {
1397   static TString s;
1398   s = "";
1399   if (code & kNoEvent)    s.Append("NOEVENT ");
1400   if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1401   if (code & kNoSPD)      s.Append("NOSPD ");
1402   if (code & kNoFMD)      s.Append("NOFMD ");
1403   if (code & kNoVertex)   s.Append("NOVERTEX ");
1404   if (code & kBadVertex)  s.Append("BADVERTEX ");
1405   return s.Data();
1406 }
1407 #define PF(N,V,...)                                     \
1408   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1409 #define PFB(N,FLAG)                             \
1410   do {                                                                  \
1411     AliForwardUtil::PrintName(N);                                       \
1412     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1413   } while(false)
1414 #define PFV(N,VALUE)                                    \
1415   do {                                                  \
1416     AliForwardUtil::PrintName(N);                       \
1417     std::cout << (VALUE) << std::endl; } while(false)
1418
1419 //____________________________________________________________________
1420 void
1421 AliFMDEventInspector::Print(Option_t*) const
1422 {
1423   //
1424   // Print information
1425   //
1426   //   option Not used
1427   //
1428   AliForwardUtil::PrintTask(*this);
1429   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
1430   sNN.Strip(TString::kBoth, '0');
1431   sNN.ReplaceAll("GeV", " GeV");
1432   TString field(AliForwardUtil::MagneticFieldString(fField));
1433   field.ReplaceAll("p",  "+");
1434   field.ReplaceAll("m",  "-");
1435   field.ReplaceAll("kG", " kG");
1436
1437   gROOT->IncreaseDirLevel();
1438   PFV("Production", "");
1439   gROOT->IncreaseDirLevel();
1440   PF("Period", "LHC%02d%c", fProdYear, fProdLetter);
1441   PFB("MC anchor", fProdMC);
1442   PFV("AliROOT Revision", fProdSVN);
1443   gROOT->DecreaseDirLevel();
1444   PFV("Vertex bins", fVtxAxis.GetNbins());
1445   PF("Vertex Range", "[%+6.1f,%+6.1f]", fVtxAxis.GetXmin(), fVtxAxis.GetXmax());
1446   PFV("Low flux cut",           fLowFluxCut );
1447   PFV("Max(delta v_z)",         fMaxVzErr );
1448   PFV("Min(nContrib_pileup)",   fMinPileupContrib );
1449   PFV("Min(v-pileup)",          fMinPileupDistance );
1450   PFV("System", AliForwardUtil::CollisionSystemString(fCollisionSystem));
1451   PFV("CMS energy per nucleon", sNN);
1452   PFV("Field",                  field);
1453   PFB("Satellite events",       fUseDisplacedVertices);
1454   PFB("Use 2012 pA vertex",     fUsepA2012Vertex );
1455   PFB("Use PWG-UD vertex",      fUseFirstPhysicsVertex);
1456   PFB("Simulation input",       fMC );
1457   PFV("Centrality method",      fCentMethod);
1458   PFV("Centrality axis",        (!fCentAxis ? "none" : ""));
1459   if (!fCentAxis) {return; }
1460   Int_t nBin = fCentAxis->GetNbins();
1461   for (Int_t i = 0; i < nBin; i++) {
1462     if (i != 0 && (i % 10) == 0) std::cout << '\n';
1463     if ((i % 10) == 0)           std::cout << "    ";
1464     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1465   }
1466   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1467   gROOT->DecreaseDirLevel();
1468 }
1469
1470
1471 //
1472 // EOF
1473 //
1474