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