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