[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / InterpolateRatiosAndYields.C

#if !defined (__CINT__) || (defined(__MAKECINT__))
#include <iostream>
#include "TClonesArray.h"
#include "AliParticleYield.h"

#endif
TClonesArray * arr =0;

void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1="V0M0005", TString centr2="V0M0510", TString centrfinal="V0M0010") ;
void Interpolate0010(Int_t pdg) ;
void Interpolate6080(Int_t pdg) ;

void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest);
Int_t collSystem = -1;
Float_t energy = -1;

void InterpolateRatiosAndYields() {
#if !(!defined (__CINT__) || (defined(__MAKECINT__)))
  LoadLibs();
#endif
  collSystem = 2; energy =2760;
  // *************** pi, K, pr *****************
  arr=   AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt");
  // Interpolate0010(211);
  // Interpolate0010(-211);
  // Interpolate0010(321);
  // Interpolate0010(-321);
  // Interpolate0010(2212);
  // Interpolate0010(-2212);
  //  InterpolateRatios(2212,211);  
  //  InterpolateRatios(321,211);  
  // Interpolate6080(211);
  // Interpolate6080(-211);
  // Interpolate6080(2212);
  // Interpolate6080(-2212);
  Interpolate6080(321);
  Interpolate6080(-321);
  // InterpolateRatios(2212,211, "V0M6070", "V0M7080", "V0M6080");  
  // InterpolateRatios(321,211, "V0M6070", "V0M7080", "V0M6080");    

  // *************** Lambda and K0 *****************
  // arr=   AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt");
  // Interpolate0010(3122);
  // Interpolate0010(310);
  // *************** Helium 3 *****************
  // arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
  // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
  // ExtrapolateWithConstantRatioToPions(1000020030, "V0M0020", "V0M0010");
  // *************** Kstar *****************
  // arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt");
  // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
  // ExtrapolateWithConstantRatioToPions(313, "V0M0020", "V0M0010");

  // *************** pPb, deuteron *********************
  //  collSystem = 1; energy = 5020;
  // 1. Average pions
  //arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt");
  //  Interpolate0010(211);
  //  Interpolate0010(-211);
  // 2. Extrapolate the deuteron
  // arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
  // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
  //  ExtrapolateWithConstantRatioToPions(1000010020, "V0A0010", "V0A0005");
  // ExtrapolateWithConstantRatioToPions(-1000010020, "V0A0010", "V0A0005");


}

void DUMP(){
  AliParticleYield * p0 = AliParticleYield::FindParticle(arr, 211,  2, 2760., "V0M0005");
  //  p0 = AliParticleYield::Add(p0, AliParticleYield::FindParticle(arr, -211,  2, 2760., "V0M0005"));
  AliParticleYield * p2 = AliParticleYield::FindParticle(arr, 2212, 2, 2760., "V0M0005");
  //  p2 = AliParticleYield::Add(p2,AliParticleYield::FindParticle(arr, -2212, 2, 2760., "V0M0005"));
  AliParticleYield *pratio = AliParticleYield::Divide(p2,p0);
  pratio->Print();
  AliParticleYield::FindRatio(arr, 2212, 211, 2, 2760, "V0M0005")->Print();

}

void Interpolate0010(Int_t pdg) {

  TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";

  AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0005");
  AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0510");
  p0->Scale(0.5);
  p1->Scale(0.5);
  AliParticleYield * pav = AliParticleYield::Add(p0,p1);
  std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;


} 
void Interpolate6080(Int_t pdg) {

  TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";

  AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"6070");
  AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"7080");
  p0->Scale(0.5);
  p1->Scale(0.5);
  AliParticleYield * pav = AliParticleYield::Add(p0,p1);
  std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;


} 

void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1, TString centr2, TString centrfinal) {

  AliParticleYield * ratio[2];
  ratio[0] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr1, 1);
  ratio[1] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr2, 1);
  AliParticleYield * average = AliParticleYield::Add(ratio[0], ratio[1]);
  average->Scale(0.5);
  AliParticleYield * pi[2];
  pi[0] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr1, 0);
  pi[0] = AliParticleYield::Add(pi[0],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr1, 0));
  pi[1] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr2, 0);
  pi[1] = AliParticleYield::Add(pi[1],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr2, 0));
  

  // Scale to get protons with errors corresponding to the ratio
  ratio[0]->Scale(pi[0]->GetYield()) ;
  ratio[1]->Scale(pi[1]->GetYield()) ;

  ratio[0]->Add(ratio[0], ratio[1]);
  pi[0]->Add(pi[0],pi[1]);
  pi[0]->SetNormError(0);
  pi[0]->SetStatError(0);
  pi[0]->SetSystError(0);
  
  ratio[0]->Scale(1./pi[0]->GetYield());
  ratio[0]->SetCentr(centrfinal);

  ratio[0]->Print();
  //  average->Dump();
  average->Print();

    
}

