B2 analysis code
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / macros / DrawPt.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // Draw corrected pt
17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19 #include <Riostream.h>
20 #include <TROOT.h>
21 #include <TFile.h>
22 #include <TH1D.h>
23 #include <TString.h>
24 #include <TCanvas.h>
25 #include <TGraphErrors.h>
26 #include <TLegend.h>
27
28 #include <RooWorkspace.h>
29 #include <RooMsgService.h>
30 #include <RooPlot.h>
31 #include <RooRealVar.h>
32 #include <RooAbsPdf.h>
33 #include <RooAbsData.h>
34
35 #include "B2.h"
36
37 void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="AntiDeuteron", Bool_t m2pid=0, Int_t lowm2bin=9, Int_t him2bin=17)
38 {
39 //
40 // Draw corrected pt for debugging
41 //
42         Double_t xmin = 0;
43         Double_t xmax = 3.5;
44         
45         TFile* finput = new TFile(inputFile.Data());
46         if (finput->IsZombie()) exit(1);
47         
48         // m2 data fitted models
49         
50         if(m2pid)
51         {
52                 TH1D* hPidPt = (TH1D*)FindObj(finput, tag, Form("%s_PID_Pt",particle.Data()));
53                 Double_t binwidth = hPidPt->GetBinWidth(0);
54                 
55                 using namespace RooFit;
56                 
57                 // disable verbose in RooFit
58                 RooMsgService::instance().setSilentMode(1);
59                 RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
60                 
61                 TCanvas* c0 = new TCanvas(Form("%s.M2",particle.Data()), Form("M2 for %ss",particle.Data()));
62                 c0->Divide(3,3);
63                 
64                 TGraphErrors* grFitSigmaPt = new TGraphErrors();
65                 grFitSigmaPt->SetName(Form("%s_Fit_Sigma_Pt", particle.Data()));
66                 
67                 TGraphErrors* grFitM2Pt = new TGraphErrors();
68                 grFitM2Pt->SetName(Form("%s_Fit_M2_Pt", particle.Data()));
69                 
70                 for(Int_t i=lowm2bin, j=0; i<him2bin && i-lowm2bin < 9 ; ++i)
71                 {
72                         c0->cd(i-lowm2bin+1);
73                         //gPad->SetLogy(0);
74                         
75                         RooWorkspace* w= (RooWorkspace*)FindObj(finput, tag, Form("%s_M2_%02d",particle.Data(),i));
76                         
77                         RooPlot* m2frame = w->var("m2")->frame();
78                         
79                         w->data("data")->plotOn(m2frame);
80                         
81                         w->pdf("model")->plotOn(m2frame, LineWidth(2));
82                         w->pdf("model")->plotOn(m2frame, Components(*(w->pdf("Sd"))),LineWidth(1), LineColor(9));
83                         w->pdf("model")->plotOn(m2frame, Components(*(w->pdf("Bkg"))),LineWidth(1), LineColor(46),LineStyle(kDashed));
84                         
85                         Double_t ptMin = (i-1)*binwidth;
86                         Double_t ptMax = i*binwidth;
87                         Double_t pt = (ptMin+ptMax)/2.;
88                         
89                         m2frame->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", ptMin, ptMax));
90                         m2frame->SetMinimum(0.2);
91                 
92                         m2frame->Draw();
93                         
94                         grFitSigmaPt->SetPoint(j, pt, w->var("width")->getVal());
95                         grFitSigmaPt->SetPointError(j, 0, w->var("width")->getError());
96                         
97                         grFitM2Pt->SetPoint(j, pt, w->var("mean")->getVal());
98                         grFitM2Pt->SetPointError(j++, 0, w->var("mean")->getError());
99                 }
100                 
101                 // model parameters
102                 
103                 TCanvas* c1 = new TCanvas(Form("%s.M2.FitParameters",particle.Data()), Form("M2 model parameters for %ss",particle.Data()));
104                 
105                 c1->Divide(2,1);
106                 
107                 c1->cd(1);
108                 gPad->SetGridx();
109                 gPad->SetGridy();
110                 
111                 grFitSigmaPt->SetMarkerStyle(kFullCircle);
112                 grFitSigmaPt->SetMarkerColor(kBlue);
113                 grFitSigmaPt->SetLineColor(kBlue);
114                 grFitSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
115                 grFitSigmaPt->GetYaxis()->SetTitle("#sigma_{m^{2}}");
116                 grFitSigmaPt->Draw("ALP");
117                 
118                 c1->cd(2);
119                 gPad->SetGridx();
120                 gPad->SetGridy();
121                 
122                 grFitM2Pt->SetMarkerStyle(kFullCircle);
123                 grFitM2Pt->SetMarkerColor(kBlue);
124                 grFitM2Pt->SetLineColor(kBlue);
125                 grFitM2Pt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
126                 grFitM2Pt->GetYaxis()->SetTitle("m^{2} (GeV^{2}/c^{4})");
127                 grFitM2Pt->Draw("ALP");
128                 
129                 //delete grFitSigmaPt;
130                 //delete grFitM2Pt;
131         }
132         
133         // remaining corrections
134         
135         const Int_t kNum = 6;
136         const TString kCorr[kNum] = { "EffCor", "Unfolded", "SecCor", "FakeCor", "M2Corr", "PID"};
137         const TString kLabel[kNum] = { "Efficiency", "Unfolding", "Secondaries", "Fake tracks", "PID contamination", "Raw"};
138         const Int_t kColor[kNum]  = {kGreen-3, kGreen-2, kRed, kBlue, kOrange+1, kAzure};
139         const Int_t kMarker[kNum] = {kFullSquare, kFullCircle, kFullTriangleDown, kOpenCircle, kFullTriangleUp, kOpenSquare};
140         
141         TLegend* legend = new TLegend(0.5689655,0.6355932,0.8362069,0.8326271,0,"brNDC");
142         legend->SetTextSize(0.03);
143         legend->SetFillColor(0);
144         legend->SetBorderSize(0);
145         
146         TH1D* hPt[kNum];
147         
148         for(Int_t i=0; i<kNum; ++i)
149         {
150                 hPt[i] = (TH1D*)FindObj(finput, tag, particle + "_" + kCorr[i] + "_Pt");
151                 hPt[i]->SetLineColor(kColor[i]);
152                 hPt[i]->SetMarkerColor(kColor[i]);
153                 hPt[i]->SetMarkerStyle(kMarker[i]);
154                 legend->AddEntry(hPt[i], kLabel[i], "lp");
155         }
156         
157         TCanvas* c2 = new TCanvas(Form("%s.Pt",particle.Data()), Form("Pt for %s",particle.Data()));
158         c2->SetLogy();
159         
160         hPt[0]->SetTitle(particle.Data());
161         hPt[0]->SetAxisRange(xmin,xmax,"X");
162         hPt[0]->Draw("E");
163         for(Int_t i=1; i<kNum; ++i) hPt[i]->Draw("sameE");
164         
165         legend->Draw();
166 }