Migrating PWG2/SPECTRA/Fit to new PWG structure
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / COMBINED / FitParticle.C
CommitLineData
362a1b75 1#if !defined(__CINT__) || defined(__MAKECINT__)
2#include <iostream>
3
fef6f2b3 4#include "AliPWGHistoTools.h"
5#include "AliPWGFunc.h"
362a1b75 6#include "AliLatexTable.h"
7#include "TF1.h"
8#include "TFile.h"
9#include "TH1.h"
10#include "TROOT.h"
11#include "TDatabasePDG.h"
12#include "TCanvas.h"
13#include "TMath.h"
14#include "TSystem.h"
15#include "THashList.h"
16#include "TMinuit.h"
a4097895 17#include "TLatex.h"
362a1b75 18
19using namespace std;
20#endif
21
a4097895 22enum {kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave, kFitBoseEinstein, kFitFermiDirac};
362a1b75 23
fef6f2b3 24void FitParticle(const char * file, const char * histo, const char * partName, const char * listname=0, Float_t min = 0, Float_t max =3, Float_t scaleHisto = -1., Int_t fitFunc = kFitLevi, Int_t vartype = AliPWGFunc::kdNdpt, const char * fileOut = 0, Bool_t wait = 0) {
362a1b75 25
26 // Generic Macro to fit any particle using the PWG2/SPECTRA/Fit macros
27 //
28 // You have to provide:
a4097895 29 // - file: file name. If not set, uses gFile
0e019ec6 30 // - histo: histogram name name (assumend to be in the main folder of the file)
90421222 31 // - listname: if different from 0, histo is looked in this tlist
362a1b75 32 // - partName: it is used to get the mass of the particle from
33 // TDatabasePDG. If you are not sure, just use a part of the name: a
34 // list of names matching it is printed on screen
35 //
36 // You can optionally provide:
37 // - min, max: the fit range
38 // - scaleHisto: a scaling factor for the histo. If negative, it is
39 // ignored. If the histo is scaled, the bin width is also
40 // divided).
41 // - fitFunc: id of the function, levi is the default
5e1a7931 42 // valide options: kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave
fef6f2b3 43 // - varType: the variable used in the pt spectrum (see AliPWGFunc.h)
44
45 // Author: Michele Floris, CERN
362a1b75 46
47 // load stuff
48 gSystem->Load("libTree.so");
49 gSystem->Load("libVMC.so");
50 gSystem->Load("libMinuit.so");
51 gSystem->Load("libSTEERBase.so");
52 gSystem->Load("libESD.so");
53 gSystem->Load("libAOD.so");
54 gSystem->Load("libANALYSIS.so");
55 gSystem->Load("libANALYSISalice.so");
56 gSystem->Load("libCORRFW.so");
57 gSystem->Load("libPWG2spectra.so");
58
59
60 // get histo and draw
fef6f2b3 61 AliPWGFunc * fm = new AliPWGFunc;
62 fm->SetVarType(AliPWGFunc::VarType_t(vartype));
63 // fm->SetVarType(AliPWGFunc::VarType_t(0));//FIXME
362a1b75 64 // cout << "Warning: hacked vartype" << endl;
65
66 if (!TDatabasePDG::Instance()->GetParticle(partName)) {
67 cout << "Wrong particle name " << partName << endl;
68
69 const THashList * l = TDatabasePDG::Instance()->ParticleList();
70 Int_t npart = l->GetSize();
71 for(Int_t ipart = 0; ipart < npart; ipart++){
72 TString name = l->At(ipart)->GetName();
73 if(name.Contains(partName, TString::kIgnoreCase))
74 cout << " - Did you mean [" << name.Data() << "] ?" << endl;
75 }
76
77 cout << "Try [ TDatabasePDG::Instance()->GetParticle(partName) ] by hand "<< endl;
78 return;
79
80 }
81 Double_t mass = TDatabasePDG::Instance()->GetParticle(partName)->Mass();
0e019ec6 82 cout << "Fitting ["<<partName<<"], mass: " << mass << endl;
83
a4097895 84
85 AliLatexTable * table;
86
362a1b75 87 TF1* func = 0;
88 Int_t normPar = -1;
a4097895 89 Int_t slopePar=0;
362a1b75 90 if (fitFunc == kFitLevi) {
5e1a7931 91 func = fm->GetLevi (mass, 0.4, 20,3);
a4097895 92 func->SetParLimits(1,0.0001,20000);
93 table = new AliLatexTable(9,"c|cccccccc");
94 table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
362a1b75 95 normPar = 0;
a4097895 96 slopePar = 2;
362a1b75 97 }
98 if (fitFunc == kFitExpPt) {
99 func = fm->GetPTExp(0.2, 20);
a4097895 100 table = new AliLatexTable(9,"c|cccccccc");
101 table->InsertCustomRow("Part & Func & Yield & pT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
102 normPar = 0;
103 slopePar = 1;
362a1b75 104 }
5e1a7931 105 if (fitFunc == fFitExpMt) {
106 func = fm->GetMTExp(mass,0.2, 20);
a4097895 107 table = new AliLatexTable(9,"c|cccccccc");
108 table->InsertCustomRow("Part & Func & Yield & mT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
109 normPar = 0;
110 slopePar = 1;
5e1a7931 111 }
112 if (fitFunc == kFitBoltzmann) {
113 func = fm->GetBoltzmann(mass, 0.2, 20);
a4097895 114 table = new AliLatexTable(9,"c|cccccccc");
115 table->InsertCustomRow("Part & Func & Yield & Bz Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
116 normPar = 0;
117 slopePar = 1;
118
5e1a7931 119 }
120 if (fitFunc == kFitBlastWave) {
a4097895 121 func = fm->GetBGBW(mass,0.6,0.3, 1, 1e5);// beta, T, n, norm
122 table = new AliLatexTable(11,"cc|ccccccccccc");
123 table->InsertCustomRow("Part & Func & Yield & T & beta & n & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
124 func->SetParLimits(1, 0.1, 0.99);
125 func->SetParLimits(2, 0.01, 1);
126 func->SetParLimits(3, 0.01, 2);
127 normPar = 4;
128 slopePar = 2;
129 }
130 if (fitFunc == kFitBoseEinstein) {
131 func = fm->GetBoseEinstein(mass,0.3,20);
132 table = new AliLatexTable(9,"c|cccccccc");
133 table->InsertCustomRow("Part & Func & Yield & BE Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
134 normPar = 0;
135 slopePar = 1;
136 }
137 if (fitFunc == kFitFermiDirac) {
138 func = fm->GetFermiDirac(mass,0.3,20);
139 table = new AliLatexTable(9,"c|cccccccc");
140 table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
141 normPar = 0;
142 slopePar = 1;
5e1a7931 143 }
144
a4097895 145
146 TFile * f = file ? new TFile(file) : gFile;
90421222 147 TH1 * h = 0;
148 if(listname){
149 TList * l = (TList*) gDirectory->Get(listname);
150 h = (TH1*) l->FindObject(histo);
151 }
152 else{
153 h = (TH1*) gDirectory->Get(histo); //
154 }
362a1b75 155 h->Sumw2();
156 if (scaleHisto > 0) h->Scale(scaleHisto, "width");
fef6f2b3 157// TH1 * h = AliPWGHistoTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get(histo)); // FIXME
362a1b75 158// cout << "WARNING SCALING2PI" << endl;
159// h->Scale(2*TMath::Pi());//Fixme
160
5e1a7931 161 h->Fit(func); // FIXME
162 gMinuit->Command("SET STRATEGY 2"); // FIXME
362a1b75 163
fef6f2b3 164 if (!AliPWGHistoTools::Fit(h,func,min,max)) {
362a1b75 165 cout << "Fitting error!" << endl;
5e1a7931 166 // return;
362a1b75 167 }
5e1a7931 168 cout << "Drawing" << endl;
362a1b75 169
170 TCanvas * c = new TCanvas();
5e1a7931 171 cout << "1" << endl;
362a1b75 172 c->SetLogy();
5e1a7931 173 cout << "2" << endl;
362a1b75 174 h->Draw();
5e1a7931 175 cout << "3" << endl;
362a1b75 176 func->Draw("same");
5e1a7931 177 cout << "4" << endl;
362a1b75 178
362a1b75 179 // Print results nicely
362a1b75 180 // populate table
a4097895 181
90421222 182 Double_t yield=0,yieldE=0;
5e1a7931 183 cout << "Y" << endl;
fef6f2b3 184 AliPWGHistoTools::GetYield(h,func,yield,yieldE);
5e1a7931 185 cout << "YE" << endl;
a4097895 186 TLatex * l = new TLatex(2,(h->GetMaximum()+h->GetMinimum())/2,Form("%3.3f #pm %3.3f (%s)", yield,yieldE, func->GetName()));
187 l->Draw();
188
189 if(wait) c->WaitPrimitive();
190 if(fileOut) {
191 TFile * fout = new TFile(fileOut, "update");
192 c->SetName(Form("c_%s_%s_fit_%2.2f_%2.2f",histo, func->GetName(), min, max));
193 c->Write();
194 fout->Purge(); // remove old canvases
195 fout->Close();
196 }
197
362a1b75 198// Float_t yield = func->Integral(0.45,1.05);
199// Float_t yieldE = func->IntegralError(0.45,1.05);
a4097895 200 cout << "Slope Par: "<< slopePar << endl;
201
202 Double_t tslope = func->GetParameter(slopePar);
203 Double_t tslopeE = func->GetParError(slopePar);
204
205 table->SetNextCol(Form("%s (%s)", partName, h->GetName()));
206 table->SetNextCol(func->GetName());
207 table->SetNextCol(yield,yieldE,-4);
208 table->SetNextCol(tslope,tslopeE,-4);
209 if(fitFunc == kFitBlastWave) {
210 table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
211 table->SetNextCol(func->GetParameter(3),func->GetParError(3),-4);
212 }
213 // table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
214 table->SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
fef6f2b3 215 Float_t lowestPoint = TMath::Max(AliPWGHistoTools::GetLowestNotEmptyBinEdge(h),min);
a4097895 216 // Float_t yieldAbove = func->Integral(lowestPoint,100);
217 Float_t yieldBelow = func->Integral(0,lowestPoint);
218 table->SetNextCol(lowestPoint,-3);
219 table->SetNextCol(min,-2);
220 table->SetNextCol(max,-2);
221 table->SetNextCol(yieldBelow/yield,-3);
222
223 // Float_t mean=0, meane=0;
224 // Float_t mean2=0, mean2e=0;
fef6f2b3 225 // AliPWGHistoTools::GetMean (func, mean, meane , 0.,100., normPar);
226 // AliPWGHistoTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
a4097895 227 // table->SetNextCol(mean, meane ,-4);
228 // table->SetNextCol(mean2, mean2e,-4);
362a1b75 229 // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
a4097895 230 table->InsertRow();
231 table->PrintTable("ASCII");
362a1b75 232
233
234
235}