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