PWGLF/FORWARD/analysis2/AliFMDCorrAcceptance.cxx

   1 //
   2 // This class contains the acceptance correction due to dead channels
   3 //
   4 //
   5 #include "AliFMDCorrAcceptance.h"
   6 #include <TBrowser.h>
   7 #include <TH2D.h>
   8 #include <AliLog.h>
   9 #include <AliForwardUtil.h>
  10 #include <TPad.h>
  11 #include <TCanvas.h>
  12 #include <TLatex.h>
  13 #include <TMath.h>
  14 #include <THStack.h>
  15 #include <TROOT.h>
  16 #include <iostream>
  17
  18 //____________________________________________________________________
  19 AliFMDCorrAcceptance::AliFMDCorrAcceptance()
  20   : fRingArray(),
  21     fCache(0),
  22     fVertexAxis(0,0,0),
  23     fHasOverflow(false)
  24 {
  25   //
  26   // Default constructor
  27   //
  28   fRingArray.SetOwner(kTRUE);
  29   fRingArray.SetName("rings");
  30   fVertexAxis.SetName("vtxAxis");
  31   fVertexAxis.SetTitle("v_{z} [cm]");
  32
  33 }
  34 //____________________________________________________________________
  35 AliFMDCorrAcceptance::AliFMDCorrAcceptance(const
  36                                                AliFMDCorrAcceptance& o)
  37   : TObject(o),
  38     fRingArray(o.fRingArray),
  39     fCache(o.fCache),
  40     fVertexAxis(o.fVertexAxis.GetNbins(), o.fVertexAxis.GetXmin(),
  41                 o.fVertexAxis.GetXmax()),
  42     fHasOverflow(o.fHasOverflow)
  43 {
  44   //
  45   // Copy constructor
  46   //
  47   // Parameters:
  48   //    o Object to copy from
  49   //
  50   fVertexAxis.SetName("vtxAxis");
  51   fVertexAxis.SetTitle("v_{z} [cm]");
  52 }
  53 //____________________________________________________________________
  54 AliFMDCorrAcceptance::~AliFMDCorrAcceptance()
  55 {
  56   //
  57   // Destructor
  58   //
  59   //
  60   fRingArray.Clear();
  61   if (fCache) fCache->Clear();
  62 }
  63 //____________________________________________________________________
  64 AliFMDCorrAcceptance&
  65 AliFMDCorrAcceptance::operator=(const AliFMDCorrAcceptance& o)
  66 {
  67   //
  68   // Assignment operator
  69   //
  70   // Parameters:
  71   //    o Object to assign from
  72   //
  73   // Return:
  74   //    Reference to this object
  75   //
  76   if (&o == this) return *this;
  77
  78   fRingArray        = o.fRingArray;
  79   fCache            = o.fCache;
  80   fHasOverflow      = o.fHasOverflow;
  81   SetVertexAxis(o.fVertexAxis);
  82
  83   return *this;
  84 }
  85 //____________________________________________________________________
  86 TH2D*
  87 AliFMDCorrAcceptance::GetCorrection(UShort_t d, Char_t r, Double_t v) const
  88 {
  89   //
  90   // Get the acceptance correction @f$ a_{r,v}@f$
  91   //
  92   // Parameters:
  93   //    d  Detector number (1-3)
  94   //    r  Ring identifier (I or O)
  95   //    v  Primary interaction point @f$z@f$ coordinate
  96   //
  97   // Return:
  98   //    The correction @f$ a_{r,v}@f$
  99   //
 100   Int_t b = FindVertexBin(v);
 101   if (b <= 0) return 0;
 102   return GetCorrection(d, r, UShort_t(b));
 103 }
 104 //____________________________________________________________________
 105 TH2D*
 106 AliFMDCorrAcceptance::GetCorrection(UShort_t d, Char_t r, UShort_t b) const
 107 {
 108   //
 109   // Get the acceptance correction @f$ a_{r,v}@f$
 110   //
 111   // Parameters:
 112   //    d  Detector number (1-3)
 113   //    r  Ring identifier (I or O)
 114   //    b  Bin corresponding to the primary interaction point
 115   //           @f$z@f$ coordinate (1 based)
 116   //
 117   // Return:
 118   //    The correction @f$ a_{r,v}@f$
 119   //
 120   return static_cast<TH2D*>(GetObject(fRingArray, d, r, b));
 121 }
 122 //____________________________________________________________________
 123 TH1D*
 124 AliFMDCorrAcceptance::GetPhiAcceptance(UShort_t d, Char_t r, Double_t v) const
 125 {
 126   //
 127   // Get the acceptance correction @f$ a_{r,v}@f$
 128   //
 129   // Parameters:
 130   //    d  Detector number (1-3)
 131   //    r  Ring identifier (I or O)
 132   //    v  Primary interaction point @f$z@f$ coordinate
 133   //
 134   // Return:
 135   //    The correction @f$ a_{r,v}@f$
 136   //
 137   Int_t b = FindVertexBin(v);
 138   if (b <= 0) return 0;
 139   return GetPhiAcceptance(d, r, UShort_t(b));
 140 }
 141 //____________________________________________________________________
 142 TH1D*
 143 AliFMDCorrAcceptance::GetPhiAcceptance(UShort_t d, Char_t r, UShort_t b) const
 144 {
 145   //
 146   // Get the acceptance correction @f$ a_{r,v}@f$
 147   //
 148   // Parameters:
 149   //    d  Detector number (1-3)
 150   //    r  Ring identifier (I or O)
 151   //    b  Bin corresponding to the primary interaction point
 152   //           @f$z@f$ coordinate (1 based)
 153   //
 154   // Return:
 155   //    The correction @f$ a_{r,v}@f$
 156   //
 157   if (!fHasOverflow) return 0;
 158   if (!fCache) FillCache();
 159   return static_cast<TH1D*>(GetObject(*fCache, d, r, b));
 160 }
 161
 162 //____________________________________________________________________
 163 Int_t
 164 AliFMDCorrAcceptance::FindVertexBin(Double_t v) const
 165 {
 166   //
 167   // Find the vertex bin that corresponds to the passed vertex
 168   //
 169   // Parameters:
 170   //    vertex The interaction points @f$z@f$-coordinate
 171   //
 172   // Return:
 173   //    Vertex bin in @f$[1,N_{\mbox{vertex}}]@f$ or negative if
 174   // out of range
 175   //
 176   if (fVertexAxis.GetNbins() <= 0) {
 177     AliWarning("No vertex array defined");
 178     return 0;
 179   }
 180   Int_t bin = const_cast<TAxis&>(fVertexAxis).FindBin(v);
 181   if (bin <= 0 || bin > fVertexAxis.GetNbins()) {
 182     AliWarning(Form("vertex %+8.4f out of range [%+8.4f,%+8.4f]",
 183                     v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
 184     return 0;
 185   }
 186   return bin;
 187 }
 188 //____________________________________________________________________
 189 Int_t
 190 AliFMDCorrAcceptance::GetRingIndex(UShort_t d, Char_t r) const
 191 {
 192   // Get the index corresponding to the given ring
 193   //
 194   // Parameters:
 195   //    d Detector
 196   //    r Ring
 197   //
 198   // Return:
 199   //    Index (0 based) or negative in case of errors
 200   //
 201   switch (d) {
 202   case 1:  return 0;
 203   case 2:  return (r == 'I' || r == 'i' ? 1 : 2); break;
 204   case 3:  return (r == 'I' || r == 'i' ? 3 : 4); break;
 205
 206   }
 207   AliWarning(Form("Index for FMD%d%c not found", d, r));
 208   return -1;
 209 }
 210 //____________________________________________________________________
 211 TObject*
 212 AliFMDCorrAcceptance::GetObject(const TObjArray& m, UShort_t d,
 213                                 Char_t r, UShort_t b) const
 214 {
 215   //
 216   // Get the object @f$ a_{r,v}@f$
 217   //
 218   // Parameters:
 219   //    m  Mother list
 220   //    d  Detector number (1-3)
 221   //    r  Ring identifier (I or O)
 222   //    b  Bin corresponding to the primary interaction point
 223   //           @f$z@f$ coordinate (1 based)
 224   //
 225   // Return:
 226   //    The correction @f$ a_{r,v}@f$
 227   //
 228   TObjArray* ringArray = GetRingArray(m, d, r);
 229   if (!ringArray) return 0;
 230
 231   if (b <= 0 || b > ringArray->GetEntriesFast()) {
 232     AliWarning(Form("vertex bin %d out of range [1,%d]",
 233                     b, ringArray->GetEntriesFast()));
 234     return 0;
 235   }
 236
 237   TObject* o = ringArray->At(b-1);
 238   if (o) return o;
 239
 240   AliWarning(Form("No dead channels map found for FMD%d%c in vertex bin %d",
 241                     d,r,b));
 242   return 0;
 243 }
 244 //____________________________________________________________________
 245 TObjArray*
 246 AliFMDCorrAcceptance::GetRingArray(const TObjArray& m,
 247                                    UShort_t d, Char_t r) const
 248 {
 249   //
 250   // Get the ring array corresponding to the specified ring
 251   //
 252   // Parameters:
 253   //    d Detector
 254   //    r Ring
 255   //
 256   // Return:
 257   //    Pointer to ring array, or null in case of problems
 258   //
 259   Int_t idx = GetRingIndex(d,r);
 260   if (idx < 0) return 0;
 261
 262   TObject* o = m.At(idx);
 263   if (!o) {
 264     AliWarning(Form("No array found for FMD%d%c", d, r));
 265     return 0;
 266   }
 267
 268   return static_cast<TObjArray*>(o);
 269 }
 270 //____________________________________________________________________
 271 TObjArray*
 272 AliFMDCorrAcceptance::GetOrMakeRingArray(TObjArray& m,
 273                                          UShort_t d, Char_t r) const
 274 {
 275   //
 276   // Get the ring array corresponding to the specified ring
 277   //
 278   // Parameters:
 279   //    d Detector
 280   //    r Ring
 281   //
 282   // Return:
 283   //    Pointer to ring array, or newly created container
 284   //
 285   Int_t idx = GetRingIndex(d,r);
 286   if (idx < 0) return 0;
 287
 288   TObject* o = m.At(idx);
 289   if (!o) {
 290     TObjArray* a = new TObjArray(fVertexAxis.GetNbins());
 291     a->SetName(Form("FMD%d%c", d, r));
 292     a->SetOwner(kTRUE);
 293     m.AddAtAndExpand(a, idx);
 294     return a;
 295   }
 296
 297   return static_cast<TObjArray*>(m.At(idx));
 298 }
 299
 300 //____________________________________________________________________
 301 Bool_t
 302 AliFMDCorrAcceptance::SetCorrection(UShort_t d, Char_t r,
 303                                     UShort_t b, TH2D*  h)
 304 {
 305   //
 306   // Set the acceptance correction @f$ a_{r,v}(\eta)@f$
 307   // Note, that the object takes ownership of the passed pointer.
 308   //
 309   // Parameters:
 310   //    d    Detector number (1-3)
 311   //    r    Ring identifier (I or O)
 312   //    b    Bin corresponding to the primary interaction point
 313   //             @f$z@f$ coordinate  (1 based)
 314   //    h    @f$ a_{r,v}(\eta)@f$
 315   //
 316   // Return:
 317   //    true if operation succeeded
 318   //
 319   TObjArray* ringArray = GetOrMakeRingArray(fRingArray, d, r);
 320   if (!ringArray) return false;
 321
 322   if (b <= 0 || b > fVertexAxis.GetNbins()) {
 323     AliWarning(Form("Vertex bin %3d out of range [1,%3d]",
 324                     b, fVertexAxis.GetNbins()));
 325     return false;
 326   }
 327   h->SetName(Form("FMD%d%c_vtxbin%03d", d, r, b));
 328   h->SetTitle(Form("Acceptance correction for FMD%d%c "
 329                    "in vertex bin %d [%+8.4f,%+8.4f]",
 330                    d, r, b, fVertexAxis.GetBinLowEdge(b),
 331                    fVertexAxis.GetBinUpEdge(b)));
 332   h->SetXTitle("#eta");
 333   h->SetYTitle("N_{strips,OK}/N_{strips}");
 334   h->SetFillStyle(3001);
 335   h->SetDirectory(0);
 336   h->SetStats(0);
 337   ringArray->AddAtAndExpand(h, b-1);
 338   return kTRUE;
 339 }
 340 //____________________________________________________________________
 341 Bool_t
 342 AliFMDCorrAcceptance::SetCorrection(UShort_t d, Char_t r,
 343                                     Double_t v, TH2D*  h)
 344 {
 345   //
 346   // Set the acceptance correction @f$ a_{r,v}(\eta)@f$.
 347   // Note, that the object takes ownership of the passed pointer.
 348   //
 349   // Parameters:
 350   //    d    Detector number (1-3)
 351   //    r    Ring identifier (I or O)
 352   //    v    Primary interaction point @f$z@f$ coordinate
 353   //    h    @f$ a_{r,v}(\eta)@f$
 354   //
 355   // Return:
 356   //    true if operation succeeded
 357   //
 358   Int_t b = FindVertexBin(v);
 359   if (b <= 0 || b > fVertexAxis.GetNbins()) {
 360     AliWarning(Form("Vertex %+8.4f out of range [%+8.4f,%+8.4f]",
 361                     v, fVertexAxis.GetXmin(), fVertexAxis.GetXmax()));
 362     return false;
 363   }
 364   return SetCorrection(d, r, UShort_t(b), h);
 365 }
 366
 367 //____________________________________________________________________
 368 void
 369 AliFMDCorrAcceptance::FillCache() const
 370 {
 371   if (fCache) return;
 372
 373   fCache = new TObjArray;
 374   fCache->SetOwner(kTRUE);
 375   fCache->SetName("cache");
 376
 377   Int_t nV = fVertexAxis.GetNbins();
 378   for (UShort_t v = 1; v <= nV; v++) {
 379     for(UShort_t d = 1; d <= 3;d++) {
 380       UShort_t nR = (d == 1 ? 1 : 2);
 381       for (UShort_t q = 0; q < nR; q++) {
 382         Char_t   r  = (q == 0 ? 'I' : 'O');
 383
 384         TObjArray* a = GetOrMakeRingArray(*fCache, d, r);
 385
 386         TH2D* corr = GetCorrection(d, r, v);
 387         if (!corr) continue;
 388
 389         Int_t nY = corr->GetNbinsY();
 390         TH1D* h  = corr->ProjectionX("tmp", nY+1, nY+1, "");
 391         h->SetName(Form("FMD%d%c_vtxbin%03d", d, r, v));
 392         h->SetTitle(Form("#phi acceptance correction for FMD%d%c "
 393                            "in vertex bin %d [%+8.4f,%+8.4f]",
 394                            d, r, v, fVertexAxis.GetBinLowEdge(v),
 395                            fVertexAxis.GetBinUpEdge(v)));
 396         h->SetXTitle("#eta");
 397         h->SetYTitle("N_{strips}/N_{strips,OK}");
 398         h->SetFillStyle(3001);
 399         h->SetDirectory(0);
 400         h->SetStats(0);
 401         a->AddAtAndExpand(h,v-1);
 402
 403         if (fHasOverflow) continue;
 404
 405         // Construct the overflow bin from
 406         Int_t nX = corr->GetNbinsX();
 407         for (Int_t eta = 1; eta <= nX; eta++) {
 408           Double_t sum = 0;
 409           for (Int_t phi = 1; phi <= nY; phi++)
 410             sum += corr->GetBinContent(eta, phi);
 411           if (nY <= 0) continue;
 412           h->SetBinContent(eta, nY/sum);
 413         } // for eta
 414       } // for q
 415     } // for d
 416   } // for v
 417 }
 418 //____________________________________________________________________
 419 void
 420 AliFMDCorrAcceptance::Browse(TBrowser* b)
 421 {
 422   //
 423   // Browse this object in the browser
 424   //
 425   // Parameters:
 426   //    b
 427   //
 428   b->Add(&fRingArray);
 429   if (fCache) b->Add(fCache);
 430   b->Add(&fVertexAxis);
 431 }
 432 //____________________________________________________________________
 433 void
 434 AliFMDCorrAcceptance::Print(Option_t* option) const
 435 {
 436   //
 437   // Print this object
 438   //
 439   // Parameters:
 440   //    option
 441   //
 442   std::cout << "Acceptance correction due to dead channels" << std::endl;
 443   fRingArray.Print(option);
 444   fVertexAxis.Print(option);
 445 }
 446 //____________________________________________________________________
 447 void
 448 AliFMDCorrAcceptance::ls(Option_t* option) const
 449 {
 450   //
 451   // Print this object
 452   //
 453   // Parameters:
 454   //    option
 455   //
 456   TObject::ls(option);
 457   gROOT->IncreaseDirLevel();
 458   fVertexAxis.ls(option);
 459   fRingArray.ls(option);
 460   gROOT->DecreaseDirLevel();
 461 }
 462
 463 #if 0
 464 namespace {
 465   void ClearCanvas(TVirtualPad* c)
 466   {
 467     c->SetLeftMargin(.1);
 468     c->SetRightMargin(.05);
 469     c->SetBottomMargin(.1);
 470     c->SetTopMargin(.05);
 471     c->Clear();
 472   }
 473 }
 474
 475 //____________________________________________________________________
 476 void
 477 AliFMDCorrAcceptance::SaveAs(const Char_t* filename, Option_t* option) const
 478 {
 479   //
 480   // Override to allow saving to a PDF
 481   //
 482   TString fileName(filename);
 483   if (!fileName.EndsWith(".pdf")) {
 484     TObject::SaveAs(fileName, option);
 485     return;
 486   }
 487
 488   TVirtualPad* c = new TCanvas(filename, GetTitle(), 800/TMath::Sqrt(2), 800);
 489   c->SetFillColor(0);
 490   c->SetBorderSize(0);
 491   c->SetBorderMode(0);
 492   c->Print(Form("%s[", filename));
 493
 494   //__________________________________________________________________
 495   // Create a title page
 496   TLatex* ll = new TLatex(.5,.8, filename);
 497   ll->SetTextAlign(22);
 498   ll->SetTextSize(0.03);
 499   ll->SetNDC();
 500   ll->Draw();
 501
 502   TLatex* l = new TLatex(.5,.8, filename);
 503   l->SetNDC();
 504   l->SetTextSize(0.03);
 505   l->SetTextFont(132);
 506   l->SetTextAlign(12);
 507   l->DrawLatex(0.2, 0.70, "Acceptance due to dead channels");
 508   l->SetTextAlign(22);
 509   l->DrawLatex(0.5, 0.60, "c_{v,r}(#eta,#phi)=#frac{"
 510                "#sum active strips#in(#eta,#phi)}{"
 511                "#sum strips#in(#eta,#phi)}");
 512
 513   c->Print(filename, "Title:Title page");
 514
 515   //__________________________________________________________________
 516   // Draw all corrections
 517   const TAxis& vtxAxis = GetVertexAxis();
 518   Int_t        nVtx    = vtxAxis.GetNbins();
 519
 520   // --- Loop over detectors -----------------------------------------
 521   for (UShort_t d = 1; d <= 3; d++) {
 522     UShort_t     nQ = (d == 1 ? 1 : 2);
 523     for (UShort_t q = 0; q < nQ; q++) {
 524       Char_t r = (q == 0 ? 'I' : 'O');
 525
 526       ClearCanvas(c);
 527       c->Divide(2, (nVtx+1)/2);
 528       for (UShort_t v=1; v <= nVtx; v++) {
 529         TVirtualPad* p = c->cd(v);
 530         p->SetFillColor(kWhite);
 531
 532         TH2* h2 = GetCorrection(d, r, v);
 533         if (!h2) {
 534           Warning("DrawCorrAcc", "No correction for r=%c, v=%d", r, v);
 535           continue;
 536         }
 537         h2->Draw(option);
 538       }
 539       c->Print(filename, Form("Title:FMD%d%c", d, r));
 540     }
 541   }
 542   if (HasOverflow()){
 543     const_cast<AliFMDCorrAcceptance*>(this)->Draw(Form("%s phi", option));
 544     c->Print(filename, "Title:Phi Acceptance");
 545   }
 546   const_cast<AliFMDCorrAcceptance*>(this)->Draw(option);
 547   c->Print(filename, "Title:Summary");
 548   c->Print(Form("%s]", filename));
 549 }
 550 //____________________________________________________________________
 551 void
 552 AliFMDCorrAcceptance::Draw(Option_t* option)
 553 {
 554   //
 555   // Draw this object
 556   //
 557   // Parameters:
 558   //   option
 559   //
 560   TString opt(option);
 561   opt.ToLower();
 562   Bool_t over = opt.Contains("phi");
 563   opt.ReplaceAll("phi", "");
 564
 565   TVirtualPad* c = gPad;
 566   if (!c) c = new TCanvas(GetName(), GetTitle());
 567
 568   const TAxis& vtxAxis = fVertexAxis;
 569   Int_t        nVtx    = vtxAxis.GetNbins();
 570   Int_t        ipad    = 0;
 571   c->SetLeftMargin(.1);
 572   c->SetRightMargin(.05);
 573   c->SetBottomMargin(.1);
 574   c->SetTopMargin(.05);
 575   c->Clear();
 576   c->Divide((nVtx+2)/3, 3, 0, 0);
 577
 578   // Draw all corrections
 579   for (UShort_t v = 1; v <= nVtx; v++) {
 580     ipad++;
 581     if (ipad == 1 || ipad == 12) ipad++;
 582
 583     TVirtualPad* p = c->cd(ipad);
 584     p->SetFillColor(kWhite);
 585
 586     THStack* stack = new THStack(Form("vtx%02d", v),
 587                                  Form("%+5.1f<v_{z}<%+5.1f",
 588                                       vtxAxis.GetBinLowEdge(v),
 589                                       vtxAxis.GetBinUpEdge(v)));
 590     for (UShort_t d = 1; d <= 3; d++) {
 591       UShort_t     nQ = (d == 1 ? 1 : 2);
 592       for (UShort_t q = 0; q < nQ; q++) {
 593         Char_t r = (q == 0 ? 'I' : 'O');
 594
 595         if (over) {
 596           TH1* hp = GetPhiAcceptance(d, r, v);
 597           if (!hp) {
 598             Error("", "No phi acceptance at v=%d", v-1);
 599             continue;
 600           }
 601           hp->SetDirectory(0);
 602           hp->SetMarkerColor(AliForwardUtil::RingColor(d, r));
 603           hp->SetLineColor(AliForwardUtil::RingColor(d, r));
 604           hp->SetFillColor(AliForwardUtil::RingColor(d, r));
 605           hp->SetFillStyle(3001);
 606           // Info("", "Adding phi acceptance plot %d", int(hp->GetEntries()));
 607           stack->Add(hp);
 608           continue;
 609         }
 610
 611         TH2* h1 = GetCorrection(d, r, v);
 612         if (!h1) {
 613           Warning("Draw", "No correction for r=%c, v=%d", r, v);
 614           continue;
 615         }
 616         Int_t nY = h1->GetNbinsY();
 617         TH1* hh = h1->ProjectionX(Form("FMD%d%c", d, r), 1, nY);
 618         hh->Scale(1. / nY);
 619         hh->SetDirectory(0);
 620         hh->SetMarkerColor(AliForwardUtil::RingColor(d, r));
 621         hh->SetLineColor(AliForwardUtil::RingColor(d, r));
 622         hh->SetFillColor(AliForwardUtil::RingColor(d, r));
 623         hh->SetFillStyle(3004);
 624
 625         stack->Add(hh);
 626       }
 627     }
 628     stack->SetMaximum(1.2);
 629     stack->Draw(Form("nostack %s", opt.Data()));
 630   }
 631 }
 632 #endif
 633
 634 //____________________________________________________________________
 635 //
 636 // EOF
 637 //