PWGLF/SPECTRA/PiKaPr/COMBINED/FitParticle.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <iostream>
   3
   4 #include "AliPWGHistoTools.h"
   5 #include "AliPWGFunc.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 = AliPWGFunc::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 AliPWGFunc.h)
  44
  45   // Author: Michele Floris, CERN
  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
  61   AliPWGFunc * fm = new AliPWGFunc;
  62   fm->SetVarType(AliPWGFunc::VarType_t(vartype));
  63   //  fm->SetVarType(AliPWGFunc::VarType_t(0));//FIXME
  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();
  82   cout << "Fitting ["<<partName<<"], mass: " << mass << endl;
  83
  84
  85   AliLatexTable * table;
  86
  87   TF1* func = 0;
  88   Int_t normPar = -1;
  89   Int_t slopePar=0;
  90   if (fitFunc == kFitLevi)   {
  91     func = fm->GetLevi (mass, 0.4, 20,3);
  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 ");
  95     normPar = 0;
  96     slopePar = 2;
  97   }
  98   if (fitFunc == kFitExpPt)  {
  99     func = fm->GetPTExp(0.2, 20);
 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;
 104   }
 105   if (fitFunc == fFitExpMt)  {
 106     func = fm->GetMTExp(mass,0.2, 20);
 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;
 111   }
 112   if (fitFunc == kFitBoltzmann)  {
 113     func = fm->GetBoltzmann(mass, 0.2, 20);
 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
 119   }
 120   if (fitFunc == kFitBlastWave)  {
 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;
 143   }
 144
 145
 146   TFile * f = file ? new TFile(file) : gFile;
 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   }
 155   h->Sumw2();
 156   if (scaleHisto > 0) h->Scale(scaleHisto, "width");
 157 //   TH1 * h = AliPWGHistoTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get(histo)); // FIXME
 158 //   cout << "WARNING SCALING2PI" << endl;
 159 //   h->Scale(2*TMath::Pi());//Fixme
 160
 161    h->Fit(func); // FIXME
 162    gMinuit->Command("SET STRATEGY 2"); // FIXME
 163
 164   if (!AliPWGHistoTools::Fit(h,func,min,max)) {
 165     cout << "Fitting error!" << endl;
 166     //    return;
 167   }
 168   cout << "Drawing" << endl;
 169
 170   TCanvas * c = new TCanvas();
 171   cout << "1" << endl;
 172   c->SetLogy();
 173   cout << "2" << endl;
 174   h->Draw();
 175   cout << "3" << endl;
 176   func->Draw("same");
 177   cout << "4" << endl;
 178
 179   // Print results nicely
 180   // populate table
 181
 182   Double_t yield=0,yieldE=0;
 183   cout << "Y" << endl;
 184   AliPWGHistoTools::GetYield(h,func,yield,yieldE);
 185   cout << "YE" << endl;
 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
 198 //   Float_t yield  = func->Integral(0.45,1.05);
 199 //   Float_t yieldE = func->IntegralError(0.45,1.05);
 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()));
 215   Float_t lowestPoint = TMath::Max(AliPWGHistoTools::GetLowestNotEmptyBinEdge(h),min);
 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;
 225   // AliPWGHistoTools::GetMean      (func, mean,  meane , 0.,100., normPar);
 226   // AliPWGHistoTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
 227   // table->SetNextCol(mean,  meane ,-4);
 228   // table->SetNextCol(mean2, mean2e,-4);
 229   //                     fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
 230   table->InsertRow();
 231   table->PrintTable("ASCII");
 232
 233
 234
 235 }