PWGLF/FORWARD/analysis2/corrs/Outliers.C

   1 void
   2 TestSPD(const TString& which, Double_t nVar=2)
   3 {
   4   TFile* file = TFile::Open("forward.root", "READ");
   5   if (!file) return;
   6
   7   Bool_t spd = which.EqualTo("spd", TString::kIgnoreCase);
   8
   9   TList* l = 0;
  10   if (spd) l = static_cast<TList*>(file->Get("CentralSums"));
  11   else     l = static_cast<TList*>(file->Get("ForwardSums"));
  12   if (!l) {
  13     Warning("", "%sSums not found", spd ? "Central" : "Forward");
  14     return;
  15   }
  16
  17   TList* ei = static_cast<TList*>(l->FindObject("fmdEventInspector"));
  18   if (!l) {
  19     Warning("", "fmdEventInspector not found");
  20     return;
  21   }
  22
  23   TObject* run = ei->FindObject("runNo");
  24   if (!run)
  25     Warning("", "No run number found");
  26   ULong_t runNo = run ? run->GetUniqueID() : 0;
  27
  28   TH2* h = 0;
  29   if (spd) h = static_cast<TH2*>(l->FindObject("nClusterVsnTracklet"));
  30   else {
  31     TList* den = static_cast<TList*>(l->FindObject("fmdDensityCalculator"));
  32     if (!den) {
  33       Error("", "fmdDensityCalculator not found");
  34       return;
  35     }
  36     TList* rng = static_cast<TList*>(den->FindObject(which));
  37     if (!rng) {
  38       Error("", "%s not found", which.Data());
  39       return;
  40     }
  41     h = static_cast<TH2*>(rng->FindObject("elossVsPoisson"));
  42   }
  43   if (!h) {
  44     Warning("", "%s not found", spd ? nClusterVsnTracklet : "elossVsPoisson");
  45     return;
  46   }
  47
  48   gStyle->SetOptFit(1111);
  49   gStyle->SetOptStat(0);
  50   TCanvas* c = new TCanvas("c", Form("Run %u", runNo));
  51   c->Divide(2,2);
  52
  53   TVirtualPad* p = c->cd(1);
  54   if (spd) {
  55     p->SetLogx();
  56     p->SetLogy();
  57   }
  58   p->SetLogz();
  59   h->Draw("colz");
  60
  61   TObjArray* fits = new TObjArray;
  62   h->FitSlicesY(0, 1, -1, 0, "QN", fits);
  63
  64   TF1* mean = new TF1("mean", "pol1");
  65   TF1* var  = new TF1("var", "pol1");
  66   // mean->FixParameter(0, 0);
  67   // var->FixParameter(0, 0);
  68   for (Int_t i = 0; i < 3; i++) {
  69     p = c->cd(2+i);
  70     if (spd) {
  71       p->SetLogx();
  72       p->SetLogy();
  73     }
  74     TH1* hh = static_cast<TH1*>(fits->At(i));
  75     hh->Draw();
  76
  77     if (i == 0) continue;
  78
  79     hh->Fit((i == 1? mean : var), "+Q");
  80
  81   }
  82
  83   TGraphErrors* g1 = new TGraphErrors(h->GetNbinsX());
  84   g1->SetFillColor(kBlue-10);
  85   g1->SetFillStyle(3001);
  86   g1->SetLineStyle(1);
  87   TGraph* u1 = new TGraph(h->GetNbinsX());
  88   TGraph* l1 = new TGraph(h->GetNbinsX());
  89   u1->SetLineColor(kBlue+1);
  90   l1->SetLineColor(kBlue+1);
  91   u1->SetName("u1");
  92   l1->SetName("l1");
  93   TGraphErrors* g2 = new TGraphErrors(h->GetNbinsX());
  94   g2->SetFillColor(kRed-10);
  95   g2->SetFillStyle(3001);
  96   g2->SetLineStyle(2);
  97   TGraph* u2 = new TGraph(h->GetNbinsX());
  98   TGraph* l2 = new TGraph(h->GetNbinsX());
  99   u2->SetLineColor(kRed+1);
 100   l2->SetLineColor(kRed+1);
 101   u2->SetName("u2");
 102   l2->SetName("l2");
 103   for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 104     Double_t x  = hh->GetXaxis()->GetBinCenter(i);
 105     Double_t y  = mean->Eval(x);
 106     Double_t e  = var->Eval(y);
 107     Double_t e1 = nVar * e;
 108     if (spd) e1 *= TMath::Log10(e);
 109     // Printf("%10e -> %10e +/- %10e", x, y, ee);
 110     g1->SetPoint(i-1, x, y);
 111     g1->SetPointError(i-1, 0, e1);
 112     u1->SetPoint(i-1, x, y+e1);
 113     l1->SetPoint(i-1, x, y-e1);
 114     // Printf("%3d: %f -> %f +/- %f", i, x, y, ee);
 115
 116     Double_t e2 = nVar*0.05*x;
 117     g2->SetPoint(i-1, x, x);
 118     g2->SetPointError(i-1, 0, e2);
 119     u2->SetPoint(i-1, x, x+e2);
 120     l2->SetPoint(i-1, x, x-e2);
 121   }
 122
 123   p = c->cd(1);
 124   c->Clear();
 125   c->cd();
 126   c->SetLogz();
 127   h->Draw("colz");
 128   g1->Draw("3 same");
 129   u1->Draw("l same");
 130   l1->Draw("l same");
 131   g2->Draw("3 same");
 132   u2->Draw("l same");
 133   l2->Draw("l same");
 134
 135   Double_t ly = 0.9;
 136   Double_t dy = 0.06;
 137   TLatex* ltx = new TLatex(0.15, ly, Form("#LTy#GT = %f + %f x",
 138                                            mean->GetParameter(0),
 139                                            mean->GetParameter(1)));
 140   ltx->SetNDC();
 141   ltx->SetTextSize(dy);
 142   ltx->SetTextAlign(13);
 143   ltx->Draw();
 144
 145   ly -= dy + 0.01;
 146   ltx->DrawLatex(0.15, ly, Form("#sigma_{y} = %f + %f x",
 147                                 var->GetParameter(0),
 148                                 var->GetParameter(1)));
 149
 150   ly -= dy + 0.01;
 151   ltx->DrawLatex(0.15, ly, Form("#delta = %f #sigma %s",
 152                                 nVar, (spd ? "log_{10}(#sigma" : "")));
 153
 154
 155 }
 156
 157
 158