[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, "V0M0005", "V0M0510", "V0M0010");  
  // InterpolateRatios(321,211 , "V0M0005", "V0M0510", "V0M0010");  
  // 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");
  // ExtrapolateWithConstantRatioToPions(1000010020, "V0A6000", "V0A6080");
  // ExtrapolateWithConstantRatioToPions(-1000010020, "V0A6000", "V0A6080");


}

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) {

  // Get the ratios from the DB for the correct centralities
  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] = AliParticleYield::Add(ratio[0], ratio[1]);
  pi[0]    = AliParticleYield::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);
  cout << "*** "<< ratio[0]->GetPartName() << " " <<  centrfinal << " ***"<< std::endl;
  std::cout << "RATIO OF AVERAGE: " ;
  ratio[0]->Print("justvalue");
  //  average->Dump();
  std::cout << "AVERAGE OF RATIO: " ;
  average->Print("justvalue");

    
}

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);
4c7bd9d1	30	// InterpolateRatios(2212,211, "V0M0005", "V0M0510", "V0M0010");
4c7bd9d1	31	// InterpolateRatios(321,211 , "V0M0005", "V0M0510", "V0M0010");
ad431d97	32	// Interpolate6080(211);
	33	// Interpolate6080(-211);
	34	// Interpolate6080(2212);
	35	// Interpolate6080(-2212);
4c7bd9d1	36	// Interpolate6080(321);
4c7bd9d1	37	// Interpolate6080(-321);
ad431d97	38	// InterpolateRatios(2212,211, "V0M6070", "V0M7080", "V0M6080");
4c7bd9d1	39	// InterpolateRatios(321,211, "V0M6070", "V0M7080", "V0M6080");
ad431d97	40
9d84731e	41	// ************* Lambda and K0 ***************
	42	// arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt");
	43	// Interpolate0010(3122);
	44	// Interpolate0010(310);
	45	// ************* Helium 3 ***************
4c7bd9d1	46	arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
	47	arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
	48	ExtrapolateWithConstantRatioToPions(1000020030, "V0M0020", "V0M0010");
9d84731e	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");
4c7bd9d1	65	// ExtrapolateWithConstantRatioToPions(1000010020, "V0A6000", "V0A6080");
4c7bd9d1	66	// ExtrapolateWithConstantRatioToPions(-1000010020, "V0A6000", "V0A6080");
9d84731e	67
	68
	69	}
	70
	71	void DUMP(){
	72	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, 211, 2, 2760., "V0M0005");
	73	// p0 = AliParticleYield::Add(p0, AliParticleYield::FindParticle(arr, -211, 2, 2760., "V0M0005"));
	74	AliParticleYield * p2 = AliParticleYield::FindParticle(arr, 2212, 2, 2760., "V0M0005");
	75	// p2 = AliParticleYield::Add(p2,AliParticleYield::FindParticle(arr, -2212, 2, 2760., "V0M0005"));
	76	AliParticleYield *pratio = AliParticleYield::Divide(p2,p0);
	77	pratio->Print();
	78	AliParticleYield::FindRatio(arr, 2212, 211, 2, 2760, "V0M0005")->Print();
	79
	80	}
	81
	82	void Interpolate0010(Int_t pdg) {
	83
ad431d97	84	TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";
	85
	86	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0005");
	87	AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0510");
	88	p0->Scale(0.5);
	89	p1->Scale(0.5);
	90	AliParticleYield * pav = AliParticleYield::Add(p0,p1);
	91	std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;
	92
	93
	94	}
	95	void Interpolate6080(Int_t pdg) {
	96
	97	TString centrPrefix = collSystem == 2 ? "V0M" : "V0A";
	98
	99	AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"6070");
	100	AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"7080");
9d84731e	101	p0->Scale(0.5);
	102	p1->Scale(0.5);
	103	AliParticleYield * pav = AliParticleYield::Add(p0,p1);
	104	std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;
	105
	106
	107	}
	108
ad431d97	109	void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1, TString centr2, TString centrfinal) {
9d84731e	110
4c7bd9d1	111	// Get the ratios from the DB for the correct centralities
9d84731e	112	AliParticleYield * ratio[2];
ad431d97	113	ratio[0] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr1, 1);
ad431d97	114	ratio[1] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr2, 1);
4c7bd9d1	115
9d84731e	116	AliParticleYield * average = AliParticleYield::Add(ratio[0], ratio[1]);
9d84731e	117	average->Scale(0.5);
4c7bd9d1	118
9d84731e	119	AliParticleYield * pi[2];
ad431d97	120	pi[0] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr1, 0);
	121	pi[0] = AliParticleYield::Add(pi[0],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr1, 0));
	122	pi[1] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr2, 0);
	123	pi[1] = AliParticleYield::Add(pi[1],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr2, 0));
9d84731e	124
	125
	126	// Scale to get protons with errors corresponding to the ratio
	127	ratio[0]->Scale(pi[0]->GetYield()) ;
	128	ratio[1]->Scale(pi[1]->GetYield()) ;
	129
4c7bd9d1	130	ratio[0] = AliParticleYield::Add(ratio[0], ratio[1]);
4c7bd9d1	131	pi[0] = AliParticleYield::Add(pi[0],pi[1]);
9d84731e	132	pi[0]->SetNormError(0);
	133	pi[0]->SetStatError(0);
	134	pi[0]->SetSystError(0);
	135
	136	ratio[0]->Scale(1./pi[0]->GetYield());
ad431d97	137	ratio[0]->SetCentr(centrfinal);
4c7bd9d1	138	cout << "* "<< ratio[0]->GetPartName() << " " << centrfinal << " *"<< std::endl;
	139	std::cout << "RATIO OF AVERAGE: " ;
	140	ratio[0]->Print("justvalue");
9d84731e	141	// average->Dump();
4c7bd9d1	142	std::cout << "AVERAGE OF RATIO: " ;
4c7bd9d1	143	average->Print("justvalue");
9d84731e	144
	145
	146	}
	147
	148	void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest) {
	149
ad431d97	150	AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg,collSystem, energy, centrOrigin);
	151	AliParticleYield * pionOrigin = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrOrigin);
	152	AliParticleYield * pionDest = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrDest);
9d84731e	153	if(!part \|\| !pionOrigin \| !pionDest) {
	154	return;
	155	}
	156	// The pion yield is takes
	157	part->Scale (1./pionOrigin->GetYield());
	158	part->Scale (pionDest->GetYield());
	159	part->SetCentr(centrDest);
	160	part->SetMeasurementType(part->GetMeasurementType() \| AliParticleYield::kTypeExtrPionRatio);
	161	part->Print();
	162	// std::cout << "1" << std::endl;
	163	TClonesArray * arrOut = new TClonesArray("AliParticleYield");
	164	// std::cout << "2" << std::endl;
	165	new((arrOut)[0]) AliParticleYield(part) ;
	166
	167	// std::cout << "3" << std::endl;
	168	// std::cout << "4" << std::endl;
	169	AliParticleYield::SaveAsASCIIFile(arrOut, "temp.txt");
	170	}