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