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