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 TMath::Min(Int_t(fMaxParticles), 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 << ClassName() << ": " << GetName() << '\n'
 718             << ind << " Multiplicity cut:       " << fMultCut << '\n'
 719             << ind << " Max(particles):         " << fMaxParticles << '\n'
 720             << ind << " Eta lumping:            " << fEtaLumping << '\n'
 721             << ind << " Phi lumping:            " << fPhiLumping << '\n'
 722             << std::flush;
 723   TString opt(option);
 724   opt.ToLower();
 725   if (opt.Contains("nomax")) return;
 726
 727   std::cout << ind << " Max weights:\n";
 728
 729   for (UShort_t d=1; d<=3; d++) {
 730     UShort_t nr = (d == 1 ? 1 : 2);
 731     for (UShort_t q=0; q<nr; q++) {
 732       ind[gROOT->GetDirLevel()]   = ' ';
 733       ind[gROOT->GetDirLevel()+1] = '\0';
 734       Char_t      r = (q == 0 ? 'I' : 'O');
 735       std::cout << ind << " FMD" << d << r << ":";
 736       ind[gROOT->GetDirLevel()+1] = ' ';
 737       ind[gROOT->GetDirLevel()+2] = '\0';
 738
 739       const TArrayI& a = (d == 1 ? fFMD1iMax :
 740                           (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
 741                            (r == 'I' ? fFMD3iMax : fFMD3oMax)));
 742       Int_t j = 0;
 743       for (Int_t i = 0; i < a.fN; i++) {
 744         if (a.fArray[i] < 1) continue;
 745         if (j % 6 == 0)      std::cout << "\n " << ind;
 746         j++;
 747         std::cout << "  " << std::setw(3) << i << ": " << a.fArray[i];
 748       }
 749       std::cout << std::endl;
 750     }
 751   }
 752 }
 753
 754 //====================================================================
 755 AliFMDDensityCalculator::RingHistos::RingHistos()
 756   : AliForwardUtil::RingHistos(),
 757     fEvsN(0),
 758     fEvsM(0),
 759     fEtaVsN(0),
 760     fEtaVsM(0),
 761     fCorr(0),
 762     fDensity(0),
 763     fELossVsPoisson(0),
 764     fTotalStrips(0),
 765     fEmptyStrips(0),
 766     fBasicHits(0),
 767     fEmptyVsTotal(0)
 768 {
 769   //
 770   // Default CTOR
 771   //
 772 }
 773
 774 //____________________________________________________________________
 775 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
 776   : AliForwardUtil::RingHistos(d,r),
 777     fEvsN(0),
 778     fEvsM(0),
 779     fEtaVsN(0),
 780     fEtaVsM(0),
 781     fCorr(0),
 782     fDensity(0),
 783     fELossVsPoisson(0),
 784     fTotalStrips(0),
 785     fEmptyStrips(0),
 786     fBasicHits(0),
 787     fEmptyVsTotal(0)
 788 {
 789   //
 790   // Constructor
 791   //
 792   // Parameters:
 793   //    d detector
 794   //    r ring
 795   //
 796   fEvsN = new TH2D(Form("%s_Eloss_N_nocorr", fName.Data()),
 797                    Form("#Delta E/#Delta E_{mip} vs uncorrected inclusive "
 798                         "N_{ch} in %s", fName.Data()),
 799                    2500, -.5, 24.5, 2500, -.5, 24.5);
 800   fEvsM = new TH2D(Form("%s_Eloss_N_corr", fName.Data()),
 801                    Form("#Delta E/#Delta E_{mip} vs corrected inclusive "
 802                         "N_{ch} in %s", fName.Data()),
 803                    2500, -.5, 24.5, 2500, -.5, 24.5);
 804   fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
 805   fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
 806   fEvsN->Sumw2();
 807   fEvsN->SetDirectory(0);
 808   fEvsM->SetXTitle("#Delta E/#Delta E_{mip}");
 809   fEvsM->SetYTitle("Inclusive N_{ch} (corrected)");
 810   fEvsM->Sumw2();
 811   fEvsM->SetDirectory(0);
 812
 813   fEtaVsN = new TProfile(Form("%s_Eta_N_nocorr", fName.Data()),
 814                          Form("Average inclusive N_{ch} vs #eta (uncorrected) "
 815                               "in %s", fName.Data()), 200, -4, 6);
 816   fEtaVsM = new TProfile(Form("%s_Eta_N_corr", fName.Data()),
 817                          Form("Average inclusive N_{ch} vs #eta (corrected) "
 818                               "in %s", fName.Data()), 200, -4, 6);
 819   fEtaVsN->SetXTitle("#eta");
 820   fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
 821   fEtaVsN->SetDirectory(0);
 822   fEtaVsN->SetLineColor(Color());
 823   fEtaVsN->SetFillColor(Color());
 824   fEtaVsM->SetXTitle("#eta");
 825   fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
 826   fEtaVsM->SetDirectory(0);
 827   fEtaVsM->SetLineColor(Color());
 828   fEtaVsM->SetFillColor(Color());
 829
 830
 831   fCorr = new TProfile(Form("%s_corr", fName.Data()),
 832                          Form("Average correction in %s", fName.Data()),
 833                        200, -4, 6);
 834   fCorr->SetXTitle("#eta");
 835   fCorr->SetYTitle("#LT correction#GT");
 836   fCorr->SetDirectory(0);
 837   fCorr->SetLineColor(Color());
 838   fCorr->SetFillColor(Color());
 839
 840   fDensity = new TH2D(Form("%s_Incl_Density", fName.Data()),
 841                       Form("Inclusive N_{ch} density in %s", fName.Data()),
 842                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
 843                       0, 2*TMath::Pi());
 844   fDensity->SetDirectory(0);
 845   fDensity->SetXTitle("#eta");
 846   fDensity->SetYTitle("#phi [radians]");
 847   fDensity->SetZTitle("Inclusive N_{ch} density");
 848
 849   fELossVsPoisson = new TH2D(Form("%s_eloss_vs_poisson", fName.Data()),
 850                              Form("N_{ch} from energy loss vs from Poission %s",
 851                                   fName.Data()), 100, 0, 20, 100, 0, 20);
 852   fELossVsPoisson->SetDirectory(0);
 853   fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
 854   fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
 855   fELossVsPoisson->SetZTitle("Correlation");
 856
 857   fEmptyVsTotal = new TH2D(Form("%s_empty_vs_total", fName.Data()),
 858                            Form("# of empty strips vs. total @ # strips in %s",
 859                                 fName.Data()), 21, -.5, 20.5, 21, -0.5, 20.5);
 860   fEmptyVsTotal->SetDirectory(0);
 861   fEmptyVsTotal->SetXTitle("Total # strips");
 862   fEmptyVsTotal->SetYTitle("Empty # strips");
 863   fEmptyVsTotal->SetZTitle("Correlation");
 864 }
 865 //____________________________________________________________________
 866 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
 867   : AliForwardUtil::RingHistos(o),
 868     fEvsN(o.fEvsN),
 869     fEvsM(o.fEvsM),
 870     fEtaVsN(o.fEtaVsN),
 871     fEtaVsM(o.fEtaVsM),
 872     fCorr(o.fCorr),
 873     fDensity(o.fDensity),
 874     fELossVsPoisson(o.fELossVsPoisson),
 875     fTotalStrips(o.fTotalStrips),
 876     fEmptyStrips(o.fEmptyStrips),
 877     fBasicHits(o.fBasicHits),
 878     fEmptyVsTotal(o.fEmptyVsTotal)
 879 {
 880   //
 881   // Copy constructor
 882   //
 883   // Parameters:
 884   //    o Object to copy from
 885   //
 886 }
 887
 888 //____________________________________________________________________
 889 AliFMDDensityCalculator::RingHistos&
 890 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
 891 {
 892   //
 893   // Assignment operator
 894   //
 895   // Parameters:
 896   //    o Object to assign from
 897   //
 898   // Return:
 899   //    Reference to this
 900   //
 901   AliForwardUtil::RingHistos::operator=(o);
 902
 903   if (fEvsN)           delete fEvsN;
 904   if (fEvsM)           delete fEvsM;
 905   if (fEtaVsN)         delete fEtaVsN;
 906   if (fEtaVsM)         delete fEtaVsM;
 907   if (fCorr)           delete fCorr;
 908   if (fDensity)        delete fDensity;
 909   if (fELossVsPoisson) delete fELossVsPoisson;
 910   if (fTotalStrips)    delete fTotalStrips;
 911   if (fEmptyStrips)    delete fEmptyStrips;
 912
 913   fEvsN           = static_cast<TH2D*>(o.fEvsN->Clone());
 914   fEvsM           = static_cast<TH2D*>(o.fEvsM->Clone());
 915   fEtaVsN         = static_cast<TProfile*>(o.fEtaVsN->Clone());
 916   fEtaVsM         = static_cast<TProfile*>(o.fEtaVsM->Clone());
 917   fCorr           = static_cast<TProfile*>(o.fCorr->Clone());
 918   fDensity        = static_cast<TH2D*>(o.fDensity->Clone());
 919   fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson);
 920   fTotalStrips    = static_cast<TH2D*>(o.fTotalStrips);
 921   fEmptyStrips    = static_cast<TH2D*>(o.fEmptyStrips);
 922
 923   return *this;
 924 }
 925 //____________________________________________________________________
 926 AliFMDDensityCalculator::RingHistos::~RingHistos()
 927 {
 928   //
 929   // Destructor
 930   //
 931   if (fEvsN)           delete fEvsN;
 932   if (fEvsM)           delete fEvsM;
 933   if (fEtaVsN)         delete fEtaVsN;
 934   if (fEtaVsM)         delete fEtaVsM;
 935   if (fCorr)           delete fCorr;
 936   if (fDensity)        delete fDensity;
 937   if (fELossVsPoisson) delete fELossVsPoisson;
 938   if (fTotalStrips)    delete fTotalStrips;
 939   if (fEmptyStrips)    delete fEmptyStrips;
 940 }
 941
 942 //____________________________________________________________________
 943 void
 944 AliFMDDensityCalculator::RingHistos::ResetPoissonHistos(const TH2D* h,
 945                                                         Int_t etaLumping,
 946                                                         Int_t phiLumping)
 947 {
 948   if (fTotalStrips) {
 949     fTotalStrips->Reset();
 950     fEmptyStrips->Reset();
 951     fBasicHits->Reset();
 952     return;
 953   }
 954   //Inserted by HHD
 955
 956   Int_t    nEta    = h->GetNbinsX() / etaLumping;
 957   Int_t    nEtaF   = h->GetNbinsX();
 958   Double_t etaMin  = h->GetXaxis()->GetXmin();
 959   Double_t etaMax  = h->GetXaxis()->GetXmax();
 960   Int_t    nPhi    = h->GetNbinsY() / phiLumping;
 961   Int_t    nPhiF   = h->GetNbinsY();
 962   Double_t phiMin  = h->GetYaxis()->GetXmin();
 963   Double_t phiMax  = h->GetYaxis()->GetXmax();
 964
 965   fTotalStrips = new TH2D(Form("total%s", fName.Data()),
 966                           Form("Total number of strips in %s", fName.Data()),
 967                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
 968   fEmptyStrips = new TH2D(Form("empty%s", fName.Data()),
 969                           Form("Empty number of strips in %s", fName.Data()),
 970                           nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
 971   fBasicHits   = new TH2D(Form("basichits%s", fName.Data()),
 972                           Form("Basic number of hits in %s", fName.Data()),
 973                           nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
 974
 975   fTotalStrips->SetDirectory(0);
 976   fEmptyStrips->SetDirectory(0);
 977   fBasicHits->SetDirectory(0);
 978   fTotalStrips->SetXTitle("#eta");
 979   fEmptyStrips->SetXTitle("#eta");
 980   fBasicHits->SetXTitle("#eta");
 981   fTotalStrips->SetYTitle("#varphi [radians]");
 982   fEmptyStrips->SetYTitle("#varphi [radians]");
 983   fBasicHits->SetYTitle("#varphi [radians]");
 984   fTotalStrips->Sumw2();
 985   fEmptyStrips->Sumw2();
 986   fBasicHits->Sumw2();
 987 }
 988
 989 //____________________________________________________________________
 990 void
 991 AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
 992 {
 993 }
 994
 995 //____________________________________________________________________
 996 void
 997 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
 998 {
 999   //
1000   // Make output
1001   //
1002   // Parameters:
1003   //    dir Where to put it
1004   //
1005   TList* d = DefineOutputList(dir);
1006   d->Add(fEvsN);
1007   d->Add(fEvsM);
1008   d->Add(fEtaVsN);
1009   d->Add(fEtaVsM);
1010   d->Add(fCorr);
1011   d->Add(fDensity);
1012   d->Add(fELossVsPoisson);
1013   d->Add(fEmptyVsTotal);
1014   d->Add(fTotalStrips);
1015   d->Add(fEmptyStrips);
1016   d->Add(fBasicHits);
1017
1018
1019 }
1020
1021 //____________________________________________________________________
1022 void
1023 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1024 {
1025   //
1026   // Scale the histograms to the total number of events
1027   //
1028   // Parameters:
1029   //    dir     Where the output is
1030   //    nEvents Number of events
1031   //
1032   TList* l = GetOutputList(dir);
1033   if (!l) return;
1034
1035   TH1* density = GetOutputHist(l,Form("%s_Incl_Density", fName.Data()));
1036   if (density) density->Scale(1./nEvents);
1037 }
1038
1039 //____________________________________________________________________
1040 //
1041 // EOF
1042 //
1043
1044
1045