[u/mrichter/AliRoot.git] / PWGLF / totEt / macros / calculateCorrections.C

#include "TTree.h"
#include "TFile.h"
#include <TList.h>
#include <Rtypes.h>
#include "TCanvas.h"
#include <TH2I.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include "TF1.h"
#include <iostream>
#include <AliAnalysisEtTrackMatchCorrections.h>
#include <AliAnalysisEtRecEffCorrection.h>

TCanvas *c1 = 0;

TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
int createDummy(Double_t p0 = 0.366, Double_t p1 = 0.0)
{
  TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
  TF1 fitneutral("fitneutral","0", 0, 100);
  TF1 fitcharged("fitcharged","0", 0, 100);
  TF1 fitgamma("fitgamma","0", 0, 100);
  TF1 fitsecondary("fitsecondary","0", 0, 100);
  AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections("TmCorrectionsPhos", fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
  TF1 *func = generateRecEffFunction(p0, p1);
  AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection("ReCorrectionsPhos", *func, 1000);

  cor->Write();
  recor->Write();
  outfile->Close();


}
int calculateCorrections(TString filename="Et.ESD.simPbPb.PHOS.root", Double_t p0 = 0.366, Double_t p1 = 0.0)
{
  c1 = new TCanvas;
  
  TFile *f = TFile::Open(filename, "READ");
  
  TList *l = (TList*)f->Get("out1");
  
  TTree *primaryTree = (TTree*)l->FindObject("fPrimaryTreePhosMC");
  TTree *recTree = (TTree*)l->FindObject("fEventSummaryTreePhosRec");
  TTree *mcTree = (TTree*)l->FindObject("fEventSummaryTreePhosMC");
  std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl;
 
  Int_t clusterMult = 0;
  Int_t nChargedNotRemoved = 0;
  Int_t nNeutralNotRemoved = 0;
  Int_t nGammaRemoved = 0;
  Int_t nSecNotRemoved = 0;
  
  Double_t emEtRec = 0;
  Double_t emEtMc = 0;

  recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult);
  mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved);
  mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved);
  mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved);
  mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved);
  
  TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
  TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
  TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
  TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
  
  Int_t nEvents = mcTree->GetEntriesFast();
    for(Int_t i = 0; i < nEvents; i++)
    {
      mcTree->GetEvent(i);
      recTree->GetEvent(i);
      hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved);
      hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved);
      hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved);
      hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved);
    }
    
  c1->Divide(2,2);
  c1->cd(1);
  TString title = "Charged particles not removed";
  TString xtitle = "Cluster Multiplicity";
  TString ytitle = "N_{ch}";
  hChargedVsClusterMult->SetTitle(title);
  hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hChargedVsClusterMult->SetStats(0);
  hChargedVsClusterMult->Draw();
  TProfile *chProf = hChargedVsClusterMult->ProfileX();
  chProf->SetStats(0);
  chProf->Draw("SAME");
  //TF1 fitcharged("fitcharged","([0] + [1]*x)*(0.48/([2] + [3]*[2]))", 0, 100);
  TF1 fitcharged("fitcharged","pol1", 0, 100);
