PWGLF/FORWARD/analysis2/AliFMDEventInspector.cxx

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