B2 analysis code
authoreserradi <eserradi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 Jan 2013 19:25:32 +0000 (19:25 +0000)
committereserradi <eserradi@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 11 Jan 2013 19:25:32 +0000 (19:25 +0000)
40 files changed:
PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/B2.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/Config.h [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/B2.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPCTOF_LHC10x.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/run.C [new file with mode: 0644]
PWGLF/SPECTRA/Nuclei/B2/macros/runLocalExample.C [new file with mode: 0644]

diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx
new file mode 100644 (file)
index 0000000..15a9d61
--- /dev/null
@@ -0,0 +1,236 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// coalescence parameter
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TMath.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+
+#include "AliLnB2.h"
+#include "B2.h"
+
+ClassImp(AliLnB2)
+
+AliLnB2::AliLnB2(const TString& protonSpectra, const TString& protonTag, const TString& nucleusSpectra, const TString& nucleusTag, const TString& outputFilename, const TString& otag, Int_t a, Int_t z)
+: TObject()
+, fProtonSpectra(protonSpectra)
+, fProtonTag(protonTag)
+, fNucleusSpectra(nucleusSpectra)
+, fNucleusTag(nucleusTag)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+, fA(2)
+, fZ(1)
+, fNucleusName("Deuteron")
+, fCd(0.15)
+{
+//
+// constructor
+//
+       this->SetNucleus(a,z);
+}
+
+AliLnB2::~AliLnB2()
+{
+//
+// destructor
+//
+}
+
+void AliLnB2::SetNucleus(Int_t a, Int_t z)
+{
+//
+// set nucleus mass and name
+//
+       fA = a;
+       fZ = z;
+       
+       Int_t zz = TMath::Abs(z);
+       
+       if(a==1 && zz==1)      fNucleusName = "Proton";
+       else if(a==2 && zz==1) fNucleusName = "Deuteron";
+       else if(a==3 && zz==1) fNucleusName = "Triton";
+       else if(a==3 && zz==2) fNucleusName = "He3";
+       else if(a==4 && zz==2) fNucleusName = "Alpha";
+       else
+       {
+               this->Warning("SetNucleus", "unknown nucleus A = %d, Z = %d", a, z);
+               fNucleusName = "Unknown";
+       }
+}
+
+Int_t AliLnB2::Run()
+{
+//
+// coalescence parameter
+//
+       TFile* finput1 = new TFile(fProtonSpectra.Data(), "read");
+       if (finput1->IsZombie()) exit(1);
+       
+       TFile* finputA = new TFile(fNucleusSpectra.Data(), "read");
+       if (finputA->IsZombie()) exit(1);
+       
+       TString prefix = "";
+       TString suffix = "";
+       if(fZ < 0)
+       {
+               prefix = "Anti";
+               suffix = "bar";
+       }
+       
+       // 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");
+       
+       TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
+       
+       foutput->mkdir(fOutputTag.Data());
+       foutput->cd(fOutputTag.Data());
+       
+       // coalescence parameter
+       
+       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
+       
+       delete grB2Pt;
+       delete grSysErrB2Pt;
+       delete grR3Pt;
+       delete grSysErrR3Pt;
+       
+       delete foutput;
+       delete finput1;
+       delete finputA;
+       
+       return 0;
+}
+
+TGraphErrors* AliLnB2::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name, Double_t errPt) const
+{
+//
+// coalescence parameter
+//
+       TGraphErrors* grBAPt = new TGraphErrors();
+       grBAPt->SetName(name.Data());
+       
+       Int_t offset = 0;
+       Double_t pt, yPrt, ptNuc, yNuc;
+       
+       grPrtInvDYieldPt->GetPoint(0,pt,yPrt);
+       
+       // find offset
+       for(Int_t i=0; i < grNucInvDYieldPt->GetN(); ++i)
+       {
+               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(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 eNuc = grNucInvDYieldPt->GetErrorY(i);
+               
+               Double_t errBA = bA*TMath::Sqrt(TMath::Power(eNuc/yNuc,2) + TMath::Power(fA*ePrt/yPrt,2));
+               
+               grBAPt->SetPoint(j, pt, bA);
+               grBAPt->SetPointError(j++, errPt, errBA);
+       }
+       
+       return grBAPt;
+}
+
+Double_t AliLnB2::Rside2Rlong(Double_t pt, Double_t B2, Double_t Cd) const
+{
+//
+// Rside^2*Rlong from B2 value
+// Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
+// (for pp ignore the exponential term)
+//
+       Double_t hbarc = 0.1973269631; // GeV*fm
+       Double_t m = 0.938272013;  // GeV/c^2
+       Double_t pi = TMath::Pi();
+       
+       if(B2==0) return 0;
+       
+       Double_t mt = TMath::Sqrt(pt*pt + m*m);
+       Double_t r3 = 3.*TMath::Power(pi,3./2.)*TMath::Power(hbarc,3.)*Cd/(2.*mt*B2);
+       
+       return r3;
+}
+
+TGraphErrors* AliLnB2::Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd, Double_t errPt) const
+{
+//
+// Rside^2*Rlong from B2 value
+// Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
+// (for pp ignore the exponential term)
+//
+       TGraphErrors* grR3 = new TGraphErrors(*grB2);
+       grR3->SetName(name.Data());
+       
+       Double_t pt, b2;
+       
+       for(Int_t i=0; i < grB2->GetN(); ++i)
+       {
+               grB2->GetPoint(i, pt, b2);
+               
+               if(b2==0) continue;
+               Double_t r3 = Rside2Rlong(pt, b2, Cd);
+               
+               grR3->SetPoint(i, pt, r3);
+               
+               // only statistical error propagation of B2
+               Double_t errB2 = grR3->GetErrorY(i);
+               Double_t errR = r3*errB2/b2;
+               
+               grR3->SetPointError(i, errPt, errR);
+       }
+       
+       return grR3;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.h
new file mode 100644 (file)
index 0000000..dee42d2
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALILNB2_H
+#define ALILNB2_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// coalescence parameter
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+
+class TGraphErrors;
+
+class AliLnB2: public TObject
+{
+  public:
+       
+       AliLnB2(const TString& protonSpectra, const TString& nucleusSpectra, const TString& protonTag, const TString& nucleusTag, const TString& outputFilename, const TString& otag, Int_t a, Int_t z);
+       
+       virtual ~AliLnB2();
+       
+       TGraphErrors* GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name, Double_t errPt=0) 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;
+       
+       Int_t Run();
+       
+       void SetNucleus(Int_t a, Int_t z);
+       void SetCd(Double_t Cd) { fCd = Cd; }
+       
+  private:
+       AliLnB2(const AliLnB2& other);
+       AliLnB2& operator=(const AliLnB2& other);
+       
+  private:
+       TString fProtonSpectra; // proton spectra filename
+       TString fProtonTag; // tag for proton file
+       TString fNucleusSpectra; // nucleus spectra filename
+       TString fNucleusTag; // tag for nucleus file
+       
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for output file
+       
+       Int_t fA; // nucleus mass
+       Int_t fZ; // nucleus charge
+       TString fNucleusName; // nucleus name
+       
+       Double_t fCd; // correction factor for homogeneity volume
+       
+       ClassDef(AliLnB2,1)
+};
+
+#endif // ALILNB2_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx
new file mode 100644 (file)
index 0000000..e812533
--- /dev/null
@@ -0,0 +1,92 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// driver for rebuilding the correction file
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TSystem.h>
+#include <TString.h>
+#include <TROOT.h>
+#include <TFileMerger.h>
+
+#include "AliLnUnfolding.h"
+#include "AliLnFakeTracks.h"
+#include "AliLnSecondaries.h"
+#include "AliLnEfficiency.h"
+#include "AliLnCorr.h"
+
+ClassImp(AliLnCorr)
+
+AliLnCorr::AliLnCorr(const TString& particle, const TString& dataFilename, const TString& simuFilename, const TString& simuFixFilename, const TString& outputFilename, const TString& otag)
+: TObject()
+, fOutputFilename(outputFilename)
+{
+//
+// constructor
+//
+
+       fUnfolding   = new AliLnUnfolding(particle, simuFilename, Form("Unfolding-%s",outputFilename.Data()), otag);
+       fFakeTracks  = new AliLnFakeTracks(particle, simuFilename, Form("FakeTracks-%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
+}
+
+AliLnCorr::~AliLnCorr()
+{
+//
+// destructor
+//
+       delete fUnfolding;
+       delete fFakeTracks;
+       delete fSecondaries;
+       delete fEfficiency;
+}
+
+Int_t AliLnCorr::Exec()
+{
+//
+// rebuild correction file
+//
+       fUnfolding->Exec();
+       fFakeTracks->Exec();
+       fSecondaries->Exec();
+       fEfficiency->Exec();
+       
+       // merge the root files
+       
+       TString output1 = *fUnfolding->GetOutputFilename();
+       TString output2 = *fFakeTracks->GetOutputFilename();
+       TString output3 = *fSecondaries->GetOutputFilename();
+       TString output4 = *fEfficiency->GetOutputFilename();
+       
+       TFileMerger m;
+       
+       m.AddFile(output1.Data(),0);
+       m.AddFile(output2.Data(),0);
+       m.AddFile(output3.Data(),0);
+       m.AddFile(output4.Data(),0);
+       
+       m.OutputFile(fOutputFilename.Data());
+       
+       m.Merge();
+       
+       // remove tmp files
+       gSystem->Exec(Form("rm -f %s %s %s %s", output1.Data(), output2.Data(),  output3.Data(),  output4.Data()));
+       
+       gSystem->Exec(Form("mv debug-%s debug-%s",output3.Data(),fOutputFilename.Data()));
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.h
new file mode 100644 (file)
index 0000000..d183159
--- /dev/null
@@ -0,0 +1,49 @@
+#ifndef ALILNCORR_H
+#define ALILNCORR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// driver for rebuilding the correction file
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class AliLnUnfolding;
+class AliLnFakeTracks;
+class AliLnSecondaries;
+class AliLnEfficiency;
+
+class AliLnCorr: public TObject
+{
+  public:
+       
+       AliLnCorr(const TString& particle, const TString& dataFilename, const TString& simuFilename, const TString& simuFixFilename, const TString& outputFilename, const TString& otag);
+       virtual ~AliLnCorr();
+       
+       Int_t Exec();
+       
+       AliLnUnfolding*   GetLnUnfolding() { return fUnfolding; }
+       AliLnFakeTracks*  GetLnFakeTracks() { return fFakeTracks; }
+       AliLnSecondaries* GetLnSecondaries() { return fSecondaries; }
+       AliLnEfficiency*  GetLnEfficiency() { return fEfficiency; }
+       
+  private:
+       AliLnCorr(const AliLnCorr& other);
+       AliLnCorr& operator=(const AliLnCorr& other);
+       
+  private:
+       TString fOutputFilename; // output filename
+       
+       AliLnUnfolding*   fUnfolding; // unfolding correction
+       AliLnFakeTracks*  fFakeTracks; // fake track correction
+       AliLnSecondaries* fSecondaries; // secondaries correction
+       AliLnEfficiency*  fEfficiency; // eficiency correction
+       
+       ClassDef(AliLnCorr,1)
+};
+
+#endif // ALILNCORR_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx
new file mode 100644 (file)
index 0000000..c7f3b43
--- /dev/null
@@ -0,0 +1,334 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// driver for computing the pt and spectra
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TObject.h>
+#include <TSystem.h>
+#include <TFileMerger.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TError.h>
+
+#include "AliLnCorr.h"
+#include "AliLnPt.h"
+#include "AliLnSecondaries.h"
+#include "AliLnEfficiency.h"
+#include "AliLnRatio.h"
+#include "AliLnSpectra.h"
+#include "AliLnB2.h"
+#include "AliLnDriver.h"
+#include "B2.h"
+
+ClassImp(AliLnDriver)
+
+AliLnDriver::AliLnDriver()
+: TObject()
+, fSpecies("Deuteron")
+, fOutputTag("test")
+, fOutputCorTag("")
+, fIsOnlyGen(0)
+, fNormToInel(1)
+, fMakeCorr(1)
+, fMakePt(1)
+, fMakeRatio(1)
+, fMakeSpectra(1)
+, fMakeStats(1)
+, fLowPtBin(3)
+, fHighPtBin(15)
+, fLowM2Bin(9)
+, fHighM2Bin(17)
+, fUnfolding(0)
+, fNIter(4)
+, fFakeTracks(1)
+, fSecondaries(1)
+, fSecProd(AliLnSecondaries::kTFractionFitter)
+, fMatDCAxyMod(AliLnSecondaries::kGeantDCAxy)
+, fNbin(1)
+, fYMin(-0.5)
+, fYMax(0.5)
+, fMinDCAxy(-1.5)
+, fMaxDCAxy(1.5)
+, fMinM2Bkg(2.2)
+, fMaxM2Bkg(5.)
+, fMinM2tpc(2.)
+, fMaxM2tpc(6.5)
+, fEfficiency(1)
+, fG3Fluka(0)
+, fScMat(1)
+, fScFd(1)
+, fSysPos(1)
+, fSysNeg(1)
+, fInputData()
+, fInputSimu()
+, fInputSimuFix()
+, fOutputPtCorr()
+, fOutputPt()
+, fOutputRatio()
+, fOutputSpectra()
+, fOutputPtCorrDebug()
+, fOutputPtDebug()
+, fFitFrac(1)
+, fSameFdwn(0)
+, fVtxCorr(0)
+, fVtxCorrVal(1)
+{
+//
+// constructor
+//
+       for(Int_t i=0; i<3; ++i)
+       {
+               fTrigEff[i] = 1;
+               fXsec[i] = 1;
+       }
+       
+       gErrorIgnoreLevel = kWarning;
+       // gErrorIgnoreLevel = kPrint, kInfo, kWarning, kError, kBreak, kSysError, kFatal;
+}
+
+AliLnDriver::~AliLnDriver()
+{
+//
+// destructor
+//
+}
+
+void AliLnDriver::SetInputFilenames(const TString& data, const TString& simu, const TString& simuFix, const TString& ptcorr)
+{
+//
+// input filenames to get the pt
+//
+       fInputData = data;
+       fInputSimu = simu;
+       fInputSimuFix = simuFix;
+       fOutputPtCorr = ptcorr;
+       
+       fOutputPtCorrDebug = ptcorr;
+       fOutputPtCorrDebug.Replace(ptcorr.Length()-5,5,"-debug.root");
+}
+
+void AliLnDriver::SetOutputFilenames(const TString& pt, const TString& ratio, const TString& spectra)
+{
+//
+// output filenames
+//
+       fOutputPt = pt;
+       fOutputRatio=ratio;
+       fOutputSpectra=spectra;
+       
+       fOutputPtDebug = pt;
+       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])
+{
+//
+// set trigger efficiency
+//
+       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();
+       
+       if(fMakePt) this->MakePt();
+       
+       if(fMakeRatio) this->MakeRatio();
+       
+       if(fMakeSpectra) this->MakeSpectra();
+       
+       return 0;
+}
+
+void AliLnDriver::MakePtCorr() const
+{
+//
+// make pt correction file
+//
+       const Int_t kNpart = 2;
+       const TString kParticle[kNpart] = {fSpecies, Form("Anti%s",fSpecies.Data())};
+       
+       TFileMerger m1, m2;
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               TString outputfile1 = kParticle[i] + "-corrections.root";
+               TString outputfile2 = Form("debug-%s", outputfile1.Data());
+               
+               AliLnCorr lncorr(kParticle[i], fInputData, fInputSimu, fInputSimuFix, outputfile1, fOutputCorTag);
+               
+               lncorr.GetLnSecondaries()->SetCorBins(fLowPtBin, fHighPtBin);
+               lncorr.GetLnSecondaries()->SetProcedure(fSecProd);
+               lncorr.GetLnSecondaries()->SetMatDCAxyModel(fMatDCAxyMod);
+               lncorr.GetLnSecondaries()->SetDCAxyInterval(fMinDCAxy, fMaxDCAxy);
+               lncorr.GetLnSecondaries()->SetNBin(fNbin);
+               lncorr.GetLnSecondaries()->SetScalingFactors(fScMat, fScFd);
+               
+               lncorr.GetLnEfficiency()->SetG3Fluka(fG3Fluka);
+               
+               lncorr.Exec();
+               
+               m1.AddFile(outputfile1.Data(),0);
+               m2.AddFile(outputfile2.Data(),0);
+       }
+       
+       // merge and remove tmp files
+       
+       m1.OutputFile(fOutputPtCorr.Data());
+       m1.Merge();
+       
+       gSystem->Exec(Form("rm -f %s-corrections.root Anti%s-corrections.root", fSpecies.Data(), fSpecies.Data()));
+       
+       m2.OutputFile(fOutputPtCorrDebug.Data());
+       m2.Merge();
+       
+       gSystem->Exec(Form("rm -f debug-%s-corrections.root debug-Anti%s-corrections.root", fSpecies.Data(), fSpecies.Data()));
+}
+
+void AliLnDriver::MakePt() const
+{
+//
+// make pt from data and correction file
+//
+       const Int_t kNpart = 2;
+       
+       const TString kParticle[kNpart] = { fSpecies, Form("Anti%s",fSpecies.Data())};
+       
+       TFileMerger m, m2;
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               TString ptfile   = kParticle[i] + "-Pt.root";
+               
+               AliLnPt lnpt(kParticle[i], fTrigEff[0], fInputData, ptfile, fOutputTag, fOutputPtCorr, fOutputCorTag);
+               
+               lnpt.SetOnlyGeneration(fIsOnlyGen);
+               lnpt.SetRapidityInterval(fYMin, fYMax);
+               
+               lnpt.SetPtBinInterval(fLowPtBin, fHighPtBin);
+               lnpt.SetM2BinInterval(fLowM2Bin, fHighM2Bin);
+               lnpt.SetM2BkgInterval(fMinM2Bkg, fMaxM2Bkg);
+               lnpt.SetM2TPCInterval(fMinM2tpc, fMaxM2tpc);
+               lnpt.SetPidM2(fPidM2);
+               lnpt.SetUnfolding(fUnfolding, fNIter);
+               lnpt.SetFakeTracks(fFakeTracks);
+               lnpt.SetSecondaries(fSecondaries);
+               lnpt.SetEfficiency(fEfficiency);
+               lnpt.SetMakeStats(fMakeStats);
+               lnpt.SetVertexCorrection(fVtxCorr, fVtxCorrVal);
+               lnpt.SetFitFractionCorr(fFitFrac);
+               lnpt.SetSameFeedDownCorr(fSameFdwn);
+               
+               lnpt.Exec();
+               
+               m.AddFile(ptfile.Data(),0);
+               if(!fIsOnlyGen) m2.AddFile(Form("debug-%s",ptfile.Data()),0);
+       }
+       
+       // merge and remove tmp files
+       
+       m.OutputFile(fOutputPt.Data());
+       m.Merge();
+       
+       gSystem->Exec(Form("rm -f %s-Pt.root Anti%s-Pt.root", fSpecies.Data(), fSpecies.Data()));
+       
+       if(!fIsOnlyGen)
+       {
+               m2.OutputFile(fOutputPtDebug.Data());
+               m2.Merge();
+               
+               gSystem->Exec(Form("rm -f debug-%s-Pt.root debug-Anti%s-Pt.root", fSpecies.Data(),fSpecies.Data()));
+       }
+}
+
+void AliLnDriver::MakeRatio() const
+{
+//
+// make antiparticle/particle ratio
+//
+       AliLnRatio lnr(fSpecies, fOutputPt, fOutputTag, fOutputRatio, fOutputTag);
+       
+       lnr.Exec();
+}
+
+void AliLnDriver::MakeSpectra() const
+{
+//
+// make differential yields
+//
+       TFileMerger m;
+       
+       const TString kParticle[]  = { fSpecies, Form("Anti%s",fSpecies.Data())};
+       const Double_t kSysErr[] = { fSysPos, fSysNeg };
+       
+       for(Int_t i=0; i<2; ++i)
+       {
+               TString spectrafile = kParticle[i] + "-Spectra.root";
+               
+               AliLnSpectra lnspectra(kParticle[i], fOutputPt, fOutputTag, spectrafile, fOutputTag, fXsec);
+               
+               lnspectra.SetRapidityInterval(fYMin, fYMax);
+               lnspectra.SetNormalizeToINEL(fNormToInel);
+               lnspectra.SetOnlyGeneration(fIsOnlyGen);
+               lnspectra.SetScalingFactor(kSysErr[i]);
+               
+               lnspectra.Exec();
+               
+               m.AddFile(spectrafile.Data(),0);
+       }
+       
+       // merge and remove tmp files
+       
+       m.OutputFile(fOutputSpectra.Data());
+       m.Merge();
+       
+       gSystem->Exec(Form("rm -f %s-Spectra.root Anti%s-Spectra.root", fSpecies.Data(), fSpecies.Data()));
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.h
new file mode 100644 (file)
index 0000000..cb0291b
--- /dev/null
@@ -0,0 +1,152 @@
+#ifndef ALILNDRIVER_H
+#define ALILNDRIVER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// driver for computing the pt and spectra
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class TString;
+class TFileMerger;
+
+class AliLnDriver: public TObject
+{
+  public:
+       
+       AliLnDriver();
+       virtual ~AliLnDriver();
+       
+       void MakePtCorr() const;
+       void MakePt() const;
+       void MakeRatio() const;
+       void MakeSpectra() const;
+       
+       void PrintFilenames() const;
+       
+       Int_t Run() const;
+       
+       TString GetSpecies() const { return fSpecies; }
+       TString GetOutputTag() const { return fOutputTag; }
+       TString GetOutputCorrTag() const { return fOutputCorTag; }
+       
+       TString GetPtCorrDebugFilename() const { return fOutputPtCorrDebug; }
+       TString GetPtDebugFilename() const { return fOutputPtDebug; }
+       
+       void SetInputFilenames(const TString& data, const TString& simu, const TString& simuFix, const TString& ptcorr);
+       void SetOutputFilenames(const TString& pt, const TString& ratio, const TString& spectra);
+       void SetDebugFilenames(const TString& corrdebug, const TString& ptdebug) { fOutputPtCorrDebug = corrdebug; fOutputPtDebug = ptdebug; }
+       
+       void SetSpecies(const TString& species) { fSpecies = species; }
+       
+       void SetOutputTag(const TString& tag) { fOutputTag = tag; }
+       void SetOutputCorTag(const TString& tag) { fOutputCorTag = tag; }
+       
+       void SetTriggerEfficiency(Double_t eff[3]);
+       void SetInelXSection(Double_t xsec[3]);
+       void SetNormalizeToINEL(Bool_t flag=1) { fNormToInel = flag; }
+       void SetOnlyGeneration(Bool_t flag=1) { fIsOnlyGen = flag; }
+       
+       void SetMakeCorrections(Bool_t flag=1) { fMakeCorr = flag; }
+       void SetMakePt(Bool_t flag=1) { fMakePt = flag; }
+       void SetMakeRatio(Bool_t flag=1) { fMakeRatio = flag; }
+       void SetMakeSpectra(Bool_t flag=1) { fMakeSpectra = flag; }
+       void SetMakeStats(Bool_t flag=1) { fMakeStats = flag; }
+       
+       void SetRapidityInterval(Double_t ymin, Double_t ymax) { fYMin = ymin; fYMax = ymax; }
+       
+       void SetPtBinInterval(Int_t lowbin, Int_t hibin) { fLowPtBin = lowbin; fHighPtBin = 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 SetFakeTracks(Bool_t flag=1) { fFakeTracks = flag; }
+       void SetSecondaries(Bool_t flag=1) { fSecondaries = flag; }
+       void SetSecProd(Int_t prod) { fSecProd = prod; }
+       void SetMatDCAxyModel(Int_t model=1) { fMatDCAxyMod = model; }
+       void SetNBin(Int_t nbin) { fNbin = nbin; }
+       void SetDCAxyInterval(Double_t dcamin, Double_t dcamax) { fMinDCAxy = dcamin; fMaxDCAxy = dcamax; }
+       void SetEfficiency(Bool_t flag=1, Bool_t g3Fluka=0) { fEfficiency = flag; fG3Fluka=g3Fluka; }
+       void SetScalingFactors(Double_t mat, Double_t fd) { fScMat=mat; fScFd=fd; }
+       
+       void SetVertexCorrection(Bool_t flag=1, Double_t val=1) { fVtxCorr = flag; fVtxCorrVal=val; }
+       void SetFitFractionCorr(Bool_t flag=1) { fFitFrac=flag; }
+       void SetSameFeedDownCorr(Bool_t flag=1) { fSameFdwn = flag; }
+       
+       void SetSysErr( Double_t pos, Double_t neg) { fSysPos = pos; fSysNeg = neg; }
+       
+  private:
+       AliLnDriver(const AliLnDriver& other);
+       AliLnDriver& operator=(const AliLnDriver& other);
+       
+  private:
+       TString fSpecies;       // particle species
+       TString fOutputTag;     // tag for output file
+       TString fOutputCorTag;        // tag for correction file
+       
+       Double_t fTrigEff[3];   // trigger efficiency, stat. and syst. errors
+       Double_t fXsec[3];      // total inelastic cross section, stat. and syst. errors
+       Bool_t  fIsOnlyGen;     // if it is only generation
+       Bool_t  fNormToInel;    // normalize to inelastic events
+       
+       Bool_t  fMakeCorr;      // make corrections
+       Bool_t  fMakePt;        // make pt
+       Bool_t  fMakeRatio;     // make antiparticle/particle ratio
+       Bool_t  fMakeSpectra;   // make spectra
+       Bool_t  fMakeStats;     // make event stats
+       
+       Int_t    fLowPtBin;     // low pt bin
+       Int_t    fHighPtBin;    // 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
+       Bool_t   fFakeTracks;   // fake tracks correction
+       Bool_t   fSecondaries;  // correction of secondaries
+       Int_t    fSecProd;      // procedure for estimating fractions
+       Int_t    fMatDCAxyMod;  // DCAxy model for correction of secondaries
+       Int_t    fNbin;         // rebin of DCAxy distribution
+       Double_t fYMin;         // min rapidity
+       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
+       Bool_t   fEfficiency;   // efficiency correction
+       Bool_t   fG3Fluka;      // enable G3/Fluka correction
+       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
+       
+       TString fOutputPtCorr;  // output filename for pt corrections
+       TString fOutputPt;      // output filename for pt
+       TString fOutputRatio;   // output filename for antiparticle/particle ratio
+       TString fOutputSpectra; // output filename for differential yields
+       
+       TString fOutputPtCorrDebug; // output filename for debugging pt corrections
+       TString fOutputPtDebug;     // output filename for debugging pt
+       
+       Bool_t fFitFrac;        // fit for fraction of secondaries
+       Bool_t fSameFdwn;       // same feed-down correction for positives and negatives
+       Bool_t fVtxCorr;        // enable vertex correction from simulation
+       Double_t fVtxCorrVal;   // vertex correction value
+       
+       ClassDef(AliLnDriver,3)
+};
+
+#endif // ALILNDRIVER_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx
new file mode 100644 (file)
index 0000000..ad93b99
--- /dev/null
@@ -0,0 +1,152 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// efficiency correction
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TString.h>
+#include <TF1.h>
+
+#include "AliLnEfficiency.h"
+#include "B2.h"
+
+ClassImp(AliLnEfficiency)
+
+AliLnEfficiency::AliLnEfficiency(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag)
+: TObject()
+, fParticle(particle)
+, fSimuFilename(simuFilename)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+, fG3Fluka(0)
+{
+//
+// constructor
+//
+}
+
+AliLnEfficiency::~AliLnEfficiency()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnEfficiency::Exec()
+{
+//
+// extract the efficiency for the given simulation
+//
+       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());
+       }
+       
+       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* hEffTrigPt   = Divide(hGenTrigPt, hGenPhSpPt, fParticle + "_Eff_Trig_Pt");
+       TH1D* hEffVtxPt    = Divide(hGenVtxPt, hGenTrigPt, fParticle + "_Eff_Vtx_Pt");
+       TH1D* hEffAccPt    = Divide(hGenAccPt, hGenVtxPt, fParticle + "_Eff_Acc_Pt");
+       TH1D* hEffTrkPt    = Divide(hPrimRecPt, hGenAccPt, fParticle + "_Eff_Trk_Pt");
+       TH1D* hEffAccTrkPt = Divide(hPrimRecPt, hGenVtxPt, fParticle + "_Eff_AccTrk_Pt");
+       
+       if(fG3Fluka)
+       {
+               TF1* fncG3Fluka = 0;
+               if(fParticle == "Proton") fncG3Fluka = GetProtonG3FlukaCor("Proton_G3FLUKA_Corr_Pt");
+               if(fParticle == "AntiProton") fncG3Fluka = GetAntiProtonG3FlukaCor("AntiProton_G3FLUKA_Corr_Pt");
+               
+               if(fncG3Fluka != 0)
+               {
+                       hEffTrkPt->Divide(fncG3Fluka);
+                       hEffAccTrkPt->Divide(fncG3Fluka);
+               }
+               
+               delete fncG3Fluka;
+       }
+       
+       hEffAccTrkPt->SetTitle(fParticle.Data());
+       hEffTrkPt->SetTitle(fParticle.Data());
+       hEffAccPt->SetTitle(fParticle.Data());
+       hEffVtxPt->SetTitle(fParticle.Data());
+       hEffTrigPt->SetTitle(fParticle.Data());
+       
+       hEffAccTrkPt->SetYTitle("Acceptance #times Track Efficiency");
+       hEffAccPt->SetYTitle("Acceptance Efficiency");
+       hEffTrkPt->SetYTitle("Track Efficiency");
+       hEffVtxPt->SetYTitle("Vertex Efficiency");
+       hEffTrigPt->SetYTitle("Trigger Efficiency");
+       
+       hEffAccTrkPt->Write();
+       hEffAccPt->Write();
+       hEffTrkPt->Write();
+       hEffVtxPt->Write();
+       hEffTrigPt->Write();
+       
+       delete hEffAccTrkPt;
+       delete hEffAccPt;
+       delete hEffTrkPt;
+       delete hEffVtxPt;
+       delete hEffTrigPt;
+       
+       delete foutput;
+       delete fsimu;
+       
+       return 0;
+}
+
+TF1* AliLnEfficiency::GetProtonG3FlukaCor(const TString& name) const
+{
+//
+// GEANT3/FLUKA correction for protons (Baryon Twiki)
+// (TPC tracks)
+//
+       TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]", 0.001, 10);
+       fnc->SetParameters(4113, -32.9, -0.009383);
+       Double_t err[] = {1640.4, 1.5, 0.000917};
+       fnc->SetParErrors(err);
+       
+       return fnc;
+}
+
+TF1* AliLnEfficiency::GetAntiProtonG3FlukaCor(const TString& name) const
+{
+//
+// GEANT3/FLUKA correction for antiprotons (Baryon Twiki)
+// (TPC tracks)
+//
+       TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]+[3]*log(x)/pow(x,0.2)", 0.001, 10);
+       fnc->SetParameters(177.4, -22.03, -0.06538, 0.04431 );
+       Double_t err[] = {73., 1.5, 0.00221, 0.00555};
+       fnc->SetParErrors(err);
+       
+       return fnc;
+}
+
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.h
new file mode 100644 (file)
index 0000000..77e0713
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef ALILNEFFICIENCY_H
+#define ALILNEFFICIENCY_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// efficiency correction
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class TString;
+class TF1;
+
+class AliLnEfficiency: public TObject
+{
+  public:
+       
+       AliLnEfficiency(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag);
+       virtual ~AliLnEfficiency();
+       
+       TF1* GetProtonG3FlukaCor(const TString& name) const;
+       TF1* GetAntiProtonG3FlukaCor(const TString& name) const;
+       
+       const TString* GetOutputFilename() const { return &fOutputFilename; }
+       
+       Int_t Exec();
+       
+       void SetParticle(const TString& particle) { fParticle = particle; }
+       
+       void SetOutputTag(const TString& tag) { fOutputTag = tag; }
+       
+       void SetG3Fluka(Bool_t flag=1) { fG3Fluka = flag; }
+       
+  private:
+       AliLnEfficiency(const AliLnEfficiency& other);
+       AliLnEfficiency& operator=(const AliLnEfficiency& other);
+       
+  private:
+       TString fParticle; // particle name
+       
+       TString fSimuFilename; // simulation filename
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for the ouput file
+       
+       Bool_t fG3Fluka; // use G3/FLUKA correction for (anti)protons
+       
+       ClassDef(AliLnEfficiency,1)
+};
+
+#endif // ALILNEFFICIENCY_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.cxx
new file mode 100644 (file)
index 0000000..8ecd6e3
--- /dev/null
@@ -0,0 +1,82 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// fake track fraction
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TString.h>
+#include <TF1.h>
+
+#include "AliLnFakeTracks.h"
+#include "B2.h"
+
+ClassImp(AliLnFakeTracks)
+
+AliLnFakeTracks::AliLnFakeTracks(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag)
+: TObject()
+, fParticle(particle)
+, fSimuFilename(simuFilename)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+{
+//
+// constructor
+//
+}
+
+AliLnFakeTracks::~AliLnFakeTracks()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnFakeTracks::Exec()
+{
+//
+// extract the fraction of fake tracks for the given fParticle.Data()
+//
+       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());
+       }
+       
+       TH1D* hFakeTracks = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Fake_Pt");
+       TH1D* hAllTracks  = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Pt");
+       
+       TH1D* hFracFakePt = Divide(hFakeTracks, hAllTracks, fParticle + "_Frac_Fake_Pt");
+       
+       hFracFakePt->SetYTitle("Fake Tracks/All Tracks");
+       hFracFakePt->Write();
+       
+       delete hFracFakePt;
+       
+       delete foutput;
+       delete fsimu;
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnFakeTracks.h
new file mode 100644 (file)
index 0000000..d75b833
--- /dev/null
@@ -0,0 +1,47 @@
+#ifndef ALILNFAKETRACKS_H
+#define ALILNFAKETRACKS_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// fake track fraction
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class TString;
+
+class AliLnFakeTracks: public TObject
+{
+  public:
+       
+       AliLnFakeTracks(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag);
+       virtual ~AliLnFakeTracks();
+       
+       const TString* GetOutputFilename() const { return &fOutputFilename; }
+       
+       Int_t Exec();
+       
+       void SetParticle(const TString& particle) { fParticle = particle; }
+       
+       void SetOutputTag(const TString& tag) { fOutputTag = tag; }
+       
+  private:
+       AliLnFakeTracks(const AliLnFakeTracks& other);
+       AliLnFakeTracks& operator=(const AliLnFakeTracks& other);
+       
+  private:
+       TString fParticle; // particle
+       
+       TString fSimuFilename; // simulation filename
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for the ouput
+       
+       
+       ClassDef(AliLnFakeTracks,1)
+};
+
+#endif // ALILNFAKETRACKS_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx
new file mode 100644 (file)
index 0000000..6a07e3d
--- /dev/null
@@ -0,0 +1,566 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// reconstructed pt
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TF1.h>
+#include <TString.h>
+#include <TROOT.h>
+#include <TString.h>
+#include <RooWorkspace.h>
+#include <RooRealVar.h>
+#include <RooDataSet.h>
+#include <RooDataHist.h>
+#include <RooHistPdf.h>
+#include <RooMsgService.h>
+
+#include "AliLnPt.h"
+#include "B2.h"
+
+#ifdef HAVE_ROOUNFOLD
+  #include "/opt/RooUnfold/src/RooUnfoldResponse.h"
+  #include "/opt/RooUnfold/src/RooUnfoldBayes.h"
+#endif
+
+ClassImp(AliLnPt)
+
+AliLnPt::AliLnPt(const TString& particle, Double_t trigEff, const TString& inputFilename, const TString& outputFilename, const TString& otag, const TString& corrFilename, const TString& corrtag)
+: TObject()
+, fParticle(particle)
+, fTrigEff(trigEff)
+, fInputFilename(inputFilename)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+, fCorrFilename(corrFilename)
+, fCorrTag(corrtag)
+, fLowPtBin(1)
+, fHiPtBin(19)
+, fPidM2(0)
+, fUnfolding(0)
+, fNIter(4)
+, fFakeTracks(1)
+, 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.5)
+, fMakeStats(1)
+, fVtxCorr(0)
+, fVtxFactor(1)
+, fFitFrac(1)
+, fSameFdwn(0)
+{
+//
+// constructor
+//
+       TH1::SetDefaultSumw2(); // switch on histogram errors
+       
+       // disable verbose in RooFit
+       RooMsgService::instance().setSilentMode(1);
+       RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+}
+
+AliLnPt::~AliLnPt()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnPt::Exec()
+{
+//
+// Get the true pt distribution
+//
+       using namespace std;
+       
+       TFile* finput = new TFile(fInputFilename.Data(), "read");
+       if (finput->IsZombie()) exit(1);
+       
+       TFile* foutput = new TFile(fOutputFilename.Data(), "recreate");
+       if (foutput->IsZombie()) exit(1);
+       
+       TFile* fcorr = 0;
+       TFile* fdebug = 0;
+       
+       TString species = fParticle;
+       species.ReplaceAll("Anti","");
+       
+       TH1D* hStats = (TH1D*)FindObj(finput, species + "_Stats");
+       
+       TH1D* hCorPt = 0;
+       
+       if(fIsOnlyGen)
+       {
+               hCorPt = (TH1D*)((TH1D*)FindObj(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data()));
+       }
+       else // correct the reconstructed pt
+       {
+               fcorr = new TFile(fCorrFilename.Data(),"read");
+               if (fcorr->IsZombie()) exit(1);
+               
+               fdebug =  new TFile(Form("debug-%s",fOutputFilename.Data()),"recreate");
+               if (fdebug->IsZombie()) exit(1);
+               
+               fdebug->mkdir(fOutputTag.Data());
+               
+               // -------------- get pointers to the required histograms ---------
+               
+               // pt distributions
+               TH1D* hOrigPt = (TH1D*)FindObj(finput, fParticle + "_PID_Pt");
+               TH2D* hM2Pt   = (TH2D*)FindObj(finput, fParticle + "_PID_M2_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
+               
+               TH1D* hFracFakePt  = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fake_Pt");
+               
+               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");
+               
+               if(fSameFdwn)
+               {
+                       hFracFdwnPt = (TH1D*)FindObj(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
+                       fncFracFdwnPt = (TF1*)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");
+               
+               // -------------- apply corrections -----------
+               
+               fdebug->cd(fOutputTag.Data());
+               
+               TH1D* hPidPt = 0;
+               if(fPidM2)
+               {
+                       hPidPt = this->PID(hOrigPt, hM2Pt, fLowM2Bin, fHiM2Bin, fParticle + "_M2Corr_Pt");
+               }
+               else
+               {
+                       hPidPt = (TH1D*)hOrigPt->Clone(Form("%s_M2Corr_Pt", fParticle.Data()));
+               }
+               
+               TH1D* hFakeCorPt = 0;
+               if(fFakeTracks)
+               {
+                       hFakeCorPt = this->FakeTracks(hPidPt, hFracFakePt, fParticle + "_FakeCor_Pt");
+               }
+               else
+               {
+                       hFakeCorPt = (TH1D*)hPidPt->Clone(Form("%s_FakeCor_Pt", fParticle.Data()));
+               }
+               
+               TH1D* hSecCorPt = 0;
+               if(fSecondaries)
+               {
+                       if(fFitFrac)
+                       {
+                               hSecCorPt = this->Secondaries(hFakeCorPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCor_Pt");
+                       }
+                       else
+                       {
+                               hSecCorPt = this->Secondaries(hFakeCorPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCor_Pt");
+                       }
+               }
+               else
+               {
+                       hSecCorPt = (TH1D*)hFakeCorPt->Clone(Form("%s_SecCor_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;
+               if(fEfficiency)
+               {
+                       hEffCorPt = this->Efficiency(hUnfoldedPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCor_Pt");
+               }
+               else
+               {
+                       hEffCorPt = (TH1D*)hUnfoldedPt->Clone(Form("%s_EffCor_Pt", fParticle.Data()));
+               }
+               
+               // --------- corrected pt --------
+               
+               hCorPt = (TH1D*) hEffCorPt->Clone(Form("%s_Cor_Pt", fParticle.Data()));
+               
+               // ---------- debug -----
+               
+               fdebug->cd(fOutputTag.Data());
+               
+               hOrigPt->Write(); // original distribution
+               hPidPt->Write();
+               hUnfoldedPt->Write();
+               hFakeCorPt->Write();
+               hSecCorPt->Write();
+               hEffCorPt->Write();
+               hCorPt->Write();
+               
+               // ------------- clean --------------
+               
+               delete hPidPt;
+               delete hUnfoldedPt;
+               delete hFakeCorPt;
+               delete hSecCorPt;
+               delete hEffCorPt;
+       }
+       
+       // ---------- save the final pt ----------
+       
+       TH1D* hPt = (TH1D*) hCorPt->Clone(Form("%s_Pt", fParticle.Data()));
+       hPt->SetTitle(fParticle.Data());
+       hPt->Reset();
+       
+       for(Int_t i=fLowPtBin; i<fHiPtBin; ++i)
+       {
+               hPt->SetBinContent(i, hCorPt->GetBinContent(i));
+               hPt->SetBinError(i, hCorPt->GetBinError(i));
+       }
+       
+       foutput->cd();
+       
+       hPt->Write();
+       
+       delete hPt;
+       delete hCorPt;
+       
+       delete fdebug;
+       delete fcorr;
+       
+       // --------- stats -------------------
+       
+       if(fMakeStats)
+       {
+               TH1D* hStatsFix = new TH1D(Form("%s_Stats", fParticle.Data()), "Stats", 6, 0, 6);
+               hStatsFix->GetXaxis()->SetBinLabel(1,"Events");
+               hStatsFix->GetXaxis()->SetBinLabel(2,"Trig");
+               hStatsFix->GetXaxis()->SetBinLabel(3,"Ana");
+               hStatsFix->GetXaxis()->SetBinLabel(4,"Inel");
+               hStatsFix->GetXaxis()->SetBinLabel(5,"Vtx");
+               hStatsFix->GetXaxis()->SetBinLabel(6,"Vz");
+               hStatsFix->SetStats(0);
+               
+               // Fill
+               hStatsFix->Fill("Events",hStats->Integral(1,1));
+               hStatsFix->Fill("Trig",hStats->Integral(2,2));
+               if(fIsOnlyGen)
+               {
+                       hStatsFix->Fill("Ana",hStats->Integral(2,2));
+                       hStatsFix->Fill("Inel",hStats->Integral(2,2));
+               }
+               else
+               {
+                       Double_t nTrig   = hStats->Integral(2,2);
+                       Double_t nAna    = hStats->Integral(3,3);
+                       Double_t nVtx    = hStats->Integral(4,4);
+                       Double_t nVz     = hStats->Integral(5,5);
+                       Double_t nInel   = nTrig/fTrigEff;
+                       
+                       if(!fVtxCorr)
+                       {
+                               // Add events without vertex and not triggered
+                               //Double_t nNoVtx = (nVz/nVtx)*(nTrig-nVtx);
+                               Double_t nNoVtx = (nAna/nVtx)*(nTrig-nVtx);
+                               nInel = (nAna + nNoVtx)/fTrigEff;
+                               nAna  = nAna + nNoVtx;
+                               // Nvz/Ninel * (Ninel - Nvtx)
+                       }
+                       
+                       hStatsFix->Fill("Ana",nAna);
+                       hStatsFix->Fill("Inel",nInel);
+                       hStatsFix->Fill("Vtx",nVtx);
+                       hStatsFix->Fill("Vz",nVz);
+               }
+               
+               foutput->cd();
+               hStatsFix->Write();
+               delete hStatsFix;
+       }
+       
+       delete foutput;
+       delete finput;
+       
+       return 0;
+}
+
+Double_t AliLnPt::GetM2Width(Double_t pt, Double_t m) const
+{
+//
+// parameterization of the expected M2 width
+// as a function of pt from LHC12a5b
+//
+       Double_t c0 = 1.90212e-03;
+       Double_t c1 = 1.29814e-02;
+       
+       TString species = fParticle;
+       species.ReplaceAll("Anti","");
+       
+       if(species == "Deuteron")
+       {
+               c0 = 1.37785e-02;
+               c1 = 6.78951e-02;
+       }
+       else if(species == "Triton")
+       {
+               c0 = 3.79471e-02;
+               c1 = 1.87478e-01;
+       }
+       else if(species == "He3")
+       {
+               c0 = 2.05887e-01;
+               c1 = 1.17813e-01;
+       }
+       else if(species == "Alpha")
+       {
+               c0 = 4.13329e-01;
+               c1 = 2.27076e-01;
+       }
+       
+       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
+{
+//
+// deuteron m2 model for TPC pid
+// gaussian + exponential background of pi, K and p
+//
+       using namespace RooFit;
+       
+       RooWorkspace* w = new RooWorkspace(name.Data());
+       
+       RooRealVar m2("m2", "m^{2} (GeV^{2}/c^{4})", fMinM2Bkg, fMaxM2Bkg);
+       w->import(m2);
+       
+       Double_t sigma = GetM2Width(pt, m);
+       
+       // 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));
+       
+       // 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("SUM::model( Nd[0,%f]*Sd, Nbkg[0,%f]*Bkg )",max,max));
+       
+       return w;
+}
+
+void AliLnPt::GetPtFromM2(TH1D* hPt, Int_t lowbin, Int_t hibin, const TH2D* hM2Pt, Double_t m2min, Double_t m2max) const
+{
+//
+// Fill the pt from lowbin to hibin integrating
+// the m2-pt distribution from m2min to m2max
+//
+       for(Int_t i=lowbin; i<hibin; ++i)
+       {
+               TH1D* hM2 = hM2Pt->ProjectionY(Form("_M2_%02d",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
+               
+               Double_t err = 0;
+               Double_t n = hM2->IntegralAndError(m2minbin, m2maxbin, err);
+               
+               hPt->SetBinContent(i, n);
+               hPt->SetBinError(i, err);
+               
+               delete hM2;
+       }
+}
+
+TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hM2Pt, Int_t lowbin, Int_t hibin, const TString& name)
+{
+//
+// Remove contamination in the m^2 distribution
+//
+       const Int_t kNBin = 2;
+       Double_t m = GetMass(fParticle);
+       
+       TH1D* hPidPt = (TH1D*) hPt->Clone(name.Data());
+       hPidPt->Reset();
+       
+       // first bins are not contaminated since the identification is clear
+       // integrate around the m2 value
+       
+       this->GetPtFromM2(hPidPt, 3, lowbin, hM2Pt, fMinM2tpc, fMaxM2tpc);
+       
+       // the remaining bins are contaminated
+       // slice the m2-pt and fit each slice with a gaussian + exponentials bkg
+       
+       using namespace RooFit;
+       
+       for(Int_t i=lowbin; i<hibin; ++i)
+       {
+               TH1D* hM2 = hM2Pt->ProjectionY(Form("%s_PID_M2_%02d",fParticle.Data(),i),i,i);
+               hM2->Rebin(kNBin);
+               
+               RooWorkspace* w = GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hM2->Integral());
+               
+               RooRealVar* m2 = w->var("m2");
+               RooDataHist data("data", "data to fit", *m2, hM2);
+               
+               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());
+               
+               w->Write();
+               
+               delete w;
+               delete hM2;
+       }
+       
+       return hPidPt;
+}
+
+#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::FakeTracks(const TH1D* hPt, const TH1D* hFracFakePt, const TString& name) const
+{
+//
+// remove fake tracks
+// N_true = N - N_fake = N - f_fake*N
+//
+       TH1D* hFakeTracks = (TH1D*)hPt->Clone("_fake_tracks_");
+       hFakeTracks->Multiply(hFracFakePt);
+       
+       TH1D* hFakeCorPt = (TH1D*)hPt->Clone(name.Data());
+       hFakeCorPt->Add(hFakeTracks,-1.);
+       
+       delete hFakeTracks;
+       
+       return hFakeCorPt;
+}
+
+TH1D* AliLnPt::Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* hFracFdwnPt, const TString& name) const
+{
+//
+// remove contamination of secondaries using fractions
+// N_prim = N/(1+fmat+fdwn)
+//
+       TF1* one = new TF1("one","1",0,10);
+       
+       TH1D* xfactor = (TH1D*)hFracFdwnPt->Clone(Form("%s_fdwn_fmat_1_",fParticle.Data()));
+       xfactor->Add(hFracMatPt);
+       xfactor->Add(one);
+       
+       TH1D* hSecCorPt = (TH1D*)hPt->Clone(name.Data());
+       
+       hSecCorPt->Divide(xfactor);
+       
+       delete xfactor;
+       delete one;
+       
+       return hSecCorPt;
+}
+
+TH1D* AliLnPt::Secondaries(const TH1D* hPt, TF1* fncMatPt, TF1* fncFdwnPt, const TString& name) const
+{
+//
+// remove contamination of secondaries using fitted 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);
+       
+       Double_t param[6];
+       fncMatPt->GetParameters(param);
+       fncFdwnPt->GetParameters(&param[3]);
+       
+       Double_t err[6];
+       for(Int_t i=0; i<3; ++i) err[i] = fncMatPt->GetParError(i);
+       for(Int_t i=3; i<6; ++i) err[i] = fncFdwnPt->GetParError(i);
+       
+       fnc->SetParameters(param);
+       fnc->SetParErrors(err);
+       
+       TH1D* hSecCorPt = (TH1D*)hPt->Clone(name.Data());
+       hSecCorPt->Divide(fnc);
+       
+       delete fnc;
+       
+       return hSecCorPt;
+}
+
+TH1D* AliLnPt::Efficiency(const TH1D* hPt, const TH1D* hEffVtxPt, const TH1D* hEffAccTrkPt, const TString& name) const
+{
+//
+// correct by efficiency 1/(eff_vtx x eff_acc x eff_trk)
+//
+       TH1D* hEffCorPt = (TH1D*)hPt->Clone(name.Data());
+       hEffCorPt->Divide(hEffAccTrkPt);
+       
+       if(fVtxCorr)
+       {
+               hEffCorPt->Divide(hEffVtxPt);
+               hEffCorPt->Scale(fVtxFactor);
+       }
+       
+       return hEffCorPt;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.h
new file mode 100644 (file)
index 0000000..299620e
--- /dev/null
@@ -0,0 +1,126 @@
+#ifndef ALILNPT_H
+#define ALILNPT_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// reconstructed pt
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+#ifndef HAVE_ROOUNFOLD
+//#define HAVE_ROOUNFOLD
+#endif
+
+class TString;
+class TH1D;
+class TH2D;
+class RooWorkspace;
+class TF1;
+
+class AliLnPt: public TObject
+{
+  public:
+       
+       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* FakeTracks(const TH1D* hPt, const TH1D* hFracFakePt, const TString& name) const;
+       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;
+       
+       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 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 SetFakeTracks(Bool_t flag=1) { fFakeTracks = flag; }
+       void SetSecondaries(Bool_t flag=1) { fSecondaries = flag; }
+       void SetEfficiency(Bool_t flag=1) { fEfficiency = flag; }
+       
+       void SetOnlyGeneration(Bool_t flag=1) { fIsOnlyGen = flag; }
+       
+       void SetMakeStats(Bool_t flag=1) { fMakeStats = flag; }
+       
+       void SetVertexCorrection(Bool_t flag=1, Double_t factor=1) { fVtxCorr=flag; fVtxFactor=factor; }
+       void SetFitFractionCorr(Bool_t flag=1) { fFitFrac=flag; }
+       void SetSameFeedDownCorr(Bool_t flag=1) { fSameFdwn=flag; }
+       
+  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;
+       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;
+       
+  private:
+       
+       TString fParticle; // particle name
+       
+       Double_t fTrigEff; // trigger efficiency
+       
+       TString fInputFilename; // input filename
+       
+       TString fOutputFilename; // output for the result
+       TString fOutputTag; // tag for the ouput
+       
+       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
+       
+       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 fFakeTracks; // remove faketracks
+       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
+       
+       Int_t fLowM2Bin; // low pt bin for pid correction
+       Int_t fHiM2Bin; // high pt bin for pid correction
+       
+       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
+       
+       Bool_t fMakeStats; // make event stats
+       
+       Bool_t fVtxCorr; // use simulation to correct from non reconstructed vertex events
+       Double_t fVtxFactor; // correction factor between simulation and data
+       
+       Bool_t fFitFrac; // fit for fraction of secondaries
+       Bool_t fSameFdwn; // same feed-down correction for positives and negatives
+       
+       ClassDef(AliLnPt,1)
+};
+
+#endif // ALILNPT_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx
new file mode 100644 (file)
index 0000000..da37f36
--- /dev/null
@@ -0,0 +1,84 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// antiparticle / particle ratio
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TROOT.h>
+
+#include "AliLnRatio.h"
+#include "B2.h"
+
+ClassImp(AliLnRatio)
+
+AliLnRatio::AliLnRatio(const TString& species, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag)
+: TObject()
+, fSpecies(species)
+, fPtFilename(ptFilename)
+, fPtTag(tag)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+{
+//
+// constructor
+//
+}
+
+AliLnRatio::~AliLnRatio()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnRatio::Exec()
+{
+//
+// antiparticle/particle ratio
+//
+       using namespace std;
+       
+       TFile* finput = new TFile(fPtFilename.Data(), "read");
+       if (finput->IsZombie()) exit(1);
+       
+       const TString kPrefix[] = { "", "Anti" };
+       
+       TH1D* hPt[2];
+       
+       for(Int_t i=0; i<2; ++i)
+       {
+               hPt[i] = (TH1D*)FindObj(finput, kPrefix[i] + fSpecies + "_Pt");
+       }
+       
+       TH1D* hRatioPt = Divide(hPt[1],hPt[0],Form("Anti%s%s_Ratio_Pt",fSpecies.Data(),fSpecies.Data()));
+       
+       TFile* foutput = new TFile(fOutputFilename.Data(), "recreate");
+       
+       foutput->mkdir(fOutputTag.Data());
+       foutput->cd(fOutputTag.Data());
+       
+       hRatioPt->Write();
+       
+       delete hRatioPt;
+       delete foutput;
+       delete finput;
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.h
new file mode 100644 (file)
index 0000000..ad5e4cc
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef ALILNRATIO_H
+#define ALILNRATIO_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// antiparticle / particle ratio
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+
+class AliLnRatio: public TObject
+{
+  public:
+       
+       AliLnRatio(const TString& species, const TString& ptFilename, const TString& ptTag, const TString& outputFilename, const TString& otag);
+       
+       virtual ~AliLnRatio();
+       
+       Int_t Exec();
+       
+  private:
+       AliLnRatio(const AliLnRatio& other);
+       AliLnRatio& operator=(const AliLnRatio& other);
+       
+  private:
+       TString fSpecies; // particle species
+       TString fPtFilename; // pt filename
+       TString fPtTag; // tag for pt
+       
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for the output
+       
+       ClassDef(AliLnRatio,1)
+};
+
+#endif // ALILNRATIO_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx
new file mode 100644 (file)
index 0000000..d783e05
--- /dev/null
@@ -0,0 +1,700 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// removal of secondaries using DCA distributions
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TFractionFitter.h>
+#include <TObjArray.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <RooRealVar.h>
+#include <RooDataSet.h>
+#include <RooDataHist.h>
+#include <RooHistPdf.h>
+#include <RooWorkspace.h>
+#include <RooMsgService.h>
+
+#include "AliLnSecondaries.h"
+#include "B2.h"
+
+ClassImp(AliLnSecondaries)
+
+AliLnSecondaries::AliLnSecondaries(const TString& particle, const TString& dataFilename, const TString& simuFilename, const TString& outputFilename, const TString& otag)
+: TObject()
+, fParticle(particle)
+, fDataFilename(dataFilename)
+, fSimuFilename(simuFilename)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+, fLowBin(3)
+, fHiBin(15)
+, fNbin(1)
+, fMinDCAxy(-1.5)
+, fMaxDCAxy(1.5)
+, fFracProc(0)
+, fMatDCAxyMod(AliLnSecondaries::kFlatDCAxy)
+, fScMat(1)
+, fScFd(1)
+{
+//
+// constructor
+//
+       TH1::SetDefaultSumw2(); // switch on histogram errors
+       
+       // disable verbose in RooFit
+       RooMsgService::instance().setSilentMode(1);
+       RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+}
+
+AliLnSecondaries::~AliLnSecondaries()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnSecondaries::Exec()
+{
+//
+// Fractions of primaries and secondaries using templates from simulations
+// Secondaries are from feed-down and/or from materials
+//
+       using namespace std;
+       
+       TFile* fdata = new TFile(fDataFilename.Data(),"read");
+       if(fdata->IsZombie()) exit(1);
+       
+       TFile* fsimu = new TFile(fSimuFilename.Data(),"read");
+       if(fsimu->IsZombie()) exit(1);
+       
+       TFile* fdebug =  new TFile(Form("debug-%s",fOutputFilename.Data()),"recreate");
+       
+       if(fOutputTag != "") fdebug->mkdir(fOutputTag.Data());
+       if(fOutputTag != "") fdebug->cd(fOutputTag.Data());
+       
+       // --------- ideal values: primaries + no secondaries --------
+       
+       TH1D* hPidPt = (TH1D*)FindObj(fdata, fParticle + "_PID_Pt");
+       hPidPt->Write();
+       
+       const char* contrib[] = { "prim", "mat", "fdwn" };
+       
+       TH1D* hFracPt[3];
+       for(Int_t i=0; i<3; ++i)
+       {
+               hFracPt[i] = this->ZeroClone(hPidPt, Form("%s_P%s_Pt",fParticle.Data(),contrib[i]));
+               hFracPt[i]->SetYTitle(Form("P_{%s}",contrib[i]));
+       }
+       
+       // fill
+       
+       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");
+               
+               delete hFracPt[0];
+               delete hFracPt[1];
+               delete hFracPt[2];
+               
+               hFracPt[0] = Divide(hPrim, hAll, fParticle + "_Pprim_Pt");
+               hFracPt[1] = Divide(hMat,  hAll, fParticle + "_Pmat_Pt");
+               hFracPt[2] = Divide(hFdwn, hAll, fParticle + "_Pfdwn_Pt");
+               
+               hFracPt[1]->Scale(fScMat);
+               hFracPt[2]->Scale(fScFd);
+       }
+       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");
+               
+               hDataDCAxyPt->Rebin2D(1,fNbin);
+               hDataMCDCAxyPt->Rebin2D(1,fNbin);
+               hPrimDCAxyPt->Rebin2D(1,fNbin);
+               hFdwnDCAxyPt->Rebin2D(1,fNbin);
+               
+               this->GetFraction(hFracPt[0], hFracPt[2], hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hFdwnDCAxyPt, "Fdwn");
+       }
+       else if(fParticle.Contains("Anti"))
+       {
+               // assume there are no secondaries for antinuclei
+               this->GetFraction(hFracPt[0]);
+       }
+       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");
+               
+               if(fMatDCAxyMod == kFlatDCAxy)
+               {
+                       hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
+               }
+               
+               hDataDCAxyPt->Rebin2D(1,fNbin);
+               hDataMCDCAxyPt->Rebin2D(1,fNbin);
+               hPrimDCAxyPt->Rebin2D(1,fNbin);
+               hMatDCAxyPt->Rebin2D(1,fNbin);
+               hFdwnDCAxyPt->Rebin2D(1,fNbin);
+               
+               this->GetFraction(hFracPt, hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hMatDCAxyPt, hFdwnDCAxyPt);
+       }
+       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"); // for debugging
+               TH2D* hPrimDCAxyPt   = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
+               // use the anti-ion as template for primaries
+               //TH2D* hPrimDCAxyPt   = (TH2D*)FindObj(fdata, Form("Anti%s_PID_DCAxy_Pt",fParticle.Data()));
+               TH2D* hMatDCAxyPt    = (TH2D*)FindObj(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
+               
+               if(fMatDCAxyMod == kFlatDCAxy)
+               {
+                       hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
+               }
+               
+               hDataDCAxyPt->Rebin2D(1,fNbin);
+               hDataMCDCAxyPt->Rebin2D(1,fNbin);
+               hPrimDCAxyPt->Rebin2D(1,fNbin);
+               hMatDCAxyPt->Rebin2D(1,fNbin);
+               
+               this->GetFraction(hFracPt[0], hFracPt[1], hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hMatDCAxyPt, "Mat");
+       }
+       
+       // --------- save --------------
+       
+       TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
+       if(fOutputTag != "")
+       {
+               foutput->mkdir(fOutputTag.Data());
+               foutput->cd(fOutputTag.Data());
+       }
+       
+       // ----- fractions w.r.t. primaries ------------
+       
+       TH1D* hMatFracPt = (TH1D*)hFracPt[1]->Clone(Form("%s_Frac_Mat_Pt", fParticle.Data()));
+       hMatFracPt->SetYTitle("P_{1} / P_{0}");
+       hMatFracPt->Divide(hFracPt[0]);
+       hMatFracPt->Write();
+       
+       TH1D* hFdwnFracPt = (TH1D*)hFracPt[2]->Clone( Form("%s_Frac_Fdwn_Pt",fParticle.Data()));
+       hFdwnFracPt->SetYTitle("P_{2} / P_{0}");
+       hFdwnFracPt->Divide(hFracPt[0]);
+       hFdwnFracPt->Write();
+       
+       // --------- optionally fit fractions to extrapolate at high pt -----------
+       
+       TF1* fncMat = this->GetMatFraction(fParticle + "_Frac_Mat_Fit_Pt");
+       if(fParticle.Contains("Anti"))
+       {
+               fncMat->SetParameters(0,0,0);
+       }
+       else if(fParticle=="Proton")
+       {
+               hMatFracPt->Fit(fncMat, "NQ", "", 0.,2);
+       }
+       else
+       {
+               hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.4);
+       }
+       
+       TF1* fncFdwn = this->GetFdwnFraction(fParticle + "_Frac_Fdwn_Fit_Pt");
+       if(!fParticle.Contains("Proton"))
+       {
+               fncFdwn->SetParameters(0,0,0);
+       }
+       else
+       {
+               hFdwnFracPt->Fit(fncFdwn, "NQ", "", 0.5,3.);
+       }
+       
+       fncFdwn->Write();
+       fncMat->Write();
+       
+       // ---------- clean ---------
+       
+       delete fncFdwn;
+       delete fncMat;
+       
+       delete hMatFracPt;
+       delete hFdwnFracPt;
+       
+       for(Int_t i=0; i<3; ++i)
+       {
+               hFracPt[i]->Write();
+               delete hFracPt[i];
+       }
+       
+       delete foutput;
+       delete fdebug;
+       delete fsimu;
+       delete fdata;
+       
+       return 0;
+}
+
+void AliLnSecondaries::GetFraction(TH1D* hPrimPt) const
+{
+//
+// No secondaries only primaries
+//
+       TF1* one = new TF1("one", "1", hPrimPt->GetXaxis()->GetXmin(), hPrimPt->GetXaxis()->GetXmax());
+       hPrimPt->Reset();
+       hPrimPt->Add(one);
+       delete one;
+}
+
+void AliLnSecondaries::GetFraction(TH1D* hPrimPt, TH1D* hSecPt, const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hSecDCAxyPt, const TString& sec) const
+{
+//
+// slice the DCA distributions and compute the fractions
+// 2 contributions (primaries and secondaries)
+//
+       TString nosec = (sec == "Fdwn") ? "Mat" : "Fdwn";
+       
+       for(Int_t i=fLowBin; i<fHiBin; ++i)
+       {
+               TH1D* hDCAxy      = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hMCDCAxy    = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hPrimDCAxy  = hPrimDCAxyPt->ProjectionY(Form("%s_Prim_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hSecDCAxy   = hSecDCAxyPt->ProjectionY(Form("%s_%s_DCAxy_%02d",fParticle.Data(),sec.Data(),i),i,i);
+               
+               // fractions
+               Double_t frac[2] = {0,0};
+               Double_t err[2]  = {0,0};
+               
+               if(fFracProc == kTFractionFitter)
+               {
+                       this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
+               }
+               else
+               {
+                       this->GetRooFitFractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
+               }
+               
+               // fill
+               TH1D* hFracPt[] = { hPrimPt, hSecPt};
+               for(Int_t j=0; j<2; ++j)
+               {
+                       hFracPt[j]->SetBinContent(i, frac[j]);
+                       hFracPt[j]->SetBinError(i, err[j]);
+               }
+               
+               // ------------- debug histograms ------------
+               
+               hDCAxy->Write();
+               hMCDCAxy->Write();
+               hPrimDCAxy->Write();
+               hSecDCAxy->Write();
+               
+               TH1D* hNoSecDCAxy = this->ZeroClone(hMCDCAxy, Form("%s_%s_DCAxy_%02d",fParticle.Data(),nosec.Data(),i));
+               TH1D* hNoSecFit   = this->ZeroClone(hMCDCAxy, Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),nosec.Data(),i));
+               
+               hNoSecDCAxy->Write();
+               hNoSecFit->Write();
+               
+               delete hNoSecDCAxy;
+               delete hNoSecFit;
+               
+               // ----------------------- end debug ----------------------
+               
+               delete hDCAxy;
+               delete hMCDCAxy;
+               delete hPrimDCAxy;
+               delete hSecDCAxy;
+       }
+}
+
+void AliLnSecondaries::GetFraction(TH1D* hFracPt[3], const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hMatDCAxyPt, const TH2D* hFdwnDCAxyPt) const
+{
+//
+// slice the DCA distribution and get the fractions for each pt bin
+// (3 contributions)
+//
+       for(Int_t i=fLowBin; i<fHiBin; ++i)
+       {
+               // slices
+               TH1D* hDCAxy     = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hMCDCAxy   = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hPrimDCAxy = hPrimDCAxyPt->ProjectionY(Form("%s_Prim_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hMatDCAxy  = hMatDCAxyPt->ProjectionY(Form("%s_Mat_DCAxy_%02d",fParticle.Data(),i),i,i);
+               TH1D* hFdwnDCAxy = hFdwnDCAxyPt->ProjectionY(Form("%s_Fdwn_DCAxy_%02d",fParticle.Data(),i),i,i);
+               
+               // fractions
+               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);
+               }
+               
+               // fill
+               for(Int_t k=0; k<3; ++k)
+               {
+                       hFracPt[k]->SetBinContent(i, frac[k]);
+                       hFracPt[k]->SetBinError(i, err[k]);
+               }
+               
+               // -------------- debug ---------------------------
+               
+               hDCAxy->Write();
+               hMCDCAxy->Write();
+               hPrimDCAxy->Write();
+               hMatDCAxy->Write();
+               hFdwnDCAxy->Write();
+               
+               // --------------------- end debug ----------------
+               
+               delete hDCAxy;
+               delete hMCDCAxy;
+               delete hPrimDCAxy;
+               delete hMatDCAxy;
+               delete hFdwnDCAxy;
+       }
+}
+
+Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hSec, Int_t ibin, const char* sec) const
+{
+//
+// find the fractions of 2 contributions
+//
+       TObjArray* mc = new TObjArray(2);
+       
+       mc->Add(hPrim);
+       mc->Add(hSec);
+       
+       TFractionFitter* fit = new TFractionFitter(hData, mc, "Q"); // initialise
+       fit->Constrain(1,0.0001,1.);
+       fit->Constrain(2,0.0001,1.);
+       
+       Int_t status = fit->Fit();
+       
+       if (status == 0) // get fractions
+       {
+               fit->GetResult(0, frac[0], err[0]);
+               fit->GetResult(1, frac[1], err[1]);
+       }
+       
+       const char* contrib[] = {"Prim", sec };
+       this->WriteTFFdebug(hData, fit, status, ibin, contrib, frac, 2);
+       
+       delete mc;
+       delete fit;
+       
+       return status;
+}
+
+Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hMat, TH1D* hFdwn, Int_t ibin) const
+{
+//
+// find the fractions of 3 contributions
+//
+       TObjArray* mc = new TObjArray(3);
+       
+       mc->Add(hPrim);
+       mc->Add(hMat);
+       mc->Add(hFdwn);
+       
+       TFractionFitter* fit = new TFractionFitter(hData, mc, "Q");
+       fit->Constrain(1,0.0001,1.);
+       fit->Constrain(2,0.0001,1.);
+       fit->Constrain(3,0.0001,1.);
+       //fit->SetRangeX(1,50);
+       
+       Int_t status = fit->Fit();
+       
+       if (status == 0) // get fractions
+       {
+               fit->GetResult(0, frac[0], err[0]);
+               fit->GetResult(1, frac[1], err[1]);
+               fit->GetResult(2, frac[2], err[2]);
+       }
+       
+       const char* contrib[] = {"Prim", "Mat", "Fdwn"};
+       this->WriteTFFdebug(hData, fit, status, ibin, contrib, frac, 3);
+       
+       delete mc;
+       delete fit;
+       
+       return status;
+}
+
+void AliLnSecondaries::WriteTFFdebug(TH1D* hData, TFractionFitter* fit, Int_t status, Int_t ibin, const char* contrib[], Double_t* frac, Int_t kmax) const
+{
+//
+// Write TFractionFitter debug histograms
+//
+       if (status == 0)
+       {
+               TH1D* hResult = (TH1D*) fit->GetPlot();
+               
+               hResult->SetName(Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
+               hResult->Write();
+               
+               for(Int_t k=0; k<kmax; ++k)
+               {
+                       TH1D* hMCpred = (TH1D*) fit->GetMCPrediction(k);
+                       hMCpred->SetName(Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),contrib[k],ibin));
+                       
+                       hMCpred->Sumw2();
+                       hMCpred->Scale(frac[k]*hResult->Integral()/hMCpred->Integral());
+                       
+                       hMCpred->Write();
+               }
+       }
+       else // write void histograms
+       {
+               TH1D* hZero = this->ZeroClone(hData,Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
+               hZero->Write();
+               
+               for(Int_t k=0; k<kmax; ++k)
+               {
+                       hZero->SetName(Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),contrib[k],ibin));
+                       hZero->Write();
+               }
+               delete hZero;
+       }
+}
+
+TH1D* AliLnSecondaries::ZeroClone(TH1D* h, const char* name) const
+{
+//
+// clone histogram and reset to zero
+//
+       TH1D* clone = (TH1D*)h->Clone(name);
+       clone->Reset();
+       
+       return clone;
+}
+
+void AliLnSecondaries::GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hSec, Int_t ibin, const char* 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 =  (TH1D*)hData->Clone(Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
+       TH1D* hDebugPrim = (TH1D*)hData->Clone(Form("%s_Fit_Prim_DCAxy_%02d",fParticle.Data(),ibin));
+       TH1D* hDebugSec =  (TH1D*)hData->Clone(Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),secname,ibin));
+       
+       hDebugFit->Reset();
+       hDebugPrim->Reset();
+       hDebugSec->Reset();
+       
+       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 =  (TH1D*)hData->Clone(Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
+       TH1D* hDebugPrim = (TH1D*)hData->Clone(Form("%s_Fit_Prim_DCAxy_%02d",fParticle.Data(),ibin));
+       TH1D* hDebugMat =  (TH1D*)hData->Clone(Form("%s_Fit_Mat_DCAxy_%02d",fParticle.Data(),ibin));
+       TH1D* hDebugFdwn = (TH1D*)hData->Clone(Form("%s_Fit_Fdwn_DCAxy_%02d",fParticle.Data(),ibin));
+       
+       hDebugFit->Reset();
+       hDebugPrim->Reset();
+       hDebugMat->Reset();
+       hDebugFdwn->Reset();
+       
+       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
+{
+//
+// generate a flat dcaxy-pt distribution
+//
+       TF1* fnc = new TF1("gen","1", fMinDCAxy, fMaxDCAxy);
+       TH2D* h = (TH2D*)hDCAxyPt->Clone(name.Data());
+       h->Reset();
+       
+       TAxis* xAxis = h->GetXaxis();
+       TAxis* yAxis = h->GetYaxis();
+       
+       for(Int_t i=1; i<xAxis->GetNbins()+1; ++i)
+       {
+               Double_t pt = xAxis->GetBinCenter(i);
+               
+               for(Int_t j=1; j<yAxis->GetNbins()+1; ++j)
+               {
+                       for(Int_t k=0; k<max; ++k) h->Fill(pt,fnc->GetRandom());
+               }
+       }
+       
+       delete fnc;
+       return h;
+}
+
+TF1* AliLnSecondaries::GetMatFraction(const TString& name) const
+{
+//
+// model the material fraction as an exponential + a constant
+//
+       TF1* fnc = new TF1(name.Data(), "[0]*exp(-[1]*x)+[2]",0,10);
+       fnc->SetParameters(1.7,5.5,0.01);
+       
+       fnc->SetParLimits(2,0,0.02); // asymptotic value
+       
+       return fnc;
+}
+
+TF1* AliLnSecondaries::GetFdwnFraction(const TString& name) const
+{
+//
+// model the feed-down fraction as an exponential + constant
+//
+       TF1* fnc = new TF1(name, "[0]*exp(-[1]*x)+[2]",0,10);
+       fnc->SetParameters(0.5,1,0.02);
+       
+       fnc->SetParLimits(2,0,0.3); // asymptotic value
+       
+       return fnc;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.h
new file mode 100644 (file)
index 0000000..941be4b
--- /dev/null
@@ -0,0 +1,99 @@
+#ifndef ALILNSECONDARIES_H
+#define ALILNSECONDARIES_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// removal of secondaries using DCA distributions
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class TString;
+class TH1D;
+class TH2D;
+class RooWorkspace;
+class RooDataSet;
+class RooRealVar;
+class TF1;
+class TFractionFitter;
+
+class AliLnSecondaries: public TObject
+{
+  public:
+       
+       AliLnSecondaries(const TString& particle, const TString& dataFilename, const TString& simuFilename, const TString& outputFilename, const TString& otag);
+       virtual ~AliLnSecondaries();
+       
+       const TString* GetOutputFilename() const { return &fOutputFilename; }
+       
+       Int_t Exec();
+       
+       void SetParticle(const TString& particle) { fParticle = particle; }
+       
+       void SetOutputTag(const TString& tag) { fOutputTag = tag; }
+       
+       void SetCorBins(Int_t lowbin, Int_t hibin) { fLowBin = lowbin; fHiBin = hibin; }
+       void SetDCAxyInterval(Double_t lowdca, Double_t hidca) { fMinDCAxy = lowdca; fMaxDCAxy = hidca; }
+       
+       void SetNBin(Int_t nbin) { fNbin = nbin; }
+       
+       void SetProcedure(Int_t prod) { fFracProc=prod; }
+       void SetMatDCAxyModel(Int_t model) { fMatDCAxyMod = model; }
+       void SetScalingFactors(Double_t mat, Double_t fd) { fScMat=mat; fScFd=fd; }
+       
+       enum { kTFractionFitter=0, kRooFit, kMonteCarlo };
+       enum { kGeantDCAxy=0, kFlatDCAxy, kGaussianDCAxy};
+       
+  private:
+       AliLnSecondaries(const AliLnSecondaries& other);
+       AliLnSecondaries& operator=(const AliLnSecondaries& other);
+       
+       void GetFraction(TH1D* hPrimPt) const;
+       void GetFraction(TH1D* hPrimFracPt, TH1D* hSecFracPt, const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hSecDCAxyPt, const TString& secName) const;
+       void GetFraction(TH1D* hFracPt[3], const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hMatDCAxyPt, const TH2D* hFdwnDCAxyPt) const;
+       
+       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 char* secName) const;
+       
+       void GetRooFitFractions(Double_t* frac, Double_t* err, const TH1D* hData, const TH1D* hPrim, const TH1D* hSec, Int_t ibin, const char* sec) 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;
+       //TH2D* GetGaussianDCAxyPt(Double_t norm, Double_t mean, Double_t sigma, const TH2D* hDCAxyPt, const char* name) const;
+       
+       TF1* GetMatFraction(const TString& name) const;
+       TF1* GetFdwnFraction(const TString& name) const;
+       
+       TH1D* ZeroClone(TH1D* h, const char* name) const;
+       void WriteTFFdebug(TH1D* hData, TFractionFitter* fit, Int_t status, Int_t ibin, const char* contrib[], Double_t* frac, Int_t kmax) const;
+       
+  private:
+  
+       TString fParticle; // particle
+       
+       TString fDataFilename; // data filename
+       TString fSimuFilename; // simulation filename
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for the ouput
+       
+       Int_t fLowBin; // low bin for the corrections
+       Int_t fHiBin ; // high bin for the corrections
+       
+       Int_t fNbin; // for rebinning DCA distributions
+       
+       Double_t fMinDCAxy; // low DCAxy value
+       Double_t fMaxDCAxy; // high DCAxy value
+       
+       Int_t fFracProc;  // procedure for estimating the fractions
+       Int_t fMatDCAxyMod; // DCAxy model for secondaries from materials
+       
+       Double_t fScMat; // scaling factor for material fraction
+       Double_t fScFd; // scaling factor for feed-down fraction
+       
+       ClassDef(AliLnSecondaries,1)
+};
+
+#endif // ALILNSECONDARIES_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx
new file mode 100644 (file)
index 0000000..f9ae536
--- /dev/null
@@ -0,0 +1,313 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// invariant differential yields and cross sections
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+
+#include "AliLnSpectra.h"
+#include "B2.h"
+
+ClassImp(AliLnSpectra)
+
+AliLnSpectra::AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag, const Double_t xsec[3])
+: TObject()
+, fParticle(particle)
+, fPtFilename(ptFilename)
+, fTag(tag)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+, fYMin(-0.5)
+, fYMax(0.5)
+, fNormToInel(1)
+, fIsOnlyGen(0)
+, fSysErr(1)
+{
+//
+// constructor
+//
+       for(Int_t i=0; i<3; ++i) fInelXsec[i] = xsec[i];
+}
+
+AliLnSpectra::~AliLnSpectra()
+{
+//
+// destructor
+//
+}
+
+Int_t AliLnSpectra::Exec()
+{
+//
+// (invariant) differential yield projection in pt for |y| < 0.5
+//
+       using namespace std;
+       
+       TFile* finput = new TFile(fPtFilename.Data(), "read");
+       if(finput->IsZombie()) exit(1);
+       
+       // pt and number of events
+       
+       TH1D* hPt = (TH1D*)FindObj(finput, fParticle + "_Pt");
+       TH1D* hStats = (TH1D*)FindObj(finput, fParticle + "_Stats");
+       
+       Double_t nEvent = hStats->Integral(3,3);
+       if(fNormToInel) nEvent = hStats->Integral(4,4);
+       
+       // ouputfile
+       
+       TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
+       
+       foutput->mkdir(fOutputTag.Data());
+       foutput->cd(fOutputTag.Data());
+       
+       // differential yield
+       
+       TGraphErrors* grDYieldPt = this->GetDiffYieldPt(hPt, nEvent, fParticle + "_DiffYield_Pt");
+       grDYieldPt->Write();
+       
+       TGraphErrors* grSysErrDYieldPt = this->AddSystError(grDYieldPt, fSysErr, fParticle + "_SysErr_DiffYield_Pt");
+       grSysErrDYieldPt->Write();
+       
+       TF1* fncTsallis0 = this->TsallisDiffYield(GetMass(fParticle), fParticle + "_Fit_DiffYield_Pt");
+       fncTsallis0->SetParLimits(0, 0, 1);
+       fncTsallis0->SetParLimits(1, 4, 50);
+       fncTsallis0->SetParLimits(2, 0.01, 10);
+       
+       grSysErrDYieldPt->Fit(fncTsallis0,"RNQ");
+       fncTsallis0->Write();
+       
+       // invariant differential yield
+       
+       TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, fParticle + "_InvDiffYield_Pt");
+       grInvDYieldPt->Write();
+       
+       TGraphErrors* grSysErrInvDYieldPt = this->AddSystError(grInvDYieldPt, fSysErr, fParticle + "_SysErr_InvDiffYield_Pt");
+       grSysErrInvDYieldPt->Write();
+       
+       TF1* fncTsallis1 = this->Tsallis(GetMass(fParticle), fParticle + "_Fit_InvDiffYield_Pt");
+       fncTsallis1->SetParLimits(0, 0, 1);
+       fncTsallis1->SetParLimits(1, 4, 50);
+       fncTsallis1->SetParLimits(2, 0.01, 10);
+       
+       grSysErrInvDYieldPt->Fit(fncTsallis1,"RNQ");
+       fncTsallis1->Write();
+       
+       // invariant differential cross section
+       
+       TGraphErrors* grInvXsectPt = GetInvDiffXsectionPt(grInvDYieldPt, fInelXsec, fParticle + "_InvDiffXSection_Pt");
+       grInvXsectPt->Write();
+       
+       TGraphErrors* grSysErrInvXsectPt = GetInvDiffXsectionPt(grSysErrInvDYieldPt, fInelXsec, fParticle + "_SysErr_InvDiffXSection_Pt");
+       grSysErrInvXsectPt->Write();
+       
+       TF1* fncTsallis2 = this->Tsallis(GetMass(fParticle), fInelXsec[0], fParticle + "_Fit_InvDiffXSection_Pt");
+       fncTsallis2->SetParLimits(0, 0, 1);
+       fncTsallis2->SetParLimits(1, 4, 50);
+       fncTsallis2->SetParLimits(2, 0.01, 10);
+       
+       grSysErrInvXsectPt->Fit(fncTsallis2,"RNQ");
+       fncTsallis2->Write();
+       
+       // clean
+       
+       delete fncTsallis0;
+       delete fncTsallis1;
+       delete fncTsallis2;
+       
+       delete grDYieldPt;
+       delete grInvDYieldPt;
+       delete grSysErrDYieldPt;
+       delete grSysErrInvDYieldPt;
+       delete grInvXsectPt;
+       
+       delete foutput;
+       delete finput;
+       
+       return 0;
+}
+
+TGraphErrors* AliLnSpectra::GetDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const
+{
+//
+// projection of the differential yield in pt
+// only statistical error of N is taken into account
+//
+       TGraphErrors* grDYieldPt = new TGraphErrors();
+       grDYieldPt->SetName(name.Data());
+       
+       Double_t dy = fYMax - fYMin;
+       
+       TAxis* xaxis = hPt->GetXaxis();
+       
+       for(Int_t i=1, j=0; i < xaxis->GetNbins()+1; ++i)
+       {
+               Double_t pt  = xaxis->GetBinCenter(i);
+               Double_t dpt = xaxis->GetBinWidth(i);
+               
+               Double_t n = hPt->GetBinContent(i);
+               Double_t eN = hPt->GetBinError(i);
+               
+               if(n==0) continue;
+               
+               Double_t yield = n/(nEvent*dpt*dy);
+               Double_t yError = yield*eN/n;
+               
+               grDYieldPt->SetPoint(j,pt,yield);
+               grDYieldPt->SetPointError(j++,0,yError);
+       }
+       
+       return grDYieldPt;
+}
+
+TGraphErrors* AliLnSpectra::GetInvDiffYieldPt(const TGraphErrors* grDYieldPt, const TString& name) const
+{
+//
+// projection of the invariant differential yield in pt
+// only statistical error of N is taken into account
+//
+       TGraphErrors* grInvDYieldPt = new TGraphErrors();
+       grInvDYieldPt->SetName(name.Data());
+       
+       for(Int_t i=0; i < grDYieldPt->GetN(); ++i)
+       {
+               Double_t pt, y;
+               grDYieldPt->GetPoint(i,pt,y);
+               
+               Double_t errpt = grDYieldPt->GetErrorX(i);
+               Double_t erry = grDYieldPt->GetErrorY(i);
+               
+               y = y/(2.*TMath::Pi()*pt);
+               erry = erry / (2.*TMath::Pi()*pt);
+               
+               grInvDYieldPt->SetPoint(i, pt, y);
+               grInvDYieldPt->SetPointError(i, errpt, erry);
+       }
+       
+       return grInvDYieldPt;
+}
+
+TGraphErrors* AliLnSpectra::GetInvDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const
+{
+//
+// projection of the invariant differential yield in pt
+// only statistical error of N is taken into account
+//
+       TGraphErrors* grDYieldPt = this->GetDiffYieldPt(hPt, nEvent, "_tmp_diff_yield_");
+       TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, name);
+       
+       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;
+}
+
+TF1* AliLnSpectra::Tsallis(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
+{
+//
+// Tsallis distribution
+// Phys. Rev. C 83, 064903 (2011)
+// Phys. Rev. C 75, 064901 (2007)
+//
+       TF1* fnc = new TF1(name.Data(), Form("[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
+       fnc->SetParNames("dN/dy","n","C");
+       fnc->SetParameters(0.1, 7, 0.2);
+       
+       return fnc;
+}
+
+TF1* AliLnSpectra::Tsallis(Double_t m0, Double_t xsect, const TString& name, Double_t xmin, Double_t xmax) const
+{
+//
+// Tsallis distribution to fit to invariant cross section
+// Phys. Rev. C 83, 064903 (2011)
+// Phys. Rev. C 75, 064901 (2007)
+//
+       TF1* fnc = new TF1(name.Data(), Form("%f*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))",xsect, m0, m0, m0, m0), xmin, xmax);
+       fnc->SetParNames("dN/dy","n","C");
+       fnc->SetParameters(0.1, 7, 0.2);
+       
+       return fnc;
+}
+
+TF1* AliLnSpectra::TsallisDiffYield(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
+{
+//
+// Tsallis distribution to fit differential yield
+// Phys. Rev. C 83, 064903 (2011)
+// Phys. Rev. C 75, 064901 (2007)
+//
+       TF1* fnc = new TF1(name.Data(), Form("x*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/([1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
+       fnc->SetParNames("dN/dy","n","C");
+       fnc->SetParameters(0.1, 7, 0.2);
+       
+       return fnc;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h b/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.h
new file mode 100644 (file)
index 0000000..1c84c2e
--- /dev/null
@@ -0,0 +1,72 @@
+#ifndef ALILNSPECTRA_H
+#define ALILNSPECTRA_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// invariant differential yields and cross sections
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+
+class TGraphErrors;
+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]);
+       
+       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 SetNormalizeToINEL(Bool_t flag=1) { fNormToInel = 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; }
+       
+       TF1* Tsallis(Double_t m0, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
+       TF1* Tsallis(Double_t m0, Double_t xsect, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
+       TF1* TsallisDiffYield(Double_t m0, const TString& name, Double_t xmin=0., Double_t xmax=10.) const;
+       
+  private:
+       AliLnSpectra(const AliLnSpectra& other);
+       AliLnSpectra& operator=(const AliLnSpectra& other);
+       
+  private:
+       TString fParticle; // particle name
+       
+       TString fPtFilename; // pt filename
+       TString fTag; // tag for pt file
+       
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for output file
+       
+       Double_t fYMin; // rapidity low limit
+       Double_t fYMax; // rapidity high limit
+       
+       Bool_t fNormToInel; // normalize 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)
+};
+
+#endif // ALILNSPECTRA_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx b/PWGLF/SPECTRA/Nuclei/B2/AliLnUnfolding.cxx
new file mode 100644 (file)
index 0000000..6dceed7
--- /dev/null
@@ -0,0 +1,86 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// unfolding correction for the pt
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+//#include <Riostream.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TString.h>
+#include <TF1.h>
+
+#include "AliLnUnfolding.h"
+#include "B2.h"
+
+ClassImp(AliLnUnfolding)
+
+AliLnUnfolding::AliLnUnfolding(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag)
+: TObject()
+, fParticle(particle)
+, fSimuFilename(simuFilename)
+, fOutputFilename(outputFilename)
+, fOutputTag(otag)
+{
+//
+// constructor
+//
+}
+
+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
new file mode 100644 (file)
index 0000000..96bc032
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef ALILNUNFOLDING_H
+#define ALILNUNFOLDING_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// unfolding correction for the pt
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TString.h>
+
+class TString;
+
+class AliLnUnfolding: public TObject
+{
+  public:
+       
+       AliLnUnfolding(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag);
+       virtual ~AliLnUnfolding();
+       
+       const TString* GetOutputFilename() const { return &fOutputFilename; }
+       
+       Int_t Exec();
+       
+       void SetParticle(const TString& particle) { fParticle = particle; }
+       
+       void SetOutputTag(const TString& tag) { fOutputTag = tag; }
+       
+  private:
+       AliLnUnfolding(const AliLnUnfolding& other);
+       AliLnUnfolding& operator=(const AliLnUnfolding& other);
+       
+  private:
+
+       TString fParticle; // particle
+       
+       TString fSimuFilename; // simulation filename
+       TString fOutputFilename; // output filename
+       TString fOutputTag; // tag for the ouput
+       
+       ClassDef(AliLnUnfolding,1)
+};
+
+#endif // ALILNUNFOLDING_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/B2.h b/PWGLF/SPECTRA/Nuclei/B2/B2.h
new file mode 100644 (file)
index 0000000..21e424a
--- /dev/null
@@ -0,0 +1,99 @@
+#ifndef B2_H
+#define B2_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// some common functions
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TObject.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TString.h>
+#include <TH1D.h>
+#include <cstdlib>
+
+TObject* FindObj(TFile* f, const TString& name)
+{
+//
+// check if the object exists
+//
+       TObject* obj = f->Get(name.Data());
+       if(obj == 0)
+       {
+               f->Error("Get","%s not found",name.Data());
+               exit(1);
+       }
+       
+       return obj;
+}
+
+TObject* 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;
+}
+
+TObject* FindObj(const TList* l, const TString& name)
+{
+//
+// check if the object exists
+//
+       TObject* obj = l->FindObject(name.Data());
+       if(obj == 0)
+       {
+               l->Error("FindObject","%s not found",name.Data());
+               exit(1);
+       }
+       
+       return obj;
+}
+
+TH1D* Divide(const TH1* hX, const TH1* hY, const TString& name)
+{
+//
+// clone and divide
+//
+       TH1D* q = (TH1D*)hX->Clone(name.Data());
+       q->Sumw2();
+       q->Divide(hY);
+       return q;
+}
+
+Double_t GetMass(const TString& name)
+{
+//
+// return particle mass
+//
+       TString particle=name;
+       particle.ToLower();
+       
+       if(particle.Contains("electron")) return 0.000511;
+       if(particle.Contains("muon")) return 0.10566;
+       if(particle.Contains("pion")) return 0.13957;
+       if(particle.Contains("kaon")) return 0.49368;
+       if(particle.Contains("proton")) return 0.93827;
+       if(particle.Contains("deuteron")) return 1.87561;
+       if(particle.Contains("triton")) return 2.80925;
+       if(particle.Contains("he3")) return 2.80923;
+       if(particle.Contains("alpha")) return 3.727417;
+       
+       //cerr << "Warning unknown particle " << particle << endl;
+       
+       return 0;
+}
+
+#endif // B2_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/Config.h b/PWGLF/SPECTRA/Nuclei/B2/Config.h
new file mode 100644 (file)
index 0000000..fcf75a6
--- /dev/null
@@ -0,0 +1,339 @@
+#ifndef CONFIG_H
+#define CONFIG_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// functions for configs
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TString.h>
+#include <TStyle.h>
+#include <cstdlib>
+
+namespace B2mult
+{
+//
+// multiplicity classes
+//
+       const Int_t kNmult = 6;
+       const TString kMultClass[kNmult]   = { "nch0002", "nch0204", "nch0408", "nch0811", "nch1120", "nch20xx" };
+       const Double_t kKNOmult[kNmult]    = { 0.20, 0.60, 1.01, 1.60, 2.60, 4.35 };
+       const Double_t kKNOmultErr[kNmult] = { 0, 0, 0, 0, 0, 0 };
+};
+
+Double_t GetVertexCorrection(const TString& period)
+{
+//
+// vertex correction factor for |Vz| < 10 cm and LHC12a5x simulations
+//
+       TString name = period;
+       name.ToLower();
+       
+       if(name == "lhc10c900")   return 1.0142;
+       if(name == "lhc10b")      return 1.0006;
+       if(name == "lhc10c")      return 1.0007;
+       if(name == "lhc10d")      return 1.0224;
+       if(name == "lhc10e")      return 1.0423;
+       if(name == "lhc11a_wsdd") return 1.0728;
+       
+       return 1;
+}
+
+TString GetCollSystem(const TString& period)
+{
+//
+// translate period name into colliding system name
+//
+       TString name = period;
+       name.ToLower();
+       
+       if(name == "lhc10c900")    return "pp0.9TeV";
+       if(name == "lhc10b_pass2") return "pp7TeV";
+       if(name == "lhc10c_pass2") return "pp7TeV";
+       if(name == "lhc10b")       return "pp7TeV";
+       if(name == "lhc10c")       return "pp7TeV";
+       if(name == "lhc10d")       return "pp7TeV";
+       if(name == "lhc10e")       return "pp7TeV";
+       if(name.Contains("lhc10")) return "pp7TeV";
+       if(name == "lhc11a_wsdd")  return "pp2.76TeV";
+       if(name == "lhc11a_wosdd") return "pp2.76TeV";
+       if(name.Contains("900"))   return "pp0.9TeV";
+       if(name.Contains("7000"))  return "pp7TeV";
+       if(name.Contains("2760"))  return "pp2.76TeV";
+       
+       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
+// http://arxiv.org/abs/1208.4968
+//
+       TString collsystem = GetCollSystem(period);
+       TString trigger = trigname;
+       trigger.ToLower();
+       
+       if( trigger == "mbor"  && collsystem == "pp0.9TeV" )
+       {
+               trigeff[0] = 0.910;
+               trigeff[1] = 0; // stat. error
+               trigeff[2] = 0.032; // syst. error
+       }
+       else if( trigger == "mbor"  && collsystem == "pp2.76TeV")
+       {
+               trigeff[0] = 0.881;
+               trigeff[1] = 0;
+               trigeff[2] = 0.059;
+       }
+       else if( trigger == "mbor"  && collsystem == "pp7TeV")
+       {
+               trigeff[0] = 0.852;
+               trigeff[1] = 0;
+               trigeff[2] = 0.062;
+       }
+       else if( trigger == "mband" && collsystem == "pp0.9TeV" )
+       {
+               trigeff[0] = 0.763;
+               trigeff[1] = 0;
+               trigeff[2] = 0.022;
+       }
+       else if( trigger == "mband" && collsystem == "pp2.76TeV")
+       {
+               trigeff[0] = 0.760;
+               trigeff[1] = 0;
+               trigeff[2] = 0.052;
+       }
+       else if( trigger == "mband" && collsystem == "pp7TeV")
+       {
+               trigeff[0] = 0.742;
+               trigeff[1] = 0;
+               trigeff[2] = 0.050;
+       }
+       else
+       {
+               std::cerr << "Warning: unknown trigger/colliding system " << trigname << "/" << collsystem << " for period " << period << std::endl;
+       }
+}
+
+TString GetSimuPeriod(const TString& period)
+{
+//
+// simulation code for the given period
+//
+       TString name = period;
+       name.ToLower();
+       
+       if(period=="lhc10c900")   return "lhc10e13";
+       if(period=="lhc10b")      return "lhc10d1";
+       if(period=="lhc10c")      return "lhc10d4";
+       //if(period=="lhc10b")      return "lhc12a5bb";
+       //if(period=="lhc10c")      return "lhc12a5bc";
+       if(period=="lhc10d")      return "lhc10f6a";
+       if(period=="lhc10e")      return "lhc10e21";
+       if(period=="lhc11a_wsdd") return "lhc11e3a_wsdd";
+       
+       return "";
+}
+
+TString GetSimuFixPeriod(const TString& period)
+{
+//
+// simulation fix code for the given period
+//
+       TString name = period;
+       name.ToLower();
+       
+       if(period=="lhc10c900")   return "lhc12a5a";
+       if(period=="lhc10b")      return "lhc12a5bb";
+       if(period=="lhc10c")      return "lhc12a5bc";
+       if(period=="lhc10d")      return "lhc12a5bd";
+       if(period=="lhc10e")      return "lhc12a5be";
+       if(period=="lhc10bcde")   return "lhc12a5b";
+       if(period=="lhc11a_wsdd") return "lhc12a5c_wsdd";
+       
+       return "";
+}
+
+TString MakeInputName(const TString& species, const TString& period, const TString& trksel)
+{
+//
+// make input name for data
+//
+       return species + "-" + period + "-" + trksel;
+}
+
+TString MakeSimuName(const TString& species, const TString& period, const TString& trksel)
+{
+//
+// make input name for simulation
+//
+       TString simu = (species == "Proton") ? GetSimuPeriod(period) : GetSimuFixPeriod(period);
+       
+       if(simu=="")
+       {
+               std::cerr << "No simulation for period: " << period << std::endl;
+               exit(1);
+       }
+       
+       return species + "-" + simu + "-" + trksel;
+}
+
+TString MakeSimuFixName(const TString& species, const TString& period, const TString& trksel, Bool_t g3Fluka=0)
+{
+//
+// make input name for simulation fix
+//
+       TString simufix = (species == "Proton" && g3Fluka) ? GetSimuPeriod(period) : GetSimuFixPeriod(period);
+       
+       if(simufix=="")
+       {
+               std::cerr << "No simulation fix for period: " << period << std::endl;
+               exit(1);
+       }
+       
+       return species + "-" + simufix + "-" + trksel;
+}
+
+TString MakeOutputName(const TString& species, const TString& outputTag)
+{
+//
+// make output name
+//
+       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="")
+{
+//
+// draw corrections for the given species
+//
+       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)
+{
+//
+// draw correction for secondaries
+//
+       const TString kParticle[] = { species, Form("Anti%s", species.Data())};
+       
+       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));
+       }
+}
+
+void DrawPtDebug(const TString& pt, const TString& tag, const TString& species, Bool_t m2pid=0, Int_t lowm2bin=9, Int_t him2bin=17)
+{
+//
+// draw pt debug for the particle species
+//
+       const TString kParticle[] = { species, Form("Anti%s", species.Data())};
+       
+       for(Int_t i=0; i<2; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawPt.C+g(\"%s\",\"%s\",\"%s\",%d,%d, %d)", pt.Data(), tag.Data(), kParticle[i].Data(), m2pid, lowm2bin, him2bin));
+       }
+}
+
+void DrawOutputRatio(const TString& ratio, const TString& tag, const TString& species)
+{
+//
+// draw ratio antiparticle/particle
+//
+       gROOT->ProcessLine(Form(".x DrawRatio.C+g(\"%s\",\"%s\",\"%s\")", ratio.Data(), tag.Data(), species.Data()));
+}
+
+void DrawOutputSpectra(const TString& spectra, const TString& tag, const TString& species)
+{
+//
+// draw spectra for the particle species
+//
+       const TString kParticle[] = { species, Form("Anti%s", species.Data())};
+       
+       for(Int_t i=0; i<2; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawSpectra.C+g(\"%s\",\"%s\",\"%s\")", spectra.Data(), tag.Data(), kParticle[i].Data()));
+       }
+}
+
+void DrawOutputSpectraMult(const TString& spectra, const TString& species, Double_t ymin, Double_t ymax, Int_t option=0, const TString& refdir="")
+{
+//
+// draw spectra for each multiplicity class
+//
+       const TString kParticle[] = { species, Form("Anti%s", species.Data())};
+       
+       for(Int_t i=0; i<2; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_InvDiffYield_Pt\",\"%s\",0,4.5, %g, %g, \"p_{T} (GeV/c)\", \"#frac{1}{2#piN_{inel}} #frac{d^{2}N}{p_{T}dp_{T}dy} (GeV^{-2}c^{3})\", %d, \"c%d\",\"%s\")", spectra.Data(), kParticle[i].Data(), refdir.Data(),ymin, ymax, option, i, kParticle[i].Data()));
+       }
+}
+
+#endif // CONFIG_H
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C b/PWGLF/SPECTRA/Nuclei/B2/macros/B2.C
new file mode 100644 (file)
index 0000000..e248955
--- /dev/null
@@ -0,0 +1,66 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// coalescence parameter
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TSystem.h>
+#include <TROOT.h>
+#include <TFileMerger.h>
+#include <TString.h>
+
+#include "AliLnB2.h"
+
+Int_t B2(const TString& pSpectra   = "~/alice/output/Proton-lhc10d-Spectra.root",
+         const TString& ptag       = "lhc10d",
+         const TString& dSpectra   = "~/alice/output/Deuteron-lhc10d-Spectra.root",
+         const TString& dtag       = "lhc10d",
+         const TString& outputfile = "~/alice/output/lhc10d-B2.root",
+         const TString& otag       = "lhc10d")
+{
+//
+// coalescence parameter
+//
+       const Int_t kNpart      = 2;
+       const TString kPrefix[] = { "", "Anti" };
+       const Int_t kCharge[]   = { 1, -1 };
+       
+       TFileMerger m;
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               TString b2file = kPrefix[i] + "B2.root";
+               
+               AliLnB2 b2(pSpectra, ptag, dSpectra, dtag, b2file, otag, 2, kCharge[i]);
+               
+               b2.Run();
+               
+               m.AddFile(b2file.Data(),0);
+       }
+       
+       m.OutputFile(outputfile.Data());
+       m.Merge();
+       
+       gSystem->Exec("rm -f B2.root AntiB2.root");
+       
+       // draw
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawB2.C+g(\"%s\",\"%s\",\"%s\")", outputfile.Data(), otag.Data(), kPrefix[i].Data()));
+       }
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C
new file mode 100644 (file)
index 0000000..87dc522
--- /dev/null
@@ -0,0 +1,213 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// B2 as a function of multiplicity
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TSystem.h>
+#include <TROOT.h>
+#include <TFileMerger.h>
+#include <TString.h>
+#include <TFile.h>
+#include <TGraphErrors.h>
+
+#include "AliLnB2.h"
+#include "B2.h"
+#include "Config.h"
+
+Double_t GetCd(Double_t z)
+{
+//
+// parameterization of <Cd> as a function of multiplicity
+// from ALICE Rlong and Rside measurements
+//
+       return 0.046133 + 0.0484458*z;
+}
+
+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-Mult-Pt.root",
+             const TString& outputPtMult = "~/alice/output/B2-Pt-Mult.root",
+             const TString& otag         = "lhc10d")
+{
+//
+// B2 as a function of multiplicity
+//
+       using namespace B2mult;
+       
+       const Int_t kNpart = 2;
+       
+       const Int_t kNMinPt = 0;
+       const Int_t kNMaxPt = 6;
+       
+       const TString kPrefix[] = { "", "Anti"};
+       const TString kSuffix[] = { "", "bar" };
+       Int_t kCharge[]         = {1, -1};
+       
+       // B2 as a function of pt for each multiplicity class
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               TFileMerger m;
+               
+               for(Int_t j=0; j<kNpart; ++j)
+               {
+                       TString b2file = kPrefix[j] + "B2.root";
+                       
+                       AliLnB2 b2(pSpectra, ptag + "-" + kMultClass[i], dSpectra, dtag + "-" + kMultClass[i], b2file, otag + "-" + kMultClass[i], 2, kCharge[j]);
+                       
+                       b2.SetCd(GetCd(kKNOmult[i]));
+                       
+                       b2.Run();
+                       
+                       m.AddFile(b2file.Data(),0);
+               }
+               
+               // merge B2 and B2bar
+               
+               TString outputfile = otag + "-" + kMultClass[i] + "-B2.root";
+               
+               m.OutputFile(outputfile.Data());
+               m.Merge();
+               
+               gSystem->Exec("rm -f B2.root AntiB2.root");
+       }
+       
+       // merge multiplicity classes
+       
+       TFileMerger m;
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               TString b2 = otag + "-" + kMultClass[i] + "-B2.root";
+               m.AddFile(b2.Data(),0);
+       }
+       
+       m.OutputFile(outputMultPt.Data());
+       m.Merge();
+       
+       // delete tmp files
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               gSystem->Exec(Form("rm -f %s-%s-B2.root",otag.Data(),kMultClass[i].Data()));
+       }
+       
+       // B2 as a function of multiplicity for each pt
+       
+       TFile* finput = new TFile(outputMultPt.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       TGraphErrors* grB2pt[kNpart][kNmult];
+       TGraphErrors* grR3pt[kNpart][kNmult];
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               for(Int_t j=0; j<kNmult; ++j)
+               {
+                       grB2pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("B2%s_Pt", kSuffix[i].Data()));
+                       grR3pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("R3%s_Pt", kSuffix[i].Data()));
+               }
+       }
+       
+       TFile* foutput = new TFile(outputPtMult.Data(),"recreate");
+       
+       Double_t* pt = grB2pt[0][0]->GetX();
+       TString ptLabel[kNMaxPt];
+       
+       if(kNMaxPt > grB2pt[0][0]->GetN())
+       {
+               std::cerr << "max pt too big" << std::endl;
+               exit(1);
+       }
+       
+       for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
+       {
+               ptLabel[i] = Form("pt%.02fA",pt[i]);
+               foutput->mkdir(ptLabel[i].Data());
+       }
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               for(Int_t j=kNMinPt; j<kNMaxPt; ++j)
+               {
+                       Double_t B2[kNmult];
+                       Double_t B2StatErr[kNmult];
+                       
+                       Double_t R3[kNmult];
+                       Double_t R3StatErr[kNmult];
+                       
+                       for(Int_t k=0; k<kNmult; ++k)
+                       {
+                               Double_t x, y, ey;
+                               grB2pt[i][k]->GetPoint(j,x,y);
+                               ey = grB2pt[i][k]->GetErrorY(j);
+                               
+                               B2[k] = y;
+                               B2StatErr[k] = ey;
+                               
+                               grR3pt[i][k]->GetPoint(j,x,y);
+                               ey = grR3pt[i][k]->GetErrorY(j);
+                               
+                               R3[k] = y;
+                               R3StatErr[k] = ey;
+                       }
+                       
+                       TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
+                       grB2Mult->SetName(Form("B2%s_Zmult", kSuffix[i].Data()));
+                       
+                       TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
+                       grR3Mult->SetName(Form("R3%s_Zmult", kSuffix[i].Data()));
+                       
+                       foutput->cd(ptLabel[j].Data());
+                       
+                       grB2Mult->Write();
+                       grR3Mult->Write();
+                       
+                       delete grB2Mult;
+                       delete grR3Mult;
+               }
+       }
+       
+       delete foutput;
+       delete finput;
+       
+       // particle ratios
+       
+       gROOT->ProcessLine(Form(".x RatioMult.C+g(\"%s\",\"%s\",\"%s\",\"%s\")",pSpectra.Data(), ptag.Data(),dSpectra.Data(), dtag.Data()));
+       
+       // draw B2 as a function of pt
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Pt\",\"\",0,2, 1.e-3, 7.e-2,\"p_{T}/A (GeV/c)\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2pt\",\"B2%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
+               
+               gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Pt\",\"\",0,2, 0, 1.7,\"p_{T}/A (GeV/c)\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3pt\",\"R3%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
+       }
+       
+       // draw B2 as a function of z
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Zmult\",\"\",0,7, 1.e-3, 7.e-2,\"z\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2z\",\"B2%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
+               
+               gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Zmult\",\"\",0,8, 0, 4,\"z\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3z\",\"R3%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
+       }
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Deuteron_TOF_LHC10x.C
new file mode 100644 (file)
index 0000000..0b2b4e8
--- /dev/null
@@ -0,0 +1,124 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+// LHC10x config for deuterons and antideuterons\r
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>\r
+\r
+#include <TROOT.h>\r
+#include <TString.h>\r
+#include "AliLnDriver.h"\r
+#include "Config.h"\r
+\r
+Int_t Config_Deuteron_TOF_LHC10x(const TString& inputDir   = "~/alice/input",\r
+                                 const TString& outputDir  = "~/alice/output",\r
+                                 const TString& period     = "lhc10d",\r
+                                 const TString& outputTag  = "lhc10d",\r
+                                 const TString& multTag    = "",\r
+                                 const TString& multCorTag = "",\r
+                                 Bool_t normToInel         = 1,  // for mult\r
+                                 Bool_t drawOutput         = 1)  // for batch\r
+{\r
+//\r
+// lhc10b, lhc10c, lhc10d, lhc10e config for deuterons and antideuterons\r
+//\r
+       const TString  kSpecies     = "Deuteron";\r
+       const TString  kTrkSel      = "its_tpc_tof_dca-tpc3";\r
+       const TString  kTrigName    = "mbor";\r
+       const Bool_t   kVtxCorr     = 0;\r
+       const Double_t kVtxCorrVal  = GetVertexCorrection(period);\r
+       const Int_t    kPtBin[2]    = {4,13};\r
+       const Int_t    kM2Bin[2]    = {8,18};\r
+       const Bool_t   kPidM2       = 1;\r
+       const Bool_t   kUnfolding   = 1;\r
+       const Int_t    kIter        = 5;\r
+       const Bool_t   kFakeTracks  = 0;\r
+       const Bool_t   kSecondaries = 1;\r
+       const Int_t    kSecProd     = 0; // 0 tff, 1 roofit, 2 mc\r
+       const Int_t    kMatDCAxyMod = 1; // 0 geant, 1 flat\r
+       const Int_t    kNbin        = 5;\r
+       const Double_t kDCAxy[2]    = {-0.2,0.2};\r
+       const Double_t kM2Bkg[2]    = {2.2,5.};\r
+       const Double_t kM2tpc[2]    = {2.,6.5};\r
+       const Bool_t   kEfficiency  = 1;\r
+       const Bool_t   kFitFrac     = 1;\r
+       const Double_t kSysErr[2]   = {0.10, 0.11} ;\r
+       \r
+       Double_t xsec[3];\r
+       GetInelXSection(xsec, period);\r
+       \r
+       Double_t trigEff[3];\r
+       GetTriggerEfficiency(trigEff, kTrigName, period);\r
+       \r
+       // input and output filenames\r
+       \r
+       TString inputData     = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + ".root";\r
+       TString inputSimu     = inputDir + "/" + period + "/" + MakeSimuName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
+       TString inputSimuFix  = inputDir + "/" + period + "/" + MakeSimuFixName(kSpecies, period, kTrkSel+multCorTag) + ".root";\r
+       TString inputCorr     = inputDir + "/" + period + "/" + MakeInputName(kSpecies, period, kTrkSel+multTag) + "-corr.root";\r
+       \r
+       TString outputPt      = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";\r
+       TString outputRatio   = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";\r
+       TString outputSpectra = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Spectra.root";\r
+       \r
+       // configure the driver and run\r
+       \r
+       AliLnDriver driver;\r
+       \r
+       driver.SetSpecies(kSpecies);\r
+       \r
+       driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);\r
+       driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);\r
+       \r
+       driver.SetOutputTag(outputTag);\r
+       driver.SetTriggerEfficiency(trigEff);\r
+       driver.SetInelXSection(xsec);\r
+       driver.SetNormalizeToINEL(normToInel);\r
+       driver.SetVertexCorrection(kVtxCorr, kVtxCorrVal);\r
+       driver.SetPtBinInterval(kPtBin[0], kPtBin[1]);\r
+       driver.SetPidM2(kPidM2);\r
+       driver.SetM2BinInterval(kM2Bin[0], kM2Bin[1]);\r
+       driver.SetM2BkgInterval(kM2Bkg[0], kM2Bkg[1]);\r
+       driver.SetM2TPCInterval(kM2tpc[0], kM2tpc[1]);\r
+       driver.SetUnfolding(kUnfolding, kIter);\r
+       driver.SetFakeTracks(kFakeTracks);\r
+       driver.SetSecondaries(kSecondaries);\r
+       driver.SetSecProd(kSecProd);\r
+       driver.SetMatDCAxyModel(kMatDCAxyMod);\r
+       driver.SetNBin(kNbin);\r
+       driver.SetDCAxyInterval(kDCAxy[0], kDCAxy[1]);\r
+       driver.SetEfficiency(kEfficiency,0);\r
+       driver.SetFitFractionCorr(kFitFrac);\r
+       driver.SetSysErr(kSysErr[0],kSysErr[1]);\r
+       \r
+       driver.Run();\r
+       \r
+       // draw output\r
+       \r
+       if(!drawOutput) return 0;\r
+       \r
+       TStyle* st = GetDrawingStyle();\r
+       st->cd();\r
+       gROOT->ForceStyle();\r
+       \r
+       DrawOutputCorr(kSpecies,inputCorr);\r
+       \r
+       if(kSecProd != 2) gROOT->ProcessLine(Form(".x DrawSec.C+g(\"%s\",\"\",\"Deuteron\", %d, %d, %f, %f)", driver.GetPtCorrDebugFilename().Data(), kPtBin[0], kPtBin[1], kDCAxy[0], kDCAxy[1]));\r
+       \r
+       DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, kPidM2, kM2Bin[0], kM2Bin[1]);\r
+       DrawOutputRatio(outputRatio, outputTag, kSpecies);\r
+       DrawOutputSpectra(outputSpectra, outputTag, kSpecies);\r
+       \r
+       return 0;\r
+}\r
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TOF_LHC10x.C
new file mode 100644 (file)
index 0000000..31c7c2d
--- /dev/null
@@ -0,0 +1,125 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// LHC10x config for protons and antiprotons
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TString.h>
+#include "AliLnDriver.h"
+#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 normToInel         = 1,  // for mult
+                               Bool_t drawOutput         = 1)  // for batch
+{
+//
+// 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  kTrigName    = "mbor";
+       const Bool_t   kVtxCorr     = 0;
+       const Double_t kVtxCorrVal  = GetVertexCorrection(period);
+       const Int_t    kPtBin[2]    = {11,36};
+       const Bool_t   kUnfolding   = 0;
+       const Bool_t   kFakeTracks  = 0;
+       const Bool_t   kSecondaries = 1;
+       const Int_t    kSecProd     = 0; // 0 tff, 1 roofit, 2 mc
+       const Int_t    kMatDCAxyMod = 1; // 0 geant, 1 flat
+       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 Bool_t   kSameFdwn    = 1;
+       const Double_t kSysErr[2]   = {0.08,0.08} ;
+       
+       Double_t xsec[3];
+       GetInelXSection(xsec, period);
+       
+       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 outputPt      = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";
+       TString outputRatio   = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";
+       TString outputSpectra = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Spectra.root";
+       
+       // configure the driver and run
+       
+       AliLnDriver driver;
+       
+       driver.SetSpecies(kSpecies);
+       
+       driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);
+       driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
+       
+       driver.SetOutputTag(outputTag);
+       driver.SetTriggerEfficiency(trigEff);
+       driver.SetInelXSection(xsec);
+       driver.SetNormalizeToINEL(normToInel);
+       driver.SetVertexCorrection(kVtxCorr, kVtxCorrVal);
+       
+       driver.SetPtBinInterval(kPtBin[0], kPtBin[1]);
+       driver.SetPidM2(0);
+       driver.SetUnfolding(kUnfolding);
+       driver.SetFakeTracks(kFakeTracks);
+       driver.SetSecondaries(kSecondaries);
+       driver.SetSecProd(kSecProd);
+       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.Run();
+       
+       // draw output
+       
+       if(!drawOutput) return 0;
+       
+       TStyle* st = GetDrawingStyle();
+       st->cd();
+       gROOT->ForceStyle();
+       
+       DrawOutputCorr(kSpecies,inputCorr);
+       
+       if(kSecProd != 2) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, kPtBin[0], kPtBin[1], kDCAxy[0], kDCAxy[1]);
+       
+       DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0);
+       DrawOutputRatio(outputRatio, outputTag, kSpecies);
+       DrawOutputSpectra(outputSpectra, outputTag, kSpecies);
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPCTOF_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPCTOF_LHC10x.C
new file mode 100644 (file)
index 0000000..d216f3c
--- /dev/null
@@ -0,0 +1,224 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// LHC10x config for protons and antiprotons
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TSystem.h>
+#include <TString.h>
+#include <TFileMerger.h>
+
+#include "AliLnDriver.h"
+#include "Config.h"
+
+Int_t Config_Proton_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 normToInel         = 1,  // for mult
+                                  Bool_t drawOutput         = 1)  // for batch
+{
+//
+// lhc10b, lhc10c, lhc10d, lhc10e config for protons and antiprotons
+// (combine TPC and TOF)
+//
+       const TString  kSpecies      = "Proton";
+       
+       const TString  kTrkSelTPC    = "its_tpc_dca_spd-bayes";
+       const TString  kTrkSelTOF    = "its_tpc_tof_dca_spd-bayes";
+       
+       const TString  kOutputTagTPC = outputTag + "-tpc";
+       const TString  kOutputTagTOF = outputTag + "-tof";
+       
+       const Int_t    kJointBin     = 11;
+       const Int_t    kHiPtBin      = 36;
+       
+       const TString  kTrigName     = "mbor";
+       const Bool_t   kVtxCorr      = 0;
+       const Double_t kVtxCorrVal   = GetVertexCorrection(period);
+       
+       const Bool_t   kUnfolding    = 0;
+       const Bool_t   kFakeTracks   = 0;
+       const Bool_t   kSecondaries  = 1;
+       const Int_t    kSecProd      = 0; // 0 tff, 1 roofit, 2 mc
+       const Int_t    kMatDCAxyMod  = 1; // 0 geant, 1 flat
+       const Int_t    kNbin         = 10;
+       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 Bool_t   kSameFdwn     = 1;
+       
+       const Double_t kSysErr[2]    = {0.08, 0.08} ;
+       
+       Double_t xsec[3];
+       GetInelXSection(xsec, period);
+       
+       Double_t trigEff[3];
+       GetTriggerEfficiency(trigEff, kTrigName, period);
+       
+       // common options
+       
+       AliLnDriver driver;
+       
+       driver.SetSpecies(kSpecies);
+       
+       driver.SetTriggerEfficiency(trigEff);
+       driver.SetInelXSection(xsec);
+       driver.SetNormalizeToINEL(normToInel);
+       driver.SetVertexCorrection(kVtxCorr, kVtxCorrVal);
+       
+       driver.SetPidM2(0);
+       driver.SetUnfolding(kUnfolding);
+       driver.SetFakeTracks(kFakeTracks);
+       driver.SetSecondaries(kSecondaries);
+       driver.SetSecProd(kSecProd);
+       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]);
+       
+       // get the pt with TPC up to kJointBin
+       
+       TString inputDataTPC     = inputDir + "/" + period + "/"
+                                + MakeInputName(kSpecies, period, kTrkSelTPC+multTag) + ".root";
+       
+       TString inputSimuTPC     = inputDir + "/" + period + "/"
+                                + MakeSimuName(kSpecies, period, kTrkSelTPC+multCorTag) + ".root";
+       
+       TString inputSimuFixTPC  = inputDir + "/" + period + "/"
+                                + MakeSimuFixName(kSpecies, period, kTrkSelTPC+multCorTag, kG3Fluka) + ".root";
+       
+       TString inputCorrTPC     = inputDir + "/" + period + "/"
+                                + MakeInputName(kSpecies, period, kTrkSelTPC+multTag) + "-corr.root";
+       
+       TString outputPtTPC      = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTPC) + "-Pt.root";
+       TString outputRatioTPC   = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTPC) + "-Ratio.root";
+       TString outputSpectraTPC = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTPC) + "-Spectra.root";
+       
+       
+       driver.SetInputFilenames(inputDataTPC, inputSimuTPC, inputSimuFixTPC, inputCorrTPC);
+       driver.SetOutputFilenames(outputPtTPC, outputRatioTPC, outputSpectraTPC);
+       
+       driver.SetOutputTag(kOutputTagTPC);
+       
+       driver.SetMakeStats(0);
+       driver.SetMakeCorrections(1);
+       driver.SetMakePt(1);
+       driver.SetMakeRatio(0);
+       driver.SetMakeSpectra(0);
+       
+       driver.SetPtBinInterval(5, kJointBin);
+       
+       TString outputCorrDebugTPC = driver.GetPtCorrDebugFilename();
+       TString outputPtDebugTPC   = driver.GetPtDebugFilename();
+       
+       driver.Run();
+       
+       // get the pt with TOF from kJointBin to kHiPtBin
+       
+       TString inputDataTOF     = inputDir + "/" + period + "/"
+                                + MakeInputName(kSpecies, period, kTrkSelTOF+multTag) + ".root";
+       
+       TString inputSimuTOF     = inputDir + "/" + period + "/"
+                                + MakeSimuName(kSpecies, period, kTrkSelTOF+multCorTag) + ".root";
+       
+       TString inputSimuFixTOF  = inputDir + "/" + period + "/"
+                                + MakeSimuFixName(kSpecies, period, kTrkSelTOF+multCorTag) + ".root";
+       
+       TString inputCorrTOF     = inputDir + "/" + period + "/"
+                                + MakeInputName(kSpecies, period, kTrkSelTOF+multTag) + "-corr.root";
+       
+       TString outputPtTOF      = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTOF) + "-Pt.root";
+       TString outputRatioTOF   = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTOF) + "-Ratio.root";
+       TString outputSpectraTOF = outputDir + "/" + MakeOutputName(kSpecies, kOutputTagTOF) + "-Spectra.root";
+       
+       
+       driver.SetInputFilenames(inputDataTOF, inputSimuTOF, inputSimuFixTOF, inputCorrTOF);
+       driver.SetOutputFilenames(outputPtTOF, outputRatioTOF, outputSpectraTOF);
+       
+       driver.SetOutputTag(kOutputTagTOF);
+       
+       driver.SetMakeStats(1);
+       driver.SetMakeCorrections(1);
+       driver.SetMakePt(1);
+       driver.SetMakeRatio(0);
+       driver.SetMakeSpectra(0);
+       
+       driver.SetPtBinInterval(kJointBin, kHiPtBin);
+       driver.SetUnfolding(0);
+       driver.SetEfficiency(kEfficiency,0);
+       
+       TString outputCorrDebugTOF = driver.GetPtCorrDebugFilename();
+       TString outputPtDebugTOF   = driver.GetPtDebugFilename();
+       
+       driver.Run();
+       
+       // combine TPC and TOF pt
+       
+       TString outputPt      = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";
+       TString outputRatio   = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";
+       TString outputSpectra = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Spectra.root";
+       
+       TFileMerger m;
+       
+       m.AddFile(outputPtTPC.Data(),0);
+       m.AddFile(outputPtTOF.Data(),0);
+       
+       m.OutputFile(outputPt.Data());
+       
+       m.Merge();
+       
+       // make ratio and spectra
+       
+       driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
+       
+       driver.SetOutputTag(outputTag);
+       
+       driver.SetMakeCorrections(0);
+       driver.SetMakePt(0);
+       driver.SetMakeRatio(1);
+       driver.SetMakeSpectra(1);
+       
+       driver.Run();
+       
+       // delete tmp files
+       
+       gSystem->Exec(Form("rm -f %s %s %s %s", inputCorrTPC.Data(), outputPtTPC.Data(), outputCorrDebugTPC.Data(), outputPtDebugTPC.Data()));
+       gSystem->Exec(Form("rm -f %s %s %s %s", inputCorrTOF.Data(), outputPtTOF.Data(), outputCorrDebugTOF.Data(), outputPtDebugTOF.Data()));
+       
+       // draw output
+       
+       if(!drawOutput) return 0;
+       
+       TStyle* st = GetDrawingStyle();
+       st->cd();
+       gROOT->ForceStyle();
+       
+       DrawOutputRatio(outputRatio, outputTag, kSpecies);
+       DrawOutputSpectra(outputSpectra, outputTag, kSpecies);
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C b/PWGLF/SPECTRA/Nuclei/B2/macros/Config_Proton_TPC_LHC10x.C
new file mode 100644 (file)
index 0000000..2590623
--- /dev/null
@@ -0,0 +1,123 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// LHC10x config for protons and antiprotons
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TString.h>
+#include "AliLnDriver.h"
+#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 normToInel         = 1,  // for mult
+                               Bool_t drawOutput         = 1)  // for batch
+{
+//
+// lhc10b, lhc10c, lhc10d, lhc10e config for protons and antiprotons
+// (TPC)
+//
+       const TString  kSpecies     = "Proton";
+       const TString  kTrkSel      = "its_tpc_dca_spd-bayes";
+       const TString  kTrigName    = "mbor";
+       const Bool_t   kVtxCorr     = 0;
+       const Double_t kVtxCorrVal  = GetVertexCorrection(period);
+       const Int_t    kPtBin[2]    = {5,14};
+       const Bool_t   kUnfolding   = 0;
+       const Int_t    kIter        = 5;
+       const Bool_t   kFakeTracks  = 0;
+       const Bool_t   kSecondaries = 1;
+       const Int_t    kSecProd     = 0; // 0 tff, 1 roofit, 2 mc
+       const Int_t    kMatDCAxyMod = 1; // 0 geant, 1 flat
+       const Int_t    kNbin        = 10;
+       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 Bool_t   kSameFdwn    = 1;
+       const Double_t kSysErr[2]   = {0.08,0.08} ;
+       
+       Double_t xsec[3];
+       GetInelXSection(xsec, period);
+       
+       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 outputPt      = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Pt.root";
+       TString outputRatio   = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Ratio.root";
+       TString outputSpectra = outputDir + "/" + MakeOutputName(kSpecies, outputTag) + "-Spectra.root";
+       
+       // configure the driver and run
+       
+       AliLnDriver driver;
+       
+       driver.SetSpecies(kSpecies);
+       
+       driver.SetInputFilenames(inputData, inputSimu, inputSimuFix, inputCorr);
+       driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);
+       
+       driver.SetOutputTag(outputTag);
+       driver.SetTriggerEfficiency(trigEff);
+       driver.SetInelXSection(xsec);
+       driver.SetNormalizeToINEL(normToInel);
+       driver.SetVertexCorrection(kVtxCorr, kVtxCorrVal);
+       driver.SetPtBinInterval(kPtBin[0], kPtBin[1]);
+       driver.SetPidM2(0);
+       driver.SetUnfolding(kUnfolding, kIter);
+       driver.SetFakeTracks(kFakeTracks);
+       driver.SetSecondaries(kSecondaries);
+       driver.SetSecProd(kSecProd);
+       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.Run();
+       
+       // draw output
+       
+       if(!drawOutput) return 0;
+       
+       TStyle* st = GetDrawingStyle();
+       st->cd();
+       gROOT->ForceStyle();
+       
+       DrawOutputCorr(kSpecies,inputCorr);
+       
+       if(kSecProd != 2) DrawCorrDebug(driver.GetPtCorrDebugFilename(), driver.GetOutputCorrTag(), kSpecies, kPtBin[0], kPtBin[1], kDCAxy[0], kDCAxy[1]);
+       
+       DrawPtDebug(driver.GetPtDebugFilename(), outputTag, kSpecies, 0);
+       DrawOutputRatio(outputRatio, outputTag, kSpecies);
+       DrawOutputSpectra(outputSpectra, outputTag, kSpecies);
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawB2.C
new file mode 100644 (file)
index 0000000..5966121
--- /dev/null
@@ -0,0 +1,98 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw B2
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TGraphErrors.h>
+#include <TStyle.h>
+
+#include "B2.h"
+
+void DrawB2(const TString& inputFile="b2.root", const TString& tag="test", const TString& prefix="", const TString& species="Deuteron")
+{
+//
+// 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();
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       TString nucleus = prefix + species;
+       TString suffix = "";
+       if(prefix=="Anti") suffix="bar";
+       
+       TCanvas* c0 = new TCanvas(Form("%s.B2", nucleus.Data()), Form("Coalescence parameter for %s", nucleus.Data()));
+       
+       c0->Divide(2,1);
+       
+       // B2
+       
+       TGraphErrors* grB2Pt = (TGraphErrors*)FindObj(finput, tag, Form("B2%s_Pt",suffix.Data()));
+       TGraphErrors* grSysErrB2Pt = (TGraphErrors*)FindObj(finput, tag, Form("B2%s_SysErr_Pt",suffix.Data()));
+       
+       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->SetLineColor(kRed);
+       grB2Pt->SetMarkerStyle(kFullCircle);
+       grB2Pt->SetMarkerColor(kRed);
+       grB2Pt->Draw("P");
+       
+       // 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()));
+       
+       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->SetLineColor(kRed);
+       grR3Pt->SetMarkerStyle(kFullCircle);
+       grR3Pt->SetMarkerColor(kRed);
+       grR3Pt->Draw("P");
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawCorr.C
new file mode 100644 (file)
index 0000000..643e48e
--- /dev/null
@@ -0,0 +1,164 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw pt corrections for debugging
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TF1.h>
+
+#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 DrawCorr(const TString& species="Deuteron", const TString& inputFile="corrections.root", const TString& tag="")
+{
+//
+// Draw pt corrections for debugging
+//
+       Double_t xmin = 0;
+       Double_t xmax = 3.5;
+       
+       const Int_t kNpart = 2;
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       const TString kPrefix[] = {"", "Anti"};
+       
+       // Response matrix
+       
+       TCanvas* c0 = new TCanvas(Form("%s.Unfolding", species.Data()), Form("Response matrix for (Anti)%ss", species.Data()));
+       
+       c0->Divide(2,1);
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               c0->cd(i+1);
+               
+               TH2D* hResponseMtx = (TH2D*)FindObj(finput, tag, Form("%s%s_Response_Matrix",kPrefix[i].Data(), species.Data()));
+               hResponseMtx->SetAxisRange(xmin, xmax, "X");
+               hResponseMtx->SetAxisRange(xmin, xmax, "Y");
+               hResponseMtx->GetYaxis()->SetTitleOffset(1.3);
+               hResponseMtx->DrawCopy("cont0colZ");
+       }
+       
+       // Fake tracks
+       
+       TH1D* hFracFakeTracks[kNpart];
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               hFracFakeTracks[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Fake_Pt");
+       }
+       
+       TCanvas* c1 = new TCanvas(Form("%s.FakeTracks",species.Data()), Form("Fake Tracks for (Anti)%ss",species.Data()));
+       c1->cd();
+       
+       DrawPair(hFracFakeTracks[0], hFracFakeTracks[1], xmin, xmax, -0.01, 0.2);
+       
+       // Reconstruction efficiency
+       
+       TCanvas* c2 = new TCanvas(Form("%s.Efficiency",species.Data()), Form("Reconstruction Efficiency for (Anti)%ss",species.Data()));
+       c2->Divide(2,2);
+       
+       TH1D* hEffVtxPt[kNpart];
+       TH1D* hEffAccTrkPt[kNpart];
+       TH1D* hEffTrigPt[kNpart];
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               hEffTrigPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_Trig_Pt");
+               hEffVtxPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_Vtx_Pt");
+               hEffAccTrkPt[i] = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Eff_AccTrk_Pt");
+       }
+       
+       c2->cd(1);
+       DrawPair(hEffTrigPt[0], hEffTrigPt[1], xmin, xmax, 0, 1.1);
+       
+       c2->cd(2);
+       DrawPair(hEffVtxPt[0], hEffVtxPt[1], xmin, xmax, 0, 1.1);
+       
+       c2->cd(3);
+       DrawPair(hEffAccTrkPt[0], hEffAccTrkPt[1], xmin, xmax, 0, 1.1);
+       
+       // Secondaries
+       
+       TCanvas* c3 = new TCanvas(Form("%s.Secondaries",species.Data()), Form("Fraction of secondaries for (Anti)%ss",species.Data()));
+       
+       c3->Divide(2,2);
+       
+       for(Int_t i=0; i<kNpart; ++i)
+       {
+               TH1D* hFracFdwnPt = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Fdwn_Pt");
+               TH1D* hFracMatPt  = (TH1D*)FindObj(finput, tag, kPrefix[i] + species + "_Frac_Mat_Pt");
+               
+               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");
+               
+               c3->cd(2*i+1);
+               hFracFdwnPt->SetAxisRange(xmin, xmax, "X");
+               hFracFdwnPt->SetAxisRange(0., 0.6, "Y");
+               hFracFdwnPt->GetYaxis()->SetTitleOffset(1.4);
+               
+               hFracFdwnPt->SetMarkerColor(kBlue);
+               hFracFdwnPt->SetLineColor(kBlue);
+               hFracFdwnPt->SetMarkerStyle(kFullCircle);
+               hFracFdwnPt->DrawCopy("E");
+               
+               fncFracFdwnPt->SetLineColor(kRed);
+               fncFracFdwnPt->SetLineWidth(1);
+               fncFracFdwnPt->Draw("same");
+               
+               c3->cd(2*i+2);
+               hFracMatPt->SetAxisRange(xmin, xmax, "X");
+               hFracMatPt->SetAxisRange(0.,0.6, "Y");
+               hFracMatPt->GetYaxis()->SetTitleOffset(1.4);
+               
+               hFracMatPt->SetMarkerColor(kBlue);
+               hFracMatPt->SetLineColor(kBlue);
+               hFracMatPt->SetMarkerStyle(kFullCircle);
+               hFracMatPt->DrawCopy("E");
+               
+               fncFracMatPt->SetLineColor(kRed);
+               fncFracMatPt->SetLineWidth(1);
+               fncFracMatPt->Draw("same");
+       }
+}
+
+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)
+{
+//
+// Draw a pair of histograms in the current pad
+//
+       hX->SetTitle(title);
+       hX->SetAxisRange(xmin, xmax, "X");
+       hX->SetAxisRange(ymin, ymax, "Y");
+       hX->SetMarkerColor(xColor);
+       hX->SetLineColor(xColor);
+       hX->SetMarkerStyle(xMarker);
+       
+       hY->SetMarkerColor(yColor);
+       hY->SetLineColor(yColor);
+       hY->SetMarkerStyle(yMarker);
+       
+       hX->DrawCopy(option);
+       hY->DrawCopy(Form("same%s",option));
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawDir.C
new file mode 100644 (file)
index 0000000..82c57b7
--- /dev/null
@@ -0,0 +1,314 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw histograms/graphs with same name located in different directories
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TLegend.h>
+#include <TH1.h>
+#include <TGraphErrors.h>
+#include <TKey.h>
+#include <TLine.h>
+#include <TStyle.h>
+
+TGraphErrors* Divide(const TGraphErrors* grX, const TGraphErrors* grY, const TString& name);
+
+void DrawDir(const TString& inputFile,
+             const TString& graphName,
+             const TString& refdir="",
+             Double_t xmin = 0,
+             Double_t xmax = 3.5,
+             Double_t ymin = 0,
+             Double_t ymax = 1.,
+             const TString& xtitle = "",
+             const TString& ytitle = "",
+             Int_t option = 1,
+             const TString& canvasName = "c0",
+             const TString& canvasTitle = "DrawDir")
+{
+//
+// draw histograms/graphs with same name located in different directories
+// if option = 0 no comparison,
+// if option = 1 draw a bottom pad with the comparison, and
+// if option = 2 draw the comparasion in a different canvas
+//
+       const Int_t kMax = 10; // maximum number of subdirectories
+       
+       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 };
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       // iterate over the list of keys
+       
+       TObject* obj[kMax];
+       TString* subdir[kMax];
+       obj[0] = 0; // reference object
+       Int_t nDir = 0;
+       if(option>0) nDir = 1;
+       
+       TIter next(finput->GetListOfKeys());
+       TKey* key = 0;
+       while( (key = (TKey*)next()) && (nDir < kMax) )
+       {
+               TString classname = key->GetClassName();
+               if(classname == "TDirectoryFile")
+               {
+                       TObject* i = 0;
+                       finput->GetObject(Form("%s/%s;1", key->GetName(), graphName.Data()),i);
+                       if(i == 0)
+                       {
+                               finput->Error("GetObject","%s/%s not found", key->GetName(), graphName.Data());
+                               exit(1);
+                       }
+                       else if(i->InheritsFrom("TH1") || i->InheritsFrom("TGraph"))
+                       {
+                               Int_t j = 0;
+                               
+                               if(option==0) j = nDir++;
+                               else if(TString(key->GetName()) != refdir) j = nDir++;
+                               
+                               obj[j] = i;
+                               (dynamic_cast<TAttLine*>(i))->SetLineColor(kColor[j]);
+                               (dynamic_cast<TAttMarker*>(i))->SetMarkerColor(kColor[j]);
+                               (dynamic_cast<TAttMarker*>(i))->SetMarkerStyle(kMarker[j]);
+                               subdir[j] = new TString(key->GetName());
+                       }
+                       else
+                       {
+                               finput->Warning("GetObject", "%s does not contain %s", key->GetName(), graphName.Data());
+                       }
+               }
+       }
+       
+       // compare w.r.t. reference data
+       
+       TGraphErrors* grDiv[kMax];
+       
+       if(option > 0)
+       {
+               if(obj[0] == 0)
+               {
+                       finput->Error("GetObject", "reference directory %s not found", refdir.Data());
+                       for(Int_t i=0; i < nDir; ++i) delete subdir[i];
+                       exit(1);
+               }
+               
+               TGraphErrors* grY = 0;
+               
+               if(obj[0]->InheritsFrom("TH1"))
+               {
+                       grY = new TGraphErrors(dynamic_cast<TH1*>(obj[0]));
+               }
+               else if(obj[0]->InheritsFrom("TGraph"))
+               {
+                       grY = new TGraphErrors(*dynamic_cast<TGraphErrors*>(obj[0]));
+               }
+               
+               TGraphErrors* grX = 0;
+               
+               for(Int_t i=1; i < nDir; ++i)
+               {
+                       if(obj[i]->InheritsFrom("TH1"))
+                       {
+                               grX = new TGraphErrors(dynamic_cast<TH1*>(obj[i]));
+                       }
+                       else if(obj[i]->InheritsFrom("TGraph"))
+                       {
+                               grX = new TGraphErrors(*dynamic_cast<TGraphErrors*>(obj[i]));
+                       }
+                       
+                       grDiv[i] = Divide(grX, grY, Form("%s_Ratio", obj[i]->GetName()));
+                       
+                       grDiv[i]->SetLineColor(kColor[i]);
+                       grDiv[i]->SetMarkerColor(kColor[i]);
+                       grDiv[i]->SetMarkerStyle(kMarker[i]);
+                       
+                       delete grX;
+               }
+               
+               delete grY;
+       }
+       
+       // 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();
+       
+       TCanvas* c0 = new TCanvas(canvasName.Data(), canvasTitle.Data());
+       
+       if(option == 1) // create a top pad
+       {
+               TPad* top = new TPad("top", "Variable", 0, 0.25, 1, 1, 0, 0, 0);
+               
+               top->SetBottomMargin(0.);
+               top->Draw();
+               
+               top->cd();
+               
+               TH1F* frm = top->DrawFrame(xmin, ymin , xmax, ymax);
+               frm->GetYaxis()->SetTitle(ytitle);
+       }
+       else // only draw the frame
+       {
+               TH1F* frm = c0->DrawFrame(xmin, ymin , xmax, ymax);
+               frm->GetXaxis()->SetTitle(xtitle);
+               frm->GetYaxis()->SetTitle(ytitle);
+       }
+       
+       // draw objects in current pad
+       for(Int_t i=0; i < nDir; ++i)
+       {
+               if(obj[i]->InheritsFrom("TH1"))
+               {
+                       obj[i]->Draw("same");
+               }
+               else if(obj[i]->InheritsFrom("TGraph"))
+               {
+                       obj[i]->Draw("PZ");
+               }
+       }
+       
+       // build a legend
+       TLegend* legend = new TLegend(0.5718391,0.6991525,0.8390805,0.8368644,0,"brNDC");
+       legend->SetTextSize(0.03);
+       legend->SetFillColor(0);
+       legend->SetBorderSize(0);
+       
+       for(Int_t i=0; i < nDir; ++i)
+       {
+               legend->AddEntry(obj[i], subdir[i]->Data(), "lp");
+       }
+       
+       legend->Draw();
+       
+       if(option == 1) // create a bottom pad
+       {
+               c0->cd();
+               
+               TPad* bottom = new TPad("bottom", "ratio", 0, 0, 1, 0.25, 0, 0, 0);
+               
+               bottom->SetBottomMargin(0.3);
+               bottom->SetTopMargin(0.);
+               bottom->Draw();
+               
+               bottom->cd();
+               
+               TH1F* frm = bottom->DrawFrame(xmin, 0.7 , xmax, 1.3);
+               
+               frm->GetXaxis()->SetLabelSize(0.13);
+               frm->GetXaxis()->SetTitle(xtitle);
+               frm->GetXaxis()->SetTitleSize(0.13);
+               frm->GetYaxis()->SetTitle("Ratio");
+               frm->GetYaxis()->SetLabelSize(0.1);
+               frm->GetYaxis()->SetTitleSize(0.12);
+               frm->GetYaxis()->SetTitleOffset(0.3);
+       }
+       else if(option == 2) // create a new canvas
+       {
+               TCanvas* c1 = new TCanvas(Form("%s.Ratio",canvasName.Data()), Form("%s ratio",canvasTitle.Data()));
+               
+               TH1F* frm = c1->DrawFrame(xmin, 0.5 ,xmax, 1.5);
+               frm->GetXaxis()->SetTitle(xtitle);
+               frm->GetYaxis()->SetTitle("Ratio");
+       }
+       
+       // draw comparison
+       if(nDir>0 && option > 0)
+       {
+               for(Int_t i=1; i<nDir; ++i)
+               {
+                       grDiv[i]->Draw("PZ");
+               }
+               
+               // draw a red line for the reference
+               TLine* ref = new TLine(xmin,1,xmax,1);
+               ref->SetLineWidth(1);
+               ref->SetLineColor(kColor[0]);
+               ref->Draw();
+               
+               if(option == 2)
+               {
+                       TLegend* legendRatio = new TLegend(0.5718391,0.6991525,0.8390805,0.8368644,0,"brNDC");
+                       legendRatio->SetTextSize(0.03);
+                       legendRatio->SetFillColor(0);
+                       legendRatio->SetBorderSize(0);
+                       
+                       legendRatio->AddEntry(ref, subdir[0]->Data(), "l");
+                       for(Int_t i=1; i < nDir; ++i)
+                       {
+                               legendRatio->AddEntry(grDiv[i], subdir[i]->Data(), "lp");
+                       }
+                       
+                       legendRatio->Draw();
+               }
+       }
+       
+       // release memory
+}
+
+TGraphErrors* Divide(const TGraphErrors* grX, const TGraphErrors* grY, const TString& name)
+{
+//
+// Divide two TGraphErrors using first one as reference
+//
+       TGraphErrors* grQ = new TGraphErrors();
+       grQ->SetName(name.Data());
+       
+       for(Int_t i=0; i < grX->GetN(); ++i)
+       {
+               Double_t x, y1;
+               grX->GetPoint(i, x, y1);
+               Double_t y2 = grY->Eval(x);
+               
+               if(y1 == 0 || y2 == 0) continue;
+               
+               Double_t r = y1/y2;
+               
+               Double_t erry1 = grX->GetErrorY(i);
+               Double_t erry2 = grY->GetErrorY(i); // guess error
+               Double_t err = r*TMath::Sqrt(TMath::Power(erry1/y1,2)+TMath::Power(erry2/y2,2));
+               
+               grQ->SetPoint(i, x, r);
+               grQ->SetPointError(i, grX->GetErrorX(i), err);
+       }
+       
+       return grQ;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawPt.C
new file mode 100644 (file)
index 0000000..a326902
--- /dev/null
@@ -0,0 +1,166 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw corrected pt
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TGraphErrors.h>
+#include <TLegend.h>
+
+#include <RooWorkspace.h>
+#include <RooMsgService.h>
+#include <RooPlot.h>
+#include <RooRealVar.h>
+#include <RooAbsPdf.h>
+#include <RooAbsData.h>
+
+#include "B2.h"
+
+void DrawPt(const TString& inputFile="debug.root", const TString& tag="test", const TString& particle="AntiDeuteron", Bool_t m2pid=0, Int_t lowm2bin=9, Int_t him2bin=17)
+{
+//
+// Draw corrected pt for debugging
+//
+       Double_t xmin = 0;
+       Double_t xmax = 3.5;
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       // m2 data fitted models
+       
+       if(m2pid)
+       {
+               TH1D* hPidPt = (TH1D*)FindObj(finput, tag, Form("%s_PID_Pt",particle.Data()));
+               Double_t binwidth = hPidPt->GetBinWidth(0);
+               
+               using namespace RooFit;
+               
+               // disable verbose in RooFit
+               RooMsgService::instance().setSilentMode(1);
+               RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+               
+               TCanvas* c0 = new TCanvas(Form("%s.M2",particle.Data()), Form("M2 for %ss",particle.Data()));
+               c0->Divide(3,3);
+               
+               TGraphErrors* grFitSigmaPt = new TGraphErrors();
+               grFitSigmaPt->SetName(Form("%s_Fit_Sigma_Pt", particle.Data()));
+               
+               TGraphErrors* grFitM2Pt = new TGraphErrors();
+               grFitM2Pt->SetName(Form("%s_Fit_M2_Pt", particle.Data()));
+               
+               for(Int_t i=lowm2bin, j=0; i<him2bin && i-lowm2bin < 9 ; ++i)
+               {
+                       c0->cd(i-lowm2bin+1);
+                       //gPad->SetLogy(0);
+                       
+                       RooWorkspace* w= (RooWorkspace*)FindObj(finput, tag, Form("%s_M2_%02d",particle.Data(),i));
+                       
+                       RooPlot* m2frame = w->var("m2")->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("Bkg"))),LineWidth(1), LineColor(46),LineStyle(kDashed));
+                       
+                       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->SetMinimum(0.2);
+               
+                       m2frame->Draw();
+                       
+                       grFitSigmaPt->SetPoint(j, pt, w->var("width")->getVal());
+                       grFitSigmaPt->SetPointError(j, 0, w->var("width")->getError());
+                       
+                       grFitM2Pt->SetPoint(j, pt, w->var("mean")->getVal());
+                       grFitM2Pt->SetPointError(j++, 0, w->var("mean")->getError());
+               }
+               
+               // model parameters
+               
+               TCanvas* c1 = new TCanvas(Form("%s.M2.FitParameters",particle.Data()), Form("M2 model parameters for %ss",particle.Data()));
+               
+               c1->Divide(2,1);
+               
+               c1->cd(1);
+               gPad->SetGridx();
+               gPad->SetGridy();
+               
+               grFitSigmaPt->SetMarkerStyle(kFullCircle);
+               grFitSigmaPt->SetMarkerColor(kBlue);
+               grFitSigmaPt->SetLineColor(kBlue);
+               grFitSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+               grFitSigmaPt->GetYaxis()->SetTitle("#sigma_{m^{2}}");
+               grFitSigmaPt->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->Draw("ALP");
+               
+               //delete grFitSigmaPt;
+               //delete grFitM2Pt;
+       }
+       
+       // remaining corrections
+       
+       const Int_t kNum = 6;
+       const TString kCorr[kNum] = { "EffCor", "Unfolded", "SecCor", "FakeCor", "M2Corr", "PID"};
+       const TString kLabel[kNum] = { "Efficiency", "Unfolding", "Secondaries", "Fake tracks", "PID contamination", "Raw"};
+       const Int_t kColor[kNum]  = {kGreen-3, kGreen-2, kRed, kBlue, kOrange+1, kAzure};
+       const Int_t kMarker[kNum] = {kFullSquare, kFullCircle, kFullTriangleDown, kOpenCircle, kFullTriangleUp, kOpenSquare};
+       
+       TLegend* legend = new TLegend(0.5689655,0.6355932,0.8362069,0.8326271,0,"brNDC");
+       legend->SetTextSize(0.03);
+       legend->SetFillColor(0);
+       legend->SetBorderSize(0);
+       
+       TH1D* hPt[kNum];
+       
+       for(Int_t i=0; i<kNum; ++i)
+       {
+               hPt[i] = (TH1D*)FindObj(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[0]->SetTitle(particle.Data());
+       hPt[0]->SetAxisRange(xmin,xmax,"X");
+       hPt[0]->Draw("E");
+       for(Int_t i=1; i<kNum; ++i) hPt[i]->Draw("sameE");
+       
+       legend->Draw();
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawRatio.C
new file mode 100644 (file)
index 0000000..ef34ff5
--- /dev/null
@@ -0,0 +1,63 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw antiparticle/particle ratio
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TCanvas.h>
+
+#include "B2.h"
+
+void DrawRatio(const TString& inputFile="pt.root", const TString& tag="test", const TString& species="Deuteron")
+{
+//
+// Draw antiparticle/particle ratio
+//
+       Double_t xmin = 0;
+       Double_t xmax = 3.5;
+       Double_t ymin = 0.73;
+       Double_t ymax = 1.14;
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       // antiparticle/particle ratio
+       
+       TString ytitle = "Negative/Positive";
+       if(species=="Proton")        ytitle = "#bar{p}/p";
+       else if(species=="Deuteron") ytitle = "#bar{d}/d";
+       else if(species=="Triton")   ytitle = "#bar{t}/t";
+       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()));
+       
+       TCanvas* c0 = new TCanvas(Form("c0.Ratio%s",species.Data()), Form("Anti%s/%s ratio", species.Data(), species.Data()));
+       c0->cd();
+       
+       hRatioPt->SetTitle("");
+       hRatioPt->SetAxisRange(xmin,xmax,"X");
+       hRatioPt->SetAxisRange(ymin,ymax,"Y");
+       hRatioPt->SetLineColor(kRed);
+       hRatioPt->SetMarkerStyle(kFullCircle);
+       hRatioPt->SetMarkerColor(kRed);
+       hRatioPt->SetYTitle(ytitle);
+       hRatioPt->DrawCopy();
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSec.C
new file mode 100644 (file)
index 0000000..cc17cb7
--- /dev/null
@@ -0,0 +1,263 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw corrections of secondaries
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TString.h>
+#include <TMath.h>
+#include <TCanvas.h>
+
+#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* 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* 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);
+
+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)
+{
+//
+// draw corrections of secondaries
+//
+       const Int_t kNBin = hibin;
+       
+       TFile* finput = new TFile(inputFile.Data());
+       if (finput->IsZombie()) exit(1);
+       
+       // data and templates
+       
+       TH1D* hData[kNBin];
+       TH1D* hDataMC[kNBin];
+       TH1D* hPrim[kNBin];
+       TH1D* hMat[kNBin];
+       TH1D* hFdwn[kNBin];
+       
+       // fitted data
+       
+       TH1D* hDataFit[kNBin];
+       TH1D* hPrimFit[kNBin];
+       TH1D* hMatFit[kNBin];
+       TH1D* hFdwnFit[kNBin];
+       
+       for(Int_t i=lowbin; i<hibin && i<kNBin; ++i)
+       {
+               // MC
+               hData[i]   = (TH1D*)FindObj(finput, tag, particle + Form("_Data_DCAxy_%02d",i));
+               hDataMC[i] = (TH1D*)FindObj(finput, tag, particle + Form("_SimData_DCAxy_%02d",i));
+               hPrim[i]   = (TH1D*)FindObj(finput, tag, particle + Form("_Prim_DCAxy_%02d",i));
+               hMat[i]    = (TH1D*)FindObj(finput, tag, particle + Form("_Mat_DCAxy_%02d",i));
+               
+               if(particle.Contains("Proton"))
+               {
+                       hFdwn[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fdwn_DCAxy_%02d",i));
+               }
+               else
+               {
+                       hFdwn[i] = 0;
+               }
+               
+               // 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));
+               
+               if(particle.Contains("Proton"))
+               {
+                       hFdwnFit[i] = (TH1D*)FindObj(finput, tag, particle + Form("_Fit_Fdwn_DCAxy_%02d",i));
+               }
+               else
+               {
+                       hFdwnFit[i] = 0;
+               }
+       }
+       
+       // 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);
+       }
+       
+       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);
+       }
+       
+       TCanvas* c2;
+       
+       if(hDataFit[lowbin]!=0)
+       {
+               c2 = DrawGoF(hData, hDataFit, lowbin, hibin, particle + ".GoF", Form("GoF for %s",particle.Data()), binwidth, dcaxyMin, dcaxyMax);
+       }
+}
+
+Int_t GetNumberCols(Int_t lowbin, Int_t hibin)
+{
+//
+// number of rows and cols for plotting the slices
+// square matrix (cols = rows)
+//
+       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)
+{
+//
+// Draw MC DCA distributions for each pt bin
+//
+       Int_t ncol = GetNumberCols(lowbin, hibin);
+       
+       TCanvas* c0 = new TCanvas(name.Data(), title.Data());
+       c0->Divide(ncol, ncol);
+       
+       for(Int_t i=lowbin; i < hibin; ++i)
+       {
+               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]->SetMinimum(0.2);
+               
+               hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
+               
+               hData[i]->SetMarkerStyle(kFullCircle);
+               hData[i]->Draw("E");
+               
+               if(hPrim[i] != 0)
+               {
+                       hPrim[i]->SetLineColor(9);
+                       hPrim[i]->Draw("sameHist");
+               }
+               
+               if(hMat[i] != 0)
+               {
+                       hMat[i]->SetLineColor(46);
+                       hMat[i]->Draw("sameHist");
+               }
+               
+               if(hFdwn != 0)
+               if(hFdwn[i] != 0)
+               {
+                       hFdwn[i]->SetLineColor(8);
+                       hFdwn[i]->Draw("sameHist");
+               }
+       }
+       
+       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)
+{
+//
+// Draw DCA fitted model for each pt bin in current canvas
+//
+       Int_t ncol = GetNumberCols(lowbin, hibin);
+       
+       TCanvas* c0 = new TCanvas(name.Data(), title.Data());
+       c0->Divide(ncol, ncol);
+       
+       for(Int_t i=lowbin; i < hibin; ++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]->SetMinimum(0.2);
+               
+               hData[i]->SetAxisRange(dcaxyMin, dcaxyMax, "X");
+               
+               hData[i]->SetMarkerStyle(kFullCircle);
+               hData[i]->Draw("E");
+               
+               hDataFit[i]->SetLineColor(2);
+               hDataFit[i]->Draw("sameHist");
+               
+               if(hPrim[i] != 0)
+               {
+                       hPrim[i]->SetLineColor(9);
+                       hPrim[i]->Draw("sameHist");
+               }
+               
+               if(hMat[i] != 0)
+               {
+                       hMat[i]->SetLineColor(46);
+                       hMat[i]->Draw("sameHist");
+               }
+               
+               if(hFdwn != 0)
+               if(hFdwn[i] != 0)
+               {
+                       hFdwn[i]->SetLineColor(8);
+                       hFdwn[i]->Draw("sameHist");
+               }
+       }
+       
+       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)
+{
+//
+// Draw a goodness of fit for each pt bin consisting of data/fit
+//
+       Int_t ncol = GetNumberCols(lowbin, hibin);
+       
+       TCanvas* c0 = new TCanvas(name.Data(), title.Data());
+       c0->Divide(ncol, ncol);
+       
+       for(Int_t i=lowbin; i < hibin; ++i)
+       {
+               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->SetAxisRange(dcaxyMin, dcaxyMax, "X");
+               hDiv->SetAxisRange(0,3,"Y");
+               hDiv->SetMarkerStyle(kFullCircle);
+               hDiv->DrawCopy("E");
+               
+               delete hDiv;
+       }
+       
+       return c0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C b/PWGLF/SPECTRA/Nuclei/B2/macros/DrawSpectra.C
new file mode 100644 (file)
index 0000000..4c75bf7
--- /dev/null
@@ -0,0 +1,114 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// Draw differential yields
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <Riostream.h>
+#include <TROOT.h>
+#include <TStyle.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+
+#include "B2.h"
+
+void DrawSpectra(const TString& inputFile="spectra.root", const TString& tag="lhc10d", const TString& particle="Deuteron")
+{
+//
+// Draw differential yields
+//
+       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);
+       
+       // differential yields
+       
+       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");
+       
+       grDYieldPt->SetMarkerStyle(kFullCircle);
+       grDYieldPt->SetMarkerColor(kRed);
+       grDYieldPt->SetLineColor(kRed);
+       grDYieldPt->Draw("P");
+       
+       // fitting function
+       TF1* fit0 = (TF1*)FindObj(finput, tag, particle + "_Fit_DiffYield_Pt");
+       fit0->SetLineWidth(1);
+       fit0->SetLineColor(kRed);
+       fit0->Draw("same");
+       
+       // invariant differential yields
+       
+       c0->cd(2);
+       
+       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");
+       
+       grInvDYieldPt->SetMarkerStyle(kFullCircle);
+       grInvDYieldPt->SetMarkerColor(kRed);
+       grInvDYieldPt->SetLineColor(kRed);
+       grInvDYieldPt->Draw("P");
+       
+       // fitting function
+       TF1* fit1 = (TF1*)FindObj(finput, tag, particle + "_Fit_InvDiffYield_Pt");
+       fit1->SetLineWidth(1);
+       fit1->SetLineColor(kRed);
+       fit1->Draw("same");
+       
+       // invariant differential cross section
+       
+       c0->cd(3);
+       
+       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");
+       
+       // fitting function
+       TF1* fit2 = (TF1*)FindObj(finput, tag, particle + "_Fit_InvDiffXSection_Pt");
+       fit2->SetLineWidth(1);
+       fit2->SetLineColor(kRed);
+       fit2->Draw("same");
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10bcde.C
new file mode 100644 (file)
index 0000000..aa186d6
--- /dev/null
@@ -0,0 +1,150 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+// call Config_XXX for each period\r
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>\r
+\r
+#include <Riostream.h>\r
+#include <TROOT.h>\r
+#include <TSystem.h>\r
+#include <TString.h>\r
+#include <TFileMerger.h>\r
+\r
+#include "AliLnDriver.h"\r
+#include "Config.h"\r
+\r
+Int_t LHC10bcde(const TString& species    = "Deuteron",\r
+                const TString& inputDir   = "~/alice/input",\r
+                const TString& outputDir  = "~/alice/output",\r
+                const TString& outputTag  = "lhc10bcde",\r
+                const TString& multTag    = "",\r
+                const TString& multCorTag = "",\r
+                Bool_t normToInel         = 1,  // for mult\r
+                Bool_t drawOutput         = 1,  // for batch\r
+                Int_t option              = 0)\r
+{\r
+//\r
+// call Config_XXX for each period, merge the corrected pt and then get the results\r
+//\r
+       using namespace std;\r
+       \r
+       const Int_t kNper = 4;\r
+       const TString kPeriod[kNper]    = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
+       const TString kOutputTag[kNper] = { "lhc10b", "lhc10c", "lhc10d", "lhc10e" };\r
+       \r
+       Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;\r
+       Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;\r
+       \r
+       TFileMerger m;\r
+       \r
+       for(Int_t i=0; i<kNper; ++i)\r
+       {\r
+               cout << endl << "Period: " << kPeriod[i] << endl;\r
+               \r
+               TString arg =           inputDir       + "\","\r
+                             + "\""  + outputDir      + "\","\r
+                             + "\""  + kPeriod[i]     + "\","\r
+                             + "\""  + kOutputTag[i]  + "\","\r
+                             + "\""  + multTag        + "\","\r
+                             + "\""  + multCorTag;\r
+                       \r
+               if(species=="Proton" && option==0)\r
+               {\r
+                       cout << "Config_Proton_TPC_LHC10x.C" << endl << endl;\r
+                       gROOT->ProcessLine(Form(".x Config_Proton_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), normToInel));\r
+               }\r
+               else if(species=="Proton" && option==1)\r
+               {\r
+                       cout << "Config_Proton_TOF_LHC10x.C" << endl << endl;\r
+                       gROOT->ProcessLine(Form(".x Config_Proton_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), normToInel));\r
+               }\r
+               else if(species=="Proton" && option==2)\r
+               {\r
+                       cout << "Config_Proton_TPCTOF_LHC10x.C" << endl << endl;\r
+                       gROOT->ProcessLine(Form(".x Config_Proton_TPCTOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), normToInel));\r
+               }\r
+               else if(species=="Deuteron")\r
+               {\r
+                       cout << "Config_Deuteron_TOF_LHC10x.C" << endl << endl;\r
+                       gROOT->ProcessLine(Form(".x Config_Deuteron_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), normToInel));\r
+               }\r
+               else\r
+               {\r
+                       cerr << "unknown species/option: " << species << "/" << option << endl;\r
+                       cerr << "usage: Proton/0, Proton/1, Proton/2 or Deuteron" << endl;\r
+                       exit(1);\r
+               }\r
+               \r
+               TString ptfile = outputDir + "/" + species + "-" + kOutputTag[i] + "-Pt.root";\r
+               m.AddFile(ptfile,0);\r
+       }\r
+       \r
+       TString outputPt      = outputDir + "/" + species + "-" + outputTag + "-Pt.root";\r
+       TString outputRatio   = outputDir + "/" + species + "-" + outputTag + "-Ratio.root";\r
+       TString outputSpectra = outputDir + "/" + species + "-" + outputTag + "-Spectra.root";\r
+       \r
+       // pt\r
+       \r
+       m.OutputFile(outputPt.Data());\r
+       m.Merge();\r
+       \r
+       // spectra\r
+       \r
+       AliLnDriver driver;\r
+       \r
+       driver.SetSpecies(species);\r
+       \r
+       driver.SetOutputFilenames(outputPt, outputRatio, outputSpectra);\r
+       driver.SetOutputTag(outputTag);\r
+       \r
+       driver.SetMakeCorrections(0);\r
+       driver.SetMakePt(0);\r
+       driver.SetMakeRatio(1);\r
+       driver.SetMakeSpectra(1);\r
+       \r
+       driver.Run();\r
+       \r
+       // merge all results for comparison\r
+       \r
+       TFileMerger m2, m3;\r
+       \r
+       for(Int_t i=0; i<kNper; ++i)\r
+       {\r
+               m2.AddFile((outputDir + "/" + species + "-" + kPeriod[i] + "-Ratio.root").Data(),0);\r
+               m3.AddFile((outputDir + "/" + species + "-" + kPeriod[i] + "-Spectra.root").Data(),0);\r
+       }\r
+       \r
+       m2.AddFile(outputRatio.Data(),0);\r
+       m3.AddFile(outputSpectra.Data(),0);\r
+       \r
+       TString allRatios = outputDir + "/" + species + "-lhc10bcde1-Ratio.root";\r
+       TString allSpectra = outputDir + "/" + species + "-lhc10bcde1-Spectra.root";\r
+       \r
+       m2.OutputFile(allRatios.Data());\r
+       m3.OutputFile(allSpectra.Data());\r
+       \r
+       m2.Merge();\r
+       m3.Merge();\r
+       \r
+       // compare periods\r
+       \r
+       if(!drawOutput) return 0;\r
+       \r
+       gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"Anti%s%s_Ratio_Pt\",\"lhc10bcde\",0,4.5, 0., 1.8, \"p_{T} (GeV/c)\", \"neg/pos\", 1, \"cRatio\",\"Particle ratio\")", allRatios.Data(), species.Data(), species.Data()));\r
+       \r
+       DrawOutputSpectraMult(allSpectra, species, ymin, ymax, 1, "lhc10bcde");\r
+       \r
+       return 0;\r
+}\r
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/LHC10xMult.C
new file mode 100644 (file)
index 0000000..5c81f2f
--- /dev/null
@@ -0,0 +1,129 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// LHC10x multiplicity config for protons and deuterons
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TSystem.h>
+#include <TROOT.h>
+#include <TString.h>
+#include <TFileMerger.h>
+
+#include "AliLnDriver.h"
+#include "Config.h"
+
+Int_t LHC10xMult(const TString& species   = "Deuteron",
+                 const TString& inputDir  = "~/alice/input",
+                 const TString& outputDir = "~/alice/output",
+                 const TString& period    = "lhc10d",
+                 Int_t option = 2)
+{
+//
+// lhc10x multiplicity config
+// call Config_XXX for each multiplicity class
+//
+// if Proton and option = 0 then use Config_Proton_TPC
+// if Proton and option = 1 then use Config_Proton_TOF
+// if Proton and option = 2 then use Config_Proton_TPCTOF
+//
+       using namespace B2mult;
+       using namespace std;
+       
+       const Bool_t  kNormToInel[kNmult] = { 1, 0 }; // only normalize first bin
+       
+       Double_t ymin = (species=="Proton") ? 1.1e-6 : 1.1e-8;
+       Double_t ymax = (species=="Proton") ? 4.e-1 : 4.e-4;
+       
+       TFileMerger m1,m2;
+       
+       TString ratio[kNmult];
+       TString spectra[kNmult];
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               cout << endl;
+               cout << "Multiplicity class : " << kMultClass[i] << endl;
+               cout << "Period             : " << period << endl;
+               
+               TString outputTag  = period + "-" + kMultClass[i];
+               
+               TString arg =          inputDir            + "\","
+                             + "\"" + outputDir           + "\","
+                             + "\"" + period              + "\","
+                             + "\"" + outputTag           + "\","
+                             + "\"" + "-" + kMultClass[i] + "\"," // data
+                             + "\"" + "";                         // same simulations for all mult
+                       
+               if(species=="Proton" && option==0)
+               {
+                       cout << "Config_Proton_TPC_LHC10x.C" << endl << endl;
+                       gROOT->ProcessLine(Form(".x Config_Proton_TPC_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
+               }
+               else if(species=="Proton" && option==1)
+               {
+                       cout << "Config_Proton_TOF_LHC10x.C" << endl << endl;
+                       gROOT->ProcessLine(Form(".x Config_Proton_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
+               }
+               else if(species=="Proton" && option==2)
+               {
+                       cout << "Config_Proton_TPCTOF_LHC10x.C" << endl << endl;
+                       gROOT->ProcessLine(Form(".x Config_Proton_TPCTOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
+               }
+               else if(species=="Deuteron")
+               {
+                       cout << "Config_Deuteron_TOF_LHC10x.C" << endl << endl;
+                       gROOT->ProcessLine(Form(".x Config_Deuteron_TOF_LHC10x.C+g(\"%s\",%d,0)", arg.Data(), kNormToInel[i]));
+               }
+               else
+               {
+                       cerr << "unknown species/option: " << species << "/" << option << endl;
+                       cerr << "usage: Proton/0, Proton/1, Proton/2 or Deuteron" << endl;
+                       exit(1);
+               }
+               
+               ratio[i]   = outputDir + "/" + species + "-" + outputTag + "-Ratio.root";
+               spectra[i] = outputDir + "/" + species + "-" + outputTag + "-Spectra.root";
+               
+               m1.AddFile(ratio[i].Data(),0);
+               m2.AddFile(spectra[i].Data(),0);
+       }
+       
+       // merge
+       
+       TString allRatios  = outputDir + "/" + species + "-" + period + "-Mult-Ratio.root";
+       TString allSpectra = outputDir + "/" + species + "-" + period + "-Mult-Spectra.root";
+       
+       m1.OutputFile(allRatios.Data());
+       m2.OutputFile(allSpectra.Data());
+       
+       m1.Merge();
+       m2.Merge();
+       
+       // delete tmp files
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               gSystem->Exec(Form("rm -f %s", ratio[i].Data()));
+               gSystem->Exec(Form("rm -f %s", spectra[i].Data()));
+       }
+       
+       // 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()));
+       
+       DrawOutputSpectraMult(allSpectra, species, ymin, ymax);
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C b/PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C
new file mode 100644 (file)
index 0000000..d066cd8
--- /dev/null
@@ -0,0 +1,254 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// particle ratios as a function of multiplicity
+// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TStyle.h>
+
+#include "B2.h"
+#include "Config.h"
+
+void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr);
+void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle);
+
+void RatioMult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
+               const TString& ptag     = "lhc10d",
+               const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
+               const TString& dtag     = "lhc10d")
+{
+//
+// particle ratios as a function of multiplicity (from fitting distributions)
+//
+       using namespace B2mult;
+       
+       const Int_t kNpart = 2;
+       
+       const TString kProton[kNpart] = { "Proton", "AntiProton" };
+       const TString kDeuteron[kNpart] = { "Deuteron", "AntiDeuteron" };
+       
+       // open files
+       
+       TFile* fproton = new TFile(pSpectra.Data());
+       if(fproton->IsZombie()) exit(1);
+       
+       TFile* fdeuteron = new TFile(dSpectra.Data());
+       if(fdeuteron->IsZombie()) exit(1);
+       
+       // particle ratios
+       TGraphErrors*  grRatio[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grMixRatio[kNpart] = { new TGraphErrors(), new TGraphErrors() };
+       
+       // model parameters
+       TGraphErrors*  grProtondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grProtonN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grProtonC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grProtonChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
+       //TGraphErrors*  grProtonProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
+       
+       TGraphErrors*  grDeuterondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grDeuteronN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grDeuteronC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
+       TGraphErrors*  grDeuteronChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
+       //TGraphErrors*  grDeuteronProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
+       
+       for(Int_t i=0; i<kNmult; ++i)
+       {
+               TF1* fncProton[kNpart];
+               TF1* fncDeuteron[kNpart];
+               
+               for(Int_t j=0; j<kNpart; ++j)
+               {
+                       fncProton[j] = (TF1*)FindObj(fproton, ptag + "-" + kMultClass[i], kProton[j] + "_Fit_InvDiffYield_Pt");
+                       fncDeuteron[j] = (TF1*)FindObj(fdeuteron, dtag + "-" + kMultClass[i], kDeuteron[j] + "_Fit_InvDiffYield_Pt");
+               }
+               
+               // integrated ratios
+               
+               Double_t dNdy[kNpart] = {fncProton[0]->GetParameter(0), fncDeuteron[0]->GetParameter(0)};
+               Double_t dNdyErr[kNpart] = {fncProton[0]->GetParError(0), fncDeuteron[0]->GetParError(0)};
+               
+               Double_t dNdyBar[kNpart] = {fncProton[1]->GetParameter(0), fncDeuteron[1]->GetParameter(0)};
+               Double_t dNdyBarErr[kNpart] = {fncProton[1]->GetParError(0), fncDeuteron[1]->GetParError(0)};
+               
+               // ratios
+               Double_t ratio[kNpart], ratioErr[kNpart];
+               for(Int_t j=0; j<kNpart; ++j)
+               {
+                       GetRatio(dNdyBar[j], dNdy[j], dNdyBarErr[j], dNdyErr[j], ratio[j], ratioErr[j]);
+                       grRatio[j]->SetPoint(i, kKNOmult[i], ratio[j]);
+                       grRatio[j]->SetPointError(i, 0, ratioErr[j]);
+               }
+               
+               // mixed ratios
+               Double_t mixRatio[kNpart], mixRatioErr[kNpart];
+               
+               GetRatio(dNdy[1], dNdy[0], dNdyErr[1], dNdyErr[0], mixRatio[0], mixRatioErr[0]);
+               GetRatio(dNdyBar[1], dNdyBar[0], dNdyBarErr[1], dNdyBarErr[0], mixRatio[1], mixRatioErr[1]);
+               
+               for(Int_t j=0; j<kNpart; ++j)
+               {
+                       grMixRatio[j]->SetPoint(i, kKNOmult[i], mixRatio[j]);
+                       grMixRatio[j]->SetPointError(i, 0, mixRatioErr[j]);
+               }
+               
+               // model parameters
+               for(Int_t j=0; j<kNpart; ++j)
+               {
+                       grProtondNdy[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(0));
+                       grProtondNdy[j]->SetPointError(i, 0, fncProton[j]->GetParError(0));
+                       
+                       grProtonN[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(1));
+                       grProtonN[j]->SetPointError(i, 0, fncProton[j]->GetParError(1));
+                       
+                       grProtonC[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(2));
+                       grProtonC[j]->SetPointError(i, 0, fncProton[j]->GetParError(2));
+                       
+                       grProtonChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetChisquare()/fncProton[j]->GetNDF());
+                       //grProtonProb[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetProb());
+                       
+                       // deuteron
+                       grDeuterondNdy[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(0));
+                       grDeuterondNdy[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(0));
+                       
+                       grDeuteronN[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(1));
+                       grDeuteronN[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(1));
+                       
+                       grDeuteronC[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(2));
+                       grDeuteronC[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(2));
+                       
+                       grDeuteronChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetChisquare()/fncDeuteron[j]->GetNDF());
+                       //grDeuteronProb[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetProb());
+               }
+       }
+       
+       delete fproton;
+       delete fdeuteron;
+       
+       // 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");
+       
+       st->cd();
+       gROOT->ForceStyle();
+       
+       const Int_t kNCol = 4;
+       
+       TCanvas* c0 = new TCanvas("c0", "proton model parameters");
+       c0->Divide(kNCol,2);
+       
+       for(Int_t j=0; j<kNpart; ++j)
+       {
+               c0->cd(kNCol*j+1);
+               Draw(grProtondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
+               
+               c0->cd(kNCol*j+2);
+               Draw(grProtonN[j], kFullCircle, kBlue, "z", "n");
+               
+               c0->cd(kNCol*j+3);
+               Draw(grProtonC[j], kFullCircle, kBlue, "z", "C");
+               
+               c0->cd(kNCol*j+4);
+               Draw(grProtonChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
+               
+       /*      c0->cd(kNCol*j+5);
+               Draw(grProtonProb[j], kFullCircle, kBlue, "z", "Prob");
+       */
+       }
+       
+       TCanvas* c1 = new TCanvas("c1", "deuteron model parameters");
+       c1->Divide(kNCol,2);
+       
+       for(Int_t j=0; j<kNpart; ++j)
+       {
+               c1->cd(kNCol*j+1);
+               Draw(grDeuterondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
+               
+               c1->cd(kNCol*j+2);
+               Draw(grDeuteronN[j], kFullCircle, kBlue, "z", "n");
+               
+               c1->cd(kNCol*j+3);
+               Draw(grDeuteronC[j], kFullCircle, kBlue, "z", "C");
+               
+               c1->cd(kNCol*j+4);
+               Draw(grDeuteronChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
+               
+       /*      c1->cd(kNCol*j+5);
+               Draw(grDeuteronProb[j], kFullCircle, kBlue, "z", "Prob");
+       */
+       }
+       
+       TCanvas* c2 = new TCanvas("c2","particle ratios");
+       c2->Divide(2,2);
+       
+       c2->cd(1);
+       Draw(grRatio[0], kFullCircle, kRed, "z", "#bar{p}/p");
+       
+       c2->cd(2);
+       Draw(grRatio[1], kFullCircle, kRed, "z", "#bar{d}/d");
+       
+       c2->cd(3);
+       Draw(grMixRatio[0], kFullCircle, kRed, "z", "d/p");
+       
+       c2->cd(4);
+       Draw(grMixRatio[1], kFullCircle, kRed, "z", "#bar{d}/#bar{p}");
+}
+
+void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr)
+{
+//
+// get the ratio of x/y, its error and save in r and rerr
+//
+       if(x == 0 || y == 0)
+       {
+               r=0;
+               rerr=1e+16;
+               
+               return;
+       }
+       
+       r = x/y;
+       rerr = r*TMath::Sqrt(TMath::Power(errx/x,2)+TMath::Power(erry/y,2));
+}
+
+void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle)
+{
+//
+// draw a graph in current canvas
+//
+       gr->SetMarkerStyle(marker);
+       gr->SetMarkerColor(color);
+       gr->SetLineColor(color);
+       gr->GetXaxis()->SetTitle(xtitle.Data());
+       gr->GetYaxis()->SetTitle(ytitle.Data());
+       gr->Draw("AP");
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/run.C b/PWGLF/SPECTRA/Nuclei/B2/macros/run.C
new file mode 100644 (file)
index 0000000..a6cf8de
--- /dev/null
@@ -0,0 +1,34 @@
+Int_t run(const TString& config)
+{
+//
+// load libraries for the analysis
+// and run the config macro
+//
+       TStopwatch timer;
+       timer.Start();
+       
+       gROOT->SetMacroPath("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/");
+       
+       // 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/AliLnFakeTracks.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnCorr.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnRatio.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnSpectra.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnB2.cxx+g");
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/AliLnDriver.cxx+g");
+       
+       gROOT->ProcessLine(Form(".x %s", config.Data()));
+       
+       timer.Stop();
+       timer.Print();
+       
+       return 0;
+}
diff --git a/PWGLF/SPECTRA/Nuclei/B2/macros/runLocalExample.C b/PWGLF/SPECTRA/Nuclei/B2/macros/runLocalExample.C
new file mode 100644 (file)
index 0000000..48ce485
--- /dev/null
@@ -0,0 +1,91 @@
+void runLocalExample(const TString& dir="./")
+{
+//
+// example for testing the analysis task
+//
+       TStopwatch timer;
+       timer.Start();
+       
+       // Global configuration flags
+       
+       const TString kPeriod   = "lhc10d";
+       // with period name the next flags can be worked out
+       const Bool_t kSimulation = 0;
+       const Bool_t kHeavyIons  = 0;
+       
+       // Load common libraries
+       
+       gSystem->Load("libTree");
+       gSystem->Load("libGeom");
+       gSystem->Load("libVMC");
+       gSystem->Load("libPhysics");
+       gSystem->Load("libMinuit.so");
+       
+       gSystem->Load("libSTEERBase");
+       gSystem->Load("libESD");
+       gSystem->Load("libAOD");
+       
+       // Load analysis framework
+       
+       gSystem->Load("libANALYSIS.so");
+       gSystem->Load("libANALYSISalice.so"); // alice specific, e.g. alien plugin
+       
+       // Load task
+       
+       gSystem->Load("libPWGLFspectra.so");
+       
+       // Get input data
+       
+       gROOT->LoadMacro("$ALICE_ROOT/PWGUD/macros/CreateESDChain.C");
+       TChain* chain = CreateESDChain(dir,900);
+       
+       // Create the analysis manager
+       
+       AliAnalysisManager* mgr = new AliAnalysisManager("B2");
+       
+       // Input event handler (input container created automatically)
+       AliESDInputHandler* esdH = new AliESDInputHandler();
+       esdH->SetReadFriends(kFALSE);
+       mgr->SetInputEventHandler(esdH);
+       
+       // MonteCarlo handler
+       if(kSimulation)
+       {
+               AliMCEventHandler* mcHandler = new AliMCEventHandler();
+               mcHandler->SetReadTR(kFALSE); // read AliTrackReferences?
+               mgr->SetMCtruthEventHandler(mcHandler);
+       }
+       
+       // Create and add the task(s)
+       
+       // PhysicsSelection
+       gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+       
+       AliPhysicsSelectionTask* phySelectionTask = AddTaskPhysicsSelection(kSimulation);
+       
+       // Centrality selection
+       if(kHeavyIons)
+       {
+               gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+               AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
+       }
+       
+       // B2: Proton, Deuteron, Triton, ...
+       gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/Nuclei/B2/macros/AddTaskB2.C");
+       
+       // see TrackCuts.C macro for track selection criteria names
+       AliAnalysisTaskB2* proton = AddTaskB2("Proton","ProtonTest.root", "its_tpc_dca_spd", AliLnID::kBayes, kPeriod, kSimulation, kHeavyIons, 1);
+       
+       AliAnalysisTaskB2* deuteron = AddTaskB2("Deuteron","DeuteronTest.root", "its_tpc_tof_dca", AliLnID::kTPC, kPeriod, kSimulation, kHeavyIons, 0.2);
+       
+       // Run the analysis
+       
+       if (mgr->InitAnalysis())
+       {
+               mgr->PrintStatus();
+               mgr->StartAnalysis("local", chain);
+       }
+       
+       timer.Stop();
+       timer.Print();
+}