1 #if !defined(__CINT__) || defined(__MAKECINT__)
4 #include "AliPWGHistoTools.h"
5 #include "AliPWGFunc.h"
6 #include "AliLatexTable.h"
11 #include "TDatabasePDG.h"
15 #include "THashList.h"
22 enum {kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave, kFitBoseEinstein, kFitFermiDirac};
25 void FitParticle(TH1 * h, const char * partName, 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, Float_t meanMin = 0., Float_t meanMax = 100) ;
26 void 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);
28 void FitParticle(const char * file, const char * histo, const char * partName, const char * listname, Float_t min, Float_t max, Float_t scaleHisto, Int_t fitFunc, Int_t vartype, const char * fileOut, Bool_t wait) {
30 // Generic Macro to fit any particle using the PWG2/SPECTRA/Fit macros
32 // You have to provide:
33 // - file: file name. If not set, uses gFile
34 // - histo: histogram name name (assumend to be in the main folder of the file)
35 // - listname: if different from 0, histo is looked in this tlist
36 // - partName: it is used to get the mass of the particle from
37 // TDatabasePDG. If you are not sure, just use a part of the name: a
38 // list of names matching it is printed on screen
40 // You can optionally provide:
41 // - min, max: the fit range
42 // - scaleHisto: a scaling factor for the histo. If negative, it is
43 // ignored. If the histo is scaled, the bin width is also
45 // - fitFunc: id of the function, levi is the default
46 // valide options: kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave
47 // - varType: the variable used in the pt spectrum (see AliPWGFunc.h)
49 // Author: Michele Floris, CERN
52 gSystem->Load("libTree");
53 gSystem->Load("libVMC");
54 gSystem->Load("libMinuit");
55 gSystem->Load("libSTEERBase");
56 gSystem->Load("libESD");
57 gSystem->Load("libAOD");
58 gSystem->Load("libANALYSIS");
59 gSystem->Load("libANALYSISalice");
60 gSystem->Load("libCORRFW");
61 gSystem->Load("libPWG2spectra");
65 TFile * f = file ? new TFile(file) : gFile;
68 TList * l = (TList*) gDirectory->Get(listname);
69 h = (TH1*) l->FindObject(histo);
72 h = (TH1*) gDirectory->Get(histo); //
75 FitParticle(h, partName, min, max, scaleHisto, fitFunc, vartype, fileOut, wait);
81 void FitParticle(TH1 * h, const char * partName, Float_t min , Float_t max, Float_t scaleHisto, Int_t fitFunc, Int_t vartype, const char * fileOut, Bool_t wait, Float_t meanMin, Float_t meanMax) {
85 AliPWGFunc * fm = new AliPWGFunc;
86 fm->SetVarType(AliPWGFunc::VarType_t(vartype));
87 // fm->SetVarType(AliPWGFunc::VarType_t(0));//FIXME
88 // cout << "Warning: hacked vartype" << endl;
90 if (!TDatabasePDG::Instance()->GetParticle(partName)) {
91 cout << "Wrong particle name " << partName << endl;
93 const THashList * l = TDatabasePDG::Instance()->ParticleList();
94 Int_t npart = l->GetSize();
95 for(Int_t ipart = 0; ipart < npart; ipart++){
96 TString name = l->At(ipart)->GetName();
97 if(name.Contains(partName, TString::kIgnoreCase))
98 cout << " - Did you mean [" << name.Data() << "] ?" << endl;
101 cout << "Try [ TDatabasePDG::Instance()->GetParticle(partName) ] by hand "<< endl;
105 Double_t mass = TDatabasePDG::Instance()->GetParticle(partName)->Mass();
106 cout << "Fitting ["<<partName<<"], mass: " << mass << endl;
109 AliLatexTable * table;
114 if (fitFunc == kFitLevi) {
115 func = fm->GetLevi (mass, 0.4, 750,3);
116 func->SetParLimits(1,0.0001,20000);
117 table = new AliLatexTable(10,"c|ccccccccc");
118 table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
122 if (fitFunc == kFitExpPt) {
123 func = fm->GetPTExp(0.2, 20);
124 table = new AliLatexTable(10,"c|ccccccccc");
125 table->InsertCustomRow("Part & Func & Yield & pT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
129 if (fitFunc == fFitExpMt) {
130 func = fm->GetMTExp(mass,0.2, 20);
131 table = new AliLatexTable(10,"c|ccccccccc");
132 table->InsertCustomRow("Part & Func & Yield & mT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
136 if (fitFunc == kFitBoltzmann) {
137 func = fm->GetBoltzmann(mass, 0.2, 20);
138 table = new AliLatexTable(10,"c|ccccccccc");
139 table->InsertCustomRow("Part & Func & Yield & Bz Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
144 if (fitFunc == kFitBlastWave) {
145 func = fm->GetBGBW(mass,0.6,0.3, 1, 1e5);// beta, T, n, norm
146 table = new AliLatexTable(12,"cc|cccccccccccc");
147 table->InsertCustomRow("Part & Func & Yield & T & beta & n & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
148 func->SetParLimits(1, 0.1, 0.99);
149 func->SetParLimits(2, 0.01, 1);
150 func->SetParLimits(3, 0.01, 2);
154 if (fitFunc == kFitBoseEinstein) {
155 func = fm->GetBoseEinstein(mass,0.3,20);
156 table = new AliLatexTable(10,"c|ccccccccc");
157 table->InsertCustomRow("Part & Func & Yield & BE Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
161 if (fitFunc == kFitFermiDirac) {
162 func = fm->GetFermiDirac(mass,0.3,20);
163 table = new AliLatexTable(10,"c|ccccccccc");
164 table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below & Mean ");
170 if (scaleHisto > 0) h->Scale(scaleHisto, "width");
171 // TH1 * h = AliPWGHistoTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get(histo)); // FIXME
172 // cout << "WARNING SCALING2PI" << endl;
173 // h->Scale(2*TMath::Pi());//Fixme
175 h->Fit(func); // FIXME
176 h->Fit(func); // FIXME
177 h->Fit(func); // FIXME
178 h->Fit(func,"IME0","",min,max);
179 h->Fit(func,"IME0","",min,max);
181 // gMinuit->Command("SET STRATEGY 2"); // FIXME
183 // if (!AliPWGHistoTools::Fit(h,func,min,max)) {
184 // cout << "Fitting error!" << endl;
187 cout << "Drawing" << endl;
189 TCanvas * c = new TCanvas();
198 // Print results nicely
201 Double_t yield=0,yieldE=0;
203 AliPWGHistoTools::GetYield(h,func,yield,yieldE);
204 cout << "YE" << endl;
205 TLatex * l = new TLatex(2,(h->GetMaximum()+h->GetMinimum())/2,Form("%3.3f #pm %3.3f (%s)", yield,yieldE, func->GetName()));
208 // Float_t yield = func->Integral(0.45,1.05);
209 // Float_t yieldE = func->IntegralError(0.45,1.05);
210 cout << "Slope Par: "<< slopePar << endl;
212 Double_t tslope = func->GetParameter(slopePar);
213 Double_t tslopeE = func->GetParError(slopePar);
215 table->SetNextCol(Form("%s (%s)", partName, h->GetName()));
216 table->SetNextCol(func->GetName());
217 table->SetNextCol(yield,yieldE,-4);
218 table->SetNextCol(tslope,tslopeE,-4);
219 if(fitFunc == kFitBlastWave) {
220 table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
221 table->SetNextCol(func->GetParameter(3),func->GetParError(3),-4);
223 // table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
224 table->SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
225 Float_t lowestPoint = TMath::Max(AliPWGHistoTools::GetLowestNotEmptyBinEdge(h),min);
226 // Float_t yieldAbove = func->Integral(lowestPoint,100);
227 Float_t yieldBelow = func->Integral(0,lowestPoint);
228 table->SetNextCol(lowestPoint,-3);
229 table->SetNextCol(min,-2);
230 table->SetNextCol(max,-2);
231 table->SetNextCol(yieldBelow/yield,-3);
233 Double_t mean=0, meane=0;
234 // Float_t mean2=0, mean2e=0;
235 // AliPWGHistoTools::GetMean (func, mean, meane , 0.,100., normPar);
236 AliPWGHistoTools::GetMeanDataAndExtrapolation (h, func, mean, meane , meanMin, meanMax);
237 // AliPWGHistoTools::GetMeanDataAndExtrapolation (h, func, mean, meane , 0.,4.5);
238 // AliPWGHistoTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
239 if(skipMean) table->SetNextCol("N/A");
240 else table->SetNextCol(mean, meane ,-4);
241 // table->SetNextCol(mean2, mean2e,-4);
242 // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
244 table->PrintTable("ASCII");
247 TFile * fout = new TFile(fileOut, "update");
248 c->SetName(Form("c_%s_%s_fit_%2.2f_%2.2f", h->GetName(), func->GetName(), min, max));
250 TH1F * hFit = new TH1F ("hFitResults", "hFitResults", 2,0,1);
251 hFit->SetBinContent(1, yield); hFit->SetBinError(1, yieldE);
252 hFit->SetBinContent(2, mean ); hFit->SetBinError(2, meane);
253 hFit->GetXaxis()->SetBinLabel(1, "yield");
254 hFit->GetXaxis()->SetBinLabel(2, "<pt>");
255 hFit->Write(Form("hFit_%s_%s_fit_%2.2f_%2.2f", h->GetName(), func->GetName(), min, max));
256 fout->Purge(); // remove old canvases
261 if(wait) c->WaitPrimitive();
266 // Double_t MicheleFunction(Double_t *x, Double_t *p) {
268 // Double_t xloc = x[0];
271 // y = xloc*TMath::Exp(-xloc/p[2]);
272 // } else if (xloc > p[0] && xloc < p[1]) {