[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / macros / B2Mult.C

/**************************************************************************
 * 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>

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TSystem.h>
#include <TROOT.h>
#include <TFileMerger.h>
#include <TString.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include "AliLnB2.h"
#endif

#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;
}

extern void RatioMult(  const TString&, const TString&, const TString&, const TString&, const Int_t, const TString*, const Double_t*, const Double_t*, const TString&, const Bool_t, const TString& , const TString& );

Int_t B2Mult(  const TString& pSpectra     = "~/alice/output/Proton-lhc10d-Mult-Spectra.root"
             , const TString& ptag         = "lhc10d"
             , const TString& dSpectra     = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root"
             , const TString& dtag         = "lhc10d"
             , const TString& outputMultPt = "~/alice/output/B2-lhc10d-MultPt.root"
             , const TString& outputPtMult = "~/alice/output/B2-lhc10d-PtMult.root"
             , const TString& otag         = "lhc10d"
             , const Bool_t qTsallis       = 0
             , const TString& tsallisTag   = "Tsallis"
             , const TString& oratio       = "~/alice/output/Particle-Ratios-lhc10d-Tsallis.root")
{
//
// B2 as a function of multiplicity
//
	using namespace B2mult;
	using namespace std;
	
	const Int_t kNpart = 2;
	
	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";
			
			cout << kMultTag[i] << endl;
			
			AliLnB2 b2(pSpectra, ptag + kMultTag[i], dSpectra, dtag + kMultTag[i], b2file, otag + kMultTag[i], 2, kCharge[j]);
			
			b2.SetCd(GetCd(kKNOmult[i]));
			
			b2.Run();
			
			m.AddFile(b2file.Data(),0);
		}
		
		// merge B2 and B2bar
		
		TString outputfile = otag + kMultTag[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 + kMultTag[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(),kMultTag[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];
	TGraphErrors* grSysB2pt[kNpart][kNmult];
	TGraphErrors* grSysR3pt[kNpart][kNmult];
	
	for(Int_t i=0; i<kNpart; ++i)
	{
		for(Int_t j=0; j<kNmult; ++j)
		{
			grB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("B2%s_Pt", kSuffix[i].Data()));
			grR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("R3%s_Pt", kSuffix[i].Data()));
			
			grSysB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_B2%s_Pt", kSuffix[i].Data()));
			grSysR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_R3%s_Pt", kSuffix[i].Data()));
		}
	}
	
	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 B2SystErr[kNmult];
			
			Double_t R3[kNmult];
			Double_t R3StatErr[kNmult];
			Double_t R3SystErr[kNmult];
			
			for(Int_t k=0; k<kNmult; ++k)
			{
				Double_t x, y, staterr, systerr;
				grB2pt[i][k]->GetPoint(j,x,y);
				staterr = grB2pt[i][k]->GetErrorY(j);
				systerr = grSysB2pt[i][k]->GetErrorY(j);
				
				B2[k] = y;
				B2StatErr[k] = staterr;
				B2SystErr[k] = systerr;
				
				grR3pt[i][k]->GetPoint(j,x,y);
				staterr = grR3pt[i][k]->GetErrorY(j);
				systerr = grSysR3pt[i][k]->GetErrorY(j);
				
				R3[k] = y;
				R3StatErr[k] = staterr;
				R3SystErr[k] = systerr;
			}
			
			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()));
			
			Double_t zMultSystErr[kNmult];
			for(Int_t k=0; k<kNmult; ++k) zMultSystErr[k]=0.07;
			
			TGraphErrors* grSysB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, zMultSystErr, B2SystErr);
			grSysB2Mult->SetName(Form("SystErr_B2%s_Zmult", kSuffix[i].Data()));
			
			TGraphErrors* grSysR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, zMultSystErr, R3SystErr);
			grSysR3Mult->SetName(Form("SystErr_R3%s_Zmult", kSuffix[i].Data()));
			
			foutput->cd(ptLabel[j].Data());
			
			grB2Mult->Write();
			grR3Mult->Write();
			grSysB2Mult->Write();
			grSysR3Mult->Write();
			
			delete grB2Mult;
			delete grR3Mult;
			delete grSysB2Mult;
			delete grSysR3Mult;
		}
	}
	
	delete foutput;
	delete finput;
	
	//
	// Particle ratios
	// -----------------------------------
	
	RatioMult(pSpectra, dSpectra, ptag, dtag, kNmult, kMultTag, kKNOmult, kKNOmultErr, kKNOmultName, qTsallis, oratio, tsallisTag);
	
	
	// 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,5, 3.e-3, 6.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,5, 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;
}
Commit	Line	Data
2403d402	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// B2 as a function of multiplicity
	17	// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
	18
aa54def0	19	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
2403d402	20	#include <Riostream.h>
	21	#include <TSystem.h>
	22	#include <TROOT.h>
	23	#include <TFileMerger.h>
	24	#include <TString.h>
	25	#include <TFile.h>
	26	#include <TGraphErrors.h>
2403d402	27	#include "AliLnB2.h"
aa54def0	28	#endif
aa54def0	29
2403d402	30	#include "B2.h"
	31	#include "Config.h"
	32
	33	Double_t GetCd(Double_t z)
	34	{
	35	//
	36	// parameterization of <Cd> as a function of multiplicity
	37	// from ALICE Rlong and Rside measurements
	38	//
	39	return 0.046133 + 0.0484458*z;
	40	}
	41
aa54def0	42	extern void RatioMult( const TString&, const TString&, const TString&, const TString&, const Int_t, const TString, const Double_t, const Double_t*, const TString&, const Bool_t, const TString& , const TString& );
	43
	44	Int_t B2Mult( const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root"
	45	, const TString& ptag = "lhc10d"
	46	, const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root"
	47	, const TString& dtag = "lhc10d"
	48	, const TString& outputMultPt = "~/alice/output/B2-lhc10d-MultPt.root"
	49	, const TString& outputPtMult = "~/alice/output/B2-lhc10d-PtMult.root"
	50	, const TString& otag = "lhc10d"
	51	, const Bool_t qTsallis = 0
	52	, const TString& tsallisTag = "Tsallis"
	53	, const TString& oratio = "~/alice/output/Particle-Ratios-lhc10d-Tsallis.root")
2403d402	54	{
	55	//
	56	// B2 as a function of multiplicity
	57	//
	58	using namespace B2mult;
aa54def0	59	using namespace std;
2403d402	60
	61	const Int_t kNpart = 2;
	62
	63	const Int_t kNMinPt = 0;
	64	const Int_t kNMaxPt = 6;
	65
	66	const TString kPrefix[] = { "", "Anti"};
	67	const TString kSuffix[] = { "", "bar" };
	68	Int_t kCharge[] = {1, -1};
	69
	70	// B2 as a function of pt for each multiplicity class
	71
	72	for(Int_t i=0; i<kNmult; ++i)
	73	{
	74	TFileMerger m;
	75
	76	for(Int_t j=0; j<kNpart; ++j)
	77	{
	78	TString b2file = kPrefix[j] + "B2.root";
	79
aa54def0	80	cout << kMultTag[i] << endl;
	81
	82	AliLnB2 b2(pSpectra, ptag + kMultTag[i], dSpectra, dtag + kMultTag[i], b2file, otag + kMultTag[i], 2, kCharge[j]);
2403d402	83
	84	b2.SetCd(GetCd(kKNOmult[i]));
	85
	86	b2.Run();
	87
	88	m.AddFile(b2file.Data(),0);
	89	}
	90
	91	// merge B2 and B2bar
	92
aa54def0	93	TString outputfile = otag + kMultTag[i] + "-B2.root";
2403d402	94
	95	m.OutputFile(outputfile.Data());
	96	m.Merge();
	97
	98	gSystem->Exec("rm -f B2.root AntiB2.root");
	99	}
	100
	101	// merge multiplicity classes
	102
	103	TFileMerger m;
	104
	105	for(Int_t i=0; i<kNmult; ++i)
	106	{
aa54def0	107	TString b2 = otag + kMultTag[i] + "-B2.root";
2403d402	108	m.AddFile(b2.Data(),0);
	109	}
	110
	111	m.OutputFile(outputMultPt.Data());
	112	m.Merge();
	113
	114	// delete tmp files
	115
	116	for(Int_t i=0; i<kNmult; ++i)
	117	{
aa54def0	118	gSystem->Exec(Form("rm -f %s%s-B2.root",otag.Data(),kMultTag[i].Data()));
2403d402	119	}
	120
	121	// B2 as a function of multiplicity for each pt
	122
	123	TFile* finput = new TFile(outputMultPt.Data());
	124	if (finput->IsZombie()) exit(1);
	125
	126	TGraphErrors* grB2pt[kNpart][kNmult];
	127	TGraphErrors* grR3pt[kNpart][kNmult];
aa54def0	128	TGraphErrors* grSysB2pt[kNpart][kNmult];
aa54def0	129	TGraphErrors* grSysR3pt[kNpart][kNmult];
2403d402	130
	131	for(Int_t i=0; i<kNpart; ++i)
	132	{
	133	for(Int_t j=0; j<kNmult; ++j)
	134	{
aa54def0	135	grB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("B2%s_Pt", kSuffix[i].Data()));
	136	grR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("R3%s_Pt", kSuffix[i].Data()));
	137
	138	grSysB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_B2%s_Pt", kSuffix[i].Data()));
	139	grSysR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], Form("SystErr_R3%s_Pt", kSuffix[i].Data()));
2403d402	140	}
	141	}
	142
	143	TFile* foutput = new TFile(outputPtMult.Data(),"recreate");
	144
	145	Double_t* pt = grB2pt[0][0]->GetX();
	146	TString ptLabel[kNMaxPt];
	147
	148	if(kNMaxPt > grB2pt[0][0]->GetN())
	149	{
	150	std::cerr << "max pt too big" << std::endl;
	151	exit(1);
	152	}
	153
	154	for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
	155	{
aa54def0	156	ptLabel[i] = Form("pT %.02fA",pt[i]);
2403d402	157	foutput->mkdir(ptLabel[i].Data());
	158	}
	159
	160	for(Int_t i=0; i<kNpart; ++i)
	161	{
	162	for(Int_t j=kNMinPt; j<kNMaxPt; ++j)
	163	{
	164	Double_t B2[kNmult];
	165	Double_t B2StatErr[kNmult];
aa54def0	166	Double_t B2SystErr[kNmult];
2403d402	167
	168	Double_t R3[kNmult];
	169	Double_t R3StatErr[kNmult];
aa54def0	170	Double_t R3SystErr[kNmult];
2403d402	171
	172	for(Int_t k=0; k<kNmult; ++k)
	173	{
aa54def0	174	Double_t x, y, staterr, systerr;
2403d402	175	grB2pt[i][k]->GetPoint(j,x,y);
aa54def0	176	staterr = grB2pt[i][k]->GetErrorY(j);
aa54def0	177	systerr = grSysB2pt[i][k]->GetErrorY(j);
2403d402	178
2403d402	179	B2[k] = y;
aa54def0	180	B2StatErr[k] = staterr;
aa54def0	181	B2SystErr[k] = systerr;
2403d402	182
2403d402	183	grR3pt[i][k]->GetPoint(j,x,y);
aa54def0	184	staterr = grR3pt[i][k]->GetErrorY(j);
aa54def0	185	systerr = grSysR3pt[i][k]->GetErrorY(j);
2403d402	186
2403d402	187	R3[k] = y;
aa54def0	188	R3StatErr[k] = staterr;
aa54def0	189	R3SystErr[k] = systerr;
2403d402	190	}
	191
	192	TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
	193	grB2Mult->SetName(Form("B2%s_Zmult", kSuffix[i].Data()));
	194
	195	TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
	196	grR3Mult->SetName(Form("R3%s_Zmult", kSuffix[i].Data()));
	197
aa54def0	198	Double_t zMultSystErr[kNmult];
	199	for(Int_t k=0; k<kNmult; ++k) zMultSystErr[k]=0.07;
	200
	201	TGraphErrors* grSysB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, zMultSystErr, B2SystErr);
	202	grSysB2Mult->SetName(Form("SystErr_B2%s_Zmult", kSuffix[i].Data()));
	203
	204	TGraphErrors* grSysR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, zMultSystErr, R3SystErr);
	205	grSysR3Mult->SetName(Form("SystErr_R3%s_Zmult", kSuffix[i].Data()));
	206
2403d402	207	foutput->cd(ptLabel[j].Data());
	208
	209	grB2Mult->Write();
	210	grR3Mult->Write();
aa54def0	211	grSysB2Mult->Write();
aa54def0	212	grSysR3Mult->Write();
2403d402	213
	214	delete grB2Mult;
	215	delete grR3Mult;
aa54def0	216	delete grSysB2Mult;
aa54def0	217	delete grSysR3Mult;
2403d402	218	}
	219	}
	220
	221	delete foutput;
	222	delete finput;
	223
aa54def0	224	//
	225	// Particle ratios
	226	// -----------------------------------
	227
	228	RatioMult(pSpectra, dSpectra, ptag, dtag, kNmult, kMultTag, kKNOmult, kKNOmultErr, kKNOmultName, qTsallis, oratio, tsallisTag);
2403d402	229
2403d402	230
	231	// draw B2 as a function of pt
	232
	233	for(Int_t i=0; i<kNpart; ++i)
	234	{
	235	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()));
	236
	237	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()));
	238	}
	239
	240	// draw B2 as a function of z
	241
	242	for(Int_t i=0; i<kNpart; ++i)
	243	{
18bd338a	244	gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Zmult\",\"\",0,5, 3.e-3, 6.e-2,\"z\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2z\",\"B2%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
2403d402	245
18bd338a	246	gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Zmult\",\"\",0,5, 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()));
2403d402	247	}
	248
	249	return 0;
	250	}