]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliForwardMCMultiplicityTask.cxx
index 778432f5daa3610418e0a1c81ae49f91488eece6..2d52956590f553b812bae399810d24826f923231 100644 (file)
 //====================================================================
 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
   : AliForwardMultiplicityBase(),
-    fHData(0),
     fESDFMD(),
-    fHistos(),
-    fAODFMD(),
-    fAODEP(),
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(),
-    fRingSums(),
     fMCRingSums(),
     fPrimary(0),
     fEventInspector(),
@@ -47,8 +42,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
     fDensityCalculator(),
     fCorrections(),
     fHistCollector(),
-    fEventPlaneFinder(),
-    fList(0)
+    fEventPlaneFinder()
 {
   // 
   // Constructor
@@ -58,15 +52,10 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
 //____________________________________________________________________
 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
   : AliForwardMultiplicityBase(name), 
-    fHData(0),
     fESDFMD(),
-    fHistos(),
-    fAODFMD(kFALSE),
-    fAODEP(kFALSE),
     fMCESDFMD(),
     fMCHistos(),
     fMCAODFMD(kTRUE),
-    fRingSums(),
     fMCRingSums(),
     fPrimary(0),
     fEventInspector("event"),
@@ -74,8 +63,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
     fDensityCalculator("density"),
     fCorrections("corrections"),
     fHistCollector("collector"),
-    fEventPlaneFinder("eventplane"),
-    fList(0)
+    fEventPlaneFinder("eventplane")
 {
   // 
   // Constructor 
@@ -83,143 +71,64 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
   // Parameters:
   //    name Name of task 
   //
-  DefineOutput(1, TList::Class());
 }
 
 //____________________________________________________________________
-AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
-  : AliForwardMultiplicityBase(o),
-    fHData(o.fHData),
-    fESDFMD(o.fESDFMD),
-    fHistos(o.fHistos),
-    fAODFMD(o.fAODFMD),
-    fAODEP(o.fAODEP),
-    fMCESDFMD(o.fMCESDFMD),
-    fMCHistos(o.fMCHistos),
-    fMCAODFMD(o.fMCAODFMD),
-    fRingSums(o.fRingSums),
-    fMCRingSums(o.fMCRingSums),
-    fPrimary(o.fPrimary),
-    fEventInspector(o.fEventInspector),
-    fSharingFilter(o.fSharingFilter),
-    fDensityCalculator(o.fDensityCalculator),
-    fCorrections(o.fCorrections),
-    fHistCollector(o.fHistCollector),
-    fEventPlaneFinder(o.fEventPlaneFinder),
-    fList(o.fList) 
+void
+AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
 {
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //    o Object to copy from 
-  //
-  DefineOutput(1, TList::Class());
+  fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
+  fCorrections.SetSecondaryForMC(!use);
 }
 
 //____________________________________________________________________
-AliForwardMCMultiplicityTask&
-AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
+void
+AliForwardMCMultiplicityTask::CreateBranches(AliAODHandler* ah)
 {
   // 
-  // Assignment operator 
-  // 
-  // Parameters:
-  //    o Object to assign from 
+  // Create output objects 
   // 
-  // Return:
-  //    Reference to this object 
   //
-  if (&o == this) return *this;
-  AliForwardMultiplicityBase::operator=(o);
+  AliForwardMultiplicityBase::CreateBranches(ah);
 
-  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;
-  fMCHistos          = o.fMCHistos;
-  fMCAODFMD          = o.fMCAODFMD;
-  fRingSums          = o.fRingSums;
-  fMCRingSums        = o.fMCRingSums;
-  fPrimary           = o.fPrimary;
-  fList              = o.fList;
+  TObject* mcobj = &fMCAODFMD;
+  ah->AddBranch("AliAODForwardMult", &mcobj);
 
-  return *this;
+  fPrimary = new TH2D("primary", "MC Primaries", 1,0,1,20,0,TMath::TwoPi());
+  fPrimary->SetXTitle("#eta");
+  fPrimary->SetYTitle("#varphi [radians]");
+  fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
+  fPrimary->Sumw2();
+  fPrimary->SetStats(0);
+  fPrimary->SetDirectory(0);
+    
+  ah->AddBranch("TH2D", &fPrimary);
 }
 
-//____________________________________________________________________
-void
-AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
-{
-  // 
-  // 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);
-}
-//____________________________________________________________________
-void
-AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
-{
-  fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
-  fCorrections.SetSecondaryForMC(!use);
-}
 
 //____________________________________________________________________
-Bool_t
-AliForwardMCMultiplicityTask::InitializeSubs()
+void
+AliForwardMCMultiplicityTask::InitMembers(const TAxis& eta, const TAxis& vertex)
 {
   // 
   // Initialise the sub objects and stuff.  Called on first event 
   // 
   //
-  const TAxis* pe = 0;
-  const TAxis* pv = 0;
-
-  if (!ReadCorrections(pe,pv,true)) return false;
+  AliForwardMultiplicityBase::InitMembers(eta, vertex);
 
-  fHistos.Init(*pe);
-  fAODFMD.Init(*pe);
-  fAODEP.Init(*pe);
-  fMCHistos.Init(*pe);
-  fMCAODFMD.Init(*pe);
-  fRingSums.Init(*pe);
-  fMCRingSums.Init(*pe);
+  fMCHistos.Init(eta);
+  fMCAODFMD.Init(eta);
+  fMCRingSums.Init(eta);
 
-  fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
-  fHData->SetStats(0);
-  fHData->SetDirectory(0);
+  AliForwardUtil::Histos::RebinEta(fPrimary, eta);
+  DMSG(fDebug,0,"Primary histogram rebinned to %d,%f,%f eta axis %d,%f,%f", 
+       fPrimary->GetXaxis()->GetNbins(), 
+       fPrimary->GetXaxis()->GetXmin(),
+       fPrimary->GetXaxis()->GetXmax(),
+       eta.GetNbins(), 
+       eta.GetXmin(),
+       eta.GetXmax());
 
-  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'));
 
   TList* mcRings = new TList;
   mcRings->SetName("mcRingSums");
@@ -236,69 +145,28 @@ AliForwardMCMultiplicityTask::InitializeSubs()
   fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
   fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
   fMCRingSums.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();
-
-  return true;
 }
 
 //____________________________________________________________________
-void
-AliForwardMCMultiplicityTask::UserCreateOutputObjects()
+Bool_t
+AliForwardMCMultiplicityTask::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* mcobj = &fMCAODFMD;
-  ah->AddBranch("AliAODForwardMult", &mcobj);
-
-  TObject* epobj = &fAODEP;
-  ah->AddBranch("AliAODForwardEP", &epobj);
-  
-  fPrimary = new TH2D("primary", "MC Primaries", 
-                     200, -4, 6, 20, 0, 2*TMath::Pi());
-  fPrimary->SetXTitle("#eta");
-  fPrimary->SetYTitle("#varphi [radians]");
-  fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
-  fPrimary->Sumw2();
-  fPrimary->SetStats(0);
-  fPrimary->SetDirectory(0);
-  ah->AddBranch("TH2D", &fPrimary);
-
-  fEventInspector.DefineOutput(fList);
-  fSharingFilter.DefineOutput(fList);
-  fDensityCalculator.DefineOutput(fList);
-  fCorrections.DefineOutput(fList);
-  fHistCollector.DefineOutput(fList);
-  fEventPlaneFinder.DefineOutput(fList);
-
-  PostData(1, fList);
+ if (fFirstEvent) 
+    fEventInspector.ReadProductionDetails(MCEvent());
+  // Clear stuff 
+  fHistos.Clear();
+  fESDFMD.Clear();
+  fAODFMD.Clear();
+  fAODEP.Clear();
+  fMCHistos.Clear();
+  fMCESDFMD.Clear();
+  fMCAODFMD.Clear();
+  fPrimary->Reset();
+  return true;
 }
 //____________________________________________________________________
