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