PWGLF/FORWARD/analysis/scripts/ViewEDCor.C

   1 //____________________________________________________________________
   2 AliFMDAnaCalibEnergyDistribution*
   3 ViewEDCor(const char* filename)
   4 {
   5   // gSystem->Load("libANALYSIS");
   6   // gSystem->Load("libANALYSISalice");
   7   // gSystem->Load("libPWG0base");
   8   // gSystem->Load("libPWGLFforward");
   9
  10   AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();
  11   param->SetEnergy(900.F);
  12   param->SetCollisionSystem("pp");
  13   param->SetMagField(+5.F);
  14   param->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
  15
  16   TFile* file = TFile::Open(filename, "READ");
  17   if (!file) {
  18     Error("ViewEDCor", "Cannot open file %s", filename);
  19     return 0;
  20   }
  21
  22   AliFMDAnaCalibEnergyDistribution* bg =
  23     static_cast<AliFMDAnaCalibEnergyDistribution*>(file->Get("energydistributions"));
  24
  25   return bg;
  26 }
  27
  28
  29 //____________________________________________________________________
  30 struct RingHistos
  31 {
  32   TH1D* MakeHisto(const char* name, const char* title,
  33                   const TAxis* etaAxis)
  34   {
  35     TH1D* ret = new TH1D(Form(name, fD, fR),
  36                          Form(title, fD, fR),
  37                          etaAxis->GetNbins(),
  38                          etaAxis->GetXmin(),
  39                          etaAxis->GetXmax());
  40     ret->SetDirectory(0);
  41     ret->SetLineColor(fC);
  42     ret->SetMarkerColor(fC);
  43     ret->SetMarkerStyle(20+2*fD+(fR == 'I' || fR == 'i' ? 0 : 1));
  44     ret->SetFillColor(fC);
  45     ret->SetFillStyle(3001);
  46     ret->SetXTitle("#eta");
  47     ret->SetYTitle(title);
  48     ret->SetStats(0);
  49     return ret;
  50   }
  51
  52
  53   RingHistos(UShort_t d, Char_t r, const TAxis* etaAxis, Int_t c)
  54     : fD(d), fR(r), fC(c)
  55   {
  56     fHists.AddAt(MakeHisto("c%d%c",     "Constant FMD%d%c",      etaAxis),0);
  57     fHists.AddAt(MakeHisto("mpv%d%c",   "MPV FMD%d%c",           etaAxis),1);
  58     fHists.AddAt(MakeHisto("w%d%c",     "Width FMD%d%c",         etaAxis),2);
  59     fHists.AddAt(MakeHisto("alpha%d%c", "#alpha FMD%d%c",        etaAxis),3);
  60     fHists.AddAt(MakeHisto("beta%d%c",  "#beta FMD%d%c",         etaAxis),4);
  61     fHists.AddAt(MakeHisto("chi%d%c",   "#chi^{2}/#nu FMD%d%c",  etaAxis),5);
  62     fHists.AddAt(MakeHisto("acc%d%c",   "Acceptance in FMD%d%c", etaAxis),6);
  63   }
  64   void Fill(Int_t iEta, Int_t p, Double_t v, Double_t e=0)
  65   {
  66     TH1D* h = Get(p);
  67     if (!h) return;
  68     h->SetBinContent(iEta, v);
  69     h->SetBinError(iEta, e);
  70   }
  71   TH1D* Get(Int_t p) const
  72   {
  73     return static_cast<TH1D*>(fHists.At(p));
  74   }
  75   UShort_t fD;
  76   Char_t   fR;
  77   Int_t    fC;
  78   TObjArray fHists;
  79 };
  80
  81 //____________________________________________________________________
  82 struct All
  83 {
  84   All(const TAxis* etaAxis)
  85   {
  86     fFMD1i = new RingHistos(1,'i', etaAxis, kRed+1);
  87     fFMD2i = new RingHistos(2,'i', etaAxis, kGreen+1);
  88     fFMD2o = new RingHistos(2,'o', etaAxis, kGreen+7);
  89     fFMD3i = new RingHistos(3,'i', etaAxis, kBlue+1);
  90     fFMD3o = new RingHistos(3,'o', etaAxis, kBlue+7);
  91   }
  92   RingHistos* Get(UShort_t d, Char_t r) const
  93   {
  94     switch (d) {
  95     case 1:    return fFMD1i;
  96     case 2:    return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
  97     case 3:    return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
  98     }
  99     return 0;
 100   }
 101   RingHistos* fFMD1i;
 102   RingHistos* fFMD2i;
 103   RingHistos* fFMD2o;
 104   RingHistos* fFMD3i;
 105   RingHistos* fFMD3o;
 106 };
 107
 108 //____________________________________________________________________
 109 PlotEDCor(const char* filename)
 110 {
 111   AliFMDAnaCalibEnergyDistribution* ed = ViewEDCor(filename);
 112
 113   AliFMDAnaParameters* param = AliFMDAnaParameters::Instance();
 114
 115   TAxis* etaAxis = param->GetBackgroundCorrection(1,'I',4)->GetXaxis();
 116   All* all = new All(etaAxis);
 117
 118   for (Int_t iEta = 1; iEta < etaAxis->GetNbins(); iEta++) {
 119     Float_t eta = etaAxis->GetBinCenter(iEta);
 120     Info("PlotEDCor", "eta bin %3d, eta=%f", iEta, eta);
 121
 122     for (UShort_t d=1; d <= 3; d++) {
 123       UShort_t nq = d == 1 ? 1 : 2;
 124       for (UShort_t q = 0; q < nq; q++) {
 125         Char_t r = q == 0 ? 'I' : 'O';
 126         RingHistos* hists = all->Get(d, r);
 127         if (!hists) continue;
 128
 129         TH1F* hed = ed->GetEnergyDistribution(d, r, eta);
 130         if (!hed) continue;
 131
 132         if (hed->GetEntries() > 0)
 133           hists->Fill(iEta, 6, 1);
 134
 135         TF1* fed = hed->GetFunction("FMDfitFunc");
 136         if (!fed) continue;
 137
 138         Int_t npar = fed->GetNpar();
 139         for (Int_t i = 0; i < npar; i++) {
 140           Double_t par = fed->GetParameter(i);
 141           Double_t err = fed->GetParError(i);
 142           hists->Fill(iEta, i, par, err);
 143         }
 144         Double_t chi2 = fed->GetChisquare();
 145         Double_t ndf  = fed->GetNDF();
 146         if (ndf != 0)
 147           hists->Fill(iEta, 5, chi2/ndf);
 148       }
 149     }
 150   }
 151
 152   TCanvas* c = new TCanvas("c", "c", 800, 800);
 153   c->SetFillColor(0);
 154   c->SetBorderMode(0);
 155   c->SetBorderSize(0);
 156   c->Divide(1,7,0,0);
 157
 158   THStack* stack = 0;
 159
 160   const char* yTitles[] = { "Constant", "#Delta_{p}", "w",
 161                             "#alpha", "#beta", "#chi^2/NDF",
 162                             "Acceptance" };
 163   for (Int_t i = 0; i < 7; i++) {
 164     TVirtualPad* p = c->cd(i+1);
 165     p->SetFillColor(0);
 166     p->SetFillStyle(0);
 167     THStack* stack = new THStack;
 168     stack->Add(all->Get(1,'i')->Get(i));
 169     stack->Add(all->Get(2,'i')->Get(i));
 170     stack->Add(all->Get(2,'o')->Get(i));
 171     stack->Add(all->Get(3,'i')->Get(i));
 172     stack->Add(all->Get(3,'o')->Get(i));
 173     stack->Draw("nostack");
 174     stack->GetHistogram()->SetYTitle(yTitles[i]);
 175     stack->GetHistogram()->SetXTitle("#eta");
 176     TAxis* yaxis = stack->GetHistogram()->GetYaxis();
 177     yaxis->SetTitleSize(0.15);
 178     yaxis->SetLabelSize(0.08);
 179     yaxis->SetTitleOffset(0.35);
 180     yaxis->SetNdivisions(10);
 181     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 182     xaxis->SetTitleSize(0.15);
 183     xaxis->SetLabelSize(0.08);
 184     xaxis->SetTitleOffset(0.35);
 185     xaxis->SetNdivisions(320);
 186     p->cd();
 187   }
 188   c->cd();
 189   c->Print("energyDist.png");
 190 }
 191 //____________________________________________________________________
 192 //
 193 // EOF
 194 //
 195