-void
-AliForwardMCMultiplicityTask::UserExec(Option_t*)
+Bool_t
+AliForwardMCMultiplicityTask::Event(AliESDEvent& esd)
 {
   // 
   // Process each event 
@@ -308,42 +176,30 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   //  
 
   // Read production details 
-  if (fFirstEvent) 
-    fEventInspector.ReadProductionDetails(MCEvent());
     
   // Get the input data 
-  AliESDEvent* esd     = GetESDEvent();
   AliMCEvent*  mcEvent = MCEvent();
-  if (!esd || !mcEvent) return;
-
-  // Clear stuff 
-  fHistos.Clear();
-  fESDFMD.Clear();
-  fAODFMD.Clear();
-  fAODEP.Clear();
-  fMCHistos.Clear();
-  fMCESDFMD.Clear();
-  fMCAODFMD.Clear();
-  fPrimary->Reset();
+  if (!mcEvent) return false;
 
   Bool_t   lowFlux   = kFALSE;
   UInt_t   triggers  = 0;
   UShort_t ivz       = 0;
-  Double_t vz        = 0;
+  TVector3 ip(1024, 1024, 0);
   Double_t cent      = -1;
   UShort_t nClusters = 0;
-  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
-                                              ivz, vz, cent, nClusters);
+  UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux, 
+                                              ivz, ip, cent, nClusters);
   UShort_t ivzMC    = 0;
   Double_t vzMC     = 0;
   Double_t phiR     = 0;
   Double_t b        = 0;
+  Double_t cMC      = 0;
   Int_t    npart    = 0;
   Int_t    nbin     = 0;
   // UInt_t   foundMC  = 
-  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, 
+  fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
                            npart, nbin, phiR);
