[u/mrichter/AliRoot.git] / PWG0 / dNdEta / runMultiplicitySelector.C

/* $Id$ */

//
// script to run the AliMultiplicityESDSelector
//

#include "../CreateESDChain.C"
#include "../PWG0Helper.C"

void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
{
  if (aProof)
    connectProof("jgrosseo@lxb6046");

  TString libraries("libEG;libGeom;libESD;libPWG0base");
  TString packages("PWG0base");

  if (aMC != kFALSE)
  {
    libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6";
    packages += ";PWG0dep";
  }

  if (!prepareQuery(libraries, packages, kTRUE))
    return;

  gProof->SetParallel(0);
  gProof->SetLogLevel(4);
  gProof->SetParallel(9999);

  gROOT->ProcessLine(".L CreateCuts.C");
  gROOT->ProcessLine(".L drawPlots.C");

  // selection of esd tracks
  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
  if (!esdTrackCuts)
  {
    printf("ERROR: esdTrackCuts could not be created\n");
    return;
  }

  TList inputList;
  inputList.Add(esdTrackCuts);

  TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);

  TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
  AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);

  selectorName += ".cxx+";

  if (aDebug != kFALSE)
    selectorName += "g";

  Int_t result = executeQuery(chain, &inputList, selectorName, option);

  if (result != 0)
  {
    printf("ERROR: Executing process failed with %d.\n", result);
    return;
  }
}

void draw(const char* fileName = "multiplicityMC.root")
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");

  TFile::Open(fileName);
  mult->LoadHistograms("Multiplicity");

  mult->DrawHistograms();
}

void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");

  TFile::Open(fileName);
  mult->LoadHistograms("Multiplicity");

  mult->ApplyBayesianMethod(hist, kFALSE, eventType);
  mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());

  //mult->ApplyMinuitFit(hist, kFALSE);
  //mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());

  //mult->ApplyGaussianMethod(hist, kFALSE);

  return mult;
}

void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2)
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");

  TFile::Open(fileNameMC);
  mult->LoadHistograms("Multiplicity");

  TFile::Open(fileNameESD);
  TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));

  mult->SetMultiplicityESD(histID, hist);

  //mult->ApplyGaussianMethod(histID, kFALSE);
  //for (Float_t f=0; f<0.1; f+=0.01)
  //mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
  mult->ApplyMinuitFit(histID, kFALSE);
  mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, hist2->ProjectionY());

  //mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);

  return mult;
}

const char* GetRegName(Int_t type)
{
  switch (type)
  {
    case AliMultiplicityCorrection::kNone:      return "None"; break;
    case AliMultiplicityCorrection::kPol0:      return "Pol0"; break;
    case AliMultiplicityCorrection::kPol1:      return "Pol1"; break;
    case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
    case AliMultiplicityCorrection::kEntropy:   return "Reduced cross-entropy"; break;
  }
  return 0;
}

void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", Int_t histID = 2)
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");

  TFile::Open(fileNameMC);
  mult->LoadHistograms("Multiplicity");

  TFile::Open(fileNameESD);
  TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
  TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));

  mult->SetMultiplicityESD(histID, hist);

  TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
  TLegend* legend = new TLegend(0.6, 0.1, 0.98, 0.4);
  legend->SetFillColor(0);

  Float_t min = 1e10;
  Float_t max = 0;

  TGraph* first = 0;

  for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
  {
    TGraph* fitResultsMC = new TGraph;
    fitResultsMC->SetTitle(";Weight Parameter");
    TGraph* fitResultsRes = new TGraph;
    fitResultsRes->SetTitle(";Weight Parameter");

    fitResultsMC->SetFillColor(0);
    fitResultsRes->SetFillColor(0);
    fitResultsMC->SetMarkerStyle(19+type);
    fitResultsRes->SetMarkerStyle(19+type);
    fitResultsRes->SetMarkerColor(kRed);
    fitResultsRes->SetLineColor(kRed);

    legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
    legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));

    if (first == 0)
      first = fitResultsMC;

    for (Float_t weight = 1e-2; weight < 2e4; weight *= 1e2)
    {
      Float_t chi2MC = 0;
      Float_t residuals = 0;

      mult->SetRegularizationParameters(type, weight);
      mult->ApplyMinuitFit(histID, kFALSE);
      mult->DrawComparison(Form("MinuitChi2_%d_%f", type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
      mult->GetComparisonResults(chi2MC, residuals);

      fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
      fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);

      min = TMath::Min(min, TMath::Min(chi2MC, residuals));
      max = TMath::Max(max, TMath::Max(chi2MC, residuals));
    }

    fitResultsMC->Print();
    fitResultsRes->Print();

    canvas->cd();
    fitResultsMC->Draw(Form("%s CP", (type == AliMultiplicityCorrection::kPol0) ? "A" : "SAME"));
    fitResultsRes->Draw("SAME CP");
  }

  gPad->SetLogx();
  gPad->SetLogy();
  printf("min = %f, max = %f\n", min, max);
  if (min <= 0)
    min = 1e-5;
  first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);

  legend->Draw();

  canvas->SaveAs(Form("%s.gif", canvas->GetName()));
}

