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) ;
  13 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
  14
  15 // Plots ratios for QM and saves input files for thermal models
  16
  17 TClonesArray * PlotRatiosForQM14() {
  18   LoadLibs();
  19   TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
  20   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
  21   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
  22   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
  23   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
  24   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
  25   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
  26   arr->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
  27
  28   TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
  29
  30   // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
  31   // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
  32
  33   TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
  34   c1->SetLogy();
  35   //  GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010")->Draw();
  36   GetHistoRatios(arr, AliParticleYield::kCSpp,   7000, ".*")->Draw("");
  37   return arr;
  38 }
  39
  40 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) {
  41
  42   const Int_t nratio = 10;
  43   Int_t num  [nratio] = {kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
  44   Int_t denum[nratio] = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGProton   , kPDGDeuteron , kPDGPi          , kPDGK   , -kPDGK};
  45   Int_t isSum[nratio] = {1      ,1           ,0           ,1        ,1          ,0             ,0             ,1                ,0        ,1};
  46   TH1F * h = new TH1F("hRatio", "hRatio", nratio, 1, nratio+1);
  47
  48   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
  49   for(Int_t iratio = 1; iratio <= nratio; iratio++){
  50     AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
  51     std::cout << num[iratio-1] << " " <<  denum[iratio-1]<< " " ;
  52     if(ratio)ratio->Print("short");
  53
  54
  55     if(!ratio) {
  56       // If the ratio is not found, try to build it!
  57       AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,  isSum[iratio-1]);
  58       if(!part1) {// Try with the !sum, if part 1 is not found
  59         part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
  60       }
  61       AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
  62       if(!part2) {// Try with the !sum, if part 2 is not found
  63         part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
  64       }
  65       ratio = AliParticleYield::Divide(part1, part2);
  66       if(ratio) {
  67         std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
  68         ratio->Print("short");
  69       }
  70     }
  71     if(ratio){
  72       h->SetBinContent(iratio, ratio->GetYield());
  73       h->SetBinError  (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
  74       h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",ratio->GetCentr().Data() , ratio->GetLatexName()));
  75     }
  76     else {
  77       h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
  78
  79     }
  80   }
  81
  82   return h;
  83
  84
  85
  86 }
  87
  88 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges)  {
  89   // If "Separate charges" is true, tries to dump both charges are dumped
  90   TClonesArray * arrOut = new TClonesArray("AliParticleYield");
  91   const Int_t npart = 12;
  92   Int_t particles  [npart] = {kPDGPi ,kPDGK   ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3    };
  93   Int_t isSum[npart]       = {1      ,1       ,0      , 1        , 0      , 1           ,0           ,1        ,1          ,0           , 1              , 0          };
  94
  95   Int_t ipartOut = 0; // Index for the array
  96   for(Int_t ipart = 0; ipart < npart; ipart++){
  97     if(!separateCharges) {
  98       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  isSum[ipart]);
  99       if(!part && isSum[ipart]) {
 100         //Could not find the particle, but the sum was requested: build the sum!
 101         part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 102         AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 103         if(part2 && part) part = AliParticleYield::Add(part, part2);
 104         else part = 0;
 105       }
 106       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 107     }
 108     else {
 109       // ignore isSum and try to find both particles
 110       Bool_t notFound = 0;
 111       AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality,  0);
 112       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 113       else notFound=1;
 114       // Try to find antiparticle (-pdg code)
 115       part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality,  0);
 116       if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 117       else if (notFound) {
 118         // If neither charge was found, check if we at least have the sum
 119         part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality,  1);
 120         if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 121       }
 122
 123     }
 124   }
 125   std::cout << "Particles for thermal model fits:" << std::endl;
 126   arrOut->Print("short");
 127   std::cout << "" << std::endl;
 128   // Write GSI input file
 129   TIter it(arrOut);
 130   AliParticleYield * part = 0;
 131   ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 132   while ((part = (AliParticleYield*) it.Next())){
 133     fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
 134   }
 135   fout.close();
 136   // Write thermus file
 137   AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 138 }