PWGLF/FORWARD/analysis2/corrs/MultCutDrawer.C

   1 #ifndef __CINT__
   2 # include "SummaryDrawer.C"
   3 # include "AliFMDCorrAcceptance.h"
   4 # include "AliFMDCorrSecondaryMap.h"
   5 # include "AliFMDCorrELossFit.h"
   6 # include "AliFMDMultCuts.h"
   7 # include "AliForwardUtil.h"
   8 # include "AliForwardCorrectionManager.h"
   9 # include "AliLog.h"
  10 # include <TString.h>
  11 # include <TError.h>
  12 #else
  13 class SummaryDrawer;
  14 class TObject;
  15 class AliFMDMultCuts;
  16 class THStack;
  17 class TH1;
  18 #include <TString.h>
  19 #endif
  20
  21 struct MultCutDrawer : public SummaryDrawer
  22 {
  23   UShort_t fMinQuality;
  24   TList    fCuts;
  25   TList    fStacks;
  26   Bool_t   fMC;
  27   //__________________________________________________________________
  28   /**
  29    * Constructor
  30    */
  31   MultCutDrawer()
  32     : fMinQuality(8),
  33       fCuts(),
  34       fStacks(),
  35       fMC(false)
  36   {
  37     // Rough equvilance:
  38     //
  39     //   Cut name  |           Parameter values
  40     //  -----------+----------------------------------------
  41     //    mpv      | 0.85     0.7      0.4     0.15
  42     //    xi       | 1        2.5      4.5     6.8
  43     //    sig      | .5       1        2       2.9
  44     //    prob     | 1e-1     2.5e-2   5e-4    2.5e-6
  45     //  -----------+----------------------------------------
  46     //   Cut name  |           Mean values
  47     //  -----------+----------------------------------------
  48     //   mpv       | 0.43    0.36      0.20    0.08
  49     //   xi        | 0.49    0.36      0.23    0.09
  50     //   sig       | 0.44    0.37      0.22    0.10
  51     //   prob      | 0.43    0.35      0.21    0.09
  52     //
  53     fCuts.Add(new TNamed("mpv",  "0.85  0.7    0.4  0.15"));
  54     fCuts.Add(new TNamed("xi",   "1     2.5    4.5  6.8"));
  55     fCuts.Add(new TNamed("sig",  ".5    1      2    2.9"));
  56     fCuts.Add(new TNamed("prob", "1e-1  2.5e-2 5e-4 2.5e-6"));
  57     // fCuts.Add(new TNamed("prob", "1e-2 1e-3 1e-5 1e-7"));
  58   }
  59   //__________________________________________________________________
  60   /**
  61    * Copy constructor
  62    *
  63    * @param o Object to copy from
  64    */
  65   MultCutDrawer(const MultCutDrawer& o)
  66     : fMinQuality(o.fMinQuality),
  67       fCuts(),
  68       fStacks(),
  69       fMC(false)
  70   {}
  71   //__________________________________________________________________
  72   /**
  73    * Assignment operator
  74    *
  75    * @param o Obejct to assign from
  76    *
  77    * @return Reference to this
  78    */
  79   MultCutDrawer& operator=(const MultCutDrawer& o)
  80   {
  81     if (&o == this) return *this;
  82     fMinQuality = o.fMinQuality;
  83     fCuts.AddAll(&o.fCuts);
  84     fStacks.AddAll(&o.fStacks);
  85     return *this;
  86   }
  87   //__________________________________________________________________
  88   /**
  89    * Destructor
  90    */
  91   ~MultCutDrawer()
  92   {
  93     CloseCanvas();
  94   }
  95   //__________________________________________________________________
  96   /**
  97    * Run the class
  98    *
  99    * @param runNo   Run number (or 999 for don't care)
 100    * @param sys     System (or 0 for don't care)
 101    * @param sNN     Collision energy in GeV (or 0 for don't care)
 102    * @param field   L3 Field in kG (or 999 for don't care)
 103    * @param mc      True of MC
 104    * @param local   Possible local database
 105    */
 106   void Run(ULong_t       runNo=999,
 107            UShort_t      sys=0,
 108            UShort_t      sNN=0,
 109            Short_t       field=999,
 110            Bool_t        mc=false,
 111            const Char_t* local=0)
 112   {
 113     Bool_t sat = false;
 114     if (!Init(runNo, sys, sNN, field, mc, sat, local)) return;
 115
 116     Double_t savX = fParVal->GetX();
 117     Double_t savY = fParVal->GetY();
 118     fParVal->SetX(.4);
 119     fParVal->SetY(.4);
 120     // fPause = true;
 121
 122     TIter    iCut(&fCuts);
 123     TObject* pCut = 0;
 124     while ((pCut = iCut())) {
 125       TString                 method(pCut->GetName());
 126       TString                 sP(pCut->GetTitle());
 127       TObjArray*              aP = sP.Tokenize(" ");
 128       TIter                   iP(aP);
 129       TObjString*             pP = 0;
 130       TString                 tM;
 131       fBody->SetBottomMargin(0.20);
 132       fBody->SetLeftMargin(0.06);
 133       fBody->Divide(1, aP->GetEntries(), 0, 0);
 134       Int_t iPad = 1;
 135       while ((pP = static_cast<TObjString*>(iP()))) {
 136         THStack* all   = AllStack(iPad-1);
 137         Double_t p     = pP->String().Atof();
 138         Double_t vP[]  = { p, p, p, p, p };
 139         THStack* stack = CutStack(method, vP, all);
 140         if (tM.IsNull()) tM = stack->GetTitle();
 141         // Kill title on all but first sub-panel
 142         stack->SetTitle("");
 143         DrawInPad(fBody, iPad, stack, "nostack p");
 144         stack->GetYaxis()->SetTitleSize(0.12);
 145         stack->GetYaxis()->SetTitleOffset(0.2);
 146         stack->GetYaxis()->SetLabelSize(0.07);
 147         if (iPad == 1) stack->GetYaxis()->SetTitle(tM);
 148         stack->GetXaxis()->SetTitle("#eta");
 149         stack->GetXaxis()->SetTitleSize(0.12);
 150         stack->GetXaxis()->SetTitleOffset(0.6);
 151         stack->GetXaxis()->SetTitleColor(kBlack);
 152         stack->GetXaxis()->SetLabelSize(0.07);
 153         stack->GetXaxis()->Dump();
 154
 155
 156         if (iPad == 1) {
 157           Color_t    col   = kBlack;
 158           Double_t   hLtx  = 0.07;
 159           Double_t   yLtx  = 7*(hLtx+.005)+0.01;
 160           TLatex*    nLtx  = new TLatex(-0.75, yLtx, "Ring");
 161           TLatex*    pLtx  = new TLatex(-0.7,  yLtx, "Param");
 162           TLatex*    vLtx  = new TLatex(+0, yLtx, "Mean#pmVar_{min}^{max}");
 163           nLtx->SetTextAlign(31);pLtx->SetTextAlign(11);
 164           nLtx->SetTextSize(hLtx);
 165           pLtx->SetTextSize(hLtx);
 166           vLtx->SetTextSize(hLtx);
 167           nLtx->SetTextColor(col);
 168           pLtx->SetTextColor(col);
 169           vLtx->SetTextColor(col);
 170           nLtx->Draw();
 171           pLtx->Draw();
 172           vLtx->Draw();
 173         }
 174         // if (iPad == 1) {
 175         //   fBody->cd(1);
 176         //   DrawRingLegend(0.4, 0.4, 0.7, 0.9);
 177         // }
 178         iPad++;
 179       }
 180       PrintCanvas(Form("%s   X={%s}", tM.Data(), sP.Data()));
 181     }
 182
 183     Int_t nAll = fStacks.GetEntries();
 184     fBody->SetBottomMargin(0.20);
 185     fBody->SetLeftMargin(0.06);
 186     fBody->Divide(1, nAll, 0, 0);
 187     for (Int_t iAll = 0; iAll < nAll; iAll++) {
 188       THStack* all = AllStack(iAll);
 189       DrawInPad(fBody, iAll+1, all, "nostack hist p");
 190       all->GetYaxis()->SetTitleSize(0.12);
 191       all->GetYaxis()->SetTitleOffset(0.2);
 192       all->GetYaxis()->SetLabelSize(0.07);
 193       if (iAll == 0) all->GetYaxis()->SetTitle("c");
 194       all->GetXaxis()->SetTitle("#eta");
 195       all->GetXaxis()->SetTitleSize(0.12);
 196       all->GetXaxis()->SetTitleOffset(0.6);
 197       all->GetXaxis()->SetTitleColor(kBlack);
 198       all->GetXaxis()->SetLabelSize(0.07);
 199
 200       TVirtualPad* p = fBody->GetPad(iAll+1);
 201       p->cd();
 202       Double_t yT = 1-p->GetTopMargin();
 203       if      (iAll == 0) DrawRingLegend(p, kNorth|kCenter);
 204       else if (iAll == 1) DrawMethodLegend(0.35, 0.4, 0.55,yT);
 205
 206       Double_t y1 = ((iAll + 2 >= nAll) ? yT - .3 : p->GetBottomMargin());
 207       Double_t y2 = ((iAll + 2 >= nAll) ? yT      : 0.3);
 208       DrawValueLegend(all, 0.2, y1, 0.9, y2);
 209
 210     }
 211     PrintCanvas("Comparisons");
 212
 213     fParVal->SetX(savX);
 214     fParVal->SetY(savY);
 215
 216     CloseCanvas();
 217   }
 218   //__________________________________________________________________
 219   /**
 220    * Draw legend of methods
 221    *
 222    * @param x1    Left x
 223    * @param y1    Lower y
 224    * @param x2    Right x
 225    * @param y2    Upper y
 226    */
 227   void DrawMethodLegend(Double_t x1, Double_t y1,
 228                         Double_t x2, Double_t y2)
 229   {
 230     TLegend* l = new TLegend(x1,y1,x2,y2);
 231     l->SetBorderSize(0);
 232     l->SetFillColor(0);
 233     l->SetFillStyle(0);
 234
 235     TIter    iCut(&fCuts);
 236     TObject* pCut = 0;
 237     while ((pCut = iCut())) {
 238       TString                 method(pCut->GetName());
 239       AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(method);
 240       TString                 title(AliFMDMultCuts::Method2String(m,true));
 241       Style_t                 style = CutStyle(m);
 242       TLegendEntry*           e = l->AddEntry("dummy", title, "p");
 243       e->SetMarkerStyle(style);
 244       e->SetMarkerColor(kBlack);
 245     }
 246     l->Draw();
 247   }
 248   //__________________________________________________________________
 249   /**
 250    * Draw a value legend
 251    *
 252    * @param stack Stack to take values from
 253    * @param x1    Left x
 254    * @param y1    Lower y
 255    * @param x2    Right x
 256    * @param y2    Upper y
 257    */
 258   void DrawValueLegend(THStack* stack,
 259                        Double_t x1, Double_t y1,
 260                        Double_t x2, Double_t y2)
 261   {
 262     TLegend* l = new TLegend(x1,y1,x2,y2);
 263     l->SetBorderSize(0);
 264     l->SetFillColor(0);
 265     l->SetFillStyle(0);
 266
 267     TString seen;
 268     TIter   iHist(stack->GetHists());
 269     TH1*    pHist = 0;
 270     Int_t   nHist = 0;
 271     while ((pHist = static_cast<TH1*>(iHist()))) {
 272       TString name(pHist->GetName());
 273       if (seen.Contains(name)) continue;
 274       seen.Append(Form(" %s", name.Data()));
 275       nHist++;
 276
 277       AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(name);
 278       Style_t                 s = CutStyle(m);
 279       TLegendEntry*           e = l->AddEntry("dummy", pHist->GetTitle(), "p");
 280       e->SetMarkerStyle(s);
 281       e->SetMarkerColor(kBlack);
 282     }
 283     if (nHist < 5) l->SetNColumns(nHist);
 284     else           l->SetNColumns(nHist/2);
 285
 286     l->Draw();
 287   }
 288   //__________________________________________________________________
 289   /**
 290    * Get stack at @a i.  If the stack doesn't exist, make it
 291    *
 292    * @param i Location (0-based)
 293    *
 294    * @return Stack
 295    */
 296   THStack* AllStack(Int_t i)
 297   {
 298     TObject* o = fStacks.At(i);
 299     if (o) return static_cast<THStack*>(o);
 300     THStack* s = new THStack(Form("all%02d", i), "");
 301     fStacks.AddAt(s, i);
 302     return s;
 303   }
 304   //__________________________________________________________________
 305   /**
 306    * Get the marker styoe associated with a cut
 307    *
 308    * @param m Cut identifier
 309    *
 310    * @return Marker style
 311    */
 312   static Style_t CutStyle(UShort_t m)
 313   {
 314     switch (m) {
 315     case AliFMDMultCuts::kFixed:            return kFullStar;
 316     case AliFMDMultCuts::kMPVFraction:      return kOpenCircle;
 317     case AliFMDMultCuts::kFitRange:         return 33; // Diamond
 318     case AliFMDMultCuts::kLandauWidth:      return 34; // Cross
 319     case AliFMDMultCuts::kLandauSigmaWidth: return kOpenSquare;
 320     case AliFMDMultCuts::kProbability:      return kFullTriangleDown;
 321     }
 322     return kFullDotMedium;
 323   }
 324   //__________________________________________________________________
 325   /**
 326    * Update statistics
 327    *
 328    * @param y     Current value
 329    * @param cnt   Current count
 330    * @param mean  Current mean
 331    * @param var   Current variance
 332    */
 333   static void Statistics(Double_t  y,
 334                          Int_t&    cnt,
 335                          Double_t& mean,
 336                          Double_t& var)
 337   {
 338     cnt        += 1;
 339     mean       += (y - mean) / cnt;
 340     var        += (cnt > 1 ? (TMath::Power(y-mean,2)/(cnt-1)-var/cnt) : 0);
 341   }
 342   //__________________________________________________________________
 343   /**
 344    * Calculate statistics for a histogram
 345    *
 346    * @param h     Histogram
 347    * @param mean  On return, the mean of y
 348    * @param var   On return, the variance in y
 349    * @param min   On return, the least y
 350    * @param max   On return, the largest y
 351    * @param rCnt  In/out: Current count
 352    * @param rMean In/out: Current mean
 353    * @param rVar  In/out: Current variance
 354    */
 355   static void HistStatistics(const TH1* h,
 356                              Double_t& mean,
 357                              Double_t& var,
 358                              Double_t& min,
 359                              Double_t& max,
 360                              Int_t&    rCnt,
 361                              Double_t& rMean,
 362                              Double_t& rVar)
 363   {
 364     mean      = 0;
 365     var       = 0;
 366     min       = +100000;
 367     max       = -100000;
 368     Int_t cnt = 0;
 369     for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 370       Double_t y = h->GetBinContent(i);
 371       if (TMath::Abs(y) <= 1e-9) continue;
 372       min        =  TMath::Min(min, y);
 373       max        =  TMath::Max(max, y);
 374       Statistics(y, cnt, mean, var);
 375       Statistics(y, rCnt, rMean, rVar);
 376     }
 377     // Info("", "Stats for %s:  mean=%f +/- %f [%f,%f]",
 378     //      h->GetTitle(), mean, var, min, max);
 379   }
 380   //__________________________________________________________________
 381   /**
 382    * Create a stack from cuts
 383    *
 384    * @param method Method to use
 385    * @param param  Parameters
 386    * @param all    Stack to add for this set of parameters
 387    *
 388    * @return Newly created stack
 389    */
 390   THStack* CutStack(const TString& method, Double_t* param, THStack* all)
 391   {
 392     AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(method);
 393     Info("CutStack", "Method %s -> %d", method.Data(), m);
 394     AliFMDMultCuts* cut = new AliFMDMultCuts(m,
 395                                              param[0],
 396                                              param[1],
 397                                              param[2],
 398                                              param[3],
 399                                              param[4]);
 400     // cut->Print();
 401
 402     TH2* hist = new TH2D("cut", cut->GetMethodString(true),
 403                          200, -4, 6, 5, 0.5, 5.5);
 404     hist->GetYaxis()->SetBinLabel(1, "FMD1i");
 405     hist->GetYaxis()->SetBinLabel(2, "FMD2i");
 406     hist->GetYaxis()->SetBinLabel(3, "FMD2o");
 407     hist->GetYaxis()->SetBinLabel(4, "FMD3o");
 408     hist->GetYaxis()->SetBinLabel(5, "FMD3i");
 409
 410     // Info("DrawMultCuts", "Filling histogram");
 411     cut->FillHistogram(hist);
 412     // Info("DrawMultCuts", "Done filling");
 413
 414     Style_t  style = CutStyle(m);
 415     THStack* stack = new THStack(hist, "x");
 416     TList*   hists = stack->GetHists();
 417     Double_t rMin  = +1000000;
 418     Double_t rMax  = -1000000;
 419     Double_t rAvg  = 0;
 420     Double_t rVar  = 0;
 421     Int_t    rCnt  = 0;
 422     printf(" %6s %7.4f | ", method.Data(), param[0]);
 423     TH1*     first = 0;
 424     for (Int_t i = 1; i <= 5; i++) {
 425       TH1*     h = static_cast<TH1*>(hists->At(i-1));
 426       TString  n(hist->GetYaxis()->GetBinLabel(i));
 427       TString  nn(n); nn.Remove(0,3);
 428       UShort_t det = nn.Atoi();
 429       Char_t   rng = nn[1];
 430       Color_t  col = AliForwardUtil::RingColor(det, rng);
 431       if (!first) first = h;
 432       h->SetName(method);
 433       h->SetTitle(Form("%f", param[i-1]));
 434       h->SetYTitle(cut->GetMethodString(true));
 435       h->SetXTitle("#eta");
 436       h->SetMarkerColor(col);
 437       h->SetFillColor(col);
 438       h->SetLineColor(col);
 439       h->SetMarkerStyle(style);
 440       h->SetFillStyle(0);
 441       Double_t avg, var, min, max;
 442
 443       HistStatistics(h, avg, var, min, max, rCnt, rAvg, rVar);
 444       rMin = TMath::Min(min, rMin);
 445       rMax = TMath::Max(max, rMax);
 446       all->Add(h);
 447       Double_t   hLtx  = 0.07;
 448       Double_t   yLtx  = i*(hLtx+.005)+0.01;
 449       TObjArray* lines = new TObjArray(3);
 450       TLatex*    nLtx  = new TLatex(-0.75, yLtx, n);
 451       TLatex*    pLtx  = new TLatex(-0.7,  yLtx, Form("X=%g", param[i-1]));
 452       TLatex*    vLtx  = new TLatex(+0, yLtx,
 453                                     Form("%5.3f#pm%6.4f_{%6.4f}^{%6.4f}",
 454                                          avg, var, max-avg, avg-min));
 455       nLtx->SetTextAlign(31);pLtx->SetTextAlign(11);
 456       nLtx->SetTextSize(hLtx);pLtx->SetTextSize(hLtx),vLtx->SetTextSize(hLtx);
 457       nLtx->SetTextColor(col);pLtx->SetTextColor(col);vLtx->SetTextColor(col);
 458       lines->Add(nLtx);lines->Add(pLtx);lines->Add(vLtx);
 459       h->GetListOfFunctions()->Add(lines);
 460       printf("%5.3f+/-%6.4f ", avg, var);
 461     }
 462     TLatex* rLtx = new TLatex(6, fMC ? 0.65 : 0.55,
 463                               Form("All: %5.3f#pm%6.4f_{%6.4f}^{%6.4f}",
 464                                    rAvg, rVar, rMin, rMax));
 465     rLtx->SetTextSize(0.05);
 466     rLtx->SetTextAlign(31);
 467     first->GetListOfFunctions()->Add(rLtx);
 468     Printf("-> %5.3f+/-%6.4f", rAvg, rVar);
 469     stack->SetTitle(cut->GetMethodString(true)); // hist->GetTitle());
 470     stack->SetMinimum(0); // 0.98*min);
 471     stack->SetMaximum(fMC ? 0.7 : 0.6); // 1.02*max);
 472     all->SetMinimum(0);
 473     all->SetMaximum(fMC ? 0.7 : 0.6);
 474
 475     delete hist;
 476     return stack;
 477   }
 478   //__________________________________________________________________
 479   /**
 480    * Initialize
 481    *
 482    * @param runNo   Run number (or 999 for don't care)
 483    * @param sys     System (or 0 for don't care)
 484    * @param sNN     Collision energy in GeV (or 0 for don't care)
 485    * @param field   L3 Field in kG (or 999 for don't care)
 486    * @param mc      True of MC
 487    * @param sat     True for including satellite collisions
 488    * @param local   Possible local database
 489    *
 490    * @return true on sucess
 491    */
 492   Bool_t Init(ULong_t       runNo=999,
 493               UShort_t      sys=0,
 494               UShort_t      sNN=0,
 495               Short_t       field=999,
 496               Bool_t        mc=false,
 497               Bool_t        sat=false,
 498               const Char_t* local=0)
 499   {
 500     fMC = mc;
 501     AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
 502     mgr.SetDebug(true);
 503     UShort_t flags = AliForwardCorrectionManager::kELossFits;
 504
 505     if (local && local[0] != '\0') mgr.SetELossFitsPath(local);
 506
 507     if (!mgr.Init(runNo, sys, sNN, field, mc, false, flags, true)) {
 508       Error("DrawMultCuts", "Failed to initialize for flags=0x%02x, "
 509             "run=%lu, sys=%hu, sNN=%hu, field=%hd, mc=%s, sat=%s",
 510             flags, runNo, sys, sNN, field, mc ? "true" : "false", "false");
 511       return false;
 512     }
 513     const AliFMDCorrELossFit* cFit = mgr.GetELossFit();
 514     AliFMDCorrELossFit*       fit  = const_cast<AliFMDCorrELossFit*>(cFit);
 515     fit->CacheBins(8);
 516
 517     CreateCanvas("multCuts.pdf", true);
 518
 519     fBody->cd();
 520
 521     Double_t y = .85;
 522     TLatex* title = new TLatex(.5, y, "#Delta Cuts");
 523     title->SetTextAlign(23);
 524     title->SetTextFont(42);
 525     title->SetTextSize(0.1);
 526     title->Draw();
 527
 528     y -= 0.11;
 529     DrawParameter(y, "Run #", Form("%lu", runNo));
 530     DrawParameter(y, "System", AliForwardUtil::CollisionSystemString(sys));
 531     DrawParameter(y, "#sqrt{s_{NN}}",
 532                   AliForwardUtil::CenterOfMassEnergyString(sNN));
 533     DrawParameter(y, "L3 field", AliForwardUtil::MagneticFieldString(field));
 534     DrawParameter(y, "Simulation", Form("%s", mc ? "yes" : "no"));
 535     DrawParameter(y, "Satellite", Form("%s", sat ? "yes" : "no"));
 536     PrintCanvas("Delta cuts");
 537
 538     return true;
 539   }
 540 };