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