void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest) {

  AliParticleYield * part       =  AliParticleYield::FindParticle(arr, pdg,collSystem, energy, centrOrigin);
  AliParticleYield * pionOrigin =  AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrOrigin);
  AliParticleYield * pionDest   =  AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrDest);
  if(!part || !pionOrigin | !pionDest) {
    return;
  }
  // The pion yield is takes
  part->Scale (1./pionOrigin->GetYield());
  part->Scale (pionDest->GetYield());
  part->SetCentr(centrDest);
  part->SetMeasurementType(part->GetMeasurementType() | AliParticleYield::kTypeExtrPionRatio);
  part->Print();
  //  std::cout << "1" << std::endl;
  TClonesArray * arrOut = new TClonesArray("AliParticleYield");
  //  std::cout << "2" << std::endl;
  new((*arrOut)[0]) AliParticleYield(*part) ;

  //  std::cout << "3" << std::endl;
  //  std::cout << "4" << std::endl;
  AliParticleYield::SaveAsASCIIFile(arrOut, "temp.txt");
}
Commit	Line	Data
9d84731e	1	#if !defined (__CINT__) \|\| (defined(__MAKECINT__))
	2	#include <iostream>
	3	#include "TClonesArray.h"
	4	#include "AliParticleYield.h"
	5
	6	#endif
	7	TClonesArray * arr =0;
	8
ad431d97	9	void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1="V0M0005", TString centr2="V0M0510", TString centrfinal="V0M0010") ;
9d84731e	10	void Interpolate0010(Int_t pdg) ;
ad431d97	11	void Interpolate6080(Int_t pdg) ;
9d84731e	12
ad431d97	13	void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest);
	14	Int_t collSystem = -1;
	15	Float_t energy = -1;
9d84731e	16
	17	void InterpolateRatiosAndYields() {
	18	#if !(!defined (__CINT__) \|\| (defined(__MAKECINT__)))
	19	LoadLibs();
	20	#endif
ad431d97	21	collSystem = 2; energy =2760;
9d84731e	22	// ************* pi, K, pr ***************
ad431d97	23	arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt");
9d84731e	24	// Interpolate0010(211);
	25	// Interpolate0010(-211);
	26	// Interpolate0010(321);
	27	// Interpolate0010(-321);
	28	// Interpolate0010(2212);
	29	// Interpolate0010(-2212);
	30	// InterpolateRatios(2212,211);
	31	// InterpolateRatios(321,211);
ad431d97	32	// Interpolate6080(211);
	33	// Interpolate6080(-211);
	34	// Interpolate6080(2212);
	35	// Interpolate6080(-2212);
	36	Interpolate6080(321);
	37	Interpolate6080(-321);
	38	// InterpolateRatios(2212,211, "V0M6070", "V0M7080", "V0M6080");
	39	// InterpolateRatios(321,211, "V0M6070", "V0M7080", "V0M6080");
	40
9d84731e	41	// ************* Lambda and K0 ***************
	42	// arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt");
	43	// Interpolate0010(3122);
	44	// Interpolate0010(310);
	45	// ************* Helium 3 ***************
	46	// arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
	47	// arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
	48	// ExtrapolateWithConstantRatioToPions(1000020030, "V0M0020", "V0M0010");
	49	// ************* Kstar ***************
ad431d97	50	// arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt");
	51	// arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
	52	// ExtrapolateWithConstantRatioToPions(313, "V0M0020", "V0M0010");
	53
	54	// ************* pPb, deuteron *******************
	55	// collSystem = 1; energy = 5020;
	56	// 1. Average pions
	57	//arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt");
	58	// Interpolate0010(211);
	59	// Interpolate0010(-211);
	60	// 2. Extrapolate the deuteron
	61	// arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
	62	// arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
	63	// ExtrapolateWithConstantRatioToPions(1000010020, "V0A0010", "V0A0005");
	64	// ExtrapolateWithConstantRatioToPions(-1000010020, "V0A0010", "V0A0005");
9d84731e	65
	66
	67	}
	68
	69	void DUMP(){
	70	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, 211, 2, 2760., "V0M0005");
	71	// p0 = AliParticleYield::Add(p0, AliParticleYield::FindParticle(arr, -211, 2, 2760., "V0M0005"));
	72	AliParticleYield * p2 = AliParticleYield::FindParticle(arr, 2212, 2, 2760., "V0M0005");
	73	// p2 = AliParticleYield::Add(p2,AliParticleYield::FindParticle(arr, -2212, 2, 2760., "V0M0005"));
	74	AliParticleYield *pratio = AliParticleYield::Divide(p2,p0);
	75	pratio->Print();
	76	AliParticleYield::FindRatio(arr, 2212, 211, 2, 2760, "V0M0005")->Print();
	77
	78	}
	79
	80	void Interpolate0010(Int_t pdg) {
	81
ad431d97	82	TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";
	83
	84	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0005");
	85	AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0510");
	86	p0->Scale(0.5);
	87	p1->Scale(0.5);
	88	AliParticleYield * pav = AliParticleYield::Add(p0,p1);
	89	std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;
	90
	91
	92	}
	93	void Interpolate6080(Int_t pdg) {
	94
	95	TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";
	96
	97	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"6070");
	98	AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"7080");
