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 <fstream>
   8 #endif
   9
  10 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};
  11
  12 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ;
  13 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ;
  14 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
  15
  16 // Plots ratios for QM and saves input files for thermal models
  17
  18 TClonesArray * PlotRatiosForQM14() {
  19 #if !(!defined (__CINT__) || (defined(__MAKECINT__)))
  20    LoadLibs();
  21 #endif
  22   TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
  23   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
  24   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
  25   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
  26   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
  27   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
  28   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
  29   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
  30
  31   TClonesArray * arrpp7   = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
  32
  33   TClonesArray * arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt");
  34   TClonesArray * arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt");
  35
  36   TClonesArray * arrpPb   = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt");
  37   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
  38   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"));
  39
  40   TClonesArray * arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt");
  41
  42   // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges0);
  43   // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges1);
  44   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0);
  45   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/1);
  46
  47   TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
  48   c1->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368);
  49
  50   //  c1->SetLogy();
  51   // CENTRAL
  52   TH1 * h =  GetHistoRatios(arr,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%");
  53   h->GetYaxis()->SetRangeUser(0, 0.4);
  54   h->Draw();
  55   // //GetHistoRatios(arrThermus,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
  56   // GetHistoRatios(arrpp7, AliParticleYield::kCSpp,   7000, ""       , "pp #sqrt{s} = 7 TeV"                   )->Draw("same");
  57   // GetHistoRatios(arrpPb,    AliParticleYield::kCSpPb,  5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same");
  58   // // GetHistoRatios(arrpp276, AliParticleYield::kCSpp,   2760, ""       , "pp #sqrt{s} = 2.76 TeV"                   )->Draw("same");
  59   // // GetHistoRatios(arrpp900, AliParticleYield::kCSpp,   900, ""       , "pp #sqrt{s} = 0.9 TeV"                   )->Draw("same");
  60   // NewLegend("", "lp",0,1,1);
  61
  62   // Peripheral
  63   GetHistoRatios(arr,       AliParticleYield::kCSPbPb, 2760, "V0M6080", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 60-80%")->Draw("same");
  64   //GetHistoRatios(arrThermus,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
  65   // GetHistoRatios(arrpp7, AliParticleYield::kCSpp,   7000, ""       , "pp #sqrt{s} = 7 TeV"                   )->Draw("same");
  66   GetHistoRatios(arrpPb,    AliParticleYield::kCSpPb,  5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same");
  67   // GetHistoRatios(arrpp276, AliParticleYield::kCSpp,   2760, ""       , "pp #sqrt{s} = 2.76 TeV"                   )->Draw("same");
  68   // GetHistoRatios(arrpp900, AliParticleYield::kCSpp,   900, ""       , "pp #sqrt{s} = 0.9 TeV"                   )->Draw("same");
  69   NewLegend("", "lp",0,1,1);
  70
  71   //return;
  72   TCanvas * c2 = new TCanvas("Yields", "Yields", 1400, 600);
  73   c2->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368);
  74
  75   GetHistoYields(arr,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%")->Draw();
  76   GetHistoYields(arrThermus,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
  77   NewLegend("", "lp",0,1,1);
  78   return arr;
  79 }
  80
  81 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) {
  82   // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY
  83
  84   const Int_t nratio = 10;
  85   Int_t num  [nratio]    = {kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
  86   //  Int_t denum[nratio]    = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGPi       , kPDGDeuteron , kPDGPi          , kPDGK   , -kPDGK};
  87   Int_t denum[nratio]    = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGProton   , kPDGDeuteron , kPDGPi          , kPDGK   , -kPDGK};
  88   Int_t isSum[nratio]    = {1      ,1           ,0           ,1        ,1          ,0             ,0             ,1                ,0        ,1      };
  89   Double_t scale[nratio] = {1      ,3           ,1           ,30       ,250         ,50            ,100           ,4e5              ,2       ,2      };
  90   TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, nratio, 1, nratio+1);
  91
  92   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
  93   for(Int_t iratio = 1; iratio <= nratio; iratio++){
  94     AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
  95     std::cout << num[iratio-1] << " " <<  denum[iratio-1]<< " " ;
  96     if(ratio)ratio->Print("short");
  97
  98
  99     if(!ratio) {
 100       // If the ratio is not found, try to build it!
 101       AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,  isSum[iratio-1]);
 102       // Try with the !sum, if part 1 is not found
 103       if(!part1) {
 104         part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 105         // If the sum was requested, try to recover it!
 106         if(isSum[iratio-1]) {
 107           AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0);
 108           if(part1 && part1_bar) {
 109             std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl;
 110
 111             part1 = AliParticleYield::Add(part1, part1_bar);
 112
 113           }
 114         } else if(part1) {
 115           // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 116           part1->Scale(0.5);
 117         }
 118
 119       }
 120       AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
 121       if(!part2) {// Try with the !sum, if part 2 is not found
 122         part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 123         if(isSum[iratio-1]) {
 124           AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0);
 125           if(part2 && part2_bar) part2 = AliParticleYield::Add(part2, part2_bar);
 126         } else if(part2){
 127           // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 128           part2->Scale(0.5);
 129         }
 130
 131       }
 132       ratio = AliParticleYield::Divide(part1, part2, 0, "YQ");
 133       if(ratio) {
 134         std::cout << "" << std::endl;
 135         std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
 136         ratio->Print("");
 137         std::cout << "" << std::endl;
 138       }
 139     }
 140     if(ratio){
 141       ratio->Scale(scale[iratio-1]);
 142       h->SetBinContent(iratio, ratio->GetYield());
 143       h->SetBinError  (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
 144       h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f",  scale[iratio-1]), ratio->GetLatexName()));
 145     }
 146     else {
 147       h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
 148
 149     }
 150   }
 151
 152   return h;
 153
 154
 155
 156 }
 157
 158 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges)  {
 159   // If "Separate charges" is true, tries to dump both charges are dumped
 160   TClonesArray * arrOut = new TClonesArray("AliParticleYield");
 161   const Int_t npart = 12;
 162   Int_t particles  [npart] = {kPDGPi ,kPDGK   ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3    };
 163   Int_t isSum[npart]       = {1      ,1       ,0      , 1        , 0      , 1           ,0           ,1        ,1          ,0           , 1              , 0          };
 164
 165   Int_t ipartOut = 0; // Index for the array
 166   for(Int_t ipart = 0; ipart < npart; ipart++){
 167     if(!separateCharges) {
 168       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  isSum[ipart]);
 169       if(!part && isSum[ipart]) {
 170         //Could not find the particle, but the sum was requested: build the sum!
 171         part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 172         AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 173         if(part2 && part) part = AliParticleYield::Add(part, part2);
 174         else part = 0;
 175       }
 176       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 177     }
 178     else {
 179       // ignore isSum and try to find both particles
 180       Bool_t notFound = 0;
 181       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 182       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 183       else notFound=1;
 184       // Try to find antiparticle (-pdg code)
 185       part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 186       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 187       else if (notFound) {
 188         // If neither charge was found, check if we at least have the sum
 189         part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality,  1);
 190         if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 191       }
 192
 193     }
 194   }
 195   std::cout << "Particles for thermal model fits:" << std::endl;
 196   arrOut->Print("short");
 197   std::cout << "" << std::endl;
 198   // Write GSI input file
 199   TIter it(arrOut);
 200   AliParticleYield * part = 0;
 201   ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 202   while ((part = (AliParticleYield*) it.Next())){
 203     fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
 204   }
 205   fout.close();
 206   // Write thermus file
 207   AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 208 }
 209
 210
 211 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) {
 212
 213   const Int_t npart = 11;
 214   Int_t pdg  [npart]    = {kPDGPi, kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
 215   //  Int_t isSum[npart]    = {1      ,1      ,1           ,0           ,1        ,1          ,0             ,0             ,1                ,0        ,1      };
 216   Int_t isSum[npart]    = {0,0,0,0,0,0,0,0,0,0,1};
 217   //  Double_t scale[npart] = {1      ,1      ,3           ,1           ,30       ,250         ,50            ,10            ,4e5              ,2       ,2      };
 218   Double_t scale[npart] = {1,5,30,30,200,1000,4000,2e6,2e6,20,20,};
 219   TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1, npart+1);
 220
 221   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
 222   for(Int_t ipart = 1; ipart <= npart; ipart++){
 223     AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg[ipart-1], system, energy, centrality,isSum[ipart-1]);
 224     if(!part) continue;
 225     part->Scale(scale[ipart-1]);
 226     h->SetBinContent(ipart, part->GetYield());
 227     h->SetBinError  (ipart, part->GetTotalError(0/* 0 = no normalization error */));
 228     h->GetXaxis()->SetBinLabel(ipart, Form("#splitline{%s}{%s}",Form("#times%2.0g",  scale[ipart-1]), part->GetLatexName()));
 229
 230   }
 231   return h;
 232 }