//   fitcharged.FixParameter(2, p0);
//   fitcharged.FixParameter(3, p1);
  TFitResultPtr chRes = chProf->Fit(&fitcharged,"S");
  TArrayD ch;
  if(!chRes)
  {
    ch = TArrayD(chRes->NPar(),chRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  c1->cd(2);
  title = "Neutral particles not removed";
  ytitle = "N_{neutral}";
  hNeutralVsClusterMult->SetTitle(title);
  hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hNeutralVsClusterMult->SetStats(0);
  hNeutralVsClusterMult->Draw();
  TProfile *neuProf = hNeutralVsClusterMult->ProfileX();
  neuProf->SetStats(0);
  neuProf->Draw("SAME");
  TF1 fitneutral("fitneutral","pol1", 0, 100);
  TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S");
  TArrayD neu;
  if(!neuRes)
  {
    neu = TArrayD(neuRes->NPar(),neuRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  c1->cd(3);
  
  title = "Gammas removed";
  ytitle = "N_{#gamma}";
  hGammaVsClusterMult->SetTitle(title);
  hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hGammaVsClusterMult->SetStats(0);
  hGammaVsClusterMult->Draw();
  TProfile *gamProf = hGammaVsClusterMult->ProfileX();
  gamProf->SetStats(0);
  gamProf->Draw("SAME");
  TF1 fitgamma("fitgamma","pol1", 0, 100);
  TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S");
  TArrayD gamma;
  if(!gammaRes)
  {
    gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  
  c1->cd(4);
  
  title = "Secondaries not removed";
  ytitle = "N_{sec}";
  hSecVsClusterMult->SetTitle(title);
  hSecVsClusterMult->GetXaxis()->SetTitle(xtitle);
  hSecVsClusterMult->GetYaxis()->SetTitle(ytitle);
  hSecVsClusterMult->SetStats(0);
  hSecVsClusterMult->Draw();
  TProfile *secProf = hSecVsClusterMult->ProfileX();
  secProf->SetStats(0);
  secProf->Draw("SAME");
  TF1 fitsecondary("fitsecondary","pol1", 0, 100);
  TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S");
  TArrayD sec;
  if(!secRes)
  {
    sec = TArrayD(secRes->NPar(),secRes->GetParams());
  }
  else
  {
    std::cout << "Could not extract charged contribution params" << std::endl;
  }
  Double_t meanCharged = 0.48/(p0 + p1*0.48);
  Double_t meanNeutral = 0.53/(p0 + p1*0.53);
  Double_t meanGamma = 0.51/(p0 + p1*0.51);
  Double_t meanSecondary = meanGamma; 
  
  AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections("TmCorrectionsPhos", fitcharged, fitneutral, fitgamma, fitsecondary, 
										   meanCharged, meanNeutral, meanGamma, meanSecondary );
  
  TF1 *func = generateRecEffFunction(p0, p1);
  AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection("ReCorrectionsPhos", *func, 1000);
  TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
  cor->Write();
  recor->Write();
  return 0;
}


TF1* generateRecEffFunction(Double_t p0, Double_t p1)
{
  TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200);
  
  Double_t params[2] = {p0, p1};
  f->SetParameters(params);
  return f;
}
Commit	Line	Data
943dc129	1	#include "TTree.h"
	2	#include "TFile.h"
	3	#include <TList.h>
	4	#include <Rtypes.h>
	5	#include "TCanvas.h"
	6	#include <TH2I.h>
	7	#include <TH2F.h>
	8	#include <TProfile.h>
	9	#include <TFitResultPtr.h>
	10	#include <TFitResult.h>
	11	#include "TF1.h"
	12	#include <iostream>
	13	#include <AliAnalysisEtTrackMatchCorrections.h>
	14	#include <AliAnalysisEtRecEffCorrection.h>
	15
	16	TCanvas *c1 = 0;
	17
	18	TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
02c62614	19	int createDummy(Double_t p0 = 0.366, Double_t p1 = 0.0)
	20	{
	21	TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
	22	TF1 fitneutral("fitneutral","0", 0, 100);
	23	TF1 fitcharged("fitcharged","0", 0, 100);
	24	TF1 fitgamma("fitgamma","0", 0, 100);
	25	TF1 fitsecondary("fitsecondary","0", 0, 100);
	26	AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections("TmCorrectionsPhos", fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
	27	TF1 *func = generateRecEffFunction(p0, p1);
	28	AliAnalysisEtRecEffCorrection recor = new AliAnalysisEtRecEffCorrection("ReCorrectionsPhos", func, 1000);
	29
	30	cor->Write();
	31	recor->Write();
	32	outfile->Close();
	33
943dc129	34
02c62614	35	}
943dc129	36	int calculateCorrections(TString filename="Et.ESD.simPbPb.PHOS.root", Double_t p0 = 0.366, Double_t p1 = 0.0)
	37	{
	38	c1 = new TCanvas;
	39
	40	TFile *f = TFile::Open(filename, "READ");
	41
	42	TList l = (TList)f->Get("out1");
	43
	44	TTree primaryTree = (TTree)l->FindObject("fPrimaryTreePhosMC");
	45	TTree recTree = (TTree)l->FindObject("fEventSummaryTreePhosRec");
	46	TTree mcTree = (TTree)l->FindObject("fEventSummaryTreePhosMC");
	47	std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl;
	48
	49	Int_t clusterMult = 0;
	50	Int_t nChargedNotRemoved = 0;
	51	Int_t nNeutralNotRemoved = 0;
	52	Int_t nGammaRemoved = 0;
	53	Int_t nSecNotRemoved = 0;
	54
	55	Double_t emEtRec = 0;
	56	Double_t emEtMc = 0;
	57
	58	recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult);
	59	mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved);
	60	mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved);
	61	mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved);
	62	mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved);
	63
	64	TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
	65	TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
	66	TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
	67	TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", 100, -0.5, 99.5, 100, -0.5, 99.5);
	68
	69	Int_t nEvents = mcTree->GetEntriesFast();
	70	for(Int_t i = 0; i < nEvents; i++)
	71	{
	72	mcTree->GetEvent(i);
	73	recTree->GetEvent(i);
	74	hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved);
	75	hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved);
	76	hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved);
	77	hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved);
	78	}
	79
	80	c1->Divide(2,2);
	81	c1->cd(1);
	82	TString title = "Charged particles not removed";
	83	TString xtitle = "Cluster Multiplicity";
	84	TString ytitle = "N_{ch}";
	85	hChargedVsClusterMult->SetTitle(title);
	86	hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle);
	87	hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle);
	88	hChargedVsClusterMult->SetStats(0);
	89	hChargedVsClusterMult->Draw();
	90	TProfile *chProf = hChargedVsClusterMult->ProfileX();
	91	chProf->SetStats(0);
	92	chProf->Draw("SAME");
	93	//TF1 fitcharged("fitcharged","([0] + [1]x)(0.48/([2] + [3]*[2]))", 0, 100);
	94	TF1 fitcharged("fitcharged","pol1", 0, 100);
	95	// fitcharged.FixParameter(2, p0);
	96	// fitcharged.FixParameter(3, p1);
	97	TFitResultPtr chRes = chProf->Fit(&fitcharged,"S");
	98	TArrayD ch;
	99	if(!chRes)