9d84731e	99	p0->Scale(0.5);
	100	p1->Scale(0.5);
	101	AliParticleYield * pav = AliParticleYield::Add(p0,p1);
	102	std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;
	103
	104
	105	}
	106
ad431d97	107	void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1, TString centr2, TString centrfinal) {
9d84731e	108
9d84731e	109	AliParticleYield * ratio[2];
ad431d97	110	ratio[0] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr1, 1);
ad431d97	111	ratio[1] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr2, 1);
9d84731e	112	AliParticleYield * average = AliParticleYield::Add(ratio[0], ratio[1]);
	113	average->Scale(0.5);
	114	AliParticleYield * pi[2];
ad431d97	115	pi[0] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr1, 0);
	116	pi[0] = AliParticleYield::Add(pi[0],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr1, 0));
	117	pi[1] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr2, 0);
	118	pi[1] = AliParticleYield::Add(pi[1],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr2, 0));
9d84731e	119
	120
	121	// Scale to get protons with errors corresponding to the ratio
	122	ratio[0]->Scale(pi[0]->GetYield()) ;
	123	ratio[1]->Scale(pi[1]->GetYield()) ;
	124
	125	ratio[0]->Add(ratio[0], ratio[1]);
	126	pi[0]->Add(pi[0],pi[1]);
	127	pi[0]->SetNormError(0);
	128	pi[0]->SetStatError(0);
	129	pi[0]->SetSystError(0);
	130
	131	ratio[0]->Scale(1./pi[0]->GetYield());
ad431d97	132	ratio[0]->SetCentr(centrfinal);
9d84731e	133
	134	ratio[0]->Print();
	135	// average->Dump();
	136	average->Print();
	137
	138
	139	}
	140
	141	void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest) {
	142
ad431d97	143	AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg,collSystem, energy, centrOrigin);
	144	AliParticleYield * pionOrigin = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrOrigin);
	145	AliParticleYield * pionDest = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrDest);
9d84731e	146	if(!part \|\| !pionOrigin \| !pionDest) {
	147	return;
	148	}
	149	// The pion yield is takes
	150	part->Scale (1./pionOrigin->GetYield());
	151	part->Scale (pionDest->GetYield());
	152	part->SetCentr(centrDest);
	153	part->SetMeasurementType(part->GetMeasurementType() \| AliParticleYield::kTypeExtrPionRatio);
	154	part->Print();
	155	// std::cout << "1" << std::endl;
	156	TClonesArray * arrOut = new TClonesArray("AliParticleYield");
	157	// std::cout << "2" << std::endl;
	158	new((arrOut)[0]) AliParticleYield(part) ;
	159
	160	// std::cout << "3" << std::endl;
	161	// std::cout << "4" << std::endl;
	162	AliParticleYield::SaveAsASCIIFile(arrOut, "temp.txt");
	163	}