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   LoadStuff(0); // True PbPb, False pPb
  32
  33   //  PrintYieldAndError(hSpectraCentr_sys[0][kMyProton][0], hSpectraCentr_stat[0][kMyProton][0], 0, mass[2]);
  34
  35   //  PrintYieldAndError(hLambdaSyst[0], hLambdaStat[0], 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
  36   // icentr, ipart, icharge
  37
  38   //  AverageAndExtrapolate("pions_pos_0020");
  39   //  AverageAndExtrapolate("pions_pos_0010");
  40   // AverageAndExtrapolate("pions_neg_0010");
  41   AverageAndExtrapolate("pions_sum_0010");
  42   // AverageAndExtrapolate("kaons_pos_0010");
  43   // AverageAndExtrapolate("kaons_neg_0010");
  44   // AverageAndExtrapolate("protons_pos_0010");
  45   // AverageAndExtrapolate("protons_neg_0010");
  46   //  AverageAndExtrapolate("lambda_0010");
  47   //  AverageAndExtrapolate("k0s_0010");
  48   // AverageAndExtrapolate("pions_pos_6080");
  49   // AverageAndExtrapolate("pions_neg_6080");
  50   // AverageAndExtrapolate("kaons_pos_6080");
  51   // AverageAndExtrapolate("kaons_neg_6080");
  52   // AverageAndExtrapolate("protons_pos_6080");
  53   // AverageAndExtrapolate("protons_neg_6080");
  54
  55 }
  56
  57 TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc, Double_t massParticle) {
  58
  59   TF1 * f = 0;
  60   if (ifunc == 0) {
  61     f= BGBlastWave("fBlastWave",  massParticle);
  62     f->SetParameters(massParticle, 0.640, 0.097, 0.73, 100); // FIXME
  63   }
  64   else if (ifunc == 1) {
  65     f  = fm.GetTsallis(massParticle, 0.2, 11, 800, "fTsallisPion");
  66   }
  67
  68   //  TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100, 0.01, 0.1, "");
  69   TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100);
  70   //  std::cout << "H " << h << std::endl;
  71   std::cout << "" << std::endl;
  72   std::cout << h->GetBinContent(1) << ", " << h->GetBinContent(2) << ", " << h->GetBinContent(3) << std::endl;
  73   std::cout << "" << std::endl;
  74   return h;
  75
  76
  77 }
  78
  79
  80 TH1* ReweightSpectra(TObjArray * histos, Double_t * weights, Int_t isLinear, TString suffix) {
  81
  82   // sums a number of spectra with a given weight. Used to combine
  83   // several centrality bins.  The weights should add up to 1 and be
  84   // proportional to the width of the centrality bin for each
  85   // histogram
  86   // Suffix is added to the name of the histogram.
  87   // if linear = 1 errors are summed linearly rather than in quadrature
  88
  89
  90   TIter it(histos);
  91   TH1 * h = 0;
  92   TH1 * hsum = 0;
  93   Int_t ihisto = 0;
  94   while ((h = dynamic_cast<TH1*>(it.Next()))) {
  95     if(!h) {
  96       std::cout << "ERROR cannot get one of the histos!" << std::endl;
  97       return 0;
  98     }
  99
 100     if (!hsum) {
 101       // First histogram, clone it
 102       hsum =  (TH1D*) h->Clone(TString(h->GetName())+suffix);
 103       hsum->Scale(weights[ihisto++]);
 104     }else {
 105       AddHistograms(hsum, h, weights[ihisto++], isLinear);
 106     }
 107   }
 108   return hsum;
 109 }
 110
 111 void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale, Int_t isLinear) {
 112   // THis method assumes that the 2 histos are consistent!!
 113   TH1 * hsourceLoc = (TH1*) hsource->Clone("hsourceLoc");
 114   hsourceLoc->Scale(scale);
 115   if(!isLinear) {
 116     hdest->Add(hsourceLoc);
 117   }
 118   else {
 119     Int_t nbin = hsourceLoc->GetNbinsX();
 120     for(Int_t ibin = 0; ibin < nbin; ibin++){
 121       Double_t content = hdest->GetBinContent(ibin) + hsourceLoc->GetBinContent(ibin);
 122       Double_t error   = hdest->GetBinError  (ibin) + hsourceLoc->GetBinError  (ibin);
 123       hdest->SetBinContent(ibin,content);
 124       hdest->SetBinError(ibin, error);
 125     }
 126
 127
 128   }
 129
 130   delete hsourceLoc;
 131
 132
 133
 134 }
 135
 136
 137 void LoadStuff(Bool_t loadPbPb) {
 138   if(loadPbPb){
 139     gROOT->LoadMacro("LoadSpectraPiKPPbPb.C");
 140     LoadSpectraPiKPPbPb();
 141   }
 142   else {
 143     gROOT->LoadMacro("LoadSpectraPiKPProtonLead.C");
 144     LoadSpectraPiKPProtonLead();
 145   }
 146   LoadLibs();
 147   gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/YieldMean.C");
 148   gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/SpectraUtils.C");
 149
 150   // Load Lambdas and K0s
 151   TFile * f = new TFile("k0s_lambda_final_spectra.root");
 152   const char * multTags[] = {  "0005",  "0510",  "1020",  "2040",  "4060",  "6080", "8090"};
 153   for (Int_t icentr = 0; icentr<7; icentr++) {
 154     hLambdaStat[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
 155     hLambdaSyst[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
 156     hK0SStat[icentr]    = (TH1*) f->Get(Form("systonly_cent%s_K0s",multTags[icentr]));
 157     hK0SSyst[icentr]    = (TH1*) f->Get(Form("statonly_cent%s_K0s",multTags[icentr]));
 158
 159     // The bin 300-400 MeV was not used in the analysis
 160     hK0SStat[icentr]->SetBinContent(4,0);
 161     hK0SSyst[icentr]->SetBinContent(4,0);
 162     hK0SStat[icentr]->SetBinError(4,0);
 163     hK0SSyst[icentr]->SetBinError(4,0);
 164   }
 165
 166 }
 167
 168 void AverageAndExtrapolate(TString what) {
 169
 170   TH1 * hstat =0;
 171   TH1 * hsyst =0;
 172   TH1 * hyieldmean =0;
 173
 174   if(what == "pions_pos_0020"){
 175
 176     TObjArray * arrStat = new TObjArray();
 177     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
 178     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
 179     arrStat->Add(hSpectraCentr_stat[2][kMyPion][kMyPos]);
 180
 181     TObjArray * arrSyst = new TObjArray();
 182     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
 183     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
 184     arrSyst->Add(hSpectraCentr_sys [2][kMyPion][kMyPos]);
 185
 186     Double_t weights[] = {0.25, 0.25, 0.5};
 187     //Double_t weights[] = {0.5, 0.5};
 188
 189     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0020");
 190     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0020");
 191
 192     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 193
 194
 195   }
 196   if(what == "pions_pos_0010"){
 197
 198     TObjArray * arrStat = new TObjArray();
 199     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
 200     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
 201
 202     TObjArray * arrSyst = new TObjArray();
 203     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
 204     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
 205
 206     Double_t weights[] = {0.5, 0.5};
 207
 208     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 209     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 210
 211     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 212
 213
 214   }
 215   if(what == "pions_sum_0010"){
 216
 217     TObjArray * arrStat = new TObjArray();
 218     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
 219     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
 220     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
 221     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
 222
 223     TObjArray * arrSyst = new TObjArray();
 224     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
 225     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
 226     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
 227     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
 228
 229     Double_t weights[] = {0.5, 0.5, 0.5, 0.5};
 230
 231     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 232     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 233
 234     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 235
 236
 237   }
 238   if(what == "pions_neg_0010"){
 239
 240     TObjArray * arrStat = new TObjArray();
 241     arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
 242     arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
 243
 244     TObjArray * arrSyst = new TObjArray();
 245     arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
 246     arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
 247
 248     Double_t weights[] = {0.5, 0.5};
 249
 250     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 251     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 252
 253     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 254
 255
 256   }
 257   if(what == "protons_pos_0010"){
 258
 259     TObjArray * arrStat = new TObjArray();
 260     arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyPos]);
 261     arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyPos]);
 262
 263     TObjArray * arrSyst = new TObjArray();
 264     arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyPos]);
 265     arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyPos]);
 266
 267     Double_t weights[] = {0.5, 0.5};
 268
 269     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 270     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 271
 272     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
 273
 274
 275   }
 276   if(what == "protons_neg_0010"){
 277
 278     TObjArray * arrStat = new TObjArray();
 279     arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyNeg]);
 280     arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyNeg]);
 281
 282     TObjArray * arrSyst = new TObjArray();
 283     arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyNeg]);
 284     arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyNeg]);
 285
 286     Double_t weights[] = {0.5, 0.5};
 287
 288     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 289     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 290
 291     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
 292
 293
 294   }
 295     if(what == "kaons_pos_0010"){
 296
 297     TObjArray * arrStat = new TObjArray();
 298     arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyPos]);
 299     arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyPos]);
 300
 301     TObjArray * arrSyst = new TObjArray();
 302     arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyPos]);
 303     arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyPos]);
 304
 305     Double_t weights[] = {0.5, 0.5};
 306
 307     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 308     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 309
 310     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
 311
 312
 313   }
 314   if(what == "kaons_neg_0010"){
 315
 316     TObjArray * arrStat = new TObjArray();
 317     arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyNeg]);
 318     arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyNeg]);
 319
 320     TObjArray * arrSyst = new TObjArray();
 321     arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyNeg]);
 322     arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyNeg]);
 323
 324     Double_t weights[] = {0.5, 0.5};
 325
 326     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 327     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 328
 329     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
 330
 331
 332   }
 333
 334   if(what == "pions_pos_6080"){
 335
 336     TObjArray * arrStat = new TObjArray();
 337     arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyPos]);
 338     arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyPos]);
 339
 340     TObjArray * arrSyst = new TObjArray();
 341     arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyPos]);
 342     arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyPos]);
 343
 344     Double_t weights[] = {0.5, 0.5};
 345
 346     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 347     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 348
 349     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 350
 351
 352   }
 353   if(what == "pions_neg_6080"){
 354
 355     TObjArray * arrStat = new TObjArray();
 356     arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyNeg]);
 357     arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyNeg]);
 358
 359     TObjArray * arrSyst = new TObjArray();
 360     arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyNeg]);
 361     arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyNeg]);
 362
 363     Double_t weights[] = {0.5, 0.5};
 364
 365     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 366     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 367
 368     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
 369
 370
 371   }
 372   if(what == "protons_pos_6080"){
 373
 374     TObjArray * arrStat = new TObjArray();
 375     arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyPos]);
 376     arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyPos]);
 377
 378     TObjArray * arrSyst = new TObjArray();
 379     arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyPos]);
 380     arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyPos]);
 381
 382     Double_t weights[] = {0.5, 0.5};
 383
 384     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 385     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 386
 387     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
 388
 389
 390   }
 391   if(what == "protons_neg_6080"){
 392
 393     TObjArray * arrStat = new TObjArray();
 394     arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyNeg]);
 395     arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyNeg]);
 396
 397     TObjArray * arrSyst = new TObjArray();
 398     arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyNeg]);
 399     arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyNeg]);
 400
 401     Double_t weights[] = {0.5, 0.5};
 402
 403     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 404     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 405
 406     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
 407
 408
 409   }
 410     if(what == "kaons_pos_6080"){
 411
 412     TObjArray * arrStat = new TObjArray();
 413     arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyPos]);
 414     arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyPos]);
 415
 416     TObjArray * arrSyst = new TObjArray();
 417     arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyPos]);
 418     arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyPos]);
 419
 420     Double_t weights[] = {0.5, 0.5};
 421
 422     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 423     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 424
 425     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
 426
 427
 428   }
 429   if(what == "kaons_neg_6080"){
 430
 431     TObjArray * arrStat = new TObjArray();
 432     arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyNeg]);
 433     arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyNeg]);
 434
 435     TObjArray * arrSyst = new TObjArray();
 436     arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyNeg]);
 437     arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyNeg]);
 438
 439     Double_t weights[] = {0.5, 0.5};
 440
 441     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
 442     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
 443
 444     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
 445
 446
 447   }
 448
 449
 450   if(what == "lambda_0010"){
 451
 452     TObjArray * arrStat = new TObjArray();
 453     arrStat->Add(hLambdaStat[0]);
 454     arrStat->Add(hLambdaStat[1]);
 455
 456     TObjArray * arrSyst = new TObjArray();
 457     arrSyst->Add(hLambdaSyst[0]);
 458     arrSyst->Add(hLambdaSyst[1]);
 459
 460     Double_t weights[] = {0.5, 0.5};
 461
 462     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 463     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 464
 465     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
 466
 467
 468   }
 469
 470   if(what == "k0s_0010"){
 471
 472     TObjArray * arrStat = new TObjArray();
 473     arrStat->Add(hK0SStat[0]);
 474     arrStat->Add(hK0SStat[1]);
 475
 476     TObjArray * arrSyst = new TObjArray();
 477     arrSyst->Add(hK0SSyst[0]);
 478     arrSyst->Add(hK0SSyst[1]);
 479
 480     Double_t weights[] = {0.5, 0.5};
 481
 482     TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
 483     TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
 484
 485     hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("K_S0")->Mass());
 486
 487
 488   }
 489
 490
 491
 492
 493   TCanvas * c1 = new TCanvas("Averaging", "Averaging");
 494   c1->Divide(2,1);
 495   c1->cd(1);
 496   hstat->Draw();
 497   hsyst->Draw("same,e3");
 498   c1->cd(2);
 499   hyieldmean->Draw();
 500
 501
 502 }
 503
 504