PWG2/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 <iostream>
  18 #include <iomanip>
  19
  20 ClassImp(AliFMDDensityCalculator)
  21 #if 0
  22 ; // For Emacs
  23 #endif
  24
  25 //____________________________________________________________________
  26 AliFMDDensityCalculator::AliFMDDensityCalculator()
  27   : TNamed(),
  28     fRingHistos(),
  29     fMultCut(0),
  30     fNXi(1),
  31     fIncludeSigma(true),
  32     fSumOfWeights(0),
  33     fWeightedSum(0),
  34     fCorrections(0),
  35     fMaxParticles(5),
  36     fUsePoisson(false),
  37     fUsePhiAcceptance(false),
  38     fAccI(0),
  39     fAccO(0),
  40     fFMD1iMax(0),
  41     fFMD2iMax(0),
  42     fFMD2oMax(0),
  43     fFMD3iMax(0),
  44     fFMD3oMax(0),
  45     fMaxWeights(0),
  46     fLowCuts(0),
  47     fEtaLumping(5),
  48     fPhiLumping(5),
  49     fDebug(0)
  50 {
  51   //
  52   // Constructor
  53   //
  54
  55 }
  56
  57 //____________________________________________________________________
  58 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
  59   : TNamed("fmdDensityCalculator", title),
  60     fRingHistos(),
  61     fMultCut(0),
  62     fNXi(1),
  63     fIncludeSigma(true),
  64     fSumOfWeights(0),
  65     fWeightedSum(0),
  66     fCorrections(0),
  67     fMaxParticles(5),
  68     fUsePoisson(false),
  69     fUsePhiAcceptance(false),
  70     fAccI(0),
  71     fAccO(0),
  72     fFMD1iMax(0),
  73     fFMD2iMax(0),
  74     fFMD2oMax(0),
  75     fFMD3iMax(0),
  76     fFMD3oMax(0),
  77     fMaxWeights(0),
  78     fLowCuts(0),
  79     fEtaLumping(5),
  80     fPhiLumping(5),
  81     fDebug(0)
  82 {
  83   //
  84   // Constructor
  85   //
  86   // Parameters:
  87   //    name Name of object
  88   //
  89   fRingHistos.SetName(GetName());
  90   fRingHistos.SetOwner();
  91   fRingHistos.Add(new RingHistos(1, 'I'));
  92   fRingHistos.Add(new RingHistos(2, 'I'));
  93   fRingHistos.Add(new RingHistos(2, 'O'));
  94   fRingHistos.Add(new RingHistos(3, 'I'));
  95   fRingHistos.Add(new RingHistos(3, 'O'));
  96   fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
  97                            200, 0, 20);
  98   fSumOfWeights->SetFillColor(kRed+1);
  99   fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
 100   fWeightedSum  = new TH1D("weightedSum", "Weighted sum of Landau propability",
 101                            200, 0, 20);
 102   fWeightedSum->SetFillColor(kBlue+1);
 103   fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
 104   fCorrections  = new TH1D("corrections", "Distribution of corrections",
 105                            100, 0, 10);
 106   fCorrections->SetFillColor(kBlue+1);
 107   fCorrections->SetXTitle("correction");
 108
 109   fAccI = GenerateAcceptanceCorrection('I');
 110   fAccO = GenerateAcceptanceCorrection('O');
 111
 112   fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
 113                          1, 0, 1, 1, 0, 1);
 114   fMaxWeights->SetXTitle("#eta");
 115   fMaxWeights->SetDirectory(0);
 116
 117   fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
 118   fLowCuts->SetXTitle("#eta");
 119   fLowCuts->SetDirectory(0);
 120
 121   for (Int_t i = 0; i < 5; i++) fMultCuts[i] = 0;
 122 }
 123
 124 //____________________________________________________________________
 125 AliFMDDensityCalculator::AliFMDDensityCalculator(const
 126                                                  AliFMDDensityCalculator& o)
 127   : TNamed(o),
 128     fRingHistos(),
 129     fMultCut(o.fMultCut),
 130     fNXi(o.fNXi),
 131     fIncludeSigma(o.fIncludeSigma),
 132     fSumOfWeights(o.fSumOfWeights),
 133     fWeightedSum(o.fWeightedSum),
 134     fCorrections(o.fCorrections),
 135     fMaxParticles(o.fMaxParticles),
 136     fUsePoisson(o.fUsePoisson),
 137     fUsePhiAcceptance(o.fUsePhiAcceptance),
 138     fAccI(o.fAccI),
 139     fAccO(o.fAccO),
 140     fFMD1iMax(o.fFMD1iMax),
 141     fFMD2iMax(o.fFMD2iMax),
 142     fFMD2oMax(o.fFMD2oMax),
 143     fFMD3iMax(o.fFMD3iMax),
 144     fFMD3oMax(o.fFMD3oMax),
 145     fMaxWeights(o.fMaxWeights),
 146     fLowCuts(o.fLowCuts),
 147     fEtaLumping(o.fEtaLumping),
 148     fPhiLumping(o.fPhiLumping),
 149     fDebug(o.fDebug)
 150 {
 151   //
 152   // Copy constructor
 153   //
 154   // Parameters:
 155   //    o Object to copy from
 156   //
 157   TIter    next(&o.fRingHistos);
 158   TObject* obj = 0;
 159   while ((obj = next())) fRingHistos.Add(obj);
 160 }
 161
 162 //____________________________________________________________________
 163 AliFMDDensityCalculator::~AliFMDDensityCalculator()
 164 {
 165   //
 166   // Destructor
 167   //
 168   fRingHistos.Delete();
 169 }
 170
 171 //____________________________________________________________________
 172 AliFMDDensityCalculator&
 173 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
 174 {
 175   //
 176   // Assignement operator
 177   //
 178   // Parameters:
 179   //    o Object to assign from
 180   //
 181   // Return:
 182   //    Reference to this object
 183   //
 184   TNamed::operator=(o);
 185
 186   fMultCut      = o.fMultCut;
 187   fNXi          = o.fNXi;
 188   fIncludeSigma = o.fIncludeSigma;
 189   fDebug        = o.fDebug;
 190   fMaxParticles = o.fMaxParticles;
 191   fUsePoisson   = o.fUsePoisson;
 192   fUsePhiAcceptance  = o.fUsePhiAcceptance;
 193   fAccI         = o.fAccI;
 194   fAccO         = o.fAccO;
 195   fFMD1iMax     = o.fFMD1iMax;
 196   fFMD2iMax     = o.fFMD2iMax;
 197   fFMD2oMax     = o.fFMD2oMax;
 198   fFMD3iMax     = o.fFMD3iMax;
 199   fFMD3oMax     = o.fFMD3oMax;
 200   fMaxWeights   = o.fMaxWeights;
 201   fLowCuts      = o.fLowCuts;
 202   fEtaLumping   = o.fEtaLumping;
 203   fPhiLumping   = o.fPhiLumping;
 204   fRingHistos.Delete();
 205   TIter    next(&o.fRingHistos);
 206   TObject* obj = 0;
 207   while ((obj = next())) fRingHistos.Add(obj);
 208
 209   return *this;
 210 }
 211
 212 //____________________________________________________________________
 213 void
 214 AliFMDDensityCalculator::Init(const TAxis& axis)
 215 {
 216   // Intialize this sub-algorithm
 217   //
 218   // Parameters:
 219   //   etaAxis   Not used
 220   CacheMaxWeights();
 221
 222   TIter    next(&fRingHistos);
 223   RingHistos* o = 0;
 224   while ((o = static_cast<RingHistos*>(next())))
 225     o->Init(axis);
 226 }
 227
 228 //____________________________________________________________________
 229 AliFMDDensityCalculator::RingHistos*
 230 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
 231 {
 232   //
 233   // Get the ring histogram container
 234   //
 235   // Parameters:
 236   //    d Detector
 237   //    r Ring
 238   //
 239   // Return:
 240   //    Ring histogram container
 241   //
 242   Int_t idx = -1;
 243   switch (d) {
 244   case 1: idx = 0; break;
 245   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 246   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 247   }
 248   if (idx < 0 || idx >= fRingHistos.GetEntries()) {
 249     AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
 250     return 0;
 251   }
 252
 253   return static_cast<RingHistos*>(fRingHistos.At(idx));
 254 }
 255
 256 //____________________________________________________________________
 257 void
 258 AliFMDDensityCalculator::SetMultCuts(Double_t fmd1i,
 259                                      Double_t fmd2i,
 260                                      Double_t fmd2o,
 261                                      Double_t fmd3i,
 262                                      Double_t fmd3o)
 263 {
 264   fMultCuts[0] = fmd1i;
 265   fMultCuts[1] = fmd2i;
 266   fMultCuts[2] = fmd2o;
 267   fMultCuts[3] = fmd3i;
 268   fMultCuts[4] = fmd3o;
 269   fMultCut = (fmd1i+fmd2i+fmd2o+fmd3i+fmd3o) / 5;
 270 }
 271
 272 //____________________________________________________________________
 273 Double_t
 274 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
 275                                     Bool_t errors) const
 276 {
 277   //
 278   // Get the multiplicity cut.  If the user has set fMultCut (via
 279   // SetMultCut) then that value is used.  If not, then the lower
 280   // value of the fit range for the enery loss fits is returned.
 281   //
 282   // Return:
 283   //    Lower cut on multiplicity
 284   //
 285   Int_t    idx = (d == 1 ? 0 : 2*(d - 2) + 1 + ((r=='I' || r=='i') ? 0 : 1));
 286   if (fMultCuts[idx] > 0) return fMultCuts[idx];
 287   if (fMultCut > 0) return fMultCut;
 288
 289   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 290   AliFMDCorrELossFit* fits = fcm.GetELossFit();
 291   if (fNXi < 0) return fits->GetLowCut();
 292
 293   return fits->GetLowerBound(d, r, ieta, fNXi, errors, fIncludeSigma);
 294 }
 295
 296 //____________________________________________________________________
 297 Double_t
 298 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
 299                                     Bool_t errors) const
 300 {
 301   //
 302   // Get the multiplicity cut.  If the user has set fMultCut (via
 303   // SetMultCut) then that value is used.  If not, then the lower
 304   // value of the fit range for the enery loss fits is returned.
 305   //
 306   // Return:
 307   //    Lower cut on multiplicity
 308   //
 309   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 310   AliFMDCorrELossFit* fits = fcm.GetELossFit();
 311   Int_t iEta = fits->FindEtaBin(eta);
 312
 313   return GetMultCut(d, r, iEta, errors);
 314 }
 315
 316 //____________________________________________________________________
 317 Bool_t
 318 AliFMDDensityCalculator::Calculate(const AliESDFMD&        fmd,
 319                                    AliForwardUtil::Histos& hists,
 320                                    UShort_t                vtxbin,
 321                                    Bool_t                  lowFlux)
 322 {
 323   //
 324   // Do the calculations
 325   //
 326   // Parameters:
 327   //    fmd      AliESDFMD object (possibly) corrected for sharing
 328   //    hists    Histogram cache
 329   //    vtxBin   Vertex bin
 330   //    lowFlux  Low flux flag.
 331   //
 332   // Return:
 333   //    true on successs
 334   //
 335   for (UShort_t d=1; d<=3; d++) {
 336     UShort_t nr = (d == 1 ? 1 : 2);
 337     for (UShort_t q=0; q<nr; q++) {
 338       Char_t      r = (q == 0 ? 'I' : 'O');
 339       UShort_t    ns= (q == 0 ?  20 :  40);
 340       UShort_t    nt= (q == 0 ? 512 : 256);
 341       TH2D*       h = hists.Get(d,r);
 342       RingHistos* rh= GetRingHistos(d,r);
 343       if (!rh) {
 344         AliError(Form("No ring histogram found for FMD%d%c", d, r));
 345         fRingHistos.ls();
 346         return false;
 347       }
 348       rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
 349
 350       for (UShort_t s=0; s<ns; s++) {
 351         for (UShort_t t=0; t<nt; t++) {
 352           Float_t  mult = fmd.Multiplicity(d,r,s,t);
 353           Float_t  phi  = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
 354           Float_t  eta  = fmd.Eta(d,r,s,t);
 355           rh->fTotalStrips->Fill(eta, phi);
 356
 357           if (mult == AliESDFMD::kInvalidMult || mult > 20) {
 358             rh->fEmptyStrips->Fill(eta,phi);
 359             rh->fEvsM->Fill(mult,0);
 360             continue;
 361           }
 362
 363           Double_t cut  = 1024;
 364           if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
 365
 366           Double_t n   = 0;
 367           if (cut > 0 && mult > cut)
 368             n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
 369
 370           rh->fEvsN->Fill(mult,n);
 371           rh->fEtaVsN->Fill(eta, n);
 372
 373           Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
 374           fCorrections->Fill(c);
 375           if (c > 0) n /= c;
 376           rh->fEvsM->Fill(mult,n);
 377           rh->fEtaVsM->Fill(eta, n);
 378           rh->fCorr->Fill(eta, c);
 379
 380           if (n > 0.9 && c > 0) rh->fBasicHits->Fill(eta,phi, 1./c);
 381           else                  rh->fEmptyStrips->Fill(eta,phi);
 382
 383           h->Fill(eta,phi,n);
 384           if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
 385         } // for t
 386       } // for s
 387
 388
 389       // --- Loop over poisson histograms ----------------------------
 390       for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
 391         for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
 392           Double_t eLossV   = h->GetBinContent(ieta, iphi);
 393           Double_t eta      = h->GetXaxis()->GetBinCenter(ieta);
 394           Double_t phi      = h->GetYaxis()->GetBinCenter(iphi);
 395           Int_t    jeta     = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
 396           Int_t    jphi     = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
 397           Double_t empty    = rh->fEmptyStrips->GetBinContent(jeta, jphi);
 398           Double_t total    = rh->fTotalStrips->GetBinContent(jeta, jphi);
 399           Double_t hits     = rh->fBasicHits->GetBinContent(ieta,iphi);
 400
 401           // Mean in region of interest
 402           Double_t poissonM = (total <= 0 || empty <= 0 ? 0 :
 403                                -TMath::Log(empty / total));
 404
 405           // Note, that given filled=total-empty, and
 406           //
 407           //     m = -log(empty/total)
 408           //       = -log(1 - filled/total)
 409           //
 410           //     v = m / (1 - exp(-m))
 411           //       = -total/filled * (log(total-filled)-log(total))
 412           //       = -total / (total-empty) * log(empty/total)
 413           //       = total (log(total)-log(empty)) / (total-empty)
 414           //
 415           Double_t poissonV = hits;
 416           if(poissonM > 0)
 417             // Correct for counting statistics and weight by counts
 418             poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
 419           Double_t poissonE = TMath::Sqrt(hits);
 420           if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
 421
 422           rh->fELossVsPoisson->Fill(eLossV, poissonV);
 423           rh->fEmptyVsTotal->Fill(total, empty);
 424           if (fUsePoisson) {
 425             h->SetBinContent(ieta,iphi,poissonV);
 426             h->SetBinError(ieta,iphi,poissonE);
 427             rh->fDensity->Fill(eta, phi, poissonV);
 428           }
 429         }
 430       }
 431
 432     } // for q
 433   } // for d
 434
 435   return kTRUE;
 436 }
 437
 438 //_____________________________________________________________________
 439 Int_t
 440 AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
 441                                        UShort_t d, Char_t r, Int_t iEta) const
 442 {
 443   //
 444   // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
 445   //
 446   // Parameters:
 447   //    cor   Correction
 448   //    d     Detector
 449   //    r     Ring
 450   //    iEta  Eta bin
 451   //
 452   AliFMDCorrELossFit::ELossFit* fit = cor->GetFit(d,r,iEta);
 453   if (!fit) {
 454     // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
 455     return -1;
 456   }
 457   return TMath::Min(Int_t(fMaxParticles), fit->FindMaxWeight());
 458 }
 459
 460 //_____________________________________________________________________
 461 void
 462 AliFMDDensityCalculator::CacheMaxWeights()
 463 {
 464   //
 465   // Find the max weights and cache them
 466   //
 467   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 468   AliFMDCorrELossFit*           cor = fcm.GetELossFit();
 469   const TAxis&                  eta = cor->GetEtaAxis();
 470
 471   Int_t nEta = eta.GetNbins();
 472   fFMD1iMax.Set(nEta);
 473   fFMD2iMax.Set(nEta);
 474   fFMD2oMax.Set(nEta);
 475   fFMD3iMax.Set(nEta);
 476   fFMD3oMax.Set(nEta);
 477
 478   fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
 479   fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
 480   fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
 481   fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
 482   fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
 483   fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
 484
 485   AliInfo(Form("Get eta axis with %d bins from %f to %f",
 486                nEta, eta.GetXmin(), eta.GetXmax()));
 487   fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
 488   fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
 489   fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
 490   fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
 491   fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
 492   fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
 493
 494   for (Int_t i = 0; i < nEta; i++) {
 495     Double_t w[5];
 496     w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
 497     w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
 498     w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
 499     w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
 500     w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
 501     Double_t l[5];
 502     l[0] = GetMultCut(1, 'I', i+1, false);
 503     l[1] = GetMultCut(2, 'I', i+1, false);
 504     l[2] = GetMultCut(2, 'O', i+1, false);
 505     l[3] = GetMultCut(3, 'I', i+1, false);
 506     l[4] = GetMultCut(3, 'O', i+1, false);
 507     for (Int_t j = 0; j < 5; j++) {
 508       if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
 509       if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
 510     }
 511   }
 512 }
 513
 514 //_____________________________________________________________________
 515 Int_t
 516 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
 517 {
 518   //
 519   // Find the (cached) maximum weight for FMD<i>dr</i> in
 520   // @f$\eta@f$ bin @a iEta
 521   //
 522   // Parameters:
 523   //    d     Detector
 524   //    r     Ring
 525   //    iEta  Eta bin
 526   //
 527   // Return:
 528   //    max weight or <= 0 in case of problems
 529   //
 530   if (iEta < 0) return -1;
 531
 532   const TArrayI* max  = 0;
 533   switch (d) {
 534   case 1:  max = &fFMD1iMax;                                       break;
 535   case 2:  max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
 536   case 3:  max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
 537   }
 538   if (!max) {
 539     AliWarning(Form("No array for FMD%d%c", d, r));
 540     return -1;
 541   }
 542
 543   if (iEta >= max->fN) {
 544     AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
 545                     iEta, max->fN-1));
 546     return -1;
 547   }
 548
 549   AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
 550                    max->At(iEta)));
 551   return max->At(iEta);
 552 }
 553
 554 //_____________________________________________________________________
 555 Int_t
 556 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
 557 {
 558   //
 559   // Find the (cached) maximum weight for FMD<i>dr</i> iat
 560   // @f$\eta@f$
 561   //
 562   // Parameters:
 563   //    d     Detector
 564   //    r     Ring
 565   //    eta   Eta bin
 566   //
 567   // Return:
 568   //    max weight or <= 0 in case of problems
 569   //
 570   AliForwardCorrectionManager&  fcm  = AliForwardCorrectionManager::Instance();
 571   Int_t                         iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
 572
 573   return GetMaxWeight(d, r, iEta);
 574 }
 575
 576 //_____________________________________________________________________
 577 Float_t
 578 AliFMDDensityCalculator::NParticles(Float_t  mult,
 579                                     UShort_t d,
 580                                     Char_t   r,
 581                                     UShort_t /*s*/,
 582                                     UShort_t /*t*/,
 583                                     UShort_t /*v*/,
 584                                     Float_t  eta,
 585                                     Bool_t   lowFlux) const
 586 {
 587   //
 588   // Get the number of particles corresponding to the signal mult
 589   //
 590   // Parameters:
 591   //    mult     Signal
 592   //    d        Detector
 593   //    r        Ring
 594   //    s        Sector
 595   //    t        Strip (not used)
 596   //    v        Vertex bin
 597   //    eta      Pseudo-rapidity
 598   //    lowFlux  Low-flux flag
 599   //
 600   // Return:
 601   //    The number of particles
 602   //
 603   // if (mult <= GetMultCut()) return 0;
 604   if (lowFlux) return 1;
 605
 606   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 607   AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
 608   if (!fit) {
 609     AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
 610     return 0;
 611   }
 612
 613   Int_t    m   = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
 614   if (m < 1) {
 615     AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
 616     return 0;
 617   }
 618   UShort_t n   = TMath::Min(fMaxParticles, UShort_t(m));
 619   Double_t ret = fit->EvaluateWeighted(mult, n);
 620
 621   if (fDebug > 10) {
 622     AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
 623   }
 624
 625   fWeightedSum->Fill(ret);
 626   fSumOfWeights->Fill(ret);
 627
 628   return ret;
 629 }
 630
 631 //_____________________________________________________________________
 632 Float_t
 633 AliFMDDensityCalculator::Correction(UShort_t d,
 634                                     Char_t   r,
 635                                     UShort_t /*s*/,
 636                                     UShort_t t,
 637                                     UShort_t /*v*/,
 638                                     Float_t  eta,
 639                                     Bool_t   lowFlux) const
 640 {
 641   //
 642   // Get the inverse correction factor.  This consist of
 643   //
 644   // - acceptance correction (corners of sensors)
 645   // - double hit correction (for low-flux events)
 646   // - dead strip correction
 647   //
 648   // Parameters:
 649   //    d        Detector
 650   //    r        Ring
 651   //    s        Sector
 652   //    t        Strip (not used)
 653   //    v        Vertex bin
 654   //    eta      Pseudo-rapidity
 655   //    lowFlux  Low-flux flag
 656   //
 657   // Return:
 658   //
 659   //
 660   AliForwardCorrectionManager&  fcm = AliForwardCorrectionManager::Instance();
 661
 662   Float_t correction = 1;
 663   if (fUsePhiAcceptance) correction = AcceptanceCorrection(r,t);
 664   if (lowFlux) {
 665     TH1D* dblHitCor = 0;
 666     if (fcm.GetDoubleHit())
 667       dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
 668
 669     if (dblHitCor) {
 670       Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
 671       if (dblC > 0) correction *= dblC;
 672     }
 673     else {
 674       AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
 675     }
 676   }
 677   return correction;
 678 }
 679
 680 //_____________________________________________________________________
 681 TH1D*
 682 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
 683 {
 684   //
 685   // Generate the acceptance corrections
 686   //
 687   // Parameters:
 688   //    r Ring to generate for
 689   //
 690   // Return:
 691   //    Newly allocated histogram of acceptance corrections
 692   //
 693   const Double_t ic1[] = { 4.9895, 15.3560 };
 694   const Double_t ic2[] = { 1.8007, 17.2000 };
 695   const Double_t oc1[] = { 4.2231, 26.6638 };
 696   const Double_t oc2[] = { 1.8357, 27.9500 };
 697   const Double_t* c1   = (r == 'I' || r == 'i' ? ic1      : oc1);
 698   const Double_t* c2   = (r == 'I' || r == 'i' ? ic2      : oc2);
 699   Double_t  minR       = (r == 'I' || r == 'i' ?   4.5213 :  15.4);
 700   Double_t  maxR       = (r == 'I' || r == 'i' ?  17.2    :  28.0);
 701   Int_t     nStrips    = (r == 'I' || r == 'i' ? 512      : 256);
 702   Int_t     nSec       = (r == 'I' || r == 'i' ?  20      :  40);
 703   Float_t   basearc    = 2 * TMath::Pi() / nSec;
 704   Double_t  rad        = maxR - minR;
 705   Float_t   segment    = rad / nStrips;
 706   Float_t   cr         = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
 707
 708   // Numbers used to find end-point of strip.
 709   // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
 710   Float_t D            = c1[0] * c2[1] - c1[1] * c2[0];
 711   Float_t dx           = c2[0] - c1[0];
 712   Float_t dy           = c2[1] - c1[1];
 713   Float_t dr           = TMath::Sqrt(dx*dx+dy*dy);
 714
 715   TH1D* ret = new TH1D(Form("acc%c", r),
 716                        Form("Acceptance correction for FMDx%c", r),
 717                        nStrips, -.5, nStrips-.5);
 718   ret->SetXTitle("Strip");
 719   ret->SetYTitle("#varphi acceptance");
 720   ret->SetDirectory(0);
 721   ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
 722   ret->SetFillStyle(3001);
 723
 724   for (Int_t t = 0; t < nStrips; t++) {
 725     Float_t   radius     = minR + t * segment;
 726
 727     // If the radius of the strip is smaller than the radius corresponding
 728     // to the first corner we have a full strip length
 729     if (radius <= cr) {
 730       ret->SetBinContent(t+1, 1);
 731       continue;
 732     }
 733
 734     // Next, we should find the end-point of the strip - that is,
 735     // the coordinates where the line from c1 to c2 intersects a circle
 736     // with radius given by the strip.
 737     // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
 738     // Calculate the determinant
 739     Float_t det = radius * radius * dr * dr - D*D;
 740
 741     if (det <=  0) {
 742       // <0 means No intersection
 743       // =0 means Exactly tangent
 744       ret->SetBinContent(t+1, 1);
 745       continue;
 746     }
 747
 748     // Calculate end-point and the corresponding opening angle
 749     Float_t x   = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
 750     Float_t y   = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
 751     Float_t th  = TMath::ATan2(x, y);
 752
 753     ret->SetBinContent(t+1, th / basearc);
 754   }
 755   return ret;
 756 }
 757
 758 //_____________________________________________________________________
 759 Float_t
 760 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
 761 {
 762   //
 763   // Get the acceptance correction for strip @a t in an ring of type @a r
 764   //
 765   // Parameters:
 766   //    r  Ring type ('I' or 'O')
 767   //    t  Strip number
 768   //
 769   // Return:
 770   //    Inverse acceptance correction
 771   //
 772   TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
 773   return acc->GetBinContent(t+1);
 774 }
 775
 776 //____________________________________________________________________
 777 void
 778 AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents)
 779 {
 780   //
 781   // Scale the histograms to the total number of events
 782   //
 783   // Parameters:
 784   //    dir     where to put the output
 785   //    nEvents Number of events
 786   //
 787   if (nEvents <= 0) return;
 788   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 789   if (!d) return;
 790
 791   TIter    next(&fRingHistos);
 792   RingHistos* o = 0;
 793   THStack* sums = new THStack("sums", "sums of ring signals");
 794   while ((o = static_cast<RingHistos*>(next()))) {
 795     o->ScaleHistograms(d, nEvents);
 796     TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
 797                                          o->fDensity->GetNbinsY(),"e");
 798     sum->Scale(1., "width");
 799     sum->SetTitle(o->GetName());
 800     sum->SetDirectory(0);
 801     sum->SetYTitle("#sum N_{ch,incl}");
 802     sums->Add(sum);
 803   }
 804   d->Add(sums);
 805 }
 806
 807 //____________________________________________________________________
 808 void
 809 AliFMDDensityCalculator::DefineOutput(TList* dir)
 810 {
 811   //
 812   // Output diagnostic histograms to directory
 813   //
 814   // Parameters:
 815   //    dir List to write in
 816   //
 817   TList* d = new TList;
 818   d->SetOwner();
 819   d->SetName(GetName());
 820   dir->Add(d);
 821   d->Add(fWeightedSum);
 822   d->Add(fSumOfWeights);
 823   d->Add(fCorrections);
 824   d->Add(fAccI);
 825   d->Add(fAccO);
 826   d->Add(fMaxWeights);
 827   d->Add(fLowCuts);
 828
 829   TIter    next(&fRingHistos);
 830   RingHistos* o = 0;
 831   while ((o = static_cast<RingHistos*>(next()))) {
 832     o->Output(d);
 833   }
 834 }
 835 //____________________________________________________________________
 836 void
 837 AliFMDDensityCalculator::Print(Option_t* option) const
 838 {
 839   //
 840   // Print information
 841   //
 842   // Parameters:
 843   //    option Not used
 844   //
 845   char ind[gROOT->GetDirLevel()+3];
 846   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 847   ind[gROOT->GetDirLevel()] = '\0';
 848   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 849             << std::boolalpha
 850             << ind << " Multiplicity cut:       " << fMultCut << '\n'
 851             << ind << " # of (xi+sigma) factor: " << fNXi << '\n'
 852             << ind << " Include sigma in cut:   " << fIncludeSigma << '\n'
 853             << ind << " Low cut method:         "
 854             << (fMultCut > 0 ? "fixed" :
 855                 (fNXi >= 0 ? "xi+sigma" : "fit range")) << '\n'
 856             << ind << " Max(particles):         " << fMaxParticles << '\n'
 857             << ind << " Poisson method:         " << fUsePoisson << '\n'
 858             << ind << " Use phi acceptance:     " << fUsePhiAcceptance << '\n'
 859             << ind << " Eta lumping:            " << fEtaLumping << '\n'
 860             << ind << " Phi lumping:            " << fPhiLumping << '\n'
 861             << std::noboolalpha
 862             << std::flush;
 863   TString opt(option);
 864   opt.ToLower();
 865   if (opt.Contains("nomax")) return;
 866
 867   std::cout << ind << " Max weights:\n";
 868
 869   for (UShort_t d=1; d<=3; d++) {
 870     UShort_t nr = (d == 1 ? 1 : 2);
 871     for (UShort_t q=0; q<nr; q++) {
 872       ind[gROOT->GetDirLevel()]   = ' ';
 873       ind[gROOT->GetDirLevel()+1] = '\0';
 874       Char_t      r = (q == 0 ? 'I' : 'O');
 875       std::cout << ind << " FMD" << d << r << ":";
 876       ind[gROOT->GetDirLevel()+1] = ' ';
 877       ind[gROOT->GetDirLevel()+2] = '\0';
 878
 879       const TArrayI& a = (d == 1 ? fFMD1iMax :
 880                           (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
 881                            (r == 'I' ? fFMD3iMax : fFMD3oMax)));
 882       Int_t j = 0;
 883       for (Int_t i = 0; i < a.fN; i++) {
 884         if (a.fArray[i] < 1) continue;
 885         if (j % 6 == 0)      std::cout << "\n " << ind;
 886         j++;
 887         std::cout << "  " << std::setw(3) << i << ": " << a.fArray[i];
 888       }
 889       std::cout << std::endl;
 890     }
 891   }
 892 }
 893
 894 //====================================================================
 895 AliFMDDensityCalculator::RingHistos::RingHistos()
 896   : AliForwardUtil::RingHistos(),
 897     fEvsN(0),
 898     fEvsM(0),
 899     fEtaVsN(0),
 900     fEtaVsM(0),
 901     fCorr(0),
 902     fDensity(0),
 903     fELossVsPoisson(0),
 904     fTotalStrips(0),
 905     fEmptyStrips(0),
 906     fBasicHits(0),
 907     fEmptyVsTotal(0)
 908 {
 909   //
 910   // Default CTOR
 911   //
 912 }
 913
 914 //____________________________________________________________________
 915 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
 916   : AliForwardUtil::RingHistos(d,r),
 917     fEvsN(0),
 918     fEvsM(0),
 919     fEtaVsN(0),
 920     fEtaVsM(0),
 921     fCorr(0),
 922     fDensity(0),
 923     fELossVsPoisson(0),
 924     fTotalStrips(0),
 925     fEmptyStrips(0),
 926     fBasicHits(0),
 927     fEmptyVsTotal(0)
 928 {
 929   //
 930   // Constructor
 931   //
 932   // Parameters:
 933   //    d detector
 934   //    r ring
 935   //
 936   fEvsN = new TH2D("elossVsNnocorr",
 937                    "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
 938                    250, -.5, 24.5, 250, -.5, 24.5);
 939   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
 940   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
 941   fEvsN->Sumw2();
 942   fEvsN->SetDirectory(0);
 943
 944   fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
 945   fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
 946   fEvsM->SetDirectory(0);
 947
 948   fEtaVsN = new TProfile("etaVsNnocorr",
 949                          "Average inclusive N_{ch} vs #eta (uncorrected)",
 950                          200, -4, 6);
 951   fEtaVsN->SetXTitle("#eta");
 952   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
 953   fEtaVsN->SetDirectory(0);
 954   fEtaVsN->SetLineColor(Color());
 955   fEtaVsN->SetFillColor(Color());
 956
 957   fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
 958   fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
 959   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
 960   fEtaVsM->SetDirectory(0);
 961
 962
 963   fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
 964   fCorr->SetXTitle("#eta");
 965   fCorr->SetYTitle("#LT correction#GT");
 966   fCorr->SetDirectory(0);
 967   fCorr->SetLineColor(Color());
 968   fCorr->SetFillColor(Color());
 969
 970   fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
 971                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
 972                       0, 2*TMath::Pi());
 973   fDensity->SetDirectory(0);
 974   fDensity->Sumw2();
 975   fDensity->SetMarkerColor(Color());
 976   fDensity->SetXTitle("#eta");
 977   fDensity->SetYTitle("#phi [radians]");
 978   fDensity->SetZTitle("Inclusive N_{ch} density");
 979
 980   fELossVsPoisson = new TH2D("elossVsPoisson",
 981                              "N_{ch} from energy loss vs from Poission",
 982                              150, 0, 30, 150, 0, 30);
 983   fELossVsPoisson->SetDirectory(0);
 984   fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
 985   fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
 986   fELossVsPoisson->SetZTitle("Correlation");
 987
 988   Int_t nStrips = 100;
 989   fEmptyVsTotal = new TH2D("emptyVsTotal",
 990                            "# of empty strips vs. total # strips",
 991                            nStrips+1, -.5, nStrips+.5,
 992                            nStrips+1, -.5, nStrips+.5);
 993   fEmptyVsTotal->SetDirectory(0);
 994   fEmptyVsTotal->SetXTitle("Total # strips");
 995   fEmptyVsTotal->SetYTitle("Empty # strips");
 996   fEmptyVsTotal->SetZTitle("Correlation");
 997 }
 998 //____________________________________________________________________
 999 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
