PWG2/FORWARD/analysis2/AliFMDHistCollector.cxx

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