[u/mrichter/AliRoot.git] / PWG4 / totEt / macros / CorrNeutralLevyFit.C

//Christine Nattrass, University of Tennessee at Knoxville
//This macro is to calculate the contributions to Et from various particles based on Levy fits to ALICE data at 900 GeV
//It uses d^2N/dpTdy from the papers, weighs it by Et, and integrates over all pT. 
//A=0 for mesons
//A=1 for antinucleons since they would anihilate in the calorimeter if they were observed by a calorimeter
//A=-1 for nucleons since they would not anihilate and their mass would not be measured
//At the end this is used to calculate the neutral energy correction
class AliAnalysisLevyPtModified{
public:
  virtual ~AliAnalysisLevyPtModified(){;};
  

  Double_t Evaluate(Double_t *pt, Double_t *par)
  {
    Double_t ldNdy  = par[0];
    Double_t l2pi   = 2*TMath::Pi();
    Double_t lTemp = par[1];
    Double_t lPower = par[2];
    Double_t lMass  = par[3];
    Double_t A = par[4];
    Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)+A*lMass;
    
    Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
    Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
    
    return ldNdy *Et* pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
  }
  ClassDef(AliAnalysisLevyPtModified, 1);
};

void CorrNeutralLevyFit(){

// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// $\pi^{+}+\pi^{-}$ & 2.977  $\pm$ 0.15  & 0.126 $\pm$ 0.001 & 7.82 $\pm$ 0.1  &                  &   &     \\
  AliAnalysisLevyPtModified *function = new AliAnalysisLevyPtModified();
  fPion = new TF1("fPion",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fPion->SetParameter(0,2.977);//dN/dy
  fPion->SetParameter(1,0.126);//T
  fPion->SetParameter(2,7.82);//n
  fPion->SetParameter(3,0.13957);//mass
  fPion->SetParameter(4,0);//A=0 for pions
  float integralPion = fPion->Integral(0,50);
  float integralErrPion = integralPion*0.15/2.977;
  float myerrorPionT = 0.0;
  float myerrorPionn = 0.0;
  float tmpPion;
  float T = 0.126;
  float Terr = 0.001;
  float n = 7.82;
  float nerr = 0.1;
  fPion->SetParameter(1,T+Terr);//T
  tmpPion = fPion->Integral(0,50);
  myerrorPionT = TMath::Abs(integralPion-tmpPion);
  fPion->SetParameter(1,T-Terr);//T
  tmpPion = fPion->Integral(0,50);
  if(TMath::Abs(integralPion-tmpPion)>myerrorPionT) myerrorPionT = TMath::Abs(integralPion-tmpPion);
  fPion->SetParameter(1,T);//T
  fPion->SetParameter(2,n+nerr);//n
  tmpPion = fPion->Integral(0,50);
  cout<<"integral "<<tmpPion<<endl;
  myerrorPionn = TMath::Abs(integralPion-tmpPion);
  fPion->SetParameter(2,n-nerr);//n
  tmpPion = fPion->Integral(0,50);
  cout<<"integral "<<tmpPion<<endl;
  if(TMath::Abs(integralPion-tmpPion)>myerrorPionn) myerrorPionn = TMath::Abs(integralPion-tmpPion);
  cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<"$\\pm$"<<myerrorPionT<<"$\\pm$"<<myerrorPionn<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrPion = TMath::Sqrt(TMath::Power(integralErrPion,2)+TMath::Power(myerrorPionT,2)+TMath::Power(myerrorPionn,2));
  cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<endl;

// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// $K^{+}+K^{-}$     & 0.366 $\pm$ 0.03   & 0.160 $\pm$ 0.006 & 6.087 $\pm$ 0.4 &                  &   &     \\
  fKaon = new TF1("fKaon",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fKaon->SetParameter(0,0.366);//dN/dy
  fKaon->SetParameter(1,0.160);//T
  fKaon->SetParameter(2,6.087);//n
  fKaon->SetParameter(3,0.493677);//mass
  fKaon->SetParameter(4,0);//A=0 for kaons
  float integralKaon = fKaon->Integral(0,50);
  float integralErrKaon = integralKaon*0.03/0.366;
  float myerrorKaonT = 0.0;
  float myerrorKaonn = 0.0;
  float tmpKaon;
  float T = 0.160;
  float Terr = 0.006;
  float n = 6.087;
  float nerr = 0.4;
  fKaon->SetParameter(1,T+Terr);//T
  tmpKaon = fKaon->Integral(0,50);
  myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
  fKaon->SetParameter(1,T-Terr);//T
  tmpKaon = fKaon->Integral(0,50);
  if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonT) myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
  fKaon->SetParameter(1,T);//T
  fKaon->SetParameter(2,n+nerr);//n
  tmpKaon = fKaon->Integral(0,50);
  cout<<"integral "<<tmpKaon<<endl;
  myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
  fKaon->SetParameter(2,n-nerr);//n
  tmpKaon = fKaon->Integral(0,50);
  cout<<"integral "<<tmpKaon<<endl;
  if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonn) myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
  cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<"$\\pm$"<<myerrorKaonT<<"$\\pm$"<<myerrorKaonn<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrKaon = TMath::Sqrt(TMath::Power(integralErrKaon,2)+TMath::Power(myerrorKaonT,2)+TMath::Power(myerrorKaonn,2));
  cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<endl;

// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// $p + \bar{p}$     &  0.157 $\pm$ 0.012 & 0.196 $\pm$ 0.009 & 8.6 $\pm$ 1.1   &                  &   &     \\
  fProton = new TF1("fProton",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fProton->SetParameter(0,0.157/2.0);//dN/dy
  fProton->SetParameter(1,0.196);//T
  fProton->SetParameter(2,8.6);//n
  fProton->SetParameter(3,0.938272);//mass
  fProton->SetParameter(4,-1);//A=0 for protons
  float integralProton = fProton->Integral(0,50);
  float integralErrProton = integralProton*0.012/0.157;
  float myerrorProtonT = 0.0;
  float myerrorProtonn = 0.0;
  float tmpProton;
  float T = 0.196;
  float Terr = 0.009;
  float n = 8.6;
  float nerr = 1.1;
  fProton->SetParameter(1,T+Terr);//T
  tmpProton = fProton->Integral(0,50);
  myerrorProtonT = TMath::Abs(integralProton-tmpProton);
  fProton->SetParameter(1,T-Terr);//T
  tmpProton = fProton->Integral(0,50);
  if(TMath::Abs(integralProton-tmpProton)>myerrorProtonT) myerrorProtonT = TMath::Abs(integralProton-tmpProton);
  fProton->SetParameter(1,T);//T
  fProton->SetParameter(2,n+nerr);//n
  tmpProton = fProton->Integral(0,50);
  cout<<"integral "<<tmpProton<<endl;
  myerrorProtonn = TMath::Abs(integralProton-tmpProton);
  fProton->SetParameter(2,n-nerr);//n
  tmpProton = fProton->Integral(0,50);
  cout<<"integral "<<tmpProton<<endl;
  if(TMath::Abs(integralProton-tmpProton)>myerrorProtonn) myerrorProtonn = TMath::Abs(integralProton-tmpProton);
  cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<"$\\pm$"<<myerrorProtonT<<"$\\pm$"<<myerrorProtonn<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrProton = TMath::Sqrt(TMath::Power(integralErrProton,2)+TMath::Power(myerrorProtonT,2)+TMath::Power(myerrorProtonn,2));
  cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<endl;


  //Antiprotons...
  fProton->SetParameter(2,n);//n
  fProton->SetParameter(4,1);//A=0 for protons
  float integralAntiProton = fProton->Integral(0,50);
  float integralErrAntiProton = integralAntiProton*0.012/0.157;
  float myerrorAntiProtonT = 0.0;
  float myerrorAntiProtonn = 0.0;
  float tmpAntiProton;
  float T = 0.196;
  float Terr = 0.009;
  float n = 8.6;
  float nerr = 1.1;
  fProton->SetParameter(1,T+Terr);//T
  tmpAntiProton = fProton->Integral(0,50);
  myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
  fProton->SetParameter(1,T-Terr);//T
  tmpAntiProton = fProton->Integral(0,50);
  if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonT) myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
  fProton->SetParameter(1,T);//T
  fProton->SetParameter(2,n+nerr);//n
  tmpAntiProton = fProton->Integral(0,50);
  cout<<"integral "<<tmpAntiProton<<endl;
  myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
  fProton->SetParameter(2,n-nerr);//n
  tmpAntiProton = fProton->Integral(0,50);
  cout<<"integral "<<tmpAntiProton<<endl;
  if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonn) myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
  cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<"$\\pm$"<<myerrorAntiProtonT<<"$\\pm$"<<myerrorAntiProtonn<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrAntiProton = TMath::Sqrt(TMath::Power(integralErrAntiProton,2)+TMath::Power(myerrorAntiProtonT,2)+TMath::Power(myerrorAntiProtonn,2));
  cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<endl;


// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// \kzeroshort       &0.184 $\pm$ 0.006   & 0.168 $\pm$ 0.005 & 6.6 $\pm$ 0.3   &                  &   &     \\
  fK0S = new TF1("fK0S",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fK0S->SetParameter(0,0.184);//dN/dy
  fK0S->SetParameter(1,0.168);//T
  fK0S->SetParameter(2,6.6);//n
  fK0S->SetParameter(3,.497614);//mass
  fK0S->SetParameter(4,0);//A=0 for kaons
  float integralK0S = fK0S->Integral(0,50);
  float integralErrK0S = integralK0S*0.006/0.184;
  float myerrorK0ST = 0.0;
  float myerrorK0Sn = 0.0;
  float tmpK0S;
  float T = 0.184;
  float Terr = 0.006;
  float n = 6.6;
  float nerr = 0.3;
  fK0S->SetParameter(1,T+Terr);//T
  tmpK0S = fK0S->Integral(0,50);
  myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
  fK0S->SetParameter(1,T-Terr);//T
  tmpK0S = fK0S->Integral(0,50);
  if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0ST) myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
  fK0S->SetParameter(1,T);//T
  fK0S->SetParameter(2,n+nerr);//n
  tmpK0S = fK0S->Integral(0,50);
  cout<<"integral "<<tmpK0S<<endl;
  myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
  fK0S->SetParameter(2,n-nerr);//n
  tmpK0S = fK0S->Integral(0,50);
  cout<<"integral "<<tmpK0S<<endl;
  if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0Sn) myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
  cout<<"K0S Et = "<<integralK0S<<"$\\pm$"<<integralErrK0S<<"$\\pm$"<<myerrorK0ST<<"$\\pm$"<<myerrorK0Sn<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrK0S = TMath::Sqrt(TMath::Power(integralErrK0S,2)+TMath::Power(myerrorK0ST,2)+TMath::Power(myerrorK0Sn,2));
  cout<<"K0S Et = "<<integralK0S<<"$\\pm$"<<integralErrK0S<<endl;

// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// \lam              &0.048 $\pm$ 0.004   & 0.229 $\pm$ 0.015 & 10.8 $\pm$ 2.0  &                  &   &     \\
  fLambda = new TF1("fLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fLambda->SetParameter(0,0.048);//dN/dy
  fLambda->SetParameter(1,0.229);//T
  fLambda->SetParameter(2,10.8);//n
  fLambda->SetParameter(3,1.115683);//mass
  fLambda->SetParameter(4,0);//A=0 for kaons
  float integralLambda = fLambda->Integral(0,50);
  float integralErrLambda = integralLambda*0.004/0.048;
  float myerrorLambdaT = 0.0;
  float myerrorLambdan = 0.0;
  float tmpLambda;
  float T = 0.229;
  float Terr = 0.015;
  float n = 10.8;
  float nerr = 2.0;
  fLambda->SetParameter(1,T+Terr);//T
  tmpLambda = fLambda->Integral(0,50);
  myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
  fLambda->SetParameter(1,T-Terr);//T
  tmpLambda = fLambda->Integral(0,50);
  if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdaT) myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
  fLambda->SetParameter(1,T);//T
  fLambda->SetParameter(2,n+nerr);//n
  tmpLambda = fLambda->Integral(0,50);
  cout<<"integral "<<tmpLambda<<endl;
  myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
  fLambda->SetParameter(2,n-nerr);//n
  tmpLambda = fLambda->Integral(0,50);
  cout<<"integral "<<tmpLambda<<endl;
  if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdan) myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
  cout<<"Lambda Et = "<<integralLambda<<"$\\pm$"<<integralErrLambda<<"$\\pm$"<<myerrorLambdaT<<"$\\pm$"<<myerrorLambdan<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrLambda = TMath::Sqrt(TMath::Power(integralErrLambda,2)+TMath::Power(myerrorLambdaT,2)+TMath::Power(myerrorLambdan,2));
  cout<<"Lambda Et = "<<integralLambda<<"$\\pm$"<<integralErrLambda<<endl;


// particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
// \alam             & 0.047 $\pm$ 0.005  & 0.210 $\pm$ 0.015 & 9.2 $\pm$ 1.4   &                  &   &     \\                  &   &     \\
  fAntiLambda = new TF1("fAntiLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
  fAntiLambda->SetParameter(0,0.047);//dN/dy
  fAntiLambda->SetParameter(1,0.210);//T
  fAntiLambda->SetParameter(2,9.2);//n
  fAntiLambda->SetParameter(3,1.115683);//mass
  fAntiLambda->SetParameter(4,0);//A=0 for kaons
  float integralAntiLambda = fAntiLambda->Integral(0,50);
  float integralErrAntiLambda = integralAntiLambda*0.005/0.047;
  float myerrorAntiLambdaT = 0.0;
  float myerrorAntiLambdan = 0.0;
  float tmpAntiLambda;
  float T = 0.210;
  float Terr = 0.015;
  float n = 9.2;
  float nerr = 1.4;
  fAntiLambda->SetParameter(1,T+Terr);//T
  tmpAntiLambda = fAntiLambda->Integral(0,50);
  myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
  fAntiLambda->SetParameter(1,T-Terr);//T
  tmpAntiLambda = fAntiLambda->Integral(0,50);
  if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdaT) myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
  fAntiLambda->SetParameter(1,T);//T
  fAntiLambda->SetParameter(2,n+nerr);//n
  tmpAntiLambda = fAntiLambda->Integral(0,50);
  cout<<"integral "<<tmpAntiLambda<<endl;
  myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
  fAntiLambda->SetParameter(2,n-nerr);//n
  tmpAntiLambda = fAntiLambda->Integral(0,50);
  cout<<"integral "<<tmpAntiLambda<<endl;
  if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdan) myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
  cout<<"AntiLambda Et = "<<integralAntiLambda<<"$\\pm$"<<integralErrAntiLambda<<"$\\pm$"<<myerrorAntiLambdaT<<"$\\pm$"<<myerrorAntiLambdan<<endl;
  //This isn't strictly correct because the errors on the parameters should be correlated but it's close
  integralErrAntiLambda = TMath::Sqrt(TMath::Power(integralErrAntiLambda,2)+TMath::Power(myerrorAntiLambdaT,2)+TMath::Power(myerrorAntiLambdan,2));
  cout<<"AntiLambda Et = "<<integralAntiLambda<<"$\\pm$"<<integralErrAntiLambda<<endl;

  //now we apply various assumptions
  float integralK0L = integralK0S;
  float integralErrK0L = integralErrK0S;
  float integralNeutron = integralProton;
  float integralErrNeutron = integralErrProton;

  float totalEt = integralPion+integralKaon+2.0*integralProton+2.0*integralAntiProton+2.0*integralK0S+integralLambda+integralAntiLambda;
  float measuredEt = integralPion+integralKaon+1.0*integralProton+1.0*integralAntiProton;
  float fneutral = measuredEt/totalEt;
  cout<<"fneutral = "<<fneutral<<endl;
  //this is from ugly derivative taking for error propagation
  //form 1:  pions, kaons, protons
  //for f=(xa+A)/(ya+B)
  //df/da=(xf-yf^2)/(xa+A)
  //x=y=1; xa+A=measuredEt
  float errPion = (fneutral-fneutral*fneutral)/measuredEt*integralErrPion;
  float errKaon = (fneutral-fneutral*fneutral)/measuredEt*integralErrKaon;
  //x=1,y=2
  float errProton = (fneutral-2.0*fneutral*fneutral)/measuredEt*integralErrProton;
  float errAntiProton = (fneutral-2.0*fneutral*fneutral)/measuredEt*integralErrAntiProton;
  //form 2: K0S, lambda, antilambda
  //for f=A/(B+xa)
  //df/da=-xf^2/A
  //x=1,A=measuredEt
  float errLambda = fneutral*fneutral/measuredEt*integralErrLambda;
  float errAntiLambda = fneutral*fneutral/measuredEt*integralErrAntiLambda;
  //x=2,A=measuredEt
  float errK0S = 2.0*fneutral*fneutral/measuredEt*integralErrK0S;
  cout<<"Errors Pion "<<errPion<<" Kaon "<<errKaon<<" Proton "<<errProton<<" Lambda "<<errLambda<<" AntiLambda "<<errAntiLambda<<" K0S "<<errK0S<<endl;
  float totalErr = TMath::Sqrt( TMath::Power(errPion,2) + TMath::Power(errKaon,2) + TMath::Power(errProton,2)+ TMath::Power(errAntiProton,2) + TMath::Power(errLambda,2) + TMath::Power(errAntiLambda,2) + TMath::Power(errK0S,2) );
  
  cout<<"fneutral = "<<fneutral<<"$\\pm$"<<totalErr<<endl;
  cout<<"1/fneutral = "<<1.0/fneutral<<"$\\pm$"<<1.0/fneutral*totalErr/fneutral<<endl;
  cout<<"Percentage error "<<totalErr/fneutral*100.0<<endl;
}
Commit	Line	Data
bc92a3aa	1	//Christine Nattrass, University of Tennessee at Knoxville
	2	//This macro is to calculate the contributions to Et from various particles based on Levy fits to ALICE data at 900 GeV
	3	//It uses d^2N/dpTdy from the papers, weighs it by Et, and integrates over all pT.
	4	//A=0 for mesons
	5	//A=1 for antinucleons since they would anihilate in the calorimeter if they were observed by a calorimeter
	6	//A=-1 for nucleons since they would not anihilate and their mass would not be measured
	7	//At the end this is used to calculate the neutral energy correction
	8	class AliAnalysisLevyPtModified{
	9	public:
	10	virtual ~AliAnalysisLevyPtModified(){;};
	11
	12
	13	Double_t Evaluate(Double_t pt, Double_t par)
	14	{
	15	Double_t ldNdy = par[0];
	16	Double_t l2pi = 2*TMath::Pi();
	17	Double_t lTemp = par[1];
	18	Double_t lPower = par[2];
	19	Double_t lMass = par[3];
	20	Double_t A = par[4];
	21	Double_t Et = TMath::Sqrt(pt[0]pt[0]+lMasslMass)+A*lMass;
	22
	23	Double_t lBigCoef = ((lPower-1)(lPower-2)) / (l2pilPowerlTemp(lPowerlTemp+lMass(lPower-2)));
	24	Double_t lInPower = 1 + (TMath::Sqrt(pt[0]pt[0]+lMasslMass)-lMass) / (lPower*lTemp);
	25
	26	return ldNdy Et pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
	27	}
	28	ClassDef(AliAnalysisLevyPtModified, 1);
	29	};
	30
	31	void CorrNeutralLevyFit(){
	32
	33	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
	34	// $\pi^{+}+\pi^{-}$ & 2.977 $\pm$ 0.15 & 0.126 $\pm$ 0.001 & 7.82 $\pm$ 0.1 & & & \\
	35	AliAnalysisLevyPtModified *function = new AliAnalysisLevyPtModified();
	36	fPion = new TF1("fPion",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
	37	fPion->SetParameter(0,2.977);//dN/dy
	38	fPion->SetParameter(1,0.126);//T
	39	fPion->SetParameter(2,7.82);//n
	40	fPion->SetParameter(3,0.13957);//mass
	41	fPion->SetParameter(4,0);//A=0 for pions
	42	float integralPion = fPion->Integral(0,50);
	43	float integralErrPion = integralPion*0.15/2.977;
	44	float myerrorPionT = 0.0;
	45	float myerrorPionn = 0.0;
	46	float tmpPion;
	47	float T = 0.126;
	48	float Terr = 0.001;
	49	float n = 7.82;
	50	float nerr = 0.1;
	51	fPion->SetParameter(1,T+Terr);//T
	52	tmpPion = fPion->Integral(0,50);
	53	myerrorPionT = TMath::Abs(integralPion-tmpPion);
	54	fPion->SetParameter(1,T-Terr);//T
	55	tmpPion = fPion->Integral(0,50);
	56	if(TMath::Abs(integralPion-tmpPion)>myerrorPionT) myerrorPionT = TMath::Abs(integralPion-tmpPion);
	57	fPion->SetParameter(1,T);//T
	58	fPion->SetParameter(2,n+nerr);//n
	59	tmpPion = fPion->Integral(0,50);
	60	cout<<"integral "<<tmpPion<<endl;
	61	myerrorPionn = TMath::Abs(integralPion-tmpPion);
	62	fPion->SetParameter(2,n-nerr);//n
	63	tmpPion = fPion->Integral(0,50);
	64	cout<<"integral "<<tmpPion<<endl;
65	if(TMath::Abs(integralPion-tmpPion)>myerrorPionn) myerrorPionn = TMath::Abs(integralPion-tmpPion);
66	cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<"$\\pm$"<<myerrorPionT<<"$\\pm$"<<myerrorPionn<<endl;
67	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
68	integralErrPion = TMath::Sqrt(TMath::Power(integralErrPion,2)+TMath::Power(myerrorPionT,2)+TMath::Power(myerrorPionn,2));
69	cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<endl;
70
71	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
72	// $K^{+}+K^{-}$ & 0.366 $\pm$ 0.03 & 0.160 $\pm$ 0.006 & 6.087 $\pm$ 0.4 & & & \\
73	fKaon = new TF1("fKaon",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
74	fKaon->SetParameter(0,0.366);//dN/dy
75	fKaon->SetParameter(1,0.160);//T
76	fKaon->SetParameter(2,6.087);//n
77	fKaon->SetParameter(3,0.493677);//mass
78	fKaon->SetParameter(4,0);//A=0 for kaons
79	float integralKaon = fKaon->Integral(0,50);
80	float integralErrKaon = integralKaon*0.03/0.366;
81	float myerrorKaonT = 0.0;
82	float myerrorKaonn = 0.0;
83	float tmpKaon;
84	float T = 0.160;
85	float Terr = 0.006;
86	float n = 6.087;
87	float nerr = 0.4;
88	fKaon->SetParameter(1,T+Terr);//T
89	tmpKaon = fKaon->Integral(0,50);
90	myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
91	fKaon->SetParameter(1,T-Terr);//T
92	tmpKaon = fKaon->Integral(0,50);
93	if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonT) myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
94	fKaon->SetParameter(1,T);//T
95	fKaon->SetParameter(2,n+nerr);//n
96	tmpKaon = fKaon->Integral(0,50);
97	cout<<"integral "<<tmpKaon<<endl;
98	myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
99	fKaon->SetParameter(2,n-nerr);//n
100	tmpKaon = fKaon->Integral(0,50);
101	cout<<"integral "<<tmpKaon<<endl;
102	if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonn) myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
103	cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<"$\\pm$"<<myerrorKaonT<<"$\\pm$"<<myerrorKaonn<<endl;
104	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
105	integralErrKaon = TMath::Sqrt(TMath::Power(integralErrKaon,2)+TMath::Power(myerrorKaonT,2)+TMath::Power(myerrorKaonn,2));
106	cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<endl;
107
108	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
109	// $p + \bar{p}$ & 0.157 $\pm$ 0.012 & 0.196 $\pm$ 0.009 & 8.6 $\pm$ 1.1 & & & \\
110	fProton = new TF1("fProton",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
111	fProton->SetParameter(0,0.157/2.0);//dN/dy
112	fProton->SetParameter(1,0.196);//T
113	fProton->SetParameter(2,8.6);//n
114	fProton->SetParameter(3,0.938272);//mass
115	fProton->SetParameter(4,-1);//A=0 for protons
116	float integralProton = fProton->Integral(0,50);
117	float integralErrProton = integralProton*0.012/0.157;
118	float myerrorProtonT = 0.0;
119	float myerrorProtonn = 0.0;
120	float tmpProton;
121	float T = 0.196;
122	float Terr = 0.009;
123	float n = 8.6;
124	float nerr = 1.1;
125	fProton->SetParameter(1,T+Terr);//T
126	tmpProton = fProton->Integral(0,50);
127	myerrorProtonT = TMath::Abs(integralProton-tmpProton);
128	fProton->SetParameter(1,T-Terr);//T
129	tmpProton = fProton->Integral(0,50);
130	if(TMath::Abs(integralProton-tmpProton)>myerrorProtonT) myerrorProtonT = TMath::Abs(integralProton-tmpProton);
131	fProton->SetParameter(1,T);//T
132	fProton->SetParameter(2,n+nerr);//n
133	tmpProton = fProton->Integral(0,50);
134	cout<<"integral "<<tmpProton<<endl;
135	myerrorProtonn = TMath::Abs(integralProton-tmpProton);
136	fProton->SetParameter(2,n-nerr);//n
137	tmpProton = fProton->Integral(0,50);
138	cout<<"integral "<<tmpProton<<endl;
139	if(TMath::Abs(integralProton-tmpProton)>myerrorProtonn) myerrorProtonn = TMath::Abs(integralProton-tmpProton);
140	cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<"$\\pm$"<<myerrorProtonT<<"$\\pm$"<<myerrorProtonn<<endl;
141	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
142	integralErrProton = TMath::Sqrt(TMath::Power(integralErrProton,2)+TMath::Power(myerrorProtonT,2)+TMath::Power(myerrorProtonn,2));
143	cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<endl;
144
145
146
147	//Antiprotons...
148	fProton->SetParameter(2,n);//n
149	fProton->SetParameter(4,1);//A=0 for protons
150	float integralAntiProton = fProton->Integral(0,50);
151	float integralErrAntiProton = integralAntiProton*0.012/0.157;
152	float myerrorAntiProtonT = 0.0;
153	float myerrorAntiProtonn = 0.0;
154	float tmpAntiProton;
155	float T = 0.196;
156	float Terr = 0.009;
157	float n = 8.6;
158	float nerr = 1.1;
159	fProton->SetParameter(1,T+Terr);//T
160	tmpAntiProton = fProton->Integral(0,50);
161	myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
162	fProton->SetParameter(1,T-Terr);//T
163	tmpAntiProton = fProton->Integral(0,50);
164	if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonT) myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
165	fProton->SetParameter(1,T);//T
166	fProton->SetParameter(2,n+nerr);//n
167	tmpAntiProton = fProton->Integral(0,50);
168	cout<<"integral "<<tmpAntiProton<<endl;
169	myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
170	fProton->SetParameter(2,n-nerr);//n
171	tmpAntiProton = fProton->Integral(0,50);
172	cout<<"integral "<<tmpAntiProton<<endl;
173	if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonn) myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
174	cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<"$\\pm$"<<myerrorAntiProtonT<<"$\\pm$"<<myerrorAntiProtonn<<endl;
175	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
176	integralErrAntiProton = TMath::Sqrt(TMath::Power(integralErrAntiProton,2)+TMath::Power(myerrorAntiProtonT,2)+TMath::Power(myerrorAntiProtonn,2));
177	cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<endl;
178
179
180	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
181	// \kzeroshort &0.184 $\pm$ 0.006 & 0.168 $\pm$ 0.005 & 6.6 $\pm$ 0.3 & & & \\
182	fK0S = new TF1("fK0S",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
183	fK0S->SetParameter(0,0.184);//dN/dy
184	fK0S->SetParameter(1,0.168);//T
185	fK0S->SetParameter(2,6.6);//n
186	fK0S->SetParameter(3,.497614);//mass
187	fK0S->SetParameter(4,0);//A=0 for kaons
188	float integralK0S = fK0S->Integral(0,50);
189	float integralErrK0S = integralK0S*0.006/0.184;
190	float myerrorK0ST = 0.0;
191	float myerrorK0Sn = 0.0;
192	float tmpK0S;
193	float T = 0.184;
194	float Terr = 0.006;
195	float n = 6.6;
196	float nerr = 0.3;
197	fK0S->SetParameter(1,T+Terr);//T
198	tmpK0S = fK0S->Integral(0,50);
199	myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
200	fK0S->SetParameter(1,T-Terr);//T
201	tmpK0S = fK0S->Integral(0,50);
202	if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0ST) myerrorK0ST = TMath::Abs(integralK0S-tmpK0S);
203	fK0S->SetParameter(1,T);//T
204	fK0S->SetParameter(2,n+nerr);//n
205	tmpK0S = fK0S->Integral(0,50);
206	cout<<"integral "<<tmpK0S<<endl;
207	myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
208	fK0S->SetParameter(2,n-nerr);//n
209	tmpK0S = fK0S->Integral(0,50);
210	cout<<"integral "<<tmpK0S<<endl;
211	if(TMath::Abs(integralK0S-tmpK0S)>myerrorK0Sn) myerrorK0Sn = TMath::Abs(integralK0S-tmpK0S);
212	cout<<"K0S Et = "<<integralK0S<<"$\\pm$"<<integralErrK0S<<"$\\pm$"<<myerrorK0ST<<"$\\pm$"<<myerrorK0Sn<<endl;
213	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
214	integralErrK0S = TMath::Sqrt(TMath::Power(integralErrK0S,2)+TMath::Power(myerrorK0ST,2)+TMath::Power(myerrorK0Sn,2));
215	cout<<"K0S Et = "<<integralK0S<<"$\\pm$"<<integralErrK0S<<endl;
216
217	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
218	// \lam &0.048 $\pm$ 0.004 & 0.229 $\pm$ 0.015 & 10.8 $\pm$ 2.0 & & & \\
219	fLambda = new TF1("fLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
220	fLambda->SetParameter(0,0.048);//dN/dy
221	fLambda->SetParameter(1,0.229);//T
222	fLambda->SetParameter(2,10.8);//n
223	fLambda->SetParameter(3,1.115683);//mass
224	fLambda->SetParameter(4,0);//A=0 for kaons
225	float integralLambda = fLambda->Integral(0,50);
226	float integralErrLambda = integralLambda*0.004/0.048;
227	float myerrorLambdaT = 0.0;
228	float myerrorLambdan = 0.0;
229	float tmpLambda;
230	float T = 0.229;
231	float Terr = 0.015;
232	float n = 10.8;
233	float nerr = 2.0;
234	fLambda->SetParameter(1,T+Terr);//T
235	tmpLambda = fLambda->Integral(0,50);
236	myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
237	fLambda->SetParameter(1,T-Terr);//T
238	tmpLambda = fLambda->Integral(0,50);
239	if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdaT) myerrorLambdaT = TMath::Abs(integralLambda-tmpLambda);
240	fLambda->SetParameter(1,T);//T
241	fLambda->SetParameter(2,n+nerr);//n
242	tmpLambda = fLambda->Integral(0,50);
243	cout<<"integral "<<tmpLambda<<endl;
244	myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
245	fLambda->SetParameter(2,n-nerr);//n
246	tmpLambda = fLambda->Integral(0,50);
247	cout<<"integral "<<tmpLambda<<endl;
248	if(TMath::Abs(integralLambda-tmpLambda)>myerrorLambdan) myerrorLambdan = TMath::Abs(integralLambda-tmpLambda);
249	cout<<"Lambda Et = "<<integralLambda<<"$\\pm$"<<integralErrLambda<<"$\\pm$"<<myerrorLambdaT<<"$\\pm$"<<myerrorLambdan<<endl;
250	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
251	integralErrLambda = TMath::Sqrt(TMath::Power(integralErrLambda,2)+TMath::Power(myerrorLambdaT,2)+TMath::Power(myerrorLambdan,2));
252	cout<<"Lambda Et = "<<integralLambda<<"$\\pm$"<<integralErrLambda<<endl;
253
254
255	// particle & $\frac{dN}{dy}$ & T (GeV) & n & $\frac{dE_T}{dy}$& a &\ET \\ \hline
256	// \alam & 0.047 $\pm$ 0.005 & 0.210 $\pm$ 0.015 & 9.2 $\pm$ 1.4 & & & \\ & & \\
257	fAntiLambda = new TF1("fAntiLambda",function, &AliAnalysisLevyPtModified::Evaluate,0,50,5,"AliAnalysisLevyPtModified","Evaluate");
258	fAntiLambda->SetParameter(0,0.047);//dN/dy
259	fAntiLambda->SetParameter(1,0.210);//T
260	fAntiLambda->SetParameter(2,9.2);//n
261	fAntiLambda->SetParameter(3,1.115683);//mass
262	fAntiLambda->SetParameter(4,0);//A=0 for kaons
263	float integralAntiLambda = fAntiLambda->Integral(0,50);
264	float integralErrAntiLambda = integralAntiLambda*0.005/0.047;
265	float myerrorAntiLambdaT = 0.0;
266	float myerrorAntiLambdan = 0.0;
267	float tmpAntiLambda;
268	float T = 0.210;
269	float Terr = 0.015;
270	float n = 9.2;
271	float nerr = 1.4;
272	fAntiLambda->SetParameter(1,T+Terr);//T
273	tmpAntiLambda = fAntiLambda->Integral(0,50);
274	myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
275	fAntiLambda->SetParameter(1,T-Terr);//T
276	tmpAntiLambda = fAntiLambda->Integral(0,50);
277	if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdaT) myerrorAntiLambdaT = TMath::Abs(integralAntiLambda-tmpAntiLambda);
278	fAntiLambda->SetParameter(1,T);//T
279	fAntiLambda->SetParameter(2,n+nerr);//n
280	tmpAntiLambda = fAntiLambda->Integral(0,50);
281	cout<<"integral "<<tmpAntiLambda<<endl;
282	myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
283	fAntiLambda->SetParameter(2,n-nerr);//n
284	tmpAntiLambda = fAntiLambda->Integral(0,50);
285	cout<<"integral "<<tmpAntiLambda<<endl;
286	if(TMath::Abs(integralAntiLambda-tmpAntiLambda)>myerrorAntiLambdan) myerrorAntiLambdan = TMath::Abs(integralAntiLambda-tmpAntiLambda);
287	cout<<"AntiLambda Et = "<<integralAntiLambda<<"$\\pm$"<<integralErrAntiLambda<<"$\\pm$"<<myerrorAntiLambdaT<<"$\\pm$"<<myerrorAntiLambdan<<endl;
288	//This isn't strictly correct because the errors on the parameters should be correlated but it's close
289	integralErrAntiLambda = TMath::Sqrt(TMath::Power(integralErrAntiLambda,2)+TMath::Power(myerrorAntiLambdaT,2)+TMath::Power(myerrorAntiLambdan,2));
290	cout<<"AntiLambda Et = "<<integralAntiLambda<<"$\\pm$"<<integralErrAntiLambda<<endl;
291
292	//now we apply various assumptions
293	float integralK0L = integralK0S;
294	float integralErrK0L = integralErrK0S;
295	float integralNeutron = integralProton;
296	float integralErrNeutron = integralErrProton;
297
298	float totalEt = integralPion+integralKaon+2.0integralProton+2.0integralAntiProton+2.0*integralK0S+integralLambda+integralAntiLambda;
299	float measuredEt = integralPion+integralKaon+1.0integralProton+1.0integralAntiProton;
300	float fneutral = measuredEt/totalEt;
301	cout<<"fneutral = "<<fneutral<<endl;
302	//this is from ugly derivative taking for error propagation
303	//form 1: pions, kaons, protons
304	//for f=(xa+A)/(ya+B)
305	//df/da=(xf-yf^2)/(xa+A)
306	//x=y=1; xa+A=measuredEt
307	float errPion = (fneutral-fneutralfneutral)/measuredEtintegralErrPion;
308	float errKaon = (fneutral-fneutralfneutral)/measuredEtintegralErrKaon;
309	//x=1,y=2
310	float errProton = (fneutral-2.0fneutralfneutral)/measuredEt*integralErrProton;
311	float errAntiProton = (fneutral-2.0fneutralfneutral)/measuredEt*integralErrAntiProton;
312	//form 2: K0S, lambda, antilambda
313	//for f=A/(B+xa)
314	//df/da=-xf^2/A
315	//x=1,A=measuredEt
316	float errLambda = fneutralfneutral/measuredEtintegralErrLambda;
317	float errAntiLambda = fneutralfneutral/measuredEtintegralErrAntiLambda;
318	//x=2,A=measuredEt
319	float errK0S = 2.0fneutralfneutral/measuredEt*integralErrK0S;
320	cout<<"Errors Pion "<<errPion<<" Kaon "<<errKaon<<" Proton "<<errProton<<" Lambda "<<errLambda<<" AntiLambda "<<errAntiLambda<<" K0S "<<errK0S<<endl;
321	float totalErr = TMath::Sqrt( TMath::Power(errPion,2) + TMath::Power(errKaon,2) + TMath::Power(errProton,2)+ TMath::Power(errAntiProton,2) + TMath::Power(errLambda,2) + TMath::Power(errAntiLambda,2) + TMath::Power(errK0S,2) );
322
323	cout<<"fneutral = "<<fneutral<<"$\\pm$"<<totalErr<<endl;
324	cout<<"1/fneutral = "<<1.0/fneutral<<"$\\pm$"<<1.0/fneutral*totalErr/fneutral<<endl;
325	cout<<"Percentage error "<<totalErr/fneutral*100.0<<endl;
326	}
327