void Merge(Int_t n, const char** files, const char* output)
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
  TList list;
  for (Int_t i=0; i<n; ++i)
  {
    TString name("Multiplicity");
    if (i > 0)
      name.Form("Multiplicity%d", i);

    TFile::Open(files[i]);
    data[i] = new AliMultiplicityCorrection(name, name);
    data[i]->LoadHistograms("Multiplicity");
    if (i > 0)
      list.Add(data[i]);
  }

  data[0]->Merge(&list);

  data[0]->DrawHistograms();

  TFile::Open(output, "RECREATE");
  data[0]->SaveHistograms();
  gFile->Close();
}

void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
{
  gSystem->Load("libPWG0base");

  AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");

  TFile::Open(fileName);
  mult->LoadHistograms("Multiplicity");

  TF1* func = 0;

  if (caseNo >= 4)
  {
    func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
    func->SetParNames("scaling", "averagen", "k");
  }

  switch (caseNo)
  {
    case 0: func = new TF1("flat", "1"); break;
    case 1: func = new TF1("flat", "501-x"); break;
    case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
    case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
    case 4: func->SetParameters(1000, 10, 2); break;
    case 5: func->SetParameters(1000, 20, 3); break;
    case 6: func->SetParameters(1000, 30, 4); break;
    case 7: func->SetParameters(1000, 40, 5); break;

    default: return;
  }

  mult->SetGenMeasFromFunc(func, 2);

  mult->ApplyBayesianMethod(2, kFALSE);
  mult->ApplyMinuitFit(2, kFALSE);
  //mult->ApplyGaussianMethod(2, kFALSE);
}
Commit	Line	Data
0ab29cfa	1	/* $Id$ */
	2
	3	//
	4	// script to run the AliMultiplicityESDSelector
	5	//
	6
	7	#include "../CreateESDChain.C"
0bd1f8a0	8	#include "../PWG0Helper.C"
0ab29cfa	9
03244034	10	void runMultiplicitySelector(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aMC = kFALSE, Bool_t aDebug = kFALSE, Bool_t aProof = kFALSE, const char* option = "")
0ab29cfa	11	{
0bd1f8a0	12	if (aProof)
cfc19dd5	13	connectProof("jgrosseo@lxb6046");
0bd1f8a0	14
	15	TString libraries("libEG;libGeom;libESD;libPWG0base");
	16	TString packages("PWG0base");
10ebe68d	17
0ab29cfa	18	if (aMC != kFALSE)
0bd1f8a0	19	{
	20	libraries += ";libVMC;libMinuit;libSTEER;libPWG0dep;libEVGEN;libFASTSIM;libmicrocern;libpdf;libpythia6;libEGPythia6;libAliPythia6";
	21	packages += ";PWG0dep";
	22	}
	23
	24	if (!prepareQuery(libraries, packages, kTRUE))
	25	return;
0ab29cfa	26
cfc19dd5	27	gProof->SetParallel(0);
	28	gProof->SetLogLevel(4);
	29	gProof->SetParallel(9999);
	30
0ab29cfa	31	gROOT->ProcessLine(".L CreateCuts.C");
	32	gROOT->ProcessLine(".L drawPlots.C");
	33
0ab29cfa	34	// selection of esd tracks
	35	AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
	36	if (!esdTrackCuts)
	37	{
	38	printf("ERROR: esdTrackCuts could not be created\n");
	39	return;
	40	}
	41
0bd1f8a0	42	TList inputList;
	43	inputList.Add(esdTrackCuts);
	44
cfc19dd5	45	TChain* chain = CreateESDChain(data, nRuns, offset, kFALSE, kFALSE);
0ab29cfa	46
	47	TString selectorName = ((aMC == kFALSE) ? "AliMultiplicityESDSelector" : "AliMultiplicityMCSelector");
	48	AliLog::SetClassDebugLevel(selectorName, AliLog::kInfo);
	49
0bd1f8a0	50	selectorName += ".cxx+";
0ab29cfa	51
	52	if (aDebug != kFALSE)
	53	selectorName += "g";
	54
03244034	55	Int_t result = executeQuery(chain, &inputList, selectorName, option);
0ab29cfa	56
	57	if (result != 0)
	58	{
	59	printf("ERROR: Executing process failed with %d.\n", result);
	60	return;
	61	}
0a173978	62	}
0ab29cfa	63
0a173978	64	void draw(const char* fileName = "multiplicityMC.root")
	65	{
	66	gSystem->Load("libPWG0base");
	67
	68	AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
	69
	70	TFile::Open(fileName);
	71	mult->LoadHistograms("Multiplicity");
	72
	73	mult->DrawHistograms();
0ab29cfa	74	}
0ab29cfa	75
cfc19dd5	76	void* fit(const char* fileName = "multiplicityMC.root", Int_t hist = 2, Int_t eventType = 0)
0a173978	77	{
	78	gSystem->Load("libPWG0base");
	79
	80	AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
	81
	82	TFile::Open(fileName);
	83	mult->LoadHistograms("Multiplicity");
	84
cfc19dd5	85	mult->ApplyBayesianMethod(hist, kFALSE, eventType);
	86	mult->DrawComparison("Bayesian", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, eventType)->ProjectionY());
	87
	88	//mult->ApplyMinuitFit(hist, kFALSE);
	89	//mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, mult->GetMultiplicityMC(hist, AliMultiplicityCorrection::kTrVtx)->ProjectionY());
	90
	91	//mult->ApplyGaussianMethod(hist, kFALSE);
