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 Bool_t skipMean = 0;
  24
  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) ;
  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);
  27
  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) {
  29
  30   // Generic Macro to fit any particle using the PWG2/SPECTRA/Fit macros
  31   //
  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
  39   //
  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
  44   //   divided).
  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)
  48
  49   // Author: Michele Floris, CERN
  50
  51   // load stuff
  52   gSystem->Load("libTree.so");
  53   gSystem->Load("libVMC.so");
  54   gSystem->Load("libMinuit.so");
  55   gSystem->Load("libSTEERBase.so");
  56   gSystem->Load("libESD.so");
  57   gSystem->Load("libAOD.so");
  58   gSystem->Load("libANALYSIS.so");
  59   gSystem->Load("libANALYSISalice.so");
  60   gSystem->Load("libCORRFW.so");
  61   gSystem->Load("libPWG2spectra.so");
  62
  63
  64   // open file
  65   TFile * f = file ? new TFile(file) : gFile;
  66   TH1 * h = 0;
  67   if(listname){
  68     TList * l = (TList*) gDirectory->Get(listname);
  69     h = (TH1*) l->FindObject(histo);
  70   }
  71   else{
  72     h = (TH1*) gDirectory->Get(histo); //
  73   }
  74
  75   FitParticle(h, partName, min, max, scaleHisto, fitFunc, vartype, fileOut, wait);
  76
  77
  78
  79 }
  80
  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) {
  82
  83
  84   // get histo and draw
  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;
  89
  90   if (!TDatabasePDG::Instance()->GetParticle(partName)) {
  91     cout << "Wrong particle name " << partName << endl;
  92
  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;
  99     }
 100
 101     cout << "Try [ TDatabasePDG::Instance()->GetParticle(partName) ] by hand "<< endl;
 102     return;
 103
 104   }
 105   Double_t mass = TDatabasePDG::Instance()->GetParticle(partName)->Mass();
 106   cout << "Fitting ["<<partName<<"], mass: " << mass << endl;
 107
 108
 109   AliLatexTable * table;
 110
 111   TF1* func = 0;
 112   Int_t normPar = -1;
 113   Int_t slopePar=0;
 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 ");
 119     normPar = 0;
 120     slopePar = 2;
 121   }
 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 ");
 126     normPar = 0;
 127     slopePar = 1;
 128   }
 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  ");
 133     normPar = 0;
 134     slopePar = 1;
 135   }
 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  ");
 140     normPar = 0;
 141     slopePar = 1;
 142
 143   }
 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);
 151     normPar = 4;
 152     slopePar = 2;
 153   }
 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 ");
 158     normPar = 0;
 159     slopePar = 1;
 160   }
 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  ");
 165     normPar = 0;
 166     slopePar = 1;
 167   }
 168
 169     h->Sumw2();
 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
 174
 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);
 180
 181    //   gMinuit->Command("SET STRATEGY 2"); // FIXME
 182
 183   // if (!AliPWGHistoTools::Fit(h,func,min,max)) {
 184   //   cout << "Fitting error!" << endl;
 185   //   //    return;
 186   // } FIXME
 187   cout << "Drawing" << endl;
 188
 189   TCanvas * c = new TCanvas();
 190   cout << "1" << endl;
 191   c->SetLogy();
 192   cout << "2" << endl;
 193   h->Draw();
 194   cout << "3" << endl;
 195   func->Draw("same");
 196   cout << "4" << endl;
 197
 198   // Print results nicely
 199   // populate table
 200
 201   Double_t yield=0,yieldE=0;
 202   cout << "Y" << endl;
 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()));
 206   l->Draw();
 207
 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;
 211
 212   Double_t tslope   = func->GetParameter(slopePar);
 213   Double_t tslopeE = func->GetParError(slopePar);
 214
 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);
 222   }
 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);
 232
 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 , 0.,100.);
 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);
 243   table->InsertRow();
 244   table->PrintTable("ASCII");
 245
 246   if(fileOut) {
 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));
 249     c->Write();
 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
 257     fout->Close();
 258   }
 259
 260
 261   if(wait)  c->WaitPrimitive();
 262
 263
 264 }
 265
 266 // Double_t MicheleFunction(Double_t *x, Double_t *p) {
 267
 268 //   Double_t xloc = x[0];
 269 //   Double_t y = 0;
 270 //   if(xloc < p[0]) {
 271 //     y = xloc*TMath::Exp(-xloc/p[2]);
 272 //   } else if (xloc > p[0] && xloc < p[1]) {
 273
 274 //   }
 275
 276
 277 // }