PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx

   1
   2 #include "AliFMDHistCollector.h"
   3 #include <AliESDFMD.h>
   4 #include <TAxis.h>
   5 #include <TList.h>
   6 #include <TMath.h>
   7 #include "AliFMDAnaParameters.h"
   8 #include "AliLog.h"
   9 #include <TH2D.h>
  10 #include <TH1I.h>
  11 #include <TProfile.h>
  12 #include <TProfile2D.h>
  13 #include <TArrayI.h>
  14
  15 ClassImp(AliFMDHistCollector)
  16 #if 0
  17 ; // For Emacs
  18 #endif
  19
  20
  21 //____________________________________________________________________
  22 AliFMDHistCollector&
  23 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
  24 {
  25   TNamed::operator=(o);
  26
  27   fNCutBins       = o.fNCutBins;
  28   fCorrectionCut  = o.fCorrectionCut;
  29   fFirstBins      = o.fFirstBins;
  30   fLastBins       = o.fLastBins;
  31   fDebug          = o.fDebug;
  32
  33   return *this;
  34 }
  35
  36 //____________________________________________________________________
  37 void
  38 AliFMDHistCollector::Init(const TAxis& vtxAxis)
  39 {
  40   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  41
  42   Int_t nVz = vtxAxis.GetNbins();
  43   fFirstBins.Set(5*nVz);
  44   fLastBins.Set(5*nVz);
  45
  46   // Find the eta bin ranges
  47   for (Int_t iVz = 0; iVz < nVz; iVz++) {
  48     // Find the first and last eta bin to use for each ring for
  49     // each vertex bin.   This is instead of using the methods
  50     // provided by AliFMDAnaParameters
  51     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
  52       UShort_t d = 0;
  53       Char_t   r = 0;
  54       GetDetRing(iIdx, d, r);
  55
  56       // Get the background object
  57       TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
  58       Int_t nEta  = bg->GetNbinsX();
  59       Int_t first = nEta+1;
  60       Int_t last  = 0;
  61
  62       // Loop over the eta bins
  63       for (Int_t ie = 1; ie <= nEta; ie++) {
  64         if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
  65
  66         // Loop over the phi bins to make sure that we
  67         // do not have holes in the coverage
  68         bool ok = true;
  69         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
  70           if (bg->GetBinContent(ie,ip) < 0.001) {
  71             ok = false;
  72             continue;
  73           }
  74         }
  75         if (!ok) continue;
  76
  77         first = TMath::Min(ie, first);
  78         last  = TMath::Max(ie, last);
  79       }
  80
  81       // Store the result for later use
  82       fFirstBins[iVz*5+iIdx] = first;
  83       fLastBins[iVz*5+iIdx]  = last;
  84     } // for j
  85   }
  86 }
  87
  88 //____________________________________________________________________
  89 Int_t
  90 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
  91 {
  92   Int_t idx = -1;
  93   switch (d) {
  94   case 1: idx = 0; break;
  95   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  96   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  97   }
  98   return idx;
  99 }
 100 //____________________________________________________________________
 101 void
 102 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
 103 {
 104   d = 0;
 105   r = '\0';
 106   switch (idx) {
 107   case 0: d = 1; r = 'I'; break;
 108   case 1: d = 2; r = 'I'; break;
 109   case 2: d = 2; r = 'O'; break;
 110   case 3: d = 3; r = 'I'; break;
 111   case 4: d = 3; r = 'O'; break;
 112   }
 113 }
 114
 115 //____________________________________________________________________
 116 void
 117 AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin,
 118                                      Int_t& first, Int_t& last) const
 119 {
 120   first = 0;
 121   last  = 0;
 122
 123   if (idx < 0) return;
 124   idx += vtxbin * 5;
 125
 126   if (idx < 0 || idx >= fFirstBins.GetSize()) return;
 127
 128   first = fFirstBins.At(idx)+fNCutBins;
 129   last  = fLastBins.At(idx)-fNCutBins;
 130 }
 131
 132 //____________________________________________________________________
 133 Int_t
 134 AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const
 135 {
 136   Int_t f, l;
 137   GetFirstAndLast(idx,v,f,l);
 138   return f;
 139 }
 140
 141
 142 //____________________________________________________________________
 143 Int_t
 144 AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const
 145 {
 146   Int_t f, l;
 147   GetFirstAndLast(idx,v,f,l);
 148   return l;
 149 }
 150
 151 //____________________________________________________________________
 152 Int_t
 153 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
 154                                 Int_t bin, Int_t vtxbin) const
 155 {
 156   // Get the possibly overlapping histogram
 157   Int_t other = -1;
 158   if (d == 1) {
 159     if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
 160   }
 161   else if (d == 2 && r == 'I') {
 162     if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
 163     else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
 164   }
 165   else if (d == 2 && r == 'O') {
 166     if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
 167   }
 168   else if (d == 3 && r == 'O') {
 169     if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
 170   }
 171   else if (d == 3 && r == 'I') {
 172     if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
 173   }
 174   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
 175   return other;
 176 }
 177 //____________________________________________________________________
 178 Int_t
 179 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
 180 {
 181   UShort_t d = 0;
 182   Char_t   r = '\0';
 183   GetDetRing(idx, d, r);
 184   if (d == 0 || r == '\0') return 0;
 185
 186   return GetOverlap(d, r, bin, vtxbin);
 187 }
 188
 189
 190 //____________________________________________________________________
 191 Bool_t
 192 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
 193                              Int_t                   vtxbin,
 194                              TH2D&                   out)
 195 {
 196   for (UShort_t d=1; d<=3; d++) {
 197     UShort_t nr = (d == 1 ? 1 : 2);
 198     for (UShort_t q=0; q<nr; q++) {
 199       Char_t      r = (q == 0 ? 'I' : 'O');
 200       TH2D*       h = hists.Get(d,r);
 201       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
 202
 203       // Outer rings have better phi segmentation - rebin to same as inner.
 204       if (q == 1) t->RebinY(2);
 205
 206       Int_t first = 0;
 207       Int_t last  = 0;
 208       GetFirstAndLast(d, r, vtxbin, first, last);
 209
 210       // Now update profile output
 211       for (Int_t iEta = first; iEta <= last; iEta++) {
 212
 213         // Get the possibly overlapping histogram
 214         Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
 215
 216         // Fill eta acceptance for this event into the phi underlow bin
 217         Float_t ooc      = out.GetBinContent(iEta,0);
 218         Float_t noc      = overlap >= 0? 0.5 : 1;
 219         out.SetBinContent(iEta, 0, ooc + noc);
 220
 221         for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
 222           Double_t c = t->GetBinContent(iEta,iPhi);
 223           Double_t e = t->GetBinError(iEta,iPhi);
 224
 225           // If there's no signal, continue
 226           // if (e <= 0) continue;
 227           if (c <= 0 || e <= 0)     continue;
 228
 229           // If there's no overlapping histogram (ring), then
 230           // fill in data and continue to the next phi bin
 231           if (overlap < 0) {
 232             out.SetBinContent(iEta,iPhi,c);
 233             out.SetBinError(iEta,iPhi,e);
 234             continue;
 235           }
 236
 237           // Get the current bin content and error
 238           Double_t oc = out.GetBinContent(iEta,iPhi);
 239           Double_t oe = out.GetBinError(iEta,iPhi);
 240
 241 #define USE_STRAIGHT_MEAN
 242 // #define USE_STRAIGHT_MEAN_NONZERO
 243 // #define USE_WEIGHTED_MEAN
 244 // #define USE_MOST_CERTAIN
 245 #if defined(USE_STRAIGHT_MEAN)
 246           // calculate the average of old value (half the original),
 247           // and this value, as well as the summed squared errors
 248           // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
 249           // and half the error of this.
 250           //
 251           // So, on the first overlapping histogram we get
 252           //
 253           //    c = c_1 / 2
 254           //    e = sqrt((e_1 / 2)^2) = e_1/2
 255           //
 256           // On the second we get
 257           //
 258           //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
 259           //    e' = sqrt(e^2 + (e_2/2)^2)
 260           //       = sqrt(e_1^2/4 + e_2^2/4)
 261           //       = sqrt(1/4 * (e_1^2+e_2^2))
 262           //       = 1/2 * sqrt(e_1^2 + e_2^2)
 263           out.SetBinContent(iEta,iPhi,oc + c/2);
 264           out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
 265 #elif defined(USE_STRAIGHT_MEAN_NONZERO)
 266 # define ZERO_OTHER
 267           // If there's data in the overlapping histogram,
 268           // calculate the average and add the errors in
 269           // quadrature.
 270           //
 271           // If there's no data in the overlapping histogram,
 272           // then just fill in the data
 273           if (oe > 0) {
 274             out.SetBinContent(iEta,iPhi,(oc + c)/2);
 275             out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
 276           }
 277           else {
 278             out.SetBinContent(iEta,iPhi,c);
 279             out.SetBinError(iEta,iPhi,e);
 280           }
 281 #elif defined(USE_WEIGHTED_MEAN)
 282           // Calculate the weighted mean
 283           Double_t w  = 1/(e*e);
 284           Double_t sc = w * c;
 285           Double_t sw = w;
 286           if (oe > 0) {
 287             Double_t ow = 1/(oe*oe);
 288             sc          += ow * oc;
 289             sw          += ow;
 290           }
 291           Double_t nc = sc / sw;
 292           Double_t ne = TMath::Sqrt(1 / sw);
 293           out.SetBinContent(iEta,iPhi,nc,ne);
 294 #elif defined(USE_MOST_CERTAIN)
 295 # define ZERO_OTHER
 296           if (e < oe) {
 297             out.SetBinContent(iEta,iPhi,c);
 298             out.SetBinError(iEta,iPhi,e);
 299           }
 300           else {
 301             out.SetBinContent(iEta,iPhi,oc);
 302             out.SetBinError(iEta,iPhi,oe);
 303           }
 304 #else
 305 #         error No method for defining content of overlapping bins defined
 306 #endif
 307 #if defined(ZERO_OTHER)
 308           // Get the content of the overlapping histogram,
 309           // and zero the content so that we won't use it
 310           // again
 311           UShort_t od; Char_t oR;
 312           GetDetRing(overlap, od, oR);
 313           TH2D* other = hists.Get(od,oR);
 314           other->SetBinContent(iEta,iPhi,0);
 315           other->SetBinError(iEta,iPhi,0);
 316 #endif
 317         }
 318       }
 319       // Remove temporary histogram
 320       delete t;
 321     } // for r
 322   } // for d
 323   return kTRUE;
 324 }
 325
 326
 327 //____________________________________________________________________
 328 //
 329 // EOF
 330 //
 331
 332
 333