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