9ca6be09	92
	93	return mult;
	94	}
	95
cfc19dd5	96	void* fitOther(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityESD.root", Int_t histID = 2)
9ca6be09	97	{
	98	gSystem->Load("libPWG0base");
	99
	100	AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
	101
	102	TFile::Open(fileNameMC);
	103	mult->LoadHistograms("Multiplicity");
	104
	105	TFile::Open(fileNameESD);
cfc19dd5	106	TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
cfc19dd5	107	TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
9ca6be09	108
cfc19dd5	109	mult->SetMultiplicityESD(histID, hist);
9ca6be09	110
cfc19dd5	111	//mult->ApplyGaussianMethod(histID, kFALSE);
	112	//for (Float_t f=0; f<0.1; f+=0.01)
	113	//mult->ApplyBayesianMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
	114	mult->ApplyMinuitFit(histID, kFALSE);
	115	mult->DrawComparison("MinuitChi2", hist, kFALSE, kTRUE, hist2->ProjectionY());
	116
	117	//mult->ApplyLaszloMethod(histID, kFALSE, AliMultiplicityCorrection::kTrVtx);
0a173978	118
	119	return mult;
	120	}
cfc19dd5	121
	122	const char* GetRegName(Int_t type)
	123	{
	124	switch (type)
	125	{
	126	case AliMultiplicityCorrection::kNone: return "None"; break;
	127	case AliMultiplicityCorrection::kPol0: return "Pol0"; break;
	128	case AliMultiplicityCorrection::kPol1: return "Pol1"; break;
	129	case AliMultiplicityCorrection::kCurvature: return "TotalCurvature"; break;
	130	case AliMultiplicityCorrection::kEntropy: return "Reduced cross-entropy"; break;
	131	}
	132	return 0;
	133	}
	134
	135	void EvaluateChi2Method(const char* fileNameMC = "multiplicityMC.root", const char* fileNameESD = "multiplicityMC.root", Int_t histID = 2)
	136	{
	137	gSystem->Load("libPWG0base");
	138
	139	AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
	140
	141	TFile::Open(fileNameMC);
	142	mult->LoadHistograms("Multiplicity");
	143
	144	TFile::Open(fileNameESD);
	145	TH2F* hist = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityESD%d", histID));
	146	TH2F* hist2 = (TH2F*) gFile->Get(Form("Multiplicity/fMultiplicityVtx%d", histID));
	147
	148	mult->SetMultiplicityESD(histID, hist);
	149
	150	TCanvas* canvas = new TCanvas("EvaluateChi2Method", "EvaluateChi2Method", 800, 600);
	151	TLegend* legend = new TLegend(0.6, 0.1, 0.98, 0.4);
	152	legend->SetFillColor(0);
	153
	154	Float_t min = 1e10;
	155	Float_t max = 0;
	156
	157	TGraph* first = 0;
	158
	159	for (AliMultiplicityCorrection::RegularizationType type = AliMultiplicityCorrection::kPol0; type <= AliMultiplicityCorrection::kEntropy; ++type)
	160	{
	161	TGraph* fitResultsMC = new TGraph;
	162	fitResultsMC->SetTitle(";Weight Parameter");
	163	TGraph* fitResultsRes = new TGraph;
	164	fitResultsRes->SetTitle(";Weight Parameter");
	165
	166	fitResultsMC->SetFillColor(0);
	167	fitResultsRes->SetFillColor(0);
	168	fitResultsMC->SetMarkerStyle(19+type);
	169	fitResultsRes->SetMarkerStyle(19+type);
	170	fitResultsRes->SetMarkerColor(kRed);
	171	fitResultsRes->SetLineColor(kRed);
	172
	173	legend->AddEntry(fitResultsMC, Form("%s MC chi2", GetRegName(type)));
	174	legend->AddEntry(fitResultsRes, Form("%s residual chi2", GetRegName(type)));
	175
	176	if (first == 0)
	177	first = fitResultsMC;
	178
	179	for (Float_t weight = 1e-2; weight < 2e4; weight *= 1e2)
	180	{
	181	Float_t chi2MC = 0;
	182	Float_t residuals = 0;
	183
	184	mult->SetRegularizationParameters(type, weight);
185	mult->ApplyMinuitFit(histID, kFALSE);
186	mult->DrawComparison(Form("MinuitChi2_%d_%f", type, weight), histID, kFALSE, kTRUE, hist2->ProjectionY());
187	mult->GetComparisonResults(chi2MC, residuals);
188
189	fitResultsMC->SetPoint(fitResultsMC->GetN(), weight, chi2MC);
190	fitResultsRes->SetPoint(fitResultsRes->GetN(), weight, residuals);
191
192	min = TMath::Min(min, TMath::Min(chi2MC, residuals));
193	max = TMath::Max(max, TMath::Max(chi2MC, residuals));
194	}
195
196	fitResultsMC->Print();
197	fitResultsRes->Print();
198
199	canvas->cd();
200	fitResultsMC->Draw(Form("%s CP", (type == AliMultiplicityCorrection::kPol0) ? "A" : "SAME"));
201	fitResultsRes->Draw("SAME CP");
202	}
203
204	gPad->SetLogx();
205	gPad->SetLogy();
206	printf("min = %f, max = %f\n", min, max);
207	if (min <= 0)
208	min = 1e-5;
209	first->GetYaxis()->SetRangeUser(min * 0.5, max * 1.5);
210
211	legend->Draw();
212
213	canvas->SaveAs(Form("%s.gif", canvas->GetName()));
214	}
215
216	void Merge(Int_t n, const char** files, const char* output)
217	{
218	gSystem->Load("libPWG0base");
219
220	AliMultiplicityCorrection** data = new AliMultiplicityCorrection*[n];
221	TList list;
222	for (Int_t i=0; i<n; ++i)
223	{
224	TString name("Multiplicity");
225	if (i > 0)
226	name.Form("Multiplicity%d", i);
227
228	TFile::Open(files[i]);
229	data[i] = new AliMultiplicityCorrection(name, name);
230	data[i]->LoadHistograms("Multiplicity");
231	if (i > 0)
232	list.Add(data[i]);
233	}
234
235	data[0]->Merge(&list);
236
237	data[0]->DrawHistograms();
238
239	TFile::Open(output, "RECREATE");
240	data[0]->SaveHistograms();
241	gFile->Close();
242	}
243
244	void testMethod(Int_t caseNo, const char* fileName = "multiplicityMC.root")
245	{
246	gSystem->Load("libPWG0base");
247
248	AliMultiplicityCorrection* mult = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
249
250	TFile::Open(fileName);
251	mult->LoadHistograms("Multiplicity");
252
253	TF1* func = 0;
254
255	if (caseNo >= 4)
256	{
257	func = new TF1("nbd", "[0] * TMath::Binomial([2]+TMath::Nint(x)-1, [2]-1) * pow([1] / ([1]+[2]), TMath::Nint(x)) * pow(1 + [1]/[2], -[2])", 0, 50);
258	func->SetParNames("scaling", "averagen", "k");
259	}
260
261	switch (caseNo)
262	{
263	case 0: func = new TF1("flat", "1"); break;
264	case 1: func = new TF1("flat", "501-x"); break;
265	case 2: func = new TF1("flat", "1000 * 1/(x+1)"); break;
266	case 3: func = new TF1("flat", "1000 * TMath::Landau(x, 10, 5)"); break;
267	case 4: func->SetParameters(1000, 10, 2); break;
268	case 5: func->SetParameters(1000, 20, 3); break;
269	case 6: func->SetParameters(1000, 30, 4); break;
270	case 7: func->SetParameters(1000, 40, 5); break;
271
272	default: return;
273	}
274
275	mult->SetGenMeasFromFunc(func, 2);
276
277	mult->ApplyBayesianMethod(2, kFALSE);
278	mult->ApplyMinuitFit(2, kFALSE);
279	//mult->ApplyGaussianMethod(2, kFALSE);
280	}