From aa54def0e65af55682297ba5cb286d627b2582f4 Mon Sep 17 00:00:00 2001 From: eserradi Date: Tue, 10 Sep 2013 13:30:46 +0000 Subject: [PATCH] cumulative changes for root scripts and code cleanup --- PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx | 106 ++++-- PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h | 6 +- PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx | 8 +- PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h | 3 - PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx | 153 ++++---- PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h | 52 +-- PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx | 14 +- .../Nuclei/B2/AliLnGaussianExpTail.cxx | 57 +++ .../SPECTRA/Nuclei/B2/AliLnGaussianExpTail.h | 42 +++ PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx | 59 +++ PWGLF/SPECTRA/Nuclei/B2/AliLnM2.h | 44 +++ PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx | 344 ++++++++++-------- PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h | 61 ++-- PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx | 2 +- PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx | 246 +++---------- PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h | 13 +- PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx | 85 +---- PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h | 18 +- PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx | 87 ----- PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h | 46 --- PWGLF/SPECTRA/Nuclei/B2/B2.h | 104 ++++-- PWGLF/SPECTRA/Nuclei/B2/Config.h | 107 ++---- PWGLF/SPECTRA/Nuclei/B2/macros/B2.C | 3 +- PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C | 81 +++-- .../B2/macros/Config_Deuteron_TOF_LHC10x.C | 119 +++--- .../B2/macros/Config_Deuteron_TPC_LHC10x.C | 83 +++-- .../B2/macros/Config_Proton_TOF_LHC10x.C | 89 ++--- .../B2/macros/Config_Proton_TPC_LHC10x.C | 83 ++--- .../Nuclei/B2/macros/Config_TPCTOF_LHC10x.C | 69 ++-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C | 53 +-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C | 50 +-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C | 43 +-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C | 90 +++-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C | 14 +- PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C | 76 ++-- PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C | 61 +--- PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C | 61 ++-- PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C | 52 +-- PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C | 184 ++++++---- PWGLF/SPECTRA/Nuclei/B2/macros/run.C | 7 +- 40 files changed, 1434 insertions(+), 1441 deletions(-) create mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.cxx create mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.h create mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx create mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnM2.h delete mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx delete mode 100644 PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx index 15a9d61c8de..dc3468ef5b2 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx @@ -96,11 +96,8 @@ Int_t AliLnB2::Run() // 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(finput1, fProtonTag, prefix + "Proton_InvDiffYield_Pt"); + TGraphErrors* grNucInvDYieldPt = FindObj(finputA, fNucleusTag, prefix + fNucleusName + "_InvDiffYield_Pt"); TFile* foutput = new TFile(fOutputFilename.Data(),"recreate"); @@ -111,25 +108,42 @@ Int_t AliLnB2::Run() 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(finput1->Get(grPrtName.Data())); + TGraphErrors* grSystErrNucInvDYieldPt = dynamic_cast(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; @@ -138,7 +152,7 @@ Int_t AliLnB2::Run() 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 @@ -146,45 +160,64 @@ TGraphErrors* AliLnB2::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGrap 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; iGetN(); ++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 { // @@ -204,7 +237,7 @@ 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 @@ -226,7 +259,8 @@ TGraphErrors* AliLnB2::Rside2Rlong(const TGraphErrors* grB2, const TString& name 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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h index dee42d2d594..ab529912f4d 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h @@ -19,10 +19,10 @@ class AliLnB2: public TObject 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(); @@ -34,6 +34,8 @@ class AliLnB2: public TObject 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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx index db2b7a6f2e3..6bcb288f4b0 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx @@ -22,7 +22,6 @@ #include #include -#include "AliLnUnfolding.h" #include "AliLnSecondaries.h" #include "AliLnEfficiency.h" #include "AliLnCorr.h" @@ -37,7 +36,6 @@ AliLnCorr::AliLnCorr(const TString& particle, const TString& dataFilename, const // 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 } @@ -47,7 +45,6 @@ AliLnCorr::~AliLnCorr() // // destructor // - delete fUnfolding; delete fSecondaries; delete fEfficiency; } @@ -57,19 +54,16 @@ Int_t AliLnCorr::Exec() // // 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); @@ -78,7 +72,7 @@ Int_t AliLnCorr::Exec() 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())); diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h index 3a1161c0126..a0a8c3b4de7 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h @@ -10,7 +10,6 @@ #include #include -class AliLnUnfolding; class AliLnSecondaries; class AliLnEfficiency; @@ -23,7 +22,6 @@ class AliLnCorr: public TObject Int_t Exec(); - AliLnUnfolding* GetLnUnfolding() { return fUnfolding; } AliLnSecondaries* GetLnSecondaries() { return fSecondaries; } AliLnEfficiency* GetLnEfficiency() { return fEfficiency; } @@ -36,7 +34,6 @@ class AliLnCorr: public TObject TString fOutputFilename; // output filename - AliLnUnfolding* fUnfolding; // unfolding correction AliLnSecondaries* fSecondaries; // secondaries correction AliLnEfficiency* fEfficiency; // eficiency correction diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx index 309ad0ae434..62fafc20f8e 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx @@ -34,6 +34,7 @@ #include "AliLnDriver.h" #include "B2.h" + ClassImp(AliLnDriver) AliLnDriver::AliLnDriver() @@ -48,14 +49,12 @@ 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) @@ -63,16 +62,14 @@ AliLnDriver::AliLnDriver() , 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() @@ -87,6 +84,9 @@ AliLnDriver::AliLnDriver() , fSameFdwn(0) , fMCtoINEL(0) , fAddFakeTracks(1) +, fPidProc(AliLnPt::kMassSquared) +, fPidEff(1) +, fDebugLevel(0) { // // constructor @@ -94,7 +94,6 @@ AliLnDriver::AliLnDriver() for(Int_t i=0; i<3; ++i) { fTrigEff[i] = 1; - fXsec[i] = 1; } gErrorIgnoreLevel = kWarning; @@ -135,28 +134,6 @@ void AliLnDriver::SetOutputFilenames(const TString& pt, const TString& ratio, co 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]) { // @@ -165,26 +142,80 @@ 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; } @@ -206,8 +237,8 @@ void AliLnDriver::MakePtCorr() const 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); @@ -255,14 +286,13 @@ void AliLnDriver::MakePt() const 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); @@ -270,6 +300,8 @@ void AliLnDriver::MakePt() const lnpt.SetFitFractionCorr(fFitFrac); lnpt.SetFeedDownCorr(fFdwnCorr); lnpt.SetSameFeedDownCorr(fSameFdwn); + lnpt.SetPidEfficiency(fPidEff); + lnpt.SetDebugLevel(fDebugLevel); lnpt.Exec(); @@ -310,19 +342,16 @@ void AliLnDriver::MakeSpectra() const // 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(); diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h index 6e54dbd44ac..6a14cb4d084 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h @@ -25,8 +25,6 @@ class AliLnDriver: public TObject void MakeRatio() const; void MakeSpectra() const; - void PrintFilenames() const; - Int_t Run() const; TString GetSpecies() const { return fSpecies; } @@ -46,7 +44,6 @@ class AliLnDriver: public TObject 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; } @@ -58,14 +55,13 @@ class AliLnDriver: public TObject 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; } @@ -78,9 +74,13 @@ class AliLnDriver: public TObject 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: @@ -104,15 +104,12 @@ class AliLnDriver: public TObject 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 @@ -120,18 +117,15 @@ class AliLnDriver: public TObject 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 @@ -151,6 +145,12 @@ class AliLnDriver: public TObject 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) }; diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx index c7d931e0467..f5b80298410 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx @@ -67,12 +67,12 @@ Int_t AliLnEfficiency::Exec() 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(fsimu, fParticle + "_Gen_PhS_Prim_Pt"); + TH1D* hGenTrigPt = FindObj(fsimu, fParticle + "_Gen_Trig_Prim_Pt"); + TH1D* hGenVtxPt = FindObj(fsimu, fParticle + "_Gen_Vtx_Prim_Pt"); + TH1D* hGenAccPt = FindObj(fsimu, fParticle + "_Gen_Acc_Prim_Pt"); + TH1D* hPrimRecPt = FindObj(fsimu, fParticle + "_Sim_Prim_Pt"); + TH1D* hFakePrimRecPt = FindObj(fsimu, fParticle + "_Sim_Fake_Prim_Pt"); if(fAddFakeTracks) hPrimRecPt->Add(hFakePrimRecPt); @@ -99,7 +99,7 @@ Int_t AliLnEfficiency::Exec() if(!fParticle.Contains("Anti")) { - TH1D* hStats = (TH1D*)FindObj(fsimu, fParticle + "_Stats"); + TH1D* hStats = FindObj(fsimu, fParticle + "_Stats"); hStats->Write(); } diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.cxx new file mode 100644 index 00000000000..d496b215c0c --- /dev/null +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.cxx @@ -0,0 +1,57 @@ +/***************************************************************************** + * 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 +#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); + } + + + diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.h new file mode 100644 index 00000000000..f718b9eb4a8 --- /dev/null +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnGaussianExpTail.h @@ -0,0 +1,42 @@ +/***************************************************************************** + * 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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx new file mode 100644 index 00000000000..0287fe2bb0b --- /dev/null +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.cxx @@ -0,0 +1,59 @@ +/***************************************************************************** + * 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 +#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)); +} diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.h new file mode 100644 index 00000000000..dd62b442d92 --- /dev/null +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnM2.h @@ -0,0 +1,44 @@ +/***************************************************************************** + * 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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx index e5f4e96e0fe..0e13808198b 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx @@ -31,15 +31,11 @@ #include #include #include +#include #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) @@ -51,37 +47,38 @@ AliLnPt::AliLnPt(const TString& particle, Double_t trigEff, const TString& input , 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() @@ -96,8 +93,6 @@ Int_t AliLnPt::Exec() // // Get the true pt distribution // - using namespace std; - TFile* finput = new TFile(fInputFilename.Data(), "read"); if (finput->IsZombie()) exit(1); @@ -110,15 +105,15 @@ Int_t AliLnPt::Exec() TString species = fParticle; species.ReplaceAll("Anti",""); - TH1D* hStats = (TH1D*)FindObj(finput, species + "_Stats"); + TH1D* hStats = FindObj(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( (FindObj(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); @@ -128,127 +123,119 @@ Int_t AliLnPt::Exec() 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(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(finput, fParticle + "_PID_DTime_Pt"); + else if(fPidProc == kMassSquared) hPidPt = FindObj(finput, fParticle + "_PID_M2_Pt"); + else if(fPidProc == kMassDifference) hPidPt = FindObj(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(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt"); + TH1D* hFracFdwnPt = FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt"); + TF1* fncFracMatPt = FindObj(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt"); + TF1* fncFracFdwnPt = FindObj(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(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt"); + fncFracFdwnPt = FindObj(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(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt"); + TH1D* hEffAccTrkPt = FindObj(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(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(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(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(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(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(hCorrPt->Clone(Form("%s_Pt", fParticle.Data()))); hPt->SetTitle(fParticle.Data()); hPt->Reset(); hPt->Sumw2(); - for(Int_t i=fLowPtBin; iGetXaxis()->FindFixBin(fPtMin); + Int_t hibin = hPt->GetXaxis()->FindFixBin(fPtMax); + + for(Int_t i=lowbin; iSetBinContent(i, hCorPt->GetBinContent(i)); - hPt->SetBinError(i, hCorPt->GetBinError(i)); + hPt->SetBinContent(i, hCorrPt->GetBinContent(i)); + hPt->SetBinError(i, hCorrPt->GetBinError(i)); } foutput->cd(); @@ -256,12 +243,12 @@ Int_t AliLnPt::Exec() hPt->Write(); delete hPt; - delete hCorPt; + delete hCorrPt; delete fdebug; delete fcorr; - // --------- stats ------------------- + // stats if(fMakeStats) { @@ -332,8 +319,7 @@ Double_t AliLnPt::GetVertexCorrection(const TH1D* hData, const TH1D* hMC) const 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; @@ -365,130 +351,168 @@ Double_t AliLnPt::GetM2Width(Double_t pt, Double_t m) const 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 - (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; iProjectionY(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(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(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; iProjectionY(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 { // @@ -514,7 +538,7 @@ TH1D* AliLnPt::Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* 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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h index f19c80e4b2c..d6793335448 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h @@ -10,10 +10,6 @@ #include #include -#ifndef HAVE_ROOUNFOLD -//#define HAVE_ROOUNFOLD -#endif - class TString; class TH1D; class TH2D; @@ -27,33 +23,29 @@ class AliLnPt: public TObject 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; } @@ -66,15 +58,22 @@ class AliLnPt: public TObject 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: @@ -90,28 +89,22 @@ class AliLnPt: public TObject 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 @@ -122,6 +115,12 @@ class AliLnPt: public TObject 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) }; diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx index 014c0fc9c35..f8233842fba 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx @@ -65,7 +65,7 @@ Int_t AliLnRatio::Exec() for(Int_t i=0; i<2; ++i) { - hPt[i] = (TH1D*)FindObj(finput, kPrefix[i] + fSpecies + "_Pt"); + hPt[i] = FindObj(finput, kPrefix[i] + fSpecies + "_Pt"); } TH1D* hRatioPt = Divide(hPt[1],hPt[0],Form("Anti%s%s_Ratio_Pt",fSpecies.Data(),fSpecies.Data())); diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx index c0dce9689e2..3ec621b7a15 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx @@ -25,12 +25,6 @@ #include #include #include -#include -#include -#include -#include -#include -#include #include #include "AliLnSecondaries.h" @@ -45,8 +39,8 @@ AliLnSecondaries::AliLnSecondaries(const TString& particle, const TString& dataF , fSimuFilename(simuFilename) , fOutputFilename(outputFilename) , fOutputTag(otag) -, fLowPtBin(3) -, fHiPtBin(15) +, fPtMin(0.5) +, fPtMax(3.0) , fNbin(1) , fMinDCAxy(-1.5) , fMaxDCAxy(1.5) @@ -61,10 +55,6 @@ AliLnSecondaries::AliLnSecondaries(const TString& particle, const TString& dataF // constructor // TH1::SetDefaultSumw2(); // switch on histogram errors - - // disable verbose in RooFit - RooMsgService::instance().setSilentMode(1); - RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL); } AliLnSecondaries::~AliLnSecondaries() @@ -95,7 +85,7 @@ Int_t AliLnSecondaries::Exec() // --------- ideal values: primaries + no secondaries -------- - TH1D* hPidPt = (TH1D*)FindObj(fdata, fParticle + "_PID_Pt"); + TH1D* hPidPt = FindObj(fdata, fParticle + "_PID_Pt"); hPidPt->Write(); const char* contrib[] = { "prim", "mat", "fdwn" }; @@ -112,10 +102,10 @@ Int_t AliLnSecondaries::Exec() 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(fsimu, fParticle + "_Sim_Pt"); + TH1D* hPrim = FindObj(fsimu, fParticle + "_Sim_Prim_Pt"); + TH1D* hMat = FindObj(fsimu, fParticle + "_Sim_Mat_Pt"); + TH1D* hFdwn = FindObj(fsimu, fParticle + "_Sim_Fdwn_Pt"); delete hFracPt[0]; delete hFracPt[1]; @@ -131,12 +121,12 @@ Int_t AliLnSecondaries::Exec() 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(fdata, fParticle + "_PID_DCAxy_Pt"); + TH2D* hDataMCDCAxyPt = FindObj(fsimu, fParticle + "_PID_DCAxy_Pt"); + TH2D* hPrimDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt"); + TH2D* hFdwnDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt"); + TH2D* hFakePrimDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt"); + TH2D* hFakeFdwnDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"); if(fAddFakeTracks) { @@ -159,14 +149,14 @@ Int_t AliLnSecondaries::Exec() 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(fdata, fParticle + "_PID_DCAxy_Pt"); + TH2D* hDataMCDCAxyPt = FindObj(fsimu, fParticle + "_PID_DCAxy_Pt"); + TH2D* hPrimDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt"); + TH2D* hMatDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt"); + TH2D* hFdwnDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt"); + TH2D* hFakePrimDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt"); + TH2D* hFakeMatDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fake_Mat_DCAxy_Pt"); + TH2D* hFakeFdwnDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt"); if(fAddFakeTracks) { @@ -191,19 +181,19 @@ Int_t AliLnSecondaries::Exec() 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(fdata, fParticle + "_PID_DCAxy_Pt"); + TH2D* hDataMCDCAxyPt = FindObj(fsimu, fParticle + "_PID_DCAxy_Pt"); + TH2D* hPrimDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt"); + TH2D* hFakePrimDCAxyPt = FindObj(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(fdata, Form("Anti%s_PID_DCAxy_Pt",fParticle.Data())); } - TH2D* hMatDCAxyPt = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt"); + TH2D* hMatDCAxyPt = FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt"); if(fMatDCAxyMod == kFlatDCAxy) { hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt"); @@ -251,7 +241,7 @@ Int_t AliLnSecondaries::Exec() } else { - hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.4); + hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.6); } TF1* fncFdwn = this->GetFdwnFraction(fParticle + "_Frac_Fdwn_Fit_Pt"); @@ -311,7 +301,10 @@ void AliLnSecondaries::GetFraction(TH1D* hPrimPt, TH1D* hSecPt, const TH2D* hDCA // TString nosec = (sec == "Fdwn") ? "Mat" : "Fdwn"; - for(Int_t i=fLowPtBin; iGetXaxis()->FindFixBin(fPtMin); + Int_t hibin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax); + + for(Int_t i=lowbin; iProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i); TH1D* hMCDCAxy = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i); @@ -319,17 +312,10 @@ void AliLnSecondaries::GetFraction(TH1D* hPrimPt, TH1D* hSecPt, const TH2D* hDCA 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}; @@ -370,7 +356,10 @@ void AliLnSecondaries::GetFraction(TH1D* hFracPt[3], const TH2D* hDCAxyPt, const // slice the DCA distribution and get the fractions for each pt bin // (3 contributions) // - for(Int_t i=fLowPtBin; iGetXaxis()->FindFixBin(fPtMin); + Int_t hibin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax); + + for(Int_t i=lowbin; iProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i); @@ -383,14 +372,7 @@ void AliLnSecondaries::GetFraction(TH1D* hFracPt[3], const TH2D* hDCAxyPt, const 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) @@ -432,6 +414,8 @@ Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hDa fit->Constrain(2,0.0001,1.); Int_t status = fit->Fit(); + status = fit->Fit(); + status = fit->Fit(); if (status == 0) // get fractions { @@ -466,6 +450,8 @@ Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hDa //fit->SetRangeX(1,50); Int_t status = fit->Fit(); + status = fit->Fit(); + status = fit->Fit(); if (status == 0) // get fractions { @@ -525,150 +511,18 @@ TH1D* AliLnSecondaries::ZeroClone(const TH1D* h, const TString& name) const // // clone histogram and reset to zero // - TH1D* clone = (TH1D*)h->Clone(name.Data()); + TH1D* clone = dynamic_cast(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 { // diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h index eb1f8ed26c3..ce1a45d4cf1 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h @@ -34,7 +34,7 @@ class AliLnSecondaries: public TObject 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; } @@ -46,8 +46,8 @@ class AliLnSecondaries: public TObject 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: @@ -61,9 +61,6 @@ class AliLnSecondaries: public TObject 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; @@ -81,8 +78,8 @@ class AliLnSecondaries: public TObject 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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx index 54bfaf8d246..eb7fe1672ac 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx @@ -13,7 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -// invariant differential yields and cross sections +// (invariant) differential yield // author: Eulogio Serradilla #include @@ -29,7 +29,7 @@ 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) @@ -39,13 +39,10 @@ AliLnSpectra::AliLnSpectra(const TString& particle, const TString& ptFilename, c , 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() @@ -58,7 +55,7 @@ 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; @@ -67,8 +64,8 @@ Int_t AliLnSpectra::Exec() // pt and number of events - TH1D* hPt = (TH1D*)FindObj(finput, fParticle + "_Pt"); - TH1D* hStats = (TH1D*)FindObj(finput, fParticle + "_Stats"); + TH1D* hPt = FindObj(finput, fParticle + "_Pt"); + TH1D* hStats = FindObj(finput, fParticle + "_Stats"); Double_t nEvent = hStats->Integral(3,3); if(fINEL) nEvent = hStats->Integral(4,4); @@ -85,32 +82,15 @@ Int_t AliLnSpectra::Exec() 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; @@ -144,8 +124,8 @@ TGraphErrors* AliLnSpectra::GetDiffYieldPt(const TH1D* hPt, Int_t nEvent, const 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; @@ -190,54 +170,3 @@ TGraphErrors* AliLnSpectra::GetInvDiffYieldPt(const TH1D* hPt, Int_t nEvent, con 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; -} diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h index d320548bee6..d9edcbd4d79 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h +++ b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h @@ -4,40 +4,32 @@ /* 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 #include 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); @@ -57,10 +49,6 @@ class AliLnSpectra: public TObject 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) }; diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx deleted file mode 100644 index 08317ef53fc..00000000000 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx +++ /dev/null @@ -1,87 +0,0 @@ -/************************************************************************** - * 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 - -//#include -#include -#include -#include -#include -#include -#include - -#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; -} diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h deleted file mode 100644 index 96bc032d296..00000000000 --- a/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h +++ /dev/null @@ -1,46 +0,0 @@ -#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 - -#include -#include - -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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/B2.h b/PWGLF/SPECTRA/Nuclei/B2/B2.h index bfc0b184f05..c11d8169b00 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/B2.h +++ b/PWGLF/SPECTRA/Nuclei/B2/B2.h @@ -7,74 +7,78 @@ // some common functions // author: Eulogio Serradilla -#include +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include #include #include +#endif -TObject* FindObj(TFile* f, const TString& name) +template +inline T* FindObj(TFile* f, const TString& name) { // // check if the object exists // - TObject* obj = f->Get(name.Data()); + T* obj = dynamic_cast(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 +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(f,name); + return FindObj(f, dir + "/" + name + ";1"); } -TObject* FindObj(const TList* l, const TString& name) +template +inline T* FindObj(const TList* l, const TString& name) { // // check if the object exists // - TObject* obj = l->FindObject(name.Data()); + T* obj = dynamic_cast(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(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 @@ -97,7 +101,7 @@ Double_t GetMass(const TString& name) 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 @@ -109,36 +113,51 @@ TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0, Double_t xmax=10 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) // @@ -147,19 +166,17 @@ TF1* TsallisPareto(Double_t m0, const TString& name, Double_t xmin=0, Double_t x 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); @@ -173,4 +190,21 @@ TF1* TsallisParetoDYield(Double_t m0, const TString& name, Double_t xmin=0, Doub 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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/Config.h b/PWGLF/SPECTRA/Nuclei/B2/Config.h index fcb9b104db0..9cd3368cb08 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/Config.h +++ b/PWGLF/SPECTRA/Nuclei/B2/Config.h @@ -4,23 +4,25 @@ /* 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 +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include -#include #include +#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)"; }; @@ -30,21 +32,24 @@ namespace B2mult // // 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"; }; @@ -74,43 +79,10 @@ TString GetCollSystem(const TString& period) 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); @@ -191,8 +163,10 @@ TString GetSimuFixPeriod(const TString& 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"; @@ -204,6 +178,8 @@ TString MakeInputName(const TString& species, const TString& period, const TStri // // make input name for data // + if(trksel == "") return species + "-" + period; + if(trksel.BeginsWith("-")) return species + "-" + period + trksel; return species + "-" + period + "-" + trksel; } @@ -247,35 +223,6 @@ TString MakeOutputName(const TString& species, const TString& outputTag) 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="") { // @@ -284,7 +231,7 @@ void DrawOutputCorr(const TString& species, const TString& corr, const TString& 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 @@ -293,11 +240,11 @@ void DrawCorrDebug(const TString& sec, const TString& tag, const TString& specie 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 @@ -306,7 +253,7 @@ void DrawPtDebug(const TString& pt, const TString& tag, const TString& 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)); } } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C b/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C index c80ca627a52..a935e471237 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C @@ -16,12 +16,13 @@ // coalescence parameter // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include - #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" diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C index e467c3e1636..473e03646ef 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C @@ -16,6 +16,7 @@ // B2 as a function of multiplicity // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include @@ -23,8 +24,9 @@ #include #include #include - #include "AliLnB2.h" +#endif + #include "B2.h" #include "Config.h" @@ -37,18 +39,24 @@ Double_t GetCd(Double_t z) 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; @@ -69,7 +77,9 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M { 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])); @@ -80,7 +90,7 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M // merge B2 and B2bar - TString outputfile = otag + "-" + kMultTag[i] + "-B2.root"; + TString outputfile = otag + kMultTag[i] + "-B2.root"; m.OutputFile(outputfile.Data()); m.Merge(); @@ -94,7 +104,7 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M for(Int_t i=0; iExec(Form("rm -f %s%s-B2.root",otag.Data(),kMultTag[i].Data())); } // B2 as a function of multiplicity for each pt @@ -115,13 +125,18 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M TGraphErrors* grB2pt[kNpart][kNmult]; TGraphErrors* grR3pt[kNpart][kNmult]; + TGraphErrors* grSysB2pt[kNpart][kNmult]; + TGraphErrors* grSysR3pt[kNpart][kNmult]; for(Int_t i=0; i(finput, otag + kMultTag[j], Form("B2%s_Pt", kSuffix[i].Data())); + grR3pt[i][j] = FindObj(finput, otag + kMultTag[j], Form("R3%s_Pt", kSuffix[i].Data())); + + grSysB2pt[i][j] = FindObj(finput, otag + kMultTag[j], Form("SystErr_B2%s_Pt", kSuffix[i].Data())); + grSysR3pt[i][j] = FindObj(finput, otag + kMultTag[j], Form("SystErr_R3%s_Pt", kSuffix[i].Data())); } } @@ -138,7 +153,7 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M for(Int_t i=kNMinPt; imkdir(ptLabel[i].Data()); } @@ -148,24 +163,30 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M { 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; kGetPoint(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); @@ -174,22 +195,38 @@ Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd-M 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; kSetName(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 diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C index 2be6d91f20e..b088dd4e51d 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C @@ -16,62 +16,77 @@ // LHC10x config for deuterons and antideuterons // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include "AliLnDriver.h" +#endif + #include "Config.h" -Int_t Config_Deuteron_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 = 5, // 0.8 Gev/c - Int_t hiPtBin = 17, // 3.2 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_Deuteron_TOF_LHC10x( const TString& inputDir = "~/alice/input" + , const TString& outputDir = "~/alice/output" + , const TString& period = "lhc10bcde" + , const TString& outputTag = "lhc10bcde" + , const TString& trkselTag = "-tpc3-nsd-moc" + , 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.7 + , Double_t ptmax = 3.0 + , Double_t ptpid = 0.8 + , 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 deuterons and antideuterons // const TString kSpecies = "Deuteron"; - const TString kTrkSel = "its_tpc_tof_dca-tpc3"; - const TString kTrigName = "mbor"; + const TString kTrkSel = "its_tpc_tof_dca"; + const TString kTrigName = "mband"; const Bool_t kMCtoINEL = 1; - const Int_t kM2Bin[2] = {8,18}; - const Bool_t kPidM2 = 1; - const Bool_t kUnfolding = 0; - const Int_t kIter = 5; + const Bool_t kPid = 1; + const Int_t kPidProc = 0; // 0 m2, 1 dm2, 2 time + const Double_t kPidEff = 1.; 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 Bool_t kAntiNucTemplate = 0; const Int_t kNbin = 5; const Double_t kDCAxy[2] = {-0.2,0.2}; - const Double_t kM2Bkg[2] = {2.2,5.}; - const Double_t kM2tpc[2] = {2.,6.5}; const Bool_t kEfficiency = 1; - const Bool_t kFitFrac = 0; - const Double_t kSysErr[2] = {0.10, 0.11} ; + const Bool_t kFitFrac = 1; + const Int_t kDebugLevel = 1; + + Double_t bkgLimitToF[2] = {-2, 2}; + Double_t pidLimitToF[2] = {-2.,6.}; - Double_t xsec[3]; - GetInelXSection(xsec, period); + Double_t bkgLimitM2[2] = {1.8,6.0}; + Double_t pidLimitM2[2] = {1.8,6.}; + + if(kPidProc==1) + { + for(Int_t i=0; i<2; ++i) + { + bkgLimitM2[i] -= 3.51792; + pidLimitM2[i] -= 3.51792; + } + } 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"; @@ -86,27 +101,41 @@ Int_t Config_Deuteron_TOF_LHC10x(const TString& inputDir = "~/alice/input", 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(kPidM2); - driver.SetM2BinInterval(kM2Bin[0], kM2Bin[1]); - driver.SetM2BkgInterval(kM2Bkg[0], kM2Bkg[1]); - driver.SetM2TPCInterval(kM2tpc[0], kM2tpc[1]); - driver.SetUnfolding(kUnfolding, kIter); + driver.SetPtInterval(ptmin, ptmax); + driver.SetPid(kPid); + driver.SetPidProcedure(kPidProc); + driver.SetPidEfficiency(kPidEff); + driver.SetPidPt(ptpid); + + if(kPidProc==2) + { + driver.SetBkgInterval(bkgLimitToF[0], bkgLimitToF[1]); + driver.SetPidInterval(pidLimitToF[0], pidLimitToF[1]); + } + else + { + driver.SetBkgInterval(bkgLimitM2[0], bkgLimitM2[1]); + driver.SetPidInterval(pidLimitM2[0], pidLimitM2[1]); + } + driver.SetSecondaries(kSecondaries); - driver.SetSecProd(kSecProd); + driver.SetSecProcedure(kSecProc); driver.SetMatDCAxyModel(kMatDCAxyMod); driver.SetAntiNucleusAsTemplate(kAntiNucTemplate); driver.SetNBin(kNbin); driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]); driver.SetEfficiency(kEfficiency,0); driver.SetFitFractionCorr(kFitFrac); - driver.SetSysErr(kSysErr[0],kSysErr[1]); + driver.SetDebugLevel(kDebugLevel); driver.SetMakeStats(makeStats); driver.SetMakeCorrections(makeCor); @@ -120,15 +149,11 @@ Int_t Config_Deuteron_TOF_LHC10x(const TString& inputDir = "~/alice/input", if(!drawOutput) return 0; - TStyle* st = GetDrawingStyle(); - st->cd(); - gROOT->ForceStyle(); - DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag()); - 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])); + 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])); - DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, kPidM2, hiPtBin, kM2Bin[0], kM2Bin[1]); + DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, kPid, ptmax, ptpid); DrawOutputRatio(outputRatio, outputTag, kSpecies); DrawOutputSpectra(outputSpectra, outputTag, kSpecies); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TPC_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TPC_LHC10x.C index 610745f6edd..7beedb95252 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TPC_LHC10x.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TPC_LHC10x.C @@ -16,59 +16,59 @@ // LHC10x config for deuterons and antideuterons // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include "AliLnDriver.h" +#endif + #include "Config.h" -Int_t Config_Deuteron_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 = 3, // 0.4 Gev/c - Int_t hiPtBin = 6, // 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_Deuteron_TPC_LHC10x( const TString& inputDir = "~/alice/input" + , const TString& outputDir = "~/alice/output" + , const TString& period = "lhc10d" + , const TString& outputTag = "lhc10d" + , const TString& trkselTag = "-tpc3-nsd-moc" + , const TString& multTag = "" + , const TString& multCorTag = "" + , Double_t ymax = 0.5 + , Bool_t inel = 0 + , Bool_t drawOutput = 1 + , Double_t ptmin = 0.5 // GeV/c + , Double_t ptmax = 1.1 // GeV/c + , 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 deuterons and antideuterons // (TPC) // const TString kSpecies = "Deuteron"; - const TString kTrkSel = "its_tpc_dca-tpc3"; - const TString kTrigName = "mbor"; + const TString kTrkSel = "its_tpc_dca"; + const TString kTrigName = "mband"; const Bool_t kMCtoINEL = 1; - const Bool_t kUnfolding = 1; - const Int_t kIter = 7; + const Double_t kPidEff = 1.; 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 Bool_t kAntiNucTemplate = 0; const Int_t kNbin = 5; - const Double_t kDCAxy[2] = {-0.2,0.2}; + const Double_t kDCAxy[2] = {-0.2, 0.2}; const Bool_t kEfficiency = 1; - const Bool_t kFitFrac = 0; - const Double_t kSysErr[2] = {0.10, 0.11} ; - - 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) + ".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"; @@ -83,23 +83,24 @@ Int_t Config_Deuteron_TPC_LHC10x(const TString& inputDir = "~/alice/input", 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.SetPidEfficiency(kPidEff); driver.SetSecondaries(kSecondaries); - driver.SetSecProd(kSecProd); + driver.SetSecProcedure(kSecProc); driver.SetMatDCAxyModel(kMatDCAxyMod); driver.SetAntiNucleusAsTemplate(kAntiNucTemplate); driver.SetNBin(kNbin); driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]); driver.SetEfficiency(kEfficiency,0); - driver.SetFitFractionCorr(kFitFrac); - driver.SetSysErr(kSysErr[0],kSysErr[1]); driver.SetMakeStats(makeStats); driver.SetMakeCorrections(makeCor); @@ -107,19 +108,17 @@ Int_t Config_Deuteron_TPC_LHC10x(const TString& inputDir = "~/alice/input", driver.SetMakeRatio(makeRatio); driver.SetMakeSpectra(makeSpectra); + driver.SetDebugLevel(kDebugLevel); + driver.Run(); // draw output if(!drawOutput) return 0; - TStyle* st = GetDrawingStyle(); - st->cd(); - gROOT->ForceStyle(); - DrawOutputCorr(kSpecies, inputCorr, driver.GetOutputCorrTag()); - 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])); + 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])); DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0); DrawOutputRatio(outputRatio, outputTag, kSpecies); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C index cabd68f20d0..075e7a509a6 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C @@ -16,59 +16,62 @@ // LHC10x config for protons and antiprotons // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #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"; @@ -83,28 +86,30 @@ Int_t Config_Proton_TOF_LHC10x(const TString& inputDir = "~/alice/input", 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); @@ -118,13 +123,9 @@ Int_t Config_Proton_TOF_LHC10x(const TString& inputDir = "~/alice/input", 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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C index 9b48164e0d3..f61b566f8d0 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C @@ -16,61 +16,60 @@ // LHC10x config for protons and antiprotons // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #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"; @@ -85,24 +84,26 @@ Int_t Config_Proton_TPC_LHC10x(const TString& inputDir = "~/alice/input", 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); @@ -116,13 +117,9 @@ Int_t Config_Proton_TPC_LHC10x(const TString& inputDir = "~/alice/input", 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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_TPCTOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_TPCTOF_LHC10x.C index e6facb170d5..282dc039b8e 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_TPCTOF_LHC10x.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_TPCTOF_LHC10x.C @@ -16,33 +16,35 @@ // TPC+TOF config // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include - #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")) @@ -58,6 +60,7 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir = "~/alice/input", + "\"" + outputDir + "\"," + "\"" + period + "\"," + "\"" + kOutputTagTPC + "\"," + + "\"" + trkselTag + "\"," + "\"" + multTag + "\"," + "\"" + multCorTag; @@ -65,28 +68,36 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir = "~/alice/input", + "\"" + 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"; @@ -102,20 +113,18 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir = "~/alice/input", 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); @@ -130,10 +139,6 @@ Int_t Config_TPCTOF_LHC10x(const TString& inputDir = "~/alice/input", if(!drawOutput) return 0; - TStyle* st = GetDrawingStyle(); - st->cd(); - gROOT->ForceStyle(); - DrawOutputRatio(outputRatio, outputTag, species); DrawOutputSpectra(outputSpectra, outputTag, species); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C index 5966121d23e..47a2479ed0c 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C @@ -16,12 +16,13 @@ // Draw B2 // author: Eulogio Serradilla -#include +#if !defined(__CINT__) || defined(__MAKECINT__) +#include #include #include #include #include -#include +#endif #include "B2.h" @@ -30,20 +31,10 @@ void DrawB2(const TString& inputFile="b2.root", const TString& tag="test", const // // 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); @@ -58,41 +49,29 @@ void DrawB2(const TString& inputFile="b2.root", const TString& tag="test", const // 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(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(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"); } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C index 2f9e737bc9a..401edcbc1c1 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C @@ -16,24 +16,31 @@ // Draw pt corrections for debugging // author: Eulogio Serradilla -#include -#include +#if !defined(__CINT__) || defined(__MAKECINT__) +#include #include #include -#include #include #include #include +#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; @@ -44,23 +51,6 @@ void DrawCorr(const TString& species="Deuteron", const TString& inputFile="corre 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; icd(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())); @@ -72,9 +62,9 @@ void DrawCorr(const TString& species="Deuteron", const TString& inputFile="corre for(Int_t i=0; i(finput, tag, kPrefix[i] + species + "_Eff_Trig_Pt"); + hEffVtxPt[i] = FindObj(finput, tag, kPrefix[i] + species + "_Eff_Vtx_Pt"); + hEffAccTrkPt[i] = FindObj(finput, tag, kPrefix[i] + species + "_Eff_AccTrk_Pt"); } c2->cd(1); @@ -94,11 +84,11 @@ void DrawCorr(const TString& species="Deuteron", const TString& inputFile="corre for(Int_t i=0; i(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Pt"); + TH1D* hFracMatPt = FindObj(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(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Fit_Pt"); + TF1* fncFracMatPt = FindObj(finput, tag, kPrefix[i] + species + "_Frac_Mat_Fit_Pt"); c3->cd(2*i+1); hFracFdwnPt->SetAxisRange(xmin, xmax, "X"); @@ -130,12 +120,12 @@ void DrawCorr(const TString& species="Deuteron", const TString& inputFile="corre } } -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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C index 42f3e80162b..5cd24a859ab 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C @@ -16,7 +16,7 @@ // Draw histograms/graphs with same name located in different directories // author: Eulogio Serradilla -#include +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include @@ -28,6 +28,7 @@ #include #include #include +#endif TGraphErrors* Divide(const TGraphErrors* grX1, const TGraphErrors* grX2, const TString& name); @@ -54,7 +55,7 @@ void DrawDir(const TString& inputFile, 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); @@ -152,25 +153,10 @@ void DrawDir(const TString& inputFile, // 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()); @@ -202,7 +188,7 @@ void DrawDir(const TString& inputFile, } else if(obj[i]->InheritsFrom("TGraph")) { - obj[i]->Draw("PZ"); + obj[i]->Draw("zP"); } } @@ -248,22 +234,24 @@ void DrawDir(const TString& inputFile, 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; iDraw("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; iDraw("zP"); + } + if(option == 2) { TLegend* legendRatio = new TLegend(0.5718391,0.6991525,0.8390805,0.8368644,0,"brNDC"); @@ -277,6 +265,7 @@ void DrawDir(const TString& inputFile, legendRatio->AddEntry(grDiv[i], subdir[i]->Data(), "lp"); } + legendRatio->SetTextFont(42); legendRatio->Draw(); } } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C index 2fb7ffcc5e2..63e03742d62 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C @@ -16,25 +16,26 @@ // Draw corrected pt // author: Eulogio Serradilla -#include +#if !defined(__CINT__) || defined(__MAKECINT__) #include +#include #include #include #include #include #include #include - #include #include #include #include -#include #include +#include +#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 @@ -42,16 +43,26 @@ void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", co 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(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 @@ -59,7 +70,7 @@ void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", co 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())); @@ -67,37 +78,36 @@ void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", co TGraphErrors* grFitM2Pt = new TGraphErrors(); grFitM2Pt->SetName(Form("%s_Fit_M2_Pt", particle.Data())); - for(Int_t i=lowm2bin, j=0; icd(i-lowm2bin+1); - //gPad->SetLogy(0); - RooWorkspace* w= (RooWorkspace*)FindObj(finput, tag, Form("%s_M2_%02d",particle.Data(),i)); + RooWorkspace* w= FindObj(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())); @@ -105,38 +115,39 @@ void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", co 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); @@ -147,20 +158,19 @@ void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", co for(Int_t i=0; i(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; iDraw("sameE"); + for(Int_t i=0; iDraw("sameE"); legend->Draw(); + + c2->Update(); } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C index ef34ff52c17..c5ed9ea8db6 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C @@ -16,12 +16,13 @@ // Draw antiparticle/particle ratio // author: Eulogio Serradilla -#include -#include +#if !defined(__CINT__) || defined(__MAKECINT__) +#include #include #include #include #include +#endif #include "B2.h" @@ -35,6 +36,13 @@ void DrawRatio(const TString& inputFile="pt.root", const TString& tag="test", co 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); @@ -47,7 +55,7 @@ void DrawRatio(const TString& inputFile="pt.root", const TString& tag="test", co 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(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(); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C index cc17cb795ac..fd74d1fd396 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C @@ -16,32 +16,37 @@ // Draw corrections of secondaries // author: Eulogio Serradilla -#include +#if !defined(__CINT__) || defined(__MAKECINT__) #include +#include #include -#include #include #include #include #include +#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); @@ -61,17 +66,22 @@ void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", c TH1D* hMatFit[kNBin]; TH1D* hFdwnFit[kNBin]; + TH1D* hPidPt = FindObj(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(finput, tag, particle + Form("_Data_DCAxy_%02d",i)); + hDataMC[i] = FindObj(finput, tag, particle + Form("_SimData_DCAxy_%02d",i)); + hPrim[i] = FindObj(finput, tag, particle + Form("_Prim_DCAxy_%02d",i)); + hMat[i] = FindObj(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(finput, tag, particle + Form("_Fdwn_DCAxy_%02d",i)); } else { @@ -79,13 +89,13 @@ void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", c } // 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(finput, tag, particle + Form("_Fit_Data_DCAxy_%02d",i)); + hPrimFit[i] = FindObj(finput, tag, particle + Form("_Fit_Prim_DCAxy_%02d",i)); + hMatFit[i] = FindObj(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(finput, tag, particle + Form("_Fit_Fdwn_DCAxy_%02d",i)); } else { @@ -93,31 +103,30 @@ void DrawSec(const TString& inputFile="debug.root", const TString& tag="test", c } } - // 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(); } } @@ -130,7 +139,7 @@ Int_t GetNumberCols(Int_t lowbin, Int_t hibin) 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 @@ -144,10 +153,9 @@ TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t low { 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"); @@ -178,7 +186,7 @@ TCanvas* DrawMC(TH1D** hData, TH1D** hPrim, TH1D** hMat, TH1D** hFdwn, Int_t low 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 @@ -188,15 +196,13 @@ TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D* 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"); @@ -230,7 +236,7 @@ TCanvas* DrawFit(TH1D** hData, TH1D** hDataFit, TH1D** hPrim, TH1D** hMat, TH1D* 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 @@ -244,13 +250,11 @@ TCanvas* DrawGoF(TH1D** hData, TH1D** hDataFit, Int_t lowbin, Int_t hibin, const { 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); diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C index 15292846e03..dfb0eb7536d 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C @@ -16,14 +16,13 @@ // Draw differential yields // author: Eulogio Serradilla -#include -#include +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include #include -#include +#endif #include "B2.h" @@ -32,65 +31,43 @@ void DrawSpectra(const TString& inputFile="spectra.root", const TString& tag="lh // // 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(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(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"); } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C index bac11e0b38b..79bcc9d2362 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C @@ -16,24 +16,28 @@ // call Config_XXX for each period // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include #include - #include "AliLnDriver.h" +#endif + #include "Config.h" -Int_t LHC10bcde(const TString& species = "Deuteron", - const TString& inputDir = "~/alice/input", - const TString& outputDir = "~/alice/output", - const TString& outputTag = "lhc10bcde", - const TString& multTag = "", - const TString& multCorTag = "", - Bool_t inel = 1, // for mult - Bool_t drawOutput = 1, // for batch - Int_t option = 2) +Int_t LHC10bcde( const TString& species = "Deuteron" + , const TString& inputDir = "~/alice/input" + , const TString& outputDir = "~/alice/output" + , const TString& outputTag = "lhc10bcde" + , const TString& trkselTag = "-tpc3-nsd-moc" + , const TString& multTag = "" + , const TString& multCorTag = "" + , Double_t ymax = 0.5 + , Bool_t inel = 0 // for mult + , Bool_t drawOutput = 1 // for batch + , Int_t option = 2) { // // call Config_XXX for each period, merge the corrected pt and then get the results @@ -42,12 +46,10 @@ Int_t LHC10bcde(const TString& species = "Deuteron", const TString kPeriod[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" }; const TString kOutputTag[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" }; - Int_t lowbin = (species=="Proton") ? 5 : 4; - Int_t jointbin = (species=="Proton") ? 11 : 6; - Int_t hibin = (species=="Proton") ? 36 : 15; - - Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8; - Double_t ymax = (species=="Proton") ? 4.e-1 : 4.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 : 3.0; + Double_t ptpid = (species=="Proton") ? 3.5 : 1.6; using namespace std; @@ -68,6 +70,7 @@ Int_t LHC10bcde(const TString& species = "Deuteron", + "\"" + outputDir + "\"," + "\"" + kPeriod[i] + "\"," + "\"" + kOutputTag[i] + "\"," + + "\"" + trkselTag + "\"," + "\"" + multTag + "\"," + "\"" + multCorTag; @@ -75,15 +78,15 @@ Int_t LHC10bcde(const TString& species = "Deuteron", { 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(), inel)); + gROOT->ProcessLine(Form(".x Config_%s_TPC_LHC10x.C+g(\"%s\", %f, %d, 0)", species.Data(), arg.Data(), ymax, inel)); 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(), inel)); + gROOT->ProcessLine(Form(".x Config_%s_TOF_LHC10x.C+g(\"%s\", %f, %d, 0)", species.Data(), arg.Data(), ymax, inel)); 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(), inel, 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, inel, species.Data(), ptmin, ptjoint, ptmax, ptpid)); break; } @@ -91,9 +94,9 @@ Int_t LHC10bcde(const TString& species = "Deuteron", m.AddFile(ptfile,0); } - TString outputPt = outputDir + "/" + species + "-" + outputTag + "-Pt.root"; - TString outputRatio = outputDir + "/" + species + "-" + outputTag + "-Ratio.root"; - TString outputSpectra = outputDir + "/" + species + "-" + outputTag + "-Spectra.root"; + TString outputPt = outputDir + "/" + species + "-" + outputTag + multTag + "-Pt.root"; + TString outputRatio = outputDir + "/" + species + "-" + outputTag + multTag + "-Ratio.root"; + TString outputSpectra = outputDir + "/" + species + "-" + outputTag + multTag + "-Spectra.root"; // pt @@ -109,6 +112,9 @@ Int_t LHC10bcde(const TString& species = "Deuteron", driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra); driver.SetOutputTag(outputTag); + driver.SetRapidityInterval(-ymax,ymax); + driver.SetExtrapolateToINEL(inel); + driver.SetMakeCorrections(0); driver.SetMakePt(0); driver.SetMakeRatio(1); @@ -129,8 +135,8 @@ Int_t LHC10bcde(const TString& species = "Deuteron", m2.AddFile(outputRatio.Data(),0); m3.AddFile(outputSpectra.Data(),0); - TString allRatios = outputDir + "/" + species + "-lhc10bcde1-Ratio.root"; - TString allSpectra = outputDir + "/" + species + "-lhc10bcde1-Spectra.root"; + TString allRatios = outputDir + "/" + species + "-" + outputTag + "-2" + multTag + "-Ratio.root"; + TString allSpectra = outputDir + "/" + species + "-" + outputTag + "-2" + multTag + "-Spectra.root"; m2.OutputFile(allRatios.Data()); m3.OutputFile(allSpectra.Data()); @@ -142,9 +148,12 @@ Int_t LHC10bcde(const TString& species = "Deuteron", if(!drawOutput) return 0; - 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())); + 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())); + + Double_t minYield = (species=="Proton") ? 1.1e-6 : 2.e-8; + Double_t maxYield = (species=="Proton") ? 4.e-1 : 9.e-4; - DrawOutputSpectraMult(allSpectra, species, ymin, ymax, 1, "lhc10bcde"); + DrawOutputSpectraMult(allSpectra, species, minYield, maxYield, 2, outputTag); return 0; } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C index 7a7bc5d3108..5b9bd50d8b3 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C @@ -16,20 +16,24 @@ // LHC10x multiplicity config for protons and deuterons // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) #include #include #include #include - #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 @@ -44,13 +48,10 @@ Int_t LHC10xMult( const TString& species = "Proton" 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"))) { @@ -66,32 +67,36 @@ Int_t LHC10xMult( const TString& species = "Proton" for(Int_t i=0; iProcessLine(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; } @@ -123,9 +128,12 @@ Int_t LHC10xMult( const TString& species = "Proton" // 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; } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C index bad3efb8de2..9dac990cd3c 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C @@ -16,6 +16,8 @@ // d/p ratio // author: Eulogio Serradilla +#if !defined(__CINT__) || defined(__MAKECINT__) +#include #include #include #include @@ -27,6 +29,9 @@ #include #include #include +#include +#include +#endif #include "B2.h" #include "Config.h" @@ -34,19 +39,20 @@ 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 // @@ -55,6 +61,12 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- 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"); @@ -77,12 +89,16 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- {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]; @@ -91,39 +107,87 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- Double_t dNdy[kNspec][kNpart]; Double_t dNdyErr[kNspec][kNpart]; + cout << endl << xtag[i] << endl << endl; + for(Int_t j=0; j(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: %g +/- %g\n", kParticle[j][k].Data(), status, dNdy[j][k], dNdyErr[j][k], 1.-extrFrac, extrFracErr, meanPt, meanPtErr); } } @@ -160,8 +224,6 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- 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()); } } } @@ -196,31 +258,20 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- 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]; @@ -239,9 +290,6 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- 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"); } } @@ -262,21 +310,30 @@ void RatioMult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-nsd- // spectra - TCanvas* c2[kNspec]; + TPaveText* label[nx][kNspec][kNpart]; + + TCanvas* c2[kNspec][kNpart]; for(Int_t j=0; jDivide(nx,2); - for(Int_t k=0; kDivide(3,2); + for(Int_t i=0; iSetLogy(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(); } } } @@ -315,5 +372,6 @@ void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TS gr->SetLineColor(color); gr->GetXaxis()->SetTitle(xtitle.Data()); gr->GetYaxis()->SetTitle(ytitle.Data()); - gr->Draw("AP"); + gr->GetYaxis()->SetTitleOffset(1.3); + gr->Draw("zAP"); } diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/run.C b/PWGLF/SPECTRA/Nuclei/B2/macros/run.C index 617e2f6dae8..669fdd41f4f 100644 --- a/PWGLF/SPECTRA/Nuclei/B2/macros/run.C +++ b/PWGLF/SPECTRA/Nuclei/B2/macros/run.C @@ -9,12 +9,10 @@ Int_t run(const TString& config) 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"); @@ -23,6 +21,7 @@ Int_t run(const TString& config) 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())); -- 2.43.0