1000   : AliForwardUtil::RingHistos(o),
1001     fEvsN(o.fEvsN),
1002     fEvsM(o.fEvsM),
1003     fEtaVsN(o.fEtaVsN),
1004     fEtaVsM(o.fEtaVsM),
1005     fCorr(o.fCorr),
1006     fDensity(o.fDensity),
1007     fELossVsPoisson(o.fELossVsPoisson),
1008     fTotalStrips(o.fTotalStrips),
1009     fEmptyStrips(o.fEmptyStrips),
1010     fBasicHits(o.fBasicHits),
1011     fEmptyVsTotal(o.fEmptyVsTotal)
1012 {
1013   //
1014   // Copy constructor
1015   //
1016   // Parameters:
1017   //    o Object to copy from
1018   //
1019 }
1020
1021 //____________________________________________________________________
1022 AliFMDDensityCalculator::RingHistos&
1023 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1024 {
1025   //
1026   // Assignment operator
1027   //
1028   // Parameters:
1029   //    o Object to assign from
1030   //
1031   // Return:
1032   //    Reference to this
1033   //
1034   AliForwardUtil::RingHistos::operator=(o);
1035
1036   if (fEvsN)           delete fEvsN;
1037   if (fEvsM)           delete fEvsM;
1038   if (fEtaVsN)         delete fEtaVsN;
1039   if (fEtaVsM)         delete fEtaVsM;
1040   if (fCorr)           delete fCorr;
1041   if (fDensity)        delete fDensity;
1042   if (fELossVsPoisson) delete fELossVsPoisson;
1043   if (fTotalStrips)    delete fTotalStrips;
1044   if (fEmptyStrips)    delete fEmptyStrips;
1045
1046   fEvsN           = static_cast<TH2D*>(o.fEvsN->Clone());
1047   fEvsM           = static_cast<TH2D*>(o.fEvsM->Clone());
1048   fEtaVsN         = static_cast<TProfile*>(o.fEtaVsN->Clone());
1049   fEtaVsM         = static_cast<TProfile*>(o.fEtaVsM->Clone());
1050   fCorr           = static_cast<TProfile*>(o.fCorr->Clone());
1051   fDensity        = static_cast<TH2D*>(o.fDensity->Clone());
1052   fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson);
1053   fTotalStrips    = static_cast<TH2D*>(o.fTotalStrips);
1054   fEmptyStrips    = static_cast<TH2D*>(o.fEmptyStrips);
1055
1056   return *this;
1057 }
1058 //____________________________________________________________________
1059 AliFMDDensityCalculator::RingHistos::~RingHistos()
1060 {
1061   //
1062   // Destructor
1063   //
1064   if (fEvsN)           delete fEvsN;
1065   if (fEvsM)           delete fEvsM;
1066   if (fEtaVsN)         delete fEtaVsN;
1067   if (fEtaVsM)         delete fEtaVsM;
1068   if (fCorr)           delete fCorr;
1069   if (fDensity)        delete fDensity;
1070   if (fELossVsPoisson) delete fELossVsPoisson;
1071   if (fTotalStrips)    delete fTotalStrips;
1072   if (fEmptyStrips)    delete fEmptyStrips;
1073 }
1074
1075 //____________________________________________________________________
1076 void
1077 AliFMDDensityCalculator::RingHistos::ResetPoissonHistos(const TH2D* h,
1078                                                         Int_t etaLumping,
1079                                                         Int_t phiLumping)
1080 {
1081   if (fTotalStrips) {
1082     fTotalStrips->Reset();
1083     fEmptyStrips->Reset();
1084     fBasicHits->Reset();
1085     return;
1086   }
1087   //Inserted by HHD
1088
1089   Int_t    nEta    = h->GetNbinsX() / etaLumping;
1090   Int_t    nEtaF   = h->GetNbinsX();
1091   Double_t etaMin  = h->GetXaxis()->GetXmin();
1092   Double_t etaMax  = h->GetXaxis()->GetXmax();
1093   Int_t    nPhi    = h->GetNbinsY() / phiLumping;
1094   Int_t    nPhiF   = h->GetNbinsY();
1095   Double_t phiMin  = h->GetYaxis()->GetXmin();
1096   Double_t phiMax  = h->GetYaxis()->GetXmax();
1097
1098   fTotalStrips = new TH2D(Form("total%s", fName.Data()),
1099                           Form("Total number of strips in %s", fName.Data()),
1100                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
1101   fEmptyStrips = new TH2D(Form("empty%s", fName.Data()),
1102                           Form("Empty number of strips in %s", fName.Data()),
1103                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
1104   fBasicHits   = new TH2D(Form("basichits%s", fName.Data()),
1105                           Form("Basic number of hits in %s", fName.Data()),
1106                           nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
1107
1108   fTotalStrips->SetDirectory(0);
1109   fEmptyStrips->SetDirectory(0);
1110   fBasicHits->SetDirectory(0);
1111   fTotalStrips->SetXTitle("#eta");
1112   fEmptyStrips->SetXTitle("#eta");
1113   fBasicHits->SetXTitle("#eta");
1114   fTotalStrips->SetYTitle("#varphi [radians]");
1115   fEmptyStrips->SetYTitle("#varphi [radians]");
1116   fBasicHits->SetYTitle("#varphi [radians]");
1117   fTotalStrips->Sumw2();
1118   fEmptyStrips->Sumw2();
1119   fBasicHits->Sumw2();
1120 }
1121
1122 //____________________________________________________________________
1123 void
1124 AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
1125 {
1126 }
1127
1128 //____________________________________________________________________
1129 void
1130 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
1131 {
1132   //
1133   // Make output
1134   //
1135   // Parameters:
1136   //    dir Where to put it
1137   //
1138   TList* d = DefineOutputList(dir);
1139   d->Add(fEvsN);
1140   d->Add(fEvsM);
1141   d->Add(fEtaVsN);
1142   d->Add(fEtaVsM);
1143   d->Add(fCorr);
1144   d->Add(fDensity);
1145   d->Add(fELossVsPoisson);
1146   d->Add(fEmptyVsTotal);
1147   d->Add(fTotalStrips);
1148   d->Add(fEmptyStrips);
1149   d->Add(fBasicHits);
1150
1151
1152 }
1153
1154 //____________________________________________________________________
1155 void
1156 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1157 {
1158   //
1159   // Scale the histograms to the total number of events
1160   //
1161   // Parameters:
1162   //    dir     Where the output is
1163   //    nEvents Number of events
1164   //
1165   TList* l = GetOutputList(dir);
1166   if (!l) return;
1167
1168   TH1* density = GetOutputHist(l,"inclDensity");
1169   if (density) density->Scale(1./nEvents);
1170 }
1171
1172 //____________________________________________________________________
1173 //
1174 // EOF
1175 //
1176
1177
1178