rewrite d/p ratio and simplify some macros
authoreserradi <eserradi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 Feb 2013 13:09:11 +0000 (13:09 +0000)
committereserradi <eserradi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 8 Feb 2013 13:09:11 +0000 (13:09 +0000)
PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx
PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h
PWGLF/SPECTRA/Nuclei/B2/B2.h
PWGLF/SPECTRA/Nuclei/B2/Config.h
PWGLF/SPECTRA/Nuclei/B2/macros/B2.C
PWGLF/SPECTRA/Nuclei/B2/macros/Config_TPCTOF_LHC10x.C
PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C
PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C
PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C
PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C

index 20c6346..54bfaf8 100644 (file)
@@ -88,14 +88,6 @@ Int_t AliLnSpectra::Exec()
        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");
@@ -104,14 +96,6 @@ Int_t AliLnSpectra::Exec()
        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");
@@ -120,20 +104,8 @@ Int_t AliLnSpectra::Exec()
        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;
@@ -269,45 +241,3 @@ TGraphErrors* AliLnSpectra::AddSystError(const TGraphErrors* gr, Double_t percen
        
        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;
-}
index ceb5878..d320548 100644 (file)
@@ -38,10 +38,6 @@ class AliLnSpectra: public TObject
        
        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);
index 21e424a..bfc0b18 100644 (file)
@@ -12,6 +12,7 @@
 #include <TList.h>
 #include <TString.h>
 #include <TH1D.h>
+#include <TF1.h>
 #include <cstdlib>
 
 TObject* FindObj(TFile* f, const TString& name)
@@ -96,4 +97,80 @@ Double_t GetMass(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
index 3a8e3aa..aa598c1 100644 (file)
 #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)
@@ -313,7 +326,7 @@ void DrawOutputSpectraMult(const TString& spectra, const TString& species, Doubl
        
        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()));
        }
 }
 
index e248955..c80ca62 100644 (file)
 
 #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
@@ -57,6 +58,8 @@ Int_t B2(const TString& pSpectra   = "~/alice/output/Proton-lhc10d-Spectra.root"
        
        // 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()));
index c0d7f44..e6facb1 100644 (file)
@@ -32,23 +32,24 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir   = "~/alice/input",
                            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";
@@ -67,46 +68,22 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir   = "~/alice/input",
                              + "\"" + 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";
        
index 4c75bf7..1529284 100644 (file)
@@ -58,12 +58,6 @@ void DrawSpectra(const TString& inputFile="spectra.root", const TString& tag="lh
        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);
@@ -82,12 +76,6 @@ void DrawSpectra(const TString& inputFile="spectra.root", const TString& tag="lh
        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);
@@ -105,10 +93,4 @@ void DrawSpectra(const TString& inputFile="spectra.root", const TString& tag="lh
        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");
 }
index f682fc2..bac11e0 100644 (file)
@@ -38,14 +38,25 @@ Int_t LHC10bcde(const TString& species    = "Deuteron",
 //\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
@@ -60,36 +71,20 @@ Int_t LHC10bcde(const TString& species    = "Deuteron",
                              + "\""  + 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
index c60c1fa..4420c3b 100644 (file)
@@ -38,13 +38,24 @@ Int_t LHC10xMult(const TString& species   = "Proton",
 // 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;
        
@@ -66,36 +77,20 @@ Int_t LHC10xMult(const TString& species   = "Proton",
                              + "\"" + "-" + 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";
index d066cd8..cf5d8ba 100644 (file)
@@ -13,7 +13,7 @@
  * 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>
@@ -24,6 +24,9 @@
 #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();
@@ -163,64 +222,70 @@ void RatioMult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spec
        
        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)