PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx

   1 // This class calculates the inclusive charged particle density
   2 // in each for the 5 FMD rings.
   3 //
   4 #include "AliFMDDensityCalculator.h"
   5 #include <AliESDFMD.h>
   6 #include <TAxis.h>
   7 #include <TList.h>
   8 #include <TMath.h>
   9 #include "AliForwardCorrectionManager.h"
  10 #include "AliFMDCorrDoubleHit.h"
  11 #include "AliFMDCorrELossFit.h"
  12 #include "AliLog.h"
  13 #include "AliForwardUtil.h"
  14 #include <TH2D.h>
  15 #include <TProfile.h>
  16 #include <THStack.h>
  17 #include <TROOT.h>
  18 #include <TVector3.h>
  19 #include <TStopwatch.h>
  20 #include <TParameter.h>
  21 #include <iostream>
  22 #include <iomanip>
  23 #include <cstring>
  24
  25 ClassImp(AliFMDDensityCalculator)
  26 #if 0
  27 ; // For Emacs
  28 #endif
  29
  30 //____________________________________________________________________
  31 const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator";
  32
  33 //____________________________________________________________________
  34 AliFMDDensityCalculator::AliFMDDensityCalculator()
  35   : TNamed(),
  36     fRingHistos(),
  37     fSumOfWeights(0),
  38     fWeightedSum(0),
  39     fCorrections(0),
  40     fMaxParticles(5),
  41     fUsePoisson(false),
  42     fUsePhiAcceptance(kPhiCorrectNch),
  43     fAccI(0),
  44     fAccO(0),
  45     fFMD1iMax(0),
  46     fFMD2iMax(0),
  47     fFMD2oMax(0),
  48     fFMD3iMax(0),
  49     fFMD3oMax(0),
  50     fMaxWeights(0),
  51     fLowCuts(0),
  52     fEtaLumping(32),
  53     fPhiLumping(4),
  54     fDebug(0),
  55     fCuts(),
  56     fRecalculateEta(false),
  57     fRecalculatePhi(false),
  58     fMinQuality(10),
  59     fCache(),
  60     fDoTiming(false),
  61   fHTiming(0),
  62   fMaxOutliers(0.05),
  63   fOutlierCut(0.50)
  64 {
  65   //
  66   // Constructor
  67   //
  68   DGUARD(fDebug, 3, "Default CTOR of FMD density calculator");
  69 }
  70
  71 //____________________________________________________________________
  72 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
  73   : TNamed(fgkFolderName, title),
  74     fRingHistos(),
  75     fSumOfWeights(0),
  76     fWeightedSum(0),
  77     fCorrections(0),
  78     fMaxParticles(5),
  79     fUsePoisson(false),
  80     fUsePhiAcceptance(kPhiCorrectNch),
  81     fAccI(0),
  82     fAccO(0),
  83     fFMD1iMax(0),
  84     fFMD2iMax(0),
  85     fFMD2oMax(0),
  86     fFMD3iMax(0),
  87     fFMD3oMax(0),
  88     fMaxWeights(0),
  89     fLowCuts(0),
  90     fEtaLumping(32),
  91     fPhiLumping(4),
  92     fDebug(0),
  93     fCuts(),
  94     fRecalculateEta(false),
  95     fRecalculatePhi(false),
  96     fMinQuality(10),
  97     fCache(),
  98     fDoTiming(false),
  99     fHTiming(0),
 100   fMaxOutliers(0.05),
 101   fOutlierCut(0.50)
 102 {
 103   //
 104   // Constructor
 105   //
 106   // Parameters:
 107   //    name Name of object
 108   //
 109   DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title);
 110   fRingHistos.SetName(GetName());
 111   fRingHistos.SetOwner();
 112   fRingHistos.Add(new RingHistos(1, 'I'));
 113   fRingHistos.Add(new RingHistos(2, 'I'));
 114   fRingHistos.Add(new RingHistos(2, 'O'));
 115   fRingHistos.Add(new RingHistos(3, 'I'));
 116   fRingHistos.Add(new RingHistos(3, 'O'));
 117   fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
 118                            200, 0, 20);
 119   fSumOfWeights->SetFillColor(kRed+1);
 120   fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
 121   fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
 122                            200, 0, 20);
 123   fWeightedSum->SetFillColor(kBlue+1);
 124   fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
 125   fCorrections  = new TH1D("corrections", "Distribution of corrections",
 126                            100, 0, 10);
 127   fCorrections->SetFillColor(kBlue+1);
 128   fCorrections->SetXTitle("correction");
 129
 130   fAccI = GenerateAcceptanceCorrection('I');
 131   fAccO = GenerateAcceptanceCorrection('O');
 132
 133   fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
 134                          1, 0, 1, 1, 0, 1);
 135   fMaxWeights->SetXTitle("#eta");
 136   fMaxWeights->SetDirectory(0);
 137
 138   fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
 139   fLowCuts->SetXTitle("#eta");
 140   fLowCuts->SetDirectory(0);
 141
 142 }
 143
 144 //____________________________________________________________________
 145 AliFMDDensityCalculator::AliFMDDensityCalculator(const
 146                                                  AliFMDDensityCalculator& o)
 147   : TNamed(o),
 148     fRingHistos(),
 149     fSumOfWeights(o.fSumOfWeights),
 150     fWeightedSum(o.fWeightedSum),
 151     fCorrections(o.fCorrections),
 152     fMaxParticles(o.fMaxParticles),
 153     fUsePoisson(o.fUsePoisson),
 154     fUsePhiAcceptance(o.fUsePhiAcceptance),
 155     fAccI(o.fAccI),
 156     fAccO(o.fAccO),
 157     fFMD1iMax(o.fFMD1iMax),
 158     fFMD2iMax(o.fFMD2iMax),
 159     fFMD2oMax(o.fFMD2oMax),
 160     fFMD3iMax(o.fFMD3iMax),
 161     fFMD3oMax(o.fFMD3oMax),
 162     fMaxWeights(o.fMaxWeights),
 163     fLowCuts(o.fLowCuts),
 164     fEtaLumping(o.fEtaLumping),
 165     fPhiLumping(o.fPhiLumping),
 166     fDebug(o.fDebug),
 167     fCuts(o.fCuts),
 168     fRecalculateEta(o.fRecalculateEta),
 169     fRecalculatePhi(o.fRecalculatePhi),
 170     fMinQuality(o.fMinQuality),
 171     fCache(o.fCache),
 172     fDoTiming(o.fDoTiming),
 173     fHTiming(o.fHTiming),
 174   fMaxOutliers(o.fMaxOutliers),
 175   fOutlierCut(o.fOutlierCut)
 176 {
 177   //
 178   // Copy constructor
 179   //
 180   // Parameters:
 181   //    o Object to copy from
 182   //
 183   DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator");
 184   TIter    next(&o.fRingHistos);
 185   TObject* obj = 0;
 186   while ((obj = next())) fRingHistos.Add(obj);
 187 }
 188
 189 //____________________________________________________________________
 190 AliFMDDensityCalculator::~AliFMDDensityCalculator()
 191 {
 192   //
 193   // Destructor
 194   //
 195   DGUARD(fDebug, 3, "DTOR of FMD density calculator");
 196   // fRingHistos.Delete();
 197 }
 198
 199 //____________________________________________________________________
 200 AliFMDDensityCalculator&
 201 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
 202 {
 203   //
 204   // Assignement operator
 205   //
 206   // Parameters:
 207   //    o Object to assign from
 208   //
 209   // Return:
 210   //    Reference to this object
 211   //
 212   DGUARD(fDebug, 3, "Assignment of FMD density calculator");
 213   if (&o == this) return *this;
 214   TNamed::operator=(o);
 215
 216   fDebug              = o.fDebug;
 217   fMaxParticles       = o.fMaxParticles;
 218   fUsePoisson         = o.fUsePoisson;
 219   fUsePhiAcceptance   = o.fUsePhiAcceptance;
 220   fAccI               = o.fAccI;
 221   fAccO               = o.fAccO;
 222   fFMD1iMax           = o.fFMD1iMax;
 223   fFMD2iMax           = o.fFMD2iMax;
 224   fFMD2oMax           = o.fFMD2oMax;
 225   fFMD3iMax           = o.fFMD3iMax;
 226   fFMD3oMax           = o.fFMD3oMax;
 227   fMaxWeights         = o.fMaxWeights;
 228   fLowCuts            = o.fLowCuts;
 229   fEtaLumping         = o.fEtaLumping;
 230   fPhiLumping         = o.fPhiLumping;
 231   fCuts               = o.fCuts;
 232   fRecalculateEta     = o.fRecalculateEta;
 233   fRecalculatePhi     = o.fRecalculatePhi;
 234   fMinQuality         = o.fMinQuality;
 235   fCache              = o.fCache;
 236   fDoTiming           = o.fDoTiming;
 237   fHTiming            = o.fHTiming;
 238   fMaxOutliers        = o.fMaxOutliers;
 239   fOutlierCut         = o.fOutlierCut;
 240
 241   fRingHistos.Delete();
 242   TIter    next(&o.fRingHistos);
 243   TObject* obj = 0;
 244   while ((obj = next())) fRingHistos.Add(obj);
 245
 246   return *this;
 247 }
 248
 249 //____________________________________________________________________
 250 void
 251 AliFMDDensityCalculator::SetupForData(const TAxis& axis)
 252 {
 253   // Intialize this sub-algorithm
 254   //
 255   // Parameters:
 256   //   etaAxis   Eta axis
 257   DGUARD(fDebug, 1, "Initialize FMD density calculator");
 258   CacheMaxWeights(axis);
 259
 260   fCache.Init(axis);
 261
 262   TIter    next(&fRingHistos);
 263   RingHistos* o = 0;
 264   while ((o = static_cast<RingHistos*>(next()))) {
 265     o->SetupForData(axis);
 266     // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
 267     // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
 268   }
 269 }
 270
 271 //____________________________________________________________________
 272 AliFMDDensityCalculator::RingHistos*
 273 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
 274 {
 275   //
 276   // Get the ring histogram container
 277   //
 278   // Parameters:
 279   //    d Detector
 280   //    r Ring
 281   //
 282   // Return:
 283   //    Ring histogram container
 284   //
 285   Int_t idx = -1;
 286   switch (d) {
 287   case 1: idx = 0; break;
 288   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 289   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 290   }
 291   if (idx < 0 || idx >= fRingHistos.GetEntries()) {
 292     AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
 293     return 0;
 294   }
 295
 296   return static_cast<RingHistos*>(fRingHistos.At(idx));
 297 }
 298
 299 //____________________________________________________________________
 300 Double_t
 301 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
 302                                     Bool_t errors) const
 303 {
 304   //
 305   // Get the multiplicity cut.  If the user has set fMultCut (via
 306   // SetMultCut) then that value is used.  If not, then the lower
 307   // value of the fit range for the enery loss fits is returned.
 308   //
 309   // Return:
 310   //    Lower cut on multiplicity
 311   //
 312   return fCuts.GetMultCut(d,r,ieta,errors);
 313 }
 314
 315 //____________________________________________________________________
 316 Double_t
 317 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
 318                                     Bool_t errors) const
 319 {
 320   //
 321   // Get the multiplicity cut.  If the user has set fMultCut (via
 322   // SetMultCut) then that value is used.  If not, then the lower
 323   // value of the fit range for the enery loss fits is returned.
 324   //
 325   // Return:
 326   //    Lower cut on multiplicity
 327   //
 328   return fCuts.GetMultCut(d,r,eta,errors);
 329 }
 330
 331 //____________________________________________________________________
 332 Bool_t
 333 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
 334                                    AliForwardUtil::Histos& hists,
 335                                    Bool_t                  lowFlux,
 336                                    Double_t                /*cent*/,
 337                                    const TVector3&         ip)
 338 {
 339   //
 340   // Do the calculations
 341   //
 342   // Parameters:
 343   //    fmd      AliESDFMD object (possibly) corrected for sharing
 344   //    hists    Histogram cache
 345   //    vtxBin   Vertex bin
 346   //    lowFlux  Low flux flag.
 347   //
 348   // Return:
 349   //    true on successs
 350   DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
 351
 352   TStopwatch timer;
 353   TStopwatch totalT;
 354
 355   // First measurements of timing
 356   //  Re-calculation      : fraction of sum 32.0%   of total 18.1%
 357   //  N_{particle}        : fraction of sum 15.2%   of total  8.6%
 358   //  Correction          : fraction of sum 26.4%   of total 14.9%
 359   //  #phi acceptance     : fraction of sum  0.2%   of total  0.1%
 360   //  Copy to cache       : fraction of sum  3.9%   of total  2.2%
 361   //  Poisson calculation : fraction of sum 18.7%   of total 10.6%
 362   //  Diagnostics         : fraction of sum  3.7%   of total  2.1%
 363   Double_t reEtaTime  = 0;
 364   Double_t nPartTime  = 0;
 365   Double_t corrTime   = 0;
 366   Double_t rePhiTime  = 0;
 367   Double_t copyTime   = 0;
 368   Double_t poissonTime= 0;
 369   Double_t diagTime   = 0;
 370   if (fDoTiming) totalT.Start(true);
 371
 372   Double_t etaCache[20*512]; // Same number of strips per ring
 373   Double_t phiCache[20*512]; // whether it is inner our outer.
 374   // We do not use TArrayD because we do not wont a bounds check
 375   // TArrayD etaCache(20*512); // Same number of strips per ring
 376   // TArrayD phiCache(20*512); // whether it is inner our outer.
 377
 378   // --- Loop over detectors -----------------------------------------
 379   for (UShort_t d=1; d<=3; d++) {
 380     UShort_t nr = (d == 1 ? 1 : 2);
 381     for (UShort_t q=0; q<nr; q++) {
 382       Char_t      r = (q == 0 ? 'I' : 'O');
 383       UShort_t    ns= (q == 0 ?  20 :  40);
 384       UShort_t    nt= (q == 0 ? 512 : 256);
 385       TH2D*       h = hists.Get(d,r);
 386       RingHistos* rh= GetRingHistos(d,r);
 387       if (!rh) {
 388         AliError(Form("No ring histogram found for FMD%d%c", d, r));
 389         fRingHistos.ls();
 390         return false;
 391       }
 392       // rh->fPoisson.SetObject(d,r,vtxbin,cent);
 393       rh->fPoisson.Reset(0);
 394       rh->fTotal->Reset();
 395       rh->fGood->Reset();
 396       // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
 397
 398       // Reset our eta cache
 399       // for (Int_t i = 0; i < 20*512; i++)
 400       // etaCache[i] = phiCache[i] = AliESDFMD::kInvalidEta;
 401       memset(etaCache, 0, sizeof(Double_t)*20*512);
 402       memset(phiCache, 0, sizeof(Double_t)*20*512);
 403       // etaCache.Reset(AliESDFMD::kInvalidEta);
 404       // phiCache.Reset(AliESDFMD::kInvalidEta);
 405
 406       // --- Loop over sectors and strips ----------------------------
 407       for (UShort_t s=0; s<ns; s++) {
 408         for (UShort_t t=0; t<nt; t++) {
 409
 410           Float_t  mult   = fmd.Multiplicity(d,r,s,t);
 411           Float_t  phi    = fmd.Phi(d,r,s,t) * TMath::DegToRad();
 412           Float_t  eta    = fmd.Eta(d,r,s,t);
 413           Double_t oldPhi = phi;
 414
 415           // --- Re-calculate eta - needed for satelittes ------------
 416           if (fDoTiming) timer.Start(true);
 417           if (eta == AliESDFMD::kInvalidEta || fRecalculateEta)
 418             eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
 419           if (fDoTiming) reEtaTime += timer.CpuTime();
 420           etaCache[s*nt+t] = eta;
 421
 422           // --- Check this strip ------------------------------------
 423           rh->fTotal->Fill(eta);
 424           if (mult == AliESDFMD::kInvalidMult) { //  || mult > 20) {
 425             // Do not count invalid stuff
 426             rh->fELoss->Fill(-1);
 427             // rh->fEvsN->Fill(mult,-1);
 428             // rh->fEvsM->Fill(mult,-1);
 429             continue;
 430           }
 431           if (mult > 20)
 432             AliWarningF("Raw multiplicity of FMD%d%c[%02d,%03d] = %f > 20",
 433                         d, r, s, t, mult);
 434           // --- Automatic calculation of acceptance -----------------
 435           rh->fGood->Fill(eta);
 436
 437           // --- If we asked to re-calculate phi for (x,y) IP --------
 438           if (fDoTiming) timer.Start(true);
 439           if (fRecalculatePhi) {
 440             oldPhi = phi;
 441             phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
 442           }
 443           phiCache[s*nt+t] = phi;
 444           if (fDoTiming) rePhiTime += timer.CpuTime();
 445
 446           // --- Apply phi corner correction to eloss ----------------
 447           if (fUsePhiAcceptance == kPhiCorrectELoss)
 448             mult *= AcceptanceCorrection(r,t);
 449
 450           // --- Get the low multiplicity cut ------------------------
 451           Double_t cut  = 1024;
 452           if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
 453           else AliWarningF("Eta for FMD%d%c[%02d,%03d] is invalid: %f",
 454                            d, r, s, t, eta);
 455
 456           // --- Now caluculate Nch for this strip using fits --------
 457           if (fDoTiming) timer.Start(true);
 458           Double_t n   = 0;
 459           if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
 460           rh->fELoss->Fill(mult);
 461           // rh->fEvsN->Fill(mult,n);
 462           // rh->fEtaVsN->Fill(eta, n);
 463           if (fDoTiming) nPartTime += timer.CpuTime();
 464
 465           // --- Calculate correction if needed ----------------------
 466           if (fDoTiming) timer.Start(true);
 467           // Temporary stuff - remove Correction call
 468           Double_t c = 1;
 469           if (fUsePhiAcceptance == kPhiCorrectNch)
 470             c = AcceptanceCorrection(r,t);
 471           // Double_t c = Correction(d,r,t,eta,lowFlux);
 472           if (fDoTiming) corrTime += timer.CpuTime();
 473           fCorrections->Fill(c);
 474           if (c > 0) n /= c;
 475           // rh->fEvsM->Fill(mult,n);
 476           // rh->fEtaVsM->Fill(eta, n);
 477           rh->fCorr->Fill(eta, c);
 478
 479           // --- Accumulate Poisson statistics -----------------------
 480           Bool_t hit = (n > 0.9 && c > 0);
 481           if (hit) {
 482             rh->fELossUsed->Fill(mult);
 483             if (fRecalculatePhi) {
 484               rh->fPhiBefore->Fill(oldPhi);
 485               rh->fPhiAfter->Fill(phi);
 486             }
 487           }
 488           rh->fPoisson.Fill(t,s,hit,1./c);
 489           h->Fill(eta,phi,n);
 490
 491           // --- If we use ELoss fits, apply now ---------------------
 492           if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
 493         } // for t
 494       } // for s
 495
 496       // --- Automatic acceptance - Calculate as an efficiency -------
 497       // This is very fast, so we do not bother to time it
 498       rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
 499
 500       // --- Make a copy and reset as needed -------------------------
 501       if (fDoTiming) timer.Start(true);
 502       TH2D* hclone = fCache.Get(d,r);
 503       // hclone->Reset();
 504       // TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
 505       if (!fUsePoisson) hclone->Reset();
 506       else {
 507         for (Int_t i = 0; i <= h->GetNbinsX()+1; i++) {
 508           for (Int_t j = 0; j <= h->GetNbinsY()+1; j++) {
 509             hclone->SetBinContent(i,j,h->GetBinContent(i,j));
 510             hclone->SetBinError(i,j,h->GetBinError(i,j));
 511           }
 512         }
 513         // hclone->Add(h);
 514         h->Reset();
 515       }
 516       if (fDoTiming) copyTime += timer.CpuTime();
 517
 518       // --- Store Poisson result ------------------------------------
 519       if (fDoTiming) timer.Start(true);
 520       TH2D* poisson = rh->fPoisson.Result();
 521       for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
 522         for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
 523
 524           Double_t poissonV = poisson->GetBinContent(t+1,s+1);
 525           // Use cached eta - since the calls to GetEtaFromStrip and
 526           // GetPhiFromStrip are _very_ expensive
 527           Double_t  phi  = phiCache[s*nt+t];
 528           Double_t  eta  = etaCache[s*nt+t];
 529           // Double_t  phi  = fmd.Phi(d,r,s,t) * TMath::DegToRad();
 530           // Double_t  eta  = fmd.Eta(d,r,s,t);
 531           // if (fRecalculateEta)
 532           //   eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
 533           // if (fRecalculatePhi)
 534           //   phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
 535           if (fUsePoisson) {
 536             h->Fill(eta,phi,poissonV);
 537             rh->fDensity->Fill(eta, phi, poissonV);
 538           }
 539           else
 540             hclone->Fill(eta,phi,poissonV);
 541         }
 542       }
 543       if (fDoTiming) poissonTime += timer.CpuTime();
 544
 545       // --- Make diagnostics - eloss vs poisson ---------------------
 546       if (fDoTiming) timer.Start(true);
 547       Int_t nY = h->GetNbinsY();
 548       Int_t nIn  = 0; // Count non-outliers
 549       Int_t nOut = 0; // Count outliers
 550       for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
 551         // Set the overflow bin to contain the phi acceptance
 552         Double_t phiAcc  = rh->fGood->GetBinContent(ieta);
 553         Double_t phiAccE = rh->fGood->GetBinError(ieta);
 554         h->SetBinContent(ieta, nY+1, phiAcc);
 555         h->SetBinError(ieta, nY+1, phiAccE);
 556         Double_t eta     = h->GetXaxis()->GetBinCenter(ieta);
 557         rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
 558         for (Int_t iphi=1; iphi<= nY; iphi++) {
 559
 560           Double_t poissonV =  0; //h->GetBinContent(,s+1);
 561           Double_t eLossV =  0;
 562           if(fUsePoisson) {
 563             poissonV = h->GetBinContent(ieta,iphi);
 564             eLossV  = hclone->GetBinContent(ieta,iphi);
 565           }
 566           else {
 567             poissonV = hclone->GetBinContent(ieta,iphi);
 568             eLossV  = h->GetBinContent(ieta,iphi);
 569           }
 570
 571           if (poissonV < 1e-12 && eLossV < 1e-12)
 572             // we do not care about trivially empty bins
 573             continue;
 574
 575           Bool_t   outlier = CheckOutlier(eLossV, poissonV, fOutlierCut);
 576           Double_t rel     = eLossV < 1e-12 ? 0 : (poissonV - eLossV) / eLossV;
 577           if (outlier) {
 578             rh->fELossVsPoissonOut->Fill(eLossV, poissonV);
 579             rh->fDiffELossPoissonOut->Fill(rel);
 580             nOut++;
 581         }
 582           else {
 583             rh->fELossVsPoisson->Fill(eLossV, poissonV);
 584             rh->fDiffELossPoisson->Fill(rel);
 585             nIn++;
 586           } // if (outlier)
 587         } // for (iphi)
 588       } // for (ieta)
 589       Double_t outRatio = Double_t(nOut) / (nIn+nOut);
 590       rh->fOutliers->Fill(outRatio);
 591       if (outRatio < fMaxOutliers) rh->fPoisson.FillDiagnostics();
 592       else                         h->SetBit(AliForwardUtil::kSkipRing);
 593       if (fDoTiming) diagTime += timer.CpuTime();
 594       // delete hclone;
 595
 596     } // for q
 597   } // for d
 598
 599   if (fDoTiming) {
 600     fHTiming->Fill(1,reEtaTime);
 601     fHTiming->Fill(2,nPartTime);
 602     fHTiming->Fill(3,corrTime);
 603     fHTiming->Fill(4,rePhiTime);
 604     fHTiming->Fill(5,copyTime);
 605     fHTiming->Fill(6,poissonTime);
 606     fHTiming->Fill(7,diagTime);
 607     fHTiming->Fill(8,totalT.CpuTime());
 608   }
 609
 610   return kTRUE;
 611 }
 612
 613 //_____________________________________________________________________
 614 Bool_t
 615 AliFMDDensityCalculator::CheckOutlier(Double_t eloss,
 616                                       Double_t poisson,
 617                                       Double_t cut) const
 618 {
 619   if (eloss < 1e-6) return true;
 620   Double_t diff = TMath::Abs(poisson - eloss) / eloss;
 621   return diff > cut;
 622 }
 623 //_____________________________________________________________________
 624 Int_t
 625 AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
 626                                        UShort_t d, Char_t r, Int_t iEta) const
 627 {
 628   //
 629   // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
 630   //
 631   // Parameters:
 632   //    cor   Correction
 633   //    d     Detector
 634   //    r     Ring
 635   //    iEta  Eta bin
 636   //
 637   DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
 638   if(!cor) return -1;
 639
 640   AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,iEta, -1);
 641   if (!fit) {
 642     // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
 643     return -1;
 644   }
 645   return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError,
 646                             AliFMDCorrELossFit::ELossFit::fgLeastWeight,
 647                             fMaxParticles);
 648 }
 649
 650 //_____________________________________________________________________
 651 void
 652 AliFMDDensityCalculator::CacheMaxWeights(const TAxis& axis)
 653 {
 654   //
 655   // Find the max weights and cache them
 656   //
 657   DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
 658   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 659   const AliFMDCorrELossFit*     cor = fcm.GetELossFit();
 660   cor->CacheBins(fMinQuality);
 661
 662   TAxis eta(axis.GetNbins(),
 663             axis.GetXmin(),
 664             axis.GetXmax());
 665
 666   if(cor)
 667     eta.Set(cor->GetEtaAxis().GetNbins(),
 668             cor->GetEtaAxis().GetXmin(),
 669             cor->GetEtaAxis().GetXmax());
 670
 671   Int_t nEta = eta.GetNbins();
 672   fFMD1iMax.Set(nEta);
 673   fFMD2iMax.Set(nEta);
 674   fFMD2oMax.Set(nEta);
 675   fFMD3iMax.Set(nEta);
 676   fFMD3oMax.Set(nEta);
 677
 678   fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
 679   fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
 680   fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
 681   fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
 682   fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
 683   fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
 684
 685   AliInfo(Form("Get eta axis with %d bins from %f to %f",
 686                nEta, eta.GetXmin(), eta.GetXmax()));
 687   fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
 688   fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
 689   fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
 690   fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
 691   fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
 692   fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
 693
 694   for (Int_t i = 0; i < nEta; i++) {
 695     Double_t w[5];
 696     w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
 697     w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
 698     w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
 699     w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
 700     w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
 701     Double_t l[5];
 702     l[0] = GetMultCut(1, 'I', i+1, false);
 703     l[1] = GetMultCut(2, 'I', i+1, false);
 704     l[2] = GetMultCut(2, 'O', i+1, false);
 705     l[3] = GetMultCut(3, 'I', i+1, false);
 706     l[4] = GetMultCut(3, 'O', i+1, false);
 707     for (Int_t j = 0; j < 5; j++) {
 708       if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
 709       if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
 710     }
 711   }
 712 }
 713
 714 //_____________________________________________________________________
 715 Int_t
 716 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
 717 {
 718   //
 719   // Find the (cached) maximum weight for FMD<i>dr</i> in
 720   // @f$\eta@f$ bin @a iEta
 721   //
 722   // Parameters:
 723   //    d     Detector
 724   //    r     Ring
 725   //    iEta  Eta bin
 726   //
 727   // Return:
 728   //    max weight or <= 0 in case of problems
 729   //
 730   if (iEta < 0) return -1;
 731
 732   const TArrayI* max  = 0;
 733   switch (d) {
 734   case 1:  max = &fFMD1iMax;                                       break;
 735   case 2:  max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
 736   case 3:  max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
 737   }
 738   if (!max) {
 739     AliWarning(Form("No array for FMD%d%c", d, r));
 740     return -1;
 741   }
 742
 743   if (iEta >= max->fN) {
 744     AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
 745                     iEta, max->fN-1));
 746     return -1;
 747   }
 748
 749   AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
 750                    max->At(iEta)));
 751   return max->At(iEta);
 752 }
 753
 754 //_____________________________________________________________________
 755 Int_t
 756 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
 757 {
 758   //
 759   // Find the (cached) maximum weight for FMD<i>dr</i> iat
 760   // @f$\eta@f$
 761   //
 762   // Parameters:
 763   //    d     Detector
 764   //    r     Ring
 765   //    eta   Eta bin
 766   //
 767   // Return:
 768   //    max weight or <= 0 in case of problems
 769   //
 770   AliForwardCorrectionManager&  fcm  = AliForwardCorrectionManager::Instance();
 771   Int_t                         iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
 772
 773   return GetMaxWeight(d, r, iEta);
 774 }
 775
 776 //_____________________________________________________________________
 777 Float_t
 778 AliFMDDensityCalculator::NParticles(Float_t  mult,
 779                                     UShort_t d,
 780                                     Char_t   r,
 781                                     Float_t  eta,
 782                                     Bool_t   lowFlux) const
 783 {
 784   //
 785   // Get the number of particles corresponding to the signal mult
 786   //
 787   // Parameters:
 788   //    mult     Signal
 789   //    d        Detector
 790   //    r        Ring
 791   //    s        Sector
 792   //    t        Strip (not used)
 793   //    v        Vertex bin
 794   //    eta      Pseudo-rapidity
 795   //    lowFlux  Low-flux flag
 796   //
 797   // Return:
 798   //    The number of particles
 799   //
 800   // if (mult <= GetMultCut()) return 0;
 801   DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
 802   if (lowFlux) return 1;
 803
 804   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 805   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta, -1);
 806   if (!fit) {
 807     AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f qual=%d",
 808                     d, r, eta, fMinQuality));
 809     return 0;
 810   }
 811
 812   Int_t    m   = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
 813   if (m < 1) {
 814     AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
 815     return 0;
 816   }
 817
 818   UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
 819   Double_t ret = fit->EvaluateWeighted(mult, n);
 820
 821   if (fDebug > 10) {
 822     AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
 823   }
 824
 825   fWeightedSum->Fill(ret);
 826   fSumOfWeights->Fill(ret);
 827
 828   return ret;
 829 }
 830
 831 //_____________________________________________________________________
 832 Float_t
 833 AliFMDDensityCalculator::Correction(UShort_t d,
 834                                     Char_t   r,
 835                                     UShort_t t,
 836                                     Float_t  eta,
 837                                     Bool_t   lowFlux) const
 838 {
 839   //
 840   // Get the inverse correction factor.  This consist of
 841   //
 842   // - acceptance correction (corners of sensors)
 843   // - double hit correction (for low-flux events)
 844   // - dead strip correction
 845   //
 846   // Parameters:
 847   //    d        Detector
 848   //    r        Ring
 849   //    s        Sector
 850   //    t        Strip (not used)
 851   //    v        Vertex bin
 852   //    eta      Pseudo-rapidity
 853   //    lowFlux  Low-flux flag
 854   //
 855   // Return:
 856   //
 857   //
 858   DGUARD(fDebug, 10, "Apply correction in FMD density calculator");
 859   Float_t correction = 1;
 860   if (fUsePhiAcceptance == kPhiCorrectNch)
 861     correction = AcceptanceCorrection(r,t);
 862   if (lowFlux) {
 863     AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 864
 865     TH1D* dblHitCor = 0;
 866     if (fcm.GetDoubleHit())
 867       dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
 868
 869     if (dblHitCor) {
 870       Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
 871       if (dblC > 0) correction *= dblC;
 872     }
 873     // else {
 874     //   AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
 875     // }
 876   }
 877   return correction;
 878 }
 879
 880 //_____________________________________________________________________
 881 TH1D*
 882 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
 883 {
 884   //
 885   // Generate the acceptance corrections
 886   //
 887   // Parameters:
 888   //    r Ring to generate for
 889   //
 890   // Return:
 891   //    Newly allocated histogram of acceptance corrections
 892   //
 893   DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
 894   const Double_t ic1[] = { 4.9895, 15.3560 };
 895   const Double_t ic2[] = { 1.8007, 17.2000 };
 896   const Double_t oc1[] = { 4.2231, 26.6638 };
 897   const Double_t oc2[] = { 1.8357, 27.9500 };
 898   const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
 899   const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
 900   Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
 901   Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
 902   Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
 903   Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
 904   Float_t   basearc    = 2 * TMath::Pi() / nSec;
 905   Double_t  rad        = maxR - minR;
 906   Float_t   segment    = rad / nStrips;
 907   Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
 908
 909   // Numbers used to find end-point of strip.
 910   // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
 911   Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
 912   Float_t dx           = c2[0] - c1[0];
 913   Float_t dy           = c2[1] - c1[1];
 914   Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);
 915
 916   TH1D* ret = new TH1D(Form("acc%c", r),
 917                        Form("Acceptance correction for FMDx%c", r),
 918                        nStrips, -.5, nStrips-.5);
 919   ret->SetXTitle("Strip");
 920   ret->SetYTitle("#varphi acceptance");
 921   ret->SetDirectory(0);
 922   ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
 923   ret->SetFillStyle(3001);
 924
 925   for (Int_t t = 0; t < nStrips; t++) {
 926     Float_t   radius     = minR + t * segment;
 927
 928     // If the radius of the strip is smaller than the radius corresponding
 929     // to the first corner we have a full strip length
 930     if (radius <= cr) {
 931       ret->SetBinContent(t+1, 1);
 932       continue;
 933     }
 934
 935     // Next, we should find the end-point of the strip - that is,
 936     // the coordinates where the line from c1 to c2 intersects a circle
 937     // with radius given by the strip.
 938     // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
 939     // Calculate the determinant
 940     Float_t det = radius * radius * dr * dr - D*D;
 941
 942     if (det <=  0) {
 943       // <0 means No intersection
 944       // =0 means Exactly tangent
 945       ret->SetBinContent(t+1, 1);
 946       continue;
 947     }
 948
 949     // Calculate end-point and the corresponding opening angle
 950     Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
 951     Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
 952     Float_t th  = TMath::ATan2(x, y);
 953
 954     ret->SetBinContent(t+1, th / basearc);
 955   }
 956   return ret;
 957 }
 958
 959 //_____________________________________________________________________
 960 Float_t
 961 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
 962 {
 963   //
 964   // Get the acceptance correction for strip @a t in an ring of type @a r
 965   //
 966   // Parameters:
 967   //    r  Ring type ('I' or 'O')
 968   //    t  Strip number
 969   //
 970   // Return:
 971   //    Inverse acceptance correction
 972   //
 973   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
 974   return acc->GetBinContent(t+1);
 975 }
 976
 977 //____________________________________________________________________
 978 void
 979 AliFMDDensityCalculator::Terminate(const TList* dir, TList* output,
 980                                    Int_t nEvents)
 981 {
 982   //
 983   // Scale the histograms to the total number of events
 984   //
 985   // Parameters:
 986   //    dir     where to put the output
 987   //    nEvents Number of events
 988   //
 989   DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
 990   if (nEvents <= 0) return;
 991   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 992   if (!d) return;
 993
 994   TList* out = new TList;
 995   out->SetName(d->GetName());
 996   out->SetOwner();
 997
 998   TIter    next(&fRingHistos);
 999   RingHistos* o = 0;
1000   THStack* sums = new THStack("sums", "sums of ring signals");
1001   while ((o = static_cast<RingHistos*>(next()))) {
1002     o->Terminate(d, nEvents);
1003     if (!o->fDensity) {
1004       Warning("Terminate", "No density in %s", o->GetName());
1005       continue;
1006     }
1007     TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
1008                                          o->fDensity->GetNbinsY(),"e");
1009     sum->Scale(1., "width");
1010     sum->SetTitle(o->GetName());
1011     sum->SetDirectory(0);
1012     sum->SetYTitle("#sum N_{ch,incl}");
1013     sums->Add(sum);
1014   }
1015   out->Add(sums);
1016   output->Add(out);
1017 }
1018
1019 //____________________________________________________________________
1020 void
1021 AliFMDDensityCalculator::CreateOutputObjects(TList* dir)
1022 {
1023   //
1024   // Output diagnostic histograms to directory
1025   //
1026   // Parameters:
1027   //    dir List to write in
1028   //
1029   DGUARD(fDebug, 1, "Define output FMD density calculator");
1030   TList* d = new TList;
1031   d->SetOwner();
1032   d->SetName(GetName());
1033   dir->Add(d);
1034   d->Add(fWeightedSum);
1035   d->Add(fSumOfWeights);
1036   d->Add(fCorrections);
1037   d->Add(fAccI);
1038   d->Add(fAccO);
1039   d->Add(fMaxWeights);
1040   d->Add(fLowCuts);
1041
1042   TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
1043   nFiles->SetMergeMode('+');
1044
1045   // d->Add(sigma);
1046   d->Add(AliForwardUtil::MakeParameter("maxParticle",  fMaxParticles));
1047   d->Add(AliForwardUtil::MakeParameter("method",       fUsePoisson));
1048   d->Add(AliForwardUtil::MakeParameter("phiAcceptance",fUsePhiAcceptance));
1049   d->Add(AliForwardUtil::MakeParameter("etaLumping",   fEtaLumping));
1050   d->Add(AliForwardUtil::MakeParameter("phiLumping",   fPhiLumping));
1051   d->Add(AliForwardUtil::MakeParameter("recalcEta",    fRecalculateEta));
1052   d->Add(AliForwardUtil::MakeParameter("recalcPhi",    fRecalculatePhi));
1053   d->Add(AliForwardUtil::MakeParameter("maxOutliers",  fMaxOutliers));
1054   d->Add(AliForwardUtil::MakeParameter("outlierCut",   fOutlierCut));
1055   d->Add(nFiles);
1056   // d->Add(nxi);
1057   fCuts.Output(d,"lCuts");
1058
1059   TIter    next(&fRingHistos);
1060   RingHistos* o = 0;
1061   while ((o = static_cast<RingHistos*>(next()))) {
1062     o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
1063     o->CreateOutputObjects(d);
1064   }
1065
1066   if (!fDoTiming) return;
1067
1068   fHTiming = new TProfile("timing", "#LTt_{part}#GT", 8, .5, 8.5);
1069   fHTiming->SetDirectory(0);
1070   fHTiming->SetYTitle("#LTt_{part}#GT");
1071   fHTiming->SetXTitle("Part");
1072   fHTiming->SetFillColor(kRed+1);
1073   fHTiming->SetFillStyle(3001);
1074   fHTiming->SetMarkerStyle(20);
1075   fHTiming->SetMarkerColor(kBlack);
1076   fHTiming->SetLineColor(kBlack);
1077   fHTiming->SetStats(0);
1078   TAxis* xaxis = fHTiming->GetXaxis();
1079   xaxis->SetBinLabel(1, "Re-calculation of #eta");
1080   xaxis->SetBinLabel(2, "N_{particle}");
1081   xaxis->SetBinLabel(3, "Correction");
1082   xaxis->SetBinLabel(4, "Re-calculation of #phi");
1083   xaxis->SetBinLabel(5, "Copy to cache");
1084   xaxis->SetBinLabel(6, "Poisson calculation");
1085   xaxis->SetBinLabel(7, "Diagnostics");
1086   xaxis->SetBinLabel(8, "Total");
1087   d->Add(fHTiming);
1088 }
1089 #define PF(N,V,...)                                     \
1090   AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1091 #define PFB(N,FLAG)                             \
1092   do {                                                                  \
1093     AliForwardUtil::PrintName(N);                                       \
1094     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1095   } while(false)
1096 #define PFV(N,VALUE)                                    \
1097   do {                                                  \
1098     AliForwardUtil::PrintName(N);                       \
1099     std::cout << (VALUE) << std::endl; } while(false)
1100
1101 //____________________________________________________________________
1102 void
1103 AliFMDDensityCalculator::Print(Option_t* option) const
1104 {
1105   //
1106   // Print information
1107   //
1108   // Parameters:
1109   //    option Not used
1110   //
1111   AliForwardUtil::PrintTask(*this);
1112   gROOT->IncreaseDirLevel();
1113
1114   TString phiM("none");
1115   switch (fUsePhiAcceptance) {
1116   case kPhiNoCorrect:    phiM = "none"; break;
1117   case kPhiCorrectNch:   phiM = "correct Nch"; break;
1118   case kPhiCorrectELoss: phiM = "correct energy loss"; break;
1119   }
1120
1121   PFV("Max(particles)",         fMaxParticles );
1122   PFV("Poisson method",         fUsePoisson );
1123   PFV("Eta lumping",            fEtaLumping );
1124   PFV("Phi lumping",            fPhiLumping );
1125   PFV("Recalculate eta",        fRecalculateEta );
1126   PFV("Recalculate phi",        fRecalculatePhi );
1127   PFV("Use phi acceptance",     phiM);
1128   PFV("Lower cut", "");
1129   fCuts.Print();
1130
1131   TString opt(option);
1132   opt.ToLower();
1133   if (!opt.Contains("nomax")) {
1134     PFV("Max weights", "");
1135
1136     char ind[64];
1137     for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1138     ind[gROOT->GetDirLevel()] = '\0';
1139     for (UShort_t d=1; d<=3; d++) {
1140       UShort_t nr = (d == 1 ? 1 : 2);
1141       for (UShort_t q=0; q<nr; q++) {
1142         ind[gROOT->GetDirLevel()]   = ' ';
1143         ind[gROOT->GetDirLevel()+1] = '\0';
1144         Char_t      r = (q == 0 ? 'I' : 'O');
1145         std::cout << ind << " FMD" << d << r << ":";
1146         ind[gROOT->GetDirLevel()+1] = ' ';
1147         ind[gROOT->GetDirLevel()+2] = '\0';
1148
1149         const TArrayI& a = (d == 1 ? fFMD1iMax :
1150                             (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
1151                              (r == 'I' ? fFMD3iMax : fFMD3oMax)));
1152         Int_t j = 0;
1153         for (Int_t i = 0; i < a.fN; i++) {
1154           if (a.fArray[i] < 1) continue;
1155           if (j % 6 == 0)      std::cout << "\n " << ind;
1156           j++;
1157           std::cout << "  " << std::setw(3) << i << ": " << a.fArray[i];
1158         }
1159         std::cout << std::endl;
1160       }
1161     }
1162   }
1163   gROOT->DecreaseDirLevel();
1164 }
1165
1166 //====================================================================
1167 AliFMDDensityCalculator::RingHistos::RingHistos()
1168   : AliForwardUtil::RingHistos(),
1169     fList(0),
1170     // fEvsN(0),
1171     // fEvsM(0),
1172     // fEtaVsN(0),
1173     // fEtaVsM(0),
1174     fCorr(0),
1175     fDensity(0),
1176     fELossVsPoisson(0),
1177     fDiffELossPoisson(0),
1178     fELossVsPoissonOut(0),
1179     fDiffELossPoissonOut(0),
1180     fOutliers(0),
1181     fPoisson(),
1182     fELoss(0),
1183     fELossUsed(0),
1184     fMultCut(0),
1185     fTotal(0),
1186     fGood(0),
1187     fPhiAcc(0),
1188     fPhiBefore(0),
1189     fPhiAfter(0)
1190 {
1191   //
1192   // Default CTOR
1193   //
1194 }
1195
1196 //____________________________________________________________________
1197 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
1198   : AliForwardUtil::RingHistos(d,r),
1199     fList(0),
1200     // fEvsN(0),
1201     // fEvsM(0),
1202     // fEtaVsN(0),
1203     // fEtaVsM(0),
1204     fCorr(0),
1205     fDensity(0),
1206     fELossVsPoisson(0),
1207     fDiffELossPoisson(0),
1208     fELossVsPoissonOut(0),
1209     fDiffELossPoissonOut(0),
1210     fOutliers(0),
1211     fPoisson("ignored"),
1212     fELoss(0),
1213     fELossUsed(0),
1214     fMultCut(0),
1215     fTotal(0),
1216     fGood(0),
1217     fPhiAcc(0),
1218     fPhiBefore(0),
1219     fPhiAfter(0)
1220 {
1221   //
1222   // Constructor
1223   //
1224   // Parameters:
1225   //    d detector
1226   //    r ring
1227   //
1228 #if 0
1229   fEvsN = new TH2D("elossVsNnocorr",
1230                    "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
1231                    250, -.5, 24.5, 251, -1.5, 24.5);
1232   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
1233   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
1234   fEvsN->Sumw2();
1235   fEvsN->SetDirectory(0);
1236
1237   fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
1238   fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
1239   fEvsM->SetDirectory(0);
1240
1241   fEtaVsN = new TProfile("etaVsNnocorr",
1242                          "Average inclusive N_{ch} vs #eta (uncorrected)",
1243                          200, -4, 6);
1244   fEtaVsN->SetXTitle("#eta");
1245   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
1246   fEtaVsN->SetDirectory(0);
1247   fEtaVsN->SetLineColor(Color());
1248   fEtaVsN->SetFillColor(Color());
1249
1250   fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
1251   fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
1252   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
1253   fEtaVsM->SetDirectory(0);
1254 #endif
1255
1256   fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
1257   fCorr->SetXTitle("#eta");
1258   fCorr->SetYTitle("#LT correction#GT");
1259   fCorr->SetDirectory(0);
1260   fCorr->SetLineColor(Color());
1261   fCorr->SetFillColor(Color());
1262
1263   fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
1264                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
1265                       0, 2*TMath::Pi());
1266   fDensity->SetDirectory(0);
1267   fDensity->Sumw2();  fDensity->SetMarkerColor(Color());
1268   fDensity->SetXTitle("#eta");
1269   fDensity->SetYTitle("#phi [radians]");
1270   fDensity->SetZTitle("Inclusive N_{ch} density");
1271
1272   // --- Create increasing sized bins --------------------------------
1273   TArrayD bins;
1274   // bins, lowest order, higest order, return array
1275   const char* nchP = "N_{ch}^{Poisson}";
1276   const char* nchE = "N_{ch}^{#Delta}";
1277   AliForwardUtil::MakeLogScale(300, 0, 2, bins);
1278   fELossVsPoisson = new TH2D("elossVsPoisson",
1279                              "N_{ch} from energy loss vs from Poisson",
1280                              bins.GetSize()-1, bins.GetArray(),
1281                              bins.GetSize()-1, bins.GetArray());
1282   fELossVsPoisson->SetDirectory(0);
1283   fELossVsPoisson->SetXTitle(nchE);
1284   fELossVsPoisson->SetYTitle(nchP);
1285   fELossVsPoisson->SetZTitle("Correlation");
1286   fELossVsPoissonOut =
1287     static_cast<TH2D*>(fELossVsPoisson
1288                        ->Clone(Form("%sOutlier",
1289                                     fELossVsPoisson->GetName())));
1290   fELossVsPoissonOut->SetDirectory(0);
1291   fELossVsPoissonOut->SetMarkerStyle(20);
1292   fELossVsPoissonOut->SetMarkerSize(0.3);
1293   fELossVsPoissonOut->SetMarkerColor(kBlack);
1294   fELossVsPoissonOut->SetTitle(Form("%s for outliers",
1295                                     fELossVsPoisson->GetTitle()));
1296
1297   fDiffELossPoisson = new TH1D("diffElossPoisson",
1298                                Form("(%s-%s)/%s", nchP, nchE, nchE),
1299                                100, -1, 1);
1300   fDiffELossPoisson->SetDirectory(0);
1301   fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
1302   fDiffELossPoisson->SetYTitle("Frequency");
1303   fDiffELossPoisson->SetMarkerColor(Color());
1304   fDiffELossPoisson->SetFillColor(Color());
1305   fDiffELossPoisson->SetFillStyle(3001);
1306   fDiffELossPoisson->Sumw2();
1307
1308   fDiffELossPoissonOut =
1309     static_cast<TH1D*>(fDiffELossPoisson
1310                        ->Clone(Form("%sOutlier",fDiffELossPoisson->GetName())));
1311   fDiffELossPoissonOut->SetDirectory(0);
1312   fDiffELossPoissonOut->SetTitle(Form("%s for outliers",
1313                                       fDiffELossPoisson->GetTitle()));
1314   fDiffELossPoissonOut->SetMarkerColor(Color()-2);
1315   fDiffELossPoissonOut->SetFillColor(Color()-2);
1316   fDiffELossPoissonOut->SetFillStyle(3002);
1317
1318   fOutliers = new TH1D("outliers", "Fraction of outliers", 100, 0, 1);
1319   fOutliers->SetDirectory(0);
1320   fOutliers->SetXTitle("N_{outlier}/(N_{outlier}+N_{inside})");
1321   fOutliers->SetYTitle("#sum_{events}#sum_{bins}");
1322   fOutliers->SetFillColor(Color());
1323   fOutliers->SetFillStyle(3001);
1324   fOutliers->SetLineColor(kBlack);
1325
1326   fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
1327                     640, -1, 15);
1328   fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1329   fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1330   fELoss->SetFillColor(Color()-2);
1331   fELoss->SetFillStyle(3003);
1332   fELoss->SetLineColor(kBlack);
1333   fELoss->SetLineStyle(2);
1334   fELoss->SetLineWidth(1);
1335   fELoss->SetDirectory(0);
1336
1337   fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1338   fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1339   fELossUsed->SetFillStyle(3002);
1340   fELossUsed->SetLineStyle(1);
1341   fELossUsed->SetDirectory(0);
1342
1343   fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
1344                         (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1345   fPhiBefore->SetDirectory(0);
1346   fPhiBefore->SetXTitle("#phi");
1347   fPhiBefore->SetYTitle("Events");
1348   fPhiBefore->SetMarkerColor(Color());
1349   fPhiBefore->SetLineColor(Color());
1350   fPhiBefore->SetFillColor(Color());
1351   fPhiBefore->SetFillStyle(3001);
1352   fPhiBefore->SetMarkerStyle(20);
1353
1354   fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
1355   fPhiAfter->SetTitle("#phi distribution (after re-calc)");
1356   fPhiAfter->SetDirectory(0);
1357
1358 }
1359 //____________________________________________________________________
1360 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
1361   : AliForwardUtil::RingHistos(o),
1362     fList(o.fList),
1363     // fEvsN(o.fEvsN),
1364     // fEvsM(o.fEvsM),
1365     // fEtaVsN(o.fEtaVsN),
1366     // fEtaVsM(o.fEtaVsM),
1367     fCorr(o.fCorr),
1368     fDensity(o.fDensity),
1369     fELossVsPoisson(o.fELossVsPoisson),
1370     fDiffELossPoisson(o.fDiffELossPoisson),
1371     fELossVsPoissonOut(o.fELossVsPoissonOut),
1372     fDiffELossPoissonOut(o.fDiffELossPoissonOut),
1373     fOutliers(o.fOutliers),
1374     fPoisson(o.fPoisson),
1375     fELoss(o.fELoss),
1376     fELossUsed(o.fELossUsed),
1377     fMultCut(o.fMultCut),
1378     fTotal(o.fTotal),
1379     fGood(o.fGood),
1380     fPhiAcc(o.fPhiAcc),
1381     fPhiBefore(o.fPhiBefore),
1382     fPhiAfter(o.fPhiAfter)
1383 {
1384   //
1385   // Copy constructor
1386   //
1387   // Parameters:
1388   //    o Object to copy from
1389   //
1390 }
1391
1392 //____________________________________________________________________
1393 AliFMDDensityCalculator::RingHistos&
1394 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1395 {
1396   //
1397   // Assignment operator
1398   //
1399   // Parameters:
1400   //    o Object to assign from
1401   //
1402   // Return:
1403   //    Reference to this
1404   //
1405   if (&o == this) return *this;
1406   AliForwardUtil::RingHistos::operator=(o);
1407
1408   // if (fEvsN)             delete fEvsN;
1409   // if (fEvsM)             delete fEvsM;
1410   // if (fEtaVsN)           delete fEtaVsN;
1411   // if (fEtaVsM)           delete fEtaVsM;
1412   if (fCorr)             delete fCorr;
1413   if (fDensity)          delete fDensity;
1414   if (fELossVsPoisson)   delete fELossVsPoisson;
1415   if (fDiffELossPoisson) delete fDiffELossPoisson;
1416   if (fTotal)            delete fTotal;
1417   if (fGood)             delete fGood;
1418   if (fPhiAcc)           delete fPhiAcc;
1419   if (fPhiBefore)        delete fPhiBefore;
1420   if (fPhiAfter)         delete fPhiAfter;
1421
1422   // fEvsN             = static_cast<TH2D*>(o.fEvsN->Clone());
1423   // fEvsM             = static_cast<TH2D*>(o.fEvsM->Clone());
1424   // fEtaVsN           = static_cast<TProfile*>(o.fEtaVsN->Clone());
1425   // fEtaVsM           = static_cast<TProfile*>(o.fEtaVsM->Clone());
1426   fCorr             = static_cast<TProfile*>(o.fCorr->Clone());
1427   fDensity          = static_cast<TH2D*>(o.fDensity->Clone());
1428   fELossVsPoisson   = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1429   fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1430   fELossVsPoissonOut   = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1431   fDiffELossPoissonOut = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1432   fOutliers            = static_cast<TH1D*>(o.fOutliers->Clone());
1433   fPoisson          = o.fPoisson;
1434   fELoss            = static_cast<TH1D*>(o.fELoss->Clone());
1435   fELossUsed        = static_cast<TH1D*>(o.fELossUsed->Clone());
1436   fTotal            = static_cast<TH1D*>(o.fTotal->Clone());
1437   fGood             = static_cast<TH1D*>(o.fGood->Clone());
1438   fPhiAcc           = static_cast<TH2D*>(o.fPhiAcc->Clone());
1439   fPhiBefore        = static_cast<TH1D*>(o.fPhiBefore->Clone());
1440   fPhiAfter         = static_cast<TH1D*>(o.fPhiAfter->Clone());
1441   return *this;
1442 }
1443 //____________________________________________________________________
1444 AliFMDDensityCalculator::RingHistos::~RingHistos()
1445 {
1446   //
1447   // Destructor
1448   //
1449 }
1450
1451
1452 //____________________________________________________________________
1453 void
1454 AliFMDDensityCalculator::RingHistos::SetupForData(const TAxis& eAxis)
1455 {
1456   // Initialize
1457   // This is called on first event
1458   fPoisson.Init(-1,-1);
1459   fTotal = new TH1D("total", "Total # of strips per #eta",
1460                     eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
1461   fTotal->SetDirectory(0);
1462   fTotal->SetXTitle("#eta");
1463   fTotal->SetYTitle("# of strips");
1464   fGood = static_cast<TH1D*>(fTotal->Clone("good"));
1465   fGood->SetTitle("# of good strips per #eta");
1466   fGood->SetDirectory(0);
1467
1468   fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}",
1469                      eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
1470                      10, -10, 10);
1471   fPhiAcc->SetXTitle("#eta");
1472   fPhiAcc->SetYTitle("v_{z} [cm]");
1473   fPhiAcc->SetZTitle("#phi acceptance");
1474   fPhiAcc->SetDirectory(0);
1475
1476   if (fList) fList->Add(fPhiAcc);
1477 }
1478
1479 //____________________________________________________________________
1480 void
1481 AliFMDDensityCalculator::RingHistos::CreateOutputObjects(TList* dir)
1482 {
1483   //
1484   // Make output.  This is called as part of SlaveBegin
1485   //
1486   // Parameters:
1487   //    dir Where to put it
1488   //
1489   TList* d = DefineOutputList(dir);
1490   // d->Add(fEvsN);
1491   // d->Add(fEvsM);
1492   // d->Add(fEtaVsN);
1493   // d->Add(fEtaVsM);
1494   d->Add(fCorr);
1495   d->Add(fDensity);
1496   d->Add(fELossVsPoisson);
1497   d->Add(fELossVsPoissonOut);
1498   d->Add(fDiffELossPoisson);
1499   d->Add(fDiffELossPoissonOut);
1500   d->Add(fOutliers);
1501   fPoisson.Output(d);
1502   fPoisson.GetOccupancy()->SetFillColor(Color());
1503   fPoisson.GetMean()->SetFillColor(Color());
1504   fPoisson.GetOccupancy()->SetFillColor(Color());
1505   d->Add(fELoss);
1506   d->Add(fELossUsed);
1507   d->Add(fPhiBefore);
1508   d->Add(fPhiAfter);
1509
1510   TAxis x(NStrip(), -.5, NStrip()-.5);
1511   TAxis y(NSector(), -.5, NSector()-.5);
1512   x.SetTitle("strip");
1513   y.SetTitle("sector");
1514   fPoisson.Define(x, y);
1515
1516   d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
1517   fList = d;
1518 }
1519
1520 //____________________________________________________________________
1521 void
1522 AliFMDDensityCalculator::RingHistos::Terminate(TList* dir, Int_t nEvents)
1523 {
1524   //
1525   // Scale the histograms to the total number of events
1526   //
1527   // Parameters:
1528   //    dir     Where the output is
1529   //    nEvents Number of events
1530   //
1531   TList* l = GetOutputList(dir);
1532   if (!l) return;
1533
1534   TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"inclDensity"));
1535   if (density) density->Scale(1./nEvents);
1536   fDensity = density;
1537 }
1538
1539 //____________________________________________________________________
1540 //
1541 // EOF
1542 //
1543
1544
1545