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