[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / AliLnBA.cxx

/**************************************************************************
 * 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 "AliLnBA.h"
#include "B2.h"

ClassImp(AliLnBA)

AliLnBA::AliLnBA(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);
}

AliLnBA::~AliLnBA()
{
//
// destructor
//
}

void AliLnBA::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";
	}
	
	if(z<0) fNucleusName.Prepend("Anti");
}

Int_t AliLnBA::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);
	
	// invariant differential yields
	
	TString nucleonName = (fZ > 0) ? "Proton" : "AntiProton";
	
	TGraphErrors* grPrtInvDYieldPt = FindObj<TGraphErrors>(finput1, fProtonTag, nucleonName + "_InvDiffYield_Pt");
	TGraphErrors* grNucInvDYieldPt = FindObj<TGraphErrors>(finputA, fNucleusTag, fNucleusName + "_InvDiffYield_Pt");
	
	TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
	
	foutput->mkdir(fOutputTag.Data());
	foutput->cd(fOutputTag.Data());
	
	// coalescence parameter
	
	TGraphErrors* grBAPt = this->GetBAPt(grPrtInvDYieldPt, grNucInvDYieldPt, fNucleusName + Form("_B%d_Pt",fA));
	
	grBAPt->Write();
	
	// homogeneity volume
	TGraphErrors* grR3Pt = 0;
	if(fA==2)
	{
		grR3Pt = this->Rside2Rlong(grBAPt, fNucleusName + "_R3_Pt", fCd);
		grR3Pt->Write();
	}
	
	// propagate systematic errors if possible
	
	TString grPrtName = fProtonTag +  "/" + nucleonName + "_SystErr_InvDiffYield_Pt;1";
	TString grNucName = fNucleusTag + "/" + fNucleusName + "_SystErr_InvDiffYield_Pt;1";
	
	TGraphErrors* grSystErrPrtInvDYieldPt = dynamic_cast<TGraphErrors*>(finput1->Get(grPrtName.Data()));
	TGraphErrors* grSystErrNucInvDYieldPt = dynamic_cast<TGraphErrors*>(finputA->Get(grNucName.Data()));
	
	if( (grSystErrPrtInvDYieldPt != 0) && (grSystErrNucInvDYieldPt != 0) )
	{
		TGraphErrors* grSystErrBAPt = this->GetBAPt(grSystErrPrtInvDYieldPt, grSystErrNucInvDYieldPt, fNucleusName + Form("_SystErr_B%d_Pt",fA));
		
		if(fA==2)
		{
			TGraphErrors* grSystErrR3Pt = this->Rside2Rlong(grSystErrBAPt, fNucleusName + "_SystErr_R3_Pt", fCd);
			grSystErrR3Pt->Write();
			
			delete grSystErrR3Pt;
		}
		
		grSystErrBAPt->Write();
		delete grSystErrBAPt;
	}
	else
	{
		this->Warning("Run", "systematic errors are not propagated");
	}
	
	delete grBAPt;
	delete grR3Pt;
	
	delete foutput;
	delete finput1;
	delete finputA;
	
	return 0;
}

TGraphErrors* AliLnBA::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name) const
{
//
// coalescence parameter
//
	TGraphErrors* grBAPt = new TGraphErrors();
	grBAPt->SetName(name.Data());
	
	for(Int_t i=0, j=0; i < grNucInvDYieldPt->GetN(); ++i)
	{
		Double_t ptNuc, yNuc;
		
		grNucInvDYieldPt->GetPoint(i, ptNuc, yNuc);
		
		if(ptNuc/fA < 0.4) continue; // acceptance
		
		Double_t yPrt = grPrtInvDYieldPt->Eval(ptNuc/fA); // interpolate
		
		if(yPrt == 0 || yNuc == 0 ) continue;
		
		Double_t bA = yNuc/TMath::Power(yPrt,fA);
		
		// error
		Double_t ePrt = this->GetErrorY(grPrtInvDYieldPt, ptNuc/fA);
		Double_t eNuc = grNucInvDYieldPt->GetErrorY(i);
		
		Double_t errPt = grNucInvDYieldPt->GetErrorX(i)/fA;
		Double_t errBA = bA*TMath::Sqrt(TMath::Power(eNuc/yNuc,2) + TMath::Power(fA*ePrt/yPrt,2));
		
		grBAPt->SetPoint(j, ptNuc/fA, bA);
		grBAPt->SetPointError(j++, errPt, errBA);
	}
	
	return grBAPt;
}

Double_t AliLnBA::GetErrorY(const TGraphErrors* gr, Double_t x0) const
{
//
// estimate error of gr(x0) with the closest point to x0
//
	const Double_t kEpsilon  = 1.e-6;
	const Double_t kUnknownError = 1.e+6;
	
	for(Int_t i=0; i<gr->GetN(); ++i)
	{
		Double_t x = gr->GetX()[i];
		
		if( TMath::Abs(x-x0) < kEpsilon) return gr->GetErrorY(i);
		
		if( ((i == 0) && (x > x0)) || ((i == (gr->GetN()-1)) && (x < x0)) )
		{
			this->Warning("GetErrorY", "%f is out of bounds",x);
			return kUnknownError;
		}
		
		if( x > x0 )
		{
			this->Warning("GetErrorY", "Interpolating error at %f",x0);
			return (gr->GetErrorY(i)+gr->GetErrorY(i-1))/2;
		}
	}
	
	return 0;
}

Double_t AliLnBA::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* AliLnBA::Rside2Rlong(const TGraphErrors* grB2, const TString& name, 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)
//
	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 errPt = grB2->GetErrorX(i);
		Double_t errB2 = grB2->GetErrorY(i);
		Double_t errR = r3*errB2/b2;
		
		grR3->SetPointError(i, errPt, errR);
	}
	
	return grR3;
}
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	// coalescence parameter
	17	// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
	18
	19	#include <TMath.h>
	20	#include <TFile.h>
	21	#include <TString.h>
	22	#include <TGraphErrors.h>
	23	#include <TF1.h>
	24
e11ea0cb	25	#include "AliLnBA.h"
2403d402	26	#include "B2.h"
2403d402	27
e11ea0cb	28	ClassImp(AliLnBA)
2403d402	29
e11ea0cb	30	AliLnBA::AliLnBA(const TString& protonSpectra, const TString& protonTag, const TString& nucleusSpectra, const TString& nucleusTag, const TString& outputFilename, const TString& otag, Int_t a, Int_t z)
2403d402	31	: TObject()
	32	, fProtonSpectra(protonSpectra)
	33	, fProtonTag(protonTag)
	34	, fNucleusSpectra(nucleusSpectra)
	35	, fNucleusTag(nucleusTag)
	36	, fOutputFilename(outputFilename)
	37	, fOutputTag(otag)
	38	, fA(2)
	39	, fZ(1)
	40	, fNucleusName("Deuteron")
	41	, fCd(0.15)
	42	{
	43	//
	44	// constructor
	45	//
	46	this->SetNucleus(a,z);
	47	}
	48
e11ea0cb	49	AliLnBA::~AliLnBA()
2403d402	50	{
	51	//
	52	// destructor
	53	//
	54	}
	55
e11ea0cb	56	void AliLnBA::SetNucleus(Int_t a, Int_t z)
2403d402	57	{
	58	//
	59	// set nucleus mass and name
	60	//
	61	fA = a;
	62	fZ = z;
	63
	64	Int_t zz = TMath::Abs(z);
	65
	66	if(a==1 && zz==1) fNucleusName = "Proton";
	67	else if(a==2 && zz==1) fNucleusName = "Deuteron";
	68	else if(a==3 && zz==1) fNucleusName = "Triton";
	69	else if(a==3 && zz==2) fNucleusName = "He3";
	70	else if(a==4 && zz==2) fNucleusName = "Alpha";
	71	else
	72	{
	73	this->Warning("SetNucleus", "unknown nucleus A = %d, Z = %d", a, z);
	74	fNucleusName = "Unknown";
	75	}
e11ea0cb	76
e11ea0cb	77	if(z<0) fNucleusName.Prepend("Anti");
2403d402	78	}
2403d402	79
e11ea0cb	80	Int_t AliLnBA::Run()
2403d402	81	{
	82	//
	83	// coalescence parameter
	84	//
	85	TFile* finput1 = new TFile(fProtonSpectra.Data(), "read");
	86	if (finput1->IsZombie()) exit(1);
	87
	88	TFile* finputA = new TFile(fNucleusSpectra.Data(), "read");
	89	if (finputA->IsZombie()) exit(1);
	90
2403d402	91	// invariant differential yields
2403d402	92
e11ea0cb	93	TString nucleonName = (fZ > 0) ? "Proton" : "AntiProton";
	94
	95	TGraphErrors* grPrtInvDYieldPt = FindObj<TGraphErrors>(finput1, fProtonTag, nucleonName + "_InvDiffYield_Pt");
	96	TGraphErrors* grNucInvDYieldPt = FindObj<TGraphErrors>(finputA, fNucleusTag, fNucleusName + "_InvDiffYield_Pt");
2403d402	97
	98	TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
	99
	100	foutput->mkdir(fOutputTag.Data());
	101	foutput->cd(fOutputTag.Data());
	102
	103	// coalescence parameter
	104
e11ea0cb	105	TGraphErrors* grBAPt = this->GetBAPt(grPrtInvDYieldPt, grNucInvDYieldPt, fNucleusName + Form("_B%d_Pt",fA));
2403d402	106
e11ea0cb	107	grBAPt->Write();
2403d402	108
2403d402	109	// homogeneity volume
e11ea0cb	110	TGraphErrors* grR3Pt = 0;
	111	if(fA==2)
	112	{
	113	grR3Pt = this->Rside2Rlong(grBAPt, fNucleusName + "_R3_Pt", fCd);
	114	grR3Pt->Write();
	115	}
2403d402	116
aa54def0	117	// propagate systematic errors if possible
aa54def0	118
e11ea0cb	119	TString grPrtName = fProtonTag + "/" + nucleonName + "_SystErr_InvDiffYield_Pt;1";
e11ea0cb	120	TString grNucName = fNucleusTag + "/" + fNucleusName + "_SystErr_InvDiffYield_Pt;1";
aa54def0	121
aa54def0	122	TGraphErrors* grSystErrPrtInvDYieldPt = dynamic_cast<TGraphErrors*>(finput1->Get(grPrtName.Data()));
e11ea0cb	123	TGraphErrors* grSystErrNucInvDYieldPt = dynamic_cast<TGraphErrors*>(finputA->Get(grNucName.Data()));
aa54def0	124
	125	if( (grSystErrPrtInvDYieldPt != 0) && (grSystErrNucInvDYieldPt != 0) )
	126	{
e11ea0cb	127	TGraphErrors* grSystErrBAPt = this->GetBAPt(grSystErrPrtInvDYieldPt, grSystErrNucInvDYieldPt, fNucleusName + Form("_SystErr_B%d_Pt",fA));
aa54def0	128
e11ea0cb	129	if(fA==2)
	130	{
	131	TGraphErrors* grSystErrR3Pt = this->Rside2Rlong(grSystErrBAPt, fNucleusName + "_SystErr_R3_Pt", fCd);
	132	grSystErrR3Pt->Write();
	133
	134	delete grSystErrR3Pt;
	135	}
aa54def0	136
e11ea0cb	137	grSystErrBAPt->Write();
e11ea0cb	138	delete grSystErrBAPt;
aa54def0	139	}
	140	else
	141	{
	142	this->Warning("Run", "systematic errors are not propagated");
	143	}
2403d402	144
e11ea0cb	145	delete grBAPt;
2403d402	146	delete grR3Pt;
2403d402	147
	148	delete foutput;
	149	delete finput1;
	150	delete finputA;
	151
	152	return 0;
	153	}
	154
e11ea0cb	155	TGraphErrors* AliLnBA::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name) const
2403d402	156	{
	157	//
	158	// coalescence parameter
	159	//
	160	TGraphErrors* grBAPt = new TGraphErrors();
	161	grBAPt->SetName(name.Data());
	162
aa54def0	163	for(Int_t i=0, j=0; i < grNucInvDYieldPt->GetN(); ++i)
2403d402	164	{
aa54def0	165	Double_t ptNuc, yNuc;
aa54def0	166
2403d402	167	grNucInvDYieldPt->GetPoint(i, ptNuc, yNuc);
aa54def0	168
e11ea0cb	169	if(ptNuc/fA < 0.4) continue; // acceptance
aa54def0	170
aa54def0	171	Double_t yPrt = grPrtInvDYieldPt->Eval(ptNuc/fA); // interpolate
2403d402	172
2403d402	173	if(yPrt == 0 \|\| yNuc == 0 ) continue;
2403d402	174
	175	Double_t bA = yNuc/TMath::Power(yPrt,fA);
	176
	177	// error
aa54def0	178	Double_t ePrt = this->GetErrorY(grPrtInvDYieldPt, ptNuc/fA);
2403d402	179	Double_t eNuc = grNucInvDYieldPt->GetErrorY(i);
2403d402	180
aa54def0	181	Double_t errPt = grNucInvDYieldPt->GetErrorX(i)/fA;
2403d402	182	Double_t errBA = bATMath::Sqrt(TMath::Power(eNuc/yNuc,2) + TMath::Power(fAePrt/yPrt,2));
2403d402	183
aa54def0	184	grBAPt->SetPoint(j, ptNuc/fA, bA);
2403d402	185	grBAPt->SetPointError(j++, errPt, errBA);
	186	}
	187
	188	return grBAPt;
	189	}
	190
e11ea0cb	191	Double_t AliLnBA::GetErrorY(const TGraphErrors* gr, Double_t x0) const
aa54def0	192	{
	193	//
	194	// estimate error of gr(x0) with the closest point to x0
	195	//
	196	const Double_t kEpsilon = 1.e-6;
	197	const Double_t kUnknownError = 1.e+6;
	198
	199	for(Int_t i=0; i<gr->GetN(); ++i)
	200	{
	201	Double_t x = gr->GetX()[i];
	202
	203	if( TMath::Abs(x-x0) < kEpsilon) return gr->GetErrorY(i);
	204
	205	if( ((i == 0) && (x > x0)) \|\| ((i == (gr->GetN()-1)) && (x < x0)) )
	206	{
	207	this->Warning("GetErrorY", "%f is out of bounds",x);
	208	return kUnknownError;
	209	}
	210
	211	if( x > x0 )
	212	{
	213	this->Warning("GetErrorY", "Interpolating error at %f",x0);
	214	return (gr->GetErrorY(i)+gr->GetErrorY(i-1))/2;
	215	}
	216	}
	217
	218	return 0;
	219	}
	220
e11ea0cb	221	Double_t AliLnBA::Rside2Rlong(Double_t pt, Double_t B2, Double_t Cd) const
2403d402	222	{
	223	//
	224	// Rside^2*Rlong from B2 value
	225	// Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
	226	// (for pp ignore the exponential term)
	227	//
	228	Double_t hbarc = 0.1973269631; // GeV*fm
	229	Double_t m = 0.938272013; // GeV/c^2
	230	Double_t pi = TMath::Pi();
	231
	232	if(B2==0) return 0;
	233
	234	Double_t mt = TMath::Sqrt(ptpt + mm);
	235	Double_t r3 = 3.TMath::Power(pi,3./2.)TMath::Power(hbarc,3.)Cd/(2.mt*B2);
	236
	237	return r3;
	238	}
	239
e11ea0cb	240	TGraphErrors* AliLnBA::Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd) const
2403d402	241	{
	242	//
	243	// Rside^2*Rlong from B2 value
	244	// Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
	245	// (for pp ignore the exponential term)
	246	//
	247	TGraphErrors* grR3 = new TGraphErrors(*grB2);
	248	grR3->SetName(name.Data());
	249
	250	Double_t pt, b2;
	251
	252	for(Int_t i=0; i < grB2->GetN(); ++i)
	253	{
	254	grB2->GetPoint(i, pt, b2);
	255
	256	if(b2==0) continue;
	257	Double_t r3 = Rside2Rlong(pt, b2, Cd);
	258
	259	grR3->SetPoint(i, pt, r3);
	260
	261	// only statistical error propagation of B2
aa54def0	262	Double_t errPt = grB2->GetErrorX(i);
aa54def0	263	Double_t errB2 = grB2->GetErrorY(i);
2403d402	264	Double_t errR = r3*errB2/b2;
	265
	266	grR3->SetPointError(i, errPt, errR);
	267	}
	268
	269	return grR3;
	270	}