PWGLF/FORWARD/analysis2/AliForwardQATask.cxx

   1 //
   2 // Calculate the multiplicity in the forward regions event-by-event
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODForwardMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliForwardQATask.h"
  15 #include "AliForwardUtil.h"
  16 #include "AliTriggerAnalysis.h"
  17 #include "AliPhysicsSelection.h"
  18 #include "AliLog.h"
  19 #include "AliESDEvent.h"
  20 #include "AliAODHandler.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliInputEventHandler.h"
  23 #include "AliForwardCorrectionManager.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliAODForwardMult.h"
  26 #include <TH1.h>
  27 #include <TDirectory.h>
  28 #include <TTree.h>
  29 #include <TROOT.h>
  30 #include <TStopwatch.h>
  31
  32 //====================================================================
  33 AliForwardQATask::AliForwardQATask()
  34   : AliAnalysisTaskSE(),
  35     fEnableLowFlux(false),
  36     fFirstEvent(true),
  37     fCorrManager(0),
  38     fESDFMD(),
  39     fHistos(),
  40     fEventInspector(),
  41     fEnergyFitter(),
  42     fSharingFilter(),
  43     fDensityCalculator(),
  44     fList(0),
  45     fDebug(0)
  46 {
  47   //
  48   // Constructor
  49   //
  50 }
  51
  52 //____________________________________________________________________
  53 AliForwardQATask::AliForwardQATask(const char* name)
  54   : AliAnalysisTaskSE(name),
  55     fEnableLowFlux(false),
  56     fFirstEvent(true),
  57     fCorrManager(0),
  58     fESDFMD(),
  59     fHistos(),
  60     fEventInspector("event"),
  61     fEnergyFitter("energy"),
  62     fSharingFilter("sharing"),
  63     fDensityCalculator("density"),
  64     fList(0),
  65     fDebug(0)
  66 {
  67   //
  68   // Constructor
  69   //
  70   // Parameters:
  71   //    name Name of task
  72   //
  73   DefineOutput(1, TList::Class());
  74   DefineOutput(2, TList::Class());
  75   fCorrManager = &AliForwardCorrectionManager::Instance();
  76   fEnergyFitter.SetNParticles(1); // Just find the 1st peak
  77   fEnergyFitter.SetDoMakeObject(false);
  78   fEnergyFitter.SetUseIncreasingBins(true);
  79   fEnergyFitter.SetDoFits(kTRUE);
  80   fEnergyFitter.SetLowCut(0.4);
  81   fEnergyFitter.SetFitRangeBinWidth(4);
  82   fEnergyFitter.SetMinEntries(1000);
  83 }
  84
  85 //____________________________________________________________________
  86 AliForwardQATask::AliForwardQATask(const AliForwardQATask& o)
  87   : AliAnalysisTaskSE(o),
  88     fEnableLowFlux(o.fEnableLowFlux),
  89     fFirstEvent(o.fFirstEvent),
  90     fCorrManager(o.fCorrManager),
  91     fESDFMD(o.fESDFMD),
  92     fHistos(o.fHistos),
  93     fEventInspector(o.fEventInspector),
  94     fEnergyFitter(o.fEnergyFitter),
  95     fSharingFilter(o.fSharingFilter),
  96     fDensityCalculator(o.fDensityCalculator),
  97     fList(o.fList),
  98     fDebug(o.fDebug)
  99 {
 100   //
 101   // Copy constructor
 102   //
 103   // Parameters:
 104   //    o Object to copy from
 105   //
 106   DefineOutput(1, TList::Class());
 107   DefineOutput(2, TList::Class());
 108 }
 109
 110 //____________________________________________________________________
 111 AliForwardQATask&
 112 AliForwardQATask::operator=(const AliForwardQATask& o)
 113 {
 114   //
 115   // Assignment operator
 116   //
 117   // Parameters:
 118   //    o Object to assign from
 119   //
 120   // Return:
 121   //    Reference to this object
 122   //
 123   if (&o == this) return *this;
 124   AliAnalysisTaskSE::operator=(o);
 125
 126   fEnableLowFlux     = o.fEnableLowFlux;
 127   fFirstEvent        = o.fFirstEvent;
 128   fCorrManager       = o.fCorrManager;
 129   fEventInspector    = o.fEventInspector;
 130   fEnergyFitter      = o.fEnergyFitter;
 131   fSharingFilter     = o.fSharingFilter;
 132   fDensityCalculator = o.fDensityCalculator;
 133   fHistos            = o.fHistos;
 134   fList              = o.fList;
 135   fDebug             = o.fDebug;
 136
 137   return *this;
 138 }
 139
 140 //____________________________________________________________________
 141 void
 142 AliForwardQATask::SetDebug(Int_t dbg)
 143 {
 144   //
 145   // Set debug level
 146   //
 147   // Parameters:
 148   //    dbg Debug level
 149   //
 150   fDebug = dbg;
 151   fEventInspector.SetDebug(dbg);
 152   fEnergyFitter.SetDebug(dbg);
 153   fSharingFilter.SetDebug(dbg);
 154   fDensityCalculator.SetDebug(dbg);
 155 }
 156
 157 //____________________________________________________________________
 158 Bool_t
 159 AliForwardQATask::CheckCorrections(UInt_t what)
 160 {
 161   //
 162   // Check if all needed corrections are there and accounted for.  If not,
 163   // do a Fatal exit
 164   //
 165   // Parameters:
 166   //    what Which corrections is needed
 167   //
 168   // Return:
 169   //    true if all present, false otherwise
 170   //
 171
 172   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 173   // Check that we have the energy loss fits, needed by
 174   //   AliFMDSharingFilter
 175   //   AliFMDDensityCalculator
 176   if (what & AliForwardCorrectionManager::kELossFits) {
 177     if (!fcm.GetELossFit()) {
 178       AliWarning("No energy loss fits");
 179
 180       // Fall-back values if we do not have the energy loss fits
 181       AliFMDMultCuts& sfLCuts = GetSharingFilter().GetLCuts();
 182       if (sfLCuts.GetMethod() != AliFMDMultCuts::kFixed) {
 183         Double_t cut = 0.3;
 184         AliWarningF("Using fixed cut @ %f for the lower bound "
 185                     "of the sharing filter", cut);
 186         sfLCuts.SetMultCuts(cut);
 187       }
 188       AliFMDMultCuts& sfHCuts = GetSharingFilter().GetHCuts();
 189       if (sfHCuts.GetMethod() != AliFMDMultCuts::kFixed) {
 190         Double_t cut = 100;
 191         AliWarningF("Using fixed cut @ %f for the upper bound "
 192                     "of the sharing filter", cut);
 193         sfHCuts.SetMultCuts(cut);
 194       }
 195       AliFMDMultCuts& dcCuts  = GetDensityCalculator().GetCuts();
 196       if (dcCuts.GetMethod() != AliFMDMultCuts::kFixed) {
 197         Double_t cut = 0.3;
 198         AliWarningF("Using fixed cut @ %f for the lower bound "
 199                     "of the density calculator", cut);
 200         dcCuts.SetMultCuts(cut);
 201       }
 202     }
 203     else
 204       fcm.GetELossFit()->CacheBins(GetDensityCalculator().GetMinQuality());
 205   }
 206
 207   return true;
 208 }
 209
 210 //____________________________________________________________________
 211 Bool_t
 212 AliForwardQATask::ReadCorrections(const TAxis*& pe,
 213                                   const TAxis*& pv,
 214                                   Bool_t        mc,
 215                                   Bool_t        sat)
 216 {
 217   //
 218   // Read corrections
 219   //
 220   //
 221   UInt_t what = AliForwardCorrectionManager::kAll;
 222   what ^= AliForwardCorrectionManager::kDoubleHit;
 223   what ^= AliForwardCorrectionManager::kVertexBias;
 224   what ^= AliForwardCorrectionManager::kAcceptance;
 225   what ^= AliForwardCorrectionManager::kMergingEfficiency;
 226
 227   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 228   if (!fcm.Init(GetEventInspector().GetRunNumber(),
 229                 GetEventInspector().GetCollisionSystem(),
 230                 GetEventInspector().GetEnergy(),
 231                 GetEventInspector().GetField(),
 232                 mc,
 233                 sat,
 234                 what,
 235                 false)) return false;
 236   if (!CheckCorrections(what)) {
 237     return false;
 238   }
 239
 240   // Sett our persistency pointer
 241   // fCorrManager = &fcm;
 242
 243   // Get the eta axis from the secondary maps - if read in
 244   if (!pe) {
 245     pe = fcm.GetEtaAxis();
 246     if (!pe) AliFatal("No eta axis defined");
 247   }
 248   // Get the vertex axis from the secondary maps - if read in
 249   if (!pv) {
 250     pv = fcm.GetVertexAxis();
 251     if (!pv) AliFatal("No vertex axis defined");
 252   }
 253
 254   return true;
 255 }
 256
 257 //____________________________________________________________________
 258 AliESDEvent*
 259 AliForwardQATask::GetESDEvent()
 260 {
 261   //
 262   // Get the ESD event. IF this is the first event, initialise
 263   //
 264   DGUARD(fDebug,2,"Get the ESD event");
 265   if (IsZombie()) {
 266     DMSG(fDebug,3,"We're a Zombie - bailing out");
 267     return 0;
 268   }
 269   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 270   if (!esd) {
 271     AliWarning("No ESD event found for input event");
 272     return 0;
 273   }
 274
 275   // On the first event, initialize the parameters
 276   if (fFirstEvent && esd->GetESDRun()) {
 277     GetEventInspector().ReadRunDetails(esd);
 278
 279     AliInfoF("Initializing with parameters from the ESD:\n"
 280              "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 281              "         AliESDEvent::GetBeamType()     ->%s\n"
 282              "         AliESDEvent::GetCurrentL3()    ->%f\n"
 283              "         AliESDEvent::GetMagneticField()->%f\n"
 284              "         AliESDEvent::GetRunNumber()    ->%d",
 285              esd->GetBeamEnergy(),
 286              esd->GetBeamType(),
 287              esd->GetCurrentL3(),
 288              esd->GetMagneticField(),
 289              esd->GetRunNumber());
 290
 291
 292     if (!SetupForData()) {
 293       AliWarning("Initialisation of sub algorithms failed!");
 294       SetZombie(true);
 295       esd = 0;
 296       return 0;
 297     }
 298     AliInfoF("Clearing first event flag from %s to false",
 299              fFirstEvent ? "true" : "false");
 300     fFirstEvent = false;
 301   }
 302   return esd;
 303 }
 304 //____________________________________________________________________
 305 Bool_t
 306 AliForwardQATask::SetupForData()
 307 {
 308   //
 309   // Initialise the sub objects and stuff.  Called on first event
 310   //
 311   //
 312   const TAxis* pe = 0;
 313   const TAxis* pv = 0;
 314
 315   if (!ReadCorrections(pe,pv))  {
 316     AliWarning("Using default binning");
 317     pv = new TAxis(10,-10, 10);
 318     pe = new TAxis(240,-6,6);
 319   }
 320
 321   fHistos.Init(*pe);
 322
 323   fEventInspector.SetupForData(*pv);
 324   fEnergyFitter.SetupForData(*pe);
 325   fSharingFilter.SetupForData(*pe);
 326   fDensityCalculator.SetupForData(*pe);
 327
 328   this->Print();
 329
 330   return true;
 331 }
 332
 333 //____________________________________________________________________
 334 void
 335 AliForwardQATask::UserCreateOutputObjects()
 336 {
 337   //
 338   // Create output objects
 339   //
 340   //
 341   fList = new TList;
 342   fList->SetOwner();
 343
 344   fEventInspector.CreateOutputObjects(fList);
 345   fEnergyFitter.CreateOutputObjects(fList);
 346   fSharingFilter.CreateOutputObjects(fList);
 347   fDensityCalculator.CreateOutputObjects(fList);
 348
 349   PostData(1, fList);
 350 }
 351 //____________________________________________________________________
 352 void
 353 AliForwardQATask::UserExec(Option_t*)
 354 {
 355   //
 356   // Process each event
 357   //
 358   // Parameters:
 359   //    option Not used
 360   //
 361   DGUARD(fDebug,1,"Process the input event");
 362
 363   // static Int_t cnt = 0;
 364   // cnt++;
 365   // Get the input data
 366   AliESDEvent* esd = GetESDEvent();
 367   if (!esd) {
 368     AliWarning("Got no ESD event");
 369     return;
 370   }
 371   if (fFirstEvent) {
 372     // If the first event flag wasn't cleared in the above call to
 373     // GetESDEvent, we should not do anything, since nothing has been
 374     // initialised yet, so we opt out here (with a warning)
 375     AliWarning("Nothing has been initialized yet, opt'ing out");
 376     return;
 377   }
 378
 379   // Clear stuff
 380   fHistos.Clear();
 381   fESDFMD.Clear();
 382
 383   Bool_t   lowFlux   = kFALSE;
 384   UInt_t   triggers  = 0;
 385   UShort_t ivz       = 0;
 386   TVector3 ip;
 387   Double_t cent      = -1;
 388   UShort_t nClusters = 0;
 389   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 390                                                ivz, ip, cent, nClusters);
 391
 392   Bool_t ok = true;
 393   if (found & AliFMDEventInspector::kNoEvent)    ok = false;
 394   if (found & AliFMDEventInspector::kNoTriggers) ok = false;
 395   if (found & AliFMDEventInspector::kNoSPD)      ok = false;
 396   if (found & AliFMDEventInspector::kNoFMD)      ok = false;
 397   if (found & AliFMDEventInspector::kNoVertex)   ok = false;
 398   if (triggers & AliAODForwardMult::kPileUp)     ok = false;
 399   if (found & AliFMDEventInspector::kBadVertex)  ok = false;
 400   if (!ok) {
 401     DMSG(fDebug,2,"Event failed selection: %s",
 402          AliFMDEventInspector::CodeString(found));
 403     return;
 404   }
 405   DMSG(fDebug,2,"Event triggers: %s",
 406        AliAODForwardMult::GetTriggerString(triggers));
 407
 408   // We we do not want to use low flux specific code, we disable it here.
 409   if (!fEnableLowFlux) lowFlux = false;
 410
 411   // Get FMD data
 412   AliESDFMD* esdFMD = esd->GetFMDData();
 413
 414   // Run the energy loss fitter
 415   if (!fEnergyFitter.Accumulate(*esdFMD, cent,
 416                                 triggers & AliAODForwardMult::kEmpty)) {
 417     AliWarning("Energy fitter failed");
 418     return;
 419   }
 420
 421   //  // Apply the sharing filter (or hit merging or clustering if you like)
 422   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) {
 423     AliWarning("Sharing filter failed!");
 424     return;
 425   }
 426
 427   // Calculate the inclusive charged particle density
 428   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
 429     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
 430     AliWarning("Density calculator failed!");
 431     return;
 432   }
 433
 434   PostData(1, fList);
 435 }
 436
 437 //____________________________________________________________________
 438 void
 439 AliForwardQATask::Terminate(Option_t*)
 440 {
 441   //
 442   // End of job
 443   //
 444   // Parameters:
 445   //    option Not used
 446   //
 447   if (fDebug) AliInfo("In Forwards terminate");
 448   TStopwatch swt;
 449   swt.Start();
 450
 451   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 452   if (!list) {
 453     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 454     if (GetOutputData(1)) GetOutputData(1)->Print();
 455     return;
 456   }
 457
 458   // Make a deep copy and post that as output 2
 459   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
 460                                                       list->GetName())));
 461   list2->SetOwner();
 462
 463   // Get our histograms from the container
 464   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 465   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 466   TH1I* hEventsAcc   = 0;
 467   TH1I* hTriggers    = 0;
 468   if (!fEventInspector.FetchHistograms(list,
 469                                        hEventsTr,
 470                                        hEventsTrVtx,
 471                                        hEventsAcc,
 472                                        hTriggers)) {
 473     AliError(Form("Didn't get histograms from event selector "
 474                   "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p)",
 475                   hEventsTr, hEventsTrVtx,hEventsAcc));
 476     return;
 477   }
 478
 479   TStopwatch swf;
 480   swf.Start();
 481   fEnergyFitter.Fit(list2);
 482   swf.Stop();
 483   AliInfoF("Fitting took %d real-time seconds, and %f CPU seconds",
 484            Int_t(swf.RealTime()), swf.CpuTime());
 485
 486   fSharingFilter.Terminate(list,list2,Int_t(hEventsTr->Integral()));
 487   fDensityCalculator.Terminate(list,list2,Int_t(hEventsTrVtx->Integral()));
 488
 489   if (fDebug) AliInfoF("Posting post processing results to %s",
 490                        list2->GetName());
 491   PostData(2, list2);
 492
 493   swt.Stop();
 494   AliInfoF("Terminate took %d real-time seconds, and %f CPU seconds",
 495            Int_t(swt.RealTime()), swt.CpuTime());
 496
 497
 498 }
 499
 500 //____________________________________________________________________
 501 void
 502 AliForwardQATask::Print(Option_t* option) const
 503 {
 504   //
 505   // Print information
 506   //
 507   // Parameters:
 508   //    option Not used
 509   //
 510
 511   std::cout << ClassName() << ": " << GetName() << "\n"
 512             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 513             << "\n"
 514             << "  Off-line trigger mask:  0x"
 515             << std::hex     << std::setfill('0')
 516             << std::setw (8) << fOfflineTriggerMask
 517             << std::dec     << std::setfill (' ') << std::endl;
 518   gROOT->IncreaseDirLevel();
 519   if (fCorrManager) fCorrManager->Print();
 520   else
 521     std::cout << "  Correction manager not set yet" << std::endl;
 522   GetEventInspector()   .Print(option);
 523   GetEnergyFitter()     .Print(option);
 524   GetSharingFilter()    .Print(option);
 525   GetDensityCalculator().Print(option);
 526   gROOT->DecreaseDirLevel();
 527 }
 528
 529 //
 530 // EOF
 531 //