]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx
Various updates, and prep for Train base class
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMultiplicityTask.cxx
index 20b3052299ad24e21d3a8f49ca2f86fe05b9d864..7afb3ef3353256b6c2b6ad82ec3b994876a3028f 100644 (file)
 #include "AliForwardCorrectionManager.h"
 #include "AliAnalysisManager.h"
 #include <TH1.h>
+#include <TH3D.h>
 #include <TDirectory.h>
 #include <TTree.h>
 #include <TROOT.h>
-
+#include <TStopwatch.h>
+#include <TProfile.h>
+// #define ENABLE_TIMING
+#ifndef ENABLE_TIMING
+# define MAKE_SW(NAME) do {} while(false)
+# define START_SW(NAME) do {} while(false)
+# define FILL_SW(NAME,WHICH) do {} while(false)
+#else
+# define MAKE_SW(NAME) TStopwatch NAME
+# define START_SW(NAME) if (fDoTiming) NAME.Start(true)
+# define FILL_SW(NAME,WHICH)                           \
+  if (fDoTiming) fHTiming->Fill(WHICH,NAME.CpuTime())
+#endif
 
 //====================================================================
 AliForwardMultiplicityTask::AliForwardMultiplicityTask()
   : AliForwardMultiplicityBase(),
-    fHData(0),
     fESDFMD(),
-    fHistos(),
-    fAODFMD(),
-    fAODEP(),
-    fRingSums(),
     fEventInspector(),
+    fESDFixer(),
     fSharingFilter(),
     fDensityCalculator(),
     fCorrections(),
     fHistCollector(),
-    fEventPlaneFinder(),
-    fList(0)
+    fEventPlaneFinder()
 {
   // 
   // Constructor
   //
+  DGUARD(fDebug, 3,"Default CTOR of AliForwardMultiplicityTask");
 }
 
 //____________________________________________________________________
 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
   : AliForwardMultiplicityBase(name),
-    fHData(0),
     fESDFMD(),
-    fHistos(),
-    fAODFMD(false),
-    fAODEP(false),
-    fRingSums(),
     fEventInspector("event"),
+    fESDFixer("esdFizer"),
     fSharingFilter("sharing"), 
     fDensityCalculator("density"),
     fCorrections("corrections"),
     fHistCollector("collector"),
-    fEventPlaneFinder("eventplane"),
-    fList(0)
+    fEventPlaneFinder("eventplane")
 {
   // 
   // Constructor 
@@ -72,168 +76,55 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
   // Parameters:
   //    name Name of task 
   //
-  DefineOutput(1, TList::Class());
-}
-
-//____________________________________________________________________
-AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
-  : AliForwardMultiplicityBase(o),
-    fHData(o.fHData),
-    fESDFMD(o.fESDFMD),
-    fHistos(o.fHistos),
-    fAODFMD(o.fAODFMD),
-    fAODEP(o.fAODEP),
-    fRingSums(o.fRingSums),
-    fEventInspector(o.fEventInspector),
-    fSharingFilter(o.fSharingFilter),
-    fDensityCalculator(o.fDensityCalculator),
-    fCorrections(o.fCorrections),
-    fHistCollector(o.fHistCollector),
-    fEventPlaneFinder(o.fEventPlaneFinder),
-    fList(o.fList) 
-{
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //    o Object to copy from 
-  //
-  DefineOutput(1, TList::Class());
+  DGUARD(fDebug, 3,"named CTOR of AliForwardMultiplicityTask: %s", name);
 }
 
-//____________________________________________________________________
-AliForwardMultiplicityTask&
-AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
-{
-  // 
-  // Assignment operator 
-  // 
-  // Parameters:
-  //    o Object to assign from 
-  // 
-  // Return:
-  //    Reference to this object 
-  //
-  if (&o == this) return *this;
-  AliForwardMultiplicityBase::operator=(o);
-
-  fHData             = o.fHData;
-  fEventInspector    = o.fEventInspector;
-  fSharingFilter     = o.fSharingFilter;
-  fDensityCalculator = o.fDensityCalculator;
-  fCorrections       = o.fCorrections;
-  fHistCollector     = o.fHistCollector;
-  fEventPlaneFinder  = o.fEventPlaneFinder;
-  fHistos            = o.fHistos;
-  fAODFMD            = o.fAODFMD;
-  fAODEP             = o.fAODEP;
-  fRingSums          = o.fRingSums;
-  fList              = o.fList;
-
-  return *this;
-}
 
 //____________________________________________________________________
 void
