PWGLF/ThermalFits/PlotRatiosForQM14.C

   1 #if !defined (__CINT__) || (defined(__MAKECINT__))
   2 #include <iostream>
   3 #include "TClonesArray.h"
   4 #include "AliParticleYield.h"
   5 #include "TH1F.h"
   6 #include "TCanvas.h"
   7 #include "TStyle.h"
   8 #include <fstream>
   9 #include "TLatex.h"
  10 #include "TLegend.h"
  11 #include "TList.h"
  12 #include "TF1.h"
  13 #include "AliPWGHistoTools.h"
  14 #include "TGraphErrors.h"
  15 #include "TMath.h"
  16 #include "TDatabasePDG.h"
  17
  18
  19
  20 #endif
  21
  22 // Plots ratios for QM and saves input files for thermal models
  23
  24
  25 enum MyParticles { kPDGPi = 211, kPDGK = 321, kPDGProton = 2212, kPDGKS0 = 310, kPDGLambda=3122, kPDGXi=3312,kPDGOmega=3334,kPDGPhi=333,kPDGKStar=313,kPDGDeuteron=1000010020,kPDGHE3 = 1000020030, kPDGHyperTriton = 1010010030, kPDGSigmaStarPlus=3224,kPDGSigmaStarMinus=3114,kPDGXiStar=3324};
  26
  27 typedef enum {kStatError, kSystError, kTotalError} myerror_t;
  28
  29 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle, Int_t icolor, Int_t imarker = kOpenSquare, Int_t errorsType = kTotalError, Float_t shift = 0) ;
  30 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ;
  31 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
  32 void SetStyle(Bool_t graypalette=0) ;
  33 void NewLegendQM() ;
  34 void DrawRatio(TString what);
  35 void LoadArrays() ;
  36
  37 // Preferred colors and markers
  38 // const Int_t fillColors[] = {kGray+1,  kRed-10, kBlue-9, kGreen-8, kMagenta-9, kOrange-9,kCyan-8,kYellow-7, kWhite}; // for syst bands
  39 // const Int_t colors[]     = {kBlack, kRed+1 , kBlue+1, kGreen+3, kMagenta+1, kOrange-1,kCyan+2,kYellow+2  , kWhite};
  40 const Int_t markers[]    = {kFullCircle, kFullSquare,kOpenCircle,kOpenSquare,kOpenDiamond,kOpenCross,kFullCross,kFullDiamond,kFullStar,kOpenStar,0};
  41
  42 Double_t maxy = 0.4;
  43
  44 // Data arrays;
  45 TClonesArray *arrPbPb=0, *arrpp7=0, *arrpPb=0, * arrpp276=0, * arrpp900=0, * arrThermus=0;
  46 TClonesArray *arrSTARPbPb=0, *arrPHENIXPbPb=0, *arrBRAHMSPbPb=0;
  47 TClonesArray *arrSTARpp  =0, *arrPHENIXpp=0;
  48
  49 const Double_t *scaleRatios = 0;
  50 Double_t *correlatedUnc = 0;
  51 Double_t correlatedUncLocalPbPb[14] = {0.0424 , 0.0424     ,  0.041     ,  0      , 0         , 0.0424       , 0.0424       , 0               , 0.05    , 0.05   };
  52 Double_t correlatedUncLocalPP  [14] = {0.0424 , 0.0424     ,  0.041     ,  0      , 0         , 0.0424       , 0.0424       , 0               , 0.0424    , 0.0424   };
  53 Double_t correlatedUncZero[14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  54
  55
  56 TClonesArray * PlotRatiosForQM14() {
  57 #if !(!defined (__CINT__) || (defined(__MAKECINT__)))
  58   LoadLibs();
  59 #endif
  60
  61   //
  62   LoadArrays();
  63
  64   // Uncomment stuff in this section to save the inputs for thermal models
  65   //#define SAVE_INPUT_THERMAL_MODEL
  66 #ifdef SAVE_INPUT_THERMAL_MODEL
  67   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
  68   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
  69   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0);
  70   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/1);
  71
  72
  73   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", /*separateCharges*/1);
  74   //  PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A2040", /*separateCharges*/1);
  75   //PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A6080", /*separateCharges*/1);
  76   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M6080", /*separateCharges*/1);
  77   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M2030", /*separateCharges*/1);
  78
  79   return 0;
  80 #endif
  81
  82   SetStyle();
  83
  84   TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
  85   c1->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.83);
  86   //  c1->SetLogy();
  87
  88   // CENTRAL
  89   DrawRatio("frame");
  90
  91   DrawRatio("PbPb_0010");
  92   // DrawRatio("PbPbSTAR");
  93   // DrawRatio("PbPbPHENIX");
  94   // DrawRatio("PbPbBRAHMS");
  95   //  DrawRatio("PbPb_6080");
  96
  97   DrawRatio("pp7");
  98   DrawRatio("pPb0005");
  99   DrawRatio("pp276");
 100   DrawRatio("pp900");
 101   // DrawRatio("ppSTAR");
 102   // DrawRatio("ppPHENIX");
 103   //  DrawRatio("ppBRAHMS");
 104
 105   NewLegendQM();
 106
 107   return 0;
 108   //
 109   TCanvas * c2 = new TCanvas("Yields", "Yields", 1400, 600);
 110   c2->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368);
 111
 112   GetHistoYields(arrPbPb,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%")->Draw();
 113   GetHistoYields(arrThermus,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
 114   NewLegendQM();
 115   return arrPbPb;
 116 }
 117
 118 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle, Int_t icolor, Int_t imarker, Int_t errorType, Float_t shift) {
 119   // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY
 120
 121   const Int_t nratio = 10;
 122   //  Int_t denum[nratio]    = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGPi       , kPDGDeuteron , kPDGPi          , kPDGK   , -kPDGK};
 123
 124   Int_t num  [nratio]            = {kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
 125   Int_t denum[nratio]            = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGProton   , kPDGDeuteron , kPDGPi          , kPDGK   , kPDGK};
 126   Int_t isSum[nratio]            = {1      , 1          ,  1         ,   1     , 1         , 0            , 0            , 1               , 1       , 1      };
 127   static const Double_t scale[]  = {1      , 3          ,  0.5       ,  30     ,  250      , 50           , 100          , 4e5             , 2       , 2      };
 128   //  static const Double_t scale[]  = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,};
 129
 130   scaleRatios = scale;
 131   TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s_%d",system,energy,centrality.Data(),errorType), histotitle, nratio, 1+shift, nratio+1+shift);
 132
 133   TClonesArray * arrRatios = new TClonesArray ("AliParticleYield");// We save to disk the precomputed ratios
 134   Int_t iratioArr = 0;// Index used only to add particles to the array above
 135
 136   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
 137   for(Int_t iratio = 1; iratio <= nratio; iratio++){
 138     std::cout << "******** " << num[iratio-1] << "/" <<  denum[iratio-1]<< " ("<<isSum[iratio-1]<<")*******" << std::endl ;
 139
 140     AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
 141     if(ratio)
 142       {
 143         ratio = new AliParticleYield(*ratio); // We need to clone it to avoid a mess if we need to use this particle again later (e.g. double scalings)
 144         std::cout << "  " ;
 145         ratio->Print("short");
 146       }
 147
 148     if(!ratio) {
 149       // If the ratio is not found, try to build it!
 150       std::cout << "  Looking for " <<  num[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
 151       AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,  isSum[iratio-1]);
 152       if(part1) part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
 153       // Try with the !sum, if part 1 is not found
 154       if(!part1) {
 155         std::cout << "  Looking for " <<  num[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
 156         part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 157         if(part1)
 158           {
 159             part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
 160             // If the sum was requested, try to recover it!
 161             if(isSum[iratio-1]) {
 162               std::cout << "  Looking for " <<  -num[iratio-1] <<std::endl;
 163               AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0);
 164               if(part1 && part1_bar) {
 165                 std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl;
 166                 part1 = AliParticleYield::Add(part1, part1_bar);
 167
 168               } else if (TDatabasePDG::Instance()->GetParticle(-num[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
 169                 std::cout << "WARNING: Sum requested but not found, scaling x2 " << part1->GetName() << std::endl;
 170                 part1->Scale(2);
 171               }
 172             } else if(part1) { // if we are here, it means the sum was *not* requested (isSum=0), but we found something with (!isSum) = 1
 173               // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 174               std::cout << "WARNING: Using sum/2 for " << part1->GetName() << std::endl;
 175
 176               part1->Scale(0.5);
 177             }
 178         }
 179
 180       }
 181       std::cout << "  Looking for " <<  denum[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
 182       AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
 183       if(part2) part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
 184       if(!part2) {// Try with the !sum, if part 2 is not found
 185         std::cout << "  Looking for " <<  denum[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
 186         part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 187         if(part2)
 188           {
 189             part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
 190             if(isSum[iratio-1]) {
 191               std::cout << "  Looking for " <<  -denum[iratio-1] << std::endl;
 192               AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0);
 193               if(part2 && part2_bar){
 194                 std::cout << "Adding " << part2_bar->GetPartName() << " " << part2->GetPartName() << std::endl;
 195                 part2 = AliParticleYield::Add(part2, part2_bar);
 196               } else if (TDatabasePDG::Instance()->GetParticle(-denum[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
 197                 std::cout << "WARNING: Sum requested but not found, scaling x2 " << part2->GetName() << std::endl;
 198                 part2->Scale(2);
 199               }
 200             } else if(part2){
 201               // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 202               std::cout << "WARNING: Using sum/2 for " << part2->GetName() << std::endl;
 203               part2->Scale(0.5);
 204             }
 205           }
 206
 207       }
 208       ratio = AliParticleYield::Divide(part1, part2, correlatedUnc[iratio-1], "YQ"); // Assume by that the systematics of part1 and part2 are uncorrelated.
 209       if(ratio) {
 210         std::cout << "" << std::endl;
 211         std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
 212         ratio->Print("");
 213         std::cout << "" << std::endl;
 214       }
 215     }
 216     if(ratio){
 217       ratio->Scale(scale[iratio-1]);
 218       h->SetBinContent(iratio, ratio->GetYield());
 219       if(errorType == kTotalError) {
 220         h->SetBinError  (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
 221       } else if (errorType == kStatError) {
 222         h->SetBinError  (iratio, ratio->GetStatError());
 223       } else if (errorType == kSystError) {
 224         h->SetBinError  (iratio, ratio->GetSystError());
 225       } else {
 226         std::cout <<  "ERROR: Unknown Error Type " << errorType << std::endl;
 227       }
 228
 229       //      h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f",  scale[iratio-1]), ratio->GetLatexName()));
 230       h->GetXaxis()->SetBinLabel(iratio, ratio->GetLatexName());
 231       new ((*arrRatios)[iratioArr++]) AliParticleYield(*ratio);
 232     }
 233     else {
 234       h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
 235
 236     }
 237     std::cout << "*** END OF " << num[iratio-1] << "/" <<  denum[iratio-1]<< " *******" << std::endl ;
 238
 239   }
 240
 241   h->GetYaxis()->SetRangeUser(0, maxy);
 242   h->SetLineColor(icolor);
 243   h->SetMarkerColor(icolor);
 244   h->SetMarkerStyle(imarker);
 245
 246   AliParticleYield::SaveAsASCIIFile(arrRatios, TString("ratios_")+h->GetName());
 247   return h;
 248
 249
 250
 251 }
 252
 253 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges)  {
 254   // If "Separate charges" is true, tries to dump both charges are dumped
 255   TClonesArray * arrOut = new TClonesArray("AliParticleYield");
 256   const Int_t npart = 12;
 257   Int_t particles  [npart] = {kPDGPi ,kPDGK   ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3    };
 258   Int_t isSum[npart]       = {1      ,1       ,0      , 1        , 0      , 1           ,0           ,1        ,1          ,0           , 1              , 0          };
 259
 260   Int_t ipartOut = 0; // Index for the array
 261   for(Int_t ipart = 0; ipart < npart; ipart++){
 262     if(!separateCharges) {
 263       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  isSum[ipart]);
 264       if(!part && isSum[ipart]) {
 265         //Could not find the particle, but the sum was requested: build the sum!
 266         part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 267         AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 268         if(part2 && part) part = AliParticleYield::Add(part, part2);
 269         else part = 0;
 270       }
 271       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 272     }
 273     else {
 274       // ignore isSum and try to find both particles
 275       Bool_t notFound = 0;
 276       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 277       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 278       else notFound=1;
 279       // Try to find antiparticle (-pdg code)
 280       part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 281       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 282       else if (notFound) {
 283         // If neither charge was found, check if we at least have the sum
 284         part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality,  1);
 285         if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 286       }
 287
 288     }
 289   }
 290   std::cout << "Particles for thermal model fits:" << std::endl;
 291   //  arrOut->Print("short");
 292   arrOut->Print("");
 293   std::cout << "" << std::endl;
 294   // Write GSI input file
 295   TIter it(arrOut);
 296   AliParticleYield * part = 0;
 297   ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 298   while ((part = (AliParticleYield*) it.Next())){
 299     fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
 300   }
 301   fout.close();
 302   // Write thermus file
 303   AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 304 }
 305
 306
 307 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) {
 308
 309   const Int_t npart = 11;
 310   Int_t pdg  [npart]    = {kPDGPi, kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
 311   //  Int_t isSum[npart]    = {1      ,1      ,1           ,0           ,1        ,1          ,0             ,0             ,1                ,0        ,1      };
 312   Int_t isSum[npart]    = {0,0,0,0,0,0,0,0,0,0,1};
 313   //  Double_t scale[npart] = {1      ,1      ,3           ,1           ,30       ,250         ,50            ,10            ,4e5              ,2       ,2      };
 314   Double_t scale[npart] = {1,5,30,30,200,1000,4000,2e6,2e6,20,20,};
 315   TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1, npart+1);
 316
 317   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
 318   for(Int_t ipart = 1; ipart <= npart; ipart++){
 319     AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg[ipart-1], system, energy, centrality,isSum[ipart-1]);
 320     if(!part) continue;
 321     part->Scale(scale[ipart-1]);
 322     h->SetBinContent(ipart, part->GetYield());
 323     h->SetBinError  (ipart, part->GetTotalError(0/* 0 = no normalization error */));
 324     h->GetXaxis()->SetBinLabel(ipart, Form("#splitline{%s}{%s}",Form("#times%2.0g",  scale[ipart-1]), part->GetLatexName()));
 325
 326   }
 327   return h;
 328 }
 329
 330 void LoadArrays() {
 331   arrPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
 332   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
 333   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
 334   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
 335   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
 336   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
 337   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
 338   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
 339
 340   arrpp7   = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
 341
 342   arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt");
 343   arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt");
 344
 345   arrpPb   = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt");
 346   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
 347   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"));
 348   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_AveragedNumbers.txt"));
 349
 350   arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt");
 351
 352   // RHIC data
 353   arrSTARPbPb   = AliParticleYield::ReadFromASCIIFile("PbPb_200_STAR-AntonQM12.txt");
 354   arrPHENIXPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_PHENIX-AntonQM12.txt");
 355   arrBRAHMSPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_BRAHMS-AntonQM12.txt");
 356   arrSTARpp     = AliParticleYield::ReadFromASCIIFile("pp_200_STAR.txt");
 357   arrPHENIXpp   = AliParticleYield::ReadFromASCIIFile("pp_200_PHENIX.txt");
 358
 359 }
 360
 361 void SetStyle(Bool_t graypalette) {
 362   std::cout << "Setting style!" << std::endl;
 363
 364   gStyle->Reset("Plain");
 365   gStyle->SetOptTitle(0);
 366   gStyle->SetOptStat(0);
 367   if(graypalette) gStyle->SetPalette(8,0);
 368   else gStyle->SetPalette(1);
 369   gStyle->SetCanvasColor(10);
 370   gStyle->SetCanvasBorderMode(0);
 371   gStyle->SetFrameLineWidth(1);
 372   gStyle->SetFrameFillColor(kWhite);
 373   gStyle->SetPadColor(10);
 374   gStyle->SetPadTickX(1);
 375   gStyle->SetPadTickY(1);
 376   gStyle->SetPadBottomMargin(0.15);
 377   gStyle->SetPadLeftMargin(0.15);
 378   gStyle->SetHistLineWidth(1);
 379   gStyle->SetHistLineColor(kRed);
 380   gStyle->SetFuncWidth(2);
 381   gStyle->SetFuncColor(kGreen);
 382   gStyle->SetLineWidth(2);
 383   gStyle->SetLabelSize(0.045,"yz");
 384   gStyle->SetLabelSize(0.06,"x");
 385   gStyle->SetLabelOffset(0.01,"y");
 386   gStyle->SetLabelOffset(0.01,"x");
 387   gStyle->SetLabelColor(kBlack,"xyz");
 388   gStyle->SetTitleSize(0.05,"xyz");
 389   gStyle->SetTitleOffset(1.25,"y");
 390   gStyle->SetTitleOffset(1.2,"x");
 391   gStyle->SetTitleFillColor(kWhite);
 392   gStyle->SetTextSizePixels(26);
 393   gStyle->SetTextFont(42);
 394   //  gStyle->SetTickLength(0.04,"X");  gStyle->SetTickLength(0.04,"Y");
 395
 396   gStyle->SetLegendBorderSize(0);
 397   gStyle->SetLegendFillColor(kWhite);
 398   //  gStyle->SetFillColor(kWhite);
 399   gStyle->SetLegendFont(42);
 400
 401   gStyle->SetErrorX(0);
 402   gStyle->SetEndErrorSize(5);
 403 }
 404
 405 void NewLegendQM() {
 406
 407   const char * style = "lp";
 408   const char ** labels=0;
 409   //  Bool_t beautify=kFALSE;
 410   Bool_t useTitle=kTRUE;
 411
 412   TLegend * l = new TLegend(  0.0985145, 0.733119, 0.301016, 0.869775);
 413   l->SetFillColor(kWhite);
 414
 415   // const Int_t markers[] = {20,24,21,25,23,28,33,20,24,21,25,23,28,33};
 416   // const Int_t colors[]  = {1,2,3,4,6,7,8,9,10,11,1,2,3,4,6,7,8,9,10};
 417
 418   TList * list = gPad->GetListOfPrimitives();
 419   TIterator * iter = list->MakeIterator();
 420   TObject * obj = 0;
 421   Int_t ilabel = -1;
 422   while ((obj = (TObject*) iter->Next())){
 423     if (obj->InheritsFrom("TH1") || obj->InheritsFrom("TGraph") || obj->InheritsFrom("TF1")) {
 424       if( (TString(obj->GetName()) == "hframe" ) ) continue;
 425       ilabel++;
 426       if (labels != 0)
 427         l->AddEntry(obj, labels[ilabel], style);
 428       else{
 429         if (useTitle)  {
 430           if(TString(obj->GetTitle()).Contains("NoLegend")) continue;
 431           l->AddEntry(obj, obj->GetTitle(), style);
 432         }
 433         else
 434           l->AddEntry(obj, obj->GetName(), style);
 435       }
 436       // if(beautify) {
 437
 438       //   if(!obj->InheritsFrom("TF1")){
 439       //     ((TH1*)obj)->SetLineColor(colors[ilabel]);
 440       //     ((TH1*)obj)->SetMarkerColor(colors[ilabel]);
 441       //     ((TH1*)obj)->SetMarkerStyle(markers[ilabel]);
 442       //   } else {
 443       //     ((TF1*)obj)->SetLineColor(colors[ilabel]);
 444       //   }
 445
 446     }
 447   }
 448
 449   l->Draw();
 450
 451 }
 452
 453
 454 void DrawRatio(TString what) {
 455   // This is used to simplify the code above
 456   // In order to draw syst error bars, we need to convert to graphs the syst errors histos
 457
 458   // Sample colors
 459   //  const Int_t colors[]     = {kBlack, kRed+1 , kBlue+1, kGreen+3, kMagenta+1, kOrange-1,kCyan+2,kYellow+2  , kWhite};
 460
 461   TClonesArray * array = 0;
 462   Int_t system,  color, marker;
 463   Float_t energy = 0, shift = 0;
 464   TString centrality, label;
 465   // FIXME: move this in the different sections below
 466   correlatedUnc = 0;
 467   std::cout << "Plotting " << what.Data() << std::endl;
 468
 469
 470   if (what == "frame" ) {
 471     // This is a bit of an hack: since the particle labels come
 472     // directly from AliPArticleYield, and since the PbPb sample is
 473     // the only one where we have all the ratios, we draw the PbPb
 474     // ratio here and then we set lines and markers to white. We also
 475     // add the "NoLegend" flag, so that it does not show up in the legend
 476     correlatedUnc = correlatedUncZero;
 477     TH1 * h = GetHistoRatios(arrPbPb,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "NoLegend", kWhite);
 478     h->Draw();
 479     Int_t nratio = h->GetNbinsX();
 480     for(Int_t iratio = 0; iratio < nratio; iratio++){
 481       Double_t exp = TMath::Floor(TMath::Log10(TMath::Abs(scaleRatios[iratio])));
 482       Double_t man = scaleRatios[iratio] / TMath::Power(10, exp);
 483       if(exp > 2) {
 484         TLatex * scaleLabel = new TLatex(iratio+1+0.2,maxy*1.01, Form("#times %0.0f 10^{%0.0f}", man, exp));
 485         scaleLabel->Draw();
 486       } else {
 487         TLatex * scaleLabel = new TLatex(iratio+1+0.2,maxy*1.01, Form("#times %g", scaleRatios[iratio]));
 488         scaleLabel->Draw();
 489       }
 490     }
 491
 492     TLatex *   tex = new TLatex(8.8,0.037,"ALICE Preliminary");
 493     tex->SetTextFont(42);
 494     tex->SetLineWidth(2);
 495     tex->Draw();
 496
 497     h->GetYaxis()->SetDecimals(1);
 498     h->GetYaxis()->SetNdivisions(505);
 499     //    h->GetXaxis()->CenterLabels(1);
 500     gPad->SetGridx();
 501
 502   }
 503   else if (what == "PbPb_0010") {
 504     array = arrPbPb;
 505     system = 2; energy = 2760.; centrality = "V0M0010";
 506     label = "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%";
 507     color  = kRed+1;
 508     marker = kFullCircle;
 509     shift =  0;
 510     correlatedUnc = correlatedUncLocalPbPb;
 511
 512   }
 513
 514   else if (what == "PbPb_6080") {
 515     array = arrPbPb;
 516     system = 2; energy = 2760.; centrality = "V0M6080";
 517     label = "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 60-80%";
 518     color  = kBlue+1;
 519     marker = kFullCircle;
 520     shift =  0.0;
 521     correlatedUnc = correlatedUncLocalPbPb;
 522   }
 523   else if (what == "pp7") {
 524     array = arrpp7;
 525     system = 0; energy = 7000.; centrality = "";
 526     label = "pp #sqrt{s} = 7 TeV";
 527     color  = kMagenta+1;
 528     marker = kFullCircle;
 529     shift =  0.2;
 530     correlatedUnc = correlatedUncLocalPP;
 531   }
 532   else if (what == "pp900") {
 533     array = arrpp900;
 534     system = 0; energy = 900.; centrality = "";
 535     label = "pp #sqrt{s} = 0.9 TeV";
 536     color  = kCyan+2;
 537     marker = kFullCircle;
 538     shift =  -0.2;
 539     correlatedUnc = correlatedUncLocalPP;
 540   }
 541   else if (what == "pp276") {
 542     array = arrpp276;
 543     system = 0; energy = 2760.; centrality = "";
 544     label = "pp #sqrt{s} = 2.76 TeV";
 545     color  = kYellow+2;
 546     marker = kFullCircle;
 547     shift = 0;
 548     correlatedUnc = correlatedUncLocalPP;
 549   }
 550   else if (what == "pPb0005") {
 551     array = arrpPb;
 552     system = 1; energy = 5020.; centrality = "V0A0005";
 553     label = "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%";
 554     color  = kBlack;
 555     marker = kFullCircle;
 556     shift = -0.2;
 557     correlatedUnc = correlatedUncLocalPP;
 558   }
 559   else if (what == "PbPbSTAR") {
 560     array = arrSTARPbPb;
 561     system = 2; energy = 200.; centrality = "0005";
 562     label = "STAR, Pb-Pb, #sqrt{s_{NN}} = 0.2 TeV, 0-5%";
 563     color  = kBlack;
 564     marker = kOpenStar;
 565     shift = +0.2;
 566     correlatedUnc = correlatedUncZero;
 567   }
 568   else if (what == "PbPbPHENIX") {
 569     array = arrPHENIXPbPb;
 570     system = 2; energy = 200.; centrality = "0005";
 571     label = "PHENIX, Pb-Pb, #sqrt{s_{NN}} = 0.2 TeV, 0-5%";
 572     color  = kBlack;
 573     marker = kOpenSquare;
 574     shift = -0.15;
 575     correlatedUnc = correlatedUncZero;
 576   }
 577   else if (what == "PbPbBRAHMS") {
 578     array = arrBRAHMSPbPb;
 579     system = 2; energy = 200.; centrality = "0010";
 580     label = "BRAHMS, Pb-Pb, #sqrt{s_{NN}} = 0.2 TeV, 0-10%";
 581     color  = kBlack;
 582     marker = kOpenCross;
 583     shift = -0.3;
 584     correlatedUnc = correlatedUncZero;
 585   }
 586   else if (what == "ppSTAR") {
 587     array = arrSTARpp;
 588     system = 0; energy = 200.; centrality = "";
 589     label = "STAR, pp, #sqrt{s} = 0.2 TeV";
 590     color  = kBlack;
 591     marker = kOpenStar;
 592     shift = 0.;
 593     correlatedUnc = correlatedUncZero;
 594   }
 595   else if (what == "ppPHENIX") {
 596     array = arrPHENIXpp;
 597     system = 0; energy = 200.; centrality = "";
 598     label = "PHENIX, pp, #sqrt{s} = 0.2 TeV";
 599     color  = kBlack;
 600     marker = kOpenSquare;
 601     shift = -0.2;
 602     correlatedUnc = correlatedUncZero;
 603   }
 604
 605
 606   else {
 607     std::cout << "Unknown ratio " << what.Data() << std::endl;
 608   }
 609
 610   if(!correlatedUnc) {
 611     std::cout << "correlatedUnc not set!" << std::endl;
 612
 613   }
 614   std::cout << "CORR: " << correlatedUnc[1] << std::endl;
 615
 616   if(array) {
 617     AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array,  system,  energy, centrality, label, color, marker, kStatError, shift)
 618                                         ,0)->Draw("PZ");
 619     AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array,  system,  energy, centrality, label+"NoLegend", color, marker, kSystError, shift)
 620                                         ,0)->Draw("[]");
 621   }
 622
 623
 624 }