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