-AliForwardMultiplicityTask::SetDebug(Int_t dbg)
+AliForwardMultiplicityTask::SetDoTiming(Bool_t enable)
 {
-  // 
-  // Set debug level 
-  // 
-  // Parameters:
-  //    dbg Debug level
-  //
-  fEventInspector.SetDebug(dbg);
-  fSharingFilter.SetDebug(dbg);
-  fDensityCalculator.SetDebug(dbg);
-  fCorrections.SetDebug(dbg);
-  fHistCollector.SetDebug(dbg);
-  fEventPlaneFinder.SetDebug(dbg);
+#ifndef ENABLE_TIMING
+  if (enable) 
+    AliWarning("Timing of task explicitly disabled in compilation");
+#else 
+  fDoTiming = enable;
+#endif
 }
-
+      
 //____________________________________________________________________
 void
-AliForwardMultiplicityTask::InitializeSubs()
+AliForwardMultiplicityTask::PreCorrections(const AliESDEvent* esd)
 {
-  // 
-  // Initialise the sub objects and stuff.  Called on first event 
-  // 
-  //
-  const TAxis* pe = 0;
-  const TAxis* pv = 0;
-
-  if (!ReadCorrections(pe,pv)) return;
-
-  fHistos.Init(*pe);
-  fAODFMD.Init(*pe);
-  fAODEP.Init(*pe);
-  fRingSums.Init(*pe);
-
-  fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
-  fHData->SetStats(0);
-  fHData->SetDirectory(0);
-  fList->Add(fHData);
-
-  TList* rings = new TList;
-  rings->SetName("ringSums");
-  rings->SetOwner();
-  fList->Add(rings);
-
-  rings->Add(fRingSums.Get(1, 'I'));
-  rings->Add(fRingSums.Get(2, 'I'));
-  rings->Add(fRingSums.Get(2, 'O'));
-  rings->Add(fRingSums.Get(3, 'I'));
-  rings->Add(fRingSums.Get(3, 'O'));
-  fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
-  fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
-  fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
-  fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
-  fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
-
-  fEventInspector.Init(*pv);
-  fSharingFilter.Init();
-  fDensityCalculator.Init(*pe);
-  fCorrections.Init(*pe);
-  fHistCollector.Init(*pv,*pe);
-  fEventPlaneFinder.Init(*pe);
-
-  this->Print();
+  if (!esd) return; 
+  
+  AliESDFMD* esdFMD = esd->GetFMDData();  
+  if (!esdFMD) return;
+
+  Int_t tgt = GetESDFixer().FindTargetNoiseFactor(*esdFMD, false);
+  if (tgt <= 0) {
+    // If the target noise factor is 0 or less, disable the noise/gain
+    // correction.
+    GetESDFixer().SetRecoNoiseFactor(4);
+    fNeededCorrections ^= AliForwardCorrectionManager::kNoiseGain;
+  }
+  else 
+    AliWarning("The noise corrector has been enabled!");
 }
-
 //____________________________________________________________________
-void
-AliForwardMultiplicityTask::UserCreateOutputObjects()
+Bool_t
+AliForwardMultiplicityTask::PreEvent()
 {
-  // 
-  // Create output objects 
-  // 
-  //
-  fList = new TList;
-  fList->SetOwner();
-
-  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
-  AliAODHandler*      ah = 
-    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
-  if (!ah) AliFatal("No AOD output handler set in analysis manager");
-    
-    
-  TObject* obj = &fAODFMD;
-  ah->AddBranch("AliAODForwardMult", &obj);
-  TObject* epobj = &fAODEP;
-  ah->AddBranch("AliAODForwardEP", &epobj);
-
-  fEventInspector.DefineOutput(fList);
-  fSharingFilter.DefineOutput(fList);
-  fDensityCalculator.DefineOutput(fList);
-  fCorrections.DefineOutput(fList);
-  fHistCollector.DefineOutput(fList);
-  fEventPlaneFinder.DefineOutput(fList);
-
-  PostData(1, fList);
+  // Clear stuff 
+  fHistos.Clear();
+  fESDFMD.Clear();
+  fAODFMD.Clear();
+  fAODEP.Clear();
+  return true;
 }
 //____________________________________________________________________