-  fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
+  fEventInspector.CompareResults(ip.Z(), vzMC, cent, cMC, b, npart, nbin);
   
   //Store all events
   MarkEventForStore();
@@ -367,17 +223,14 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   fMCAODFMD.SetCentrality(cent);
   fMCAODFMD.SetNClusters(nClusters);
   
-  //All events should be stored - HHD
-  //if (isAccepted) MarkEventForStore();
-
   // Disable this check on SPD - will bias data 
   // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
 
   if (isAccepted) {
-    fAODFMD.SetIpZ(vz);
-    fMCAODFMD.SetIpZ(vz);
+    fAODFMD.SetIpZ(ip.Z());
+    fMCAODFMD.SetIpZ(ip.Z());
   }
   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
 
@@ -387,126 +240,98 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*)
   
 
   // Get FMD data 
-  AliESDFMD*  esdFMD  = esd->GetFMDData();
+  AliESDFMD*  esdFMD  = esd.GetFMDData();
 
   // Apply the sharing filter (or hit merging or clustering if you like)
-  if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { 
+  if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD,ip.Z())){
     AliWarning("Sharing filter failed!");
-    return;
+    return false;
   }
-  if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) { 
+  if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, ip.Z(),fMCESDFMD,fPrimary)){
     AliWarning("MC Sharing filter failed!");
-    return;
+    return false;
   }
-  if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
-  // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
+
+  // Store some MC parameters in corners of histogram :-)
+  fPrimary->SetBinContent(0,                      0,                    vzMC);
+  fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,0,                    phiR);
+  fPrimary->SetBinContent(fPrimary->GetNbinsX()+1,fPrimary->GetNbinsY(),cMC);
+  
+
+  if (!isAccepted) 
+    // Exit on MC event w/o trigger, vertex, data - since there's no more 
+    // to be done for MC 
+    return false; 
   
   //MarkEventForStore();
   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
 
   // Calculate the inclusive charged particle density 
-  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) { 
+  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) { 
     AliWarning("Density calculator failed!");
-    return;
+    return false;
   }
   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) { 
     AliWarning("MC Density calculator failed!");
-    return;
+    return false;
   }
   fDensityCalculator.CompareResults(fHistos, fMCHistos);
   
   if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
-    if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
+    if (!fEventPlaneFinder.FindEventplane(&esd, fAODEP, 
+                                         &(fAODFMD.GetHistogram()), &fHistos))
       AliWarning("Eventplane finder failed!");
   }
 
   // Do the secondary and other corrections. 
   if (!fCorrections.Correct(fHistos, ivz)) { 
     AliWarning("Corrections failed");
-    return;
+    return false;
   }
   if (!fCorrections.CorrectMC(fMCHistos, ivz)) { 
     AliWarning("MC Corrections failed");
-    return;
+    return false;
   }
   fCorrections.CompareResults(fHistos, fMCHistos);
     
   if (!fHistCollector.Collect(fHistos, fRingSums, 
-                             ivz, fAODFMD.GetHistogram())) {
+                             ivz, fAODFMD.GetHistogram(),
+                             fAODFMD.GetCentrality())) {
     AliWarning("Histogram collector failed");
-    return;
+    return false;
   }
   if (!fHistCollector.Collect(fMCHistos, fMCRingSums, 
-                             ivz, fMCAODFMD.GetHistogram())) {
+                             ivz, fMCAODFMD.GetHistogram(), -1, true)) {
     AliWarning("MC Histogram collector failed");
-    return;
+    return false;
+  }
+#if 0
+  // Copy underflow bins to overflow bins - always full phi coverage 
+  TH2&  hMC  = fMCAODFMD.GetHistogram();
+  Int_t nEta = hMC.GetNbinsX();
+  Int_t nY   = hMC.GetNbinsY();
+  for (Int_t iEta = 1; iEta <= nEta; iEta++) {
+    hMC.SetBinContent(iEta, nY+1, hMC.GetBinContent(iEta, 0));
   }
+#endif
 
   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
     fHData->Add(&(fAODFMD.GetHistogram()));
 
-  PostData(1, fList);
+  return true;
 }
 
 //____________________________________________________________________
 void
-AliForwardMCMultiplicityTask::Terminate(Option_t*)
+AliForwardMCMultiplicityTask::EstimatedNdeta(const TList* input, 
+                                            TList*       output) const
 {
-  // 
-  // 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;
-  TH1I* hEventsTrVtx = 0;
-  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;
-  }
-
-  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;
-  }
-  
-  // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
-  TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
-  TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
-  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");
-  MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
-
-  fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
-  fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
-  fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
+  AliForwardMultiplicityBase::EstimatedNdeta(input, output);
+  MakeRingdNdeta(input, "mcRingSums", output, "mcRingResults", 24);
 }
 
 
+
 //
 // EOF
 //