[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / corrs / MultCutDrawer.C

#ifndef __CINT__
# include "SummaryDrawer.C"
# include "AliFMDCorrAcceptance.h"
# include "AliFMDCorrSecondaryMap.h"
# include "AliFMDCorrELossFit.h"
# include "AliFMDMultCuts.h"
# include "AliForwardUtil.h"
# include "AliForwardCorrectionManager.h"
# include "AliLog.h"
# include <TString.h>
# include <TError.h>
#else
class SummaryDrawer;
class TObject;
class AliFMDMultCuts;
class THStack;
class TH1;
#include <TString.h>
#endif

struct MultCutDrawer : public SummaryDrawer
{
  UShort_t fMinQuality;
  TList    fCuts;
  TList    fStacks;
  Bool_t   fMC;
  //__________________________________________________________________
  /** 
   * Constructor
   */
  MultCutDrawer() 
    : fMinQuality(8),
      fCuts(),
      fStacks(),
      fMC(false)
  {
    // Rough equvilance: 
    // 
    //   Cut name  |           Parameter values   
    //  -----------+----------------------------------------
    //    mpv      | 0.85     0.7      0.4     0.15
    //    xi       | 1        2.5      4.5     6.8
    //    sig      | .5       1        2       2.9
    //    prob     | 1e-1     2.5e-2   5e-4    2.5e-6
    //  -----------+----------------------------------------
    //   Cut name  |           Mean values   
    //  -----------+----------------------------------------
    //   mpv       | 0.43    0.36      0.20    0.08
    //   xi        | 0.49    0.36      0.23    0.09
    //   sig       | 0.44    0.37      0.22    0.10
    //   prob      | 0.43    0.35      0.21    0.09
    // 
    fCuts.Add(new TNamed("mpv",  "0.85  0.7    0.4  0.15"));
    fCuts.Add(new TNamed("xi",   "1     2.5    4.5  6.8"));
    fCuts.Add(new TNamed("sig",  ".5    1      2    2.9"));
    fCuts.Add(new TNamed("prob", "1e-1  2.5e-2 5e-4 2.5e-6"));
    // fCuts.Add(new TNamed("prob", "1e-2 1e-3 1e-5 1e-7"));
  }
  //__________________________________________________________________
  /** 
   * Copy constructor 
   * 
   * @param o Object to copy from 
   */
  MultCutDrawer(const MultCutDrawer& o) 
    : fMinQuality(o.fMinQuality),
      fCuts(),
      fStacks(),
      fMC(false)
  {}
  //__________________________________________________________________
  /** 
   * Assignment operator 
   * 
   * @param o Obejct to assign from
   * 
   * @return Reference to this 
   */
  MultCutDrawer& operator=(const MultCutDrawer& o) 
  {
    if (&o == this) return *this;
    fMinQuality = o.fMinQuality;
    fCuts.AddAll(&o.fCuts);
    fStacks.AddAll(&o.fStacks);
    return *this;
  }
  //__________________________________________________________________
  /** 
   * Destructor
   */
  ~MultCutDrawer()
  {
    CloseCanvas();
  }
  //__________________________________________________________________
  /** 
   * Run the class
   * 
   * @param runNo   Run number (or 999 for don't care)
   * @param sys     System (or 0 for don't care)
   * @param sNN     Collision energy in GeV (or 0 for don't care)
   * @param field   L3 Field in kG (or 999 for don't care)
   * @param mc      True of MC
   * @param local   Possible local database 
   */  
  void Run(ULong_t       runNo=999, 
	   UShort_t      sys=0, 
	   UShort_t      sNN=0, 
	   Short_t       field=999, 
	   Bool_t        mc=false, 
	   const Char_t* local=0)
  {
    Bool_t sat = false;
    if (!Init(runNo, sys, sNN, field, mc, sat, local)) return;

    Double_t savX = fParVal->GetX();
    Double_t savY = fParVal->GetY();
    fParVal->SetX(.4);
    fParVal->SetY(.4);
    // fPause = true;

    TIter    iCut(&fCuts);
    TObject* pCut = 0;
    while ((pCut = iCut())) { 
      TString                 method(pCut->GetName());
      TString                 sP(pCut->GetTitle());
      TObjArray*              aP = sP.Tokenize(" ");
      TIter                   iP(aP);
      TObjString*             pP = 0;
      TString                 tM;
      fBody->SetBottomMargin(0.20);
      fBody->SetLeftMargin(0.06);
      fBody->Divide(1, aP->GetEntries(), 0, 0);
      Int_t iPad = 1;
      while ((pP = static_cast<TObjString*>(iP()))) {
	THStack* all   = AllStack(iPad-1);
	Double_t p     = pP->String().Atof();
	Double_t vP[]  = { p, p, p, p, p };
	THStack* stack = CutStack(method, vP, all);
	if (tM.IsNull()) tM = stack->GetTitle();
	// Kill title on all but first sub-panel
	stack->SetTitle("");
	DrawInPad(fBody, iPad, stack, "nostack p");
	stack->GetYaxis()->SetTitleSize(0.12);
	stack->GetYaxis()->SetTitleOffset(0.2);
	stack->GetYaxis()->SetLabelSize(0.07);
	if (iPad == 1) stack->GetYaxis()->SetTitle(tM);
	stack->GetXaxis()->SetTitle("#eta");
	stack->GetXaxis()->SetTitleSize(0.12);
	stack->GetXaxis()->SetTitleOffset(0.6);
	stack->GetXaxis()->SetTitleColor(kBlack);
	stack->GetXaxis()->SetLabelSize(0.07);
	stack->GetXaxis()->Dump();


	if (iPad == 1) {
	  Color_t    col   = kBlack;
	  Double_t   hLtx  = 0.07;
	  Double_t   yLtx  = 7*(hLtx+.005)+0.01;
	  TLatex*    nLtx  = new TLatex(-0.75, yLtx, "Ring");
	  TLatex*    pLtx  = new TLatex(-0.7,  yLtx, "Param");
	  TLatex*    vLtx  = new TLatex(+0, yLtx, "Mean#pmVar_{min}^{max}");
	  nLtx->SetTextAlign(31);pLtx->SetTextAlign(11);
	  nLtx->SetTextSize(hLtx);
	  pLtx->SetTextSize(hLtx);
	  vLtx->SetTextSize(hLtx);
	  nLtx->SetTextColor(col);
	  pLtx->SetTextColor(col);
	  vLtx->SetTextColor(col);
	  nLtx->Draw();
	  pLtx->Draw();
	  vLtx->Draw();
	}
	// if (iPad == 1) { 
	//   fBody->cd(1);
	//   DrawRingLegend(0.4, 0.4, 0.7, 0.9);
	// }
	iPad++;
      }
      PrintCanvas(Form("%s   X={%s}", tM.Data(), sP.Data()));
    }

    Int_t nAll = fStacks.GetEntries();
    fBody->SetBottomMargin(0.20);
    fBody->SetLeftMargin(0.06);
    fBody->Divide(1, nAll, 0, 0);
    for (Int_t iAll = 0; iAll < nAll; iAll++) {
      THStack* all = AllStack(iAll);
      DrawInPad(fBody, iAll+1, all, "nostack hist p");
      all->GetYaxis()->SetTitleSize(0.12);
      all->GetYaxis()->SetTitleOffset(0.2);
      all->GetYaxis()->SetLabelSize(0.07);
      if (iAll == 0) all->GetYaxis()->SetTitle("c");
      all->GetXaxis()->SetTitle("#eta");
      all->GetXaxis()->SetTitleSize(0.12);
      all->GetXaxis()->SetTitleOffset(0.6);
      all->GetXaxis()->SetTitleColor(kBlack);
      all->GetXaxis()->SetLabelSize(0.07);
      
      TVirtualPad* p = fBody->GetPad(iAll+1);
      p->cd();
      Double_t yT = 1-p->GetTopMargin();
      if      (iAll == 0) DrawRingLegend(p, kNorth|kCenter); 
      else if (iAll == 1) DrawMethodLegend(0.35, 0.4, 0.55,yT);
      
      Double_t y1 = ((iAll + 2 >= nAll) ? yT - .3 : p->GetBottomMargin());
      Double_t y2 = ((iAll + 2 >= nAll) ? yT      : 0.3);
      DrawValueLegend(all, 0.2, y1, 0.9, y2);
       
    }
    PrintCanvas("Comparisons");

    fParVal->SetX(savX);
    fParVal->SetY(savY);
    
    CloseCanvas();
  }
  //__________________________________________________________________
  /** 
   * Draw legend of methods 
   * 
   * @param x1    Left x
   * @param y1    Lower y
   * @param x2    Right x
   * @param y2    Upper y
   */
  void DrawMethodLegend(Double_t x1, Double_t y1, 
			Double_t x2, Double_t y2) 
  {
    TLegend* l = new TLegend(x1,y1,x2,y2);
    l->SetBorderSize(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);

    TIter    iCut(&fCuts);
    TObject* pCut = 0;
    while ((pCut = iCut())) { 
      TString                 method(pCut->GetName());
      AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(method);
      TString                 title(AliFMDMultCuts::Method2String(m,true));
      Style_t                 style = CutStyle(m);
      TLegendEntry*           e = l->AddEntry("dummy", title, "p");
      e->SetMarkerStyle(style);
      e->SetMarkerColor(kBlack);
    }
    l->Draw();
  }
  //__________________________________________________________________
  /** 
   * Draw a value legend 
   * 
   * @param stack Stack to take values from 
   * @param x1    Left x
   * @param y1    Lower y
   * @param x2    Right x
   * @param y2    Upper y
   */
  void DrawValueLegend(THStack* stack, 
		       Double_t x1, Double_t y1, 
		       Double_t x2, Double_t y2) 
  {
    TLegend* l = new TLegend(x1,y1,x2,y2);
    l->SetBorderSize(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);

    TString seen;
    TIter   iHist(stack->GetHists());
    TH1*    pHist = 0;
    Int_t   nHist = 0;
    while ((pHist = static_cast<TH1*>(iHist()))) {
      TString name(pHist->GetName());
      if (seen.Contains(name)) continue; 
      seen.Append(Form(" %s", name.Data()));
      nHist++;
      
      AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(name);
      Style_t                 s = CutStyle(m);
      TLegendEntry*           e = l->AddEntry("dummy", pHist->GetTitle(), "p");
      e->SetMarkerStyle(s);
      e->SetMarkerColor(kBlack);
    }
    if (nHist < 5) l->SetNColumns(nHist);
    else           l->SetNColumns(nHist/2);

    l->Draw();
  }
  //__________________________________________________________________
  /** 
   * Get stack at @a i.  If the stack doesn't exist, make it
   * 
   * @param i Location (0-based)
   * 
   * @return Stack 
   */
  THStack* AllStack(Int_t i) 
  {
    TObject* o = fStacks.At(i);
    if (o) return static_cast<THStack*>(o);
    THStack* s = new THStack(Form("all%02d", i), "");
    fStacks.AddAt(s, i);
    return s;
  }
  //__________________________________________________________________
  /** 
   * Get the marker styoe associated with a cut
   * 
   * @param m Cut identifier 
   * 
   * @return Marker style 
   */
  static Style_t CutStyle(UShort_t m) 
  {
    switch (m) { 
    case AliFMDMultCuts::kFixed:            return kFullStar;
    case AliFMDMultCuts::kMPVFraction:      return kOpenCircle;
    case AliFMDMultCuts::kFitRange:         return 33; // Diamond
    case AliFMDMultCuts::kLandauWidth:      return 34; // Cross
    case AliFMDMultCuts::kLandauSigmaWidth: return kOpenSquare;
    case AliFMDMultCuts::kProbability:      return kFullTriangleDown;
    }
    return kFullDotMedium;
  }
  //__________________________________________________________________
  /** 
   * Update statistics
   * 
   * @param y     Current value 
   * @param cnt   Current count
   * @param mean  Current mean
   * @param var   Current variance 
   */
  static void Statistics(Double_t  y,
			 Int_t&    cnt,
			 Double_t& mean, 
			 Double_t& var) 
  {
    cnt        += 1;
    mean       += (y - mean) / cnt;
    var        += (cnt > 1 ? (TMath::Power(y-mean,2)/(cnt-1)-var/cnt) : 0);
  }
  //__________________________________________________________________
  /** 
   * Calculate statistics for a histogram
   * 
   * @param h     Histogram
   * @param mean  On return, the mean of y
   * @param var   On return, the variance in y
   * @param min   On return, the least y
   * @param max   On return, the largest y
   * @param rCnt  In/out: Current count
   * @param rMean In/out: Current mean 
   * @param rVar  In/out: Current variance 
   */
  static void HistStatistics(const TH1* h, 
			     Double_t& mean, 
			     Double_t& var, 
			     Double_t& min, 
			     Double_t& max,
			     Int_t&    rCnt,
			     Double_t& rMean,
			     Double_t& rVar)
  {
    mean      = 0;
    var       = 0;
    min       = +100000;
    max       = -100000;
    Int_t cnt = 0;
    for (Int_t i = 1; i <= h->GetNbinsX(); i++) { 
      Double_t y = h->GetBinContent(i);
      if (TMath::Abs(y) <= 1e-9) continue;
      min        =  TMath::Min(min, y);
      max        =  TMath::Max(max, y);
      Statistics(y, cnt, mean, var);
      Statistics(y, rCnt, rMean, rVar);
    }
    // Info("", "Stats for %s:  mean=%f +/- %f [%f,%f]",
    //      h->GetTitle(), mean, var, min, max);
  }
  //__________________________________________________________________
  /** 
   * Create a stack from cuts
   * 
   * @param method Method to use 
   * @param param  Parameters 
   * @param all    Stack to add for this set of parameters
   * 
   * @return Newly created stack 
   */    
  THStack* CutStack(const TString& method, Double_t* param, THStack* all)
  {
    AliFMDMultCuts::EMethod m = AliFMDMultCuts::String2Method(method);
    Info("CutStack", "Method %s -> %d", method.Data(), m);
    AliFMDMultCuts* cut = new AliFMDMultCuts(m, 
					     param[0],
					     param[1], 
					     param[2],
					     param[3],
					     param[4]);
    // cut->Print();
    
    TH2* hist = new TH2D("cut", cut->GetMethodString(true),
			 200, -4, 6, 5, 0.5, 5.5);
    hist->GetYaxis()->SetBinLabel(1, "FMD1i");
    hist->GetYaxis()->SetBinLabel(2, "FMD2i");
    hist->GetYaxis()->SetBinLabel(3, "FMD2o");
    hist->GetYaxis()->SetBinLabel(4, "FMD3o");
    hist->GetYaxis()->SetBinLabel(5, "FMD3i");
    
    // Info("DrawMultCuts", "Filling histogram");
    cut->FillHistogram(hist);
    // Info("DrawMultCuts", "Done filling");

    Style_t  style = CutStyle(m);
    THStack* stack = new THStack(hist, "x");
    TList*   hists = stack->GetHists();
    Double_t rMin  = +1000000;
    Double_t rMax  = -1000000;
    Double_t rAvg  = 0;
    Double_t rVar  = 0;
    Int_t    rCnt  = 0;
    printf(" %6s %7.4f | ", method.Data(), param[0]);
    TH1*     first = 0;
    for (Int_t i = 1; i <= 5; i++) { 
      TH1*     h = static_cast<TH1*>(hists->At(i-1));
      TString  n(hist->GetYaxis()->GetBinLabel(i)); 
      TString  nn(n); nn.Remove(0,3);
      UShort_t det = nn.Atoi();
      Char_t   rng = nn[1];
      Color_t  col = AliForwardUtil::RingColor(det, rng);
      if (!first) first = h;
      h->SetName(method);
      h->SetTitle(Form("%f", param[i-1]));
      h->SetYTitle(cut->GetMethodString(true));
      h->SetXTitle("#eta");
      h->SetMarkerColor(col);
      h->SetFillColor(col);
      h->SetLineColor(col);
      h->SetMarkerStyle(style);
      h->SetFillStyle(0);
      Double_t avg, var, min, max;
      
      HistStatistics(h, avg, var, min, max, rCnt, rAvg, rVar);
      rMin = TMath::Min(min, rMin);
      rMax = TMath::Max(max, rMax);
      all->Add(h);
      Double_t   hLtx  = 0.07;
      Double_t   yLtx  = i*(hLtx+.005)+0.01;
      TObjArray* lines = new TObjArray(3);
      TLatex*    nLtx  = new TLatex(-0.75, yLtx, n);
      TLatex*    pLtx  = new TLatex(-0.7,  yLtx, Form("X=%g", param[i-1]));
      TLatex*    vLtx  = new TLatex(+0, yLtx, 
				    Form("%5.3f#pm%6.4f_{%6.4f}^{%6.4f}",
					 avg, var, max-avg, avg-min));
      nLtx->SetTextAlign(31);pLtx->SetTextAlign(11);
      nLtx->SetTextSize(hLtx);pLtx->SetTextSize(hLtx),vLtx->SetTextSize(hLtx);
      nLtx->SetTextColor(col);pLtx->SetTextColor(col);vLtx->SetTextColor(col);
      lines->Add(nLtx);lines->Add(pLtx);lines->Add(vLtx);
      h->GetListOfFunctions()->Add(lines);
      printf("%5.3f+/-%6.4f ", avg, var);
    }
    TLatex* rLtx = new TLatex(6, fMC ? 0.65 : 0.55, 
			      Form("All: %5.3f#pm%6.4f_{%6.4f}^{%6.4f}",
				   rAvg, rVar, rMin, rMax));
    rLtx->SetTextSize(0.05);
    rLtx->SetTextAlign(31);
    first->GetListOfFunctions()->Add(rLtx);
    Printf("-> %5.3f+/-%6.4f", rAvg, rVar);
    stack->SetTitle(cut->GetMethodString(true)); // hist->GetTitle());
    stack->SetMinimum(0); // 0.98*min);
    stack->SetMaximum(fMC ? 0.7 : 0.6); // 1.02*max);
    all->SetMinimum(0);
    all->SetMaximum(fMC ? 0.7 : 0.6);

    delete hist;
    return stack;
  }
  //__________________________________________________________________
  /** 
   * Initialize 
   * 
   * @param runNo   Run number (or 999 for don't care)
   * @param sys     System (or 0 for don't care)
   * @param sNN     Collision energy in GeV (or 0 for don't care)
   * @param field   L3 Field in kG (or 999 for don't care)
   * @param mc      True of MC
   * @param sat     True for including satellite collisions 
   * @param local   Possible local database 
   * 
   * @return true on sucess 
   */
  Bool_t Init(ULong_t       runNo=999, 
	      UShort_t      sys=0, 
	      UShort_t      sNN=0, 
	      Short_t       field=999, 
	      Bool_t        mc=false, 
	      Bool_t        sat=false,
	      const Char_t* local=0)
  {
    fMC = mc;
    AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
    mgr.SetDebug(true);
    UShort_t flags = AliForwardCorrectionManager::kELossFits;
  
    if (local && local[0] != '\0') mgr.SetELossFitsPath(local);
  
    if (!mgr.Init(runNo, sys, sNN, field, mc, false, flags, true)) {
      Error("DrawMultCuts", "Failed to initialize for flags=0x%02x, "
	    "run=%lu, sys=%hu, sNN=%hu, field=%hd, mc=%s, sat=%s",
	    flags, runNo, sys, sNN, field, mc ? "true" : "false", "false");
      return false;
    }
    const AliFMDCorrELossFit* cFit = mgr.GetELossFit();
    AliFMDCorrELossFit*       fit  = const_cast<AliFMDCorrELossFit*>(cFit);
    fit->CacheBins(8);

    CreateCanvas("multCuts.pdf", true);

    fBody->cd();
    
    Double_t y = .85;
    TLatex* title = new TLatex(.5, y, "#Delta Cuts");
    title->SetTextAlign(23);
    title->SetTextFont(42);
    title->SetTextSize(0.1);
    title->Draw();
    
    y -= 0.11;
    DrawParameter(y, "Run #", Form("%lu", runNo));
    DrawParameter(y, "System", AliForwardUtil::CollisionSystemString(sys));
    DrawParameter(y, "#sqrt{s_{NN}}", 
		  AliForwardUtil::CenterOfMassEnergyString(sNN));
    DrawParameter(y, "L3 field", AliForwardUtil::MagneticFieldString(field));
    DrawParameter(y, "Simulation", Form("%s", mc ? "yes" : "no"));
    DrawParameter(y, "Satellite", Form("%s", sat ? "yes" : "no"));
    PrintCanvas("Delta cuts");

    return true;
  }
};
Commit	Line	Data
f37328ed	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);
a19faec0	120	// fPause = true;
a19faec0	121
f37328ed	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;
a19faec0	130	TString tM;
	131	fBody->SetBottomMargin(0.20);
	132	fBody->SetLeftMargin(0.06);