-void
-AliForwardMultiplicityTask::UserExec(Option_t*)
+Bool_t
+AliForwardMultiplicityTask::Event(AliESDEvent& esd)
 {
   // 
   // Process each event 
@@ -241,29 +132,32 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   // Parameters:
   //    option Not used
   //  
-
-  // static Int_t cnt = 0;
-  // cnt++;
-  // Get the input data 
-  AliESDEvent* esd = GetESDEvent();
-
-  // Clear stuff 
-  fHistos.Clear();
-  fESDFMD.Clear();
-  fAODFMD.Clear();
-  fAODEP.Clear();
+  MAKE_SW(total);
+  MAKE_SW(individual);
+  START_SW(total);
   
+  DGUARD(fDebug,1,"Process the input event");
+
+  // Inspect the event
+  START_SW(individual);
   Bool_t   lowFlux   = kFALSE;
   UInt_t   triggers  = 0;
   UShort_t ivz       = 0;
-  Double_t vz        = 0;
+  TVector3 ip;
   Double_t cent      = -1;
   UShort_t nClusters = 0;
-  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
-                                              ivz, vz, cent, nClusters);
-  
-  if (found & AliFMDEventInspector::kNoEvent)    return;
-  if (found & AliFMDEventInspector::kNoTriggers) return;
+  UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux, 
+                                              ivz, ip, cent, nClusters);
+  FILL_SW(individual,kTimingEventInspector);
+
+  if (found & AliFMDEventInspector::kNoEvent)    { 
+    fHStatus->Fill(1);
+    return false;
+  }
+  if (found & AliFMDEventInspector::kNoTriggers) {
+    fHStatus->Fill(2);
+    return false;
+  } 
 
   // Set trigger bits, and mark this event for storage 
   fAODFMD.SetTriggerBits(triggers);
@@ -273,115 +167,124 @@ AliForwardMultiplicityTask::UserExec(Option_t*)
   fAODFMD.SetNClusters(nClusters);
   MarkEventForStore();
  
-  if (found & AliFMDEventInspector::kNoSPD)      return;
-  if (found & AliFMDEventInspector::kNoFMD)      return;
-  if (found & AliFMDEventInspector::kNoVertex)   return;
-  
-  if (triggers & AliAODForwardMult::kPileUp) return;
-  
-  fAODFMD.SetIpZ(vz);
-
-  if (found & AliFMDEventInspector::kBadVertex) return;
+  // Do not check if SPD data is there - potential bias 
+  // if (found & AliFMDEventInspector::kNoSPD) {
+  //   fHStatus->Fill(3);
+  //   return false;
+  // }
+  if (found    & AliFMDEventInspector::kNoFMD) {
+    fHStatus->Fill(4);
+    return false;
+  }
+  if (found    & AliFMDEventInspector::kNoVertex) {
+    fHStatus->Fill(5);
+    return false;
+  }
+  // Also analyse pile-up events - we'll remove them in later steps. 
+  if (triggers & AliAODForwardMult::kPileUp) {
+    fHStatus->Fill(6);
+    return false;
+  }
+  fAODFMD.SetIpZ(ip.Z());
+  if (found & AliFMDEventInspector::kBadVertex) {
+    fHStatus->Fill(7);
+    return false;
+  }
 
   // We we do not want to use low flux specific code, we disable it here. 
   if (!fEnableLowFlux) lowFlux = false;
 
   // Get FMD data 
-  AliESDFMD* esdFMD = esd->GetFMDData();
-  //  // Apply the sharing filter (or hit merging or clustering if you like)
-  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
+  AliESDFMD* esdFMD = esd.GetFMDData();  
+
+  // Fix up the the ESD 
+  GetESDFixer().Fix(*esdFMD, ip.Z());
+
+  // Apply the sharing filter (or hit merging or clustering if you like)
+  START_SW(individual);
+  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) { 
     AliWarning("Sharing filter failed!");
-    return;
+    fHStatus->Fill(8);
+    return false;
   }
+  FILL_SW(individual,kTimingSharingFilter);
   
   // Calculate the inclusive charged particle density 
