]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliFMDCorrector.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDCorrector.cxx
index 4a6b3fee7c8bc2b32086591949286d96dbcba5bf..aa29928058ba2d044e3b586b4f877472f9e16700 100644 (file)
@@ -8,7 +8,7 @@
 #include <TList.h>
 #include <TMath.h>
 #include "AliForwardCorrectionManager.h"
-// #include "AliFMDCorrDoubleHit.h"
+#include "AliFMDCorrSecondaryMap.h"
 #include "AliFMDCorrVertexBias.h"
 #include "AliFMDCorrMergingEfficiency.h"
 #include "AliFMDCorrAcceptance.h"
@@ -30,12 +30,12 @@ AliFMDCorrector::AliFMDCorrector()
     fRingHistos(),
     fUseSecondaryMap(true),
     fUseVertexBias(false),
-    fUseAcceptance(true),
+    fUseAcceptance(false),
     fUseMergingEfficiency(false),
     fDebug(0)
 {
   // Constructor
-  DGUARD(fDebug, 0, "Default CTOR of AliFMDCorrector");
+  DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
 }
 
 //____________________________________________________________________
@@ -44,7 +44,7 @@ AliFMDCorrector::AliFMDCorrector(const char* title)
     fRingHistos(), 
     fUseSecondaryMap(true),
     fUseVertexBias(false),
-    fUseAcceptance(true),
+    fUseAcceptance(false),
     fUseMergingEfficiency(false),
     fDebug(0)
 {
@@ -52,7 +52,7 @@ AliFMDCorrector::AliFMDCorrector(const char* title)
   // 
   // Parameters: 
   //   title   Title
-  DGUARD(fDebug, 0, "Named CTOR of AliFMDCorrector: %s", title);
+  DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
   fRingHistos.SetName(GetName());
   fRingHistos.Add(new RingHistos(1, 'I'));
   fRingHistos.Add(new RingHistos(2, 'I'));
@@ -75,7 +75,7 @@ AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
   // 
   // Parameters: 
   //  o  Object to copy from 
-  DGUARD(fDebug, 0, "Copy CTOR of AliFMDCorrector");
+  DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
   TIter    next(&o.fRingHistos);
   TObject* obj = 0;
   while ((obj = next())) fRingHistos.Add(obj);
@@ -118,7 +118,7 @@ AliFMDCorrector::operator=(const AliFMDCorrector& o)
 
 //____________________________________________________________________
 void
-AliFMDCorrector::Init(const TAxis&)
+AliFMDCorrector::SetupForData(const TAxis&)
 {
   //
   // Initialize this object
@@ -159,6 +159,43 @@ AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
   return static_cast<RingHistos*>(fRingHistos.At(idx));
 }
     
+//____________________________________________________________________
+void
+AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
+                          Bool_t alsoUnderOver) const
+{
+  // 
+  // Implement TH1::Divide but 
+  // - Assume compatible histograms 
+  // - Unless third argument is true, do not divide over/under flow bins
+  // 
+  if (!num || !denom) return;
+
+  Int_t first = (alsoUnderOver ? 0 : 1);
+  Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
+  Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
+  
+  for (Int_t ix = first; ix <= lastX; ix++) {
+    for (Int_t iy = first; iy <= lastY; iy++) { 
+      Int_t    bin = num->GetBin(ix,iy);
+      Double_t c0  = num->GetBinContent(bin);
+      Double_t c1  = denom->GetBinContent(bin);
+      if (!c1) { 
+       num->SetBinContent(bin,0);
+       num->SetBinError(bin, 0);
+       continue;
+      }
+      Double_t w   = c0 / c1;
+      Double_t e0  = num->GetBinError(bin);
+      Double_t e1  = denom->GetBinError(bin);
+      Double_t c12 = c1*c1;
+      Double_t e2  = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
+      
+      num->SetBinContent(bin, w);
+      num->SetBinError(bin, TMath::Sqrt(e2));
+    }
+  }
+}
 //____________________________________________________________________
 Bool_t
 AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
