// invariant differential yields
- TGraphErrors* grPrtInvDYieldPt = (TGraphErrors*)FindObj(finput1, fProtonTag, prefix + "Proton_InvDiffYield_Pt");
- TGraphErrors* grNucInvDYieldPt = (TGraphErrors*)FindObj(finputA, fNucleusTag, prefix + fNucleusName + "_InvDiffYield_Pt");
-
- TGraphErrors* grSysErrPrtInvDYieldPt = (TGraphErrors*)FindObj(finput1, fProtonTag, prefix + "Proton_SysErr_InvDiffYield_Pt");
- TGraphErrors* grSysErrNucInvDYieldPt = (TGraphErrors*)FindObj(finputA, fNucleusTag, prefix + fNucleusName + "_SysErr_InvDiffYield_Pt");
+ TGraphErrors* grPrtInvDYieldPt = FindObj<TGraphErrors>(finput1, fProtonTag, prefix + "Proton_InvDiffYield_Pt");
+ TGraphErrors* grNucInvDYieldPt = FindObj<TGraphErrors>(finputA, fNucleusTag, prefix + fNucleusName + "_InvDiffYield_Pt");
TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
TGraphErrors* grB2Pt = this->GetBAPt(grPrtInvDYieldPt, grNucInvDYieldPt, Form("B2%s_Pt",suffix.Data()));
- TGraphErrors* grSysErrB2Pt = this->GetBAPt(grSysErrPrtInvDYieldPt, grSysErrNucInvDYieldPt, Form("B2%s_SysErr_Pt",suffix.Data()), 0.025);
grB2Pt->Write();
- grSysErrB2Pt->Write();
// homogeneity volume
TGraphErrors* grR3Pt = this->Rside2Rlong(grB2Pt, Form("R3%s_Pt", suffix.Data()), fCd);
- TGraphErrors* grSysErrR3Pt = this->Rside2Rlong(grSysErrB2Pt, Form("R3%s_SysErr_Pt", suffix.Data()), fCd, 0.025);
grR3Pt->Write();
- grSysErrR3Pt->Write();
- // clean
+ // propagate systematic errors if possible
+
+ TString grPrtName = fProtonTag + "/" + prefix + "Proton_SystErr_InvDiffYield_Pt;1";
+ TString grNucName = fNucleusTag + "/" + prefix + fNucleusName + "_SystErr_InvDiffYield_Pt;1";
+
+ TGraphErrors* grSystErrPrtInvDYieldPt = dynamic_cast<TGraphErrors*>(finput1->Get(grPrtName.Data()));
+ TGraphErrors* grSystErrNucInvDYieldPt = dynamic_cast<TGraphErrors*>(finput1->Get(grNucName.Data()));
+
+ if( (grSystErrPrtInvDYieldPt != 0) && (grSystErrNucInvDYieldPt != 0) )
+ {
+ TGraphErrors* grSystErrB2Pt = this->GetBAPt(grSystErrPrtInvDYieldPt, grSystErrNucInvDYieldPt, Form("SystErr_B2%s_Pt",suffix.Data()));
+
+ grSystErrB2Pt->Write();
+
+ TGraphErrors* grSystErrR3Pt = this->Rside2Rlong(grSystErrB2Pt, Form("SystErr_R3%s_Pt", suffix.Data()), fCd);
+ grSystErrR3Pt->Write();
+
+ delete grSystErrB2Pt;
+ delete grSystErrR3Pt;
+ }
+ else
+ {
+ this->Warning("Run", "systematic errors are not propagated");
+ }
delete grB2Pt;
- delete grSysErrB2Pt;
delete grR3Pt;
- delete grSysErrR3Pt;
delete foutput;
delete finput1;
return 0;
}
-TGraphErrors* AliLnB2::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name, Double_t errPt) const
+TGraphErrors* AliLnB2::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name) const
{
//
// coalescence parameter
TGraphErrors* grBAPt = new TGraphErrors();
grBAPt->SetName(name.Data());
- Int_t offset = 0;
- Double_t pt, yPrt, ptNuc, yNuc;
-
- grPrtInvDYieldPt->GetPoint(0,pt,yPrt);
-
- // find offset
- for(Int_t i=0; i < grNucInvDYieldPt->GetN(); ++i)
+ for(Int_t i=0, j=0; i < grNucInvDYieldPt->GetN(); ++i)
{
+ Double_t ptNuc, yNuc;
+
grNucInvDYieldPt->GetPoint(i, ptNuc, yNuc);
- if(TMath::Abs(fA*pt-ptNuc)<0.001)
- {
- offset = i;
- break;
- }
- }
-
- for(Int_t i=0, j=0; i < grPrtInvDYieldPt->GetN() && i+offset < grNucInvDYieldPt->GetN(); ++i)
- {
- grPrtInvDYieldPt->GetPoint(i,pt,yPrt);
- grNucInvDYieldPt->GetPoint(i+offset,ptNuc, yNuc);
+
+ if(ptNuc<0.8) continue; // acceptance
+
+ Double_t yPrt = grPrtInvDYieldPt->Eval(ptNuc/fA); // interpolate
if(yPrt == 0 || yNuc == 0 ) continue;
- if(TMath::Abs(fA*pt-ptNuc)>0.001) continue; // compare at equal momentum per nucleon
Double_t bA = yNuc/TMath::Power(yPrt,fA);
// error
- Double_t ePrt = grPrtInvDYieldPt->GetErrorY(i);
+ Double_t ePrt = this->GetErrorY(grPrtInvDYieldPt, ptNuc/fA);
Double_t eNuc = grNucInvDYieldPt->GetErrorY(i);
+ Double_t errPt = grNucInvDYieldPt->GetErrorX(i)/fA;
Double_t errBA = bA*TMath::Sqrt(TMath::Power(eNuc/yNuc,2) + TMath::Power(fA*ePrt/yPrt,2));
- grBAPt->SetPoint(j, pt, bA);
+ grBAPt->SetPoint(j, ptNuc/fA, bA);
grBAPt->SetPointError(j++, errPt, errBA);
}
return grBAPt;
}
+Double_t AliLnB2::GetErrorY(const TGraphErrors* gr, Double_t x0) const
+{
+//
+// estimate error of gr(x0) with the closest point to x0
+//
+ const Double_t kEpsilon = 1.e-6;
+ const Double_t kUnknownError = 1.e+6;
+
+ for(Int_t i=0; i<gr->GetN(); ++i)
+ {
+ Double_t x = gr->GetX()[i];
+
+ if( TMath::Abs(x-x0) < kEpsilon) return gr->GetErrorY(i);
+
+ if( ((i == 0) && (x > x0)) || ((i == (gr->GetN()-1)) && (x < x0)) )
+ {
+ this->Warning("GetErrorY", "%f is out of bounds",x);
+ return kUnknownError;
+ }
+
+ if( x > x0 )
+ {
+ this->Warning("GetErrorY", "Interpolating error at %f",x0);
+ return (gr->GetErrorY(i)+gr->GetErrorY(i-1))/2;
+ }
+ }
+
+ return 0;
+}
+
Double_t AliLnB2::Rside2Rlong(Double_t pt, Double_t B2, Double_t Cd) const
{
//
return r3;
}
-TGraphErrors* AliLnB2::Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd, Double_t errPt) const
+TGraphErrors* AliLnB2::Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd) const
{
//
// Rside^2*Rlong from B2 value
grR3->SetPoint(i, pt, r3);
// only statistical error propagation of B2
- Double_t errB2 = grR3->GetErrorY(i);
+ Double_t errPt = grB2->GetErrorX(i);
+ Double_t errB2 = grB2->GetErrorY(i);
Double_t errR = r3*errB2/b2;
grR3->SetPointError(i, errPt, errR);
virtual ~AliLnB2();
- TGraphErrors* GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name, Double_t errPt=0) const;
+ TGraphErrors* GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name) const;
Double_t Rside2Rlong(Double_t pt, Double_t B2, Double_t Cd) const;
- TGraphErrors* Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd, Double_t errPt=0) const;
+ TGraphErrors* Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd) const;
Int_t Run();
AliLnB2(const AliLnB2& other);
AliLnB2& operator=(const AliLnB2& other);
+ Double_t GetErrorY(const TGraphErrors* gr, Double_t x0) const;
+
private:
TString fProtonSpectra; // proton spectra filename
#include <TROOT.h>
#include <TFileMerger.h>
-#include "AliLnUnfolding.h"
#include "AliLnSecondaries.h"
#include "AliLnEfficiency.h"
#include "AliLnCorr.h"
// constructor
//
- fUnfolding = new AliLnUnfolding(particle, simuFilename, Form("Unfolding-%s",outputFilename.Data()), otag);
fSecondaries = new AliLnSecondaries(particle, dataFilename, simuFilename, Form("Secondaries-%s",outputFilename.Data()), otag);
fEfficiency = new AliLnEfficiency(particle, simuFixFilename, Form("Efficiency-%s",outputFilename.Data()), otag); // simufix for efficiencies
}
//
// destructor
//
- delete fUnfolding;
delete fSecondaries;
delete fEfficiency;
}
//
// rebuild correction file
//
- fUnfolding->Exec();
fSecondaries->Exec();
fEfficiency->Exec();
// merge the root files
- TString output1 = *fUnfolding->GetOutputFilename();
TString output3 = *fSecondaries->GetOutputFilename();
TString output4 = *fEfficiency->GetOutputFilename();
TFileMerger m;
- m.AddFile(output1.Data(),0);
m.AddFile(output3.Data(),0);
m.AddFile(output4.Data(),0);
m.Merge();
// remove tmp files
- gSystem->Exec(Form("rm -f %s %s %s", output1.Data(), output3.Data(), output4.Data()));
+ gSystem->Exec(Form("rm -f %s %s", output3.Data(), output4.Data()));
gSystem->Exec(Form("mv debug-%s debug-%s",output3.Data(),fOutputFilename.Data()));
#include <TObject.h>
#include <TString.h>
-class AliLnUnfolding;
class AliLnSecondaries;
class AliLnEfficiency;
Int_t Exec();
- AliLnUnfolding* GetLnUnfolding() { return fUnfolding; }
AliLnSecondaries* GetLnSecondaries() { return fSecondaries; }
AliLnEfficiency* GetLnEfficiency() { return fEfficiency; }
TString fOutputFilename; // output filename
- AliLnUnfolding* fUnfolding; // unfolding correction
AliLnSecondaries* fSecondaries; // secondaries correction
AliLnEfficiency* fEfficiency; // eficiency correction
#include "AliLnDriver.h"
#include "B2.h"
+
ClassImp(AliLnDriver)
AliLnDriver::AliLnDriver()
, fMakeRatio(1)
, fMakeSpectra(1)
, fMakeStats(1)
-, fLowPtBin(3)
-, fHiPtBin(15)
-, fLowM2Bin(9)
-, fHighM2Bin(17)
-, fUnfolding(0)
-, fNIter(4)
+, fPtMin(3)
+, fPtMax(15)
+, fPid(0.4)
+, fPidPt(1.)
, fSecondaries(1)
-, fSecProd(AliLnSecondaries::kTFractionFitter)
+, fSecProc(AliLnSecondaries::kTFractionFitter)
, fMatDCAxyMod(AliLnSecondaries::kGeantDCAxy)
, fANucTemplate(0)
, fNbin(1)
, fYMax(0.5)
, fMinDCAxy(-1.5)
, fMaxDCAxy(1.5)
-, fMinM2Bkg(2.2)
-, fMaxM2Bkg(5.)
-, fMinM2tpc(2.)
-, fMaxM2tpc(6.5)
+, fBkgMin(2.2)
+, fBkgMax(5.)
+, fIntMin(2.)
+, fIntMax(6.5)
, fEfficiency(1)
, fG3Fluka(0)
, fScMat(1)
, fScFd(1)
-, fSysPos(1)
-, fSysNeg(1)
, fInputData()
, fInputSimu()
, fInputSimuFix()
, fSameFdwn(0)
, fMCtoINEL(0)
, fAddFakeTracks(1)
+, fPidProc(AliLnPt::kMassSquared)
+, fPidEff(1)
+, fDebugLevel(0)
{
//
// constructor
for(Int_t i=0; i<3; ++i)
{
fTrigEff[i] = 1;
- fXsec[i] = 1;
}
gErrorIgnoreLevel = kWarning;
fOutputPtDebug.Replace(pt.Length()-5,5,"-debug.root");
}
-void AliLnDriver::PrintFilenames() const
-{
-//
-// print filenames to stdout
-//
- using namespace std;
-
- cout << endl;
- cout << "input data : " << fInputData << endl;
- cout << "input simulation : " << fInputSimu << endl;
- cout << "input simu. fix : " << fInputSimuFix << endl;
-
- cout << "corrections : " << fOutputPtCorr << endl;
- cout << "corrections debug : " << fOutputPtCorrDebug << endl;
-
- cout << "output pt : " << fOutputPt << endl;
- cout << "output pt debug : " << fOutputPtDebug << endl;
- cout << "output ratio : " << fOutputRatio << endl;
- cout << "output spectra : " << fOutputSpectra << endl;
- cout << endl;
-}
-
void AliLnDriver::SetTriggerEfficiency(Double_t eff[3])
{
//
for(Int_t i=0; i<3; ++i) fTrigEff[i] = eff[i];
}
-void AliLnDriver::SetInelXSection(Double_t xsec[3])
-{
-//
-// set total cross section, stat. and syst. errors
-//
- for(Int_t i=0; i<3; ++i) fXsec[i] = xsec[i];
-}
-
Int_t AliLnDriver::Run() const
{
//
// run script
//
- if(!fIsOnlyGen && fMakeCorr) this->MakePtCorr();
+ using namespace std;
- if(fMakePt) this->MakePt();
+ if(!fIsOnlyGen && fMakeCorr)
+ {
+ if(fDebugLevel > 0 )
+ {
+ cout << endl;
+ cout << "********* Rebuild corrections ***********" << endl;
+ cout << "Species : " << fSpecies << endl;
+ cout << "Simulation file : " << fInputSimu << endl;
+ cout << "Simulation file (fix) : " << fInputSimuFix << endl;
+ cout << "Output file : " << fOutputPtCorr << endl;
+ cout << "Output file (debug) : " << fOutputPtCorrDebug << endl;
+ cout << endl;
+ }
+
+ this->MakePtCorr();
+ }
+
+ if(fMakePt)
+ {
+ if(fDebugLevel > 0 )
+ {
+ cout << endl;
+ cout << "********* Make pt ***********" << endl;
+ cout << "Species : " << fSpecies << endl;
+ cout << "pt interval (GeV/c) : (" << fPtMin << ", " << fPtMax << ")" << endl;
+ cout << "pid pt (GeV/c) : " << fPidPt << endl;
+ cout << "Input file : " << fInputData << endl;
+ cout << "Correction file : " << fOutputPtCorr << endl;
+ cout << "Output file : " << fOutputPt << endl;
+ cout << "Output file (debug) : " << fOutputPtDebug << endl;
+ cout << endl;
+ }
+
+ this->MakePt();
+ }
- if(fMakeRatio) this->MakeRatio();
+ if(fMakeRatio)
+ {
+ if(fDebugLevel > 0 )
+ {
+ cout << endl;
+ cout << "********* Make ratio ***********" << endl;
+ cout << "Species : " << fSpecies << endl;
+ cout << "Input file : " << fOutputPt << endl;
+ cout << "Output file : " << fOutputRatio << endl;
+ cout << endl;
+ }
+
+ this->MakeRatio();
+ }
- if(fMakeSpectra) this->MakeSpectra();
+ if(fMakeSpectra)
+ {
+ if(fDebugLevel > 0 )
+ {
+ cout << endl;
+ cout << "********* Make spectra ***********" << endl;
+ cout << "Species : " << fSpecies << endl;
+ cout << "INEL extrapolation : " << fINEL << endl;
+ cout << "Rapidity interval : (" << fYMin << ", " << fYMax << ")" << endl;
+ cout << "Input file : " << fOutputPt << endl;
+ cout << "Output file : " << fOutputSpectra << endl;
+ cout << endl;
+ }
+
+ this->MakeSpectra();
+ }
return 0;
}
AliLnCorr lncorr(kParticle[i], fInputData, fInputSimu, fInputSimuFix, outputfile1, fOutputCorTag);
- lncorr.GetLnSecondaries()->SetCorBins(fLowPtBin, fHiPtBin);
- lncorr.GetLnSecondaries()->SetProcedure(fSecProd);
+ lncorr.GetLnSecondaries()->SetCorBins(fPtMin, fPtMax);
+ lncorr.GetLnSecondaries()->SetProcedure(fSecProc);
lncorr.GetLnSecondaries()->SetMatDCAxyModel(fMatDCAxyMod);
lncorr.GetLnSecondaries()->SetAntiNucleusAsTemplate(fANucTemplate);
lncorr.GetLnSecondaries()->SetDCAxyInterval(fMinDCAxy, fMaxDCAxy);
AliLnPt lnpt(kParticle[i], fTrigEff[0], fInputData, ptfile, fOutputTag, fOutputPtCorr, fOutputCorTag);
lnpt.SetOnlyGeneration(fIsOnlyGen);
- lnpt.SetRapidityInterval(fYMin, fYMax);
- lnpt.SetPtBinInterval(fLowPtBin, fHiPtBin);
- lnpt.SetM2BinInterval(fLowM2Bin, fHighM2Bin);
- lnpt.SetM2BkgInterval(fMinM2Bkg, fMaxM2Bkg);
- lnpt.SetM2TPCInterval(fMinM2tpc, fMaxM2tpc);
- lnpt.SetPidM2(fPidM2);
- lnpt.SetUnfolding(fUnfolding, fNIter);
+ lnpt.SetPtInterval(fPtMin, fPtMax);
+ lnpt.SetPid(fPid);
+ lnpt.SetPidProcedure(fPidProc);
+ lnpt.SetPidPt(fPidPt);
+ lnpt.SetBkgInterval(fBkgMin, fBkgMax);
+ lnpt.SetPidInterval(fIntMin, fIntMax);
lnpt.SetSecondaries(fSecondaries);
lnpt.SetEfficiency(fEfficiency);
lnpt.SetMakeStats(fMakeStats);
lnpt.SetFitFractionCorr(fFitFrac);
lnpt.SetFeedDownCorr(fFdwnCorr);
lnpt.SetSameFeedDownCorr(fSameFdwn);
+ lnpt.SetPidEfficiency(fPidEff);
+ lnpt.SetDebugLevel(fDebugLevel);
lnpt.Exec();
//
TFileMerger m;
- const TString kParticle[] = { fSpecies, Form("Anti%s",fSpecies.Data())};
- const Double_t kSysErr[] = { fSysPos, fSysNeg };
+ const TString kParticle[2] = { fSpecies, Form("Anti%s",fSpecies.Data())};
for(Int_t i=0; i<2; ++i)
{
TString spectrafile = kParticle[i] + "-Spectra.root";
- AliLnSpectra lnspectra(kParticle[i], fOutputPt, fOutputTag, spectrafile, fOutputTag, fXsec);
+ AliLnSpectra lnspectra(kParticle[i], fOutputPt, fOutputTag, spectrafile, fOutputTag);
lnspectra.SetRapidityInterval(fYMin, fYMax);
lnspectra.SetExtrapolateToINEL(fINEL);
- lnspectra.SetOnlyGeneration(fIsOnlyGen);
- lnspectra.SetScalingFactor(kSysErr[i]);
lnspectra.Exec();
void MakeRatio() const;
void MakeSpectra() const;
- void PrintFilenames() const;
-
Int_t Run() const;
TString GetSpecies() const { return fSpecies; }
void SetOutputCorTag(const TString& tag) { fOutputCorTag = tag; }
void SetTriggerEfficiency(Double_t eff[3]);
- void SetInelXSection(Double_t xsec[3]);
void SetExtrapolateToINEL(Bool_t flag=1) { fINEL = flag; }
void SetOnlyGeneration(Bool_t flag=1) { fIsOnlyGen = flag; }
void SetRapidityInterval(Double_t ymin, Double_t ymax) { fYMin = ymin; fYMax = ymax; }
- void SetPtBinInterval(Int_t lowbin, Int_t hibin) { fLowPtBin = lowbin; fHiPtBin = hibin; }
- void SetM2BinInterval(Int_t lowbin, Int_t hibin) { fLowM2Bin = lowbin; fHighM2Bin = hibin; }
- void SetM2BkgInterval(Double_t min, Double_t max) { fMinM2Bkg = min; fMaxM2Bkg = max; }
- void SetM2TPCInterval(Double_t min, Double_t max) { fMinM2tpc = min; fMaxM2tpc = max; }
- void SetPidM2(Bool_t flag=1) { fPidM2 = flag; }
- void SetUnfolding(Bool_t flag=1, Int_t niter=4) { fUnfolding = flag; fNIter=niter; }
+ void SetPtInterval(Double_t min, Double_t max) { fPtMin = min; fPtMax = max; }
+ void SetPidPt(Double_t ptpid) { fPidPt = ptpid; }
+ void SetBkgInterval(Double_t min, Double_t max) { fBkgMin = min; fBkgMax = max; }
+ void SetPidInterval(Double_t min, Double_t max) { fIntMin = min; fIntMax = max; }
+ void SetPid(Bool_t flag=1) { fPid = flag; }
void SetSecondaries(Bool_t flag=1) { fSecondaries = flag; }
- void SetSecProd(Int_t prod) { fSecProd = prod; }
+ void SetSecProcedure(Int_t proc) { fSecProc = proc; }
void SetAntiNucleusAsTemplate(Bool_t flag=1) { fANucTemplate = flag; }
void SetMatDCAxyModel(Int_t model=1) { fMatDCAxyMod = model; }
void SetNBin(Int_t nbin) { fNbin = nbin; }
void SetFeedDownCorr(Bool_t flag=1) { fFdwnCorr=flag; }
void SetSameFeedDownCorr(Bool_t flag=1) { fSameFdwn = flag; }
+ void SetPidProcedure(Int_t proc) { fPidProc=proc; }
+
+ void SetPidEfficiency(Double_t eff) { fPidEff=eff; }
+
void SetAddFakeTracks(Bool_t flag=1) { fAddFakeTracks = flag; }
- void SetSysErr( Double_t pos, Double_t neg) { fSysPos = pos; fSysNeg = neg; }
+ void SetDebugLevel(Int_t level) { fDebugLevel = level; }
private:
Bool_t fMakeSpectra; // make spectra
Bool_t fMakeStats; // make event stats
- Int_t fLowPtBin; // low pt bin
- Int_t fHiPtBin; // high pt bin
- Bool_t fPidM2; // enable m2 pid correction
- Int_t fLowM2Bin; // low m2 bin for pid contamination
- Int_t fHighM2Bin; // high m2 bin for pid contamination
- Bool_t fUnfolding; // unfolding correction
- Int_t fNIter; // number of iterations for Bayesian unfolding
+ Double_t fPtMin; // minimum pt value
+ Double_t fPtMax; // maximum pt value
+ Bool_t fPid; // enable pid correction
+ Double_t fPidPt; // minimum pt value for pid correction
Bool_t fSecondaries; // correction of secondaries
- Int_t fSecProd; // procedure for estimating fractions
+ Int_t fSecProc; // procedure to estimate fractions
Int_t fMatDCAxyMod; // DCAxy model for correction of secondaries
Bool_t fANucTemplate; // enable antinucleus as template for primaries
Int_t fNbin; // rebin of DCAxy distribution
Double_t fYMax; // max rapidity
Double_t fMinDCAxy; // min DCAxy
Double_t fMaxDCAxy; // max DCAxy
- Double_t fMinM2Bkg; // min M2 for removing background
- Double_t fMaxM2Bkg; // max M2 for removing background
- Double_t fMinM2tpc; // min M2 for integration
- Double_t fMaxM2tpc; // max M2 for integration
+ Double_t fBkgMin; // lower limit for removing background
+ Double_t fBkgMax; // upper limit for removing background
+ Double_t fIntMin; // lower limit for integration
+ Double_t fIntMax; // upper limit for integration
Bool_t fEfficiency; // efficiency correction
Bool_t fG3Fluka; // enable G3/Fluka correction for TPC
Double_t fScMat; // scaling factor for material fraction
Double_t fScFd; // scaling factor for feed-down fraction
- Double_t fSysPos; // variation for positives
- Double_t fSysNeg; // variation for negatives
-
TString fInputData; // input data filename
TString fInputSimu; // input simulation filename
TString fInputSimuFix; // input fixed simulation filename
Bool_t fAddFakeTracks; // include fake tracks in the efficiency and templates
+ Int_t fPidProc; // pid procedure on the pt distribution
+
+ Double_t fPidEff; // pid efficiency for all pt
+
+ Int_t fDebugLevel; // 0 no verbose, > 1 verbose
+
ClassDef(AliLnDriver,3)
};
foutput->cd(fOutputTag.Data());
}
- TH1D* hGenPhSpPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_PhS_Prim_Pt");
- TH1D* hGenTrigPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Trig_Prim_Pt");
- TH1D* hGenVtxPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Vtx_Prim_Pt");
- TH1D* hGenAccPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Acc_Prim_Pt");
- TH1D* hPrimRecPt = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Prim_Pt");
- TH1D* hFakePrimRecPt = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Fake_Prim_Pt");
+ TH1D* hGenPhSpPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_PhS_Prim_Pt");
+ TH1D* hGenTrigPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_Trig_Prim_Pt");
+ TH1D* hGenVtxPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_Vtx_Prim_Pt");
+ TH1D* hGenAccPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_Acc_Prim_Pt");
+ TH1D* hPrimRecPt = FindObj<TH1D>(fsimu, fParticle + "_Sim_Prim_Pt");
+ TH1D* hFakePrimRecPt = FindObj<TH1D>(fsimu, fParticle + "_Sim_Fake_Prim_Pt");
if(fAddFakeTracks) hPrimRecPt->Add(hFakePrimRecPt);
if(!fParticle.Contains("Anti"))
{
- TH1D* hStats = (TH1D*)FindObj(fsimu, fParticle + "_Stats");
+ TH1D* hStats = FindObj<TH1D>(fsimu, fParticle + "_Stats");
hStats->Write();
}
--- /dev/null
+/*****************************************************************************
+ * Project: RooFit *
+ * *
+ * This code was autogenerated by RooClassFactory *
+ *****************************************************************************/
+
+// Your description goes here...
+
+#include "Riostream.h"
+
+#include "AliLnGaussianExpTail.h"
+#include "RooAbsReal.h"
+#include "RooAbsCategory.h"
+#include <math.h>
+#include "TMath.h"
+
+ClassImp(AliLnGaussianExpTail)
+
+ AliLnGaussianExpTail::AliLnGaussianExpTail(const char *name, const char *title,
+ RooAbsReal& _x,
+ RooAbsReal& _mu,
+ RooAbsReal& _sigma,
+ RooAbsReal& _tau) :
+ RooAbsPdf(name,title),
+ x("x","x",this,_x),
+ mu("mu","mu",this,_mu),
+ sigma("sigma","sigma",this,_sigma),
+ tau("tau","tau",this,_tau)
+ {
+ }
+
+
+ AliLnGaussianExpTail::AliLnGaussianExpTail(const AliLnGaussianExpTail& other, const char* name) :
+ RooAbsPdf(other,name),
+ x("x",this,other.x),
+ mu("mu",this,other.mu),
+ sigma("sigma",this,other.sigma),
+ tau("tau",this,other.tau)
+ {
+ }
+
+
+
+ Double_t AliLnGaussianExpTail::evaluate() const
+ {
+// from pi/K/p note
+
+ if(x <= mu + tau)
+ {
+ return exp(-0.5*(x-mu)*(x-mu)/sigma/sigma);
+ }
+
+ return exp(-0.5*tau*tau/sigma/sigma)*exp(-tau*(x-mu-tau)/sigma/sigma);
+ }
+
+
+
--- /dev/null
+/*****************************************************************************
+ * Project: RooFit *
+ * *
+ * This code was autogenerated by RooClassFactory *
+ *****************************************************************************/
+
+#ifndef ALILNGAUSSIANEXPTAIL
+#define ALILNGAUSSIANEXPTAIL
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+#include "RooCategoryProxy.h"
+#include "RooAbsReal.h"
+#include "RooAbsCategory.h"
+
+class AliLnGaussianExpTail : public RooAbsPdf {
+public:
+ AliLnGaussianExpTail() {} ;
+ AliLnGaussianExpTail(const char *name, const char *title,
+ RooAbsReal& _x,
+ RooAbsReal& _mu,
+ RooAbsReal& _sigma,
+ RooAbsReal& _tau);
+ AliLnGaussianExpTail(const AliLnGaussianExpTail& other, const char* name=0) ;
+ virtual TObject* clone(const char* newname) const { return new AliLnGaussianExpTail(*this,newname); }
+ inline virtual ~AliLnGaussianExpTail() { }
+
+protected:
+
+ RooRealProxy x ;
+ RooRealProxy mu ;
+ RooRealProxy sigma ;
+ RooRealProxy tau ;
+
+ Double_t evaluate() const ;
+
+private:
+
+ ClassDef(AliLnGaussianExpTail,1) // Your description goes here...
+};
+
+#endif
--- /dev/null
+/*****************************************************************************
+ * Project: RooFit *
+ * *
+ * This code was autogenerated by RooClassFactory *
+ *****************************************************************************/
+
+// Your description goes here...
+
+#include "Riostream.h"
+
+#include "AliLnM2.h"
+#include "RooAbsReal.h"
+#include "RooAbsCategory.h"
+#include <math.h>
+#include "TMath.h"
+
+ClassImp(AliLnM2)
+
+ AliLnM2::AliLnM2(const char *name, const char *title,
+ RooAbsReal& _x,
+ RooAbsReal& _mu,
+ RooAbsReal& _sigma,
+ RooAbsReal& _tau,
+ RooAbsReal& _p) :
+ RooAbsPdf(name,title),
+ x("x","x",this,_x),
+ mu("mu","mu",this,_mu),
+ sigma("sigma","sigma",this,_sigma),
+ tau("tau","tau",this,_tau),
+ p("p","p",this,_p)
+ {
+ }
+
+
+ AliLnM2::AliLnM2(const AliLnM2& other, const char* name) :
+ RooAbsPdf(other,name),
+ x("x",this,other.x),
+ mu("mu",this,other.mu),
+ sigma("sigma",this,other.sigma),
+ tau("tau",this,other.tau),
+ p("p",this,other.p)
+ {
+ }
+
+Double_t AliLnM2::evaluate() const
+{
+//
+// m^2 shape
+//
+ if( x <= mu + tau)
+ {
+ return TMath::Exp(-(x-mu)*(x-mu)/(2.*sigma*sigma));
+ }
+
+ Double_t lambda = 2.*tau*TMath::Sqrt(p*p+mu+tau)/(sigma*sigma);
+ Double_t a = TMath::Exp(-tau*tau/(2.*sigma*sigma)+2.*tau*(p*p+mu+tau)/(sigma*sigma));
+
+ return a*TMath::Exp(-lambda*TMath::Sqrt(p*p + x));
+}
--- /dev/null
+/*****************************************************************************
+ * Project: RooFit *
+ * *
+ * This code was autogenerated by RooClassFactory *
+ *****************************************************************************/
+
+#ifndef ALILNM2
+#define ALILNM2
+
+#include "RooAbsPdf.h"
+#include "RooRealProxy.h"
+#include "RooCategoryProxy.h"
+#include "RooAbsReal.h"
+#include "RooAbsCategory.h"
+
+class AliLnM2 : public RooAbsPdf {
+public:
+ AliLnM2() {} ;
+ AliLnM2(const char *name, const char *title,
+ RooAbsReal& _x,
+ RooAbsReal& _mu,
+ RooAbsReal& _sigma,
+ RooAbsReal& _tau,
+ RooAbsReal& _p);
+ AliLnM2(const AliLnM2& other, const char* name=0) ;
+ virtual TObject* clone(const char* newname) const { return new AliLnM2(*this,newname); }
+ inline virtual ~AliLnM2() { }
+
+protected:
+
+ RooRealProxy x ;
+ RooRealProxy mu ;
+ RooRealProxy sigma ;
+ RooRealProxy tau ;
+ RooRealProxy p ;
+
+ Double_t evaluate() const ;
+
+private:
+
+ ClassDef(AliLnM2,1) // Your description goes here...
+};
+
+#endif
#include <RooDataHist.h>
#include <RooHistPdf.h>
#include <RooMsgService.h>
+#include <TMath.h>
#include "AliLnPt.h"
#include "B2.h"
-#ifdef HAVE_ROOUNFOLD
- #include "/opt/RooUnfold/src/RooUnfoldResponse.h"
- #include "/opt/RooUnfold/src/RooUnfoldBayes.h"
-#endif
-
ClassImp(AliLnPt)
AliLnPt::AliLnPt(const TString& particle, Double_t trigEff, const TString& inputFilename, const TString& outputFilename, const TString& otag, const TString& corrFilename, const TString& corrtag)
, fOutputTag(otag)
, fCorrFilename(corrFilename)
, fCorrTag(corrtag)
-, fLowPtBin(1)
-, fHiPtBin(19)
-, fPidM2(0)
-, fUnfolding(0)
-, fNIter(4)
+, fPtMin(0.5)
+, fPtMax(3.0)
+, fPid(0)
, fSecondaries(1)
, fEfficiency(1)
, fIsOnlyGen(0)
-, fYMin(-0.5)
-, fYMax(0.5)
-, fLowM2Bin(9)
-, fHiM2Bin(17)
-, fMinM2Bkg(2.2)
-, fMaxM2Bkg(5.)
-, fMinM2tpc(2.)
-, fMaxM2tpc(6.)
+, fPtPid(1.)
+, fBkgMin(2.2)
+, fBkgMax(5.)
+, fIntMin(2.)
+, fIntMax(6.)
, fMakeStats(1)
, fMCtoINEL(0)
, fVtxFactor(1)
, fFitFrac(1)
, fFdwnCorr(1)
, fSameFdwn(0)
+, fPidProc(kMassSquared)
+, fPidEff(1)
+, fDebugLevel(0)
{
//
// constructor
//
TH1::SetDefaultSumw2(); // switch on histogram errors
- // disable verbose in RooFit
- RooMsgService::instance().setSilentMode(1);
- RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+ if(fDebugLevel < 3)
+ {
+ // disable verbose in RooFit
+ RooMsgService::instance().setSilentMode(1);
+ RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+ }
}
AliLnPt::~AliLnPt()
//
// Get the true pt distribution
//
- using namespace std;
-
TFile* finput = new TFile(fInputFilename.Data(), "read");
if (finput->IsZombie()) exit(1);
TString species = fParticle;
species.ReplaceAll("Anti","");
- TH1D* hStats = (TH1D*)FindObj(finput, species + "_Stats");
+ TH1D* hStats = FindObj<TH1D>(finput, species + "_Stats");
- TH1D* hCorPt = 0;
+ TH1D* hCorrPt = 0;
if(fIsOnlyGen)
{
- hCorPt = (TH1D*)((TH1D*)FindObj(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data()));
+ hCorrPt = dynamic_cast<TH1D*>( (FindObj<TH1D>(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data())) );
}
- else // correct the reconstructed pt
+ else
{
fcorr = new TFile(fCorrFilename.Data(),"read");
if (fcorr->IsZombie()) exit(1);
fdebug->mkdir(fOutputTag.Data());
- // -------------- get pointers to the required histograms ---------
+ // pointers to the required histograms
- // pt distributions
- TH1D* hOrigPt = (TH1D*)FindObj(finput, fParticle + "_PID_Pt");
- TH2D* hM2Pt = (TH2D*)FindObj(finput, fParticle + "_PID_M2_Pt");
+ TH1D* hOrigPt = FindObj<TH1D>(finput, fParticle + "_PID_Pt");
- // corrections
-#ifdef HAVE_ROOUNFOLD
- TH2D* hResponseMtx = (TH2D*)FindObj(fcorr, fCorrTag, fParticle + "_Response_Matrix");
- TH1D* hMeasuredPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Measured_Pt");
- TH1D* hTruePt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_True_Pt");
-#endif
+ TH2D* hPidPt = 0;
+ if(fPidProc == kTimeOfFlight) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DTime_Pt");
+ else if(fPidProc == kMassSquared) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_M2_Pt");
+ else if(fPidProc == kMassDifference) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DM2_Pt");
+ else fPid = 0;
- TH1D* hFracMatPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt");
- TH1D* hFracFdwnPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt");
- TF1* fncFracMatPt = (TF1*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt");
- TF1* fncFracFdwnPt = (TF1*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Fit_Pt");
+ // corrections
+ TH1D* hFracMatPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt");
+ TH1D* hFracFdwnPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt");
+ TF1* fncFracMatPt = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt");
+ TF1* fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Fit_Pt");
if(fSameFdwn)
{
- hFracFdwnPt = (TH1D*)FindObj(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
- fncFracFdwnPt = (TF1*)FindObj(fcorr, fCorrTag, species + "_Frac_Fdwn_Fit_Pt");
+ hFracFdwnPt = FindObj<TH1D>(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
+ fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, species + "_Frac_Fdwn_Fit_Pt");
}
- TH1D* hEffVtxPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt");
- TH1D* hEffAccTrkPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Eff_AccTrk_Pt");
+ TH1D* hEffVtxPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt");
+ TH1D* hEffAccTrkPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_AccTrk_Pt");
- // -------------- apply corrections -----------
+ // apply corrections
fdebug->cd(fOutputTag.Data());
- TH1D* hPidPt = 0;
- if(fPidM2 && (fHiPtBin > fLowM2Bin))
+ TH1D* hPidCorrPt = 0;
+ if(fPid)
{
- Int_t hibin = (fHiPtBin > fHiM2Bin) ? fHiM2Bin : fHiPtBin;
- hPidPt = this->PID(hOrigPt, hM2Pt, fLowM2Bin, hibin, fParticle + "_M2Corr_Pt");
+ if(fDebugLevel>1) std::cout << "\t\tpid correction" << std::endl;
+
+ hPidCorrPt = this->PID(hOrigPt, hPidPt, fPtPid, fPtMax, fParticle + "_PidCorr_Pt");
+ hPidCorrPt->Scale(1./fPidEff);
}
else
{
- hPidPt = (TH1D*)hOrigPt->Clone(Form("%s_M2Corr_Pt", fParticle.Data()));
+ hPidCorrPt = dynamic_cast<TH1D*>(hOrigPt->Clone(Form("%s_PidCorr_Pt", fParticle.Data())));
}
- TH1D* hSecCorPt = 0;
+ TH1D* hSecCorrPt = 0;
if(fSecondaries)
{
+ if(fDebugLevel>1) std::cout << "\t\tremoval of secondaries" << std::endl;
+
if(fFitFrac)
{
- hSecCorPt = this->Secondaries(hPidPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCor_Pt");
+ hSecCorrPt = this->Secondaries(hPidCorrPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCorr_Pt");
}
else
{
- hSecCorPt = this->Secondaries(hPidPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCor_Pt");
+ hSecCorrPt = this->Secondaries(hPidCorrPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCorr_Pt");
}
}
else
{
- hSecCorPt = (TH1D*)hPidPt->Clone(Form("%s_SecCor_Pt", fParticle.Data()));
+ hSecCorrPt = dynamic_cast<TH1D*>(hPidCorrPt->Clone(Form("%s_SecCorr_Pt", fParticle.Data())));
}
- TH1D* hUnfoldedPt = 0;
-#ifdef HAVE_ROOUNFOLD
- if(fUnfolding)
- {
- hUnfoldedPt = this->Unfolding(hSecCorPt, hMeasuredPt, hTruePt, hResponseMtx, fParticle + "_Unfolded_Pt", fNIter);
- }
- else
-#endif
- {
- hUnfoldedPt = (TH1D*)hSecCorPt->Clone(Form("%s_Unfolded_Pt", fParticle.Data()));
- }
-
- TH1D* hEffCorPt = 0;
+ TH1D* hEffCorrPt = 0;
if(fEfficiency)
{
+ if(fDebugLevel>1) std::cout << "\t\tefficiency correction" << std::endl;
+
if(fMCtoINEL)
{
- TH1D* hStatsMC = (TH1D*)FindObj(fcorr, species + "_Stats");
+ TH1D* hStatsMC = FindObj<TH1D>(fcorr, fCorrTag, species + "_Stats");
fVtxFactor = this->GetVertexCorrection(hStats, hStatsMC);
}
- hEffCorPt = this->Efficiency(hUnfoldedPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCor_Pt");
+ hEffCorrPt = this->Efficiency(hSecCorrPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCorr_Pt");
}
else
{
- hEffCorPt = (TH1D*)hUnfoldedPt->Clone(Form("%s_EffCor_Pt", fParticle.Data()));
+ hEffCorrPt = dynamic_cast<TH1D*>(hSecCorrPt->Clone(Form("%s_EffCorr_Pt", fParticle.Data())));
}
- // --------- corrected pt --------
+ // corrected pt
- hCorPt = (TH1D*) hEffCorPt->Clone(Form("%s_Cor_Pt", fParticle.Data()));
+ hCorrPt = dynamic_cast<TH1D*>(hEffCorrPt->Clone(Form("%s_Corr_Pt", fParticle.Data())));
- // ---------- debug -----
+ // debug histograms
fdebug->cd(fOutputTag.Data());
hOrigPt->Write(); // original distribution
- hPidPt->Write();
- hUnfoldedPt->Write();
- hSecCorPt->Write();
- hEffCorPt->Write();
- hCorPt->Write();
-
- // ------------- clean --------------
-
- delete hPidPt;
- delete hUnfoldedPt;
- delete hSecCorPt;
- delete hEffCorPt;
+ hPidCorrPt->Write();
+ hSecCorrPt->Write();
+ hEffCorrPt->Write();
+ hCorrPt->Write();
+
+ // free memory
+ delete hPidCorrPt;
+ delete hSecCorrPt;
+ delete hEffCorrPt;
}
- // ---------- save the final pt ----------
+ // save final pt
- TH1D* hPt = (TH1D*) hCorPt->Clone(Form("%s_Pt", fParticle.Data()));
+ TH1D* hPt = dynamic_cast<TH1D*>(hCorrPt->Clone(Form("%s_Pt", fParticle.Data())));
hPt->SetTitle(fParticle.Data());
hPt->Reset();
hPt->Sumw2();
- for(Int_t i=fLowPtBin; i<fHiPtBin; ++i)
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(fPtMin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(fPtMax);
+
+ for(Int_t i=lowbin; i<hibin; ++i)
{
- hPt->SetBinContent(i, hCorPt->GetBinContent(i));
- hPt->SetBinError(i, hCorPt->GetBinError(i));
+ hPt->SetBinContent(i, hCorrPt->GetBinContent(i));
+ hPt->SetBinError(i, hCorrPt->GetBinError(i));
}
foutput->cd();
hPt->Write();
delete hPt;
- delete hCorPt;
+ delete hCorrPt;
delete fdebug;
delete fcorr;
- // --------- stats -------------------
+ // stats
if(fMakeStats)
{
Double_t AliLnPt::GetM2Width(Double_t pt, Double_t m) const
{
//
-// parameterization of the expected M2 width
-// as a function of pt from LHC12a5b
+// expected M2 width parameterization from LHC12a5b
//
Double_t c0 = 1.90212e-03;
Double_t c1 = 1.29814e-02;
return c0/(pt*pt) + c1*(pt*pt/(m*m) + 1.);
}
-RooWorkspace* AliLnPt::GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const
+Double_t AliLnPt::GetExpectedDT(Double_t pt, Double_t m) const
{
//
-// deuteron m2 model for TPC pid
-// gaussian + exponential background of pi, K and p
+// estimate of the expected time of flight (ns) for mass hypothesis m
+//
+ Double_t l = 370.; // cm
+ Double_t c = 2.99792458e+1; // cm/ns
+
+ return (l/c)*TMath::Sqrt(1.+(m/pt)*(m/pt));
+}
+
+RooWorkspace* AliLnPt::GetToFModel(Double_t pt, const TString& name) const
+{
+//
+// model for the time of flight distribution
//
using namespace RooFit;
RooWorkspace* w = new RooWorkspace(name.Data());
- RooRealVar m2("m2", "m^{2} (GeV^{2}/c^{4})", fMinM2Bkg, fMaxM2Bkg);
- w->import(m2);
+ // signal
+ w->factory(Form("AliLnGaussianExpTail::Sd(x[%f,%f], mu[0,-0.06,0.06], sigma[0.120,0.05,0.3], tau[0.1,0.08,0.4])", fBkgMin, fBkgMax));
- Double_t sigma = GetM2Width(pt, m);
+ // background
+ Double_t expTd = this->GetExpectedDT(pt, GetMass("deuteron"));
+ Double_t expDTp = this->GetExpectedDT(pt, GetMass("proton")) - expTd;
+ Double_t expDTk = this->GetExpectedDT(pt, GetMass("kaon")) - expTd;
+ Double_t expDTpi = this->GetExpectedDT(pt, GetMass("pion")) - expTd;
+
+ w->factory(Form("AliLnGaussianExpTail::ProtonBkg(x, muP[%f,%f-0.2,%f+0.2], widthP[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTp, expDTp, expDTp));
+ w->factory(Form("AliLnGaussianExpTail::KaonBkg(x, muK[%f,%f-0.2,%f+0.2], widthK[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTk, expDTk, expDTk));
+ w->factory(Form("AliLnGaussianExpTail::PionBkg(x, muPi[%f,%f-0.2,%f+0.2], widthPi[0.120,0.05,0.3], tauPi[0.08,0.01,0.2])",expDTpi, expDTpi, expDTpi));
+
+ w->factory("SUM::Bkg( Npi[0.3,0,1.e+9]*PionBkg, Nk[0.24,0,1.e+9]*KaonBkg, Np[0.47,0,1.e+9]*ProtonBkg, Np[0.47,0,1.e+9]*ProtonBkg)");
+
+ w->factory("SUM::model( Nd[0,1.e+9]*Sd, Nbkg[0,1.e+9]*Bkg )");
+
+ w->var("x")->SetTitle("t - <t> (ns)");
+
+ return w;
+}
+
+RooWorkspace* AliLnPt::GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const
+{
+//
+// mass squared model for deuterons
+//
+ using namespace RooFit;
- // signal (m is shifted to the left)
- w->factory(Form("Gaussian::Sd(m2, mean[%f,%f,%f], width[%f,%f,%f])",m*m, m*m-0.3, m*m+0.05, sigma, sigma-0.05, sigma+0.05));
+ RooWorkspace* w = new RooWorkspace(name.Data());
// background
- w->factory("Exponential::PionBkg(m2,lambda1[-0.5,-100.,0.])");
- w->factory("Exponential::KaonBkg(m2,lambda2[-0.5,-100.,0.])");
- w->factory("Exponential::ProtonBkg(m2,lambda3[-0.5,-100.,0.])");
- w->factory(Form("SUM::Bkg( Npi[0.3,0,%f]*PionBkg, Nk[0.24,0,%f]*KaonBkg, Np[0.47,0,%f]*ProtonBkg )",max,max,max));
+ w->factory(Form("Exponential::PionBkg(x[%f,%f],lambdaPi[-0.3,-10.,0.])", fBkgMin, fBkgMax));
+ w->factory("Exponential::KaonBkg(x,lambdaK[-0.4,-10.,0.])");
+ w->factory("Exponential::ProtonBkg(x,lambdaP[-0.5,-10.,0.])");
+
+ w->factory("SUM::Bkg( Npi[0.3,0.,1.]*PionBkg, Nk[0.24,0.,1.]*KaonBkg, Np[0.47,0.,1.]*ProtonBkg )");
- w->factory(Form("SUM::model( Nd[0,%f]*Sd, Nbkg[0,%f]*Bkg )",max,max));
+ // signal
+
+ Double_t sigma = GetM2Width(pt, m);
+ Double_t mm = (fPidProc == kMassSquared) ? m*m : 0;
+
+ if(pt<2.0)
+ {
+ w->factory(Form("AliLnM2::Sd(x,mu[%f,%f,%f], sigma[%f,%f,%f],tau[0.1,0.05,0.4],p[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1, 1.16*pt,pt,1.33*pt));
+ }
+ else // neglect non-gaussian tails
+ {
+ w->factory(Form("Gaussian::Sd(x, mu[%f,%f,%f], sigma[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1));
+ }
+
+ // model
+ w->factory(Form("SUM::model( Nd[0.,%f]*Sd, Nbkg[0.,%f]*Bkg )", max, max));
+
+ w->var("x")->SetTitle("#it{m}^{2} (GeV^{2}/c^{4})");
return w;
}
-void AliLnPt::GetPtFromM2(TH1D* hPt, Int_t lowbin, Int_t hibin, const TH2D* hM2Pt, Double_t m2min, Double_t m2max) const
+void AliLnPt::GetPtFromPid(TH1D* hPt, Double_t ptmin, Double_t ptmax, const TH2D* hPidPt, Double_t pidMin, Double_t pidMax) const
{
//
-// Fill the pt from lowbin to hibin integrating
-// the m2-pt distribution from m2min to m2max
+// Integrate pid distribution
//
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(ptmax);
+
for(Int_t i=lowbin; i<hibin; ++i)
{
- TH1D* hM2 = hM2Pt->ProjectionY(Form("_M2_%02d",i),i,i);
+ TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
- Int_t m2minbin = hM2->GetXaxis()->FindFixBin(m2min);
- Int_t m2maxbin = hM2->GetXaxis()->FindFixBin(m2max); // m2 has a non-gaussian tail on the right
+ Int_t pidLowBin = hPid->GetXaxis()->FindFixBin(pidMin);
+ Int_t pidHiBin = hPid->GetXaxis()->FindFixBin(pidMax);
Double_t err = 0;
- Double_t n = hM2->IntegralAndError(m2minbin, m2maxbin, err);
+ Double_t n = hPid->IntegralAndError(pidLowBin, pidHiBin, err);
hPt->SetBinContent(i, n);
hPt->SetBinError(i, err);
- delete hM2;
+ delete hPid;
}
}
-TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hM2Pt, Int_t lowbin, Int_t hibin, const TString& name)
+TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hPidPt2, Double_t ptmin, Double_t ptmax, const TString& name)
{
//
-// Remove contamination in the m^2 distribution
+// Remove contamination in the pt distribution
//
+ TH1D* hPidCorrPt = dynamic_cast<TH1D*>(hPt->Clone(name.Data()));
+ if(hPidCorrPt == 0) return 0;
+
+ hPidCorrPt->Reset();
+ hPidCorrPt->Sumw2();
+
+ // set binning
const Int_t kNBin = 2;
- Double_t m = GetMass(fParticle);
- TH1D* hPidPt = (TH1D*) hPt->Clone(name.Data());
- hPidPt->Reset();
- hPidPt->Sumw2();
+ TH2D* hPidPt = dynamic_cast<TH2D*>(hPidPt2->Clone("_pid_pt_"));
+ hPidPt->Rebin2D(1,kNBin);
- // first bins are not contaminated since the identification is clear
- // integrate around the m2 value
+ // integrate around the mean value
- this->GetPtFromM2(hPidPt, 0, lowbin, hM2Pt, fMinM2tpc, fMaxM2tpc);
+ this->GetPtFromPid(hPidCorrPt, 0, ptmin, hPidPt, fIntMin, fIntMax);
- // the remaining bins are contaminated
- // slice the m2-pt and fit each slice with a gaussian + exponentials bkg
+ // slice the m2-pt
using namespace RooFit;
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(ptmax);
+
+ Double_t m = GetMass(fParticle);
+
for(Int_t i=lowbin; i<hibin; ++i)
{
- TH1D* hM2 = hM2Pt->ProjectionY(Form("%s_PID_M2_%02d",fParticle.Data(),i),i,i);
- hM2->Rebin(kNBin);
+ TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
- RooWorkspace* w = GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hM2->Integral());
+ RooWorkspace* w = (fPidProc==kTimeOfFlight) ? GetToFModel( hPt->GetBinCenter(i), Form("%s_M2_%02d",fParticle.Data(),i)) : GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hPid->Integral());
- RooRealVar* m2 = w->var("m2");
- RooDataHist data("data", "data to fit", *m2, hM2);
+ RooRealVar* x = w->var("x");
+ RooDataHist data("data", "data to fit", *x, hPid);
w->import(data);
w->pdf("model")->fitTo(data, Minos(kTRUE), PrintLevel(-1), Verbose(kFALSE), Warnings(kFALSE), PrintEvalErrors(-1) );
- hPidPt->SetBinContent(i, w->var("Nd")->getVal());
- hPidPt->SetBinError(i, w->var("Nd")->getError());
+ hPidCorrPt->SetBinContent(i, w->var("Nd")->getVal());
+ hPidCorrPt->SetBinError(i, w->var("Nd")->getError());
w->Write();
delete w;
- delete hM2;
+ delete hPid;
}
- return hPidPt;
+ delete hPidPt;
+ return hPidCorrPt;
}
-#ifdef HAVE_ROOUNFOLD
-
-TH1D* AliLnPt::Unfolding(const TH1D* hPt, const TH1D* hMeasuredPt, const TH1D* hTruePt, const TH2D* hResponseMtx, const TString& name, Int_t iterations) const
-{
-//
-// Bayesian unfolding
-//
- RooUnfoldResponse* matrix = new RooUnfoldResponse(hMeasuredPt, hTruePt, hResponseMtx);
-
- RooUnfoldBayes* unfold = new RooUnfoldBayes(matrix, hPt, iterations);
-
- TH1D* hUnfoldedPt = (TH1D*)unfold->Hreco();
- hUnfoldedPt->SetName(name.Data());
- hUnfoldedPt->SetXTitle("p_{T} (GeV/c)");
-
- delete unfold;
- delete matrix;
-
- return hUnfoldedPt;
-}
-
-#endif
-
TH1D* AliLnPt::Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* hFracFdwnPt, const TString& name) const
{
//
TH1D* AliLnPt::Secondaries(const TH1D* hPt, TF1* fncMatPt, TF1* fncFdwnPt, const TString& name) const
{
//
-// remove contamination of secondaries using fitted fractions
+// remove contamination of secondaries using a function fitted to the fractions
// N_prim = N/(1+fmat+fdwn)
//
TF1* fnc = new TF1("_secondaries_","1 + [0]*exp(-[1]*x)+[2] + [3]*exp(-[4]*x)+[5]", 0, 10);
#include <TObject.h>
#include <TString.h>
-#ifndef HAVE_ROOUNFOLD
-//#define HAVE_ROOUNFOLD
-#endif
-
class TString;
class TH1D;
class TH2D;
AliLnPt(const TString& particle, Double_t trigEff, const TString& inputFilename, const TString& outputFilename, const TString& otag, const TString& corrFilename, const TString& corrtag);
virtual ~AliLnPt();
- TH1D* PID(const TH1D* hPt, const TH2D* hM2Pt, Int_t lowbin, Int_t hibin, const TString& name);
-#ifdef HAVE_ROOUNFOLD
- TH1D* Unfolding(const TH1D* hPt, const TH1D* hMeasuredPt, const TH1D* hTruePt, const TH2D* ResponseMtx, const TString& name, Int_t iterations=4) const;
-#endif
+ TH1D* PID(const TH1D* hPt, const TH2D* hPidPt, Double_t ptmin, Double_t ptmax, const TString& name);
TH1D* Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* hFracFdwnPt, const TString& name) const;
TH1D* Secondaries(const TH1D* hPt, TF1* fncMatPt, TF1* fncFdwnPt, const TString& name) const;
+
TH1D* Efficiency(const TH1D* hPt, const TH1D* hEffVtxPt, const TH1D* hEffAccTrkPt, const TString& name) const;
Double_t GetVertexCorrection(const TH1D* hData, const TH1D* hMC) const;
Int_t Exec();
- void SetRapidityInterval(Double_t ymin, Double_t ymax) { fYMin = ymin; fYMax = ymax; }
-
void SetParticle(const TString& particle) { fParticle= particle; }
void SetOutputTag(const TString& tag) { fOutputTag = tag; }
void SetCorrectionTag(const TString& tag) { fCorrTag = tag; }
- void SetPtBinInterval(Int_t lowbin, Int_t hibin) { fLowPtBin = lowbin; fHiPtBin = hibin; };
- void SetM2BinInterval(Int_t lowbin, Int_t hibin) { fLowM2Bin = lowbin; fHiM2Bin = hibin; };
+ void SetPtInterval(Double_t min, Double_t max) { fPtMin = min; fPtMax = max; };
- void SetM2BkgInterval(Double_t min, Double_t max) { fMinM2Bkg = min; fMaxM2Bkg = max; };
- void SetM2TPCInterval(Double_t min, Double_t max) { fMinM2tpc = min; fMaxM2tpc = max; };
+ void SetPidProcedure(Int_t proc) { fPidProc=proc; }
+ void SetPidPt(Double_t ptpid) { fPtPid = ptpid; };
+ void SetBkgInterval(Double_t min, Double_t max) { fBkgMin = min; fBkgMax = max; };
+ void SetPidInterval(Double_t min, Double_t max) { fIntMin = min; fIntMax = max; };
- void SetPidM2(Bool_t flag=1) { fPidM2 = flag; }
- void SetUnfolding(Bool_t flag=1, Int_t niter=4) { fUnfolding = flag; fNIter=niter; }
+ void SetPid(Bool_t flag=1) { fPid = flag; }
void SetSecondaries(Bool_t flag=1) { fSecondaries = flag; }
void SetEfficiency(Bool_t flag=1) { fEfficiency = flag; }
void SetFeedDownCorr(Bool_t flag=1) { fFdwnCorr=flag; }
void SetSameFeedDownCorr(Bool_t flag=1) { fSameFdwn=flag; }
+ void SetPidEfficiency(Double_t eff) { fPidEff=eff; }
+
+ void SetDebugLevel(Int_t level) { fDebugLevel = level; }
+
+ enum { kMassSquared=0, kMassDifference, kTimeOfFlight };
private:
AliLnPt(const AliLnPt& other);
AliLnPt& operator=(const AliLnPt& other);
- void GetPtFromM2(TH1D* hPt, Int_t lowbin, Int_t hibin, const TH2D* hM2Pt, Double_t m2min, Double_t m2max) const;
+ void GetPtFromPid(TH1D* hPt, Double_t ptmin, Double_t ptmax, const TH2D* hM2Pt, Double_t m2min, Double_t m2max) const;
Double_t GetM2Width(Double_t pt, Double_t mass) const;
RooWorkspace* GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const;
+ RooWorkspace* GetToFModel(Double_t pt, const TString& name) const;
+ Double_t GetExpectedDT(Double_t pt, Double_t m) const;
private:
TString fCorrFilename; // file with the corrections
TString fCorrTag; // tag for the correction file
- Int_t fLowPtBin; // low bin for the pt in GeV/c/A
- Int_t fHiPtBin; // high bin for the pt in GeV/c/A
+ Double_t fPtMin; // minimum pt value in GeV/c
+ Double_t fPtMax; // maximum pt value in GeV/c
- Bool_t fPidM2; // correct contamination of m2 pid
- Bool_t fUnfolding; // unfold the pt correction
- Int_t fNIter; // number of iterations for Bayesian unfolding
+ Bool_t fPid; // enable additional PID
Bool_t fSecondaries; // remove secondaries
Bool_t fEfficiency; // correct by efficiency
Bool_t fIsOnlyGen; // if the simulation is only generation
- Double_t fYMin; // rapidity interval min y
- Double_t fYMax; // rapidity interval max y
+ Double_t fPtPid; // minimum pt value for pid correction
- Int_t fLowM2Bin; // low pt bin for pid correction
- Int_t fHiM2Bin; // high pt bin for pid correction
+ Double_t fBkgMin; // minimum value to remove background
+ Double_t fBkgMax; // maximum value to remove background
- Double_t fMinM2Bkg; // minimum m2 to remove background
- Double_t fMaxM2Bkg; // maximum m2 to remove background
-
- Double_t fMinM2tpc; // minimum m2 for integration
- Double_t fMaxM2tpc; // maximum m2 for integration
+ Double_t fIntMin; // minimum value for integration
+ Double_t fIntMax; // maximum value for integration
Bool_t fMakeStats; // make event stats
Bool_t fFdwnCorr; // enable feed-down correction
Bool_t fSameFdwn; // same feed-down correction for positives and negatives
+ Int_t fPidProc; // selected pid procedure on the pt distribution
+
+ Double_t fPidEff; // pid efficiency for all pt bins
+
+ Int_t fDebugLevel; // 0 no verbose, > 1 verbose
+
ClassDef(AliLnPt,1)
};
for(Int_t i=0; i<2; ++i)
{
- hPt[i] = (TH1D*)FindObj(finput, kPrefix[i] + fSpecies + "_Pt");
+ hPt[i] = FindObj<TH1D>(finput, kPrefix[i] + fSpecies + "_Pt");
}
TH1D* hRatioPt = Divide(hPt[1],hPt[0],Form("Anti%s%s_Ratio_Pt",fSpecies.Data(),fSpecies.Data()));
#include <TObjArray.h>
#include <TH1D.h>
#include <TH2D.h>
-#include <RooRealVar.h>
-#include <RooDataSet.h>
-#include <RooDataHist.h>
-#include <RooHistPdf.h>
-#include <RooWorkspace.h>
-#include <RooMsgService.h>
#include <TBackCompFitter.h>
#include "AliLnSecondaries.h"
, fSimuFilename(simuFilename)
, fOutputFilename(outputFilename)
, fOutputTag(otag)
-, fLowPtBin(3)
-, fHiPtBin(15)
+, fPtMin(0.5)
+, fPtMax(3.0)
, fNbin(1)
, fMinDCAxy(-1.5)
, fMaxDCAxy(1.5)
// constructor
//
TH1::SetDefaultSumw2(); // switch on histogram errors
-
- // disable verbose in RooFit
- RooMsgService::instance().setSilentMode(1);
- RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
}
AliLnSecondaries::~AliLnSecondaries()
// --------- ideal values: primaries + no secondaries --------
- TH1D* hPidPt = (TH1D*)FindObj(fdata, fParticle + "_PID_Pt");
+ TH1D* hPidPt = FindObj<TH1D>(fdata, fParticle + "_PID_Pt");
hPidPt->Write();
const char* contrib[] = { "prim", "mat", "fdwn" };
if(fFracProc == kMonteCarlo)
{
// compute the fractions from the montecarlo simulation
- TH1D* hAll = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Pt");
- TH1D* hPrim = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Prim_Pt");
- TH1D* hMat = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Mat_Pt");
- TH1D* hFdwn = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Fdwn_Pt");
+ TH1D* hAll = FindObj<TH1D>(fsimu, fParticle + "_Sim_Pt");
+ TH1D* hPrim = FindObj<TH1D>(fsimu, fParticle + "_Sim_Prim_Pt");
+ TH1D* hMat = FindObj<TH1D>(fsimu, fParticle + "_Sim_Mat_Pt");
+ TH1D* hFdwn = FindObj<TH1D>(fsimu, fParticle + "_Sim_Fdwn_Pt");
delete hFracPt[0];
delete hFracPt[1];
else if(fParticle == "AntiProton")
{
// assume only secondaries from feed-down
- TH2D* hDataDCAxyPt = (TH2D*)FindObj(fdata, fParticle + "_PID_DCAxy_Pt");
- TH2D* hDataMCDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_PID_DCAxy_Pt");
- TH2D* hPrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
- TH2D* hFdwnDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
- TH2D* hFakePrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
- TH2D* hFakeFdwnDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
+ TH2D* hDataDCAxyPt = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hDataMCDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hPrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
+ TH2D* hFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
+ TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
+ TH2D* hFakeFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
if(fAddFakeTracks)
{
else if(fParticle == "Proton")
{
// secondaries from materials and feed-down
- TH2D* hDataDCAxyPt = (TH2D*)FindObj(fdata, fParticle + "_PID_DCAxy_Pt");
- TH2D* hDataMCDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_PID_DCAxy_Pt");
- TH2D* hPrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
- TH2D* hMatDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
- TH2D* hFdwnDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
- TH2D* hFakePrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
- TH2D* hFakeMatDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Mat_DCAxy_Pt");
- TH2D* hFakeFdwnDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
+ TH2D* hDataDCAxyPt = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hDataMCDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hPrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
+ TH2D* hMatDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
+ TH2D* hFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
+ TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
+ TH2D* hFakeMatDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Mat_DCAxy_Pt");
+ TH2D* hFakeFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
if(fAddFakeTracks)
{
else
{
// only secondaries from materials for nuclei
- TH2D* hDataDCAxyPt = (TH2D*)FindObj(fdata, fParticle + "_PID_DCAxy_Pt");
- TH2D* hDataMCDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_PID_DCAxy_Pt");
- TH2D* hPrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
- TH2D* hFakePrimDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
+ TH2D* hDataDCAxyPt = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hDataMCDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
+ TH2D* hPrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
+ TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
if(fAddFakeTracks) hPrimDCAxyPt->Add(hFakePrimDCAxyPt);
if(fANucTemplate)
{
- hPrimDCAxyPt = (TH2D*)FindObj(fdata, Form("Anti%s_PID_DCAxy_Pt",fParticle.Data()));
+ hPrimDCAxyPt = FindObj<TH2D>(fdata, Form("Anti%s_PID_DCAxy_Pt",fParticle.Data()));
}
- TH2D* hMatDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
+ TH2D* hMatDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
if(fMatDCAxyMod == kFlatDCAxy)
{
hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
}
else
{
- hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.4);
+ hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.6);
}
TF1* fncFdwn = this->GetFdwnFraction(fParticle + "_Frac_Fdwn_Fit_Pt");
//
TString nosec = (sec == "Fdwn") ? "Mat" : "Fdwn";
- for(Int_t i=fLowPtBin; i<fHiPtBin; ++i)
+ Int_t lowbin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMin);
+ Int_t hibin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax);
+
+ for(Int_t i=lowbin; i<hibin; ++i)
{
TH1D* hDCAxy = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
TH1D* hMCDCAxy = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i);
TH1D* hSecDCAxy = hSecDCAxyPt->ProjectionY(Form("%s_%s_DCAxy_%02d",fParticle.Data(),sec.Data(),i),i,i);
// fractions
- Double_t frac[2] = {0,0};
- Double_t err[2] = {0,0};
+ Double_t frac[2] = {0};
+ Double_t err[2] = {0};
- if(fFracProc == kTFractionFitter)
- {
- this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
- }
- else
- {
- this->GetRooFitFractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
- }
+ this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
// fill
TH1D* hFracPt[] = { hPrimPt, hSecPt};
// slice the DCA distribution and get the fractions for each pt bin
// (3 contributions)
//
- for(Int_t i=fLowPtBin; i<fHiPtBin; ++i)
+ Int_t lowbin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMin);
+ Int_t hibin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax);
+
+ for(Int_t i=lowbin; i<hibin; ++i)
{
// slices
TH1D* hDCAxy = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
Double_t frac[3] = {0};
Double_t err[3] = {0};
- if(fFracProc == kTFractionFitter)
- {
- this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hMatDCAxy, hFdwnDCAxy, i);
- }
- else
- {
- this->GetRooFitFractions(frac, err, hDCAxy, hPrimDCAxy, hMatDCAxy, hFdwnDCAxy, i);
- }
+ this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hMatDCAxy, hFdwnDCAxy, i);
// fill
for(Int_t k=0; k<3; ++k)
fit->Constrain(2,0.0001,1.);
Int_t status = fit->Fit();
+ status = fit->Fit();
+ status = fit->Fit();
if (status == 0) // get fractions
{
//fit->SetRangeX(1,50);
Int_t status = fit->Fit();
+ status = fit->Fit();
+ status = fit->Fit();
if (status == 0) // get fractions
{
//
// clone histogram and reset to zero
//
- TH1D* clone = (TH1D*)h->Clone(name.Data());
+ TH1D* clone = dynamic_cast<TH1D*>(h->Clone(name.Data()));
+ if(clone == 0)
+ {
+ this->Error("ZeroClone", "could not clone %s", h->GetName());
+ exit(1);
+ }
clone->Reset();
clone->Sumw2();
return clone;
}
-void AliLnSecondaries::GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hSec, Int_t ibin, const TString& secname) const
-{
-//
-// DCAxy model 2 contributions
-// RooHistPdf class only represents the distribution as a fixed shape and
-// does not propagate the statistical uncertainty to the likelihood
-//
- using namespace RooFit;
-
- RooWorkspace* w = new RooWorkspace("RigidTemplates");
-
- RooRealVar x("x", "sign #times DCA_{xy} (cm)", fMinDCAxy, fMaxDCAxy);
- w->import(x);
-
- RooDataHist dataPrim("data_sig", "primaries", x, hPrim);
- RooHistPdf prim("Sprim", "Template for primaries", x, dataPrim);
- w->import(prim);
-
- RooDataHist dataSec("data_sec", "secondaries", x, hSec);
- RooHistPdf sec("Ssec", "Template for secondaries", x, dataSec);
- w->import(sec);
-
- w->factory(Form("SUM::model( Nprim[0,%f]*Sprim, Nsec[0,%f]*Ssec)",hData->Integral(), 0.8*hData->Integral()));
-
- // data
- RooDataHist data("data", "data to fit", x, hData);
- w->import(data);
-
- w->pdf("model")->fitTo(data, Minos(kTRUE), PrintLevel(-1), Verbose(kFALSE), Warnings(kFALSE), PrintEvalErrors(-1));
-
- // --------- get fractions --------------
-
- Double_t nevents = w->var("Nprim")->getVal()+w->var("Nsec")->getVal();
-
- frac[0] = w->var("Nprim")->getVal()/nevents;
- err[0] = w->var("Nprim")->getError()/nevents;
-
- frac[1] = w->var("Nsec")->getVal()/nevents;
- err[1] = w->var("Nsec")->getError()/nevents;
-
- // ------ debug ------------
-
- TH1D* hDebugFit = this->ZeroClone(hData, fParticle + "_Fit_Data_DCAxy_" + Form("%02d",ibin));
- TH1D* hDebugPrim = this->ZeroClone(hData, fParticle + "_Fit_Prim_DCAxy_" + Form("%02d",ibin));
- TH1D* hDebugSec = this->ZeroClone(hData, fParticle + "_Fit_" + secname + "_DCAxy_" + Form("%02d",ibin));
-
- w->pdf("model")->fillHistogram(hDebugFit,x,nevents);
- w->pdf("Sprim")->fillHistogram(hDebugPrim,x,w->var("Nprim")->getVal());
- w->pdf("Ssec")->fillHistogram(hDebugSec,x,w->var("Nsec")->getVal());
-
- hDebugFit->Write();
- hDebugPrim->Write();
- hDebugSec->Write();
-
- delete hDebugFit;
- delete hDebugPrim;
- delete hDebugSec;
-
- // ---------- end debug --------------
-
- delete w;
-}
-
-void AliLnSecondaries::GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hMat, const TH1D* hFdwn, Int_t ibin) const
-{
-//
-// DCAxy model 3 contributions
-// RooHistPdf class only represents the distribution as a fixed shape and
-// does not propagate the statistical uncertainty to the likelihood
-//
- using namespace RooFit;
-
- RooWorkspace* w = new RooWorkspace("RigidTemplates");
-
- RooRealVar x("x", "sign #times DCA_{xy} (cm)", fMinDCAxy, fMaxDCAxy);
- w->import(x);
-
- RooDataHist dataPrim("data_sig", "data for primaries", x, hPrim);
- RooHistPdf prim("Sprim", "Template for primaries", x, dataPrim);
- w->import(prim);
-
- RooDataHist dataMat("data_mat", "MC for materials", x, hMat);
- RooHistPdf mat("Smat", "Template for matondaries", x, dataMat);
- w->import(mat);
-
- RooDataHist dataFdwn("data_fdwn", "MC for feed-down", x, hFdwn);
- RooHistPdf fdwn("Sfdwn", "Template for feed-down", x, dataFdwn);
- w->import(fdwn);
-
- w->factory(Form("SUM::model( Nprim[0.,%f]*Sprim, Nmat[0.,%f]*Smat, Nfdwn[0.,%f]*Sfdwn )",hData->Integral(),0.8*hData->Integral(),0.8*hData->Integral()));
-
- // --------------- fit -----------------
- RooDataHist data("data", "data to fit", x, hData);
- w->import(data);
-
- w->pdf("model")->fitTo(data, Minos(kTRUE), PrintLevel(-1), Verbose(kFALSE), Warnings(kFALSE), PrintEvalErrors(-1));
-
- // --------- get fractions --------------
-
- Double_t nevents = w->var("Nprim")->getVal()+w->var("Nmat")->getVal()+w->var("Nfdwn")->getVal();
-
- frac[0] = w->var("Nprim")->getVal()/nevents;
- err[0] = w->var("Nprim")->getError()/nevents;
-
- frac[1] = w->var("Nmat")->getVal()/nevents;
- err[1] = w->var("Nmat")->getError()/nevents;
-
- frac[2] = w->var("Nfdwn")->getVal()/nevents;
- err[2] = w->var("Nfdwn")->getError()/nevents;
-
- // ------ debug ------------
-
- TH1D* hDebugFit = this->ZeroClone(hData, fParticle + "_Fit_Data_DCAxy_" + Form("%02d",ibin));
- TH1D* hDebugPrim = this->ZeroClone(hData, fParticle + "_Fit_Prim_DCAxy_" + Form("%02d",ibin));
- TH1D* hDebugMat = this->ZeroClone(hData, fParticle + "_Fit_Mat_DCAxy_" + Form("%02d",ibin));
- TH1D* hDebugFdwn = this->ZeroClone(hData, fParticle + "_Fit_Fdwn_DCAxy_" + Form("%02d",ibin));
-
- w->pdf("model")->fillHistogram(hDebugFit,x,hData->Integral());
- w->pdf("Sprim")->fillHistogram(hDebugPrim,x,w->var("Nprim")->getVal());
- w->pdf("Smat")->fillHistogram(hDebugMat,x,w->var("Nmat")->getVal());
- w->pdf("Sfdwn")->fillHistogram(hDebugFdwn,x,w->var("Nfdwn")->getVal());
-
- hDebugFit->Write();
- hDebugPrim->Write();
- hDebugMat->Write();
- hDebugFdwn->Write();
-
- delete hDebugFit;
- delete hDebugPrim;
- delete hDebugMat;
- delete hDebugFdwn;
-
- // ----------- end debug --------------
-
- delete w;
-}
-
TH2D* AliLnSecondaries::GetFlatDCAxyPt(Double_t max, const TH2D* hDCAxyPt, const TString& name) const
{
//
void SetOutputTag(const TString& tag) { fOutputTag = tag; }
- void SetCorBins(Int_t lowbin, Int_t hibin) { fLowPtBin = lowbin; fHiPtBin = hibin; }
+ void SetCorBins(Double_t min, Double_t max) { fPtMin = min; fPtMax = max; }
void SetDCAxyInterval(Double_t lowdca, Double_t hidca) { fMinDCAxy = lowdca; fMaxDCAxy = hidca; }
void SetNBin(Int_t nbin) { fNbin = nbin; }
void SetAddFakeTracks(Bool_t flag=1) { fAddFakeTracks = flag; }
- enum { kTFractionFitter=0, kRooFit, kMonteCarlo };
- enum { kGeantDCAxy=0, kFlatDCAxy};
+ enum { kTFractionFitter=0, kMonteCarlo };
+ enum { kGeantDCAxy=0, kFlatDCAxy };
private:
Int_t GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hMat, TH1D* hFdwn, Int_t ibin) const;
Int_t GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hSec, Int_t ibin, const TString& secName) const;
- void GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hSec, Int_t ibin, const TString& secName) const;
- void GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hMat, const TH1D* hFdwn, Int_t ibin) const;
-
TH2D* GetFlatDCAxyPt(Double_t norm, const TH2D* hDCAxyPt, const TString& name) const;
TF1* GetMatFraction(const TString& name) const;
TString fOutputFilename; // output filename
TString fOutputTag; // tag for the ouput
- Int_t fLowPtBin; // low pt bin for the corrections
- Int_t fHiPtBin ; // high pt bin for the corrections
+ Double_t fPtMin; // minimum pt value for the corrections
+ Double_t fPtMax ; // maximum pt value for the corrections
Int_t fNbin; // for rebinning DCA distributions
* provided "as is" without express or implied warranty. *
**************************************************************************/
-// invariant differential yields and cross sections
+// (invariant) differential yield
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
#include <Riostream.h>
ClassImp(AliLnSpectra)
-AliLnSpectra::AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag, const Double_t xsec[3])
+AliLnSpectra::AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag)
: TObject()
, fParticle(particle)
, fPtFilename(ptFilename)
, fYMin(-0.5)
, fYMax(0.5)
, fINEL(1)
-, fIsOnlyGen(0)
-, fSysErr(1)
{
//
// constructor
//
- for(Int_t i=0; i<3; ++i) fInelXsec[i] = xsec[i];
}
AliLnSpectra::~AliLnSpectra()
Int_t AliLnSpectra::Exec()
{
//
-// (invariant) differential yield projection in pt for |y| < 0.5
+// (invariant) differential yield projection in pt for |y| < ymax-ymin
//
using namespace std;
// pt and number of events
- TH1D* hPt = (TH1D*)FindObj(finput, fParticle + "_Pt");
- TH1D* hStats = (TH1D*)FindObj(finput, fParticle + "_Stats");
+ TH1D* hPt = FindObj<TH1D>(finput, fParticle + "_Pt");
+ TH1D* hStats = FindObj<TH1D>(finput, fParticle + "_Stats");
Double_t nEvent = hStats->Integral(3,3);
if(fINEL) nEvent = hStats->Integral(4,4);
TGraphErrors* grDYieldPt = this->GetDiffYieldPt(hPt, nEvent, fParticle + "_DiffYield_Pt");
grDYieldPt->Write();
- TGraphErrors* grSysErrDYieldPt = this->AddSystError(grDYieldPt, fSysErr, fParticle + "_SysErr_DiffYield_Pt");
- grSysErrDYieldPt->Write();
-
// invariant differential yield
TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, fParticle + "_InvDiffYield_Pt");
grInvDYieldPt->Write();
- TGraphErrors* grSysErrInvDYieldPt = this->AddSystError(grInvDYieldPt, fSysErr, fParticle + "_SysErr_InvDiffYield_Pt");
- grSysErrInvDYieldPt->Write();
-
- // invariant differential cross section
-
- TGraphErrors* grInvXsectPt = GetInvDiffXsectionPt(grInvDYieldPt, fInelXsec, fParticle + "_InvDiffXSection_Pt");
- grInvXsectPt->Write();
-
- TGraphErrors* grSysErrInvXsectPt = GetInvDiffXsectionPt(grSysErrInvDYieldPt, fInelXsec, fParticle + "_SysErr_InvDiffXSection_Pt");
- grSysErrInvXsectPt->Write();
-
// clean
delete grDYieldPt;
delete grInvDYieldPt;
- delete grSysErrDYieldPt;
- delete grSysErrInvDYieldPt;
- delete grInvXsectPt;
delete foutput;
delete finput;
Double_t yield = n/(nEvent*dpt*dy);
Double_t yError = yield*eN/n;
- grDYieldPt->SetPoint(j,pt,yield);
- grDYieldPt->SetPointError(j++,0,yError);
+ grDYieldPt->SetPoint(j, pt, yield);
+ grDYieldPt->SetPointError(j++, dpt/2., yError);
}
return grDYieldPt;
delete grDYieldPt;
return grInvDYieldPt;
}
-
-TGraphErrors* AliLnSpectra::GetInvDiffXsectionPt(const TGraphErrors* grInvDYieldPt, const Double_t* sigma, const TString& name) const
-{
-//
-// invariant differential cross section, only stat. errors
-// (multiply each value for the total cross section)
-//
- Double_t xsec = sigma[0];
- Double_t errXSec = sigma[1];
-
- TGraphErrors* grInvDXsecPt = new TGraphErrors();
- grInvDXsecPt->SetName(name.Data());
-
- for(Int_t i=0; i < grInvDYieldPt->GetN(); ++i)
- {
- Double_t pt, y;
- grInvDYieldPt->GetPoint(i,pt,y);
-
- Double_t errpt = grInvDYieldPt->GetErrorX(i);
- Double_t erry = grInvDYieldPt->GetErrorY(i);
-
- Double_t invxsec = xsec*y;
- Double_t err = invxsec*TMath::Sqrt(TMath::Power(errXSec/xsec,2) + TMath::Power(erry/y,2));
-
- grInvDXsecPt->SetPoint(i, pt, invxsec);
- grInvDXsecPt->SetPointError(i, errpt, err);
- }
-
- return grInvDXsecPt;
-}
-
-TGraphErrors* AliLnSpectra::AddSystError(const TGraphErrors* gr, Double_t percent, const TString& name) const
-{
-//
-// set the error of h as the given percent of its value
-//
- TGraphErrors* grSyst = new TGraphErrors();
- grSyst->SetName(name.Data());
-
- for(Int_t i=0; i < gr->GetN(); ++i)
- {
- Double_t x, y;
- gr->GetPoint(i,x,y);
- Double_t err = percent*y;
-
- grSyst->SetPoint(i, x, y);
- grSyst->SetPointError(i, 0.025, err);
- }
-
- return grSyst;
-}
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-// invariant differential yields and cross sections
+// (invariant) differential yield
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
#include <TObject.h>
class TGraphErrors;
+class TString;
class TF1;
class AliLnSpectra: public TObject
{
public:
- AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag, const Double_t xsec[3]);
+ AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag);
virtual ~AliLnSpectra();
TGraphErrors* GetDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const;
TGraphErrors* GetInvDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const;
TGraphErrors* GetInvDiffYieldPt(const TGraphErrors* grDYieldPt, const TString& name) const;
- TGraphErrors* GetInvDiffXsectionPt(const TGraphErrors* grInvDYieldPt, const Double_t* sigma, const TString& name) const;
-
- TGraphErrors* AddSystError(const TGraphErrors* gr, Double_t percent, const TString& name) const;
Int_t Exec();
void SetRapidityInterval(Double_t ymin, Double_t ymax) { fYMin = ymin; fYMax = ymax; }
void SetExtrapolateToINEL(Bool_t flag=1) { fINEL = flag; }
- void SetOnlyGeneration(Bool_t flag=1) { fIsOnlyGen = flag; }
-
- void SetScalingFactor(Double_t syserr=1) { fSysErr = syserr; }
-
- void SetInelXSection(Double_t xsec, Double_t statErr, Double_t systErr) { fInelXsec[0] = xsec; fInelXsec[1] = statErr; fInelXsec[2] = systErr; }
-
private:
AliLnSpectra(const AliLnSpectra& other);
Double_t fYMax; // rapidity high limit
Bool_t fINEL; // extrapolate to inelastic events
- Bool_t fIsOnlyGen; // if no need for correction
-
- Double_t fSysErr; // variation for systematics
- Double_t fInelXsec[3]; // total inelastic cross section, syst. and stat. errors (mb)
ClassDef(AliLnSpectra,1)
};
+++ /dev/null
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes is hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-
-// unfolding correction for the pt
-// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-
-//#include <Riostream.h>
-#include <TFile.h>
-#include <TSystem.h>
-#include <TH1D.h>
-#include <TH2D.h>
-#include <TString.h>
-#include <TF1.h>
-
-#include "AliLnUnfolding.h"
-#include "B2.h"
-
-ClassImp(AliLnUnfolding)
-
-AliLnUnfolding::AliLnUnfolding(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag)
-: TObject()
-, fParticle(particle)
-, fSimuFilename(simuFilename)
-, fOutputFilename(outputFilename)
-, fOutputTag(otag)
-{
-//
-// constructor
-//
- TH1::SetDefaultSumw2();
-}
-
-AliLnUnfolding::~AliLnUnfolding()
-{
-//
-// destructor
-//
-}
-
-Int_t AliLnUnfolding::Exec()
-{
-//
-// extract the detector response from the simulation for particle/antiparticle
-//
- using namespace std;
-
- TFile* fsimu = new TFile(fSimuFilename.Data(), "read");
- if(fsimu->IsZombie()) exit(1);
-
- TFile* foutput = new TFile(fOutputFilename.Data(), "recreate");
- if(fOutputTag != "")
- {
- foutput->mkdir(fOutputTag.Data());
- foutput->cd(fOutputTag.Data());
- }
-
- TH2D* hResponseMtx = (TH2D*)FindObj(fsimu, fParticle + "_Prim_Response_Matrix");
- hResponseMtx->SetName(Form("%s_Response_Matrix",fParticle.Data()));
- hResponseMtx->Write();
-
- TH1D* hMeasuredPt = (TH1D*)FindObj(fsimu, fParticle + "_PID_Pt");
- hMeasuredPt->SetName(Form("%s_Measured_Pt", fParticle.Data()));
- hMeasuredPt->Reset();
- hMeasuredPt->Write();
-
- TH1D* hTruePt = (TH1D*)FindObj(fsimu, fParticle + "_Sim_PID_Prim_Pt");
- hTruePt->SetName(Form("%s_True_Pt", fParticle.Data()));
- hTruePt->Reset();
- hTruePt->Write();
-
- delete foutput;
- delete fsimu;
-
- return 0;
-}
+++ /dev/null
-#ifndef ALILNUNFOLDING_H
-#define ALILNUNFOLDING_H
-
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice */
-
-// unfolding correction for the pt
-// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-
-#include <TObject.h>
-#include <TString.h>
-
-class TString;
-
-class AliLnUnfolding: public TObject
-{
- public:
-
- AliLnUnfolding(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag);
- virtual ~AliLnUnfolding();
-
- const TString* GetOutputFilename() const { return &fOutputFilename; }
-
- Int_t Exec();
-
- void SetParticle(const TString& particle) { fParticle = particle; }
-
- void SetOutputTag(const TString& tag) { fOutputTag = tag; }
-
- private:
-
- AliLnUnfolding(const AliLnUnfolding& other);
- AliLnUnfolding& operator=(const AliLnUnfolding& other);
-
- private:
-
- TString fParticle; // particle
-
- TString fSimuFilename; // simulation filename
- TString fOutputFilename; // output filename
- TString fOutputTag; // tag for the ouput
-
- ClassDef(AliLnUnfolding,1)
-};
-
-#endif // ALILNUNFOLDING_H
// some common functions
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <TObject.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TFile.h>
#include <TList.h>
#include <TString.h>
#include <TH1D.h>
#include <TF1.h>
#include <cstdlib>
+#endif
-TObject* FindObj(TFile* f, const TString& name)
+template <class T>
+inline T* FindObj(TFile* f, const TString& name)
{
//
// check if the object exists
//
- TObject* obj = f->Get(name.Data());
+ T* obj = dynamic_cast<T*>(f->Get(name.Data()));
if(obj == 0)
{
- f->Error("Get","%s not found",name.Data());
+ f->Error("Get","%s not found", name.Data());
exit(1);
}
return obj;
}
-TObject* FindObj(TFile* f, const TString& dir, const TString& name)
+template <class T>
+inline T* FindObj(TFile* f, const TString& dir, const TString& name)
{
//
// check if the object exists
//
- if(dir=="") return FindObj(f,name);
-
- TObject* obj;
- f->GetObject(Form("%s/%s;1",dir.Data(),name.Data()),obj);
- if(obj == 0)
- {
- f->Error("GetObject","%s/%s not found",dir.Data(),name.Data());
- exit(1);
- }
-
- return obj;
+ if(dir=="") return FindObj<T>(f,name);
+ return FindObj<T>(f, dir + "/" + name + ";1");
}
-TObject* FindObj(const TList* l, const TString& name)
+template <class T>
+inline T* FindObj(const TList* l, const TString& name)
{
//
// check if the object exists
//
- TObject* obj = l->FindObject(name.Data());
+ T* obj = dynamic_cast<T*>(l->FindObject(name.Data()));
+
if(obj == 0)
{
- l->Error("FindObject","%s not found",name.Data());
+ l->Error("FindObject","%s not found", name.Data());
exit(1);
}
return obj;
}
-TH1D* Divide(const TH1* hX, const TH1* hY, const TString& name)
+inline TH1D* Divide(const TH1* hX, const TH1* hY, const TString& name)
{
//
// clone and divide
//
- TH1D* q = (TH1D*)hX->Clone(name.Data());
+ TH1D* q = dynamic_cast<TH1D*>(hX->Clone(name.Data()));
+
+ if(q == 0)
+ {
+ hX->Warning("Clone", "could not clone %s", hX->GetName());
+ return q;
+ }
+
q->Sumw2();
q->Divide(hY);
+
return q;
}
-Double_t GetMass(const TString& name)
+inline Double_t GetMass(const TString& name)
{
//
// return particle mass
return 0;
}
-TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+inline TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
{
//
// Tsallis distribution
fnc->SetParNames("gV","q","T");
fnc->SetParameters(100., 1.1, 0.07);
- fnc->SetParLimits(0, 1., 1.e+7);
+ fnc->SetParLimits(0, 0., 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)
+inline 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(0, 0., 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)
+inline TF1* PtTsallisDYield(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
{
//
-// Tsallis-Pareto distribution
+// Tsallis distribution for mean pt
+//
+ TF1* fnc = new TF1(name.Data(), Form("[0]*x*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, 0., 1.e+7);
+ fnc->SetParLimits(1, 1.0001, 3.);
+ fnc->SetParLimits(2, 0.001, 0.3);
+
+ return fnc;
+}
+
+inline TF1* QTsallis(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// q-Tsallis distribution
// Phys. Rev. C 83, 064903 (2011)
// Phys. Rev. C 75, 064901 (2007)
//
fnc->SetParNames("dN/dy","n","C");
fnc->SetParameters(0.05, 7, 0.3);
- fnc->SetParLimits(0, 0, 1);
+ fnc->SetParLimits(0, 0, 10);
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)
+inline TF1* QTsallisDYield(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)
+// q-Tsallis distribution for differential yield
//
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);
return fnc;
}
+inline TF1* PtQTsallisDYield(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10)
+{
+//
+// q-Tsallis distribution for mean pt
+//
+ TF1* fnc = new TF1(name.Data(), Form("x*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
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* See cxx source for full Copyright notice */
-// functions for configs
+// common functions for config macros
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TROOT.h>
#include <TString.h>
-#include <TStyle.h>
#include <cstdlib>
+#endif
namespace CollSystem
{
//
-// pp collision systems
+// pp collision energies
//
const Int_t kNener = 3;
- const TString kEnergyTag[kNener] = { "900GeV", "2.76TeV", "7TeV" };
- const Double_t kEnergy[kNener] = { 0.9, 2.76, 7. };
+ const TString kEnergyTag[kNener] = { "-900GeV", "-2.76TeV", "-7TeV" };
+ const TString kEnergyLabel[kNener] = { "900 GeV", "2.76 TeV", "7 TeV" };
+ const Double_t kEnergy[kNener] = { 0.9, 2.76, 7 };
const Double_t kEnergyError[kNener] = { 0 };
const TString kEnergyName = "#sqrt{s} (TeV)";
};
//
// multiplicity classes
//
- 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 Int_t kNmult = 5;
+ const TString kMultTag[kNmult] = { "-ntrk0103", "-ntrk0405", "-ntrk0608", "-ntrk0914", "-ntrk15xx" };
+ const Double_t kKNOmult[kNmult] = { 0.35, 0.74, 1.13, 1.82, 3.14 };
+ const Double_t kDtrPtMax[kNmult] = { 1.8, 2.0, 2.2, 2.3, 2.5 };
const Double_t kKNOmultErr[kNmult] = { 0 };
const TString kKNOmultName = "z";
};
-namespace HiMult
+namespace B2HiLowMult
{
//
-// high multiplicity
+// high and low multiplicity classes
//
- const Int_t kNmult = 6;
- const TString kMultTag[kNmult] = { "kno04xx", "kno05xx", "kno06xx", "kno07xx", "kno08xx", "kno09xx"};
- const Double_t kKNOmult[kNmult] = { 4.16, 5.85, 6.70, 7.70, 8.73, 9.75 };
+ const Int_t kNmult = 2;
+ const TString kMultTag[kNmult] = { "-ntrk0108", "-ntrk09xx" };
+ const Double_t kKNOmult[kNmult] = { 0.64, 2.23 };
+ const Double_t kDtrPtMax[kNmult] = { 2.0, 2.5 };
const Double_t kKNOmultErr[kNmult] = { 0 };
const TString kKNOmultName = "z";
};
return "unknown";
}
-void GetInelXSection(Double_t xsection[3], const TString& period)
-{
-//
-// inelastic cross section in mb and largest stat. and syst error
-// measured by ALICE for the given colliding system
-// http://arxiv.org/abs/1208.4968
-//
- TString collsystem = GetCollSystem(period);
-
- if( collsystem == "pp0.9TeV" )
- {
- xsection[0] = 50.3;
- xsection[1] = 0.4;
- xsection[2] = 1.0;
- }
- else if( collsystem == "pp2.76TeV" )
- {
- xsection[0] = 62.8;
- xsection[1] = 4.0;
- xsection[2] = 4.6;
- }
- else if( collsystem == "pp7TeV" )
- {
- xsection[0] = 73.2;
- xsection[1] = 1.2;
- xsection[2] = 2.6;
- }
- else
- {
- std::cerr << "Warning: unknown colliding system " << collsystem << " for period " << period << std::endl;
- }
-}
-
void GetTriggerEfficiency(Double_t trigeff[3], const TString& trigname, const TString& period)
{
//
-// trigger efficiency for the given colliding system and largest stat. and syst. errors
+// trigger efficiency for the given colliding system and largest systematic error
// http://arxiv.org/abs/1208.4968
//
TString collsystem = GetCollSystem(period);
if(period=="lhc10c") return "lhc12a5bc";
if(period=="lhc10d") return "lhc12a5bd";
if(period=="lhc10e") return "lhc12a5be";
+ if(period=="lhc10bc") return "lhc12a5bbc";
if(period=="lhc10de") return "lhc12a5bde";
if(period=="lhc10cde") return "lhc12a5bcde";
+ if(period=="lhc10bcd") return "lhc12a5bbcd";
if(period=="lhc10bcde") return "lhc12a5bbcde";
if(period=="lhc11a_wsdd") return "lhc12a5c_wsdd";
//
// make input name for data
//
+ if(trksel == "") return species + "-" + period;
+ if(trksel.BeginsWith("-")) return species + "-" + period + trksel;
return species + "-" + period + "-" + trksel;
}
return species + "-" + outputTag;
}
-TStyle* GetDrawingStyle()
-{
-//
-// define a default style for drawing
-//
- TStyle* st = new TStyle();
-
- //st->SetPalette(51);
- st->SetPalette(1);
- st->SetPadTickX(1);
- st->SetPadTickY(1);
- st->SetPadGridX(1);
- st->SetPadGridY(1);
-
- st->SetCanvasColor(0);
- st->SetFrameBorderMode(0);
- st->SetStatBorderSize(1);
- st->SetStatColor(0);
- st->SetFrameFillColor(0);
- st->SetTitleFillColor(0);
- st->SetLabelFont(62,"XYZ");
- st->SetTitleFont(62,"XYZ");
-
- //st->SetOptTitle(1);
- st->SetOptStat(0);
-
- return st;
-}
-
void DrawOutputCorr(const TString& species, const TString& corr, const TString& tag="")
{
//
gROOT->ProcessLine(Form(".x DrawCorr.C+g(\"%s\",\"%s\",\"%s\")", species.Data(), corr.Data(), tag.Data()));
}
-void DrawCorrDebug(const TString& sec, const TString& tag, const TString& species, Int_t lowbin=1, Int_t hibin=10, Double_t dcaxyMin=-1.5, Double_t dcaxyMax=1.5)
+void DrawCorrDebug(const TString& sec, const TString& tag, const TString& species, Double_t ptmin=0.5, Double_t ptmax=3., Double_t dcaxyMin=-1.5, Double_t dcaxyMax=1.5)
{
//
// draw correction for secondaries
for(Int_t i=0; i<2; ++i)
{
- gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"%s\",\"%s\", %d, %d, %f, %f)", sec.Data(), tag.Data(), kParticle[i].Data(), lowbin, hibin, dcaxyMin, dcaxyMax));
+ gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"%s\",\"%s\", %f, %f, %f, %f)", sec.Data(), tag.Data(), kParticle[i].Data(), ptmin, ptmax, dcaxyMin, dcaxyMax));
}
}
-void DrawPtDebug(const TString& pt, const TString& tag, const TString& species, Bool_t m2pid=0, Int_t hiptbin=17, Int_t lowm2bin=9, Int_t him2bin=17)
+void DrawPtDebug(const TString& pt, const TString& tag, const TString& species, Bool_t m2pid=0, Double_t ptmax=3., Double_t ptpid=1.2)
{
//
// draw pt debug for the particle species
for(Int_t i=0; i<2; ++i)
{
- gROOT->ProcessLine(Form(".x DrawPt.C+g(\"%s\",\"%s\",\"%s\",%d, %d, %d, %d)", pt.Data(), tag.Data(), kParticle[i].Data(), hiptbin, m2pid, lowm2bin, him2bin));
+ gROOT->ProcessLine(Form(".x DrawPt.C+g(\"%s\",\"%s\",\"%s\",%f, %d, %f)", pt.Data(), tag.Data(), kParticle[i].Data(), ptmax, m2pid, ptpid));
}
}
// coalescence parameter
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TSystem.h>
#include <TROOT.h>
#include <TFileMerger.h>
#include <TString.h>
-
#include "AliLnB2.h"
+#endif
Int_t B2( const TString& pSpectra = "~/alice/output/Proton-lhc10d-Spectra.root"
, const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Spectra.root"
// B2 as a function of multiplicity
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TSystem.h>
#include <TROOT.h>
#include <TString.h>
#include <TFile.h>
#include <TGraphErrors.h>
-
#include "AliLnB2.h"
+#endif
+
#include "B2.h"
#include "Config.h"
return 0.046133 + 0.0484458*z;
}
-Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-Mult-Spectra.root"
- , const TString& ptag = "lhc10d-nsd"
- , const TString& dSpectra = "~/alice/output/Deuteron-lhc10bcde-nsd-Mult-Spectra.root"
- , const TString& dtag = "lhc10bcde-nsd"
- , const TString& outputMultPt = "~/alice/output/B2-Mult-Pt.root"
- , const TString& outputPtMult = "~/alice/output/B2-Pt-Mult.root"
- , const TString& otag = "pp-nsd")
+extern void RatioMult( const TString&, const TString&, const TString&, const TString&, const Int_t, const TString*, const Double_t*, const Double_t*, const TString&, const Bool_t, const TString& , const TString& );
+
+Int_t B2Mult( 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"
+ , const TString& outputMultPt = "~/alice/output/B2-lhc10d-MultPt.root"
+ , const TString& outputPtMult = "~/alice/output/B2-lhc10d-PtMult.root"
+ , const TString& otag = "lhc10d"
+ , const Bool_t qTsallis = 0
+ , const TString& tsallisTag = "Tsallis"
+ , const TString& oratio = "~/alice/output/Particle-Ratios-lhc10d-Tsallis.root")
{
//
// B2 as a function of multiplicity
//
using namespace B2mult;
+ using namespace std;
const Int_t kNpart = 2;
{
TString b2file = kPrefix[j] + "B2.root";
- AliLnB2 b2(pSpectra, ptag + "-" + kMultTag[i], dSpectra, dtag + "-" + kMultTag[i], b2file, otag + "-" + kMultTag[i], 2, kCharge[j]);
+ cout << kMultTag[i] << endl;
+
+ AliLnB2 b2(pSpectra, ptag + kMultTag[i], dSpectra, dtag + kMultTag[i], b2file, otag + kMultTag[i], 2, kCharge[j]);
b2.SetCd(GetCd(kKNOmult[i]));
// merge B2 and B2bar
- TString outputfile = otag + "-" + kMultTag[i] + "-B2.root";
+ TString outputfile = otag + kMultTag[i] + "-B2.root";
m.OutputFile(outputfile.Data());
m.Merge();
for(Int_t i=0; i<kNmult; ++i)
{
- TString b2 = otag + "-" + kMultTag[i] + "-B2.root";
+ TString b2 = otag + kMultTag[i] + "-B2.root";
m.AddFile(b2.Data(),0);
}
for(Int_t i=0; i<kNmult; ++i)
{
- gSystem->Exec(Form("rm -f %s-%s-B2.root",otag.Data(),kMultTag[i].Data()));
+ gSystem->Exec(Form("rm -f %s%s-B2.root",otag.Data(),kMultTag[i].Data()));
}
// B2 as a function of multiplicity for each pt
TGraphErrors* grB2pt[kNpart][kNmult];
TGraphErrors* grR3pt[kNpart][kNmult];
+ TGraphErrors* grSysB2pt[kNpart][kNmult];
+ TGraphErrors* grSysR3pt[kNpart][kNmult];
for(Int_t i=0; i<kNpart; ++i)
{
for(Int_t j=0; j<kNmult; ++j)
{
- grB2pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultTag[j], Form("B2%s_Pt", kSuffix[i].Data()));
- grR3pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultTag[j], Form("R3%s_Pt", kSuffix[i].Data()));
+ grB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("B2%s_Pt", kSuffix[i].Data()));
+ grR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("R3%s_Pt", kSuffix[i].Data()));
+
+ grSysB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_B2%s_Pt", kSuffix[i].Data()));
+ grSysR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_R3%s_Pt", kSuffix[i].Data()));
}
}
for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
{
- ptLabel[i] = Form("pT%.02fA",pt[i]);
+ ptLabel[i] = Form("pT %.02fA",pt[i]);
foutput->mkdir(ptLabel[i].Data());
}
{
Double_t B2[kNmult];
Double_t B2StatErr[kNmult];
+ Double_t B2SystErr[kNmult];
Double_t R3[kNmult];
Double_t R3StatErr[kNmult];
+ Double_t R3SystErr[kNmult];
for(Int_t k=0; k<kNmult; ++k)
{
- Double_t x, y, ey;
+ Double_t x, y, staterr, systerr;
grB2pt[i][k]->GetPoint(j,x,y);
- ey = grB2pt[i][k]->GetErrorY(j);
+ staterr = grB2pt[i][k]->GetErrorY(j);
+ systerr = grSysB2pt[i][k]->GetErrorY(j);
B2[k] = y;
- B2StatErr[k] = ey;
+ B2StatErr[k] = staterr;
+ B2SystErr[k] = systerr;
grR3pt[i][k]->GetPoint(j,x,y);
- ey = grR3pt[i][k]->GetErrorY(j);
+ staterr = grR3pt[i][k]->GetErrorY(j);
+ systerr = grSysR3pt[i][k]->GetErrorY(j);
R3[k] = y;
- R3StatErr[k] = ey;
+ R3StatErr[k] = staterr;
+ R3SystErr[k] = systerr;
}
TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
grR3Mult->SetName(Form("R3%s_Zmult", kSuffix[i].Data()));
+ Double_t zMultSystErr[kNmult];
+ for(Int_t k=0; k<kNmult; ++k) zMultSystErr[k]=0.07;
+
+ TGraphErrors* grSysB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, zMultSystErr, B2SystErr);
+ grSysB2Mult->SetName(Form("SystErr_B2%s_Zmult", kSuffix[i].Data()));
+
+ TGraphErrors* grSysR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, zMultSystErr, R3SystErr);
+ grSysR3Mult->SetName(Form("SystErr_R3%s_Zmult", kSuffix[i].Data()));
+
foutput->cd(ptLabel[j].Data());
grB2Mult->Write();
grR3Mult->Write();
+ grSysB2Mult->Write();
+ grSysR3Mult->Write();
delete grB2Mult;
delete grR3Mult;
+ delete grSysB2Mult;
+ delete grSysR3Mult;
}
}
delete foutput;
delete finput;
- // particle ratios
+ //
+ // Particle ratios
+ // -----------------------------------
+
+ RatioMult(pSpectra, dSpectra, ptag, dtag, kNmult, kMultTag, kKNOmult, kKNOmultErr, kKNOmultName, qTsallis, oratio, tsallisTag);
- gROOT->ProcessLine(Form(".x RatioMult.C+g(\"%s\",\"%s\",\"%s\",\"%s\")", pSpectra.Data(), dSpectra.Data(), ptag.Data(), dtag.Data()));
// draw B2 as a function of pt
// LHC10x config for deuterons and antideuterons\r
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>\r
\r
+#if !defined(__CINT__) || defined(__MAKECINT__)\r
#include <TROOT.h>\r
#include <TString.h>\r
#include "AliLnDriver.h"\r
+#endif\r
+\r
#include "Config.h"\r
\r
-Int_t Config_Deuteron_TOF_LHC10x(const TString& inputDir = "~/alice/input",\r
- const TString& outputDir = "~/alice/output",\r
- const TString& period = "lhc10d",\r
- const TString& outputTag = "lhc10d",\r
- const TString& multTag = "",\r
- const TString& multCorTag = "",\r
- Bool_t inel = 1, // for mult\r
- Bool_t drawOutput = 1, // for batch\r
- Int_t lowPtBin = 5, // 0.8 Gev/c\r
- Int_t hiPtBin = 17, // 3.2 GeV/c\r
- Bool_t makeStats = 1,\r
- Bool_t makeCor = 1,\r
- Bool_t makePt = 1,\r
- Bool_t makeRatio = 1,\r
- Bool_t makeSpectra = 1 )\r
+Int_t Config_Deuteron_TOF_LHC10x( const TString& inputDir = "~/alice/input"\r
+ , const TString& outputDir = "~/alice/output"\r
+ , const TString& period = "lhc10bcde"\r
+ , const TString& outputTag = "lhc10bcde"\r
+ , const TString& trkselTag = "-tpc3-nsd-moc"\r
+ , const TString& multTag = ""\r
+ , const TString& multCorTag = ""\r
+ , Double_t ymax = 0.5\r
+ , Bool_t inel = 0 // for mult\r
+ , Bool_t drawOutput = 1 // for batch\r
+ , Double_t ptmin = 0.7\r
+ , Double_t ptmax = 3.0\r
+ , Double_t ptpid = 0.8\r
+ , Bool_t makeStats = 1\r
+ , Bool_t makeCor = 1\r
+ , Bool_t makePt = 1\r
+ , Bool_t makeRatio = 1\r
+ , Bool_t makeSpectra = 1)\r
{\r
//\r
// lhc10b, lhc10c, lhc10d, lhc10e config for deuterons and antideuterons\r
//\r
const TString kSpecies = "Deuteron";\r
- const TString kTrkSel = "its_tpc_tof_dca-tpc3";\r
- const TString kTrigName = "mbor";\r
+ const TString kTrkSel = "its_tpc_tof_dca";\r
+ const TString kTrigName = "mband";\r
const Bool_t kMCtoINEL = 1;\r
- const Int_t kM2Bin[2] = {8,18};\r
- const Bool_t kPidM2 = 1;\r
- const Bool_t kUnfolding = 0;\r
- const Int_t kIter = 5;\r
+ const Bool_t kPid = 1;\r
+ const Int_t kPidProc = 0; // 0 m2, 1 dm2, 2 time\r
+ const Double_t kPidEff = 1.;\r
const Bool_t kSecondaries = 1;\r
- const Int_t kSecProd = 0; // 0 tff, 1 roofit, 2 mc\r
+ const Int_t kSecProc = 0; // 0 tff, 1 mc\r
const Int_t kMatDCAxyMod = 1; // 0 geant, 1 flat\r
const Bool_t kAntiNucTemplate = 0;\r
const Int_t kNbin = 5;\r
const Double_t kDCAxy[2] = {-0.2,0.2};\r
- const Double_t kM2Bkg[2] = {2.2,5.};\r
- const Double_t kM2tpc[2] = {2.,6.5};\r
const Bool_t kEfficiency = 1;\r
- const Bool_t kFitFrac = 0;\r
- const Double_t kSysErr[2] = {0.10, 0.11} ;\r
+ const Bool_t kFitFrac = 1;\r
+ const Int_t kDebugLevel = 1;\r
+ \r
+ Double_t bkgLimitToF[2] = {-2, 2};\r
+ Double_t pidLimitToF[2] = {-2.,6.};\r
\r
- Double_t xsec[3];\r
- GetInelXSection(xsec, period);\r
+ Double_t bkgLimitM2[2] = {1.8,6.0};\r
+ Double_t pidLimitM2[2] = {1.8,6.};\r
+ \r
+ if(kPidProc==1)\r
+ {\r
+ for(Int_t i=0; i<2; ++i)\r
+ {\r
+ bkgLimitM2[i] -= 3.51792;\r
+ pidLimitM2[i] -= 3.51792;\r
+ }\r
+ }\r
\r
Double_t trigEff[3];\r
GetTriggerEfficiency(trigEff, kTrigName, period);\r
\r
// input and output filenames\r
\r
- TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + ".root";\r
- TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
- TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
- TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + "-corr.root";\r
+ TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + ".root";\r
+ TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";\r
+ TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";\r
+ TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + "-corr.root";\r
\r
TString outputPt = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";\r
TString outputRatio = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";\r
driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);\r
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);\r
\r
+ driver.SetRapidityInterval(-ymax,ymax);\r
+ \r
driver.SetOutputTag(outputTag);\r
+ driver.SetOutputCorTag(outputTag);\r
+ \r
driver.SetTriggerEfficiency(trigEff);\r
- driver.SetInelXSection(xsec);\r
driver.SetExtrapolateToINEL(inel);\r
driver.SetMCtoINEL(kMCtoINEL);\r
- driver.SetPtBinInterval(lowPtBin, hiPtBin);\r
- driver.SetPidM2(kPidM2);\r
- driver.SetM2BinInterval(kM2Bin[0], kM2Bin[1]);\r
- driver.SetM2BkgInterval(kM2Bkg[0], kM2Bkg[1]);\r
- driver.SetM2TPCInterval(kM2tpc[0], kM2tpc[1]);\r
- driver.SetUnfolding(kUnfolding, kIter);\r
+ driver.SetPtInterval(ptmin, ptmax);\r
+ driver.SetPid(kPid);\r
+ driver.SetPidProcedure(kPidProc);\r
+ driver.SetPidEfficiency(kPidEff);\r
+ driver.SetPidPt(ptpid);\r
+ \r
+ if(kPidProc==2)\r
+ {\r
+ driver.SetBkgInterval(bkgLimitToF[0], bkgLimitToF[1]);\r
+ driver.SetPidInterval(pidLimitToF[0], pidLimitToF[1]);\r
+ }\r
+ else\r
+ {\r
+ driver.SetBkgInterval(bkgLimitM2[0], bkgLimitM2[1]);\r
+ driver.SetPidInterval(pidLimitM2[0], pidLimitM2[1]);\r
+ }\r
+ \r
driver.SetSecondaries(kSecondaries);\r
- driver.SetSecProd(kSecProd);\r
+ driver.SetSecProcedure(kSecProc);\r
driver.SetMatDCAxyModel(kMatDCAxyMod);\r
driver.SetAntiNucleusAsTemplate(kAntiNucTemplate);\r
driver.SetNBin(kNbin);\r
driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]);\r
driver.SetEfficiency(kEfficiency,0);\r
driver.SetFitFractionCorr(kFitFrac);\r
- driver.SetSysErr(kSysErr[0],kSysErr[1]);\r
\r
+ driver.SetDebugLevel(kDebugLevel);\r
\r
driver.SetMakeStats(makeStats);\r
driver.SetMakeCorrections(makeCor);\r
\r
if(!drawOutput) return 0;\r
\r
- TStyle* st = GetDrawingStyle();\r
- st->cd();\r
- gROOT->ForceStyle();\r
- \r
DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag());\r
\r
- if(kSecProd != 2) gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"\",\"Deuteron\", %d, %d, %f, %f)", driver.GetPtCorrDebugFilename().Data(), lowPtBin, hiPtBin, kDCAxy[0], kDCAxy[1]));\r
+ if(kSecProc == 0) gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"%s\",\"Deuteron\", %f, %f, %f, %f)", driver.GetPtCorrDebugFilename().Data(), driver.GetOutputCorrTag().Data(), ptmin, ptmax, kDCAxy[0], kDCAxy[1]));\r
\r
- DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, kPidM2, hiPtBin, kM2Bin[0], kM2Bin[1]);\r
+ DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, kPid, ptmax, ptpid);\r
DrawOutputRatio(outputRatio, outputTag, kSpecies);\r
DrawOutputSpectra(outputSpectra, outputTag, kSpecies);\r
\r
// LHC10x config for deuterons and antideuterons\r
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>\r
\r
+#if !defined(__CINT__) || defined(__MAKECINT__)\r
#include <TROOT.h>\r
#include <TString.h>\r
#include "AliLnDriver.h"\r
+#endif\r
+\r
#include "Config.h"\r
\r
-Int_t Config_Deuteron_TPC_LHC10x(const TString& inputDir = "~/alice/input",\r
- const TString& outputDir = "~/alice/output",\r
- const TString& period = "lhc10d",\r
- const TString& outputTag = "lhc10d",\r
- const TString& multTag = "",\r
- const TString& multCorTag = "",\r
- Bool_t inel = 1, // for mult\r
- Bool_t drawOutput = 1, // for batch\r
- Int_t lowPtBin = 3, // 0.4 Gev/c\r
- Int_t hiPtBin = 6, // 1.0 GeV/c\r
- Bool_t makeStats = 1,\r
- Bool_t makeCor = 1,\r
- Bool_t makePt = 1,\r
- Bool_t makeRatio = 1,\r
- Bool_t makeSpectra = 1 )\r
+Int_t Config_Deuteron_TPC_LHC10x( const TString& inputDir = "~/alice/input"\r
+ , const TString& outputDir = "~/alice/output"\r
+ , const TString& period = "lhc10d"\r
+ , const TString& outputTag = "lhc10d"\r
+ , const TString& trkselTag = "-tpc3-nsd-moc"\r
+ , const TString& multTag = ""\r
+ , const TString& multCorTag = ""\r
+ , Double_t ymax = 0.5\r
+ , Bool_t inel = 0\r
+ , Bool_t drawOutput = 1\r
+ , Double_t ptmin = 0.5 // GeV/c\r
+ , Double_t ptmax = 1.1 // GeV/c\r
+ , Bool_t makeStats = 1\r
+ , Bool_t makeCor = 1\r
+ , Bool_t makePt = 1\r
+ , Bool_t makeRatio = 1\r
+ , Bool_t makeSpectra = 1 )\r
{\r
//\r
// lhc10b, lhc10c, lhc10d, lhc10e config for deuterons and antideuterons\r
// (TPC)\r
//\r
const TString kSpecies = "Deuteron";\r
- const TString kTrkSel = "its_tpc_dca-tpc3";\r
- const TString kTrigName = "mbor";\r
+ const TString kTrkSel = "its_tpc_dca";\r
+ const TString kTrigName = "mband";\r
const Bool_t kMCtoINEL = 1;\r
- const Bool_t kUnfolding = 1;\r
- const Int_t kIter = 7;\r
+ const Double_t kPidEff = 1.;\r
const Bool_t kSecondaries = 1;\r
- const Int_t kSecProd = 0; // 0 tff, 1 roofit, 2 mc\r
+ const Int_t kSecProc = 0; // 0 tff, 1 mc\r
const Int_t kMatDCAxyMod = 1; // 0 geant, 1 flat\r
const Bool_t kAntiNucTemplate = 0;\r
const Int_t kNbin = 5;\r
- const Double_t kDCAxy[2] = {-0.2,0.2};\r
+ const Double_t kDCAxy[2] = {-0.2, 0.2};\r
const Bool_t kEfficiency = 1;\r
- const Bool_t kFitFrac = 0;\r
- const Double_t kSysErr[2] = {0.10, 0.11} ;\r
- \r
- Double_t xsec[3];\r
- GetInelXSection(xsec, period);\r
+ const Int_t kDebugLevel = 1;\r
\r
Double_t trigEff[3];\r
GetTriggerEfficiency(trigEff, kTrigName, period);\r
\r
// input and output filenames\r
\r
- TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + ".root";\r
- TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
- TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
- TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + "-corr.root";\r
+ TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + ".root";\r
+ TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";\r
+ TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";\r
+ TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + "-corr.root";\r
\r
TString outputPt = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";\r
TString outputRatio = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";\r
driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);\r
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);\r
\r
+ driver.SetRapidityInterval(-ymax,ymax);\r
+ \r
driver.SetOutputTag(outputTag);\r
+ driver.SetOutputCorTag(outputTag);\r
+ \r
driver.SetTriggerEfficiency(trigEff);\r
- driver.SetInelXSection(xsec);\r
driver.SetExtrapolateToINEL(inel);\r
driver.SetMCtoINEL(kMCtoINEL);\r
- driver.SetPtBinInterval(lowPtBin, hiPtBin);\r
- driver.SetPidM2(0);\r
- driver.SetUnfolding(kUnfolding, kIter);\r
+ driver.SetPtInterval(ptmin, ptmax);\r
+ driver.SetPid(0);\r
+ driver.SetPidEfficiency(kPidEff);\r
driver.SetSecondaries(kSecondaries);\r
- driver.SetSecProd(kSecProd);\r
+ driver.SetSecProcedure(kSecProc);\r
driver.SetMatDCAxyModel(kMatDCAxyMod);\r
driver.SetAntiNucleusAsTemplate(kAntiNucTemplate);\r
driver.SetNBin(kNbin);\r
driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]);\r
driver.SetEfficiency(kEfficiency,0);\r
- driver.SetFitFractionCorr(kFitFrac);\r
- driver.SetSysErr(kSysErr[0],kSysErr[1]);\r
\r
driver.SetMakeStats(makeStats);\r
driver.SetMakeCorrections(makeCor);\r
driver.SetMakeRatio(makeRatio);\r
driver.SetMakeSpectra(makeSpectra);\r
\r
+ driver.SetDebugLevel(kDebugLevel);\r
+ \r
driver.Run();\r
\r
// draw output\r
\r
if(!drawOutput) return 0;\r
\r
- TStyle* st = GetDrawingStyle();\r
- st->cd();\r
- gROOT->ForceStyle();\r
- \r
DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag());\r
\r
- if(kSecProd != 2) gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"\",\"Deuteron\", %d, %d, %f, %f)", driver.GetPtCorrDebugFilename().Data(), lowPtBin, hiPtBin, kDCAxy[0], kDCAxy[1]));\r
+ if(kSecProc == 0) gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"%s\",\"Deuteron\", %f, %f, %f, %f)", driver.GetPtCorrDebugFilename().Data(), driver.GetOutputCorrTag().Data(), ptmin, ptmax, kDCAxy[0], kDCAxy[1]));\r
\r
DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0);\r
DrawOutputRatio(outputRatio, outputTag, kSpecies);\r
// LHC10x config for protons and antiprotons
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TString.h>
#include "AliLnDriver.h"
+#endif
+
#include "Config.h"
-Int_t Config_Proton_TOF_LHC10x(const TString& inputDir = "~/alice/input",
- const TString& outputDir = "~/alice/output",
- const TString& period = "lhc10d",
- const TString& outputTag = "lhc10d",
- const TString& multTag = "",
- const TString& multCorTag = "",
- Bool_t inel = 1, // for mult
- Bool_t drawOutput = 1, // for batch
- Int_t lowPtBin = 11, // 1.0 GeV/c
- Int_t hiPtBin = 36, // 3.5 GeV/c
- Bool_t makeStats = 1,
- Bool_t makeCor = 1,
- Bool_t makePt = 1,
- Bool_t makeRatio = 1,
- Bool_t makeSpectra = 1 )
+Int_t Config_Proton_TOF_LHC10x( const TString& inputDir = "~/alice/input"
+ , const TString& outputDir = "~/alice/output"
+ , const TString& period = "lhc10d"
+ , const TString& outputTag = "lhc10d"
+ , const TString& trkselTag = "-bayes-nsd"
+ , const TString& multTag = ""
+ , const TString& multCorTag = ""
+ , Double_t ymax = 0.5
+ , Bool_t inel = 0 // for mult
+ , Bool_t drawOutput = 1 // for batch
+ , Double_t ptmin = 1.0
+ , Double_t ptmax = 3.5
+ , Double_t ptpid = 3.5
+ , Bool_t makeStats = 1
+ , Bool_t makeCor = 1
+ , Bool_t makePt = 1
+ , Bool_t makeRatio = 1
+ , Bool_t makeSpectra = 1 )
{
//
// lhc10b, lhc10c, lhc10d, lhc10e config for protons and antiprotons
// (TOF)
//
const TString kSpecies = "Proton";
- const TString kTrkSel = "its_tpc_tof_dca_spd-bayes";
+ const TString kTrkSel = "its_tpc_tof_dca_spd";
const TString kTrigName = "mbor";
const Bool_t kMCtoINEL = 1;
- const Bool_t kUnfolding = 0;
const Bool_t kSecondaries = 1;
- const Int_t kSecProd = 0; // 0 tff, 1 roofit, 2 mc
+ const Int_t kSecProc = 0; // 0 tff, 1 mc
const Int_t kMatDCAxyMod = 1; // 0 geant, 1 flat
- const Int_t kNbin = 15;
+ const Int_t kNbin = 10;
const Double_t kDCAxy[2] = {-1.,1.};
const Bool_t kEfficiency = 1;
- const Double_t kMatScaling = 1.9;
- const Double_t kFdwnScaling = 1.9;
- const Bool_t kFitFrac = 0;
+ const Double_t kMatScaling = 1;
+ const Double_t kFdwnScaling = 1;
const Bool_t kSameFdwn = 1;
- const Double_t kSysErr[2] = {0.08,0.08} ;
-
- Double_t xsec[3];
- GetInelXSection(xsec, period);
+ const Bool_t kPid = 0;
+ const Double_t kPidEff = 1;
+ const Int_t kDebugLevel = 1;
Double_t trigEff[3];
GetTriggerEfficiency(trigEff, kTrigName, period);
// input and output filenames
- TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + ".root";
- TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+multCorTag) + ".root";
- TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+multCorTag) + ".root";
- TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + "-corr.root";
+ TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + ".root";
+ TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";
+ TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";
+ TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + "-corr.root";
TString outputPt = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";
TString outputRatio = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";
driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
+ driver.SetRapidityInterval(-ymax,ymax);
+
driver.SetOutputTag(outputTag);
+ driver.SetOutputCorTag(outputTag);
+
driver.SetTriggerEfficiency(trigEff);
- driver.SetInelXSection(xsec);
driver.SetExtrapolateToINEL(inel);
driver.SetMCtoINEL(kMCtoINEL);
- driver.SetPtBinInterval(lowPtBin, hiPtBin);
- driver.SetPidM2(0);
- driver.SetUnfolding(kUnfolding);
+ driver.SetPtInterval(ptmin, ptmax);
+ driver.SetPidPt(ptpid);
+ driver.SetPid(kPid);
+ driver.SetPidEfficiency(kPidEff);
+
driver.SetSecondaries(kSecondaries);
- driver.SetSecProd(kSecProd);
+ driver.SetSecProcedure(kSecProc);
driver.SetMatDCAxyModel(kMatDCAxyMod);
driver.SetNBin(kNbin);
driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]);
driver.SetEfficiency(kEfficiency);
driver.SetScalingFactors(kMatScaling, kFdwnScaling);
- driver.SetFitFractionCorr(kFitFrac);
driver.SetSameFeedDownCorr(kSameFdwn);
- driver.SetSysErr(kSysErr[0],kSysErr[1]);
-
- //driver.PrintFilenames();
+ driver.SetDebugLevel(kDebugLevel);
driver.SetMakeStats(makeStats);
driver.SetMakeCorrections(makeCor);
if(!drawOutput) return 0;
- TStyle* st = GetDrawingStyle();
- st->cd();
- gROOT->ForceStyle();
-
DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag());
- if(kSecProd != 2) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, lowPtBin, hiPtBin, kDCAxy[0], kDCAxy[1]);
+ if(kSecProc == 0) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, ptmin, ptmax, kDCAxy[0], kDCAxy[1]);
DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0);
DrawOutputRatio(outputRatio, outputTag, kSpecies);
// LHC10x config for protons and antiprotons
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TString.h>
#include "AliLnDriver.h"
+#endif
+
#include "Config.h"
-Int_t Config_Proton_TPC_LHC10x(const TString& inputDir = "~/alice/input",
- const TString& outputDir = "~/alice/output",
- const TString& period = "lhc10d",
- const TString& outputTag = "lhc10d",
- const TString& multTag = "",
- const TString& multCorTag = "",
- Bool_t inel = 1, // for mult
- Bool_t drawOutput = 1, // for batch
- Int_t lowPtBin = 5, // 0.4 GeV/c
- Int_t hiPtBin = 11, // 1.0 GeV/c
- Bool_t makeStats = 1,
- Bool_t makeCor = 1,
- Bool_t makePt = 1,
- Bool_t makeRatio = 1,
- Bool_t makeSpectra = 1 )
+Int_t Config_Proton_TPC_LHC10x( const TString& inputDir = "~/alice/input"
+ , const TString& outputDir = "~/alice/output"
+ , const TString& period = "lhc10d"
+ , const TString& outputTag = "lhc10d"
+ , const TString& trkselTag = "-bayes-nsd"
+ , const TString& multTag = ""
+ , const TString& multCorTag = ""
+ , Double_t ymax = 0.5
+ , Bool_t inel = 0 // for mult
+ , Bool_t drawOutput = 1 // for batch
+ , Double_t ptmin = 0.4
+ , Double_t ptmax = 1.0
+ , Bool_t makeStats = 1
+ , Bool_t makeCor = 1
+ , Bool_t makePt = 1
+ , Bool_t makeRatio = 1
+ , Bool_t makeSpectra = 1 )
{
//
// lhc10b, lhc10c, lhc10d, lhc10e config for protons and antiprotons
// (TPC)
//
const TString kSpecies = "Proton";
- const TString kTrkSel = "its_tpc_dca_spd-bayes";
+ const TString& kTrkSel = "its_tpc_dca_spd";
const TString kTrigName = "mbor";
const Bool_t kMCtoINEL = 1;
- const Bool_t kUnfolding = 0;
- const Int_t kIter = 5;
const Bool_t kSecondaries = 1;
- const Int_t kSecProd = 0; // 0 tff, 1 roofit, 2 mc
+ const Int_t kSecProc = 0; // 0 tff, 1 mc
const Int_t kMatDCAxyMod = 1; // 0 geant, 1 flat
- const Int_t kNbin = 10;
+ const Int_t kNbin = 12;
const Double_t kDCAxy[2] = {-1.,1.};
const Bool_t kEfficiency = 1;
const Bool_t kG3Fluka = 0;
- const Double_t kMatScaling = 1.9;
- const Double_t kFdwnScaling = 1.9;
- const Bool_t kFitFrac = 0;
+ const Double_t kMatScaling = 1;
+ const Double_t kFdwnScaling = 1;
const Bool_t kSameFdwn = 1;
- const Double_t kSysErr[2] = {0.08,0.08} ;
-
- Double_t xsec[3];
- GetInelXSection(xsec, period);
+ const Int_t kDebugLevel = 1;
Double_t trigEff[3];
GetTriggerEfficiency(trigEff, kTrigName, period);
// input and output filenames
- TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + ".root";
- TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+multCorTag) + ".root";
- TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+multCorTag, kG3Fluka) + ".root";
- TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + "-corr.root";
+ TString inputData = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + ".root";
+ TString inputSimu = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+trkselTag+multCorTag) + ".root";
+ TString inputSimuFix = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+trkselTag+multCorTag, kG3Fluka) + ".root";
+ TString inputCorr = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+trkselTag+multTag) + "-corr.root";
TString outputPt = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";
TString outputRatio = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";
driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
+ driver.SetRapidityInterval(-ymax,ymax);
+
driver.SetOutputTag(outputTag);
+ driver.SetOutputCorTag(outputTag);
+
driver.SetTriggerEfficiency(trigEff);
- driver.SetInelXSection(xsec);
driver.SetExtrapolateToINEL(inel);
driver.SetMCtoINEL(kMCtoINEL);
- driver.SetPtBinInterval(lowPtBin, hiPtBin);
- driver.SetPidM2(0);
- driver.SetUnfolding(kUnfolding, kIter);
+ driver.SetPtInterval(ptmin, ptmax);
+ driver.SetPid(0);
driver.SetSecondaries(kSecondaries);
- driver.SetSecProd(kSecProd);
+ driver.SetSecProcedure(kSecProc);
driver.SetMatDCAxyModel(kMatDCAxyMod);
driver.SetNBin(kNbin);
driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]);
driver.SetEfficiency(kEfficiency,kG3Fluka);
driver.SetScalingFactors(kMatScaling, kFdwnScaling);
- driver.SetFitFractionCorr(kFitFrac);
driver.SetSameFeedDownCorr(kSameFdwn);
- driver.SetSysErr(kSysErr[0],kSysErr[1]);
+
+ driver.SetDebugLevel(kDebugLevel);
driver.SetMakeStats(makeStats);
driver.SetMakeCorrections(makeCor);
if(!drawOutput) return 0;
- TStyle* st = GetDrawingStyle();
- st->cd();
- gROOT->ForceStyle();
-
DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag());
- if(kSecProd != 2) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, lowPtBin, hiPtBin, kDCAxy[0], kDCAxy[1]);
+ if(kSecProc == 0) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, ptmin, ptmax, kDCAxy[0], kDCAxy[1]);
DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0);
DrawOutputRatio(outputRatio, outputTag, kSpecies);
// TPC+TOF config
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TSystem.h>
#include <TString.h>
#include <TFileMerger.h>
-
#include "AliLnDriver.h"
+#endif
+
#include "Config.h"
-Int_t Config_TPCTOF_LHC10x(const TString& inputDir = "~/alice/input",
- const TString& outputDir = "~/alice/output",
- const TString& period = "lhc10d",
- const TString& outputTag = "lhc10d",
- const TString& multTag = "",
- const TString& multCorTag = "",
- Bool_t inel = 1, // for mult
- Bool_t drawOutput = 1, // for batch
- const TString& species = "Proton",
- Int_t lowPtBin = 5,
- Int_t jointPtBin = 11,
- Int_t hiPtBin = 36)
+Int_t Config_TPCTOF_LHC10x( const TString& inputDir = "~/alice/input"
+ , const TString& outputDir = "~/alice/output"
+ , const TString& period = "lhc10d"
+ , const TString& outputTag = "lhc10d"
+ , const TString& trkselTag = "-bayes-nsd"
+ , const TString& multTag = ""
+ , const TString& multCorTag = ""
+ , Double_t ymax = 0.5
+ , Bool_t inel = 0
+ , Bool_t drawOutput = 1
+ , const TString& species = "Proton"
+ , Double_t ptmin = 0.4
+ , Double_t ptjoint = 1.0
+ , Double_t ptmax = 3.2
+ , Double_t ptpid = 4)
{
//
// combine TPC and TOF for protons and deuterons
//
- 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"))
+ "\"" + outputDir + "\","
+ "\"" + period + "\","
+ "\"" + kOutputTagTPC + "\","
+ + "\"" + trkselTag + "\","
+ "\"" + multTag + "\","
+ "\"" + multCorTag;
+ "\"" + outputDir + "\","
+ "\"" + period + "\","
+ "\"" + kOutputTagTOF + "\","
+ + "\"" + trkselTag + "\","
+ "\"" + multTag + "\","
+ "\"" + multCorTag;
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)"
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %f, %d, 0, %f, %f, 0,1,1,0,0)"
, species.Data()
, kArgTPC.Data()
+ , ymax
, inel
- , lowPtBin
- , jointPtBin));
+ , ptmin
+ , ptjoint));
+
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)"
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %f, %d, 0, %f, %f, %f, 1,1,1,0,0)"
, species.Data()
, kArgTOF.Data()
+ , ymax
, inel
- , jointPtBin
- , hiPtBin));
+ , ptjoint
+ , ptmax
+ , ptpid));
TString outputPtTPC = outputDir + "/" + MakeOutputName(species, kOutputTagTPC) + "-Pt.root";
TString outputPtTOF = outputDir + "/" + MakeOutputName(species, kOutputTagTOF) + "-Pt.root";
+ TString outputPtTPCdbg = outputDir + "/" + MakeOutputName(species, kOutputTagTPC) + "-Pt-debug.root";
+ TString outputPtTOFdbg = outputDir + "/" + MakeOutputName(species, kOutputTagTOF) + "-Pt-debug.root";
+
// combine TPC and TOF pt
TString outputPt = outputDir + "/" + MakeOutputName(species, outputTag) + "-Pt.root";
m.Merge();
+ // remove tmp files
+ gSystem->Exec(Form("rm -f %s %s %s %s", outputPtTPC.Data(), outputPtTOF.Data(), outputPtTPCdbg.Data(), outputPtTOFdbg.Data()));
+
// make ratio and spectra
AliLnDriver driver;
driver.SetSpecies(species);
- Double_t xsec[3];
- GetInelXSection(xsec, period);
+ driver.SetRapidityInterval(-ymax,ymax);
- driver.SetInelXSection(xsec);
driver.SetExtrapolateToINEL(inel);
-
- if(species == "Proton") driver.SetSysErr(kProtonSysErr[0],kProtonSysErr[1]);
- if(species == "Deuteron") driver.SetSysErr(kDeuteronSysErr[0],kDeuteronSysErr[1]);
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
if(!drawOutput) return 0;
- TStyle* st = GetDrawingStyle();
- st->cd();
- gROOT->ForceStyle();
-
DrawOutputRatio(outputRatio, outputTag, species);
DrawOutputSpectra(outputSpectra, outputTag, species);
// Draw B2
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <TROOT.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TStyle.h>
#include <TFile.h>
#include <TString.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
-#include <TStyle.h>
+#endif
#include "B2.h"
//
// Draw B2
//
- TStyle* st = new TStyle();
-
- st->SetPadTickX(1);
- st->SetPadTickY(1);
- st->SetPadGridX(1);
- st->SetPadGridY(1);
- st->SetCanvasColor(0);
- st->SetFrameBorderMode(0);
- st->SetFrameFillColor(0);
- st->SetLabelFont(62,"XYZ");
- st->SetTitleFont(62,"XYZ");
-
- st->cd();
- gROOT->ForceStyle();
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
// B2
- TGraphErrors* grB2Pt = (TGraphErrors*)FindObj(finput, tag, Form("B2%s_Pt",suffix.Data()));
- TGraphErrors* grSysErrB2Pt = (TGraphErrors*)FindObj(finput, tag, Form("B2%s_SysErr_Pt",suffix.Data()));
+ TGraphErrors* grB2Pt = FindObj<TGraphErrors>(finput, tag, Form("B2%s_Pt",suffix.Data()));
c0->cd(1);
gPad->SetLogy(1);
- grSysErrB2Pt->GetXaxis()->SetTitle("p_{T}/A (GeV/c)");
- grSysErrB2Pt->GetYaxis()->SetTitle("B_{2} (GeV^{2}c^{-3})");
- grSysErrB2Pt->SetLineColor(kRed);
- grSysErrB2Pt->SetFillStyle(1001);
- grSysErrB2Pt->SetFillColor(kRed-10);
- grSysErrB2Pt->Draw("A3");
-
+ grB2Pt->GetXaxis()->SetTitle("#it{p}_{T}/A (GeV/#it{c})");
+ grB2Pt->GetYaxis()->SetTitle("#it{B}_{2} (GeV^{2}#it{c}^{-3})");
grB2Pt->SetLineColor(kRed);
grB2Pt->SetMarkerStyle(kFullCircle);
grB2Pt->SetMarkerColor(kRed);
- grB2Pt->Draw("P");
+ grB2Pt->Draw("zAP");
// homogeneity volume
- TGraphErrors* grR3Pt = (TGraphErrors*)FindObj(finput, tag, Form("R3%s_Pt",suffix.Data()));
- TGraphErrors* grSysErrR3Pt = (TGraphErrors*)FindObj(finput, tag, Form("R3%s_SysErr_Pt",suffix.Data()));
+ TGraphErrors* grR3Pt = FindObj<TGraphErrors>(finput, tag, Form("R3%s_Pt",suffix.Data()));
c0->cd(2);
- grSysErrR3Pt->GetXaxis()->SetTitle("p_{T}/A (GeV/c)");
- grSysErrR3Pt->GetYaxis()->SetTitle("R_{side}^{2} R_{long} (fm^{3})");
- grSysErrR3Pt->GetYaxis()->SetTitleOffset(1.3);
- grSysErrR3Pt->SetLineColor(kRed);
- grSysErrR3Pt->SetFillStyle(1001);
- grSysErrR3Pt->SetFillColor(kRed-10);
- grSysErrR3Pt->Draw("A3");
-
+ grR3Pt->GetXaxis()->SetTitle("#it{p}_{T}/A (GeV/#it{c})");
+ grR3Pt->GetYaxis()->SetTitle("#it{R}_{side}^{2} #it{R}_{long} (fm^{3})");
+ grR3Pt->GetYaxis()->SetTitleOffset(1.3);
grR3Pt->SetLineColor(kRed);
grR3Pt->SetMarkerStyle(kFullCircle);
grR3Pt->SetMarkerColor(kRed);
- grR3Pt->Draw("P");
+ grR3Pt->Draw("zAP");
}
// Draw pt corrections for debugging
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <Riostream.h>
-#include <TROOT.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TStyle.h>
#include <TFile.h>
#include <TH1D.h>
-#include <TH2D.h>
#include <TString.h>
#include <TCanvas.h>
#include <TF1.h>
+#endif
#include "B2.h"
-void DrawPair(TH1* hX, TH1* hY, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, const char* title="", const char* option="E", Int_t xMarker=kFullCircle, Int_t yMarker=kFullCircle, Int_t xColor=kBlue, Int_t yColor=kRed);
+void DrawPair(TH1* hX, TH1* hY, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, const TString& title="", const char* option="E", Int_t xMarker=kFullCircle, Int_t yMarker=kFullCircle, Int_t xColor=kBlue, Int_t yColor=kRed);
void DrawCorr(const TString& species="Deuteron", const TString& inputFile="corrections.root", const TString& tag="")
{
//
// Draw pt corrections for debugging
//
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetOptStat(0);
+
Double_t xmin = 0;
Double_t xmax = 3.5;
const TString kPrefix[] = {"", "Anti"};
- // Response matrix
-
- TCanvas* c0 = new TCanvas(Form("%s.Unfolding", species.Data()), Form("Response matrix for (Anti)%ss", species.Data()));
-
- c0->Divide(2,1);
-
- for(Int_t i=0; i<kNpart; ++i)
- {
- c0->cd(i+1);
-
- TH2D* hResponseMtx = (TH2D*)FindObj(finput, tag, Form("%s%s_Response_Matrix",kPrefix[i].Data(), species.Data()));
- hResponseMtx->SetAxisRange(xmin, xmax, "X");
- hResponseMtx->SetAxisRange(xmin, xmax, "Y");
- hResponseMtx->GetYaxis()->SetTitleOffset(1.3);
- hResponseMtx->DrawCopy("cont0colZ");
- }
-
// Reconstruction efficiency
TCanvas* c2 = new TCanvas(Form("%s.Efficiency",species.Data()), Form("Reconstruction Efficiency for (Anti)%ss",species.Data()));
for(Int_t i=0; i<kNpart; ++i)
{
- hEffTrigPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_Trig_Pt");
- hEffVtxPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_Vtx_Pt");
- hEffAccTrkPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_AccTrk_Pt");
+ hEffTrigPt[i] = FindObj<TH1D>(finput, tag, kPrefix[i] + species + "_Eff_Trig_Pt");
+ hEffVtxPt[i] = FindObj<TH1D>(finput, tag, kPrefix[i] + species + "_Eff_Vtx_Pt");
+ hEffAccTrkPt[i] = FindObj<TH1D>(finput, tag, kPrefix[i] + species + "_Eff_AccTrk_Pt");
}
c2->cd(1);
for(Int_t i=0; i<kNpart; ++i)
{
- TH1D* hFracFdwnPt = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Pt");
- TH1D* hFracMatPt = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Mat_Pt");
+ TH1D* hFracFdwnPt = FindObj<TH1D>(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Pt");
+ TH1D* hFracMatPt = FindObj<TH1D>(finput, tag, kPrefix[i] + species + "_Frac_Mat_Pt");
- TF1* fncFracFdwnPt = (TF1*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Fit_Pt");
- TF1* fncFracMatPt = (TF1*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Mat_Fit_Pt");
+ TF1* fncFracFdwnPt = FindObj<TF1>(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Fit_Pt");
+ TF1* fncFracMatPt = FindObj<TF1>(finput, tag, kPrefix[i] + species + "_Frac_Mat_Fit_Pt");
c3->cd(2*i+1);
hFracFdwnPt->SetAxisRange(xmin, xmax, "X");
}
}
-void DrawPair(TH1* hX, TH1* hY, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, const char* title, const char* option, Int_t xMarker, Int_t yMarker, Int_t xColor, Int_t yColor)
+void DrawPair(TH1* hX, TH1* hY, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, const TString& title, const char* option, Int_t xMarker, Int_t yMarker, Int_t xColor, Int_t yColor)
{
//
// Draw a pair of histograms in the current pad
//
- hX->SetTitle(title);
+ hX->SetTitle(title.Data());
hX->SetAxisRange(xmin, xmax, "X");
hX->SetAxisRange(ymin, ymax, "Y");
hX->SetMarkerColor(xColor);
// Draw histograms/graphs with same name located in different directories
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <TROOT.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TFile.h>
#include <TString.h>
#include <TMath.h>
#include <TKey.h>
#include <TLine.h>
#include <TStyle.h>
+#endif
TGraphErrors* Divide(const TGraphErrors* grX1, const TGraphErrors* grX2, const TString& name);
const Int_t kColor[kMax] = { kRed, kBlue, kOrange+1, kGreen-2, kGreen+2, kAzure, kViolet+10, kAzure+2, kOrange+2, kSpring-7 };
- const Int_t kMarker[kMax] = { kFullCircle, kFullCircle, kFullCircle, kFullCircle, kOpenCircle, kOpenSquare, kOpenTriangleUp, kOpenDiamond, kOpenCross, kFullStar };
+ const Int_t kMarker[kMax] = { kFullCircle, kFullTriangleUp, kFullSquare, kFullTriangleDown, kOpenCircle, kOpenSquare, kOpenTriangleUp, kOpenDiamond, kOpenCross, kFullStar };
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
// draw
- TStyle* st = new TStyle();
-
- st->SetOptTitle(0);
- st->SetOptStat(0);
-
- st->SetPadTickX(1);
- st->SetPadTickY(1);
- st->SetPadGridX(1);
- st->SetPadGridY(1);
-
- st->SetCanvasColor(0);
- st->SetFrameBorderMode(0);
- st->SetFrameFillColor(0);
- st->SetTitleFillColor(0);
- st->SetLabelFont(62,"XYZ");
- st->SetTitleFont(62,"XYZ");
-
- st->cd();
- gROOT->ForceStyle();
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetOptStat(0);
TCanvas* c0 = new TCanvas(canvasName.Data(), canvasTitle.Data());
}
else if(obj[i]->InheritsFrom("TGraph"))
{
- obj[i]->Draw("PZ");
+ obj[i]->Draw("zP");
}
}
TH1F* frm = c1->DrawFrame(xmin, 0.5 ,xmax, 1.5);
frm->GetXaxis()->SetTitle(xtitle);
frm->GetYaxis()->SetTitle("Ratio");
+ frm->GetYaxis()->SetTitleOffset(1.20);
}
// draw comparison
if(nDir>0 && option > 0)
{
- for(Int_t i=1; i<nDir; ++i)
- {
- grDiv[i]->Draw("PZ");
- }
-
// draw a red line for the reference
TLine* ref = new TLine(xmin,1,xmax,1);
ref->SetLineWidth(1);
ref->SetLineColor(kColor[0]);
+ ref->SetLineStyle(2);
ref->Draw();
+ for(Int_t i=1; i<nDir; ++i)
+ {
+ grDiv[i]->Draw("zP");
+ }
+
if(option == 2)
{
TLegend* legendRatio = new TLegend(0.5718391,0.6991525,0.8390805,0.8368644,0,"brNDC");
legendRatio->AddEntry(grDiv[i], subdir[i]->Data(), "lp");
}
+ legendRatio->SetTextFont(42);
legendRatio->Draw();
}
}
// Draw corrected pt
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <Riostream.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TROOT.h>
+#include <TStyle.h>
#include <TFile.h>
#include <TH1D.h>
#include <TString.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TLegend.h>
-
#include <RooWorkspace.h>
#include <RooMsgService.h>
#include <RooPlot.h>
#include <RooRealVar.h>
-#include <RooAbsPdf.h>
#include <RooAbsData.h>
+#include <RooAbsPdf.h>
+#endif
#include "B2.h"
-void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="AntiDeuteron", Int_t hiptbin=17, Bool_t m2pid=0, Int_t lowm2bin=9, Int_t him2bin=17)
+void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="AntiDeuteron", Double_t ptmax=3., Bool_t m2pid=0, Double_t ptpid=1.2)
{
//
// Draw corrected pt for debugging
Double_t xmin = 0;
Double_t xmax = 3.5;
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptTitle(1);
+
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
+ TH1D* hPidPt = FindObj<TH1D>(finput, tag, Form("%s_PID_Pt",particle.Data()));
+
+ Int_t hiptbin = hPidPt->GetXaxis()->FindFixBin(ptmax);
+ Int_t lowm2bin = hPidPt->GetXaxis()->FindFixBin(ptpid);
+ Int_t him2bin = hPidPt->GetXaxis()->FindFixBin(ptmax);
+
// m2 data fitted models
if(m2pid && (hiptbin>lowm2bin))
{
- TH1D* hPidPt = (TH1D*)FindObj(finput, tag, Form("%s_PID_Pt",particle.Data()));
- Double_t binwidth = hPidPt->GetBinWidth(0);
-
using namespace RooFit;
// disable verbose in RooFit
RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
TCanvas* c0 = new TCanvas(Form("%s.M2",particle.Data()), Form("M2 for %ss",particle.Data()));
- c0->Divide(3,3);
+ c0->Divide(4,5);
TGraphErrors* grFitSigmaPt = new TGraphErrors();
grFitSigmaPt->SetName(Form("%s_Fit_Sigma_Pt", particle.Data()));
TGraphErrors* grFitM2Pt = new TGraphErrors();
grFitM2Pt->SetName(Form("%s_Fit_M2_Pt", particle.Data()));
- for(Int_t i=lowm2bin, j=0; i<him2bin && i-lowm2bin < 9 && i < hiptbin; ++i)
+ for(Int_t i=lowm2bin, j=0; i<him2bin && i-lowm2bin < 20 && i < hiptbin; ++i)
{
c0->cd(i-lowm2bin+1);
- //gPad->SetLogy(0);
- RooWorkspace* w= (RooWorkspace*)FindObj(finput, tag, Form("%s_M2_%02d",particle.Data(),i));
+ RooWorkspace* w= FindObj<RooWorkspace>(finput, tag, Form("%s_M2_%02d",particle.Data(),i));
- RooPlot* m2frame = w->var("m2")->frame();
+ RooPlot* m2frame = w->var("x")->frame();
w->data("data")->plotOn(m2frame);
- w->pdf("model")->plotOn(m2frame, LineWidth(2));
- w->pdf("model")->plotOn(m2frame, Components(*(w->pdf("Sd"))),LineWidth(1), LineColor(9));
+ w->pdf("model")->plotOn(m2frame, Components(*(w->pdf("Sd"))),LineWidth(1), LineColor(8));
w->pdf("model")->plotOn(m2frame, Components(*(w->pdf("Bkg"))),LineWidth(1), LineColor(46),LineStyle(kDashed));
+ w->pdf("model")->plotOn(m2frame, LineWidth(1));
- Double_t ptMin = (i-1)*binwidth;
- Double_t ptMax = i*binwidth;
- Double_t pt = (ptMin+ptMax)/2.;
-
- m2frame->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", ptMin, ptMax));
+ m2frame->SetTitle(Form("%0.2f < #it{p}_{T} < %0.2f GeV/#it{c}", hPidPt->GetBinLowEdge(i), hPidPt->GetBinLowEdge(i)+hPidPt->GetBinWidth(i)));
m2frame->SetMinimum(0.2);
-
+
m2frame->Draw();
- grFitSigmaPt->SetPoint(j, pt, w->var("width")->getVal());
- grFitSigmaPt->SetPointError(j, 0, w->var("width")->getError());
+ Double_t pt = hPidPt->GetBinCenter(i);
+
+ grFitSigmaPt->SetPoint(j, pt, w->var("sigma")->getVal());
+ grFitSigmaPt->SetPointError(j, 0, w->var("sigma")->getError());
- grFitM2Pt->SetPoint(j, pt, w->var("mean")->getVal());
- grFitM2Pt->SetPointError(j++, 0, w->var("mean")->getError());
+ grFitM2Pt->SetPoint(j, pt, w->var("mu")->getVal());
+ grFitM2Pt->SetPointError(j++, 0, w->var("mu")->getError());
}
+ c0->Update();
+
// model parameters
TCanvas* c1 = new TCanvas(Form("%s.M2.FitParameters",particle.Data()), Form("M2 model parameters for %ss",particle.Data()));
c1->Divide(2,1);
c1->cd(1);
- gPad->SetGridx();
- gPad->SetGridy();
grFitSigmaPt->SetMarkerStyle(kFullCircle);
grFitSigmaPt->SetMarkerColor(kBlue);
grFitSigmaPt->SetLineColor(kBlue);
grFitSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
- grFitSigmaPt->GetYaxis()->SetTitle("#sigma_{m^{2}}");
+ grFitSigmaPt->GetYaxis()->SetTitle("#sigma");
grFitSigmaPt->Draw("ALP");
c1->cd(2);
- gPad->SetGridx();
- gPad->SetGridy();
grFitM2Pt->SetMarkerStyle(kFullCircle);
grFitM2Pt->SetMarkerColor(kBlue);
grFitM2Pt->SetLineColor(kBlue);
grFitM2Pt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
- grFitM2Pt->GetYaxis()->SetTitle("m^{2} (GeV^{2}/c^{4})");
+ grFitM2Pt->GetYaxis()->SetTitle("#mu");
grFitM2Pt->Draw("ALP");
//delete grFitSigmaPt;
//delete grFitM2Pt;
+
+ c1->Update();
}
// remaining corrections
- const Int_t kNum = 5;
- const TString kCorr[kNum] = { "PID", "M2Corr", "SecCor", "Unfolded", "EffCor"};
- const TString kLabel[kNum] = { "Raw", "PID contamination", "Secondaries", "Unfolding","Efficiency" };
- const Int_t kColor[kNum] = { kRed, kAzure, kOrange+1, kGreen-2, kGreen-3};
- const Int_t kMarker[kNum] = { kFullCircle, kOpenCircle, kFullTriangleUp, kOpenTriangleUp, kFullCircle};
+ TCanvas* c2 = new TCanvas(Form("%s.Pt",particle.Data()), Form("Pt for %s",particle.Data()));
+ c2->SetLogy();
+
+ const Int_t kNum = 4;
+ const TString kCorr[kNum] = { "PID", "PidCorr", "SecCorr", "EffCorr"};
+ const TString kLabel[kNum] = { "Raw", "PID", "Secondaries","Efficiency" };
+ const Int_t kColor[] = { kRed, kAzure, kOrange+1, kGreen-3, kGreen-2};
+ const Int_t kMarker[] = { kFullCircle, kFullSquare, kFullTriangleUp, kFullTriangleDown, kFullCircle, kOpenTriangleUp};
TLegend* legend = new TLegend(0.5689655,0.6355932,0.8362069,0.8326271,0,"brNDC");
legend->SetTextSize(0.03);
for(Int_t i=0; i<kNum; ++i)
{
- hPt[i] = (TH1D*)FindObj(finput, tag, particle + "_" + kCorr[i] + "_Pt");
+ hPt[i] = FindObj<TH1D>(finput, tag, particle + "_" + kCorr[i] + "_Pt");
hPt[i]->SetLineColor(kColor[i]);
hPt[i]->SetMarkerColor(kColor[i]);
hPt[i]->SetMarkerStyle(kMarker[i]);
legend->AddEntry(hPt[i], kLabel[i], "lp");
}
- TCanvas* c2 = new TCanvas(Form("%s.Pt",particle.Data()), Form("Pt for %s",particle.Data()));
- c2->SetLogy();
-
hPt[kNum-1]->SetTitle(particle.Data());
hPt[kNum-1]->SetAxisRange(xmin,xmax,"X");
hPt[kNum-1]->Draw("E");
- for(Int_t i=0; i<kNum-1; ++i) hPt[i]->Draw("sameE");
+ for(Int_t i=0; i<kNum-1; ++i) hPt[i]->Draw("sameE");
legend->Draw();
+
+ c2->Update();
}
// Draw antiparticle/particle ratio
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <Riostream.h>
-#include <TROOT.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TStyle.h>
#include <TFile.h>
#include <TH1D.h>
#include <TString.h>
#include <TCanvas.h>
+#endif
#include "B2.h"
Double_t ymin = 0.73;
Double_t ymax = 1.14;
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetOptTitle(0);
+ gStyle->SetOptStat(0);
+
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
else if(species=="He3") ytitle = "#bar{He3}/He3";
else if(species=="Alpha") ytitle = "#bar{#alpha}/#alpha";
- TH1D* hRatioPt = (TH1D*)FindObj(finput, tag, Form("Anti%s%s_Ratio_Pt", species.Data(), species.Data()));
+ TH1D* hRatioPt = FindObj<TH1D>(finput, tag, Form("Anti%s%s_Ratio_Pt", species.Data(), species.Data()));
TCanvas* c0 = new TCanvas(Form("c0.Ratio%s",species.Data()), Form("Anti%s/%s ratio", species.Data(), species.Data()));
c0->cd();
// Draw corrections of secondaries
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <Riostream.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TROOT.h>
+#include <TStyle.h>
#include <TFile.h>
-#include <TSystem.h>
#include <TH1D.h>
#include <TString.h>
#include <TMath.h>
#include <TCanvas.h>
+#endif
#include "B2.h"
-
Int_t GetNumberCols(Int_t lowbin, Int_t hibin);
-TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
+TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax);
-TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
+TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax);
-TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax);
+TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax);
-void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="Proton", Int_t lowbin=1, Int_t hibin=10, Double_t dcaxyMin=-1.5, Double_t dcaxyMax=1.5)
+void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="Proton", Double_t ptmin=0.5, Double_t ptmax=3., Double_t dcaxyMin=-1.5, Double_t dcaxyMax=1.5)
{
//
// draw corrections of secondaries
//
- const Int_t kNBin = hibin;
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptTitle(1);
+
+ const Int_t kNBin = 50;
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
TH1D* hMatFit[kNBin];
TH1D* hFdwnFit[kNBin];
+ TH1D* hPidPt = FindObj<TH1D>(finput, tag, particle + "_PID_Pt");
+
+ Int_t lowbin = hPidPt->GetXaxis()->FindFixBin(ptmin);
+ Int_t hibin = hPidPt->GetXaxis()->FindFixBin(ptmax);
+
for(Int_t i=lowbin; i<hibin && i<kNBin; ++i)
{
// MC
- hData[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Data_DCAxy_%02d",i));
- hDataMC[i] = (TH1D*)FindObj(finput, tag, particle + Form("_SimData_DCAxy_%02d",i));
- hPrim[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Prim_DCAxy_%02d",i));
- hMat[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Mat_DCAxy_%02d",i));
+ hData[i] = FindObj<TH1D>(finput, tag, particle + Form("_Data_DCAxy_%02d",i));
+ hDataMC[i] = FindObj<TH1D>(finput, tag, particle + Form("_SimData_DCAxy_%02d",i));
+ hPrim[i] = FindObj<TH1D>(finput, tag, particle + Form("_Prim_DCAxy_%02d",i));
+ hMat[i] = FindObj<TH1D>(finput, tag, particle + Form("_Mat_DCAxy_%02d",i));
if(particle.Contains("Proton"))
{
- hFdwn[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fdwn_DCAxy_%02d",i));
+ hFdwn[i] = FindObj<TH1D>(finput, tag, particle + Form("_Fdwn_DCAxy_%02d",i));
}
else
{
}
// fitted data
- hDataFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Data_DCAxy_%02d",i));
- hPrimFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Prim_DCAxy_%02d",i));
- hMatFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Mat_DCAxy_%02d",i));
+ hDataFit[i] = FindObj<TH1D>(finput, tag, particle + Form("_Fit_Data_DCAxy_%02d",i));
+ hPrimFit[i] = FindObj<TH1D>(finput, tag, particle + Form("_Fit_Prim_DCAxy_%02d",i));
+ hMatFit[i] = FindObj<TH1D>(finput, tag, particle + Form("_Fit_Mat_DCAxy_%02d",i));
if(particle.Contains("Proton"))
{
- hFdwnFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Fdwn_DCAxy_%02d",i));
+ hFdwnFit[i] = FindObj<TH1D>(finput, tag, particle + Form("_Fit_Fdwn_DCAxy_%02d",i));
}
else
{
}
}
- // pt bin width
- TH1D* hPidPt = (TH1D*)FindObj(finput, tag, particle + "_PID_Pt");
- Double_t binwidth = hPidPt->GetBinWidth(0);
-
// draw
TCanvas* c0 = 0;
if(hDataMC[lowbin]!=0)
{
- c0 = DrawMC(hDataMC, hPrim, hMat, hFdwn, lowbin, hibin, particle + ".MC", Form("MC simulation for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
+ c0 = DrawMC(hDataMC, hPrim, hMat, hFdwn, lowbin, hibin, particle + ".MC", Form("MC simulation for %s",particle.Data()), hPidPt, dcaxyMin, dcaxyMax);
+ c0->Update();
}
TCanvas* c1;
if(hDataFit[lowbin]!=0)
{
- c1 = DrawFit(hData, hDataFit, hPrimFit, hMatFit, hFdwnFit, lowbin, hibin, particle + ".Fit", Form("Fitted DCA for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
+ c1 = DrawFit(hData, hDataFit, hPrimFit, hMatFit, hFdwnFit, lowbin, hibin, particle + ".Fit", Form("Fitted DCA for %s",particle.Data()), hPidPt, dcaxyMin, dcaxyMax);
+ c1->Update();
}
TCanvas* c2;
if(hDataFit[lowbin]!=0)
{
- c2 = DrawGoF(hData, hDataFit, lowbin, hibin, particle + ".GoF", Form("GoF for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
+ c2 = DrawGoF(hData, hDataFit, lowbin, hibin, particle + ".GoF", Form("GoF for %s",particle.Data()), hPidPt, dcaxyMin, dcaxyMax);
+ c2->Update();
}
}
return TMath::CeilNint(TMath::Sqrt(TMath::Abs(hibin - lowbin)));
}
-TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
+TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax)
{
//
// Draw MC DCA distributions for each pt bin
{
c0->cd(i+1-lowbin);
- gPad->SetGrid();
gPad->SetLogy();
- hData[i]->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
+ hData[i]->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", hPt->GetBinLowEdge(i), hPt->GetBinLowEdge(i)+hPt->GetBinWidth(i)));
hData[i]->SetMinimum(0.2);
hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
return c0;
}
-TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
+TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax)
{
//
// Draw DCA fitted model for each pt bin in current canvas
TCanvas* c0 = new TCanvas(name.Data(), title.Data());
c0->Divide(ncol, ncol);
- for(Int_t i=lowbin; i < hibin; ++i)
+ for(Int_t i=lowbin; i < hibin && i < ncol*ncol; ++i)
{
c0->cd(i+1-lowbin);
- gPad->SetGridx();
- gPad->SetGridy();
gPad->SetLogy();
- hData[i]->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
+ hData[i]->SetTitle(Form("%0.2f < #it{p}_{T} < %0.2f GeV/c", hPt->GetBinLowEdge(i), hPt->GetBinLowEdge(i)+hPt->GetBinWidth(i)));
hData[i]->SetMinimum(0.2);
hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
return c0;
}
-TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, Double_t binwidth, Double_t dcaxyMin, Double_t dcaxyMax)
+TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const TString& name, const TString& title, const TH1D* hPt, Double_t dcaxyMin, Double_t dcaxyMax)
{
//
// Draw a goodness of fit for each pt bin consisting of data/fit
{
c0->cd(i+1-lowbin);
- gPad->SetGrid();
-
TH1D* hDiv = (TH1D*)hData[i]->Clone("_Div_Data_Fit_");
hDiv->Sumw2();
hDiv->Divide(hDataFit[i]);
- hDiv->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", (i-1)*binwidth, i*binwidth));
+ hDiv->SetTitle(Form("%0.2f < p_{T} < %0.2f GeV/c", hPt->GetBinLowEdge(i), hPt->GetBinLowEdge(i)+hPt->GetBinWidth(i)));
hDiv->SetAxisRange(dcaxyMin, dcaxyMax, "X");
hDiv->SetAxisRange(0,3,"Y");
hDiv->SetMarkerStyle(kFullCircle);
// Draw differential yields
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
-#include <Riostream.h>
-#include <TROOT.h>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TStyle.h>
#include <TFile.h>
#include <TString.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
-#include <TF1.h>
+#endif
#include "B2.h"
//
// Draw differential yields
//
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
gStyle->SetOptLogy(1);
TFile* finput = new TFile(inputFile.Data());
if (finput->IsZombie()) exit(1);
TCanvas* c0 = new TCanvas(Form("%s.Spectra", particle.Data()), Form("Spectra for %ss", particle.Data()));
- c0->Divide(2,2);
+ c0->Divide(2,1);
- // differential yields
+ // differential yield
c0->cd(1);
- TGraphErrors* grDYieldPt = (TGraphErrors*)FindObj(finput, tag, particle + "_DiffYield_Pt");
- TGraphErrors* grSysErrDYieldPt = (TGraphErrors*)FindObj(finput, tag, particle + "_SysErr_DiffYield_Pt");
-
- grSysErrDYieldPt->SetLineColor(kRed);
- grSysErrDYieldPt->SetFillStyle(0);
- grSysErrDYieldPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
- grSysErrDYieldPt->GetYaxis()->SetTitle("#frac{1}{N_{inel}} #frac{d^{2}N}{dp_{T}dy} (GeV^{-1}c^{2})");
- grSysErrDYieldPt->Draw("A5");
+ TGraphErrors* grDYieldPt = FindObj<TGraphErrors>(finput, tag, particle + "_DiffYield_Pt");
grDYieldPt->SetMarkerStyle(kFullCircle);
grDYieldPt->SetMarkerColor(kRed);
grDYieldPt->SetLineColor(kRed);
- grDYieldPt->Draw("P");
+ grDYieldPt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ grDYieldPt->GetYaxis()->SetTitle("1/#it{N}_{ev} d^{2}#it{N}/(d#it{p}_{T}d#it{y}) (GeV^{-1}#it{c}^{2})");
+ grDYieldPt->GetYaxis()->SetTitleOffset(1.30);
+ grDYieldPt->Draw("zAP");
- // invariant differential yields
+ // invariant differential yield
c0->cd(2);
- TGraphErrors* grInvDYieldPt = (TGraphErrors*)FindObj(finput, tag, particle + "_InvDiffYield_Pt");
- TGraphErrors* grSysErrInvDYieldPt = (TGraphErrors*)FindObj(finput, tag, particle + "_SysErr_InvDiffYield_Pt");
-
- grSysErrInvDYieldPt->SetLineColor(kRed);
- grSysErrInvDYieldPt->SetFillStyle(0);
- grSysErrInvDYieldPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
- grSysErrInvDYieldPt->GetYaxis()->SetTitle("#frac{1}{2#piN_{inel}} #frac{d^{2}N}{p_{T}dp_{T}dy} (GeV^{-2}c^{3})");
- grSysErrInvDYieldPt->Draw("A5");
+ TGraphErrors* grInvDYieldPt = FindObj<TGraphErrors>(finput, tag, particle + "_InvDiffYield_Pt");
grInvDYieldPt->SetMarkerStyle(kFullCircle);
grInvDYieldPt->SetMarkerColor(kRed);
grInvDYieldPt->SetLineColor(kRed);
- grInvDYieldPt->Draw("P");
-
- // invariant differential cross section
-
- c0->cd(3);
-
- TGraphErrors* grInvDXsectPt = (TGraphErrors*)FindObj(finput, tag, particle + "_InvDiffXSection_Pt");
- TGraphErrors* grSysErrInvDXsectPt = (TGraphErrors*)FindObj(finput, tag, particle + "_SysErr_InvDiffXSection_Pt");
-
- grSysErrInvDXsectPt->SetLineColor(kRed);
- grSysErrInvDXsectPt->SetFillStyle(0);
- grSysErrInvDXsectPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
- grSysErrInvDXsectPt->GetYaxis()->SetTitle("E #frac{d^{3}#sigma}{dp^{3}} (mb GeV^{-2}c^{3})");
- grSysErrInvDXsectPt->Draw("A5");
-
- grInvDXsectPt->SetMarkerStyle(kFullCircle);
- grInvDXsectPt->SetMarkerColor(kRed);
- grInvDXsectPt->SetLineColor(kRed);
- grInvDXsectPt->Draw("P");
+ grInvDYieldPt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+ grInvDYieldPt->GetYaxis()->SetTitle("1/(2#pi#it{N}_{ev}) d^{2}#it{N}/(#it{p}_{T}d#it{p}_{T}d#it{y}) (GeV^{-2}#it{c}^{3})");
+ grInvDYieldPt->GetYaxis()->SetTitleOffset(1.30);
+ grInvDYieldPt->Draw("zAP");
}
// call Config_XXX for each period\r
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>\r
\r
+#if !defined(__CINT__) || defined(__MAKECINT__)\r
#include <Riostream.h>\r
#include <TROOT.h>\r
#include <TSystem.h>\r
#include <TString.h>\r
#include <TFileMerger.h>\r
-\r
#include "AliLnDriver.h"\r
+#endif\r
+\r
#include "Config.h"\r
\r
-Int_t LHC10bcde(const TString& species = "Deuteron",\r
- const TString& inputDir = "~/alice/input",\r
- const TString& outputDir = "~/alice/output",\r
- const TString& outputTag = "lhc10bcde",\r
- const TString& multTag = "",\r
- const TString& multCorTag = "",\r
- Bool_t inel = 1, // for mult\r
- Bool_t drawOutput = 1, // for batch\r
- Int_t option = 2)\r
+Int_t LHC10bcde( const TString& species = "Deuteron"\r
+ , const TString& inputDir = "~/alice/input"\r
+ , const TString& outputDir = "~/alice/output"\r
+ , const TString& outputTag = "lhc10bcde"\r
+ , const TString& trkselTag = "-tpc3-nsd-moc"\r
+ , const TString& multTag = ""\r
+ , const TString& multCorTag = ""\r
+ , Double_t ymax = 0.5\r
+ , Bool_t inel = 0 // for mult\r
+ , Bool_t drawOutput = 1 // for batch\r
+ , Int_t option = 2)\r
{\r
//\r
// call Config_XXX for each period, merge the corrected pt and then get the results\r
const TString kPeriod[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
const TString kOutputTag[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
\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
+ Double_t ptmin = (species=="Proton") ? 0.4 : 0.7;\r
+ Double_t ptjoint = (species=="Proton") ? 1.0 : 1.0;\r
+ Double_t ptmax = (species=="Proton") ? 2.0 : 3.0;\r
+ Double_t ptpid = (species=="Proton") ? 3.5 : 1.6;\r
\r
using namespace std;\r
\r
+ "\"" + outputDir + "\","\r
+ "\"" + kPeriod[i] + "\","\r
+ "\"" + kOutputTag[i] + "\","\r
+ + "\"" + trkselTag + "\","\r
+ "\"" + multTag + "\","\r
+ "\"" + multCorTag;\r
\r
{\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
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %f, %d, 0)", species.Data(), arg.Data(), ymax, 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
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %f, %d, 0)", species.Data(), arg.Data(), ymax, 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
+ gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\", %f, %d, 0, \"%s\", %f, %f, %f, %f)", arg.Data(), ymax, inel, species.Data(), ptmin, ptjoint, ptmax, ptpid));\r
break;\r
}\r
\r
m.AddFile(ptfile,0);\r
}\r
\r
- TString outputPt = outputDir + "/" + species + "-" + outputTag + "-Pt.root";\r
- TString outputRatio = outputDir + "/" + species + "-" + outputTag + "-Ratio.root";\r
- TString outputSpectra = outputDir + "/" + species + "-" + outputTag + "-Spectra.root";\r
+ TString outputPt = outputDir + "/" + species + "-" + outputTag + multTag + "-Pt.root";\r
+ TString outputRatio = outputDir + "/" + species + "-" + outputTag + multTag + "-Ratio.root";\r
+ TString outputSpectra = outputDir + "/" + species + "-" + outputTag + multTag + "-Spectra.root";\r
\r
// pt\r
\r
driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);\r
driver.SetOutputTag(outputTag);\r
\r
+ driver.SetRapidityInterval(-ymax,ymax);\r
+ driver.SetExtrapolateToINEL(inel);\r
+ \r
driver.SetMakeCorrections(0);\r
driver.SetMakePt(0);\r
driver.SetMakeRatio(1);\r
m2.AddFile(outputRatio.Data(),0);\r
m3.AddFile(outputSpectra.Data(),0);\r
\r
- TString allRatios = outputDir + "/" + species + "-lhc10bcde1-Ratio.root";\r
- TString allSpectra = outputDir + "/" + species + "-lhc10bcde1-Spectra.root";\r
+ TString allRatios = outputDir + "/" + species + "-" + outputTag + "-2" + multTag + "-Ratio.root";\r
+ TString allSpectra = outputDir + "/" + species + "-" + outputTag + "-2" + multTag + "-Spectra.root";\r
\r
m2.OutputFile(allRatios.Data());\r
m3.OutputFile(allSpectra.Data());\r
\r
if(!drawOutput) return 0;\r
\r
- gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"Anti%s%s_Ratio_Pt\",\"lhc10bcde\",0,4.5, 0., 1.8, \"p_{T} (GeV/c)\", \"neg/pos\", 1, \"cRatio\",\"Particle ratio\")", allRatios.Data(), species.Data(), species.Data()));\r
+ gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"Anti%s%s_Ratio_Pt\",\"%s\",0,4.5, 0., 1.8, \"p_{T} (GeV/c)\", \"neg/pos\", 2, \"cRatio\",\"Particle ratio\")", allRatios.Data(), species.Data(), species.Data(),outputTag.Data()));\r
+ \r
+ Double_t minYield = (species=="Proton") ? 1.1e-6 : 2.e-8;\r
+ Double_t maxYield = (species=="Proton") ? 4.e-1 : 9.e-4;\r
\r
- DrawOutputSpectraMult(allSpectra, species, ymin, ymax, 1, "lhc10bcde");\r
+ DrawOutputSpectraMult(allSpectra, species, minYield, maxYield, 2, outputTag);\r
\r
return 0;\r
}\r
// LHC10x multiplicity config for protons and deuterons
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TSystem.h>
#include <TROOT.h>
#include <TString.h>
#include <TFileMerger.h>
-
#include "AliLnDriver.h"
+#endif
+
#include "Config.h"
-Int_t LHC10xMult( const TString& species = "Proton"
+Int_t LHC10xMult( const TString& species = "Deuteron"
, const TString& inputDir = "~/alice/input"
, const TString& outputDir = "~/alice/output"
, const TString& period = "lhc10d"
, const TString& otag = "lhc10d"
- , Int_t option = 2)
+ , const TString& trkselTag = "-tpc3-nsd-moc"
+ , Double_t ymax = 0.5
+ , Int_t option = 2)
{
//
// lhc10x multiplicity config
const Bool_t kINEL[kNmult] = { 0 };
- 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 : 7.e-4;
-
+ Double_t ptmin = (species=="Proton") ? 0.4 : 0.7;
+ Double_t ptjoint = (species=="Proton") ? 1.0 : 1.0;
+ Double_t ptmax = (species=="Proton") ? 2.0 : 2.5;
+ Double_t ptpid = (species=="Proton") ? 0.4 : 1.3;
if( (option<0) || (option>2) || ((species != "Proton") && (species != "Deuteron")))
{
for(Int_t i=0; i<kNmult; ++i)
{
+ // limit the ptmax value for deuterons
+ if(species == "Deuteron") ptmax = kDtrPtMax[i];
+
cout << endl;
cout << "Multiplicity class : " << kMultTag[i] << endl;
cout << "Period : " << period << endl;
- TString outputTag = otag + "-" + kMultTag[i];
+ TString outputTag = otag + kMultTag[i];
- TString arg = inputDir + "\","
- + "\"" + outputDir + "\","
- + "\"" + period + "\","
- + "\"" + outputTag + "\","
- + "\"" + "-" + kMultTag[i] + "\"," // data
- + "\"" + ""; // same simulations for all mult
+ TString arg = inputDir + "\","
+ + "\"" + outputDir + "\","
+ + "\"" + period + "\","
+ + "\"" + outputTag + "\","
+ + "\"" + trkselTag + "\","
+ + "\"" + kMultTag[i] + "\"," // data
+ + "\"" + ""; // same simulations for all mult
switch(option)
{
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]));
+ gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %f, %d, 0, %f, %f)", species.Data(), arg.Data(), ymax, kINEL[i], ptmin, ptmax));
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]));
+ gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %f, %d, 0, %f, %f, %f)", species.Data(), arg.Data(), ymax, kINEL[i], ptmin, ptmax, ptpid));
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));
+ gROOT->ProcessLine(Form(".x Config_TPCTOF_LHC10x.C+g(\"%s\", %f, %d, 0, \"%s\", %f, %f, %f, %f)", arg.Data(), ymax, kINEL[i], species.Data(), ptmin, ptjoint, ptmax, ptpid));
break;
}
// draw output
- gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"Anti%s%s_Ratio_Pt\",\"\",0,4.5, 0., 1.8, \"p_{T} (GeV/c)\", \"#bar{p}/p\", 0, \"cRatio\",\"Particle ratio\")", allRatios.Data(),species.Data(),species.Data()));
+ gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"Anti%s%s_Ratio_Pt\",\"\",0,4.5, 0., 1.8, \"#it{p}_{T} (GeV/c)\", \"#bar{p}/p\", 0, \"cRatio\",\"Particle ratio\")", allRatios.Data(),species.Data(),species.Data()));
+
+ Double_t minYield = (species=="Proton") ? 1.1e-6 : 1.1e-8;
+ Double_t maxYield = (species=="Proton") ? 4.e-1 : 7.e-4;
- DrawOutputSpectraMult(allSpectra, species, ymin, ymax);
+ DrawOutputSpectraMult(allSpectra, species, minYield, maxYield);
return 0;
}
// d/p ratio
// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
#include <TROOT.h>
#include <TFile.h>
#include <TString.h>
#include <TFitResult.h>
#include <TFitResultPtr.h>
#include <TVirtualFitter.h>
+#include <TPaveText.h>
+#include <TAxis.h>
+#endif
#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-nsd-Mult-Spectra.root"
- , const TString& dSpectra = "~/alice/output/Deuteron-lhc10bcde-nsd-Mult-Spectra.root"
- , const TString& ptag = "lhc10d-nsd"
- , const TString& dtag = "lhc10bcde-nsd"
- , 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-Tsallis-Mult.root"
+void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root"
+ , const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-hilow-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 qTsallis = 0
+ , const TString& outputfile = "~/alice/output/Particle-Ratios-lhc10d-Tsallis.root"
, const TString& otag = "Tsallis")
{
+
//
// d/p ratio as a function of X
//
const TString kParticle[kNspec][kNpart] = { {"Proton", "AntiProton"}, { "Deuteron", "AntiDeuteron"} };
+ const Double_t xmin[kNspec] = {0., 0.};
+ const Double_t xmax[kNspec] = {100, 100};
+
+ const Double_t xminf[kNspec] = {0.4, 0.8};
+ const Double_t xmaxf[kNspec] = {2., 2.2};
+
TVirtualFitter::SetDefaultFitter("Minuit2");
TFile* fproton = new TFile(pSpectra.Data(), "read");
{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* grMeanPt[kNspec][kNpart] = {{new TGraphErrors(), new TGraphErrors()},
+ {new TGraphErrors(), new TGraphErrors()}};
+
+ using std::cout;
+ using std::endl;
+
TGraphErrors* grDYieldPt[nx][kNspec][kNpart];
TF1* fncTsallis[nx][kNspec][kNpart];
Double_t dNdy[kNspec][kNpart];
Double_t dNdyErr[kNspec][kNpart];
+ cout << endl << xtag[i] << endl << endl;
+
for(Int_t j=0; j<kNspec; ++j)
{
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");
+ // data
- 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);
- }
+ grDYieldPt[i][j][k] = FindObj<TGraphErrors>(finput[j], tag[j] + xtag[i], kParticle[j][k] + "_StatSystErr_DiffYield_Pt");
+
+ // Tsallis distribution
+
+ fncTsallis[i][j][k] = (qTsallis) ? QTsallisDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_DiffYield_Pt",xmin[j],xmax[j])
+ : TsallisDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_DiffYield_Pt",xmin[j],xmax[j]);
+
+ // fit to data
TFitResultPtr r = grDYieldPt[i][j][k]->Fit(fncTsallis[i][j][k], "RENSQ");
+ Int_t status = r;
- if(tsallis)
+ Int_t n=0;
+ while(status != 0 && n<7)
{
- 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);
+ r = grDYieldPt[i][j][k]->Fit(fncTsallis[i][j][k], "RENSQ");
+ status = r;
+ ++n;
+ printf("\t*** %s, status: %d (%d)\n", kParticle[j][k].Data(), status,n);
}
- Int_t status = r;
- printf("status: %d\tdN/dy: %g +/- %g\n",status,dNdy[j][k],dNdyErr[j][k]);
+ // integral
+
+ Double_t integral = fncTsallis[i][j][k]->Integral(xmin[j], xmax[j]);
+ Double_t intErr = fncTsallis[i][j][k]->IntegralError(xmin[j], xmax[j], r->GetParams(), r->GetCovarianceMatrix().GetMatrixArray());
+
+ // mean pt (dy cancels out)
+
+ TF1* fncPtTsallis = (qTsallis) ? PtQTsallisDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_MeanPt",xmin[j],xmax[j])
+ : PtTsallisDYield(GetMass(kParticle[j][k]), kParticle[j][k] + "_Tsallis_MeanPt",xmin[j],xmax[j]) ;
+
+ Double_t par[] = { fncTsallis[i][j][k]->GetParameter(0)
+ , fncTsallis[i][j][k]->GetParameter(1)
+ , fncTsallis[i][j][k]->GetParameter(2)
+ };
+
+ Double_t parErr[] = { fncTsallis[i][j][k]->GetParError(0)
+ , fncTsallis[i][j][k]->GetParError(1)
+ , fncTsallis[i][j][k]->GetParError(2)
+ };
+
+ fncPtTsallis->SetParameters(par);
+ fncPtTsallis->SetParErrors(parErr);
+
+ Double_t intPtTsallis = fncPtTsallis->Integral(xmin[j], xmax[j]);
+ Double_t intPtTsallisErr = fncPtTsallis->IntegralError(xmin[j], xmax[j], r->GetParams(), r->GetCovarianceMatrix().GetMatrixArray());
+
+ Double_t meanPt, meanPtErr;
+ GetRatio(intPtTsallis, integral, intPtTsallisErr, intErr, meanPt, meanPtErr);
+
+ grMeanPt[j][k]->SetPoint(i, x[i], meanPt);
+ grMeanPt[j][k]->SetPointError(i, xerr[i], meanPtErr);
+
+ delete fncPtTsallis;
+
+ // extrapolation fraction
+
+ Double_t intFrac = fncTsallis[i][j][k]->Integral(xminf[j], xmaxf[j]);
+ Double_t intFracErr = fncTsallis[i][j][k]->IntegralError(xminf[j], xmaxf[j], r->GetParams(), r->GetCovarianceMatrix().GetMatrixArray());
+
+ Double_t extrFrac, extrFracErr;
+ GetRatio(intFrac, integral, intFracErr, intErr, extrFrac, extrFracErr);
+
+ // dN/dy
+
+ dNdy[j][k] = (qTsallis) ? fncTsallis[i][j][k]->GetParameter(0) : integral;
+ dNdyErr[j][k] = (qTsallis) ? fncTsallis[i][j][k]->GetParError(0) : intErr;
grdNdy[j][k]->SetPoint(i, x[i], dNdy[j][k]);
grdNdy[j][k]->SetPointError(i, xerr[i], dNdyErr[j][k]);
+
+ // debug
+ printf("\t*** %s, status: %d\tdN/dy: %g +/- %g\t extr: %g +/- %g\t<pt>: %g +/- %g\n", kParticle[j][k].Data(), status, dNdy[j][k], dNdyErr[j][k], 1.-extrFrac, extrFracErr, meanPt, meanPtErr);
}
}
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());
}
}
}
grP2[i][j]->Write();
grdNdy[i][j]->SetName(Form("%s_dNdy", kParticle[i][j].Data()));
- grChi2Ndf[i][j]->SetName(Form("%s_Chi2Ndf", kParticle[i][j].Data()));
-
grdNdy[i][j]->Write();
- grChi2Ndf[i][j]->Write();
+
+ grMeanPt[i][j]->SetName(Form("%s_MeanPt", kParticle[i][j].Data()));
+ grMeanPt[i][j]->Write();
}
}
// draw
- TStyle* st = new TStyle();
-
- st->SetPadTickX(1);
- st->SetPadTickY(1);
- st->SetPadGridX(1);
- st->SetPadGridY(1);
- st->SetCanvasColor(0);
- st->SetFrameBorderMode(0);
- st->SetFrameFillColor(0);
- st->SetLabelFont(62,"XYZ");
- st->SetTitleFont(62,"XYZ");
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetOptStat(0);
- st->cd();
- gROOT->ForceStyle();
-
- const Int_t kNCol = 4;
+ const Int_t kNCol = 3;
TCanvas* c0[kNspec];
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");
}
}
// spectra
- TCanvas* c2[kNspec];
+ TPaveText* label[nx][kNspec][kNpart];
+
+ TCanvas* c2[kNspec][kNpart];
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)
+ c2[j][k] = new TCanvas(Form("c2.%s",kParticle[j][k].Data()), Form("%s data",kParticle[j][k].Data()));
+ c2[j][k]->Divide(3,2);
+ for(Int_t i=0; i<nx && i<6; ++i)
{
gPad->SetLogy(0);
- c2[j]->cd(nx*k+i+1);
- Draw(grDYieldPt[i][j][k], kFullCircle, kBlue, "p_{T} (GeV/c)", "DYield");
+ c2[j][k]->cd(i+1);
+ Draw(grDYieldPt[i][j][k], kFullCircle, kBlue, "p_{T} (GeV/c)", "1/N_{ev} d^{2}N/dp_{T}dy (GeV^{-1})");
+ fncTsallis[i][j][k]->SetLineWidth(1);
fncTsallis[i][j][k]->Draw("same");
+
+ label[i][j][k] = new TPaveText(0.5658525,0.7471751,0.814296,0.835452,"brNDC");
+ label[i][j][k]->AddText(Form("%s = %.2f",xname.Data(),x[i]));
+ label[i][j][k]->SetBorderSize(0);
+ label[i][j][k]->SetFillColor(0);
+ label[i][j][k]->SetTextSize(0.06);
+ label[i][j][k]->Draw();
}
}
}
gr->SetLineColor(color);
gr->GetXaxis()->SetTitle(xtitle.Data());
gr->GetYaxis()->SetTitle(ytitle.Data());
- gr->Draw("AP");
+ gr->GetYaxis()->SetTitleOffset(1.3);
+ gr->Draw("zAP");
}
gROOT->SetMacroPath(Form("%s:$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/", gROOT->GetMacroPath()));
- // Install RooUnfold and uncomment and edit AliLnPt.{h,cxx}
- //gSystem->Load("$ALICE_ROOT/../RooUnfold/libRooUnfold.so");
-
gSystem->AddIncludePath("-I\"$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/\"");
- gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx+g");
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx+g");
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx+g");
gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx+g");
+ gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C+g");
gROOT->ProcessLine(Form(".x %s", config.Data()));
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff