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