@@ -192,7 +229,7 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
           continue;
         }
         // Divide by primary/total ratio
-        h->Divide(bg);
+       DivideMap(h, bg, false);
       }
       if (fUseVertexBias) {
         TH2D*  ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
@@ -202,7 +239,7 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
           continue;
         }
         // Divide by the event selection efficiency
-        h->Divide(ef);
+       DivideMap(h, ef, false);
       }
       if (fUseAcceptance) {
         TH2D*  ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
@@ -211,8 +248,12 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
                          "vertex bin %d", d, r, uvb));
           continue;
         }
-        // Divide by the acceptance correction
-        h->Divide(ac);
+       // Fill overflow bin with ones 
+       for (Int_t i = 1; i <= h->GetNbinsX(); i++) 
+         h->SetBinContent(i, h->GetNbinsY()+1, 1);
+
+        // Divide by the acceptance correction - 
+       DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
       }
 
       if (fUseMergingEfficiency) {
@@ -231,10 +272,14 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
        for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
          Float_t c  = sf->GetBinContent(ieta);
          Float_t ec = sf->GetBinError(ieta);
-         
-         if (c == 0) continue;
-         
+                 
          for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
+           if (c == 0) {
+             h->SetBinContent(ieta,iphi,0);
+             h->SetBinError(ieta,iphi,0);
+             continue;
+           }
+
            Double_t m  = h->GetBinContent(ieta, iphi) / c;
            Double_t em = h->GetBinError(ieta, iphi);
          
@@ -245,49 +290,6 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
          }
        }
       }
-      //HHD
-      /*
-      TH2D*  bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
-      TH2D   hRing("hring","hring",bg->GetNbinsX(),
-                  bg->GetXaxis()->GetXmin(),
-                  bg->GetXaxis()->GetXmax(),
-                  bg->GetNbinsY(),
-                  bg->GetYaxis()->GetXmin(),
-                  bg->GetYaxis()->GetXmax());
-      
-      Int_t edgebin[4] = {0,0,0,0};
-      for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
-       for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
-         Float_t bgcor = bg->GetBinContent(ii,jj);
-         if(bgcor<0.1) continue;
-         if(edgebin[0] == 0) edgebin[0] = ii;
-         if(edgebin[0] == ii) continue;
-         if(edgebin[0] > 0 && edgebin[1] == 0) edgebin[1] = ii;
-         if(edgebin[0]>0 && edgebin[1]>0) break; 
-       }
-      }
-      for(Int_t ii = bg->GetNbinsX(); ii >= 1;  ii--) {
-       for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
-         Float_t bgcor = bg->GetBinContent(ii,jj);
-         if(bgcor<0.1) continue;
-         if(edgebin[2] == 0) edgebin[2] = ii;
-         if(edgebin[2] == ii) continue;
-         if(edgebin[2] > 0 && edgebin[3] == 0) edgebin[3] = ii;
-         if(edgebin[2]>0 && edgebin[3]>0) break; 
-       }
-      }
-      for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
-       for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
-         Float_t data = h->GetBinContent(ii,jj);
-         if(data <0.000001) continue;
-         if(edgebin[0] == ii || edgebin[1] == ii || edgebin[2] == ii || edgebin[3] == ii) continue;
-         hRing.SetBinContent(ii,jj,data);
-         hRing.SetBinError(ii,jj,h->GetBinError(ii,jj));
-       }
-      }
-      
-      //std::cout<<edgebin[0]<<"  "<<edgebin[1]<<"  "<<edgebin[2]<<"   "<<edgebin[3]<<std::endl;
-      */
       
       rh->fDensity->Add(h);
     }
@@ -298,7 +300,7 @@ AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
 
 //____________________________________________________________________
 void
-AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
+AliFMDCorrector::Terminate(const TList* dir, TList* output, Int_t nEvents)
 {
   // 
   // Scale the histograms to the total number of events 
@@ -311,11 +313,19 @@ AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
   if (!d) return;
 
+  TList* out = new TList;
+  out->SetName(d->GetName());
+  out->SetOwner();
+
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   THStack* sums = new THStack("sums", "Sums of ring results");
   while ((o = static_cast<RingHistos*>(next()))) {
-    o->ScaleHistograms(d, nEvents);
+    o->Terminate(d, nEvents);
+    if (!o->fDensity) { 
+      Warning("Terminate", "No density from %s", o->GetName());
+      continue;
+    }
     TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, 
                                         o->fDensity->GetNbinsY(),"e");
     sum->Scale(1., "width");
@@ -324,12 +334,12 @@ AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
     sum->SetYTitle("#sum N_{ch,primary}");
     sums->Add(sum);
   }
-  d->Add(sums);
-
+  out->Add(sums);
+  output->Add(out);
 }
 //____________________________________________________________________
 void
-AliFMDCorrector::DefineOutput(TList* dir)
+AliFMDCorrector::CreateOutputObjects(TList* dir)
 {
   // 
   // Output diagnostic histograms to directory 
@@ -351,10 +361,22 @@ AliFMDCorrector::DefineOutput(TList* dir)
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next()))) {
-    o->Output(d);
+    o->CreateOutputObjects(d);
   }
 }
 
+#define PF(N,V,...)                                    \
+  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
+#define PFB(N,FLAG)                            \
+  do {                                                                 \
+    AliForwardUtil::PrintName(N);                                      \
+    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
+  } while(false)
+#define PFV(N,VALUE)                                   \
+  do {                                                 \
+    AliForwardUtil::PrintName(N);                      \
+    std::cout << (VALUE) << std::endl; } while(false)
+
 //____________________________________________________________________
 void
 AliFMDCorrector::Print(Option_t* /* option */) const
@@ -364,16 +386,13 @@ AliFMDCorrector::Print(Option_t* /* option */) const
   // Parameters:
   //    option Not used 
   //
-  char ind[gROOT->GetDirLevel()+1];
-  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
-  ind[gROOT->GetDirLevel()] = '\0';
-  std::cout << ind << ClassName() << ": " << GetName() <<  "\n"
-            << std::boolalpha
-            << ind << " Use secondary maps:     " << fUseSecondaryMap << "\n"
-            << ind << " Use vertex bias:        " << fUseVertexBias << "\n"
-            << ind << " Use acceptance:         " << fUseAcceptance << "\n"
-            << ind << " Use merging efficiency: " << fUseMergingEfficiency
-            << std::noboolalpha << std::endl;
+  AliForwardUtil::PrintTask(*this);
+  gROOT->IncreaseDirLevel();
+  PFB("Use secondary maps",            fUseSecondaryMap );
+  PFB("Use vertex bias",               fUseVertexBias );
+  PFB("Use acceptance",                        fUseAcceptance );
+  PFB("Use merging efficiency",                fUseMergingEfficiency);
+  gROOT->DecreaseDirLevel();  
 }
 
 //====================================================================
@@ -447,12 +466,12 @@ AliFMDCorrector::RingHistos::~RingHistos()
   // 
   // Destructor 
   //
-  if (fDensity) delete fDensity;
+  // if (fDensity) delete fDensity;
 }
 
 //____________________________________________________________________
 void
-AliFMDCorrector::RingHistos::Output(TList* dir)
+AliFMDCorrector::RingHistos::CreateOutputObjects(TList* dir)
 {
   // 
   // Make output 
@@ -465,7 +484,7 @@ AliFMDCorrector::RingHistos::Output(TList* dir)
 
 //____________________________________________________________________
 void
-AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
+AliFMDCorrector::RingHistos::Terminate(TList* dir, Int_t nEvents)
 { 
   // 
   // Scale the histograms to the total number of events 
@@ -476,8 +495,9 @@ AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
   TList* l = GetOutputList(dir);
   if (!l) return; 
 
-  TH1* density = GetOutputHist(l,"primaryDensity");
+  TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
   if (density) density->Scale(1./nEvents);
+  fDensity = density;
 }
 
 //____________________________________________________________________