coverity fix and some clean up
[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / AverageAndExtrapolate.C
1 #if !defined (__CINT__) || (defined(__MAKECINT__))
2 #include <iostream>
3 #include "TClonesArray.h"
4 #include "TH1.h"
5 #include "AliPWGFunc.h"
6 //#include "LoadSpectraPiKPPbPb.C"
7 #include "AliParticleYield.h"
8 #include "TFile.h"
9 #include "TDatabasePDG.h"
10
11 #endif
12
13
14 void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale = 1. , Int_t isLinear = 0) ;
15 TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc = 0, Double_t massParticle=0.139) ;
16 TH1* ReweightSpectra(TClonesArray * histos, Double_t * weights, Int_t isLinear = 0, TString suffix = "") ;
17 void LoadStuff(Bool_t loadPbPb = kTRUE) ;
18 void AverageAndExtrapolate(TString what);
19
20
21
22 TH1 * hLambdaStat[7];
23 TH1 * hLambdaSyst[7];
24 TH1 * hK0SStat[7]   ;
25 TH1 * hK0SSyst[7]   ;
26
27
28
29 void AverageAndExtrapolate () {
30
31   // MF COMMON LIBRARIES SHOULD NOT BE LOADED FOR THIS MACRO TO WORK
32   LoadStuff(1); // True PbPb, False pPb
33
34   //  PrintYieldAndError(hSpectraCentr_sys[0][kMyProton][0], hSpectraCentr_stat[0][kMyProton][0], 0, mass[2]);
35
36   //  PrintYieldAndError(hLambdaSyst[0], hLambdaStat[0], 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
37   // icentr, ipart, icharge
38
39   //  AverageAndExtrapolate("pions_pos_0020");
40   //  AverageAndExtrapolate("pions_pos_0010");
41   // AverageAndExtrapolate("pions_neg_0010");
42   //  AverageAndExtrapolate("pions_sum_0010");
43   // AverageAndExtrapolate("kaons_pos_0010");
44   // AverageAndExtrapolate("kaons_neg_0010");
45   // AverageAndExtrapolate("protons_pos_0010");
46   // AverageAndExtrapolate("protons_neg_0010");
47   //  AverageAndExtrapolate("lambda_0010");
48   //  AverageAndExtrapolate("k0s_0010");
49   // AverageAndExtrapolate("pions_pos_6080");
50   // AverageAndExtrapolate("pions_neg_6080");
51   // AverageAndExtrapolate("kaons_pos_6080");
52   // AverageAndExtrapolate("kaons_neg_6080");
53   // AverageAndExtrapolate("protons_pos_6080");
54   // AverageAndExtrapolate("protons_neg_6080");
55
56   //  AverageAndExtrapolate("pions_pos_2040");
57   AverageAndExtrapolate("pions_neg_2040");
58   AverageAndExtrapolate("kaons_pos_2040");
59   AverageAndExtrapolate("kaons_neg_2040");
60   AverageAndExtrapolate("protons_pos_2040");
61   AverageAndExtrapolate("protons_neg_2040");
62   
63 }
64
65 TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc, Double_t massParticle) {
66
67   TF1 * f = 0;
68   if (ifunc == 0) {
69     f= BGBlastWave("fBlastWave",  massParticle);
70     f->SetParameters(massParticle, 0.640, 0.097, 0.73, 100); // FIXME
71   }
72   else if (ifunc == 1) {
73     f  = fm.GetTsallis(massParticle, 0.2, 11, 800, "fTsallisPion");
74   }
75   
76   //  TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100, 0.01, 0.1, "");
77   TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100);
78   //  std::cout << "H " << h << std::endl;
79   std::cout << "" << std::endl;
80   std::cout << h->GetBinContent(1) << ", " << h->GetBinContent(2) << ", " << h->GetBinContent(3) << std::endl;
81   std::cout << "" << std::endl;
82   return h;
83
84
85 }
86
87
88 TH1* ReweightSpectra(TObjArray * histos, Double_t * weights, Int_t isLinear, TString suffix) {
89
90   // sums a number of spectra with a given weight. Used to combine
91   // several centrality bins.  The weights should add up to 1 and be
92   // proportional to the width of the centrality bin for each
93   // histogram
94   // Suffix is added to the name of the histogram.
95   // if linear = 1 errors are summed linearly rather than in quadrature
96
97
98   TIter it(histos);
99   TH1 * h = 0;
100   TH1 * hsum = 0;
101   Int_t ihisto = 0;
102   while ((h = dynamic_cast<TH1*>(it.Next()))) {
103     if(!h) {
104       std::cout << "ERROR cannot get one of the histos!" << std::endl;
105       return 0;
106     }
107
108     if (!hsum) {
109       // First histogram, clone it
110       hsum =  (TH1D*) h->Clone(TString(h->GetName())+suffix);
111       hsum->Scale(weights[ihisto++]);
112     }else {
113       AddHistograms(hsum, h, weights[ihisto++], isLinear);
114     }
115   }
116   return hsum;
117 }
118
119 void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale, Int_t isLinear) {
120   // THis method assumes that the 2 histos are consistent!!
121   TH1 * hsourceLoc = (TH1*) hsource->Clone("hsourceLoc");
122   hsourceLoc->Scale(scale);
123   if(!isLinear) {
124     hdest->Add(hsourceLoc);
125   } 
126   else {
127     Int_t nbin = hsourceLoc->GetNbinsX();
128     for(Int_t ibin = 0; ibin < nbin; ibin++){
129       Double_t content = hdest->GetBinContent(ibin) + hsourceLoc->GetBinContent(ibin);
130       Double_t error   = hdest->GetBinError  (ibin) + hsourceLoc->GetBinError  (ibin);
131       hdest->SetBinContent(ibin,content);
132       hdest->SetBinError(ibin, error);
133     }
134     
135
136   }
137
138   delete hsourceLoc;
139   
140   
141
142 }
143
144
145 void LoadStuff(Bool_t loadPbPb) {
146   if(loadPbPb){
147     gROOT->LoadMacro("LoadSpectraPiKPPbPb.C");
148     LoadSpectraPiKPPbPb();
149   }
150   else {
151     gROOT->LoadMacro("LoadSpectraPiKPProtonLead.C");
152     LoadSpectraPiKPProtonLead();
153   }
154   LoadLibs();
155   gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/YieldMean.C");
156   gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/SpectraUtils.C");
157
158   // Load Lambdas and K0s
159   TFile * f = new TFile("k0s_lambda_final_spectra.root");
160   const char * multTags[] = {  "0005",  "0510",  "1020",  "2040",  "4060",  "6080", "8090"};
161   for (Int_t icentr = 0; icentr<7; icentr++) {
162     hLambdaStat[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
163     hLambdaSyst[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
164     hK0SStat[icentr]    = (TH1*) f->Get(Form("systonly_cent%s_K0s",multTags[icentr]));
165     hK0SSyst[icentr]    = (TH1*) f->Get(Form("statonly_cent%s_K0s",multTags[icentr]));
166
167     // The bin 300-400 MeV was not used in the analysis
168     hK0SStat[icentr]->SetBinContent(4,0);   
169     hK0SSyst[icentr]->SetBinContent(4,0);
170     hK0SStat[icentr]->SetBinError(4,0);   
171     hK0SSyst[icentr]->SetBinError(4,0);
172   }
173
174 }
175
176 void AverageAndExtrapolate(TString what) {
177
178   TH1 * hstat =0;
179   TH1 * hsyst =0;
180   TH1 * hyieldmean =0;
181
182   if(what == "pions_pos_0020"){ 
183
184     TObjArray * arrStat = new TObjArray();
185     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
186     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
187     arrStat->Add(hSpectraCentr_stat[2][kMyPion][kMyPos]);
188     
189     TObjArray * arrSyst = new TObjArray();
190     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
191     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
192     arrSyst->Add(hSpectraCentr_sys [2][kMyPion][kMyPos]);
193     
194     Double_t weights[] = {0.25, 0.25, 0.5};
195     //Double_t weights[] = {0.5, 0.5};
196     
197     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0020");
198     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0020");
199
200     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
201
202     
203   }
204   if(what == "pions_pos_0010"){ 
205
206     TObjArray * arrStat = new TObjArray();
207     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
208     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
209     
210     TObjArray * arrSyst = new TObjArray();
211     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
212     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
213     
214     Double_t weights[] = {0.5, 0.5};
215     
216     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
217     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
218
219     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
220
221     
222   }
223   if(what == "pions_sum_0010"){ 
224
225     TObjArray * arrStat = new TObjArray();
226     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
227     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
228     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
229     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
230     
231     TObjArray * arrSyst = new TObjArray();
232     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
233     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
234     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
235     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
236     
237     Double_t weights[] = {0.5, 0.5, 0.5, 0.5};
238     
239     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
240     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
241
242     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
243
244     
245   }
246   if(what == "pions_neg_0010"){ 
247
248     TObjArray * arrStat = new TObjArray();
249     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
250     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
251     
252     TObjArray * arrSyst = new TObjArray();
253     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
254     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
255     
256     Double_t weights[] = {0.5, 0.5};
257     
258     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
259     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
260
261     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
262
263     
264   }
265   if(what == "protons_pos_0010"){ 
266
267     TObjArray * arrStat = new TObjArray();
268     arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyPos]);
269     arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyPos]);
270     
271     TObjArray * arrSyst = new TObjArray();
272     arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyPos]);
273     arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyPos]);
274     
275     Double_t weights[] = {0.5, 0.5};
276     
277     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
278     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
279
280     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
281
282     
283   }
284   if(what == "protons_neg_0010"){ 
285
286     TObjArray * arrStat = new TObjArray();
287     arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyNeg]);
288     arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyNeg]);
289     
290     TObjArray * arrSyst = new TObjArray();
291     arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyNeg]);
292     arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyNeg]);
293     
294     Double_t weights[] = {0.5, 0.5};
295     
296     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
297     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
298
299     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
300
301     
302   }
303     if(what == "kaons_pos_0010"){ 
304
305     TObjArray * arrStat = new TObjArray();
306     arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyPos]);
307     arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyPos]);
308     
309     TObjArray * arrSyst = new TObjArray();
310     arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyPos]);
311     arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyPos]);
312     
313     Double_t weights[] = {0.5, 0.5};
314     
315     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
316     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
317
318     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
319
320     
321   }
322   if(what == "kaons_neg_0010"){ 
323
324     TObjArray * arrStat = new TObjArray();
325     arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyNeg]);
326     arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyNeg]);
327     
328     TObjArray * arrSyst = new TObjArray();
329     arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyNeg]);
330     arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyNeg]);
331     
332     Double_t weights[] = {0.5, 0.5};
333     
334     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
335     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
336
337     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
338
339     
340   }
341
342   if(what == "pions_pos_6080"){ 
343
344     TObjArray * arrStat = new TObjArray();
345     arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyPos]);
346     arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyPos]);
347     
348     TObjArray * arrSyst = new TObjArray();
349     arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyPos]);
350     arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyPos]);
351     
352     Double_t weights[] = {0.5, 0.5};
353     
354     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
355     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
356
357     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
358
359     
360   }
361   if(what == "pions_neg_6080"){ 
362
363     TObjArray * arrStat = new TObjArray();
364     arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyNeg]);
365     arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyNeg]);
366     
367     TObjArray * arrSyst = new TObjArray();
368     arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyNeg]);
369     arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyNeg]);
370     
371     Double_t weights[] = {0.5, 0.5};
372     
373     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
374     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
375
376     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
377
378     
379   }
380   if(what == "protons_pos_6080"){ 
381
382     TObjArray * arrStat = new TObjArray();
383     arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyPos]);
384     arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyPos]);
385     
386     TObjArray * arrSyst = new TObjArray();
387     arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyPos]);
388     arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyPos]);
389     
390     Double_t weights[] = {0.5, 0.5};
391     
392     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
393     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
394
395     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
396
397     
398   }
399   if(what == "protons_neg_6080"){ 
400
401     TObjArray * arrStat = new TObjArray();
402     arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyNeg]);
403     arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyNeg]);
404     
405     TObjArray * arrSyst = new TObjArray();
406     arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyNeg]);
407     arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyNeg]);
408     
409     Double_t weights[] = {0.5, 0.5};
410     
411     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
412     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
413
414     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
415
416     
417   }
418     if(what == "kaons_pos_6080"){ 
419
420     TObjArray * arrStat = new TObjArray();
421     arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyPos]);
422     arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyPos]);
423     
424     TObjArray * arrSyst = new TObjArray();
425     arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyPos]);
426     arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyPos]);
427     
428     Double_t weights[] = {0.5, 0.5};
429     
430     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
431     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
432
433     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
434
435     
436   }
437   if(what == "kaons_neg_6080"){ 
438
439     TObjArray * arrStat = new TObjArray();
440     arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyNeg]);
441     arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyNeg]);
442     
443     TObjArray * arrSyst = new TObjArray();
444     arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyNeg]);
445     arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyNeg]);
446     
447     Double_t weights[] = {0.5, 0.5};
448     
449     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
450     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
451
452     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
453
454     
455   }
456
457   if(what == "pions_pos_2040"){ 
458
459     TObjArray * arrStat = new TObjArray();
460     arrStat->Add(hSpectraCentr_stat[3][kMyPion][kMyPos]);
461     arrStat->Add(hSpectraCentr_stat[4][kMyPion][kMyPos]);
462     
463     TObjArray * arrSyst = new TObjArray();
464     arrSyst->Add(hSpectraCentr_sys [3][kMyPion][kMyPos]);
465     arrSyst->Add(hSpectraCentr_sys [4][kMyPion][kMyPos]);
466     
467     Double_t weights[] = {0.5, 0.5};
468     
469     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
470     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
471
472     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
473
474     
475   }
476   if(what == "pions_neg_2040"){ 
477
478     TObjArray * arrStat = new TObjArray();
479     arrStat->Add(hSpectraCentr_stat[3][kMyPion][kMyNeg]);
480     arrStat->Add(hSpectraCentr_stat[4][kMyPion][kMyNeg]);
481     
482     TObjArray * arrSyst = new TObjArray();
483     arrSyst->Add(hSpectraCentr_sys [3][kMyPion][kMyNeg]);
484     arrSyst->Add(hSpectraCentr_sys [4][kMyPion][kMyNeg]);
485     
486     Double_t weights[] = {0.5, 0.5};
487     
488     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
489     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
490
491     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
492
493     
494   }
495     if(what == "kaons_pos_2040"){ 
496
497     TObjArray * arrStat = new TObjArray();
498     arrStat->Add(hSpectraCentr_stat[3][kMyKaon][kMyPos]);
499     arrStat->Add(hSpectraCentr_stat[4][kMyKaon][kMyPos]);
500     
501     TObjArray * arrSyst = new TObjArray();
502     arrSyst->Add(hSpectraCentr_sys [3][kMyKaon][kMyPos]);
503     arrSyst->Add(hSpectraCentr_sys [4][kMyKaon][kMyPos]);
504     
505     Double_t weights[] = {0.5, 0.5};
506     
507     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
508     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
509
510     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
511
512     
513   }
514   if(what == "kaons_neg_2040"){ 
515
516     TObjArray * arrStat = new TObjArray();
517     arrStat->Add(hSpectraCentr_stat[3][kMyKaon][kMyNeg]);
518     arrStat->Add(hSpectraCentr_stat[4][kMyKaon][kMyNeg]);
519     
520     TObjArray * arrSyst = new TObjArray();
521     arrSyst->Add(hSpectraCentr_sys [3][kMyKaon][kMyNeg]);
522     arrSyst->Add(hSpectraCentr_sys [4][kMyKaon][kMyNeg]);
523     
524     Double_t weights[] = {0.5, 0.5};
525     
526     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
527     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
528
529     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
530
531     
532   }
533
534   if(what == "protons_pos_2040"){ 
535
536     TObjArray * arrStat = new TObjArray();
537     arrStat->Add(hSpectraCentr_stat[3][kMyProton][kMyPos]);
538     arrStat->Add(hSpectraCentr_stat[4][kMyProton][kMyPos]);
539     
540     TObjArray * arrSyst = new TObjArray();
541     arrSyst->Add(hSpectraCentr_sys [3][kMyProton][kMyPos]);
542     arrSyst->Add(hSpectraCentr_sys [4][kMyProton][kMyPos]);
543     
544     Double_t weights[] = {0.5, 0.5};
545     
546     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
547     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
548
549     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
550
551     
552   }
553   if(what == "protons_neg_2040"){ 
554
555     TObjArray * arrStat = new TObjArray();
556     arrStat->Add(hSpectraCentr_stat[3][kMyProton][kMyNeg]);
557     arrStat->Add(hSpectraCentr_stat[4][kMyProton][kMyNeg]);
558     
559     TObjArray * arrSyst = new TObjArray();
560     arrSyst->Add(hSpectraCentr_sys [3][kMyProton][kMyNeg]);
561     arrSyst->Add(hSpectraCentr_sys [4][kMyProton][kMyNeg]);
562     
563     Double_t weights[] = {0.5, 0.5};
564     
565     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040");
566     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040");
567
568     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
569
570     
571   }
572
573
574   if(what == "lambda_0010"){ 
575
576     TObjArray * arrStat = new TObjArray();
577     arrStat->Add(hLambdaStat[0]);
578     arrStat->Add(hLambdaStat[1]);
579     
580     TObjArray * arrSyst = new TObjArray();
581     arrSyst->Add(hLambdaSyst[0]);
582     arrSyst->Add(hLambdaSyst[1]);
583     
584     Double_t weights[] = {0.5, 0.5};
585     
586     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
587     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
588
589     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
590
591     
592   }
593
594   if(what == "k0s_0010"){ 
595
596     TObjArray * arrStat = new TObjArray();
597     arrStat->Add(hK0SStat[0]);
598     arrStat->Add(hK0SStat[1]);
599     
600     TObjArray * arrSyst = new TObjArray();
601     arrSyst->Add(hK0SSyst[0]);
602     arrSyst->Add(hK0SSyst[1]);
603     
604     Double_t weights[] = {0.5, 0.5};
605     
606     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
607     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
608
609     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("K_S0")->Mass());
610
611     
612   }
613
614
615
616   
617   TCanvas * c1 = new TCanvas("Averaging", "Averaging");
618   c1->Divide(2,1);
619   c1->cd(1);
620   hstat->Draw();
621   hsyst->Draw("same,e3");
622   c1->cd(2);
623   hyieldmean->Draw();
624   
625
626 }
627
628