PWG2/FORWARD/analysis2/AliFMDEventInspector.cxx

   1 #include "AliFMDEventInspector.h"
   2 #include "AliLog.h"
   3 #include "AliESDEvent.h"
   4 #include "AliMultiplicity.h"
   5 #include "AliAnalysisManager.h"
   6 #include "AliInputEventHandler.h"
   7 #include "AliTriggerAnalysis.h"
   8 #include "AliPhysicsSelection.h"
   9 #include "AliAODForwardMult.h"
  10 #include "AliForwardUtil.h"
  11 #include <TH1.h>
  12 #include <TList.h>
  13 #include <TDirectory.h>
  14 #include <TROOT.h>
  15 #include <iostream>
  16 #include <iomanip>
  17
  18 //====================================================================
  19 AliFMDEventInspector::AliFMDEventInspector()
  20   : TNamed(),
  21     fHEventsTr(0),
  22     fHEventsTrVtx(0),
  23     fHTriggers(0),
  24     fHType(0),
  25     fLowFluxCut(1000),
  26     fMaxVzErr(0.1),
  27     fList(0),
  28     fEnergy(0),
  29     fField(999),
  30     fCollisionSystem(kUnknown),
  31     fDebug(0)
  32 {
  33 }
  34
  35 //____________________________________________________________________
  36 AliFMDEventInspector::AliFMDEventInspector(const char* name)
  37   : TNamed("fmdEventInspector", name),
  38     fHEventsTr(0),
  39     fHEventsTrVtx(0),
  40     fHTriggers(0),
  41     fHType(0),
  42     fLowFluxCut(1000),
  43     fMaxVzErr(0.1),
  44     fList(0),
  45     fEnergy(0),
  46     fField(999),
  47     fCollisionSystem(kUnknown),
  48     fDebug(0)
  49 {
  50 }
  51
  52 //____________________________________________________________________
  53 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
  54   : TNamed(o),
  55     fHEventsTr(o.fHEventsTr),
  56     fHEventsTrVtx(o.fHEventsTrVtx),
  57     fHTriggers(o.fHTriggers),
  58     fHType(o.fHType),
  59     fLowFluxCut(o.fMaxVzErr),
  60     fMaxVzErr(o.fMaxVzErr),
  61     fList(o.fList),
  62     fEnergy(o.fEnergy),
  63     fField(o.fField),
  64     fCollisionSystem(o.fCollisionSystem),
  65     fDebug(0)
  66 {
  67 }
  68
  69 //____________________________________________________________________
  70 AliFMDEventInspector::~AliFMDEventInspector()
  71 {
  72   if (fHEventsTr)    delete fHEventsTr;
  73   if (fHEventsTrVtx) delete fHEventsTrVtx;
  74   if (fHTriggers)    delete fHTriggers;
  75   if (fHType)        delete fHType;
  76   if (fList)         delete fList;
  77 }
  78 //____________________________________________________________________
  79 AliFMDEventInspector&
  80 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
  81 {
  82   TNamed::operator=(o);
  83   fHEventsTr         = o.fHEventsTr;
  84   fHEventsTrVtx      = o.fHEventsTrVtx;
  85   fHTriggers         = o.fHTriggers;
  86   fHType             = o.fHType;
  87   fLowFluxCut        = o.fLowFluxCut;
  88   fMaxVzErr          = o.fMaxVzErr;
  89   fDebug             = o.fDebug;
  90   fList              = (o.fList ? new TList : 0);
  91   fEnergy            = o.fEnergy;
  92   fField             = o.fField;
  93   fCollisionSystem   = o.fCollisionSystem;
  94   if (fList) {
  95     fList->SetName(GetName());
  96     if (fHEventsTr)    fList->Add(fHEventsTr);
  97     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
  98     if (fHTriggers)    fList->Add(fHTriggers);
  99     if (fHType)        fList->Add(fHType);
 100   }
 101   return *this;
 102 }
 103
 104 //____________________________________________________________________
 105 Bool_t
 106 AliFMDEventInspector::FetchHistograms(TList* d,
 107                                       TH1I*& hEventsTr,
 108                                       TH1I*& hEventsTrVtx,
 109                                       TH1I*& hTriggers) const
 110 {
 111   hEventsTr    = 0;
 112   hEventsTrVtx = 0;
 113   hTriggers    = 0;
 114   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
 115   if (!dd) return kFALSE;
 116
 117   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
 118   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
 119   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
 120
 121   if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
 122   return kTRUE;
 123 }
 124 //____________________________________________________________________
 125 void
 126 AliFMDEventInspector::Init(const TAxis& vtxAxis)
 127 {
 128   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
 129                         vtxAxis.GetNbins(),
 130                         vtxAxis.GetXmin(),
 131                         vtxAxis.GetXmax());
 132   fHEventsTr->SetXTitle("v_{z} [cm]");
 133   fHEventsTr->SetYTitle("# of events");
 134   fHEventsTr->SetFillColor(kRed+1);
 135   fHEventsTr->SetFillStyle(3001);
 136   fHEventsTr->SetDirectory(0);
 137   // fHEventsTr->Sumw2();
 138   fList->Add(fHEventsTr);
 139
 140   fHEventsTrVtx = new TH1I("nEventsTrVtx",
 141                            "Number of events w/trigger and vertex",
 142                            vtxAxis.GetNbins(),
 143                            vtxAxis.GetXmin(),
 144                            vtxAxis.GetXmax());
 145   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
 146   fHEventsTrVtx->SetYTitle("# of events");
 147   fHEventsTrVtx->SetFillColor(kBlue+1);
 148   fHEventsTrVtx->SetFillStyle(3001);
 149   fHEventsTrVtx->SetDirectory(0);
 150   // fHEventsTrVtx->Sumw2();
 151   fList->Add(fHEventsTrVtx);
 152
 153
 154   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
 155   fHTriggers->SetFillColor(kRed+1);
 156   fHTriggers->SetFillStyle(3001);
 157   fHTriggers->SetStats(0);
 158   fHTriggers->SetDirectory(0);
 159   fHTriggers->GetXaxis()->SetBinLabel(kInel   +1,"INEL");
 160   fHTriggers->GetXaxis()->SetBinLabel(kInelGt0+1,"INEL>0");
 161   fHTriggers->GetXaxis()->SetBinLabel(kNSD    +1,"NSD");
 162   fHTriggers->GetXaxis()->SetBinLabel(kEmpty  +1,"Empty");
 163   fHTriggers->GetXaxis()->SetBinLabel(kA      +1,"A");
 164   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
 165   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
 166   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
 167   fHTriggers->GetXaxis()->SetBinLabel(9,         "spare1");
 168   fHTriggers->GetXaxis()->SetBinLabel(10,        "spare2");
 169   fList->Add(fHTriggers);
 170
 171   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)",
 172                                  fLowFluxCut), 2, -.5, 1.5);
 173   fHType->SetFillColor(kRed+1);
 174   fHType->SetFillStyle(3001);
 175   fHType->SetStats(0);
 176   fHType->SetDirectory(0);
 177   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
 178   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
 179   fList->Add(fHType);
 180
 181 }
 182
 183 //____________________________________________________________________
 184 void
 185 AliFMDEventInspector::DefineOutput(TList* dir)
 186 {
 187   fList = new TList;
 188   fList->SetName(GetName());
 189   dir->Add(fList);
 190 }
 191
 192 //____________________________________________________________________
 193 UInt_t
 194 AliFMDEventInspector::Process(const AliESDEvent* event,
 195                               UInt_t&            triggers,
 196                               Bool_t&            lowFlux,
 197                               UShort_t&          ivz,
 198                               Double_t&          vz)
 199 {
 200   // Check that we have an event
 201   if (!event) {
 202     AliWarning("No ESD event found for input event");
 203     return kNoEvent;
 204   }
 205
 206   // Read trigger information from the ESD and store in AOD object
 207   if (!ReadTriggers(event, triggers)) {
 208     if (fDebug > 2) {
 209       AliWarning("Failed to read triggers from ESD"); }
 210     return kNoTriggers;
 211   }
 212
 213   // Check if this is a high-flux event
 214   const AliMultiplicity* testmult = event->GetMultiplicity();
 215   if (!testmult) {
 216     if (fDebug > 3) {
 217       AliWarning("No central multiplicity object found"); }
 218     return kNoSPD;
 219   }
 220   lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
 221
 222   fHType->Fill(lowFlux ? 0 : 1);
 223
 224   // Check the FMD ESD data
 225   if (!event->GetFMDData()) {
 226     if (fDebug > 3) {
 227       AliWarning("No FMD data found in ESD"); }
 228     return kNoFMD;
 229   }
 230
 231   // Get the vertex information
 232   vz          = 0;
 233   Bool_t vzOk = ReadVertex(event, vz);
 234
 235   fHEventsTr->Fill(vz);
 236   if (!vzOk) {
 237     if (fDebug > 3) {
 238       AliWarning("Failed to read vertex from ESD"); }
 239     return kNoVertex;
 240   }
 241   fHEventsTrVtx->Fill(vz);
 242
 243   // Get the vertex bin
 244   ivz = fHEventsTr->GetXaxis()->FindBin(vz);
 245   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) {
 246     if (fDebug > 3) {
 247       AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
 248                       vz, fHEventsTr->GetXaxis()->GetXmin(),
 249                       fHEventsTr->GetXaxis()->GetXmax())); }
 250     ivz = 0;
 251     return kBadVertex;
 252   }
 253   return kOk;
 254 }
 255
 256 //____________________________________________________________________
 257 Bool_t
 258 AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
 259 {
 260   triggers = 0;
 261
 262   // Get the analysis manager - should always be there
 263   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 264   if (!am) {
 265     AliWarning("No analysis manager defined!");
 266     return kFALSE;
 267   }
 268
 269   // Get the input handler - should always be there
 270   AliInputEventHandler* ih =
 271     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
 272   if (!ih) {
 273     AliWarning("No input handler");
 274     return kFALSE;
 275   }
 276
 277   // Get the physics selection - add that by using the macro
 278   // AddTaskPhysicsSelection.C
 279   AliPhysicsSelection* ps =
 280     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
 281   if (!ps) {
 282     AliWarning("No physics selection");
 283     return kFALSE;
 284   }
 285
 286   // Check if this is a collision candidate (INEL)
 287   Bool_t inel = ps->IsCollisionCandidate(esd);
 288   if (inel) {
 289     triggers |= AliAODForwardMult::kInel;
 290     fHTriggers->Fill(kInel+0.5);
 291   }
 292
 293   // IF this is inel, see if we have a tracklet
 294   if (inel) {
 295     const AliMultiplicity* spdmult = esd->GetMultiplicity();
 296     if (!spdmult) {
 297       AliWarning("No SPD multiplicity");
 298     }
 299     else {
 300       Int_t n = spdmult->GetNumberOfTracklets();
 301       for (Int_t j = 0; j < n; j++) {
 302         if(TMath::Abs(spdmult->GetEta(j)) < 1) {
 303           triggers |= AliAODForwardMult::kInelGt0;
 304           fHTriggers->Fill(kInelGt0+.5);
 305           break;
 306         }
 307       }
 308     }
 309   }
 310
 311   // Analyse some trigger stuff
 312   AliTriggerAnalysis ta;
 313   if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
 314     triggers |= AliAODForwardMult::kNSD;
 315     fHTriggers->Fill(kNSD+.5);
 316   }
 317
 318   // Get trigger stuff
 319   TString trigStr = esd->GetFiredTriggerClasses();
 320   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
 321     triggers |= AliAODForwardMult::kEmpty;
 322     fHTriggers->Fill(kEmpty+.5);
 323   }
 324
 325   if (trigStr.Contains("CINT1A-ABCE-NOPF-ALL")) {
 326     triggers |= AliAODForwardMult::kA;
 327     fHTriggers->Fill(kA+.5);
 328   }
 329
 330   if (trigStr.Contains("CINT1B-ABCE-NOPF-ALL")) {
 331     triggers |= AliAODForwardMult::kB;
 332     fHTriggers->Fill(kB+.5);
 333   }
 334
 335
 336   if (trigStr.Contains("CINT1C-ABCE-NOPF-ALL")) {
 337     triggers |= AliAODForwardMult::kC;
 338     fHTriggers->Fill(kC+.5);
 339   }
 340
 341   if (trigStr.Contains("CINT1-E-NOPF-ALL")) {
 342     triggers |= AliAODForwardMult::kE;
 343     fHTriggers->Fill(kE+.5);
 344   }
 345
 346   return kTRUE;
 347 }
 348 //____________________________________________________________________
 349 Bool_t
 350 AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
 351 {
 352   vz = 0;
 353   // Get the vertex
 354   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
 355   if (!vertex) {
 356     if (fDebug > 2) {
 357       AliWarning("No SPD vertex found in ESD"); }
 358     return kFALSE;
 359   }
 360
 361   // Check that enough tracklets contributed
 362   if(vertex->GetNContributors() <= 0) {
 363     if (fDebug > 2) {
 364       AliWarning(Form("Number of contributors to vertex is %d<=0",
 365                       vertex->GetNContributors())); }
 366     vz = 0;
 367     return kFALSE;
 368   }
 369
 370   // Check that the uncertainty isn't too large
 371   if (vertex->GetZRes() > fMaxVzErr) {
 372     if (fDebug > 2) {
 373       AliWarning(Form("Uncertaintity in Z of vertex is too large %f > %f",
 374                       vertex->GetZRes(), fMaxVzErr)); }
 375     return kFALSE;
 376   }
 377
 378   // Get the z coordiante
 379   vz = vertex->GetZ();
 380   return kTRUE;
 381 }
 382
 383 //____________________________________________________________________
 384 Bool_t
 385 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
 386 {
 387   fCollisionSystem =
 388     AliForwardUtil::ParseCollisionSystem(esd->GetBeamType());
 389   fEnergy          =
 390     AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,
 391                                               2 * esd->GetBeamEnergy());
 392   fField           =
 393     AliForwardUtil::ParseMagneticField(esd->GetMagneticField());
 394
 395   if (fCollisionSystem   == AliForwardUtil::kUnknown ||
 396       fEnergy            <= 0                        ||
 397       TMath::Abs(fField) >  10)
 398     return kFALSE;
 399
 400   return kTRUE;
 401 }
 402
 403 //____________________________________________________________________
 404 void
 405 AliFMDEventInspector::Print(Option_t*) const
 406 {
 407   char ind[gROOT->GetDirLevel()+1];
 408   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 409   ind[gROOT->GetDirLevel()] = '\0';
 410   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
 411   sNN.Strip(TString::kBoth, '0');
 412   sNN.ReplaceAll("GeV", " GeV");
 413   TString field(AliForwardUtil::MagneticFieldString(fField));
 414   field.ReplaceAll("p",  "+");
 415   field.ReplaceAll("m",  "-");
 416   field.ReplaceAll("kG", " kG");
 417
 418   std::cout << ind << "AliFMDEventInspector: " << GetName() << '\n'
 419             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
 420             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
 421             << ind << " System:                 "
 422             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
 423             << ind << " CMS energy per nucleon: " << sNN << '\n'
 424             << ind << " Field:                  " <<  field << std::endl;
 425 }
 426
 427
 428 //
 429 // EOF
 430 //
 431