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