f37328ed	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);
a19faec0	140	if (tM.IsNull()) tM = stack->GetTitle();
f37328ed	141	// Kill title on all but first sub-panel
a19faec0	142	stack->SetTitle("");
f37328ed	143	DrawInPad(fBody, iPad, stack, "nostack p");
a19faec0	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	}
f37328ed	174	// if (iPad == 1) {
	175	// fBody->cd(1);
	176	// DrawRingLegend(0.4, 0.4, 0.7, 0.9);
	177	// }
	178	iPad++;
	179	}
a19faec0	180	PrintCanvas(Form("%s X={%s}", tM.Data(), sP.Data()));
f37328ed	181	}
a19faec0	182
f37328ed	183	Int_t nAll = fStacks.GetEntries();
a19faec0	184	fBody->SetBottomMargin(0.20);
a19faec0	185	fBody->SetLeftMargin(0.06);
f37328ed	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");
a19faec0	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
f37328ed	200	TVirtualPad* p = fBody->GetPad(iAll+1);
	201	p->cd();
	202	Double_t yT = 1-p->GetTopMargin();
a19faec0	203	if (iAll == 0) DrawRingLegend(p, kNorth\|kCenter);
f37328ed	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]));
a19faec0	434	h->SetYTitle(cut->GetMethodString(true));
a19faec0	435	h->SetXTitle("#eta");
f37328ed	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);
a19faec0	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);
f37328ed	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);
a19faec0	469	stack->SetTitle(cut->GetMethodString(true)); // hist->GetTitle());
f37328ed	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
a19faec0	517	CreateCanvas("multCuts.pdf", true);
f37328ed	518
	519	fBody->cd();
	520
a19faec0	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;
f37328ed	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"));
a19faec0	536	PrintCanvas("Delta cuts");
f37328ed	537
	538	return true;
	539	}
	540	};