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