PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx

   1 //
   2 // This class collects the event histograms into single histograms,
   3 // one for each ring in each vertex bin.
   4 //
   5 // Input:
   6 //   - AliESDFMD object possibly corrected for sharing
   7 //
   8 // Output:
   9 //   - 5 RingHistos objects - each with a number of vertex dependent
  10 //     2D histograms of the inclusive charge particle density
  11 //
  12 // HistCollector used:
  13 //   - AliFMDCorrSecondaryMap
  14 //
  15 #include "AliFMDHistCollector.h"
  16 #include <AliESDFMD.h>
  17 #include <TAxis.h>
  18 #include <TList.h>
  19 #include <TMath.h>
  20 #include "AliForwardCorrectionManager.h"
  21 #include "AliLog.h"
  22 #include <TH2D.h>
  23 #include <TH1I.h>
  24 #include <TProfile.h>
  25 #include <TProfile2D.h>
  26 #include <TArrayI.h>
  27 #include <TROOT.h>
  28 #include <iostream>
  29 #include <iomanip>
  30
  31 ClassImp(AliFMDHistCollector)
  32 #if 0
  33 ; // For Emacs
  34 #endif
  35
  36
  37 //____________________________________________________________________
  38 AliFMDHistCollector&
  39 AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
  40 {
  41   //
  42   // Assignement operator
  43   //
  44   // Parameters:
  45   //    o Object to assign from
  46   //
  47   // Return:
  48   //    Reference to this object
  49   //
  50   TNamed::operator=(o);
  51
  52   fNCutBins       = o.fNCutBins;
  53   fCorrectionCut  = o.fCorrectionCut;
  54   fFirstBins      = o.fFirstBins;
  55   fLastBins       = o.fLastBins;
  56   fDebug          = o.fDebug;
  57
  58   return *this;
  59 }
  60
  61 //____________________________________________________________________
  62 void
  63 AliFMDHistCollector::Init(const TAxis& vtxAxis)
  64 {
  65   //
  66   // Intialise
  67   //
  68   // Parameters:
  69   //    vtxAxis  Vertex axis
  70   //
  71
  72   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  73
  74   UShort_t nVz = vtxAxis.GetNbins();
  75   fFirstBins.Set(5*nVz);
  76   fLastBins.Set(5*nVz);
  77
  78   // Find the eta bin ranges
  79   for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
  80     // Find the first and last eta bin to use for each ring for
  81     // each vertex bin.   This is instead of using the methods
  82     // provided by AliFMDAnaParameters
  83     for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
  84       UShort_t d = 0;
  85       Char_t   r = 0;
  86       GetDetRing(iIdx, d, r);
  87
  88       // Get the background object
  89       // TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
  90       TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
  91       Int_t nEta  = bg->GetNbinsX();
  92       Int_t first = nEta+1;
  93       Int_t last  = 0;
  94
  95       // Loop over the eta bins
  96       for (Int_t ie = 1; ie <= nEta; ie++) {
  97         if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
  98
  99         // Loop over the phi bins to make sure that we
 100         // do not have holes in the coverage
 101         bool ok = true;
 102         for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
 103           if (bg->GetBinContent(ie,ip) < 0.001) {
 104             ok = false;
 105             continue;
 106           }
 107         }
 108         if (!ok) continue;
 109
 110         first = TMath::Min(ie, first);
 111         last  = TMath::Max(ie, last);
 112       }
 113
 114       // Store the result for later use
 115       fFirstBins[(iVz-1)*5+iIdx] = first;
 116       fLastBins[(iVz-1)*5+iIdx]  = last;
 117     } // for j
 118   }
 119 }
 120
 121 //____________________________________________________________________
 122 Int_t
 123 AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
 124 {
 125   //
 126   // Get the ring index from detector number and ring identifier
 127   //
 128   // Parameters:
 129   //    d Detector
 130   //    r Ring identifier
 131   //
 132   // Return:
 133   //    ring index or -1 in case of problems
 134   //
 135   Int_t idx = -1;
 136   switch (d) {
 137   case 1: idx = 0; break;
 138   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 139   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 140   }
 141   return idx;
 142 }
 143 //____________________________________________________________________
 144 void
 145 AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
 146 {
 147   //
 148   // Get the detector and ring from the ring index
 149   //
 150   // Parameters:
 151   //    idx Ring index
 152   //    d   On return, the detector or 0 in case of errors
 153   //    r   On return, the ring id or '0' in case of errors
 154   //
 155   d = 0;
 156   r = '\0';
 157   switch (idx) {
 158   case 0: d = 1; r = 'I'; break;
 159   case 1: d = 2; r = 'I'; break;
 160   case 2: d = 2; r = 'O'; break;
 161   case 3: d = 3; r = 'I'; break;
 162   case 4: d = 3; r = 'O'; break;
 163   }
 164 }
 165
 166 //____________________________________________________________________
 167 void
 168 AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
 169                                      Int_t& first, Int_t& last) const
 170 {
 171   //
 172   // Get the first and last eta bin to use for a given ring and vertex
 173   //
 174   // Parameters:
 175   //    idx      Ring index as given by GetIdx
 176   //    vtxBin   Vertex bin (1 based)
 177   //    first    On return, the first eta bin to use
 178   //    last     On return, the last eta bin to use
 179   //
 180   first = 0;
 181   last  = 0;
 182
 183   if (idx    <  0) return;
 184   if (vtxbin <= 0) return;
 185   idx += (vtxbin-1) * 5;
 186
 187   if (idx < 0 || idx >= fFirstBins.GetSize()) return;
 188
 189   first = fFirstBins.At(idx)+fNCutBins;
 190   last  = fLastBins.At(idx)-fNCutBins;
 191 }
 192
 193 //____________________________________________________________________
 194 Int_t
 195 AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
 196 {
 197   //
 198   // Get the first eta bin to use for a given ring and vertex
 199   //
 200   // Parameters:
 201   //    idx Ring index as given by GetIdx
 202   //    v vertex bin (1 based)
 203   //
 204   // Return:
 205   //    First eta bin to use, or -1 in case of problems
 206   //
 207   Int_t f, l;
 208   GetFirstAndLast(idx,v,f,l);
 209   return f;
 210 }
 211
 212
 213 //____________________________________________________________________
 214 Int_t
 215 AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
 216 {
 217   //
 218   // Get the last eta bin to use for a given ring and vertex
 219   //
 220   // Parameters:
 221   //    idx Ring index as given by GetIdx
 222   //    v vertex bin (1 based)
 223   //
 224   // Return:
 225   //    Last eta bin to use, or -1 in case of problems
 226   //
 227   Int_t f, l;
 228   GetFirstAndLast(idx,v,f,l);
 229   return l;
 230 }
 231
 232 //____________________________________________________________________
 233 Int_t
 234 AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
 235                                 Int_t bin,  UShort_t vtxbin) const
 236 {
 237   //
 238   // Get the possibly overlapping histogram of eta bin @a e in
 239   // detector and ring
 240   //
 241   // Parameters:
 242   //    d Detector
 243   //    r Ring
 244   //    e Eta bin
 245   //    v Vertex bin (1 based)
 246   //
 247   // Return:
 248   //    Overlapping histogram index or -1
 249   //
 250
 251   Int_t other = -1;
 252   if (d == 1) {
 253     if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
 254   }
 255   else if (d == 2 && r == 'I') {
 256     if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
 257     else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
 258   }
 259   else if (d == 2 && r == 'O') {
 260     if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
 261   }
 262   else if (d == 3 && r == 'O') {
 263     if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
 264   }
 265   else if (d == 3 && r == 'I') {
 266     if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
 267   }
 268   // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
 269   return other;
 270 }
 271 //____________________________________________________________________
 272 Int_t
 273 AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
 274 {
 275   //
 276   // Get the possibly overlapping histogram of eta bin @a e in
 277   // detector and ring
 278   //
 279   // Parameters:
 280   //    i Ring index
 281   //    e Eta bin
 282   //    v Vertex bin (1 based)
 283   //
 284   // Return:
 285   //    Overlapping histogram index or -1
 286   //
 287   UShort_t d = 0;
 288   Char_t   r = '\0';
 289   GetDetRing(idx, d, r);
 290   if (d == 0 || r == '\0') return 0;
 291
 292   return GetOverlap(d, r, bin, vtxbin);
 293 }
 294
 295
 296 //____________________________________________________________________
 297 Bool_t
 298 AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
 299                              UShort_t                vtxbin,
 300                              TH2D&                   out)
 301 {
 302   //
 303   // Do the calculations
 304   //
 305   // Parameters:
 306   //    hists    Cache of histograms
 307   //    vtxBin   Vertex bin (1 based)
 308   //    out      Output histogram
 309   //
 310   // Return:
 311   //    true on successs
 312   //
 313   for (UShort_t d=1; d<=3; d++) {
 314     UShort_t nr = (d == 1 ? 1 : 2);
 315     for (UShort_t q=0; q<nr; q++) {
 316       Char_t      r = (q == 0 ? 'I' : 'O');
 317       TH2D*       h = hists.Get(d,r);
 318       TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
 319
 320       // Outer rings have better phi segmentation - rebin to same as inner.
 321       if (q == 1) t->RebinY(2);
 322
 323       Int_t first = 0;
 324       Int_t last  = 0;
 325       GetFirstAndLast(d, r, vtxbin, first, last);
 326
 327       // Now update profile output
 328       for (Int_t iEta = first; iEta <= last; iEta++) {
 329
 330         // Get the possibly overlapping histogram
 331         Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
 332
 333         // Fill eta acceptance for this event into the phi underlow bin
 334         Float_t ooc      = out.GetBinContent(iEta,0);
 335         Float_t noc      = overlap >= 0? 0.5 : 1;
 336         out.SetBinContent(iEta, 0, ooc + noc);
 337
 338         for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
 339           Double_t c = t->GetBinContent(iEta,iPhi);
 340           Double_t e = t->GetBinError(iEta,iPhi);
 341
 342           // If there's no signal, continue
 343           // if (e <= 0) continue;
 344           if (c <= 0 || e <= 0)     continue;
 345
 346           // If there's no overlapping histogram (ring), then
 347           // fill in data and continue to the next phi bin
 348           if (overlap < 0) {
 349             out.SetBinContent(iEta,iPhi,c);
 350             out.SetBinError(iEta,iPhi,e);
 351             continue;
 352           }
 353
 354           // Get the current bin content and error
 355           Double_t oc = out.GetBinContent(iEta,iPhi);
 356           Double_t oe = out.GetBinError(iEta,iPhi);
 357
 358 #define USE_STRAIGHT_MEAN
 359 // #define USE_STRAIGHT_MEAN_NONZERO
 360 // #define USE_WEIGHTED_MEAN
 361 // #define USE_MOST_CERTAIN
 362 #if defined(USE_STRAIGHT_MEAN)
 363           // calculate the average of old value (half the original),
 364           // and this value, as well as the summed squared errors
 365           // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
 366           // and half the error of this.
 367           //
 368           // So, on the first overlapping histogram we get
 369           //
 370           //    c = c_1 / 2
 371           //    e = sqrt((e_1 / 2)^2) = e_1/2
 372           //
 373           // On the second we get
 374           //
 375           //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
 376           //    e' = sqrt(e^2 + (e_2/2)^2)
 377           //       = sqrt(e_1^2/4 + e_2^2/4)
 378           //       = sqrt(1/4 * (e_1^2+e_2^2))
 379           //       = 1/2 * sqrt(e_1^2 + e_2^2)
 380           out.SetBinContent(iEta,iPhi,oc + c/2);
 381           out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
 382 #elif defined(USE_STRAIGHT_MEAN_NONZERO)
 383 # define ZERO_OTHER
 384           // If there's data in the overlapping histogram,
 385           // calculate the average and add the errors in
 386           // quadrature.
 387           //
 388           // If there's no data in the overlapping histogram,
 389           // then just fill in the data
 390           if (oe > 0) {
 391             out.SetBinContent(iEta,iPhi,(oc + c)/2);
 392             out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
 393           }
 394           else {
 395             out.SetBinContent(iEta,iPhi,c);
 396             out.SetBinError(iEta,iPhi,e);
 397           }
 398 #elif defined(USE_WEIGHTED_MEAN)
 399           // Calculate the weighted mean
 400           Double_t w  = 1/(e*e);
 401           Double_t sc = w * c;
 402           Double_t sw = w;
 403           if (oe > 0) {
 404             Double_t ow = 1/(oe*oe);
 405             sc          += ow * oc;
 406             sw          += ow;
 407           }
 408           Double_t nc = sc / sw;
 409           Double_t ne = TMath::Sqrt(1 / sw);
 410           out.SetBinContent(iEta,iPhi,nc,ne);
 411 #elif defined(USE_MOST_CERTAIN)
 412 # define ZERO_OTHER
 413           if (e < oe) {
 414             out.SetBinContent(iEta,iPhi,c);
 415             out.SetBinError(iEta,iPhi,e);
 416           }
 417           else {
 418             out.SetBinContent(iEta,iPhi,oc);
 419             out.SetBinError(iEta,iPhi,oe);
 420           }
 421 #else
 422 #         error No method for defining content of overlapping bins defined
 423 #endif
 424 #if defined(ZERO_OTHER)
 425           // Get the content of the overlapping histogram,
 426           // and zero the content so that we won't use it
 427           // again
 428           UShort_t od; Char_t oR;
 429           GetDetRing(overlap, od, oR);
 430           TH2D* other = hists.Get(od,oR);
 431           other->SetBinContent(iEta,iPhi,0);
 432           other->SetBinError(iEta,iPhi,0);
 433 #endif
 434         }
 435       }
 436       // Remove temporary histogram
 437       delete t;
 438     } // for r
 439   } // for d
 440   return kTRUE;
 441 }
 442
 443 //____________________________________________________________________
 444 void
 445 AliFMDHistCollector::Print(Option_t* /* option */) const
 446 {
 447   //
 448   // Print information
 449   //
 450   // Parameters:
 451   //    option Not used
 452   //
 453   char ind[gROOT->GetDirLevel()+1];
 454   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 455   ind[gROOT->GetDirLevel()] = '\0';
 456   std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
 457             << ind << " # of cut bins:          " << fNCutBins << '\n'
 458             << ind << " Correction cut:         " << fCorrectionCut << '\n'
 459             << ind << " Bin ranges:\n" << ind << "  v_z bin";
 460   Int_t nVz = fFirstBins.fN / 5;
 461   for (UShort_t iVz = 1; iVz <= nVz; iVz++)
 462     std::cout << " | " << std::setw(7) << iVz;
 463   std::cout << '\n' << ind << "  / ring ";
 464   for (UShort_t iVz = 1; iVz <= nVz; iVz++)
 465     std::cout << "-+--------";
 466   std::cout << std::endl;
 467
 468   for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
 469     UShort_t d = 0;
 470     Char_t   r = 0;
 471     GetDetRing(iIdx, d, r);
 472
 473     std::cout << ind << "  FMD" << d << r << "  ";
 474     for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
 475       Int_t first, last;
 476       GetFirstAndLast(iIdx, iVz, first, last);
 477       std::cout << " | " << std::setw(3) << first << "-"
 478                 << std::setw(3) << last;
 479     }
 480     std::cout << std::endl;
 481   }
 482 }
 483
 484 //____________________________________________________________________
 485 //
 486 // EOF
 487 //
 488
 489
 490