-  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) { 
+  START_SW(individual);
+  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) { 
     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { 
     AliWarning("Density calculator failed!");
-    return;
+    fHStatus->Fill(9);
+    return false;
   }
+  FILL_SW(individual,kTimingDensityCalculator);
 
+  // Check if we should do the event plane finder
   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
-    if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos))
+    START_SW(individual);
+    if (!fEventPlaneFinder.FindEventplane(&esd, fAODEP, 
+                                         &(fAODFMD.GetHistogram()), &fHistos)){
       AliWarning("Eventplane finder failed!");
+      fHStatus->Fill(10);
+    }
+    FILL_SW(individual,kTimingEventPlaneFinder);
+  }
+  
+  // Check how many rings have been marked for skipping 
+  Int_t nSkip = 0;
+  for (UShort_t d=1; d<=3; d++) { 
+    for (UShort_t q=0; q<=(d/2); q++) { 
+      TH2D* h = fHistos.Get(d,q == 0 ? 'I' : 'O');
+      if (h && h->TestBit(AliForwardUtil::kSkipRing)) nSkip++;
+    }
+  }
+  if (nSkip > 0) {
+    // Skip the rest if we have too many outliers 
+    fHStatus->Fill(11);
+    return false;
   }
   
   // Do the secondary and other corrections. 
+  START_SW(individual);
   if (!fCorrections.Correct(fHistos, ivz)) { 
     AliWarning("Corrections failed");
-    return;
+    fHStatus->Fill(12);
+    return false;
   }
-
-  if (!fHistCollector.Collect(fHistos, fRingSums, 
-                             ivz, fAODFMD.GetHistogram())) {
+  FILL_SW(individual,kTimingCorrections);
+
+  // Check if we should add to internal caches 
+  Bool_t add = (fAODFMD.IsTriggerBits(fAddMask) && nSkip < 1);
+
+  // Collect our `super' histogram 
+  START_SW(individual);
+  if (!fHistCollector.Collect(fHistos, 
+                             fRingSums, 
+                             ivz, 
+                             fAODFMD.GetHistogram(),
+                             fAODFMD.GetCentrality(),
+                             false, 
+                             add)) {
     AliWarning("Histogram collector failed");
-    return;
+    fHStatus->Fill(13);
+    return false;
   }
+  FILL_SW(individual,kTimingHistCollector);
 
-  if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
-    fHData->Add(&(fAODFMD.GetHistogram()));
-
-  PostData(1, fList);
-}
-
-//____________________________________________________________________
-void
-AliForwardMultiplicityTask::Terminate(Option_t*)
-{
-  // 
-  // End of job
-  // 
-  // Parameters:
-  //    option Not used 
-  //
-
-  TList* list = dynamic_cast<TList*>(GetOutputData(1));
-  if (!list) {
-    AliError(Form("No output list defined (%p)", GetOutputData(1)));
-    if (GetOutputData(1)) GetOutputData(1)->Print();
-    return;
-  }
-  
-  // Get our histograms from the container 
-  TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
-  TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
-  TH1I* hTriggers    = 0;
-  if (!fEventInspector.FetchHistograms(list, hEventsTr, 
-                                      hEventsTrVtx, hTriggers)) { 
-    AliError(Form("Didn't get histograms from event selector "
-                 "(hEventsTr=%p,hEventsTrVtx=%p)", 
-                 hEventsTr, hEventsTrVtx));
-    list->ls();
-    return;
+  if (!add) {
+    fHStatus->Fill(14);
   }
-
-  TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
-  if (!hData) { 
-    AliError(Form("Couldn't get our summed histogram from output "
-                 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
-    list->ls();
-    return;
+  else {
+    // Collect rough Min. Bias result
+    fHData->Add(&(fAODFMD.GetHistogram()));
+    fHStatus->Fill(15);
   }
+  FILL_SW(total,kTimingTotal);
   
-  // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
-  TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
-  TH1D* norm   = hData->ProjectionX("norm",   0,  0,  "");
-  dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
-  dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
-  dNdeta->Divide(norm);
-  dNdeta->SetStats(0);
-  dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
-               "width");
-  list->Add(dNdeta);
-  list->Add(norm);
-
-  MakeRingdNdeta(list, "ringSums", list, "ringResults");
-
-  fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
-  fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
-  fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
+  return true;
 }
 
+
 //
 // EOF
 //