TGraphErrors* grSysErrDYieldPt = this->AddSystError(grDYieldPt, fSysErr, fParticle + "_SysErr_DiffYield_Pt");
grSysErrDYieldPt->Write();
- TF1* fncTsallis0 = this->TsallisDiffYield(GetMass(fParticle), fParticle + "_Fit_DiffYield_Pt");
- fncTsallis0->SetParLimits(0, 0, 1);
- fncTsallis0->SetParLimits(1, 4, 50);
- fncTsallis0->SetParLimits(2, 0.01, 10);
-
- grSysErrDYieldPt->Fit(fncTsallis0,"RNQ");
- fncTsallis0->Write();
-
// invariant differential yield
TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, fParticle + "_InvDiffYield_Pt");
TGraphErrors* grSysErrInvDYieldPt = this->AddSystError(grInvDYieldPt, fSysErr, fParticle + "_SysErr_InvDiffYield_Pt");
grSysErrInvDYieldPt->Write();
- TF1* fncTsallis1 = this->Tsallis(GetMass(fParticle), fParticle + "_Fit_InvDiffYield_Pt");
- fncTsallis1->SetParLimits(0, 0, 1);
- fncTsallis1->SetParLimits(1, 4, 50);
- fncTsallis1->SetParLimits(2, 0.01, 10);
-
- grSysErrInvDYieldPt->Fit(fncTsallis1,"RNQ");
- fncTsallis1->Write();
-
// invariant differential cross section
TGraphErrors* grInvXsectPt = GetInvDiffXsectionPt(grInvDYieldPt, fInelXsec, fParticle + "_InvDiffXSection_Pt");
TGraphErrors* grSysErrInvXsectPt = GetInvDiffXsectionPt(grSysErrInvDYieldPt, fInelXsec, fParticle + "_SysErr_InvDiffXSection_Pt");
grSysErrInvXsectPt->Write();
- TF1* fncTsallis2 = this->Tsallis(GetMass(fParticle), fInelXsec[0], fParticle + "_Fit_InvDiffXSection_Pt");
- fncTsallis2->SetParLimits(0, 0, 1);
- fncTsallis2->SetParLimits(1, 4, 50);
- fncTsallis2->SetParLimits(2, 0.01, 10);
-
- grSysErrInvXsectPt->Fit(fncTsallis2,"RNQ");
- fncTsallis2->Write();
-
// clean
- delete fncTsallis0;
- delete fncTsallis1;
- delete fncTsallis2;
-
delete grDYieldPt;
delete grInvDYieldPt;
delete grSysErrDYieldPt;
return grSyst;
}
-
-TF1* AliLnSpectra::Tsallis(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
-{
-//
-// Tsallis distribution
-// Phys. Rev. C 83, 064903 (2011)
-// Phys. Rev. C 75, 064901 (2007)
-//
- TF1* fnc = new TF1(name.Data(), Form("[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
- fnc->SetParNames("dN/dy","n","C");
- fnc->SetParameters(0.1, 7, 0.2);
-
- return fnc;
-}
-
-TF1* AliLnSpectra::Tsallis(Double_t m0, Double_t xsect, const TString& name, Double_t xmin, Double_t xmax) const
-{
-//
-// Tsallis distribution to fit to invariant cross section
-// Phys. Rev. C 83, 064903 (2011)
-// Phys. Rev. C 75, 064901 (2007)
-//
- TF1* fnc = new TF1(name.Data(), Form("%f*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))",xsect, m0, m0, m0, m0), xmin, xmax);
- fnc->SetParNames("dN/dy","n","C");
- fnc->SetParameters(0.1, 7, 0.2);
-
- return fnc;
-}
-
-TF1* AliLnSpectra::TsallisDiffYield(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
-{
-//
-// Tsallis distribution to fit differential yield
-// Phys. Rev. C 83, 064903 (2011)
-// Phys. Rev. C 75, 064901 (2007)
-//
- TF1* fnc = new TF1(name.Data(), Form("x*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/([1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
- fnc->SetParNames("dN/dy","n","C");
- fnc->SetParameters(0.1, 7, 0.2);
-
- return fnc;
-}
void SetInelXSection(Double_t xsec, Double_t statErr, Double_t systErr) { fInelXsec[0] = xsec; fInelXsec[1] = statErr; fInelXsec[2] = systErr; }
- TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
- TF1* Tsallis(Double_t m0, Double_t xsect, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
- TF1* TsallisDiffYield(Double_t m0, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
-
private:
AliLnSpectra(const AliLnSpectra& other);
#include <TList.h>
#include <TString.h>
#include <TH1D.h>
+#include <TF1.h>
#include <cstdlib>
TObject* FindObj(TFile* f, const TString& name)
return 0;
}
+TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// Tsallis distribution
+// Journal of Statistical Physics, Vol. 52, Nos. 1/2, 1988
+// J. Phys. G: Nucl. Part. Phys. 39 (2012)
+//
+ TF1* fnc = new TF1(name.Data(), Form("[0]*sqrt(x*x+%f)*pow(1+([1]-1)*sqrt(x*x+%f)/[2],[1]/(1-[1]))/pow(2*TMath::Pi(),3)",m0*m0,m0*m0),xmin,xmax);
+
+ fnc->SetParNames("gV","q","T");
+ fnc->SetParameters(100., 1.1, 0.07);
+
+ fnc->SetParLimits(0, 1., 1.e+7);
+ fnc->SetParLimits(1, 1.0001, 3.);
+ fnc->SetParLimits(2, 0.001, 0.3);
+
+ return fnc;
+}
+
+TF1* TsallisDYield(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// Tsallis distribution for differential yield
+// Journal of Statistical Physics, Vol. 52, Nos. 1/2, 1988
+// J. Phys. G: Nucl. Part. Phys. 39 (2012)
+//
+ TF1* fnc = new TF1(name.Data(), Form("[0]*x*sqrt(x*x+%f)*pow(1+([1]-1)*sqrt(x*x+%f)/[2],[1]/(1-[1]))/pow(2.*TMath::Pi(),2)",m0*m0,m0*m0),xmin,xmax);
+
+ fnc->SetParNames("gV","q","T");
+ fnc->SetParameters(100., 1.1, 0.07);
+
+ fnc->SetParLimits(0, 1., 1.e+7);
+ fnc->SetParLimits(1, 1.0001, 3.);
+ fnc->SetParLimits(2, 0.001, 0.3);
+
+ return fnc;
+}
+
+TF1* TsallisPareto(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// Tsallis-Pareto distribution
+// Phys. Rev. C 83, 064903 (2011)
+// Phys. Rev. C 75, 064901 (2007)
+//
+ TF1* fnc = new TF1(name.Data(), Form("[0]*([1]-1)*([1]-2)*pow(1+(sqrt(x*x+%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))", m0*m0, m0, m0), xmin, xmax);
+
+ fnc->SetParNames("dN/dy","n","C");
+ fnc->SetParameters(0.05, 7, 0.3);
+
+ fnc->SetParLimits(0, 0, 1);
+ fnc->SetParLimits(1, 4, 50);
+ fnc->SetParLimits(2, 0.01, 10);
+
+ return fnc;
+}
+
+TF1* TsallisParetoDYield(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// Tsallis-Pareto distribution for differential yield
+// Phys. Rev. C 83, 064903 (2011)
+// Phys. Rev. C 75, 064901 (2007)
+//
+ TF1* fnc = new TF1(name.Data(), Form("x*[0]*([1]-1)*([1]-2)*pow(1+(sqrt(x*x+%f)-%f)/([1]*[2]),-[1])/([1]*[2]*([1]*[2]+%f*([1]-2)))", m0*m0, m0, m0), xmin, xmax);
+
+ fnc->SetParNames("dN/dy","n","C");
+ fnc->SetParameters(0.05, 7, 0.3);
+
+ fnc->SetParLimits(0, 0, 1);
+ fnc->SetParLimits(1, 4, 50);
+ fnc->SetParLimits(2, 0.01, 10);
+
+ return fnc;
+}
+
#endif // B2_H
#include <TStyle.h>
#include <cstdlib>
+namespace CollSystem
+{
+//
+// pp collision systems
+//
+ const Int_t kNener = 3;
+ const TString kEnergyTag[kNener] = { "900GeV", "2.76TeV", "7TeV" };
+ const Double_t kEnergy[kNener] = { 0.9, 2.76, 7. };
+ const Double_t kEnergyError[kNener] = { 0 };
+ const TString kEnergyName = "#sqrt{s} (TeV)";
+};
+
namespace B2mult
{
//
// multiplicity classes
//
- const Int_t kNmult = 6;
- const TString kMultClass[kNmult] = { "ntrk0002", "ntrk0204", "ntrk0408", "ntrk0811", "ntrk1120", "ntrk20xx" };
+ const Int_t kNmult = 6;
+ const TString kMultTag[kNmult] = { "ntrk0002", "ntrk0204", "ntrk0408", "ntrk0811", "ntrk1120", "ntrk20xx" };
const Double_t kKNOmult[kNmult] = { 0.20, 0.60, 1.01, 1.60, 2.60, 4.35 };
- const Double_t kKNOmultErr[kNmult] = { 0, 0, 0, 0, 0, 0 };
+ const Double_t kKNOmultErr[kNmult] = { 0 };
+ const TString kKNOmultName = "z";
};
TString GetCollSystem(const TString& period)
for(Int_t i=0; i<2; ++i)
{
- gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_InvDiffYield_Pt\",\"%s\",0,4.5, %g, %g, \"p_{T} (GeV/c)\", \"#frac{1}{2#piN_{inel}} #frac{d^{2}N}{p_{T}dp_{T}dy} (GeV^{-2}c^{3})\", %d, \"c%d\",\"%s\")", spectra.Data(), kParticle[i].Data(), refdir.Data(),ymin, ymax, option, i, kParticle[i].Data()));
+ gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_InvDiffYield_Pt\",\"%s\",0,4.5, %g, %g, \"p_{T} (GeV/c)\", \"#frac{1}{2#piN_{ev}} #frac{d^{2}N}{p_{T}dp_{T}dy} (GeV^{-2}c^{3})\", %d, \"%s\",\"%s\")", spectra.Data(), kParticle[i].Data(), refdir.Data(),ymin, ymax, option, kParticle[i].Data(), kParticle[i].Data()));
}
}
#include "AliLnB2.h"
-Int_t B2(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Spectra.root",
- const TString& ptag = "lhc10d",
- const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Spectra.root",
- const TString& dtag = "lhc10d",
- const TString& outputfile = "~/alice/output/lhc10d-B2.root",
- const TString& otag = "lhc10d")
+Int_t B2( const TString& pSpectra = "~/alice/output/Proton-lhc10d-Spectra.root"
+ , const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Spectra.root"
+ , const TString& ptag = "lhc10d"
+ , const TString& dtag = "lhc10d"
+ , const TString& outputfile = "~/alice/output/lhc10d-B2.root"
+ , const TString& otag = "lhc10d"
+ , const Bool_t draw = 1)
{
//
// coalescence parameter
// draw
+ if(!draw) return 0;
+
for(Int_t i=0; i<kNpart; ++i)
{
gROOT->ProcessLine(Form(".x DrawB2.C+g(\"%s\",\"%s\",\"%s\")", outputfile.Data(), otag.Data(), kPrefix[i].Data()));
const TString& multCorTag = "",
Bool_t inel = 1, // for mult
Bool_t drawOutput = 1, // for batch
- const TString& species = "Proton")
+ const TString& species = "Proton",
+ Int_t lowPtBin = 5,
+ Int_t jointPtBin = 11,
+ Int_t hiPtBin = 36)
{
//
// combine TPC and TOF for protons and deuterons
//
- const Int_t kProtonLowPtBin = 5 ;
- const Int_t kProtonJointBin = 11;
- const Int_t kProtonHiPtBin = 36;
-
- const Int_t kDeuteronLowPtBin = 4;
- const Int_t kDeuteronJointBin = 6;
- const Int_t kDeuteronHiPtBin = 13;
-
const Double_t kProtonSysErr[2] = {0.08, 0.08} ;
const Double_t kDeuteronSysErr[2] = {0.10, 0.11} ;
- // -------------
+ using namespace std;
+
+ if((species != "Proton") && (species != "Deuteron"))
+ {
+ cerr << "Particle species " << species << " not implemented, only 'Proton' and 'Deuteron'." << endl;
+ exit(1);
+ }
const TString kOutputTagTPC = outputTag + "-tpc";
const TString kOutputTagTOF = outputTag + "-tof";
+ "\"" + multTag + "\","
+ "\"" + multCorTag;
- using namespace std;
-
- if(species == "Proton")
- {
- cout << "Config_Proton_TPC_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Proton_TPC_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 0,1,1,0,0)"
- , kArgTPC.Data()
- , inel
- , kProtonLowPtBin
- , kProtonJointBin));
+ cout << "Config_" << species << "_TPC_LHC10x.C" << endl << endl;
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 0,1,1,0,0)"
+ , species.Data()
+ , kArgTPC.Data()
+ , inel
+ , lowPtBin
+ , jointPtBin));
- cout << "Config_Proton_TOF_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Proton_TOF_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 1,1,1,0,0)"
- , kArgTOF.Data()
- , inel
- , kProtonJointBin
- , kProtonHiPtBin));
- }
- else if (species == "Deuteron")
- {
- cout << "Config_Deuteron_TPC_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Deuteron_TPC_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 0,1,1,0,0)"
- , kArgTPC.Data()
- , inel
- , kDeuteronLowPtBin
- , kDeuteronJointBin));
+ cout << "Config_" << species << "_TOF_LHC10x.C" << endl << endl;
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 1,1,1,0,0)"
+ , species.Data()
+ , kArgTOF.Data()
+ , inel
+ , jointPtBin
+ , hiPtBin));
- cout << "Config_Deuteron_TOF_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Deuteron_TOF_LHC10x.C+g(\"%s\", %d, 0, %d, %d, 1,1,1,0,0)"
- , kArgTOF.Data()
- , inel
- , kDeuteronJointBin
- , kDeuteronHiPtBin));
- }
- else
- {
- cerr << "Particle species " << species << " not implemented." << endl;
- exit(1);
- }
-
TString outputPtTPC = outputDir + "/" + MakeOutputName(species, kOutputTagTPC) + "-Pt.root";
TString outputPtTOF = outputDir + "/" + MakeOutputName(species, kOutputTagTOF) + "-Pt.root";
grDYieldPt->SetLineColor(kRed);
grDYieldPt->Draw("P");
- // fitting function
- TF1* fit0 = (TF1*)FindObj(finput, tag, particle + "_Fit_DiffYield_Pt");
- fit0->SetLineWidth(1);
- fit0->SetLineColor(kRed);
- fit0->Draw("same");
-
// invariant differential yields
c0->cd(2);
grInvDYieldPt->SetLineColor(kRed);
grInvDYieldPt->Draw("P");
- // fitting function
- TF1* fit1 = (TF1*)FindObj(finput, tag, particle + "_Fit_InvDiffYield_Pt");
- fit1->SetLineWidth(1);
- fit1->SetLineColor(kRed);
- fit1->Draw("same");
-
// invariant differential cross section
c0->cd(3);
grInvDXsectPt->SetMarkerColor(kRed);
grInvDXsectPt->SetLineColor(kRed);
grInvDXsectPt->Draw("P");
-
- // fitting function
- TF1* fit2 = (TF1*)FindObj(finput, tag, particle + "_Fit_InvDiffXSection_Pt");
- fit2->SetLineWidth(1);
- fit2->SetLineColor(kRed);
- fit2->Draw("same");
}
//\r
// call Config_XXX for each period, merge the corrected pt and then get the results\r
//\r
- using namespace std;\r
- \r
const Int_t kNper = 4;\r
const TString kPeriod[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
const TString kOutputTag[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
\r
- Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;\r
- Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;\r
+ Int_t lowbin = (species=="Proton") ? 5 : 4;\r
+ Int_t jointbin = (species=="Proton") ? 11 : 6;\r
+ Int_t hibin = (species=="Proton") ? 36 : 15;\r
+ \r
+ Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;\r
+ Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;\r
+ \r
+ using namespace std;\r
+ \r
+ if( (option<0) || (option>2) || ((species != "Proton") && (species != "Deuteron")))\r
+ {\r
+ cerr << "unknown species/option: " << species << "/" << option << endl;\r
+ cerr << "species: Proton or Deuteron, options: 0 (TPC), 1 (TOF), 2 (TPCTOF)" << endl;\r
+ exit(1);\r
+ }\r
\r
TFileMerger m;\r
\r
+ "\"" + multTag + "\","\r
+ "\"" + multCorTag;\r
\r
- if(species=="Proton" && option==0)\r
- {\r
- cout << "Config_Proton_TPC_LHC10x.C" << endl << endl;\r
- gROOT->ProcessLine(Form(".x Config_Proton_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), inel));\r
- }\r
- else if(species=="Proton" && option==1)\r
- {\r
- cout << "Config_Proton_TOF_LHC10x.C" << endl << endl;\r
- gROOT->ProcessLine(Form(".x Config_Proton_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), inel));\r
- }\r
- else if(species=="Deuteron" && option==0)\r
- {\r
- cout << "Config_Deuteron_TPC_LHC10x.C" << endl << endl;\r
- gROOT->ProcessLine(Form(".x Config_Deuteron_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), inel));\r
- }\r
- else if(species=="Deuteron" && option==1)\r
- {\r
- cout << "Config_Deuteron_TOF_LHC10x.C" << endl << endl;\r
- gROOT->ProcessLine(Form(".x Config_Deuteron_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), inel));\r
- }\r
- else if((species=="Proton" || species=="Deuteron") && option==2)\r
- {\r
- cout << "Config_TPCTOF_LHC10x.C" << endl << endl;\r
- gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\",%d,0,\"%s\")", arg.Data(), kNormToInel[i], species.Data()));\r
- }\r
- else\r
+ switch(option)\r
{\r
- cerr << "unknown species/option: " << species << "/" << option << endl;\r
- cerr << "usage: Proton/0, Proton/1, Proton/2 or Deuteron" << endl;\r
- exit(1);\r
+ case 0:\r
+ cout << "Config_" << species << "_TPC_LHC10x.C" << endl << endl;\r
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %d, 0)", species.Data(), arg.Data(), inel));\r
+ break;\r
+ case 1:\r
+ cout << "Config_" << species << "_LHC10x.C" << endl << endl;\r
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %d, 0)", species.Data(), arg.Data(), inel));\r
+ break;\r
+ case 2:\r
+ cout << "Config_TPCTOF_LHC10x.C" << endl << endl;\r
+ gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\", %d, 0, \"%s\", %d, %d, %d)", arg.Data(), inel, species.Data(), lowbin, jointbin, hibin));\r
+ break;\r
}\r
\r
TString ptfile = outputDir + "/" + species + "-" + kOutputTag[i] + "-Pt.root";\r
// if option = 1 then use Config_XXX_TOF
// if option = 2 then use Config_TPCTOF
//
+ const Bool_t kINEL[kNmult] = { 1, 0 }; // only normalize first bin
+
+ Int_t lowbin = (species=="Proton") ? 5 : 4;
+ Int_t jointbin = (species=="Proton") ? 11 : 6;
+ Int_t hibin = (species=="Proton") ? 31 : 12;
+
+ Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;
+ Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;
+
using namespace B2mult;
using namespace std;
- const Bool_t kNormToInel[kNmult] = { 1, 0 }; // only normalize first bin
-
- Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;
- Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;
+ if( (option<0) || (option>2) || ((species != "Proton") && (species != "Deuteron")))
+ {
+ cerr << "unknown species/option: " << species << "/" << option << endl;
+ cerr << "species: Proton or Deuteron, options: 0 (TPC), 1 (TOF), 2 (TPCTOF)" << endl;
+ exit(1);
+ }
TFileMerger m1,m2;
+ "\"" + "-" + kMultClass[i] + "\"," // data
+ "\"" + ""; // same simulations for all mult
- if(species=="Proton" && option==0)
- {
- cout << "Config_Proton_TPC_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Proton_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
- }
- else if(species=="Proton" && option==1)
- {
- cout << "Config_Proton_TOF_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Proton_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
- }
- else if(species=="Deuteron" && option==0)
- {
- cout << "Config_Deuteron_TPC_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Deuteron_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
- }
- else if(species=="Deuteron" && option==1)
- {
- cout << "Config_Deuteron_TOF_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_Deuteron_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
- }
- else if((species=="Proton" || species=="Deuteron") && option==2)
- {
- cout << "Config_TPCTOF_LHC10x.C" << endl << endl;
- gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\",%d,0,\"%s\")", arg.Data(), kNormToInel[i], species.Data()));
- }
- else
+ switch(option)
{
- cerr << "unknown species/option: " << species << "/" << option << endl;
- cerr << "usage: Proton/0, Proton/1, Proton/2 or Deuteron" << endl;
- exit(1);
+ case 0:
+ cout << "Config_" << species << "_TPC_LHC10x.C" << endl << endl;
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %d, 0)", species.Data(), arg.Data(), kINEL[i]));
+ break;
+ case 1:
+ cout << "Config_" << species << "_LHC10x.C" << endl << endl;
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %d, 0)", species.Data(), arg.Data(), kINEL[i]));
+ break;
+ case 2:
+ cout << "Config_TPCTOF_LHC10x.C" << endl << endl;
+ gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\", %d, 0, \"%s\", %d, %d, %d)", arg.Data(), kINEL[i], species.Data(), lowbin, jointbin, hibin));
+ break;
}
ratio[i] = outputDir + "/" + species + "-" + outputTag + "-Ratio.root";
* provided "as is" without express or implied warranty. *
**************************************************************************/
-// particle ratios as a function of multiplicity
+// d/p ratio
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TStyle.h>
+#include <TFitResult.h>
+#include <TFitResultPtr.h>
+#include <TVirtualFitter.h>
#include "B2.h"
#include "Config.h"
void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr);
void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle);
-void RatioMult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
- const TString& ptag = "lhc10d",
- const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
- const TString& dtag = "lhc10d")
+void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root"
+ , const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root"
+ , const TString& ptag = "lhc10d"
+ , const TString& dtag = "lhc10d"
+ , const Int_t nx = B2mult::kNmult
+ , const TString* xtag = B2mult::kMultTag
+ , const Double_t* x = B2mult::kKNOmult
+ , const Double_t* xerr = B2mult::kKNOmultErr
+ , const TString& xname = B2mult::kKNOmultName
+ , const Bool_t tsallis = 1
+ , const TString& outputfile = "~/alice/output/Ratio-Mult.root"
+ , const TString& otag = "Tsallis")
{
//
-// particle ratios as a function of multiplicity (from fitting distributions)
+// d/p ratio as a function of X
//
- using namespace B2mult;
-
+ const Int_t kNspec = 2;
const Int_t kNpart = 2;
- const TString kProton[kNpart] = { "Proton", "AntiProton" };
- const TString kDeuteron[kNpart] = { "Deuteron", "AntiDeuteron" };
+ const TString kParticle[kNspec][kNpart] = { {"Proton", "AntiProton"}, { "Deuteron", "AntiDeuteron"} };
- // open files
+ //TVirtualFitter::SetDefaultFitter("Minuit2");
- TFile* fproton = new TFile(pSpectra.Data());
+ TFile* fproton = new TFile(pSpectra.Data(), "read");
if(fproton->IsZombie()) exit(1);
- TFile* fdeuteron = new TFile(dSpectra.Data());
+ TFile* fdeuteron = new TFile(dSpectra.Data(), "read");
if(fdeuteron->IsZombie()) exit(1);
+ TFile* finput[kNspec] = { fproton, fdeuteron };
+ TString tag[kNspec] = { ptag, dtag };
+
// particle ratios
- TGraphErrors* grRatio[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grMixRatio[kNpart] = { new TGraphErrors(), new TGraphErrors() };
+ TGraphErrors* grRatio[kNspec] = { new TGraphErrors(), new TGraphErrors() };
+ TGraphErrors* grMixRatio[kNspec] = { new TGraphErrors(), new TGraphErrors() };
// model parameters
- TGraphErrors* grProtondNdy[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grProtonN[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grProtonC[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grProtonChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- //TGraphErrors* grProtonProb[kNpart] = { new TGraphErrors(), new TGraphErrors() };
-
- TGraphErrors* grDeuterondNdy[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grDeuteronN[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grDeuteronC[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- TGraphErrors* grDeuteronChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
- //TGraphErrors* grDeuteronProb[kNpart] = { new TGraphErrors(), new TGraphErrors() };
-
- for(Int_t i=0; i<kNmult; ++i)
+ TGraphErrors* grP0[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+ TGraphErrors* grP1[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+ TGraphErrors* grP2[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+ TGraphErrors* grChi2Ndf[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+
+ TGraphErrors* grdNdy[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+
+ TGraphErrors* grDYieldPt[nx][kNspec][kNpart];
+ TF1* fncTsallis[nx][kNspec][kNpart];
+
+ for(Int_t i=0; i<nx; ++i)
{
- TF1* fncProton[kNpart];
- TF1* fncDeuteron[kNpart];
+ Double_t dNdy[kNspec][kNpart];
+ Double_t dNdyErr[kNspec][kNpart];
- for(Int_t j=0; j<kNpart; ++j)
+ for(Int_t j=0; j<kNspec; ++j)
{
- fncProton[j] = (TF1*)FindObj(fproton, ptag + "-" + kMultClass[i], kProton[j] + "_Fit_InvDiffYield_Pt");
- fncDeuteron[j] = (TF1*)FindObj(fdeuteron, dtag + "-" + kMultClass[i], kDeuteron[j] + "_Fit_InvDiffYield_Pt");
+ for(Int_t k=0; k<kNpart; ++k)
+ {
+ grDYieldPt[i][j][k] = (TGraphErrors*)FindObj(finput[j], tag[j] + "-" + xtag[i], kParticle[j][k] + "_SysErr_DiffYield_Pt");
+
+ if(tsallis)
+ {
+ fncTsallis[i][j][k] = TsallisDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_DiffYield_Pt",0,10);
+ }
+ else
+ {
+ fncTsallis[i][j][k] = TsallisParetoDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_DiffYield_Pt",0,10);
+ }
+
+ TFitResultPtr r = grDYieldPt[i][j][k]->Fit(fncTsallis[i][j][k], "RENSQ");
+
+ if(tsallis)
+ {
+ dNdy[j][k] = fncTsallis[i][j][k]->Integral(0,100);
+ dNdyErr[j][k] = fncTsallis[i][j][k]->IntegralError(0,100,r->GetParams(), r->GetCovarianceMatrix().GetMatrixArray());
+ }
+ else
+ {
+ dNdy[j][k] = fncTsallis[i][j][k]->GetParameter(0);
+ dNdyErr[j][k] = fncTsallis[i][j][k]->GetParError(0);
+ }
+
+ Int_t status = r;
+ printf("status: %d\tdN/dy: %g +/- %g\n",status,dNdy[j][k],dNdyErr[j][k]);
+
+ grdNdy[j][k]->SetPoint(i, x[i], dNdy[j][k]);
+ grdNdy[j][k]->SetPointError(i, xerr[i], dNdyErr[j][k]);
+ }
}
- // integrated ratios
-
- Double_t dNdy[kNpart] = {fncProton[0]->GetParameter(0), fncDeuteron[0]->GetParameter(0)};
- Double_t dNdyErr[kNpart] = {fncProton[0]->GetParError(0), fncDeuteron[0]->GetParError(0)};
-
- Double_t dNdyBar[kNpart] = {fncProton[1]->GetParameter(0), fncDeuteron[1]->GetParameter(0)};
- Double_t dNdyBarErr[kNpart] = {fncProton[1]->GetParError(0), fncDeuteron[1]->GetParError(0)};
-
// ratios
- Double_t ratio[kNpart], ratioErr[kNpart];
- for(Int_t j=0; j<kNpart; ++j)
+ Double_t ratio[kNspec], ratioErr[kNspec];
+ for(Int_t j=0; j<kNspec; ++j)
{
- GetRatio(dNdyBar[j], dNdy[j], dNdyBarErr[j], dNdyErr[j], ratio[j], ratioErr[j]);
- grRatio[j]->SetPoint(i, kKNOmult[i], ratio[j]);
+ GetRatio(dNdy[j][1], dNdy[j][0], dNdyErr[j][1], dNdyErr[j][0], ratio[j], ratioErr[j]);
+
+ grRatio[j]->SetPoint(i, x[i], ratio[j]);
grRatio[j]->SetPointError(i, 0, ratioErr[j]);
}
// mixed ratios
Double_t mixRatio[kNpart], mixRatioErr[kNpart];
-
- GetRatio(dNdy[1], dNdy[0], dNdyErr[1], dNdyErr[0], mixRatio[0], mixRatioErr[0]);
- GetRatio(dNdyBar[1], dNdyBar[0], dNdyBarErr[1], dNdyBarErr[0], mixRatio[1], mixRatioErr[1]);
-
for(Int_t j=0; j<kNpart; ++j)
{
- grMixRatio[j]->SetPoint(i, kKNOmult[i], mixRatio[j]);
- grMixRatio[j]->SetPointError(i, 0, mixRatioErr[j]);
+ GetRatio(dNdy[1][j], dNdy[0][j], dNdyErr[1][j], dNdyErr[0][j], mixRatio[j], mixRatioErr[j]);
+
+ grMixRatio[j]->SetPoint(i, x[i], mixRatio[j]);
+ grMixRatio[j]->SetPointError(i, xerr[i], mixRatioErr[j]);
}
// model parameters
+ for(Int_t j=0; j<kNspec; ++j)
+ {
+ for(Int_t k=0; k<kNpart; ++k)
+ {
+ grP0[j][k]->SetPoint(i, x[i], fncTsallis[i][j][k]->GetParameter(0));
+ grP0[j][k]->SetPointError(i, xerr[i], fncTsallis[i][j][k]->GetParError(0));
+
+ grP1[j][k]->SetPoint(i, x[i], fncTsallis[i][j][k]->GetParameter(1));
+ grP1[j][k]->SetPointError(i, xerr[i], fncTsallis[i][j][k]->GetParError(1));
+
+ grP2[j][k]->SetPoint(i, x[i], fncTsallis[i][j][k]->GetParameter(2));
+ grP2[j][k]->SetPointError(i, xerr[i], fncTsallis[i][j][k]->GetParError(2));
+
+ grChi2Ndf[j][k]->SetPoint(i, x[i], fncTsallis[i][j][k]->GetChisquare()/fncTsallis[i][j][k]->GetNDF());
+ }
+ }
+ }
+
+ // save
+
+ TFile* foutput = new TFile(outputfile.Data(), "recreate");
+ if (foutput->IsZombie()) exit(1);
+
+ foutput->mkdir(otag.Data());
+ foutput->cd(otag.Data());
+
+ for(Int_t i=0; i<kNspec; ++i)
+ {
+ grRatio[i]->SetName(Form("%s%s_Ratio", kParticle[i][1].Data(), kParticle[i][0].Data()));
+ grRatio[i]->Write();
+
+ grMixRatio[i]->SetName(Form("%s%s_Ratio", kParticle[1][i].Data(), kParticle[0][i].Data()));
+ grMixRatio[i]->Write();
+ }
+
+ for(Int_t i=0; i<kNspec; ++i)
+ {
for(Int_t j=0; j<kNpart; ++j)
{
- grProtondNdy[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(0));
- grProtondNdy[j]->SetPointError(i, 0, fncProton[j]->GetParError(0));
-
- grProtonN[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(1));
- grProtonN[j]->SetPointError(i, 0, fncProton[j]->GetParError(1));
-
- grProtonC[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(2));
- grProtonC[j]->SetPointError(i, 0, fncProton[j]->GetParError(2));
-
- grProtonChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetChisquare()/fncProton[j]->GetNDF());
- //grProtonProb[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetProb());
+ grP0[i][j]->SetName(Form("%s_P0", kParticle[i][j].Data()));
+ grP1[i][j]->SetName(Form("%s_P1", kParticle[i][j].Data()));
+ grP2[i][j]->SetName(Form("%s_P2", kParticle[i][j].Data()));
- // deuteron
- grDeuterondNdy[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(0));
- grDeuterondNdy[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(0));
+ grP0[i][j]->Write();
+ grP1[i][j]->Write();
+ grP2[i][j]->Write();
- grDeuteronN[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(1));
- grDeuteronN[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(1));
+ grdNdy[i][j]->SetName(Form("%s_dNdy", kParticle[i][j].Data()));
+ grChi2Ndf[i][j]->SetName(Form("%s_Chi2Ndf", kParticle[i][j].Data()));
- grDeuteronC[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(2));
- grDeuteronC[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(2));
-
- grDeuteronChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetChisquare()/fncDeuteron[j]->GetNDF());
- //grDeuteronProb[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetProb());
+ grdNdy[i][j]->Write();
+ grChi2Ndf[i][j]->Write();
}
}
- delete fproton;
- delete fdeuteron;
-
// draw
TStyle* st = new TStyle();
const Int_t kNCol = 4;
- TCanvas* c0 = new TCanvas("c0", "proton model parameters");
- c0->Divide(kNCol,2);
+ TCanvas* c0[kNspec];
- for(Int_t j=0; j<kNpart; ++j)
+ for(Int_t i=0; i<kNspec; ++i)
{
- c0->cd(kNCol*j+1);
- Draw(grProtondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
+ c0[i] = new TCanvas(Form("c0.%s",kParticle[i][0].Data()), Form("%s model parameters",kParticle[i][0].Data()));
+ c0[i]->Divide(kNCol,2);
- c0->cd(kNCol*j+2);
- Draw(grProtonN[j], kFullCircle, kBlue, "z", "n");
-
- c0->cd(kNCol*j+3);
- Draw(grProtonC[j], kFullCircle, kBlue, "z", "C");
-
- c0->cd(kNCol*j+4);
- Draw(grProtonChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
-
- /* c0->cd(kNCol*j+5);
- Draw(grProtonProb[j], kFullCircle, kBlue, "z", "Prob");
- */
+ for(Int_t j=0; j<kNpart; ++j)
+ {
+ c0[i]->cd(kNCol*j+1);
+ Draw(grP0[i][j], kFullCircle, kBlue, xname, fncTsallis[0][0][0]->GetParName(0));
+
+ c0[i]->cd(kNCol*j+2);
+ Draw(grP1[i][j], kFullCircle, kBlue, xname, fncTsallis[0][0][0]->GetParName(1));
+
+ c0[i]->cd(kNCol*j+3);
+ Draw(grP2[i][j], kFullCircle, kBlue, xname, fncTsallis[0][0][0]->GetParName(2));
+
+ c0[i]->cd(kNCol*j+4);
+ Draw(grChi2Ndf[i][j], kFullCircle, kBlue, xname, "#chi^{2}/ndf");
+ }
}
- TCanvas* c1 = new TCanvas("c1", "deuteron model parameters");
- c1->Divide(kNCol,2);
+ TCanvas* c1 = new TCanvas("c1", "Particle ratios");
+ c1->Divide(2,2);
- for(Int_t j=0; j<kNpart; ++j)
- {
- c1->cd(kNCol*j+1);
- Draw(grDeuterondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
-
- c1->cd(kNCol*j+2);
- Draw(grDeuteronN[j], kFullCircle, kBlue, "z", "n");
-
- c1->cd(kNCol*j+3);
- Draw(grDeuteronC[j], kFullCircle, kBlue, "z", "C");
-
- c1->cd(kNCol*j+4);
- Draw(grDeuteronChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
-
- /* c1->cd(kNCol*j+5);
- Draw(grDeuteronProb[j], kFullCircle, kBlue, "z", "Prob");
- */
- }
+ c1->cd(1);
+ Draw(grRatio[0], kFullCircle, kRed, xname, "#bar{p}/p");
+
+ c1->cd(2);
+ Draw(grRatio[1], kFullCircle, kRed, xname, "#bar{d}/d");
+
+ c1->cd(3);
+ Draw(grMixRatio[0], kFullCircle, kRed, xname, "d/p");
+
+ c1->cd(4);
+ Draw(grMixRatio[1], kFullCircle, kRed, xname, "#bar{d}/#bar{p}");
- TCanvas* c2 = new TCanvas("c2","particle ratios");
- c2->Divide(2,2);
+ // spectra
- c2->cd(1);
- Draw(grRatio[0], kFullCircle, kRed, "z", "#bar{p}/p");
+ TCanvas* c2[kNspec];
- c2->cd(2);
- Draw(grRatio[1], kFullCircle, kRed, "z", "#bar{d}/d");
+ for(Int_t j=0; j<kNspec; ++j)
+ {
+ c2[j] = new TCanvas(Form("c2.%s",kParticle[j][0].Data()), Form("%s data",kParticle[j][0].Data()));
+ c2[j]->Divide(nx,2);
+
+ for(Int_t k=0; k<kNpart; ++k)
+ {
+ for(Int_t i=0; i<nx; ++i)
+ {
+ gPad->SetLogy(0);
+ c2[j]->cd(nx*k+i+1);
+ Draw(grDYieldPt[i][j][k], kFullCircle, kBlue, "p_{T} (GeV/c)", "DYield");
+ fncTsallis[i][j][k]->Draw("same");
+ }
+ }
+ }
- c2->cd(3);
- Draw(grMixRatio[0], kFullCircle, kRed, "z", "d/p");
+ delete foutput;
+ delete fproton;
+ delete fdeuteron;
- c2->cd(4);
- Draw(grMixRatio[1], kFullCircle, kRed, "z", "#bar{d}/#bar{p}");
+ // free memory
}
void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff