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