PWG2/FORWARD/analysis2/AliFMDDensityCalculator.cxx

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