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