100	{
101	ch = TArrayD(chRes->NPar(),chRes->GetParams());
102	}
103	else
104	{
105	std::cout << "Could not extract charged contribution params" << std::endl;
106	}
107	c1->cd(2);
108	title = "Neutral particles not removed";
109	ytitle = "N_{neutral}";
110	hNeutralVsClusterMult->SetTitle(title);
111	hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle);
112	hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle);
113	hNeutralVsClusterMult->SetStats(0);
114	hNeutralVsClusterMult->Draw();
115	TProfile *neuProf = hNeutralVsClusterMult->ProfileX();
116	neuProf->SetStats(0);
117	neuProf->Draw("SAME");
118	TF1 fitneutral("fitneutral","pol1", 0, 100);
119	TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S");
120	TArrayD neu;
121	if(!neuRes)
122	{
123	neu = TArrayD(neuRes->NPar(),neuRes->GetParams());
124	}
125	else
126	{
127	std::cout << "Could not extract charged contribution params" << std::endl;
128	}
129	c1->cd(3);
130
131	title = "Gammas removed";
132	ytitle = "N_{#gamma}";
133	hGammaVsClusterMult->SetTitle(title);
134	hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle);
135	hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle);
136	hGammaVsClusterMult->SetStats(0);
137	hGammaVsClusterMult->Draw();
138	TProfile *gamProf = hGammaVsClusterMult->ProfileX();
139	gamProf->SetStats(0);
140	gamProf->Draw("SAME");
141	TF1 fitgamma("fitgamma","pol1", 0, 100);
142	TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S");
143	TArrayD gamma;
144	if(!gammaRes)
145	{
146	gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams());
147	}
148	else
149	{
150	std::cout << "Could not extract charged contribution params" << std::endl;
151	}
152
153	c1->cd(4);
154
155	title = "Secondaries not removed";
156	ytitle = "N_{sec}";
157	hSecVsClusterMult->SetTitle(title);
158	hSecVsClusterMult->GetXaxis()->SetTitle(xtitle);
159	hSecVsClusterMult->GetYaxis()->SetTitle(ytitle);
160	hSecVsClusterMult->SetStats(0);
161	hSecVsClusterMult->Draw();
162	TProfile *secProf = hSecVsClusterMult->ProfileX();
163	secProf->SetStats(0);
164	secProf->Draw("SAME");
165	TF1 fitsecondary("fitsecondary","pol1", 0, 100);
166	TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S");
167	TArrayD sec;
168	if(!secRes)
169	{
170	sec = TArrayD(secRes->NPar(),secRes->GetParams());
171	}
172	else
173	{
174	std::cout << "Could not extract charged contribution params" << std::endl;
175	}
176	Double_t meanCharged = 0.48/(p0 + p1*0.48);
177	Double_t meanNeutral = 0.53/(p0 + p1*0.53);
178	Double_t meanGamma = 0.51/(p0 + p1*0.51);
179	Double_t meanSecondary = meanGamma;
180
181	AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections("TmCorrectionsPhos", fitcharged, fitneutral, fitgamma, fitsecondary,
182	meanCharged, meanNeutral, meanGamma, meanSecondary );
183
184	TF1 *func = generateRecEffFunction(p0, p1);
185	AliAnalysisEtRecEffCorrection recor = new AliAnalysisEtRecEffCorrection("ReCorrectionsPhos", func, 1000);
186	TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
187	cor->Write();
188	recor->Write();
189	return 0;
190	}
191
192
193	TF1* generateRecEffFunction(Double_t p0, Double_t p1)
194	{
195	TF1 f = new TF1("receff", "x/([0] + x[1])", 0, 200);
196
197	Double_t params[2] = {p0, p1};
198	f->SetParameters(params);
199	return f;
200	}
201
202