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