Excluding lowest pt point for K and p, ratios fit/data are now shown in 3 canvases
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / Fit / CombineSpectra.C
CommitLineData
5f2ca1e6 1#if !defined(__CINT__) || defined(__MAKECINT__)
4c0d7bc7 2#include <iostream>
3#include "TH1F.h"
4#include "TGraphErrors.h"
5#include "AliBWFunc.h"
6#include "AliBWTools.h"
7#include "TF1.h"
8#include "TFile.h"
9#include "TDatabasePDG.h"
10#include "TROOT.h"
11#include "TCanvas.h"
12#include "TFolder.h"
13#include "TStyle.h"
14#include "AliLatexTable.h"
15#include "TLegend.h"
16#include "TVirtualFitter.h"
17#include "TMath.h"
18#include "TH2F.h"
19#include "TSystem.h"
20#include "TLine.h"
21#include "TLatex.h"
22#include "TMath.h"
23
24#include "TASImage.h"
25#include "TPaveText.h"
4c0d7bc7 26#include "StarPPSpectra.C"
27#include "GetE735Ratios.C"
3472a634 28#include "TString.h"
3725266c 29#include "TObjString.h"
77233d47 30#include "TGraphAsymmErrors.h"
31
5f2ca1e6 32#endif
4c0d7bc7 33
3725266c 34
4c0d7bc7 35using namespace std;
36
37
38// A bunch of useful enums and constants
39enum {kPion=0,kKaon,kProton,kNPart};
ee08b77d 40enum {kTPC=0,kTOF,kITS,kITSTPC,kK0,kKinks,kCombTOFTPC,kCombAll,kNHist};// "k0" listed here as a kind of PID method...
4c0d7bc7 41const Int_t kNDet = kITS+2;
42enum {kPos=0,kNeg,kNCharge};
43enum {kPhojet=0,kPyTuneAtlasCSC, kPyTuneCMS6D6T, kPyTunePerugia0, kNTunes} ;
44enum {kFitLevi=0, kFitUA1, kFitPowerLaw,
45 kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0,
46 kNFit};
77233d47 47enum {kDoFits=0, kDoRatios, kDoSuperposition, kDoDrawWithModels, kDoCompareToStar, kDoDrawSyst, kDoCompareToAllCharged, kDoHelp};
058137a6 48enum {kStatErrors = 0, kSystErrors, kStatSystErrors}; // which errors do we put in the histo that we fit? stat,syst or stat+syst?
4c0d7bc7 49
50// flags, labels and names
51const char * partFlag[] = {"Pion", "Kaon", "Proton"};
d8e6677d 52const char * detFlag[] = {"TPC", "TOF", "ITS standalone", "ITS-TPC", "K0", "Kinks", "Combined TOF + TPC", "Combined TOF + TPC + ITS"};
4c0d7bc7 53const char * chargeFlag[] = {"Pos", "Neg"};
54const char * chargeLabel[] = {"Positive", "Negative"};
55const char * partLabel[kNPart][kNCharge] = {{"#pi^{+}", "#pi^{-}"},
56 // {"K^{+} (#times 2)", "K^{-} (#times 2)"},
57 {"K^{+}", "K^{-}"},
58 {"p" , "#bar{p}"}};
5f2ca1e6 59const char * partLatex[kNPart][kNCharge] = {{"$\\pi^{+}$", "$\\pi^{-}$"},
60 // {"K^{+} (#times 2)", "K^{-} (#times 2)"},
61 {"$K^{+}$", "$K^{-}$"},
62 {"$p$" , "$\\bar{p}$"}};
4c0d7bc7 63const char * mcTuneName[] = {"Phojet",
64 "Pythia - CSC 306",
65 "Pythia - D6T 109",
66 "Pythia - Perugia0 - 320", };
67const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
68
69// Style
70//const Int_t marker[] = {25,24,28,20,21}; // set for highlithining marek
d8e6677d 71//const Int_t marker[] = {20,24,25,28,21}; // standard set
72const Int_t marker[] = {24,25,28,27,21}; // No full symbols
4c0d7bc7 73const Int_t color [] = {1,2,4};
74
75const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
76const Int_t mcLineStyle[] = {1,2,3,4};
77
78
79// Template needed to combine different detectors
ee08b77d 80const Float_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2,2.4,2.6};
81Int_t nbinsTempl=34;
4c0d7bc7 82
83TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
84
85// Globals
86TH1F * hSpectra[kNHist][kNPart][kNCharge];
87TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
3725266c 88TH1F * hSystError[kNHist][kNPart][kNCharge];
4c0d7bc7 89Double_t mass[kNPart];
90
91// Functions:
92// Loading
93void LoadSpectra() ;
94void LoadMC() ;
95
96// Additional tasks (may be uncommented below)
3725266c 97void LoadLibs();
4c0d7bc7 98void DrawStar(Int_t icharge);
99void GetITSResiduals();
100void DrawWithModels() ;
101void DrawAllAndKaons();
102void DrawWithJacek();
103void DrawRatioToStar();
104void DrawRatios();
3725266c 105void DrawWithSyst();
5f2ca1e6 106void FitCombined();
3725266c 107void PrintCanvas(TCanvas* c,const TString formats) ;
4d4a268c 108void RebinHistosForRatios() ;
c4fe358b 109void Help();
4c0d7bc7 110
111// External stuff
5f2ca1e6 112void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance"){
113
55056a39 114 TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.72,0.89,0.89);
115 //TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.62,0.85,0.75);
5f2ca1e6 116 myPadLogo->SetFillColor(0);
117 myPadLogo->SetBorderMode(0);
118 myPadLogo->SetBorderSize(2);
119 myPadLogo->SetFrameBorderMode(0);
120 myPadLogo->SetLeftMargin(0.0);
121 myPadLogo->SetTopMargin(0.0);
122 myPadLogo->SetBottomMargin(0.0);
123 myPadLogo->SetRightMargin(0.0);
124 myPadLogo->SetFillStyle(0);
125 myPadLogo->Draw();
126 myPadLogo->cd();
127 TASImage *myAliceLogo = new TASImage("alice_logo.png");
128 myAliceLogo->Draw();
55056a39 129 c->cd();
130 TPaveText* t1=new TPaveText(0.418103, 0.837798, 0.656609, 0.888393,"NDC");
5f2ca1e6 131 t1->SetFillStyle(0);
132 t1->SetBorderSize(0);
133 t1->AddText(0.,0.,label);
134 t1->SetTextColor(kRed);
135 t1->SetTextFont(42);
55056a39 136 t1->SetTextSize(0.04);
5f2ca1e6 137 t1->Draw();
55056a39 138 TPaveText* t2=new TPaveText(0.418103, 0.80, 0.656609, 0.84,"NDC");
5f2ca1e6 139 t2->SetFillStyle(0);
140 t2->SetBorderSize(0);
55056a39 141 t2->AddText(0.,0.,"pp at #sqrt{s} = 900 GeV (2009 data)");
5f2ca1e6 142 t2->SetTextColor(kRed);
55056a39 143 t2->SetTextFont(42);
144 t2->SetTextSize(0.027);
5f2ca1e6 145 t2->Draw();
55056a39 146 TPaveText* t3=new TPaveText(0.418103, 0.76, 0.656609, 0.80,"NDC");
147 t3->SetFillStyle(0);
148 t3->SetBorderSize(0);
149 t3->AddText(0.,0.,"Statistical and systematic errors");
150 t3->SetTextColor(kRed);
151 t3->SetTextFont(42);
152 t3->SetTextSize(0.027);
153 t3->Draw();
154 // TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
155 // t2->SetFillStyle(0);
156 // t2->SetBorderSize(0);
157 // t2->SetTextColor(kRed);
158 // t2->SetTextFont(52);
159 // t2->AddText(0.,0.,today.Data());
160 // t2->Draw();
5f2ca1e6 161}
4c0d7bc7 162
163// Used to tag plots
164TDatime dt;
165TString today = "";
166
167
168
169// Switches
170Bool_t convertToMT = 0;
77233d47 171Bool_t sumCharge = 0;
78033eeb 172Int_t whatToFit = kStatErrors;
4c0d7bc7 173Bool_t scaleKaons = kFALSE;
0e6f7392 174Bool_t drawStar = kFALSE; // Overlay star when doing fits
4c0d7bc7 175Bool_t correctSecondaries = 1;
176Bool_t correctGeantFlukaXS = 1;
5f2ca1e6 177Int_t iCombInStudy = kCombAll; //kCombTOFTPC
c4fe358b 178Int_t fitFuncID = kFitLevi;
5f2ca1e6 179Bool_t showMC=kTRUE;
180Bool_t showE735=kTRUE;
77233d47 181Bool_t useSecCorrFromDCA=kTRUE;
2210b002 182Bool_t flagPreliminary=kFALSE; // Add "preliminary" flag and logo to the plots
183Bool_t doPrint = 0;
184TString printFormats = "C,eps,root"; // format in which canvases will be printed, if PrintCanvas is called (not all prints are based on printcanvas at the moment). This is a comma separated list.
22abbb27 185
4c0d7bc7 186
eea641de 187void CombineSpectra(Int_t analysisType=kDoHelp, Int_t locfitFuncID = kFitLevi) { //kDoSuperposition;//kDoDrawWithModels;// kDoFits; //kDoRatios;
4c0d7bc7 188
189 // This macro is used to combine the 900 GeV spectra from 2009
190
c4fe358b 191 fitFuncID=locfitFuncID;
192
4c0d7bc7 193 // Load stuff
3725266c 194 LoadLibs();
4c0d7bc7 195
eea641de 196 // Print Help and quit
197 if (analysisType == kDoHelp) {
198 Help();
199 return;
200 }
201
202
4c0d7bc7 203 // Set date
204 today = today + long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
205
206
207 // Set Masses
208 mass[kPion] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
209 mass[kKaon] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
210 mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
211
212 // Load histos
213 LoadSpectra();
214 LoadMC();
215
216 // Additional tasks
217 //DrawStar(icharge);
218 // GetITSResiduals();
5f2ca1e6 219 if(analysisType==kDoSuperposition) DrawAllAndKaons();
220 else if(analysisType==kDoDrawWithModels) DrawWithModels() ;
77233d47 221 else if(analysisType==kDoCompareToAllCharged) DrawWithJacek();
c4fe358b 222 else if(analysisType==kDoCompareToStar) DrawRatioToStar();
5f2ca1e6 223 else if(analysisType==kDoRatios) DrawRatios();
3725266c 224 else if(analysisType==kDoDrawSyst) DrawWithSyst();
5f2ca1e6 225 else if(analysisType==kDoFits) FitCombined();
b7018bc7 226 return;
5f2ca1e6 227}
4c0d7bc7 228
229
5f2ca1e6 230void FitCombined() {
eea641de 231 gStyle->SetOptTitle(0);
232 gStyle->SetOptStat(0);
233
4c0d7bc7 234 // Draw combined & Fit
235 AliBWFunc * fm = new AliBWFunc;
236 fm->SetVarType(AliBWFunc::kdNdpt);
237 if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
238
239 // table to print results
2f7f947e 240 AliLatexTable table(10,"c|ccccccccc");
4c0d7bc7 241 if (fitFuncID == kFitLevi) {
55056a39 242 table.InsertCustomRow("Part & Yield & Yield (FIT) & T Slope & n & $\\chi^2$/NDF & Min X & Extrap (%) & $\\langle p_{t} \\rangle (Fit)$ & $\\langle p_{t} \\rangle$ (FIT+DATA)\\\\");
4c0d7bc7 243 } else if (fitFuncID == kFitPowerLaw) {
55056a39 244 table.InsertCustomRow("Part & Yield & Norm & n & pt0 & $\\chi^2$/NDF & Min X & Frac Above & $\\langle p_{t} \\rangle (Fit)$ & $\\langle p_{t} \\rangle$ (FIT+DATA) \\\\");
4c0d7bc7 245 } else {
55056a39 246 table.InsertCustomRow("Part & Yield & Par0 & Par2 & Par1 & $\\chi^2$/NDF & Min X & Frac Above & $\\langle p_{t} \\rangle (Fit)$ & $\\langle p_{t} \\rangle$ (FIT+DATA) \\\\");
4c0d7bc7 247
248 }
2f7f947e 249 table.InsertHline();
49ad829f 250 AliLatexTable tempTable(3,"c|cc");
251 tempTable.InsertCustomRow("Part & Yield & $\\langle p_{t} \\rangle$ \\\\");
252 tempTable.InsertHline();
4c0d7bc7 253
2210b002 254 TH1F* hRatiosToFit[kNPart][kNCharge]; // Ratio data/fit function, stat error by default
255 TH1F* hRatiosToFitSyst[kNPart][kNCharge]; // Ratio data/fit, stat + syst
058137a6 256 // Fit all
257 Int_t chargeLoop = sumCharge ? 1 : 2;
57a944cd 258 Bool_t divideRatiosCanvas=kTRUE;
058137a6 259 for(Int_t icharge = 0; icharge < chargeLoop; icharge++){
5f2ca1e6 260
261 TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge],700,700);
262 c2->SetTickx();
263 c2->SetTicky();
264 c2->SetLeftMargin(0.14);
57a944cd 265 TCanvas * c2r = 0;
266 if(divideRatiosCanvas) c2r = new TCanvas(TString("cCombinedRatio")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombinedRatio")+chargeFlag[icharge],1200,500);
267 else c2r = new TCanvas(TString("cCombinedRatio")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombinedRatio")+chargeFlag[icharge],700,700);
5f2ca1e6 268 c2->cd();
269 gPad->SetLogy();
270 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
4c0d7bc7 271 hempty->SetXTitle("p_{t} (GeV/c)");
272 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
5f2ca1e6 273 hempty->GetYaxis()->SetTitleOffset(1.35);
274 hempty->GetXaxis()->SetTitleOffset(1.1);
4c0d7bc7 275 hempty->Draw();
5f2ca1e6 276 c2r->cd();
57a944cd 277 if(divideRatiosCanvas) {
278 c2r->Divide(3,1);
279 c2r->cd(1);
280 }
5f2ca1e6 281 gPad->SetGridy();
282 TH2F * hemptyR = new TH2F(TString("hemptyR")+long(icharge),"hemptyR",100,0.,2.9, 100, 0.5,1.5);
283 hemptyR->SetXTitle("p_{t} (GeV/c)");
284 hemptyR->SetYTitle("Data/Fit");
285 hemptyR->Draw();
22abbb27 286
57a944cd 287
288 TLegend * l = new TLegend(0.176724, 0.181548, 0.477011, 0.462798,chargeLabel[icharge]);
4c0d7bc7 289 l->SetFillColor(kWhite);
5f2ca1e6 290 l->SetTextSize(0.035);
55056a39 291
292
293
294 TPaveText* tf=new TPaveText(0.321839, 0.175595, 0.685345, 0.299107,"NDC");
5f2ca1e6 295 if(fitFuncID == kFitLevi){
296 tf->AddText("#frac{dN}{dp_{t}} #propto p_{t} #left(1+#frac{#sqrt{m^{2}+p_{t}^{2}} -m}{nT} #right)^{-n}");
297 // tf->SetNDC();
298 tf->SetTextFont(12);
55056a39 299 tf->SetTextSize(0.032);
5f2ca1e6 300 }
57a944cd 301
4c0d7bc7 302 for(Int_t ipart = 0; ipart < kNPart; ipart++){
5f2ca1e6 303 printf(" ----- Fit %s %s ------\n",partFlag[ipart],chargeFlag[icharge]);
4c0d7bc7 304 Float_t fitmin = 0;
305 Float_t fitmax = 3;
306
307 // Get functions
308 TF1 * func = 0;
3472a634 309 Int_t normPar = -1;
4c0d7bc7 310 if(fitFuncID == kFitLevi) {
3472a634 311 normPar = 0; // The levi is normalized by parameter 0
4c0d7bc7 312 if (ipart == kPion)
313 func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
314 if (ipart == kKaon)
315 func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
316 if (ipart == kProton)
5f2ca1e6 317 func = fm->GetLevi(mass[ipart], 0.15, 8.5, 0.09);
318 }
4c0d7bc7 319 else if(fitFuncID == kFitUA1) func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
320 else if(fitFuncID == kFitPowerLaw) {
22abbb27 321 // fitmin=0;
322 // fitmax=1;
4c0d7bc7 323 if (ipart == kPion)
324 func = fm->GetPowerLaw(1.0, 8.6, 7);
325 if (ipart == kKaon)
326 func = fm->GetPowerLaw(3.0, 12, 2.6);
327 if (ipart == kProton)
328 func = fm->GetPowerLaw(24, 72, 0.8);
329 }
330 else if(fitFuncID == kFitPhojet) func = fm->GetHistoFunc(hSpectraMC[kPhojet] [ipart][icharge]);
331 else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
332 else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
333 else if(fitFuncID == kFitCMS6D6T) func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
334 else {
335 cout << "Unknown fit ID " << fitFuncID << endl;
336 return;
337 }
338
339 if(fitFuncID >= kFitPhojet){
340 fitmin = 0.0;
341 fitmax = 1.0;
342 }
343
058137a6 344 // Temp: fit histo with sist errors
22abbb27 345 TH1F * hsyst = new TH1F(*htemplate);
346 TH1F * hsyststat = 0;
55056a39 347 hsyst->SetFillColor(kGray);
3725266c 348 AliBWTools::GetValueAndError(hsyst,hSpectra[iCombInStudy][ipart][icharge],hSystError[iCombInStudy][ipart][icharge],kTRUE);
22abbb27 349 hsyststat= new TH1F(*hsyst);
350 AliBWTools::GetHistoCombinedErrors(hsyststat,hSpectra[iCombInStudy][ipart][icharge]); // combine syst and stat
55056a39 351 hsyststat->SetFillColor(kGray);
3725266c 352
058137a6 353 TH1F * hToFit = 0;
354 if (whatToFit == kStatErrors) hToFit = hSpectra[iCombInStudy][ipart][icharge]; // Shorthand
355 if (sumCharge) hToFit->Add(hSpectra[iCombInStudy][ipart][1]);
356 if (whatToFit == kStatSystErrors) {
22abbb27 357 hToFit = hsyststat;// Shorthand
058137a6 358 }
359 if (whatToFit == kSystErrors) hToFit = hsyst;
360
3725266c 361
362 if(!AliBWTools::Fit(hToFit,func,fitmin,fitmax)) {
4c0d7bc7 363 cout << " FIT ERROR " << endl;
364 return;
365 }
3725266c 366 cout << "DRAWING" << endl;
22abbb27 367 c2->cd();
368 hsyststat->Draw("same,e2");
369 // TH1F * hsyststat2 = new TH1F(*hsyststat);
370 // hsyststat2->SetFillColor(kWhite);
371 // hsyststat2->SetLineColor(kBlack);
372 // hsyststat2->Draw("same,e2");
058137a6 373 hToFit->Draw("same");
3725266c 374 TF1* fitfunc=(TF1*)hToFit->GetListOfFunctions()->At(0);
5f2ca1e6 375 fitfunc->Draw("same");
376 fitfunc->SetRange(0,4);
377 fitfunc->SetLineColor(hSpectra[iCombInStudy][ipart][icharge]->GetLineColor());
c4fe358b 378 if(drawStar) DrawStar(icharge);
058137a6 379 hRatiosToFit[ipart][icharge]=(TH1F*)hToFit->Clone(Form("hRatio%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
2210b002 380 hRatiosToFitSyst[ipart][icharge]=(TH1F*)hsyststat->Clone(Form("hRatioSyst%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
57a944cd 381 // Syst error on ratio with no fill
382 hRatiosToFitSyst[ipart][icharge]->SetLineColor(hToFit->GetLineColor());
383 hRatiosToFitSyst[ipart][icharge]->SetFillStyle(0);
384
385
4d4a268c 386 // Compute ratio data/function integrating the function in the bin width
058137a6 387 for(Int_t iBin=1; iBin<hToFit->GetNbinsX(); iBin++){
388 Double_t lowLim=hToFit->GetBinLowEdge(iBin);
389 Double_t highLim=hToFit->GetBinLowEdge(iBin+1);
5f2ca1e6 390 Double_t contFunc=fitfunc->Integral(lowLim,highLim)/(highLim-lowLim);
058137a6 391 Double_t ratio=hToFit->GetBinContent(iBin)/contFunc;
392 Double_t eratio=hToFit->GetBinError(iBin)/contFunc;
2210b002 393 Double_t eratioSyst=hsyststat->GetBinError(iBin)/contFunc;
394 hRatiosToFitSyst[ipart][icharge]->SetBinContent(iBin,ratio);
395 hRatiosToFitSyst[ipart][icharge]->SetBinError(iBin,eratioSyst);
5f2ca1e6 396 hRatiosToFit[ipart][icharge]->SetBinContent(iBin,ratio);
397 hRatiosToFit[ipart][icharge]->SetBinError(iBin,eratio);
398 }
5f2ca1e6 399 c2->Update();
058137a6 400 l->AddEntry(hToFit,
4c0d7bc7 401 scaleKaons && ipart == kKaon ?
402 (TString(partLabel[ipart][icharge])+" #times 2").Data()
403 : partLabel[ipart][icharge]);
4c0d7bc7 404
405 // populate table
406 // Float_t yield = func->Integral(0.45,1.05);
407 // Float_t yieldE = func->IntegralError(0.45,1.05);
408
409 Float_t yield = func->Integral(0.,100);
410 //Float_t yieldE = func->IntegralError(0.,100);
411
412 Double_t yieldTools, yieldETools;
413 Double_t partialYields[3],partialYieldsErrors[3];
058137a6 414 AliBWTools::GetYield(hToFit, func, yieldTools, yieldETools,
4c0d7bc7 415 0, 100, partialYields,partialYieldsErrors);
416 Double_t tslope = func->GetParameter(2);
417 Double_t tslopeE = func->GetParError(2);
418
5f2ca1e6 419 table.SetNextCol(partLatex[ipart][icharge]);
4c0d7bc7 420 //table.SetNextCol(yield,yieldE,-4);
421 table.SetNextCol(yieldTools, yieldETools,-4);
422 table.SetNextCol(func->GetParameter(0));
423 table.SetNextCol(tslope,tslopeE,-4);
55056a39 424 table.SetNextCol(func->GetParameter(1),func->GetParError(1),-2);
4c0d7bc7 425 table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
058137a6 426 Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hToFit);
4c0d7bc7 427 //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
428 Float_t yieldAbove = func->Integral(lowestPoint,100);
429 table.SetNextCol(lowestPoint,-2);
55056a39 430 table.SetNextCol((1 - yieldAbove/yield)*100,-1);
431 Float_t mean, meane; // Mean using only fit function
432 Double_t meanDF, meanDFe; // Mean from data + fit
4c0d7bc7 433 Float_t mean2, mean2e;
55056a39 434 if (fitFuncID < kFitPhojet) {
435 AliBWTools::GetMean (func, mean, meane , 0.,100., normPar);
436 AliBWTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
437 }
ffcaf7be 438 // AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
55056a39 439 AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
55056a39 440 table.SetNextCol(mean, meane ,-4);
441 table.SetNextCol(meanDF, meanDFe,-4);
4c0d7bc7 442
443 // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
444 table.InsertRow();
445
446
447 /// TEMP TABLE
49ad829f 448 tempTable.SetNextCol(partLatex[ipart][icharge]);
4c0d7bc7 449 tempTable.SetNextCol(yieldTools, yieldETools, -4);
55056a39 450 tempTable.SetNextCol(meanDF, meanDFe ,-4);
49ad829f 451 // tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
452 // tempTable.SetNextCol(yieldAbove/yield,-2);
4c0d7bc7 453 tempTable.InsertRow();
5f2ca1e6 454 c2r->cd();
57a944cd 455 if(divideRatiosCanvas) {
456 c2r->cd(ipart+1);
457 gPad->SetGridy();
458 TH2F * heloc = new TH2F(*hemptyR);
459 heloc->Draw();
460 }
2210b002 461 hRatiosToFitSyst[ipart][icharge]->Draw("e2same");
5f2ca1e6 462 hRatiosToFit[ipart][icharge]->Draw("esame");
463
4c0d7bc7 464 }
5f2ca1e6 465 c2->cd();
2210b002 466 if (flagPreliminary) ALICEWorkInProgress(c2,"","ALICE Preliminary");
5f2ca1e6 467 l->Draw();
2210b002 468 if (flagPreliminary) tf->Draw();
5f2ca1e6 469 c2r->cd();
2210b002 470 if(flagPreliminary) ALICEWorkInProgress(c2r,"","ALICE Preliminary");
4c0d7bc7 471 l->Draw();
0e019ec6 472 PrintCanvas(c2,printFormats);
473 PrintCanvas(c2r,printFormats);
4c0d7bc7 474
5f2ca1e6 475
4c0d7bc7 476 }
477
478
57a944cd 479 // table.PrintTable("ASCII");
480 table.PrintTable("TWIKI");
4c0d7bc7 481
482 cout << "" << endl;
483 tempTable.PrintTable("ASCII");
484
5f2ca1e6 485
486
4c0d7bc7 487}
488
489void LoadSpectra() {
490
3725266c 491 // Loads spectra for different detectors and corresponding systematic errors.
492
4c0d7bc7 493 TFile * f=0;
494
3725266c 495 // Systematic errors: initialize histos
496 gROOT->cd();
497 for(Int_t ipart = 0; ipart < kNPart; ipart++){
498 for(Int_t idet = 0; idet < kNHist; idet++){
499 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
500 hSystError[idet][ipart][icharge] = new TH1F (TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
501 TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
502 nbinsTempl,templBins);
503 hSystError[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
504 hSystError[idet][ipart][icharge]->SetLineColor (color[ipart] );
505
506 }
507
508 }
509
510 }
511
4c0d7bc7 512 // Load
513
514
515 // TOF
5f2ca1e6 516 // Load Efficiencies
4d4a268c 517 f = new TFile("./Files/effTOFhistosSmooth.root");
518 TFile * ftrk = new TFile("./Files/effhistos.root");
5f2ca1e6 519 TH1D * hEffTrackTOF[kNPart][kNCharge];
520 TH1D * hEffMatchTOF[kNPart][kNCharge];
4d4a268c 521 hEffTrackTOF[kPion] [kPos] = (TH1D*) ftrk->Get("hpitrk_pos");
522 hEffTrackTOF[kKaon] [kPos] = (TH1D*) ftrk->Get("hkatrk_pos");
523 hEffTrackTOF[kProton][kPos] = (TH1D*) ftrk->Get("hprtrk_pos");
5f2ca1e6 524 hEffMatchTOF[kPion] [kPos] = (TH1D*) f->Get("hpieff_pos");
525 hEffMatchTOF[kKaon] [kPos] = (TH1D*) f->Get("hkaeff_pos");
526 hEffMatchTOF[kProton][kPos] = (TH1D*) f->Get("hpreff_pos");
4d4a268c 527 hEffTrackTOF[kPion] [kNeg] = (TH1D*) ftrk->Get("hpitrk_neg");
528 hEffTrackTOF[kKaon] [kNeg] = (TH1D*) ftrk->Get("hkatrk_neg");
529 hEffTrackTOF[kProton][kNeg] = (TH1D*) ftrk->Get("hprtrk_neg");
5f2ca1e6 530 hEffMatchTOF[kPion] [kNeg] = (TH1D*) f->Get("hpieff_neg");
531 hEffMatchTOF[kKaon] [kNeg] = (TH1D*) f->Get("hkaeff_neg");
532 hEffMatchTOF[kProton][kNeg] = (TH1D*) f->Get("hpreff_neg");
533
4c0d7bc7 534 // f = new TFile("./Files/spectra-pos-y.root");
5f2ca1e6 535 f = new TFile("./Files/spectraRaw-pos-y.root");
4c0d7bc7 536 hSpectra[kTOF][kPion] [kPos]= (TH1F*) f->Get("hpi");
537 hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
538 hSpectra[kTOF][kKaon] [kPos]= (TH1F*) f->Get("hka");
539 hSpectra[kTOF][kPion] [kPos]->SetName("hpiPos");
540 hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
541 hSpectra[kTOF][kKaon] [kPos]->SetName("hkaPos");
542 //f = new TFile("./Files/spectra-neg-y.root");
5f2ca1e6 543 f = new TFile("./Files/spectraRaw-neg-y.root");
4c0d7bc7 544 hSpectra[kTOF][kPion] [kNeg]= (TH1F*) f->Get("hpi");
545 hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
546 hSpectra[kTOF][kKaon] [kNeg]= (TH1F*) f->Get("hka");
547 hSpectra[kTOF][kPion] [kNeg]->SetName("hpiNeg");
548 hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
549 hSpectra[kTOF][kKaon] [kNeg]->SetName("hkaNeg");
550
5f2ca1e6 551 // Divide for efficiency
552 hSpectra[kTOF][kPion] [kPos]->Divide(hEffTrackTOF[kPion] [kPos]);
553 hSpectra[kTOF][kKaon] [kPos]->Divide(hEffTrackTOF[kKaon] [kPos]);
554 hSpectra[kTOF][kProton][kPos]->Divide(hEffTrackTOF[kProton][kPos]);
555 hSpectra[kTOF][kPion] [kPos]->Divide(hEffMatchTOF[kPion] [kPos]);
556 hSpectra[kTOF][kKaon] [kPos]->Divide(hEffMatchTOF[kKaon] [kPos]);
557 hSpectra[kTOF][kProton][kPos]->Divide(hEffMatchTOF[kProton][kPos]);
558 hSpectra[kTOF][kPion] [kNeg]->Divide(hEffTrackTOF[kPion] [kNeg]);
559 hSpectra[kTOF][kKaon] [kNeg]->Divide(hEffTrackTOF[kKaon] [kNeg]);
560 hSpectra[kTOF][kProton][kNeg]->Divide(hEffTrackTOF[kProton][kNeg]);
561 hSpectra[kTOF][kPion] [kNeg]->Divide(hEffMatchTOF[kPion] [kNeg]);
562 hSpectra[kTOF][kKaon] [kNeg]->Divide(hEffMatchTOF[kKaon] [kNeg]);
563 hSpectra[kTOF][kProton][kNeg]->Divide(hEffMatchTOF[kProton][kNeg]);
564
565
566 // Clean UP TOF spectra, removing unwanted points
4c0d7bc7 567 cout << "Cleaning Up TOF spectra" << endl;
568 Int_t nbin = hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
569 for(Int_t ibin = 1; ibin <= nbin; ibin++){
570 Float_t pt = hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
571 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
572 if (pt > 2.35) {
573 hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
574 hSpectra[kTOF][kKaon][icharge]->SetBinError (ibin,0);
575 hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
576 hSpectra[kTOF][kProton][icharge]->SetBinError (ibin,0);
4d4a268c 577 }
578 if (pt < 0.7) {
579 hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
580 hSpectra[kTOF][kProton][icharge]->SetBinError (ibin,0);
581 }
4c0d7bc7 582 }
583 }
584// cout << "Scaling TOF to TPC" << endl;
585// // Scale protons to TPC
586// hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
587// // Scale pbar so that pbar/p is compatible with Panos
588// hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
589
3725266c 590 // TOF: systematics
591 // Load TOF systematic errors:
592 f = new TFile ("./Files/systMatchingPos.root");
593 AliBWTools::AddHisto(hSystError[kTOF][kPion][kPos] ,(TH1*)gDirectory->Get("hErrPionMatch") );
594 AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hErrProtonMatch"));
595 AliBWTools::AddHisto(hSystError[kTOF][kKaon][kPos] ,(TH1*)gDirectory->Get("hErrKaonMatch") );
596 f = new TFile ("./Files/systMatchingNeg.root");
597 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kNeg],(TH1*)gDirectory->Get("hErrPionMatch"));
598 AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hErrProtonMatch"));
599 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kNeg],(TH1*)gDirectory->Get("hErrKaonMatch"));
600 f = new TFile ("./Files/systPIDall.root");
601 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kPos],(TH1*)gDirectory->Get("hpiCorr"));
602 AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hpCorr"));
603 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kPos],(TH1*)gDirectory->Get("hkCorr"));
604 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kNeg],(TH1*)gDirectory->Get("hpiCorr"));
605 AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hpCorr"));
606 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kNeg],(TH1*)gDirectory->Get("hkCorr"));
4c0d7bc7 607
608
5f2ca1e6 609 // ITS SA
2210b002 610 f = new TFile("./Files/ITSsaSpectraCorr-28Sep10.root");
4c0d7bc7 611 hSpectra[kITS][kPion] [kPos]= (TH1F*) f->Get("hSpectraPos0");
612 hSpectra[kITS][kKaon] [kPos]= (TH1F*) f->Get("hSpectraPos1");
613 hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
614 hSpectra[kITS][kPion] [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
615 hSpectra[kITS][kKaon] [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
616 hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
617
4d4a268c 618 // Remove unwanted bins
4c0d7bc7 619 for(Int_t ipart = 0; ipart < kNPart; ipart++){
620 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
2210b002 621 Int_t nbinITS = hSpectra[kITS][ipart][icharge]->GetNbinsX();
622 for(Int_t ibin = 1; ibin <= nbinITS; ibin++){
4c0d7bc7 623 if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
624 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
625 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
626 }
57a944cd 627 if(ipart == kProton && (ibin==9 || ibin==10)){
5f2ca1e6 628 printf("Kill bin %d (%f - %f GeV/c)for ITS protons\n",ibin,hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin),hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin)+hSpectra[kITS][ipart][icharge]->GetBinWidth(ibin));
629 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
630 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
631 }
57a944cd 632 if(ipart == kKaon && ibin==7){
633 printf("Kill bin %d (%f - %f GeV/c)for ITS kaons\n",ibin,hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin),hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin)+hSpectra[kITS][ipart][icharge]->GetBinWidth(ibin));
634 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
635 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
636 }
637
4c0d7bc7 638// if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
639// hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
640// hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
641// }
642 }
643
644 }
645 }
646
eea641de 647 if(useSecCorrFromDCA){
77233d47 648 TFile* fseccorr = new TFile("./Files/CorrFac-SecProtons3Meth-ITSsa.root");
649 TH1F* hcorrp=(TH1F*)fseccorr->Get("hSecPCorrFromDCAexpo");
650 TH1F* hcorrpbar=(TH1F*)fseccorr->Get("hSecPbarCorrFromDCAexpo");
eea641de 651 hSpectra[kITS][kProton][kPos]->Multiply(hcorrp);
652 hSpectra[kITS][kProton][kNeg]->Multiply(hcorrpbar);
653 fseccorr->Close();
654 }
4c0d7bc7 655
3725266c 656 // Load ITS sa systematics, only pt dependent ones
4d4a268c 657 // f = TFile::Open("./Files/ITSsa-systematics_20100930.root");
658 f = TFile::Open("./Files/ITSsa-systematics-20101014.root");
3725266c 659 for(Int_t ipart = 0; ipart < kNPart; ipart++){
660 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
661 AliBWTools::AddHisto(hSystError[kITS][ipart][icharge], (TH1*) gDirectory->Get(Form("hSystTot%s%s",chargeFlag[icharge],partFlag[ipart]))); // Using total error
662 }
663 }
664
4c0d7bc7 665 // ITS + TPC (Marek)
0e019ec6 666 //f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons20100720.root");
5b103cfc 667 f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons14092010new.root");
4c0d7bc7 668 TList * list = (TList*) gDirectory->Get("output");
669 hSpectra[kITSTPC][kPion] [kPos]= (TH1F*) list->FindObject("Pions");
670 hSpectra[kITSTPC][kKaon] [kPos]= (TH1F*) list->FindObject("Kaons");
671 hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
672 hSpectra[kITSTPC][kPion] [kNeg]= (TH1F*) list->FindObject("AntiPions");
673 hSpectra[kITSTPC][kKaon] [kNeg]= (TH1F*) list->FindObject("AntiKaons");
674 hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
675
676 // TPC
4c0d7bc7 677 f = new TFile("./Files/protonSpectra_20100615.root");
678 hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
679 hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
680 f = new TFile("./Files/pionSpectra_20100615.root");
681 hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
682 hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
5b103cfc 683 // f = new TFile("./Files/kaonSpectra_20100615.root");
684 f = new TFile("./Files/kaonsTPCcorr_2010_08_31.root");
4c0d7bc7 685 hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
686 hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
687
688 // Clean UP TPC spectra, removing unwanted points
689 cout << "Cleaning Up TPC spectra" << endl;
690 nbin = hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
691 for(Int_t ibin = 0; ibin < nbin; ibin++){
692 Float_t pt = hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
57a944cd 693 if (pt > 0.45 || pt<0.25) {
4c0d7bc7 694 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
57a944cd 695 hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
696 hSpectra[kTPC][kKaon][icharge]->SetBinError (ibin,0);
4c0d7bc7 697 }
698 }
57a944cd 699 // if (pt < 0.25) {
700 // hSpectra[kTPC][kKaon][kNeg]->SetBinContent(ibin,0);
701 // hSpectra[kTPC][kKaon][kNeg]->SetBinError (ibin,0);
702 // }
4c0d7bc7 703 if (pt < 0.45) {
704 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
705 hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
706 hSpectra[kTPC][kProton][icharge]->SetBinError (ibin,0);
707 }
708 }
709 }
710
3725266c 711 // Load TPC systematics
712 cout << "WARNING: TPC SYST FOR NEGATIVES TO BE CORRECTED" << endl;
713 f = TFile ::Open("./Files/pionsSystSum.root");
714 AliBWTools::AddHisto(hSystError[kTPC][kPion][kPos],(TH1*) gDirectory->Get("pionsSyst"));
715 AliBWTools::AddHisto(hSystError[kTPC][kPion][kNeg],(TH1*) gDirectory->Get("pionsSyst"));
716 f = TFile ::Open("./Files/kaonsSystSum.root");
717 AliBWTools::AddHisto(hSystError[kTPC][kKaon][kPos],(TH1*) gDirectory->Get("kaonsSyst"));
718 AliBWTools::AddHisto(hSystError[kTPC][kKaon][kNeg],(TH1*) gDirectory->Get("kaonsSyst"));
719 f = TFile ::Open("./Files/ProtonsSystSum.root");
720 AliBWTools::AddHisto(hSystError[kTPC][kProton][kPos],(TH1*) gDirectory->Get("ProtonsSyst"));
721 AliBWTools::AddHisto(hSystError[kTPC][kProton][kNeg],(TH1*) gDirectory->Get("ProtonsSyst"));
722
4c0d7bc7 723 // K0s
49ad829f 724 f = new TFile ("./Files/PtSpectraCorrectedK0sOff_20100803.root");
4c0d7bc7 725 // hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
726 hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
727 // hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
728 // hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
729 hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
730
49ad829f 731 // Kinks
4c0d7bc7 732 // f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
49ad829f 733 // f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
22abbb27 734 // f = new TFile ("./Files/KaonKinkJun16PhySel2N.root");
735 f = new TFile ("./Files/KaonKinkJun16PhySel2NRebPass6.root");
4c0d7bc7 736 hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
737 hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
49ad829f 738 // hSpectra[kKinks][kKaon][kPos]->Scale(0.5/0.7); // different rapidity range for kinks
739 // hSpectra[kKinks][kKaon][kNeg]->Scale(0.5/0.7); // different rapidity range for kinks
740 // hSpectra[kKinks][kKaon][kPos]->Scale(276004./263345.); // different N of events
741 // hSpectra[kKinks][kKaon][kNeg]->Scale(276004./263345.); // different N of events
22abbb27 742 // This normalization if for KaonKinkJun16PhySel2N
743 // hSpectra[kKinks][kKaon][kPos]->Scale(1./303214); // different N of events
744 // hSpectra[kKinks][kKaon][kNeg]->Scale(1./303214); // different N of events
745 // This normalization is for KaonKinkJun16PhySel2NRebPass6
746 hSpectra[kKinks][kKaon][kPos]->Scale(1./260432.26); // different N of events
747 hSpectra[kKinks][kKaon][kNeg]->Scale(1./260432.26); // different N of events
4c0d7bc7 748
5f2ca1e6 749 // Apply correction factors
4c0d7bc7 750 // Secondaries for protons
5f2ca1e6 751
0e019ec6 752 //f = new TFile ("./Files/corrFactorProtons_20100615.root");
2210b002 753 f = new TFile ("./Files/corrFactorProtons_2010_09_24.root");
4c0d7bc7 754 TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
755 if(correctSecondaries) {
756 cout << "CORRECTING SECONDARIES" << endl;
757
758 for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
759 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
760 Int_t ipart = kProton;
761 TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
762 if (h){
763 Int_t nbins = h->GetNbinsX();
764 for(Int_t ibin = 0; ibin < nbins; ibin++){
765 // Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
766 Float_t pt = h->GetBinCenter(ibin);
767 if (icharge == kNeg) pt = -pt;
768 Int_t binCorrection = hCorrSecondaries->FindBin(pt);
769 Float_t correction = hCorrSecondaries->GetBinContent(binCorrection);
770 // cout << pt << " " << correction << endl;
771
c4fe358b 772 if (correction) {// If the bin is empty this is a 0
4c0d7bc7 773 h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
774 h->SetBinError (ibin,h->GetBinError (ibin)/correction);
775 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
776 cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
777 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
778
779 }
780 }
781 }
782 }
783 }
784 }
785 // geant/fluka absorption
786 if(correctGeantFlukaXS) {
787 cout << "CORRECTING GEANT3/FLUKA" << endl;
eea641de 788 // common file for Kaons
789 TFile *fFlukakaons = TFile::Open("./Files/correctionForCrossSection.321.root");
790 TH1D * hCorrFlukakaon[kNCharge];
791 hCorrFlukakaon[kPos] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionParticles");
792 hCorrFlukakaon[kNeg] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionAntiParticles");
793
794 for(Int_t idet = 0; idet < kNDet; idet++){
795 if( idet != kITS) continue; // comment to use fluka for kaons for all dets
796 if (idet == kTOF) continue; // TOF already corrected
797 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
798 Int_t ipart = kKaon;
2f7f947e 799 TH1 * h = hSpectra[idet][ipart][icharge]; // only ITS sa
eea641de 800 if (h){
801 Int_t nbins = h->GetNbinsX();
802 Int_t nbinsy=hCorrFlukakaon[icharge]->GetNbinsY();
803 for(Int_t ibin = 0; ibin < nbins; ibin++){
804 Float_t pt = h->GetBinCenter(ibin);
805 Float_t minPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(1);
806 Float_t maxPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(nbinsy+1);
807 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
808 if (pt > maxPtCorrection) pt = maxPtCorrection;
809 Float_t correction = hCorrFlukakaon[icharge]->GetBinContent(hCorrFlukakaon[icharge]->GetXaxis()->FindBin(pt));
810 if (correction != 0) {// If the bin is empty this is a 0
811 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
812 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
813 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
814 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
815 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
816 }
817 }
818 }
819 }
820 }
821
2f7f947e 822 // PROTONS
823
eea641de 824 // ITS specific file for protons/antiprotons
5f2ca1e6 825 TFile* fITS = new TFile ("./Files/correctionForCrossSectionITS_20100719.root");
826 TH2D * hCorrFlukaITS[kNCharge];
827 hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
828 hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
829
830 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
831 Int_t ipart = kProton;
832 TH1 * h = hSpectra[kITS][ipart][icharge]; // only ITS sa
833 if (h){
834 Int_t nbins = h->GetNbinsX();
835 Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
836 for(Int_t ibin = 0; ibin < nbins; ibin++){
837 Float_t pt = h->GetBinCenter(ibin);
838 Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
839 Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
840 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
841 if (pt > maxPtCorrection) pt = maxPtCorrection;
842 Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
843 if (correction != 0) {// If the bin is empty this is a 0
844 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
845 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
846 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
847 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
848 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
849 }
850 }
851 }
852 }
853
4c0d7bc7 854
77233d47 855
856 //f = new TFile ("./Files/correctionForCrossSection_20100615.root");
857 f = new TFile ("./Files/correctionForCrossSection_20100920.root");
4c0d7bc7 858 TH2D * hCorrFluka[kNCharge];
859 hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
860 hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
861 for(Int_t idet = 0; idet < kNDet; idet++){
862 if (idet == kITS) continue; // skip ITS sa
863 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
864 Int_t ipart = kProton;
865 TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
866 if (h){
867 Int_t nbins = h->GetNbinsX();
868 for(Int_t ibin = 0; ibin < nbins; ibin++){
869// Float_t pt = convertToMT ?
870// TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
871// h->GetBinCenter(ibin);
872 Float_t pt = h->GetBinCenter(ibin);
873 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
2210b002 874 // hCorrFluka[icharge]->Draw();
875 TH1D * htmpFluka=hCorrFluka[icharge]->ProjectionY();
876 Float_t maxPtCorrection = AliBWTools::GetHighestNotEmptyBinEdge(htmpFluka);//->GetYaxis()->GetBinLowEdge(hCorrFluka[icharge]->GetNbinsY()+1);
4c0d7bc7 877 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
2210b002 878 if (pt > maxPtCorrection) pt = maxPtCorrection-0.0001;
4c0d7bc7 879 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
77233d47 880
5f2ca1e6 881 // already in the efficiency correction (F. Noferini)
2f7f947e 882 if (idet == kTOF) {
49ad829f 883 if (icharge == kNeg) correction = 1; // antiprotons already corrected in efficiency
884 // Scale correction for the different material budget. Recipe by Francesco Noferini
2210b002 885 else {
886 correction = TMath::Power(correction,0.07162/0.03471);
887 }
49ad829f 888 }
4c0d7bc7 889 if (correction != 0) {// If the bin is empty this is a 0
77233d47 890
4c0d7bc7 891 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
892 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
4c0d7bc7 893 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
894 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
895 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
896 }
897
898 }
899
900 }
901 }
902 }
903 }
904
905
906 // Set style and scale
907 for(Int_t idet = 0; idet < kNDet; idet++){
908 for(Int_t ipart = 0; ipart < kNPart; ipart++){
909 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
910 if (hSpectra[idet][ipart][icharge]){
911 hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
912 hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
913 hSpectra[idet][ipart][icharge]->SetLineColor (color[ipart] );
914 hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
915 hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
916 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
917 if (convertToMT) {
918 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
919 hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
920 hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
921 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
922 }
923// if (idet == kTOF || idet == kTPC) {
924// hSpectra[idet][ipart][icharge]->Scale(1./236763);
925// }
926// if (idet == kITS){
927// hSpectra[idet][ipart][icharge]->Scale(202945./236763);
928// }
929 if (scaleKaons && ipart == kKaon) {
930 hSpectra[idet][ipart][icharge]->Scale(2.);
931 }
932 } else {
933 printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
934 }
935 }
936 }
937 }
938
ee08b77d 939
940 // Create fake weights for the mean; To be update once we have syst errors
2f7f947e 941 // Using syste from table in paper. It would be better to have this as a function of pt.
942 TH1F * hWeights[3][kNPart];
943 const Double_t kWeights[3][kNPart] =
3725266c 944 // {{4, 3, 10.2}, // TPC
945 // {4.1,8.8,7.0 }, //TOF
946 // {4.5,5.6,7.0 }}; // ITS
947 {{0.1,0.1,0.1}, // TPC
948 {0.1,0.1,0.1}, //TOF
949 {0.1,0.1,0.1}}; // ITS
2f7f947e 950 for(Int_t ipart = 0; ipart <= kNPart ; ipart++){
951 for(Int_t idet = 0; idet <= kITS ; idet++){
952 hWeights[idet][ipart] = (TH1F*) hSpectra[idet][ipart][kPos]->Clone();
2210b002 953 Int_t nbinW = hWeights[idet][ipart]->GetNbinsX();
954 for(Int_t ibin = 1; ibin <= nbinW; ibin++){
2f7f947e 955 hWeights[idet][ipart]->SetBinError(ibin, kWeights[idet][ipart]);
956 }
957 }
ee08b77d 958 }
49ad829f 959 const Double_t scaleDet[] = {1.00,1.00,1.00}; // During the weekly meeting 19/08/2010 we decided not to apply any scaling.
960 // const Double_t scaleDet[] = {0.98,1,1.00}; // Scaling factor for the different detectors. Estimated from ratios, it has an estimated uncertainty of ~2%
2f7f947e 961 // const Double_t scaleDet[] = {0.88,1,0.88}; // Scaling factor for the different detectors. Estimated from ratios, it has an estimated uncertainty of ~2%
ee08b77d 962
963
4c0d7bc7 964 // Combine detectors
965 for(Int_t ipart = 0; ipart < kNPart; ipart++){
966 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
967 TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
968
969 if (convertToMT) {
970 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
971 htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
972
973 }
974 hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
975 hSpectra[kTPC][ipart][icharge],
976 htemplLocal,1.);;
3725266c 977
ee08b77d 978 hSpectra[kCombAll][ipart][icharge] =
3725266c 979 AliBWTools::Combine3HistosWithErrors(hSpectra[kITS][ipart][icharge], // Histos to combine
ee08b77d 980 hSpectra[kTPC][ipart][icharge],
981 hSpectra[kTOF][ipart][icharge],
3725266c 982 hSystError[kITS][ipart][icharge], // Errors (weights) used for the average
983 hSystError[kTPC][ipart][icharge],
984 hSystError[kTOF][ipart][icharge],
985 // hWeights[kITS][ipart],
986 // hWeights[kTPC][ipart],
987 // hWeights[kTOF][ipart],
988 htemplLocal,1, // template, take statistical error from TPC in overlap
989 scaleDet[kITS], // relative scaling
2f7f947e 990 scaleDet[kTPC],
3725266c 991 scaleDet[kTOF],
992 (TH1**)&hSystError[kCombAll][ipart][icharge], // put combined syst error here
993 1 // weights histos contain error in bin content
2f7f947e 994 );
4c0d7bc7 995// if (convertToMT) {
996// TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
997// hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
998// hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
999// hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
1000// }
1001 }
1002 }
1003
1004
1005 // Scale for the number of inelastic collisions and correct for
1006 // efficiency losses due to physics selection:
1007
1008 Double_t effPhysSel[kNPart];
1009 effPhysSel[kPion] = 1.012;
1010 effPhysSel[kKaon] = 1.013;
1011 effPhysSel[kProton] = 1.014;
1012
1013
1014 for(Int_t idet = 0; idet < kNHist; idet++){
1015 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1016 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1017 if(hSpectra[idet][ipart][icharge]) {
1018 // cout << "Scaling!" << endl;
49ad829f 1019 if(idet != kKinks && idet != kK0){ // Kinks use a different run list, k0s normalized by Helene
4d4a268c 1020 hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel efficiency (computed by Alexander)
5f2ca1e6 1021 }
4c0d7bc7 1022 }
1023 }
1024 }
1025 }
1026
1027
1028}
1029
1030void LoadMC() {
1031
1032 TFile * f = 0;
1033 const char * evClass= "INEL";
1034 const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
1035 "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
1036 "./Files/dndeta_CMSD6T109_10M_900GeV.root",
1037 "./Files/dndeta_Perugia0320_10M_900GeV.root", };
1038
1039 // Phojet
1040 for(Int_t itune = 0; itune < kNTunes; itune++){
1041 f = new TFile(files[itune]);
1042 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1043 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1044 hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
1045 }
1046 }
1047 }
1048
1049
1050 // Set style
1051 for(Int_t itune = 0; itune < kNTunes; itune++){
1052 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1053 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1054 if (hSpectraMC[itune][ipart][icharge]){
1055 hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
1056 hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
1057 hSpectraMC[itune][ipart][icharge]->SetLineColor (mcLineColor[itune] );
1058 hSpectraMC[itune][ipart][icharge]->SetLineStyle (mcLineStyle[itune] );
1059 hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
1060 if (convertToMT) {
1061 TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
1062 hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
1063 hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
1064 hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
1065 }
1066
1067 } else {
1068 printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
1069 }
1070 }
1071 }
1072 }
1073
1074}
1075
1076void DrawStar(Int_t icharge) {
1077 // cout << icharge << endl;
1078
1079 // gROOT->LoadMacro("StarPPSpectra.C");
1080 TGraphErrors ** gStar = StarPPSpectra();
1081
1082 for(Int_t istar = 0; istar < 6; istar++){
1083 gStar[istar]->SetMarkerStyle(kOpenStar);
1084 if (icharge==kPos && (istar%2) ) gStar[istar]->Draw("P");
1085 else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
1086 else cout << "Skipping Star " << istar << endl;
1087 }
1088}
1089
1090void GetITSResiduals() {
1091
1092
1093 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1094 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1095 // cout << "1 " << ipart << " " << icharge << endl;
1096
1097// gSystem->ProcessEvents();
1098// gSystem->Sleep(1000);
1099 TF1* f = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
1100 TH1F * hres1, *hres2;
1101 AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
1102
1103 TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
1104 c1->SetLogy();
1105 hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
1106 hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
1107 hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
1108 hSpectra[kITS][ipart][icharge]->Draw("same");
1109 hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
1110 TLegend* l = new TLegend( 0.182886, 0.192308, 0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
1111 l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
1112 l->AddEntry(hSpectra[kITS][ipart][icharge], "ITS");
1113 l->AddEntry(f, "Fit to TOF+TPC");
1114 l->Draw();
1115
1116 TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
1117 TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
1118 c2->SetGridy();
1119 hres2->SetMinimum(-1);
1120 hres2->SetMaximum(1);
1121 hres2->Draw();
1122 hres2->GetYaxis()->SetTitleOffset(1.2);
1123 Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
1124 TLine * line = new TLine(x,-1,x,1);
1125 line->SetLineStyle(kDashed);
1126 line->Draw("same");
1127
1128 if (doPrint) {
1129 c1->Update();
1130 c2->Update();
1131 gSystem->ProcessEvents();
1132 c1->Print(TString(c1->GetName()) + ".png");
1133 c2->Print(TString(c2->GetName()) + ".png");
1134 }
1135 }
1136 }
1137}
1138
1139void DrawWithModels() {
1140
4d4a268c 1141 Int_t chargeLoop = sumCharge ? 1 : 2;
1142 for(Int_t icharge = 0; icharge < chargeLoop; icharge++){
4c0d7bc7 1143
1144 // Draw with models
1145 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1146 // Multipad canvas
4d4a268c 1147 TString chargeFlagLocal = chargeFlag[icharge];
1148 if(sumCharge) chargeFlagLocal += chargeFlag[1];
1149 TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlagLocal,TString("cSpectra")+partFlag[ipart]+chargeFlagLocal,700,700);
1150 TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlagLocal, "p1", 0.0, 0.35, 1.0, 0.95, 0, 0, 0);
4c0d7bc7 1151 p1->SetBottomMargin(0);
1152 p1->Draw();
1153
4d4a268c 1154 TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlagLocal, "p2", 0.0, 0.05, 1.0, 0.35, 0, 0, 0);
4c0d7bc7 1155 p2->SetTopMargin(0);
5f2ca1e6 1156 p2->SetBottomMargin(0.3);
4c0d7bc7 1157 p2->Draw();
1158
1159 Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
1160
1161 // Draw spectra
1162 p1->cd();
1163 p1->SetLogy();
4d4a268c 1164 Float_t maxy = sumCharge ? 10 : 5;
1165 TH2F * hempty = new TH2F(TString("hempty")+long(icharge)+long(ipart),"hempty",100,0.,3, 100, 0.0015,maxy);
4c0d7bc7 1166 hempty->SetXTitle("p_{t} (GeV/c)");
1167 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1168 hempty->Draw();
1169 c1->SetLogy();
d8e6677d 1170
4c0d7bc7 1171
1172
d8e6677d 1173 TLegend * l =new TLegend(0.502874, 0.493056, 0.892241, 0.904762);
4d4a268c 1174 l->SetFillColor(kWhite);
1175 if(sumCharge) {
1176 hSpectra[iCombInStudy][ipart][kPos]->Add(hSpectra[iCombInStudy][ipart][kNeg]); // Draw pos+neg
1177 }
ee08b77d 1178 hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
d8e6677d 1179 l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data"));
4c0d7bc7 1180 for(Int_t itune = 0; itune < kNTunes; itune++){
4d4a268c 1181 if(sumCharge) hSpectraMC[itune][ipart][icharge]->Add(hSpectraMC[itune][ipart][kNeg]); // Draw pos+neg;
4c0d7bc7 1182 l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1183 hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
d8e6677d 1184 AliBWTools::GetGraphFromHisto(hSpectraMC[itune][ipart][icharge])->Draw("CX");
4c0d7bc7 1185 }
1186 l->Draw("same");
d8e6677d 1187
4d4a268c 1188 TLatex * tex = new TLatex(0.6712643,2.353486,sumCharge ? Form("%s+%s",partLabel[ipart][icharge],partLabel[ipart][1]) : partLabel[ipart][icharge]);
d8e6677d 1189 tex->SetTextFont(42);
1190 tex->SetTextSize(0.07936508);
1191 tex->SetLineWidth(2);
1192 tex->Draw();
4c0d7bc7 1193
1194 // Draw ratio
1195 p2->cd();
4d4a268c 1196 Float_t maxyr = sumCharge ? 3.99 : 2.99;
1197 TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(icharge)+long(ipart),"hempty",100,0.,3, 100, 0.01,maxyr);
4c0d7bc7 1198 hemptyr->SetXTitle("p_{t} (GeV/c)");
1199 hemptyr->SetYTitle("Data/MC");
1200 hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
1201 hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
1202 hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
1203 hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
1204 hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
d8e6677d 1205 hemptyr->GetXaxis()->SetTitleOffset(1.05);
4c0d7bc7 1206 hemptyr->SetTickLength(0.03*scaleFonts, "X");
1207 hemptyr->SetTickLength(0.02*scaleFonts, "Y");
1208 // hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
1209 hemptyr->GetYaxis()->SetNdivisions(5);
1210 hemptyr->Draw("");
1211
1212 AliBWFunc fm;
1213 for(Int_t itune = 0; itune < kNTunes; itune++){
1214 TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
1215
1216 // l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
ee08b77d 1217 TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[iCombInStudy][ipart][icharge],f);
4c0d7bc7 1218 hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
1219 hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
1220 hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
4d4a268c 1221 hRatio->SetMarkerStyle(0);
1222 hRatio->Draw("same");
1223 // AliBWTools::GetGraphFromHisto(hRatio)->Draw("CX");
4c0d7bc7 1224 }
1225
1226
1227 // print
4d4a268c 1228 PrintCanvas(c1,printFormats);
4c0d7bc7 1229 }
1230 }
1231
1232
1233
1234}
1235
1236void DrawAllAndKaons() {
1237
22abbb27 1238
4c0d7bc7 1239
1240 // gROOT->LoadMacro("ALICEWorkInProgress.C");
1241
eea641de 1242 // gStyle->SetOptFit(0);
1243 gStyle->SetStatX(0.9);
1244 gStyle->SetStatY(0.88);
1245 gStyle->SetStatW(0.4);
1246 gStyle->SetStatH(0.1);
4c0d7bc7 1247
ee08b77d 1248 TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1249 hKaonsAllTPCTOF->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
4c0d7bc7 1250
1251 TH1F * hK0Scaled = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
1252 hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
1253
1254 hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
5f2ca1e6 1255 hSpectra[kKinks][kKaon][kPos]->SetLineColor(4);
4c0d7bc7 1256 hSpectra[kKinks][kKaon][kPos]->SetStats(0);
1257 TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
1258 hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
1259
5f2ca1e6 1260 TCanvas * c1 = new TCanvas("cKaons","cKaons",700,700);
4c0d7bc7 1261 c1->SetLogy();
5f2ca1e6 1262 TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-3,6);
4c0d7bc7 1263 hempty->SetXTitle("p_{t} (GeV/c)");
22abbb27 1264 // hempty->SetYTitle("dN / dp_{t} (A.U.)");
1265 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
4c0d7bc7 1266 hempty->Draw();
4c0d7bc7 1267 hK0Scaled->Draw("same");
5f2ca1e6 1268 hKaonsAllTPCTOF->Draw("same");
4c0d7bc7 1269 hKinksAll->Draw("same");
22abbb27 1270
4c0d7bc7 1271
22abbb27 1272 TLegend * leg = new TLegend(0.456897, 0.71875, 0.892241, 0.936012,NULL,"brNDC");
5f2ca1e6 1273 // leg->SetBorderSize(0);
4c0d7bc7 1274// leg->SetLineColor(1);
1275// leg->SetLineStyle(1);
1276// leg->SetLineWidth(1);
1277// leg->SetFillColor(19);
5f2ca1e6 1278 leg->SetFillColor(0);
1279 TLegendEntry *entry=leg->AddEntry(hKaonsAllTPCTOF,"K^{+} + K^{-}, ITS+TPC+TOF ","lpf");
1280 entry=leg->AddEntry(hK0Scaled,"K^{0} #times 2","lpf");
1281 entry=leg->AddEntry(hKinksAll,"K^{+} + K ^{-}, Kinks","lpf");
4c0d7bc7 1282 leg->Draw();
1283
22abbb27 1284 // ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Prelimiary}{Statistical Error Only}");
49ad829f 1285 // TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
1286 // tex->SetTextColor(2);
1287 // tex->SetTextFont(42);
1288 // tex->SetTextSize(0.03496503);
1289 // tex->Draw();
4c0d7bc7 1290
1291 c1->Update();
22abbb27 1292 PrintCanvas(c1, printFormats);
4c0d7bc7 1293
1294 // Draw all "stable" hadrons
78033eeb 1295 Bool_t divideCanvas=kTRUE;
4c0d7bc7 1296 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
78033eeb 1297 TCanvas * c1h = 0;
1298 if(divideCanvas) c1h = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],1200,500);
1299 else c1h = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],700,700);
1300 // c1h->SetLogy();
2210b002 1301 c1h->SetLeftMargin(0.14);
78033eeb 1302 if(divideCanvas) {
1303 c1h->Divide(3,1);
1304 c1h->cd(1);
1305 }
2210b002 1306 TH2F * hemptyLoc = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-3,10);
1307 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1308 hemptyLoc->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1309 hemptyLoc->GetYaxis()->SetTitleOffset(1.35);
1310 hemptyLoc->GetXaxis()->SetTitleOffset(1.1);
1311 hemptyLoc->Draw();
d8e6677d 1312 leg = new TLegend(0.482759, 0.489583, 0.896552, 0.925595, NULL,"brNDC");
5f2ca1e6 1313 leg->SetFillColor(0);
4c0d7bc7 1314 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
78033eeb 1315 if(divideCanvas) {
1316 cout << "CD!!" << ipart+1 << endl;
1317 c1h->cd(ipart+1);
1318 TH2F * hClone = new TH2F(*hemptyLoc);
1319 hClone->GetXaxis()->SetRangeUser(0,1);
1320 if(ipart == kPion) hClone->GetYaxis()->SetRangeUser(0,3.5);
1321 if(ipart == kKaon) hClone->GetYaxis()->SetRangeUser(0,0.28);
1322 if(ipart == kProton) hClone->GetYaxis()->SetRangeUser(0,0.1);
1323 hClone->Draw();
1324 }
4c0d7bc7 1325 for(Int_t idet = 0; idet <= kITSTPC; idet++){
1326 // if (idet == kITS) continue;
1327 // if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
1328 hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
1329 hSpectra[idet][ipart][icharge]->Draw("same");
1330 leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet] + ")","lpf");
1331 }
1332 // leg->AddLine();
78033eeb 1333 }
1334 if(divideCanvas) c1h->cd(1);
4c0d7bc7 1335 leg->Draw();
2210b002 1336 // ALICEWorkInProgress(c1h,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1337 PrintCanvas(c1h,printFormats);
4c0d7bc7 1338 }
1339
1340
eea641de 1341 // Draw ratios
4c0d7bc7 1342
eea641de 1343 // K-/K+ in the different detectors
5f2ca1e6 1344 TCanvas * cpm=new TCanvas("cpm","Kminus/Kplus",700,700);
1345 cpm->Divide(2,2);
1346 cpm->cd(1);
1347 TH1F* hRatioKPKM_TPC=new TH1F(*(hSpectra[kTPC][kKaon][kNeg]));
1348 hRatioKPKM_TPC->SetMinimum(0.5);
1349 hRatioKPKM_TPC->SetMaximum(1.5);
1350 hRatioKPKM_TPC->Divide(hSpectra[kTPC][kKaon][kPos]);
1351 hRatioKPKM_TPC->GetYaxis()->SetTitle("K-/K+ (TPC)");
1352 hRatioKPKM_TPC->Draw();
1353 cpm->cd(2);
1354 TH1F* hRatioKPKM_ITS=new TH1F(*(hSpectra[kITS][kKaon][kNeg]));
1355 hRatioKPKM_ITS->Divide(hSpectra[kITS][kKaon][kPos]);
1356 hRatioKPKM_ITS->SetMinimum(0.5);
1357 hRatioKPKM_ITS->SetMaximum(1.5);
1358 hRatioKPKM_ITS->GetYaxis()->SetTitle("K-/K+ (ITSsa)");
1359 hRatioKPKM_ITS->Draw("");
1360 cpm->cd(3);
1361 TH1F* hRatioKPKM_TOF=new TH1F(*(hSpectra[kTOF][kKaon][kNeg]));
1362 hRatioKPKM_TOF->Divide(hSpectra[kTOF][kKaon][kPos]);
1363 hRatioKPKM_TOF->SetMinimum(0.5);
1364 hRatioKPKM_TOF->SetMaximum(1.5);
1365 hRatioKPKM_TOF->GetYaxis()->SetTitle("K-/K+ (TOF)");
1366 hRatioKPKM_TOF->Draw("");
1367 cpm->cd(4);
1368 TH1F* hRatioKPKM_ITSTPC=new TH1F(*(hSpectra[kITSTPC][kKaon][kNeg]));
1369 hRatioKPKM_ITSTPC->Divide(hSpectra[kITSTPC][kKaon][kPos]);
1370 hRatioKPKM_ITSTPC->SetMinimum(0.5);
1371 hRatioKPKM_ITSTPC->SetMaximum(1.5);
1372 hRatioKPKM_ITSTPC->GetYaxis()->SetTitle("K-/K+ (ITS Global)");
1373 hRatioKPKM_ITSTPC->Draw("");
4c0d7bc7 1374
eea641de 1375
1376 // ITS/TPC
5f2ca1e6 1377 TH1F * hRatioITSTPC[kNPart][kNCharge];
eea641de 1378 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1379 // Create canvas
2210b002 1380 TCanvas * c1h = new TCanvas(TString("cITSTPCRatio_")+chargeFlag[icharge],TString("cITSTPCRatio_")+chargeFlag[icharge],1200,500);
1381 c1h->Divide(3,1);
1382 c1h->SetGridy();
1383 TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITSsa/TPC ",100,0.,1., 100, 0.5,1.5);
1384 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1385 hemptyLoc->SetYTitle("ITSsa / TPC");
eea641de 1386 // Loop over particles
5f2ca1e6 1387 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
eea641de 1388 // Clone histo
5f2ca1e6 1389 hRatioITSTPC[ipart][icharge]=new TH1F(*hSpectra[kITS][ipart][icharge]);
1390 Int_t nBinsITS=hSpectra[kITS][ipart][icharge]->GetNbinsX();
1391 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
eea641de 1392 // Loop over ITS bins,
5f2ca1e6 1393 for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1394 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1395 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1396 Float_t lowPtITS=hSpectra[kITS][ipart][icharge]->GetBinLowEdge(iBin);
1397 Float_t binWidITS=hSpectra[kITS][ipart][icharge]->GetBinWidth(iBin);
eea641de 1398 // Loop over TPC bins and find overlapping bins to ITS
5f2ca1e6 1399 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1400 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1401 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1402 if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1403 Float_t numer=hSpectra[kITS][ipart][icharge]->GetBinContent(iBin);
1404 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1405 Float_t enumer=hSpectra[kITS][ipart][icharge]->GetBinError(iBin);
78033eeb 1406 Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
5f2ca1e6 1407 Double_t ratio=0.;
1408 Double_t eratio=0.;
1409 if(numer>0. && denom>0.){
1410 ratio=numer/denom;
1411 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1412 }
1413 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1414 hRatioITSTPC[ipart][icharge]->SetBinError(iBin,eratio);
1415 break;
1416 }
1417 }
1418 }
2210b002 1419 c1h->cd(ipart+1);
1420 // hemptyLoc->SetStats(1);
1421 // hemptyLoc->Draw();
1422 hRatioITSTPC[ipart][icharge]->SetYTitle("ITSsa/TPC");
eea641de 1423 hRatioITSTPC[ipart][icharge]->SetStats(1);
1424 hRatioITSTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1425 hRatioITSTPC[ipart][icharge]->Draw("");
1426 hRatioITSTPC[ipart][icharge]->Fit("pol0","","same");
1427
5f2ca1e6 1428 }
2210b002 1429 PrintCanvas(c1h,printFormats);
5f2ca1e6 1430 }
1431
5b103cfc 1432 // ITS/TPC
1433 TH1F * hRatioITSGlobalTPC[kNPart][kNCharge];
1434 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1435 // Create canvas
2210b002 1436 TCanvas * c1h = new TCanvas(TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],1200,500);
1437 c1h->Divide(3,1);
1438 c1h->SetGridy();
1439 TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITS Global/TPC ",100,0.,1., 100, 0.5,1.5);
1440 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1441 hemptyLoc->SetYTitle("ITS Global / TPC");
5b103cfc 1442 // Loop over particles
1443 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1444 // Clone histo
1445 hRatioITSGlobalTPC[ipart][icharge]=new TH1F(*hSpectra[kITSTPC][ipart][icharge]);
1446 Int_t nBinsITS=hSpectra[kITSTPC][ipart][icharge]->GetNbinsX();
1447 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1448 // Loop over ITS bins,
1449 for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1450 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1451 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1452 Float_t lowPtITS=hSpectra[kITSTPC][ipart][icharge]->GetBinLowEdge(iBin);
1453 Float_t binWidITS=hSpectra[kITSTPC][ipart][icharge]->GetBinWidth(iBin);
1454 // Loop over TPC bins and find overlapping bins to ITS
1455 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1456 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1457 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1458 if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1459 Float_t numer=hSpectra[kITSTPC][ipart][icharge]->GetBinContent(iBin);
1460 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1461 Float_t enumer=hSpectra[kITSTPC][ipart][icharge]->GetBinError(iBin);
78033eeb 1462 Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
5b103cfc 1463 Double_t ratio=0.;
1464 Double_t eratio=0.;
1465 if(numer>0. && denom>0.){
1466 ratio=numer/denom;
1467 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1468 }
1469 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1470 hRatioITSGlobalTPC[ipart][icharge]->SetBinError(iBin,eratio);
1471 break;
1472 }
1473 }
1474 }
2210b002 1475 c1h->cd(ipart+1);
1476 // hemptyLoc->SetStats(1);
1477 // hemptyLoc->Draw();
1478 hRatioITSGlobalTPC[ipart][icharge]->SetYTitle("ITS Global/TPC");
5b103cfc 1479 hRatioITSGlobalTPC[ipart][icharge]->SetStats(1);
1480 hRatioITSGlobalTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1481 hRatioITSGlobalTPC[ipart][icharge]->Draw("");
1482 hRatioITSGlobalTPC[ipart][icharge]->Fit("pol0","","same");
1483
1484 }
2210b002 1485 PrintCanvas(c1h,printFormats);
5b103cfc 1486 }
1487
eea641de 1488 // TOF/TPC
5f2ca1e6 1489 TH1F * hRatioTOFTPC[kNPart][kNCharge];
eea641de 1490 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1491 // create canvas
1492 TCanvas * c1t = new TCanvas(TString("cTOFTPCRatio_")+chargeFlag[icharge],TString("cTOFTPCRatio_")+chargeFlag[icharge],1200,500);
5f2ca1e6 1493 c1t->SetGridy();
eea641de 1494 c1t->Divide(3,1);
5f2ca1e6 1495 TH2F * hemptyt = new TH2F(TString("hemptyRt")+long(icharge),"TOF/TPC ",100,0.,1., 100, 0.5,1.5);
1496 hemptyt->SetXTitle("p_{t} (GeV/c)");
1497 hemptyt->SetYTitle("TOF / TPC");
eea641de 1498 // hemptyt->Draw();
1499 for(Int_t ipart = 0; ipart < kNPart; ipart++) { // loop over particles
1500 // Clone histo
5f2ca1e6 1501 hRatioTOFTPC[ipart][icharge]=new TH1F(*hSpectra[kTOF][ipart][icharge]);
1502 Int_t nBinsTOF=hSpectra[kTOF][ipart][icharge]->GetNbinsX();
1503 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
eea641de 1504 // Loop over TOF bins
5f2ca1e6 1505 for(Int_t iBin=1; iBin<=nBinsTOF; iBin++){
1506 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1507 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1508 Float_t lowPtTOF=hSpectra[kTOF][ipart][icharge]->GetBinLowEdge(iBin);
1509 Float_t binWidTOF=hSpectra[kTOF][ipart][icharge]->GetBinWidth(iBin);
eea641de 1510 // Loop over TPC bins and find overlapping bins to ITS
5f2ca1e6 1511 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1512 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1513 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1514 if(TMath::Abs(lowPtTOF-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidTOF)<0.001){
1515 Float_t numer=hSpectra[kTOF][ipart][icharge]->GetBinContent(iBin);
1516 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1517 Float_t enumer=hSpectra[kTOF][ipart][icharge]->GetBinError(iBin);
78033eeb 1518 Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
5f2ca1e6 1519 Double_t ratio=0.;
1520 Double_t eratio=0.;
1521 if(numer>0. && denom>0.){
1522 ratio=numer/denom;
1523 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1524 }
1525 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1526 hRatioTOFTPC[ipart][icharge]->SetBinError(iBin,eratio);
1527 break;
1528 }
1529 }
1530 }
eea641de 1531 c1t->cd(ipart+1);
2210b002 1532 hRatioTOFTPC[ipart][icharge]->SetYTitle("TOF/TPC");
0e6f7392 1533 hRatioTOFTPC[ipart][icharge]->SetStats(1);
1534 hRatioTOFTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1535 hRatioTOFTPC[ipart][icharge]->Draw("");
1536 hRatioTOFTPC[ipart][icharge]->Fit("pol0","","same");
5f2ca1e6 1537 }
22abbb27 1538 PrintCanvas(c1t,printFormats);
5f2ca1e6 1539 }
1540
4c0d7bc7 1541}
1542
1543void DrawWithJacek() {
1544
1545 //1. Convert spectra to dNdeta and sum
1546 TH1F * hsum = (TH1F*) htemplate->Clone();
1547 hsum->Reset();
ee08b77d 1548 Int_t idet= iCombInStudy;
4c0d7bc7 1549 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1550 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1551 TH1 * h = hSpectra[idet][ipart][icharge];
1552 Int_t nbin = h->GetNbinsX();
1553 for(Int_t ibin = 1; ibin <= nbin; ibin++){
1554 Double_t pt = h->GetBinCenter(ibin);
1555 Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
1556 Double_t jacobian = pt/mt;
1557 h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
1558 h->SetBinError (ibin,h->GetBinError (ibin)*jacobian);
1559 Int_t ibinSum = hsum->FindBin(pt);
1560 Double_t epsilon = 0.0001;
1561 if ( h->GetBinContent(ibin) > 0 &&
1562 (TMath::Abs(h->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1563 TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1564 ) {
1565 cout << "DISCREPANCY IN BIN RANGES" << endl;
1566 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
1567 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
1568 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
1569 cout << "" << endl;
1570 }
1571
1572
1573 hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
1574 hsum->SetBinError (ibinSum,0);
1575 }
1576 // hsum->Add(h);
1577 }
1578 }
1579
1580
1581 // Load Jacek and Draw both:
77233d47 1582 new TFile ("./Files/pt_allcharged_published.root");
1583 TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("dndpt");
1584 gJacek->SetMarkerStyle(20);
1585 gJacek->SetMarkerColor(kGray);
1586 gJacek->Draw("AP");
1587 gJacek->GetHistogram()->SetXTitle("p_{t} (GeV)");
1588 gJacek->GetHistogram()->SetYTitle("d^{2}N/dp_{t}d#eta (GeV^{-1})");
1589 hsum->Draw("same");
1590
1591 TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
1592 new TCanvas();
1593 gRatio->Draw("AP");
4c0d7bc7 1594
1595
1596
1597}
1598
1599
1600void DrawRatioToStar() {
1601
1602 // Star data
1603 // gROOT->LoadMacro("StarPPSpectra.C");
1604 TGraphErrors ** gStar = StarPPSpectra();
c4fe358b 1605 gStar[0]->SetMarkerStyle(kOpenStar);
1606 gStar[1]->SetMarkerStyle(kFullStar);
1607 gStar[2]->SetMarkerStyle(kOpenStar);
1608 gStar[3]->SetMarkerStyle(kFullStar);
1609 gStar[4]->SetMarkerStyle(kOpenStar);
1610 gStar[5]->SetMarkerStyle(kFullStar);
4c0d7bc7 1611
1612 // ALICE, INEL -> NSD
1613 Double_t scaleYield = 3.58/3.02; // from paper 2
1614 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1615 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
ee08b77d 1616 hSpectra[iCombInStudy][ipart][icharge]->Scale(scaleYield);
4c0d7bc7 1617 }
1618 }
1619
1620
1621 TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
1622 TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
1623 hempty->SetXTitle("p_{t} (GeV/c)");
1624 hempty->SetYTitle("ALICE/STAR (NSD)");
c4fe358b 1625 hempty->Draw();
1626
1627 TCanvas * c1Comp = new TCanvas("cCompToStarNeg","cCompToStarNeg");
1628 c1Comp->SetLogy();
1629 TH2F * hempty2 = new TH2F(TString("hemptyCompNeg"),"hemptyCompNeg",100,0.,1.5, 100, 0.001,10);
1630 hempty2->SetXTitle("p_{t} (GeV/c)");
1631 hempty2->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1} [NSD]");
1632 hempty2->Draw();
4c0d7bc7 1633
1634 TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
1635 leg->SetBorderSize(0);
1636 leg->SetLineColor(1);
1637 leg->SetLineStyle(1);
1638 leg->SetLineWidth(1);
1639 leg->SetFillColor(0);
1640 leg->SetFillStyle(1001);
1641
1642
c4fe358b 1643
1644 c1->cd();
4c0d7bc7 1645 TGraphErrors * g ;
ee08b77d 1646 g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[iCombInStudy][kPion][kNeg],1);
4c0d7bc7 1647 g->SetMarkerStyle(kFullCircle);
1648 g->SetMarkerColor(kBlack);
1649 g->Draw("p");
1650 leg->AddEntry(g,"#pi^{-}","lp");
ee08b77d 1651 g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[iCombInStudy][kKaon][kNeg],1);
4c0d7bc7 1652 g->SetMarkerStyle(kOpenCircle);
1653 g->SetMarkerColor(kRed);
1654 g->Draw("p");
1655 leg->AddEntry(g,"K^{-}","lp");
ee08b77d 1656 g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[iCombInStudy][kProton][kNeg],1);
4c0d7bc7 1657 g->SetMarkerStyle(kOpenSquare);
1658 g->SetMarkerColor(kBlue);
1659 g->Draw("p");
1660 leg->AddEntry(g,"#bar{p}","lp");
1661 leg->Draw();
1662
c4fe358b 1663 c1Comp->cd();
1664 gStar[0]->Draw("p");
1665 hSpectra[iCombInStudy][kPion][kNeg]->Draw("same");
1666 gStar[2]->Draw("p");
1667 hSpectra[iCombInStudy][kKaon][kNeg]->Draw("same");
1668 gStar[4]->Draw("p");
1669 hSpectra[iCombInStudy][kProton][kNeg]->Draw("same");
4c0d7bc7 1670
1671
1672
1673 TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
1674 hempty->Draw();
1675 leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
1676 leg->SetBorderSize(0);
1677 leg->SetLineColor(1);
1678 leg->SetLineStyle(1);
1679 leg->SetLineWidth(1);
1680 leg->SetFillColor(0);
1681 leg->SetFillStyle(1001);
c4fe358b 1682
1683 TCanvas * c2Comp = new TCanvas("cCompToStarPos","cCompToStarPos");
1684 c2Comp->SetLogy();
1685 hempty2->Draw();
1686
1687 c2->cd();
4c0d7bc7 1688 // TGraphErrors * g ;
ee08b77d 1689 g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[iCombInStudy][kPion][kPos],1);
4c0d7bc7 1690 g->SetMarkerStyle(kFullCircle);
1691 g->SetMarkerColor(kBlack);
1692 g->Draw("p");
1693 leg->AddEntry(g,"#pi^{+}","lp");
ee08b77d 1694 g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[iCombInStudy][kKaon][kPos],1);
4c0d7bc7 1695 g->SetMarkerStyle(kOpenCircle);
1696 g->SetMarkerColor(kRed);
1697 g->Draw("p");
1698 leg->AddEntry(g,"K^{+}","lp");
ee08b77d 1699 g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[iCombInStudy][kProton][kPos],1);
4c0d7bc7 1700 g->SetMarkerStyle(kOpenSquare);
1701 g->SetMarkerColor(kBlue);
1702 g->Draw("p");
1703 leg->AddEntry(g,"p","lp");
c4fe358b 1704 leg->Draw();
4c0d7bc7 1705
1706
c4fe358b 1707 c2Comp->cd();
1708 gStar[1]->Draw("p");
1709 hSpectra[iCombInStudy][kPion][kPos]->Draw("same");
1710 gStar[3]->Draw("p");
1711 hSpectra[iCombInStudy][kKaon][kPos]->Draw("same");
1712 gStar[5]->Draw("p");
1713 hSpectra[iCombInStudy][kProton][kPos]->Draw("same");
1714
4c0d7bc7 1715
1716 c1->Update();
1717 c2->Update();
1718 gSystem->ProcessEvents();
1719 c1->Print(TString(c1->GetName()) + ".eps");
1720 c2->Print(TString(c2->GetName()) + ".eps");
1721 ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1722 ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1723 c1->Update();
1724 c2->Update();
1725 c1->Print(TString(c1->GetName()) + ".png");
1726 c2->Print(TString(c2->GetName()) + ".png");
1727
1728
1729
1730
1731}
1732
1733
1734
1735void DrawRatios() {
1736
1737 // Draws ratios of combined spectra
1738
4d4a268c 1739 // Rebin histograms (hi pt was fluctuating too much in the TOF)
1740 RebinHistosForRatios();
1741
4c0d7bc7 1742 // Compute ratios
1743 TH1F * hPosNegRatio[kNPart];
1744
1745 for(Int_t ipart = 0; ipart < kNPart; ipart++){
ee08b77d 1746 hPosNegRatio[ipart] = (TH1F*) hSpectra[iCombInStudy][ipart][kPos]->Clone();
1747 hPosNegRatio[ipart]->Divide(hSpectra[iCombInStudy][ipart][kNeg]);
4c0d7bc7 1748 hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
1749 hPosNegRatio[ipart]->SetMinimum(0.5);
1750 hPosNegRatio[ipart]->SetMaximum(1.5);
1751 }
1752
4d4a268c 1753
ee08b77d 1754 TH1F * hKPiRatio = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1755 hKPiRatio->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1756 TH1F * htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1757 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
4c0d7bc7 1758 hKPiRatio->Divide(htmp);
1759 hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
4d4a268c 1760 // Get ratio of syst errors
1761 // FIXME: commented combined error
1762 TH1F * hKPiRatioSyst = new TH1F(*htemplate);;
1763 AliBWTools::GetValueAndError(hKPiRatioSyst,hSpectra[iCombInStudy][kKaon][kPos],hSystError[iCombInStudy][kKaon][kPos],kTRUE);
1764 // AliBWTools::GetHistoCombinedErrors(hKPiRatioSyst,hSpectra[iCombInStudy][kKaon][kPos]);
1765 AliBWTools::GetValueAndError(htmp,hSpectra[iCombInStudy][kKaon][kNeg],hSystError[iCombInStudy][kKaon][kNeg],kTRUE);
1766 // // AliBWTools::GetHistoCombinedErrors(htmp,hSpectra[iCombInStudy][kKaon][kNeg]);
1767 hKPiRatioSyst->Add(htmp);
1768 AliBWTools::GetValueAndError(htmp,hSpectra[iCombInStudy][kPion][kNeg],hSystError[iCombInStudy][kPion][kNeg],kTRUE);
1769 // AliBWTools::GetHistoCombinedErrors(htmp,hSpectra[iCombInStudy][kPion][kNeg]);
1770 TH1F * htmp2 = new TH1F(*htemplate);
1771 AliBWTools::GetValueAndError(htmp2,hSpectra[iCombInStudy][kPion][kPos],hSystError[iCombInStudy][kPion][kPos],kTRUE);
1772 // AliBWTools::GetHistoCombinedErrors(htmp2,hSpectra[iCombInStudy][kPion][kPos]);
1773 htmp->Add(htmp2);
1774 hKPiRatioSyst->Divide(htmp);
1775 hKPiRatioSyst->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1776 // Multiply error for sqrt(2): the syst is not reduced by combining positive and negative charges
1777 Int_t nbin = hKPiRatioSyst->GetNbinsX();
1778 for(Int_t ibin = 1; ibin <= nbin; ibin++){
1779 hKPiRatioSyst->SetBinError(ibin, hKPiRatioSyst->GetBinError(ibin)*TMath::Sqrt(2));
1780 }
1781
1782
1783
4c0d7bc7 1784
1785 TH1F * hKPiRatioMC[kNTunes];
5f2ca1e6 1786 if(showMC){
1787 for(Int_t itune = 0; itune < kNTunes; itune++){
1788 hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
1789 hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
2210b002 1790 TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1791 htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1792 hKPiRatioMC[itune]->Divide(htmpLoc);
5f2ca1e6 1793 hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1794 }
4c0d7bc7 1795 }
1796
1797
1798
ee08b77d 1799 TH1F * hPPiRatio = (TH1F*) hSpectra[iCombInStudy][kProton][kPos]->Clone();
1800 hPPiRatio->Add(hSpectra[iCombInStudy][kProton][kNeg]);
1801 htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1802 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
4c0d7bc7 1803 hPPiRatio->Divide(htmp);
1804 hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1805
55056a39 1806 TH1F * hPPiRatioMC[kNTunes];
5f2ca1e6 1807 if(showMC){
5f2ca1e6 1808 for(Int_t itune = 0; itune < kNTunes; itune++){
1809 hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
1810 hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
2210b002 1811 TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1812 htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1813 hPPiRatioMC[itune]->Divide(htmpLoc);
5f2ca1e6 1814 hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1815 }
4c0d7bc7 1816 }
1817
d8e6677d 1818 // Phenix ratio
1819 TFile * fPhenix = new TFile("./Files/phenixdata.root");
1820 TH1F * hPPiRatioPhenix = (TH1F*) fPhenix->Get("ProtonPHNX");
1821 hPPiRatioPhenix->Add((TH1F*) fPhenix->Get("PbarPHNX"));
1822 htmp = (TH1F*) fPhenix->Get("PiPlusPHNX");
1823 htmp->Add((TH1F*) fPhenix->Get("PiMinusPHNX"));
1824 hPPiRatioPhenix = AliBWTools::DivideHistosDifferentBins(hPPiRatioPhenix,htmp);
1825 hPPiRatioPhenix->SetMarkerStyle(24);
1826
4c0d7bc7 1827 // Draw
1828// TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
1829// hempty->SetXTitle("p_{t} (GeV/c)");
1830 // hempty->SetYTitle("");
1831
1832 // tmp: overlay levi fits
1833 AliBWFunc * fm2 = new AliBWFunc;
1834 fm2->SetVarType(AliBWFunc::kdNdpt);
1835 TF1 * fLevi[kNPart] [kNCharge];
1836 fLevi[kPion] [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
1837 fLevi[kKaon] [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
1838 fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
1839 fLevi[kPion] [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
1840 fLevi[kKaon] [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
1841 fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
1842 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1843 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1844 fLevi[ipart][icharge]->SetRange(0,4);
1845 }
1846 }
1847
1848
1849 for(Int_t ipart = 0; ipart < kNPart; ipart++){
b7018bc7 1850 TString detName = detFlag[iCombInStudy];
1851 detName.ReplaceAll(" ", "_");
1852 detName.ReplaceAll("+", "");
1853
5f2ca1e6 1854 TCanvas * c1 = new TCanvas(TString("cRatio_")+detName+TString("_")+partFlag[ipart], TString("cRatio_")+detName+partFlag[ipart]);
1855 c1->SetGridy();
4c0d7bc7 1856 hPosNegRatio[ipart]->Draw();
1857 TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
1858 // fRatio->Draw("same");
1859 fRatio->SetRange(0,5);
4d4a268c 1860 PrintCanvas(c1,printFormats);
4c0d7bc7 1861
1862 }
1863
1864
1865 TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
d8e6677d 1866 // c2->SetGridy();
4d4a268c 1867
1868 hKPiRatioSyst->SetMaximum(0.62);
1869 hKPiRatioSyst->SetFillColor(kGray);
1870 hKPiRatioSyst->SetMarkerStyle(24);
1871 // hKPiRatioSyst->Draw("e2");
1872 // TMP FIXME: ERROR FROM FP:
1873 TFile * f = TFile::Open("./Files/ITSsa-systematics-20101014.root");
1874 TH1F * hsystRatio = (TH1F*) f->Get("hSystRatioKpi");
1875 TH1F * hsyst = new TH1F(*htemplate);
1876 AliBWTools::GetValueAndError(hsyst,hKPiRatio,hsystRatio,1);
1877 // AliBWTools::GetHistoCombinedErrors(hsyst,hKPiRatio);
1878 // hsyst->Draw("samee1");
1879 // END TMP
1880 hKPiRatio->SetMaximum(0.62);
1881 hKPiRatio->Draw("");
4c0d7bc7 1882 TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1883 lMC->SetFillColor(kWhite);
4d4a268c 1884 TGraphErrors ** gStar = StarPPSpectra();
1885 cout << gStar[6] << endl;
1886
1887 gStar[6]->SetMarkerStyle(kOpenStar);
1888 gStar[6]->Draw("P");
1889
5f2ca1e6 1890 if(showE735){
1891 gROOT->LoadMacro("GetE735Ratios.C");
4d4a268c 1892 // GetE735Ratios(0,0)->Draw("EX0,same");
1893 // GetE735Ratios(0,1)->Draw("EX0,same");
1894 // GetE735Ratios(0,2)->Draw("EX0,same");
5f2ca1e6 1895 GetE735Ratios(0,3)->Draw("EX0,same");
1896 }
4c0d7bc7 1897 hKPiRatio->SetMarkerStyle(20);
d8e6677d 1898 hKPiRatio->SetMarkerColor(kRed);
1899 hKPiRatio->SetLineColor(kRed);
4d4a268c 1900 //hKPiRatio->Draw("same");
4c0d7bc7 1901
5f2ca1e6 1902 if(showMC){
1903 for(Int_t itune = 0; itune < kNTunes; itune++){
1904 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1905 hKPiRatioMC[itune]->SetLineWidth(2);
1906 hKPiRatioMC[itune]->Draw("same,chist");
1907 }
4c0d7bc7 1908
5f2ca1e6 1909 lMC->Draw();
1910 }
4d4a268c 1911 cout << "12" << endl;
4c0d7bc7 1912
5f2ca1e6 1913 if(showE735){
1914 TLegend * l = new TLegend( 0.1879, 0.68, 0.54,0.92);
1915 l->SetFillColor(kWhite);
1916 l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
4d4a268c 1917 l->AddEntry(gStar[6], "STAR, #sqrt{s} = 200 GeV", "lp");
1918 // l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
1919 // l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
1920 // l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
5f2ca1e6 1921 l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
1922 l->Draw();
1923 }
4c0d7bc7 1924
1925
1926 TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
4d4a268c 1927 cout << "DRAWING" << endl;
d8e6677d 1928 // c3->SetGridy();
1929 hPPiRatio->SetMarkerStyle(20);
1930 hPPiRatio->SetMarkerColor(kRed);
1931 hPPiRatio->SetLineColor(kRed);
4c0d7bc7 1932 hPPiRatio->Draw();
1933 hPPiRatio->SetMaximum(0.39);
5f2ca1e6 1934 if(showMC){
1935 lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1936 lMC->SetFillColor(kWhite);
1937
1938 for(Int_t itune = 0; itune < kNTunes; itune++){
1939 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
55056a39 1940 hPPiRatioMC[itune]->SetLineWidth(2);
1941 hPPiRatioMC[itune]->Draw("same,chist");
5f2ca1e6 1942 }
4c0d7bc7 1943
5f2ca1e6 1944 lMC->Draw();
1945 }
d8e6677d 1946 hPPiRatioPhenix->Draw("same");
1947
1948
1949 TLegend * l = new TLegend(0.186242, 0.781469, 0.538591, 0.921329);
1950 l->SetFillColor(kWhite);
1951 l->AddEntry(hPPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1952 l->AddEntry(hPPiRatioPhenix, "PHENIX, #sqrt{s} = 200 GeV","lpf");
1953 // TLegend * l = new TLegend(0.206376, 0.77972, 0.600671, 0.909091);
1954 l->Draw();
1955
4d4a268c 1956 PrintCanvas(c2,printFormats);
1957 PrintCanvas(c3,printFormats);
4c0d7bc7 1958
1959}
1960
3725266c 1961void DrawWithSyst() {
1962
1963 // Draws detector and combined with syst errors.
1964
1965 for(Int_t idet = 0; idet < kNHist; idet++){
1966 if(idet > kITS && idet != iCombInStudy) continue;
1967 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1968 TCanvas * c = new TCanvas(TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge],TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge]);
1969 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
1970 hempty->SetXTitle("p_{t} (GeV/c)");
1971 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1972 hempty->GetYaxis()->SetTitleOffset(1.35);
1973 hempty->GetXaxis()->SetTitleOffset(1.1);
1974 hempty->Draw();
1975 c->SetLogy();
1976 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1977 // cout << detFlag[idet] << " " << chargeFlag[icharge] << " " << partFlag[ipart] << endl;
1978
1979 TString opt = ipart ? "" : "same";
1980 TH1F * hsyst = new TH1F(*htemplate);
1981 AliBWTools::GetValueAndError(hsyst,hSpectra[idet][ipart][icharge],hSystError[idet][ipart][icharge],kTRUE);
55056a39 1982 hsyst->SetFillColor(kGray);
3725266c 1983 hsyst->Draw("e5,same");
1984 hSpectra[idet][ipart][icharge]->Draw("same");
22abbb27 1985 // Do not draw syst error outside of spectra range
1986 Float_t lowEdge = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1987 Float_t hiEdge = AliBWTools::GetHighestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1988 Int_t lowBin = hSystError[idet][ipart][icharge]->FindBin(lowEdge);
1989 Int_t hiBin = hSystError[idet][ipart][icharge]->FindBin(hiEdge);
1990 Int_t nbin = hSystError[idet][ipart][icharge]->GetNbinsX();
2210b002 1991 for(Int_t ibin = 1; ibin < lowBin; ibin++){
22abbb27 1992 hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1993 hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1994 }
1995 for(Int_t ibin = hiBin; ibin <= nbin; ibin++){
1996 hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1997 hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1998 }
1999 // Draw relative syst error, convert to graph to avoid vertical bars at the edges (empty bins in histo)
2000 AliBWTools::GetGraphFromHisto(hSystError[idet][ipart][icharge])->Draw("LX");
2001
3725266c 2002
2003
2004 }
22abbb27 2005 PrintCanvas(c,"png,eps");
3725266c 2006 }
2007 }
2008}
c4fe358b 2009
2010void Help() {
2011
2012 cout << "Macro: void CombineSpectra(Int_t analysisType=kDoFits, Int_t fitFuncID = kFitLevi) " << endl;
2013 cout << "" << endl;
4d4a268c 2014 cout << "WARNING: some function modify the histos upon execution. You are advised to quit root" << endl;
2015 cout << " before running a different function" << endl;
2016 cout << "" << endl;
c4fe358b 2017 cout << "Possible Arguments" << endl;
2018 cout << "- analysisType:" << endl;
eea641de 2019 cout << " kDoFits: Fit Combined Spectra " << endl;
2020 cout << " kDoRatios: Particle ratios K/pi and p/pi" << endl;
2021 cout << " kDoSuperposition: Compare different detectors (superimpose and ratios)" << endl;
2022 cout << " kDoCompareStar: Compare combined spectra to star results" << endl;
2023 cout << " kDoDrawWithModels: Compare combined spectra and models" << endl;
3725266c 2024 cout << " kDoDrawSyst: Draws spectra from individual detectors with their systematic error" << endl;
77233d47 2025 cout << " kDoCompareToAllCharged: convert to dn/deta/dpt, sum up and compare to all charged" << endl;
eea641de 2026 cout << " kDoHelp: This help" << endl;
2027 cout << "- fitFuncID, function used to extrapolate and compute yields" << endl;
2028 cout << " An analitic fit function [kFitLevi, kFitUA1, kFitPowerLaw]" << endl;
2029 cout << " Or a shape from a MC moder [kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0]" << endl;
2f7f947e 2030 cout << " Which is fitted to the data at low pt and used to extrapolate at low pt" << endl;
c4fe358b 2031
2032
2033}
3725266c 2034
2035void PrintCanvas(TCanvas* c,const TString formats) {
2036 // print a canvas in every of the given comma-separated formats
2037 // ensure the canvas is updated
22abbb27 2038 if(!doPrint) return;
3725266c 2039 c->Update();
2040 gSystem->ProcessEvents();
2041 TObjArray * arr = formats.Tokenize(",");
2042 TIterator * iter = arr->MakeIterator();
2043 TObjString * element = 0;
2044 TString name =c ->GetName();
2045 name.ReplaceAll(" ","_");
2046 name.ReplaceAll("+","Plus");
2047 name.ReplaceAll("-","");
2048 while ((element = (TObjString*) iter->Next())) {
2049 c->Print(name+ "."+element->GetString());
2050 }
2051}
2052
22abbb27 2053
3725266c 2054void LoadLibs(){
2055
2056 gSystem->Load("libTree.so");
2057 gSystem->Load("libVMC.so");
2058 gSystem->Load("libMinuit.so");
2059 gSystem->Load("libSTEERBase.so");
2060 gSystem->Load("libESD.so");
2061 gSystem->Load("libAOD.so");
2062 gSystem->Load("libANALYSIS.so");
2063 gSystem->Load("libANALYSISalice.so");
2064 gSystem->Load("libCORRFW.so");
2065 gSystem->Load("libPWG2spectra.so");
2066
2067}
22abbb27 2068
4d4a268c 2069void RebinHistosForRatios() {
2070 // rebins all spectra histograms
2071
2072 cout << "WARNING: REBINNING HISTOS BEFORE RATIOS" << endl;
2073
2074
2075 const Float_t rebinBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.6,1.8,2.4};
2076 Int_t nrebinBins=29;
2077
2078
2079
2080 Int_t idet = iCombInStudy;
2081 for(Int_t ipart = 0; ipart < kNPart; ipart++){
2082 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
2083
2084 TH1F * h = hSpectra[idet][ipart][icharge];
2085 TString hname = h->GetName();
2086 hname += "_rebin";
2087
2088 TH1F *hrebin = new TH1F(hname, hname,nrebinBins,rebinBins);
2089 hrebin->Sumw2();
2090
2091 for(Int_t i=1;i <= h->GetNbinsX();i++){
2092 Float_t x = h->GetBinCenter(i);
2093 Int_t j=0;
2094 while(x > rebinBins[j+1]) j++;
2095 j++;
2096 hrebin->AddBinContent(j,h->GetBinContent(i));
2097
2098 Float_t err = hrebin->GetBinError(j) + h->GetBinError(i)*h->GetBinError(i);
2099 hrebin->SetBinError(j,err);
2100 }
2101
2102 for(Int_t i=1;i <= hrebin->GetNbinsX();i++){
2103 hrebin->SetBinError(i,TMath::Sqrt(hrebin->GetBinError(i)));
2104 }
2105 hSpectra[idet][ipart][icharge] = hrebin;
2106 }
2107 }
2108
2109}