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 <iostream>
  36 #include <iomanip>
  37 #include "AliMCEvent.h"
  38 #include "AliHeader.h"
  39 #include "AliGenEventHeader.h"
  40 #include "AliCollisionGeometry.h"
  41 #include "AliVVZERO.h"
  42
  43 //====================================================================
  44 AliFMDEventInspector::AliFMDEventInspector()
  45   : TNamed(),
  46     fHEventsTr(0),
  47     fHEventsTrVtx(0),
  48     fHEventsAccepted(0),
  49     fHEventsAcceptedXY(0),
  50     fHTriggers(0),
  51     fHType(0),
  52   fHWords(0),
  53     fHCent(0),
  54   fHCentVsQual(0),
  55   fHStatus(0),
  56     fLowFluxCut(1000),
  57     fMaxVzErr(0.2),
  58     fList(0),
  59     fEnergy(0),
  60     fField(999),
  61     fCollisionSystem(kUnknown),
  62     fDebug(0),
  63     fCentAxis(0),
  64     fVtxAxis(10,-10,10),
  65     fUseFirstPhysicsVertex(true),
  66     fUseV0AND(false),
  67     fMinPileupContrib(3),
  68     fMinPileupDistance(0.8),
  69     fUseDisplacedVertices(false),
  70   fDisplacedVertex(),
  71   fCollWords(),
  72   fBgWords()
  73 {
  74   //
  75   // Constructor
  76   //
  77   DGUARD(fDebug,1,"Default CTOR of AliFMDEventInspector");
  78 }
  79
  80 //____________________________________________________________________
  81 AliFMDEventInspector::AliFMDEventInspector(const char* name)
  82   : TNamed("fmdEventInspector", name),
  83     fHEventsTr(0),
  84     fHEventsTrVtx(0),
  85     fHEventsAccepted(0),
  86     fHEventsAcceptedXY(0),
  87     fHTriggers(0),
  88     fHType(0),
  89     fHWords(0),
  90     fHCent(0),
  91     fHCentVsQual(0),
  92     fHStatus(0),
  93     fLowFluxCut(1000),
  94     fMaxVzErr(0.2),
  95     fList(0),
  96     fEnergy(0),
  97     fField(999),
  98     fCollisionSystem(kUnknown),
  99     fDebug(0),
 100     fCentAxis(0),
 101     fVtxAxis(10,-10,10),
 102     fUseFirstPhysicsVertex(true),
 103     fUseV0AND(false),
 104     fMinPileupContrib(3),
 105     fMinPileupDistance(0.8),
 106     fUseDisplacedVertices(false),
 107   fDisplacedVertex(),
 108   fCollWords(),
 109   fBgWords()
 110 {
 111   //
 112   // Constructor
 113   //
 114   // Parameters:
 115   //   name Name of object
 116   //
 117   DGUARD(fDebug,1,"Named CTOR of AliFMDEventInspector: %s", name);
 118 }
 119
 120 //____________________________________________________________________
 121 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
 122   : TNamed(o),
 123     fHEventsTr(o.fHEventsTr),
 124     fHEventsTrVtx(o.fHEventsTrVtx),
 125     fHEventsAccepted(o.fHEventsAccepted),
 126     fHEventsAcceptedXY(o.fHEventsAcceptedXY),
 127     fHTriggers(o.fHTriggers),
 128     fHType(o.fHType),
 129     fHWords(o.fHWords),
 130     fHCent(o.fHCent),
 131     fHCentVsQual(o.fHCentVsQual),
 132     fHStatus(o.fHStatus),
 133     fLowFluxCut(o.fLowFluxCut),
 134     fMaxVzErr(o.fMaxVzErr),
 135     fList(o.fList),
 136     fEnergy(o.fEnergy),
 137     fField(o.fField),
 138     fCollisionSystem(o.fCollisionSystem),
 139     fDebug(0),
 140     fCentAxis(0),
 141     fVtxAxis(o.fVtxAxis),
 142     fUseFirstPhysicsVertex(o.fUseFirstPhysicsVertex),
 143     fUseV0AND(o.fUseV0AND),
 144     fMinPileupContrib(o.fMinPileupContrib),
 145     fMinPileupDistance(o.fMinPileupDistance),
 146     fUseDisplacedVertices(o.fUseDisplacedVertices),
 147     fDisplacedVertex(o.fDisplacedVertex),
 148   fCollWords(),
 149   fBgWords()
 150 {
 151   //
 152   // Copy constructor
 153   //
 154   // Parameters:
 155   //   o Object to copy from
 156   //
 157   DGUARD(fDebug,1,"Copy CTOR of AliFMDEventInspector");
 158 }
 159
 160 //____________________________________________________________________
 161 AliFMDEventInspector::~AliFMDEventInspector()
 162 {
 163   //
 164   // Destructor
 165   //
 166   DGUARD(fDebug,1,"DTOR of AliFMDEventInspector");
 167   // if (fList)         delete fList;
 168 }
 169 //____________________________________________________________________
 170 AliFMDEventInspector&
 171 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
 172 {
 173   //
 174   // Assignement operator
 175   //
 176   // Parameters:
 177   //   o Object to assign from
 178   //
 179   // Return:
 180   //    Reference to this object
 181   //
 182   DGUARD(fDebug,3,"Assignment of AliFMDEventInspector");
 183   if (&o == this) return *this;
 184   TNamed::operator=(o);
 185   fHEventsTr         = o.fHEventsTr;
 186   fHEventsTrVtx      = o.fHEventsTrVtx;
 187   fHEventsAccepted   = o.fHEventsAccepted;
 188   fHEventsAcceptedXY = o.fHEventsAcceptedXY;
 189   fHTriggers         = o.fHTriggers;
 190   fHType             = o.fHType;
 191   fHWords            = o.fHWords;
 192   fHCent             = o.fHCent;
 193   fHCentVsQual       = o.fHCentVsQual;
 194   fHStatus           = o.fHStatus;
 195   fLowFluxCut        = o.fLowFluxCut;
 196   fMaxVzErr          = o.fMaxVzErr;
 197   fDebug             = o.fDebug;
 198   fList              = (o.fList ? new TList : 0);
 199   fEnergy            = o.fEnergy;
 200   fField             = o.fField;
 201   fCollisionSystem   = o.fCollisionSystem;
 202   fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(),
 203                o.fVtxAxis.GetXmax());
 204
 205   fUseFirstPhysicsVertex = o.fUseFirstPhysicsVertex;
 206   fUseV0AND              = o.fUseV0AND;
 207   fMinPileupContrib      = o.fMinPileupContrib;
 208   fMinPileupDistance     = o.fMinPileupDistance;
 209   fUseDisplacedVertices  = o.fUseDisplacedVertices;
 210   fDisplacedVertex       = o.fDisplacedVertex;
 211   if (fList) {
 212     fList->SetName(GetName());
 213     if (fHEventsTr)    fList->Add(fHEventsTr);
 214     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
 215     if (fHTriggers)    fList->Add(fHTriggers);
 216     if (fHType)        fList->Add(fHType);
 217     if (fHWords)       fList->Add(fHWords);
 218     if (fHCent)        fList->Add(fHCent);
 219     if (fHCentVsQual)  fList->Add(fHCentVsQual);
 220     if (fHStatus)      fList->Add(fHStatus);
 221   }
 222   return *this;
 223 }
 224
 225 //____________________________________________________________________
 226 Bool_t
 227 AliFMDEventInspector::FetchHistograms(const TList* d,
 228                                       TH1I*& hEventsTr,
 229                                       TH1I*& hEventsTrVtx,
 230                                       TH1I*& hTriggers) const
 231 {
 232   //
 233   // Fetch our histograms from the passed list
 234   //
 235   // Parameters:
 236   //   d             Input
 237   //   hEventsTr     On return, pointer to histogram, or null
 238   //   hEventsTrVtx  On return, pointer to histogram, or null
 239   //   hTriggers     On return, pointer to histogram, or null
 240   //
 241   // Return:
 242   //    true on success, false otherwise
 243   //
 244   DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
 245   hEventsTr    = 0;
 246   hEventsTrVtx = 0;
 247   hTriggers    = 0;
 248   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
 249   if (!dd) return kFALSE;
 250
 251   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
 252   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
 253   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
 254
 255   if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
 256   return kTRUE;
 257 }
 258 //____________________________________________________________________
 259 void
 260 AliFMDEventInspector::CacheConfiguredTriggerClasses(TList& cache,
 261                                                     const TList* classes,
 262                                                     AliOADBPhysicsSelection* o)
 263 {
 264   TIter nextClass(classes);
 265   TObjString* trigClass = 0;
 266   // Loop over all trigger classes.  Trigger classes have the format
 267   //
 268   //   class          := positive_words SPACE(s) negative_words
 269   //   positive_words :=
 270   //                  |  '+' words
 271   //   negative_words :=
 272   //                  |  '-' words
 273   //   words          := word
 274   //                  |  word ',' words
 275   //
 276   while ((trigClass = static_cast<TObjString*>(nextClass()))) {
 277     // Tokenize on space to get positive and negative parts
 278     TString     side   = o->GetBeamSide(trigClass->String());
 279     TObjArray*  parts  = trigClass->String().Tokenize(" ");
 280     TObjString* part   = 0;
 281     TIter       nextPart(parts);
 282     while ((part = static_cast<TObjString*>(nextPart()))) {
 283       // We only care about the positive ones
 284       if (part->GetName()[0] != '+') continue;
 285       part->String().Remove(0,1);
 286
 287       // Tokenize on a comma to get the words
 288       TObjArray*  words = part->String().Tokenize(",");
 289       TObjString* word  = 0;
 290       TIter       nextWord(words);
 291       while ((word = static_cast<TObjString*>(nextWord()))) {
 292         TNamed* store = new TNamed(word->String(), side);
 293         cache.Add(store);
 294         DMSG(fDebug,3,"Caching %s trigger word %s",
 295              store->GetTitle(), store->GetName());
 296       } // while (word)
 297       delete words;
 298     }
 299     delete parts;
 300   }
 301 }
 302
 303 //____________________________________________________________________
 304 void
 305 AliFMDEventInspector::Init(const TAxis& vtxAxis)
 306 {
 307   //
 308   // Initialize the object - this is called on the first seen event.
 309   //
 310   // Parameters:
 311   //   vtxAxis Vertex axis in use
 312   //
 313   DGUARD(fDebug,1,"Initialize in AliFMDEventInspector");
 314
 315   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 316
 317   // Get the input handler - should always be there
 318   AliInputEventHandler* ih =
 319     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 320   if (!ih) {
 321     AliWarning("No input handler");
 322     return;
 323   }
 324   // Get the physics selection - should always be there
 325   AliPhysicsSelection* ps =
 326     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
 327   if (!ps) {
 328     AliWarning("No physics selection");
 329     return;
 330   }
 331   // Get the configured triggers
 332   AliOADBPhysicsSelection* oadb =
 333     const_cast<AliOADBPhysicsSelection*>(ps->GetOADBPhysicsSelection());
 334   if (!oadb) {
 335     AliWarning("No OADB physics selection object");
 336     return;
 337   }
 338   // Get the configured trigger words from the physics selection
 339   const TList* collTriggClasses = ps->GetCollisionTriggerClasses();
 340   const TList* bgTriggClasses   = ps->GetBGTriggerClasses();
 341   if (!collTriggClasses) {
 342     AliWarning("No configured collision trigger classes");
 343     return;
 344   }
 345   if (!bgTriggClasses) {
 346     AliWarning("No configured background trigger classes");
 347     return;
 348   }
 349   CacheConfiguredTriggerClasses(fCollWords, collTriggClasses, oadb);
 350   CacheConfiguredTriggerClasses(fBgWords,   bgTriggClasses,   oadb);
 351   // fCollWords.ls();
 352   // fBgWords.ls();
 353
 354
 355   // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
 356   // ----- 92 number --------- ---- 1 ---
 357   TArrayD limits(93);
 358   for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
 359   limits[92] = 100.5;
 360
 361   fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
 362
 363   fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
 364   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
 365                         4*vtxAxis.GetNbins(),
 366                         2*vtxAxis.GetXmin(),
 367                         2*vtxAxis.GetXmax());
 368   fHEventsTr->SetXTitle("v_{z} [cm]");
 369   fHEventsTr->SetYTitle("# of events");
 370   fHEventsTr->SetFillColor(kRed+1);
 371   fHEventsTr->SetFillStyle(3001);
 372   fHEventsTr->SetDirectory(0);
 373   // fHEventsTr->Sumw2();
 374   fList->Add(fHEventsTr);
 375
 376   fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
 377   fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex");
 378   fHEventsTrVtx->SetFillColor(kBlue+1);
 379   fHEventsTrVtx->SetDirectory(0);
 380   // fHEventsTrVtx->Sumw2();
 381   fList->Add(fHEventsTrVtx);
 382
 383   fHEventsAccepted = new TH1I("nEventsAccepted",
 384                               "Number of events  w/trigger and vertex in range",
 385                               2*vtxAxis.GetNbins(),
 386                               2*vtxAxis.GetXmin(),
 387                               2*vtxAxis.GetXmax());
 388   fHEventsAccepted->SetXTitle("v_{z} [cm]");
 389   fHEventsAccepted->SetYTitle("# of events");
 390   fHEventsAccepted->SetFillColor(kGreen+1);
 391   fHEventsAccepted->SetFillStyle(3001);
 392   fHEventsAccepted->SetDirectory(0);
 393   // fHEventsAccepted->Sumw2();
 394   fList->Add(fHEventsAccepted);
 395
 396   fHEventsAcceptedXY = new TH2D("nEventsAcceptedXY",
 397                               "XY vertex w/trigger and Z vertex in range",
 398                                 1000,-1,1,1000,-1,1);
 399
 400   fHEventsAcceptedXY->SetXTitle("v_{x} [cm]");
 401   fHEventsAcceptedXY->SetYTitle("v_{y} [cm]");
 402   fHEventsAcceptedXY->SetDirectory(0);
 403   // fHEventsAccepted->Sumw2();
 404   fList->Add(fHEventsAcceptedXY);
 405
 406
 407   fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
 408   fHTriggers->SetFillColor(kRed+1);
 409   fHTriggers->SetFillStyle(3001);
 410   fHTriggers->SetStats(0);
 411   fHTriggers->SetDirectory(0);
 412   fHTriggers->GetXaxis()->SetBinLabel(kInel   +1,"INEL");
 413   fHTriggers->GetXaxis()->SetBinLabel(kInelGt0+1,"INEL>0");
 414   fHTriggers->GetXaxis()->SetBinLabel(kNSD    +1,"NSD");
 415   fHTriggers->GetXaxis()->SetBinLabel(kV0AND  +1,"VOAND");
 416   fHTriggers->GetXaxis()->SetBinLabel(kEmpty  +1,"Empty");
 417   fHTriggers->GetXaxis()->SetBinLabel(kA      +1,"A");
 418   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
 419   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
 420   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
 421   fHTriggers->GetXaxis()->SetBinLabel(kPileUp +1,"Pileup");
 422   fHTriggers->GetXaxis()->SetBinLabel(kMCNSD  +1,"NSD_{MC}");
 423   fHTriggers->GetXaxis()->SetBinLabel(kOffline+1,"Offline");
 424   fList->Add(fHTriggers);
 425
 426   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)",
 427                                  fLowFluxCut), 2, -.5, 1.5);
 428   fHType->SetFillColor(kRed+1);
 429   fHType->SetFillStyle(3001);
 430   fHType->SetStats(0);
 431   fHType->SetDirectory(0);
 432   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
 433   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
 434   fList->Add(fHType);
 435
 436 #if 0
 437   // This histogram disabled as it causes problems in the merge
 438   fHWords = new TH1I("words", "Trigger words seen", 1, 0, 0);
 439   fHWords->SetFillColor(kBlue+1);
 440   fHWords->SetFillStyle(3001);
 441   fHWords->SetStats(0);
 442   fHWords->SetDirectory(0);
 443   fHWords->SetBit(TH1::kCanRebin);
 444   fList->Add(fHWords);
 445 #endif
 446
 447   fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
 448   fHCent->SetFillColor(kBlue+1);
 449   fHCent->SetFillStyle(3001);
 450   fHCent->SetStats(0);
 451   fHCent->SetDirectory(0);
 452   fHCent->SetXTitle("Centrality [%]");
 453   fHCent->SetYTitle("Events");
 454   fList->Add(fHCent);
 455
 456   fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality",
 457                           5, 0, 5, limits.GetSize()-1, limits.GetArray());
 458   fHCentVsQual->SetXTitle("Quality");
 459   fHCentVsQual->SetYTitle("Centrality [%]");
 460   fHCentVsQual->SetZTitle("Events");
 461   fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
 462   fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
 463   fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
 464   fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
 465   fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
 466   fList->Add(fHCentVsQual);
 467
 468   fHStatus = new TH1I("status", "Status", 7, 1, 8);
 469   fHStatus->SetFillColor(kRed+1);
 470   fHStatus->SetFillStyle(3001);
 471   fHStatus->SetStats(0);
 472   fHStatus->SetDirectory(0);
 473   fHStatus->GetXaxis()->SetBinLabel(1, "OK");
 474   fHStatus->GetXaxis()->SetBinLabel(2, "No event");
 475   fHStatus->GetXaxis()->SetBinLabel(3, "No triggers");
 476   fHStatus->GetXaxis()->SetBinLabel(4, "No SPD");
 477   fHStatus->GetXaxis()->SetBinLabel(5, "No FMD");
 478   fHStatus->GetXaxis()->SetBinLabel(6, "No vertex");
 479   fHStatus->GetXaxis()->SetBinLabel(7, "Bad vertex");
 480   fList->Add(fHStatus);
 481 }
 482
 483 //____________________________________________________________________
 484 void
 485 AliFMDEventInspector::StoreInformation(Int_t runNo)
 486 {
 487   // Write TNamed objects to output list containing information about
 488   // the running conditions
 489   DGUARD(fDebug,2,"Store information from AliFMDEventInspector");
 490   if (!fList) return;
 491
 492
 493   fList->Add(AliForwardUtil::MakeParameter("sys", fCollisionSystem));
 494   fList->Add(AliForwardUtil::MakeParameter("sNN", fEnergy));
 495   fList->Add(AliForwardUtil::MakeParameter("field", fField));
 496   fList->Add(AliForwardUtil::MakeParameter("runNo", runNo));
 497   fList->Add(AliForwardUtil::MakeParameter("lowFlux", fLowFluxCut));
 498   fList->Add(AliForwardUtil::MakeParameter("fpVtx",fUseFirstPhysicsVertex));
 499   fList->Add(AliForwardUtil::MakeParameter("v0and",fUseV0AND));
 500   fList->Add(AliForwardUtil::MakeParameter("nPileUp", fMinPileupContrib));
 501   fList->Add(AliForwardUtil::MakeParameter("dPileup", fMinPileupDistance));
 502 }
 503
 504 //____________________________________________________________________
 505 void
 506 AliFMDEventInspector::DefineOutput(TList* dir)
 507 {
 508   //
 509   // Define the output histograms.  These are put in a sub list of the
 510   // passed list.   The histograms are merged before the parent task calls
 511   // AliAnalysisTaskSE::Terminate
 512   //
 513   //   dir Directory to add to
 514   //
 515   DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
 516   fList = new TList;
 517   fList->SetName(GetName());
 518   fList->SetOwner();
 519   dir->Add(fList);
 520 }
 521
 522 //____________________________________________________________________
 523 UInt_t
 524 AliFMDEventInspector::Process(const AliESDEvent* event,
 525                               UInt_t&            triggers,
 526                               Bool_t&            lowFlux,
 527                               UShort_t&          ivz,
 528                               Double_t&          vz,
 529                               Double_t&          cent,
 530                               UShort_t&          nClusters)
 531 {
 532   //
 533   // Process the event
 534   //
 535   // Parameters:
 536   //   event     Input event
 537   //   triggers  On return, the triggers fired
 538   //   lowFlux   On return, true if the event is considered a low-flux
 539   //                  event (according to the setting of fLowFluxCut)
 540   //   ivz       On return, the found vertex bin (1-based).  A zero
 541   //                  means outside of the defined vertex range
 542   //   vz        On return, the z position of the interaction
 543   //   cent      On return, the centrality - if not available < 0
 544   //
 545   // Return:
 546   //    0 (or kOk) on success, otherwise a bit mask of error codes
 547   //
 548   DGUARD(fDebug,1,"Process event in AliFMDEventInspector");
 549
 550   // --- Check that we have an event ---------------------------------
 551   if (!event) {
 552     AliWarning("No ESD event found for input event");
 553     fHStatus->Fill(2);
 554     return kNoEvent;
 555   }
 556
 557   // --- Read trigger information from the ESD and store in AOD object
 558   if (!ReadTriggers(*event, triggers, nClusters)) {
 559     if (fDebug > 2) {
 560       AliWarning("Failed to read triggers from ESD"); }
 561     fHStatus->Fill(3);
 562     return kNoTriggers;
 563   }
 564
 565   // --- Check if this is a high-flux event --------------------------
 566   const AliMultiplicity* testmult = event->GetMultiplicity();
 567   if (!testmult) {
 568     if (fDebug > 3) {
 569       AliWarning("No central multiplicity object found"); }
 570   }
 571   else
 572     lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 573
 574   fHType->Fill(lowFlux ? 0 : 1);
 575
 576   // --- Process satellite event information is requested ------------
 577   if (fUseDisplacedVertices) {
 578     if (!fDisplacedVertex.Process(event))
 579       AliWarning("Failed to process satellite event");
 580   }
 581
 582   // --- Read centrality information
 583   cent          = -10;
 584   UShort_t qual = 0;
 585   if (!ReadCentrality(*event, cent, qual)) {
 586     if (fDebug > 3)
 587       AliWarning("Failed to get centrality");
 588   }
 589   fHCent->Fill(cent);
 590   if (qual == 0) fHCentVsQual->Fill(0., cent);
 591   else {
 592     for (UShort_t i = 0; i < 4; i++)
 593       if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
 594   }
 595
 596   // --- Get the vertex information ----------------------------------
 597   Double_t vx = 0;
 598   Double_t vy = 0;
 599   vz          = 0;
 600
 601   Bool_t vzOk = ReadVertex(*event, vz,vx,vy);
 602
 603   fHEventsTr->Fill(vz);
 604   if (!vzOk) {
 605     if (fDebug > 3) {
 606       AliWarning("Failed to read vertex from ESD"); }
 607     fHStatus->Fill(6);
 608     return kNoVertex;
 609   }
 610   fHEventsTrVtx->Fill(vz);
 611
 612   // --- Get the vertex bin ------------------------------------------
 613   ivz = fVtxAxis.FindBin(vz);
 614   if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
 615     if (fDebug > 3) {
 616       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 617                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
 618     }
 619     ivz = 0;
 620     fHStatus->Fill(7);
 621     return kBadVertex;
 622   }
 623   fHEventsAccepted->Fill(vz);
 624   fHEventsAcceptedXY->Fill(vx,vy);
 625
 626   // --- Check the FMD ESD data --------------------------------------
 627   if (!event->GetFMDData()) {
 628     if (fDebug > 3) {
 629       AliWarning("No FMD data found in ESD"); }
 630     fHStatus->Fill(5);
 631     return kNoFMD;
 632   }
 633
 634   fHStatus->Fill(1);
 635   return kOk;
 636 }
 637
 638 //____________________________________________________________________
 639 Bool_t
 640 AliFMDEventInspector::ReadCentrality(const AliESDEvent& esd,
 641                                      Double_t& cent,
 642                                      UShort_t& qual) const
 643 {
 644   //
 645   // Read centrality from event
 646   //
 647   // Parameters:
 648   //    esd  Event
 649   //    cent On return, the centrality or negative if not found
 650   //
 651   // Return:
 652   //    False on error, true otherwise
 653   //
 654   DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
 655
 656   if(fUseDisplacedVertices) {
 657     Double_t zvtx = fDisplacedVertex.GetVertexZ();
 658     qual          = 1;
 659     if(TMath::Abs(zvtx) < 999) {
 660       cent = fDisplacedVertex.GetCentralityPercentile();
 661       qual = 0;
 662     }
 663     return true;
 664   }
 665
 666   cent = -1;
 667   qual = 0;
 668   AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
 669   if (!centObj)  return true;
 670
 671   cent = centObj->GetCentralityPercentile("V0M");
 672   qual = centObj->GetQuality();
 673
 674   return true;
 675 }
 676
 677 //____________________________________________________________________
 678 Bool_t
 679 AliFMDEventInspector::CheckpAExtraV0(const AliESDEvent& esd) const
 680 {
 681   if (fCollisionSystem != AliForwardUtil::kPPb) return true;
 682
 683    AliVVZERO* esdV0 = esd.GetVZEROData();
 684    if ((esdV0->GetV0ADecision()!=1) || (esdV0->GetV0CDecision()!=1))
 685      return false;
 686    return true;
 687 }
 688
 689 //____________________________________________________________________
 690 Bool_t
 691 AliFMDEventInspector::ReadTriggers(const AliESDEvent& esd, UInt_t& triggers,
 692                                    UShort_t& nClusters)
 693 {
 694   //
 695   // Read the trigger information from the ESD event
 696   //
 697   // Parameters:
 698   //   esd        ESD event
 699   //   triggers   On return, contains the trigger bits
 700   //
 701   // Return:
 702   //    @c true on success, @c false otherwise
 703   //
 704   DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
 705   triggers = 0;
 706
 707   // Get the analysis manager - should always be there
 708   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 709   DMSG(fDebug,10,"Got analysis manager %p", am);
 710   if (!am) {
 711     AliWarning("No analysis manager defined!");
 712     return kFALSE;
 713   }
 714
 715   // Get the input handler - should always be there
 716   AliInputEventHandler* ih =
 717     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 718   DMSG(fDebug,10,"Got input handler %p", ih);
 719   if (!ih) {
 720     AliWarning("No input handler");
 721     return kFALSE;
 722   }
 723
 724   // Check if this is a collision candidate (MB)
 725   // Note, that we should use the value cached in the input
 726   // handler rather than calling IsCollisionCandiate directly
 727   // on the AliPhysicsSelection obejct.  If we called the latter
 728   // then the AliPhysicsSelection object would overcount by a
 729   // factor of 2! :-(
 730   Bool_t  offline  = ih->IsEventSelected();
 731   Bool_t  fastonly = (ih->IsEventSelected() & AliVEvent::kFastOnly);
 732   TString trigStr  = esd.GetFiredTriggerClasses();
 733
 734   if (fHWords) fHWords->Fill(trigStr.Data(), 1);
 735
 736   if(fUseDisplacedVertices) {
 737     DMSG(fDebug,3,"Using displaced vertex stuff");
 738     if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
 739   }
 740
 741   if (CheckFastPartition(fastonly))     offline = false;
 742   if (offline && CheckCosmics(trigStr)) offline = false;
 743   // if (!CheckpAExtraV0(esd))             offline = false;
 744
 745   DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
 746
 747   if (offline) {
 748     triggers |= AliAODForwardMult::kOffline;
 749     triggers |= AliAODForwardMult::kInel;
 750     if (!fHTriggers) {
 751       AliWarning("Histogram of triggers not defined - has init been called");
 752       return false;
 753     }
 754     fHTriggers->Fill(kOffline+0.5);
 755
 756     CheckINELGT0(esd, nClusters, triggers);
 757   }
 758
 759   CheckNSD(esd,triggers);
 760   if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
 761   if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
 762
 763   CheckWords(esd, triggers);
 764
 765   // Now check - if we have a collision - for offline triggers and
 766   // fill histogram.
 767   if (triggers & AliAODForwardMult::kB) {
 768     fHTriggers->Fill(kB+.5);
 769     if (triggers & AliAODForwardMult::kInel)
 770       fHTriggers->Fill(kInel);
 771
 772     if (triggers & AliAODForwardMult::kInelGt0)
 773       fHTriggers->Fill(kInelGt0+.5);
 774
 775     if (triggers & AliAODForwardMult::kNSD)
 776       fHTriggers->Fill(kNSD+.5);
 777
 778     if (triggers & AliAODForwardMult::kV0AND)
 779       fHTriggers->Fill(kV0AND+.5);
 780   }
 781   if (triggers & AliAODForwardMult::kA) fHTriggers->Fill(kA+.5);
 782   if (triggers & AliAODForwardMult::kC) fHTriggers->Fill(kC+.5);
 783   if (triggers & AliAODForwardMult::kE) fHTriggers->Fill(kE+.5);
 784
 785   return kTRUE;
 786 }
 787
 788 //____________________________________________________________________
 789 Bool_t
 790 AliFMDEventInspector::CheckFastPartition(bool fastonly) const
 791 {
 792   // For the 2.76 TeV p+p run, the FMD ran in the slow partition
 793   // so it received no triggers from the fast partition. Therefore
 794   // the fast triggers are removed here but not for MC where all
 795   // triggers are fast.
 796   if (TMath::Abs(fEnergy - 2750.) > 20) return false;
 797   if (fCollisionSystem != AliForwardUtil::kPP) return false;
 798   if (fastonly)
 799     DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
 800
 801   return fastonly;
 802 }
 803
 804 //____________________________________________________________________
 805 Bool_t
 806 AliFMDEventInspector::CheckCosmics(const TString& trigStr) const
 807 {
 808   // MUON triggers are not strictly minimum bias (MB) so they are
 809   // removed (HHD)
 810   if(trigStr.Contains("CMUS1")) {
 811     DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
 812     return true;
 813   }
 814   return false;
 815 }
 816
 817 //____________________________________________________________________
 818 Bool_t
 819 AliFMDEventInspector::CheckINELGT0(const AliESDEvent& esd,
 820                                    UShort_t& nClusters,
 821                                    UInt_t& triggers) const
 822 {
 823   nClusters = 0;
 824
 825   // If this is inel, see if we have a tracklet
 826   const AliMultiplicity* spdmult = esd.GetMultiplicity();
 827   if (!spdmult) {
 828     AliWarning("No SPD multiplicity");
 829     return false;
 830   }
 831
 832   // Check if we have one or more tracklets
 833   // in the range -1 < eta < 1 to set the INEL>0
 834   // trigger flag.
 835   //
 836   // Also count tracklets as a single cluster
 837   Int_t n = spdmult->GetNumberOfTracklets();
 838   for (Int_t j = 0; j < n; j++) {
 839     if(TMath::Abs(spdmult->GetEta(j)) < 1) {
 840       triggers |= AliAODForwardMult::kInelGt0;
 841       nClusters++;
 842     }
 843   }
 844   n = spdmult->GetNumberOfSingleClusters();
 845   for (Int_t j = 0; j < n; j++) {
 846     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 847     if (TMath::Abs(eta) < 1) nClusters++;
 848   }
 849   if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
 850
 851   return triggers & AliAODForwardMult::kNClusterGt0;
 852 }
 853
 854 //____________________________________________________________________
 855 Bool_t
 856 AliFMDEventInspector::CheckNSD(const AliESDEvent& esd, UInt_t& triggers) const
 857 {
 858   // Analyse some trigger stuff
 859   AliTriggerAnalysis ta;
 860   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kV0AND)) {
 861     triggers |= AliAODForwardMult::kV0AND;
 862     if (fUseV0AND)
 863       triggers |= AliAODForwardMult::kNSD;
 864   }
 865   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1))
 866     triggers |= AliAODForwardMult::kNSD;
 867   return triggers & AliAODForwardMult::kNSD;
 868 }
 869 //____________________________________________________________________
 870 Bool_t
 871 AliFMDEventInspector::CheckPileup(const AliESDEvent& esd,
 872                                   UInt_t& triggers) const
 873 {
 874   // Check for multiple vertices (pile-up) with at least 3
 875   // contributors and at least 0.8cm from the primary vertex
 876   if(fCollisionSystem != AliForwardUtil::kPP) return false;
 877
 878   Bool_t pileup =  esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance);
 879   if (pileup) triggers |= AliAODForwardMult::kPileUp;
 880   return pileup;
 881 }
 882
 883 //____________________________________________________________________
 884 Bool_t
 885 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
 886 {
 887   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
 888     triggers |= AliAODForwardMult::kEmpty;
 889     return true;
 890   }
 891   return false;
 892 }
 893 //____________________________________________________________________
 894 Bool_t
 895 AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const
 896 {
 897   TObject* word = 0;
 898   TIter nextColl(&fCollWords);
 899   while ((word = nextColl())) {
 900     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
 901          word->GetTitle(), word->GetName());
 902     if (!esd.IsTriggerClassFired(word->GetName())) continue;
 903
 904     TString beamSide = word->GetTitle();
 905     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
 906
 907     if (!beamSide.EqualTo("B")) continue;
 908     triggers |= AliAODForwardMult::kB;
 909     break; // No more to do here
 910   }
 911   TIter nextBg(&fBgWords);
 912   UInt_t all = (AliAODForwardMult::kA |
 913                 AliAODForwardMult::kC |
 914                 AliAODForwardMult::kE);
 915   while ((word = nextBg())) {
 916     DMSG(fDebug,10,"Checking if %s trigger %s is fired",
 917          word->GetTitle(), word->GetName());
 918     if (!esd.IsTriggerClassFired(word->GetName())) continue;
 919
 920     TString beamSide = word->GetTitle();
 921     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
 922
 923     if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
 924     if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
 925     if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
 926
 927     if ((triggers & all) == all) break; // No more to do
 928   }
 929   return true;
 930 }
 931
 932
 933 //____________________________________________________________________
 934 Bool_t
 935 AliFMDEventInspector::ReadVertex(const AliESDEvent& esd,
 936                                  Double_t& vz,
 937                                  Double_t& vx,
 938                                  Double_t& vy)
 939 {
 940   //
 941   // Read the vertex information from the ESD event
 942   //
 943   // Parameters:
 944   //   esd  ESD event
 945   //   vz   On return, the vertex Z position
 946   //
 947   // Return:
 948   //    @c true on success, @c false otherwise
 949   //
 950   DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
 951   vz = 0;
 952   vx = 1024;
 953   vy = 1024;
 954
 955   if(fUseDisplacedVertices) {
 956     Double_t zvtx = fDisplacedVertex.GetVertexZ();
 957
 958     if(TMath::Abs(zvtx) < 999) {
 959       vz = zvtx;
 960       return true;
 961     }
 962     return false;
 963   }
 964
 965   if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, vz, vx, vy);
 966
 967
 968   return CheckVertex(esd, vz, vx, vy);
 969 }
 970
 971 //____________________________________________________________________
 972 Bool_t
 973 AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd,
 974                                        Double_t& vz,
 975                                        Double_t& vx,
 976                                        Double_t& vy) const
 977 {
 978   // This is the code used by the 1st physics people
 979   const AliESDVertex* vertex    = esd.GetPrimaryVertex();
 980   if (!vertex  || !vertex->GetStatus()) {
 981     DMSG(fDebug,2,"No primary vertex (%p) or bad status %d",
 982          vertex, (vertex ? vertex->GetStatus() : -1));
 983     return false;
 984   }
 985   const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
 986   if (!vertexSPD || !vertexSPD->GetStatus()) {
 987     DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d",
 988          vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
 989     return false;
 990   }
 991
 992   // if vertex is from SPD vertexZ, require more stringent cuts
 993   if (vertex->IsFromVertexerZ()) {
 994     if (vertex->GetDispersion() > fMaxVzErr ||
 995         vertex->GetZRes() > 1.25 * fMaxVzErr) {
 996       DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
 997            vertex->GetDispersion(), fMaxVzErr,
 998            vertex->GetZRes(), 1.25 * fMaxVzErr);
 999       return false;
1000     }
1001   }
1002   vz = vertex->GetZ();
1003
1004   if(!vertex->IsFromVertexerZ()) {
1005     vx = vertex->GetX();
1006     vy = vertex->GetY();
1007   }
1008   return true;
1009 }
1010 //____________________________________________________________________
1011 Bool_t
1012 AliFMDEventInspector::CheckVertex(const AliESDEvent& esd,
1013                                   Double_t& vz,
1014                                   Double_t& vx,
1015                                   Double_t& vy) const
1016 {
1017   // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1018   // Get the vertex
1019   const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1020   if (!vertex) {
1021     if (fDebug > 2) {
1022       AliWarning("No SPD vertex found in ESD"); }
1023     return false;
1024   }
1025
1026   // Check that enough tracklets contributed
1027   if(vertex->GetNContributors() <= 0) {
1028     DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1029          vertex->GetNContributors());
1030     vz = 0;
1031     return false;
1032   }
1033   // Check that the uncertainty isn't too large
1034   if (vertex->GetZRes() > fMaxVzErr) {
1035     DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f",
1036          vertex->GetZRes(), fMaxVzErr);
1037     return false;
1038   }
1039
1040   // Get the z coordiante
1041   vz = vertex->GetZ();
1042   const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1043
1044   if(!vertexXY->IsFromVertexerZ()) {
1045     vx = vertexXY->GetX();
1046     vy = vertexXY->GetY();
1047   }
1048   return true;
1049 }
1050
1051 //____________________________________________________________________
1052 Bool_t
1053 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
1054 {
1055   //
1056   // Read the collision system, collision energy, and L3 field setting
1057   // from the ESD
1058   //
1059   // Parameters:
1060   //   esd ESD to get information from
1061   //
1062   // Return:
1063   //    true on success, false
1064   //
1065   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1066   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(),
1067   //           esd->GetMagneticField()));
1068   DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1069   const char* sys  = esd->GetBeamType();
1070   Float_t     cms  = 2 * esd->GetBeamEnergy();
1071   Float_t     fld  = esd->GetMagneticField();
1072   fCollisionSystem = AliForwardUtil::ParseCollisionSystem(sys);
1073   fEnergy          = AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,
1074                                                              cms);
1075   fField           = AliForwardUtil::ParseMagneticField(fld);
1076
1077   StoreInformation(esd->GetRunNumber());
1078   if (fCollisionSystem   == AliForwardUtil::kUnknown) {
1079     AliWarningF("Unknown collision system: %s - please check", sys);
1080     return false;
1081   }
1082   if (fEnergy            <= 0) {
1083     AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1084     return false;
1085   }
1086   if (TMath::Abs(fField) >  10) {
1087     AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1088     return false;
1089   }
1090
1091   return true;
1092 }
1093
1094
1095 //____________________________________________________________________
1096 const Char_t*
1097 AliFMDEventInspector::CodeString(UInt_t code)
1098 {
1099   static TString s;
1100   s = "";
1101   if (code & kNoEvent)    s.Append("NOEVENT ");
1102   if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1103   if (code & kNoSPD)      s.Append("NOSPD ");
1104   if (code & kNoFMD)      s.Append("NOFMD ");
1105   if (code & kNoVertex)   s.Append("NOVERTEX ");
1106   if (code & kBadVertex)  s.Append("BADVERTEX ");
1107   return s.Data();
1108 }
1109 //____________________________________________________________________
1110 void
1111 AliFMDEventInspector::Print(Option_t*) const
1112 {
1113   //
1114   // Print information
1115   //
1116   //   option Not used
1117   //
1118   char ind[gROOT->GetDirLevel()+1];
1119   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1120   ind[gROOT->GetDirLevel()] = '\0';
1121   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
1122   sNN.Strip(TString::kBoth, '0');
1123   sNN.ReplaceAll("GeV", " GeV");
1124   TString field(AliForwardUtil::MagneticFieldString(fField));
1125   field.ReplaceAll("p",  "+");
1126   field.ReplaceAll("m",  "-");
1127   field.ReplaceAll("kG", " kG");
1128
1129   std::cout << std::boolalpha
1130             << ind << ClassName() << ": " << GetName() << '\n'
1131             << ind << " Vertex bins:            " << fVtxAxis.GetNbins() << '\n'
1132             << ind << " Vertex range:           [" << fVtxAxis.GetXmin()
1133             << "," << fVtxAxis.GetXmax() << "]\n"
1134             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
1135             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
1136             << ind << " Min(nContrib_pileup):   " << fMinPileupContrib << '\n'
1137             << ind << " Min(v-pileup):          " << fMinPileupDistance << '\n'
1138             << ind << " System:                 "
1139             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
1140             << ind << " CMS energy per nucleon: " << sNN << '\n'
1141             << ind << " Field:                  " <<  field << '\n'
1142             << ind << " Satellite events:       " << fUseDisplacedVertices
1143             << "\n" << std::noboolalpha;
1144   if (!fCentAxis) { std::cout << std::flush; return; }
1145   Int_t nBin = fCentAxis->GetNbins();
1146   std::cout << ind << " Centrality axis:        " << nBin << " bins"
1147             << std::flush;
1148   for (Int_t i = 0; i < nBin; i++) {
1149     if ((i % 10) == 0) std::cout << '\n' << ind << "  ";
1150     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1151   }
1152   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1153 }
1154
1155
1156 //
1157 // EOF
1158 //
1159