STEER/STEERBase/AliTRDNDFast.cxx

   1 // Author: Daniel.Lohner@cern.ch
   2
   3 #include "AliTRDNDFast.h"
   4 #include "TCanvas.h"
   5
   6 extern Double_t langaufunc(Double_t *x, Double_t *par) {
   7
   8     // needs protection, parameter [3]>0!!!!!
   9     // needs protection, parameter [4]>0!!!!!
  10
  11     //Fit parameters:
  12     //par[0]=Width (scale) parameter of Landau density
  13     //par[1]=Most Probable (MP, location) parameter of Landau density
  14     //par[2]=Total area (integral -inf to inf, normalization constant)
  15     //par[3]=Width (sigma) of convoluted Gaussian function
  16     //par[4]=Exponential Slope Parameter
  17     //
  18     //In the Landau distribution (represented by the CERNLIB approximation),
  19     //the maximum is located at x=-0.22278298 with the location parameter=0.
  20     //This shift is corrected within this function, so that the actual
  21     //maximum is identical to the MP parameter.
  22
  23     // Numeric constants
  24     Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  25     Double_t mpshift  = -0.22278298;       // Landau maximum location
  26
  27     // Control constants
  28     Double_t np = 100.0;      // number of convolution stpdf
  29     Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
  30
  31     // Variables
  32     Double_t xx;
  33     Double_t mpc;
  34     Double_t fland;
  35     Double_t sum = 0.0;
  36     Double_t xlow,xupp;
  37     Double_t step;
  38     Double_t i;
  39
  40     // MP shift correction
  41     mpc = par[1] - mpshift * par[0];
  42
  43     // Range of convolution integral
  44     xlow = x[0] - sc * par[3];
  45     xupp = x[0] + sc * par[3];
  46
  47     if(xlow<0)xlow=0;
  48     if(xupp<xlow)cout<<"ERROR: Convolution around Negative MPV"<<endl;
  49
  50     step = (xupp-xlow) / np;
  51
  52     // Convolution integral of Landau and Gaussian by sum
  53     for(i=1.0; i<=np/2; i++) {
  54         xx = xlow + (i-.5) * step;
  55         fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];   // mpc taken out
  56         sum += fland * TMath::Gaus(x[0],xx,par[3]);
  57
  58         xx = xupp - (i-.5) * step;
  59         fland = TMath::Landau(xx,mpc,par[0])*TMath::Exp(-par[4]*xx) / par[0];   // mpc taken out
  60         sum += fland * TMath::Gaus(x[0],xx,par[3]);
  61     }
  62
  63     return (par[2] * step * sum * invsq2pi / par[3]);
  64 }
  65
  66
  67
  68 ClassImp(AliTRDNDFast);
  69
  70 AliTRDNDFast::AliTRDNDFast(): TObject(),
  71     fNDim(0),
  72     fTitle(""),
  73     fFunc(NULL),
  74     fHistos(NULL),
  75     fPars()
  76 {
  77     // default constructor
  78 }
  79
  80 AliTRDNDFast::AliTRDNDFast(const char *name,Int_t ndim,Int_t nbins,Double_t xlow,Double_t xup): TObject(),
  81 fNDim(ndim),
  82 fTitle(name),
  83 fFunc(NULL),
  84 fHistos(NULL),
  85 fPars()
  86 {
  87     Init();
  88     for(Int_t idim=0;idim<fNDim;idim++){
  89         fHistos[idim]=new TH1F(Form("%s_axis_%d",fTitle.Data(),idim),Form("%s_axis_%d",fTitle.Data(),idim),nbins,xlow,xup);
  90         fHistos[idim]->SetDirectory(0);
  91     }
  92 }
  93
  94 AliTRDNDFast::AliTRDNDFast(const AliTRDNDFast &ref) : TObject(ref),
  95 fNDim(ref.fNDim),
  96 fTitle(ref.fTitle),
  97 fFunc(NULL),
  98 fHistos(NULL),
  99 fPars()
 100 {
 101     Cleanup();
 102     Init();
 103     for(Int_t idim=0;idim<fNDim;idim++){
 104         fHistos[idim]=(TH1F*)ref.fHistos[idim]->Clone(Form("%s_axis_%d",GetName(),idim));
 105         fHistos[idim]->SetDirectory(0);
 106         for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
 107     }
 108 }
 109
 110 AliTRDNDFast &AliTRDNDFast::operator=(const AliTRDNDFast &ref){
 111     //
 112     // Assignment operator
 113     //
 114     if(this != &ref){
 115         TObject::operator=(ref);
 116         fNDim=ref.fNDim;
 117         fTitle=ref.fTitle;
 118         fFunc = ref.fFunc;
 119         fHistos=ref.fHistos;
 120         for(Int_t idim=0;idim<fNDim;idim++){
 121             for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[idim][ipar]=ref.fPars[idim][ipar];
 122         }
 123     }
 124     return *this;
 125 }
 126
 127 AliTRDNDFast::~AliTRDNDFast(){
 128     Cleanup();
 129
 130 }
 131
 132 void AliTRDNDFast::Init(){
 133
 134     for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].Set(fNDim);
 135     fFunc=new TF1*[fNDim];
 136     fHistos=new TH1F*[fNDim];
 137     for(Int_t idim=0;idim<fNDim;idim++){
 138         fHistos[idim]=NULL;
 139         for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(0,idim);
 140         fFunc[idim]=NULL;
 141     }
 142 }
 143
 144 void AliTRDNDFast::Cleanup(){
 145     if(fHistos){
 146         for(Int_t idim=0;idim<fNDim;idim++){
 147             if(fHistos[idim]){
 148                 delete fHistos[idim];
 149                 fHistos[idim]=NULL;
 150             }
 151         }
 152         delete[] fHistos;
 153         fHistos=NULL;
 154     }
 155     if(fPars){
 156         for(Int_t ipar=0;ipar<kNpar;ipar++){
 157             fPars[ipar].Reset();
 158         }
 159     }
 160     if(fFunc){
 161         for(Int_t idim=0;idim<fNDim;idim++){
 162             if(fFunc[idim]){
 163                 delete fFunc[idim];
 164                 fFunc[idim]=NULL;
 165             }
 166         }
 167         delete[] fFunc;
 168         fFunc=NULL;
 169     }
 170 }
 171
 172 void AliTRDNDFast::PrintPars(){
 173     for(Int_t idim=0;idim<fNDim;idim++){
 174         for(Int_t ipar=0;ipar<kNpar;ipar++)cout<<idim<<" "<<ipar<<" "<<GetParam(idim,ipar)<<endl;
 175     }
 176 }
 177
 178 void AliTRDNDFast::ScaleLangauFun(TF1 *func,Double_t mpv){
 179
 180     Double_t start[kNpar],low[kNpar],high[kNpar];
 181     for(Int_t ii=0;ii<kNpar;ii++){
 182         start[ii]=func->GetParameter(ii);
 183         func->GetParLimits(ii,low[ii],high[ii]);
 184     }
 185
 186     Double_t scalefactor=mpv/start[1];
 187
 188     for(Int_t ii=0;ii<kNpar;ii++){
 189         Double_t scale=1;
 190         if(ii==0||ii==1||ii==3)scale=scalefactor;
 191         if(ii==4)scale=1./scalefactor;
 192         start[ii]*=scale;
 193         low[ii]*=scale;
 194         high[ii]*=scale;
 195         func->SetParLimits(ii, low[ii], high[ii]);
 196     }
 197     func->SetParameters(start);
 198 }
 199
 200 //---------------------------------------------------------------
 201 //---------------------------------------------------------------
 202 TF1 *AliTRDNDFast::GetLangauFun(TString funcname,Double_t range[2],Double_t scalefactor){
 203
 204     Double_t start[kNpar] = {120, 1000, 1, 100, 1e-5};
 205     Double_t low[kNpar] = {10, 300, 0.01, 1, 0.};
 206     Double_t high[kNpar] = {1000, 3000, 100, 1000, 1.};
 207
 208     TF1 *hlandauE=new TF1(funcname.Data(),langaufunc,0,range[1],kNpar);
 209     hlandauE->SetParameters(start);
 210     hlandauE->SetParNames("Width","MP","Area","GSigma","delta");
 211     for (int i=0; i<kNpar; i++) {
 212         hlandauE->SetParLimits(i, low[i], high[i]);
 213     }
 214
 215     if(scalefactor!=1){ScaleLangauFun(hlandauE,scalefactor*start[1]);}
 216
 217     return hlandauE;
 218 }
 219
 220 TF1 *AliTRDNDFast::FitLandau(TString name,TH1F *htemp,Double_t range[2],TString option){
 221
 222     TF1 *fitlandau1D=GetLangauFun(name,range);
 223     TF1 fit("land","landau");
 224     Double_t max=htemp->GetXaxis()->GetBinCenter(htemp->GetMaximumBin());
 225     fit.SetParameter(1,max);
 226     fit.SetParLimits(1,0,htemp->GetXaxis()->GetXmax());
 227     fit.SetParameter(2,0.3*max); // MPV may never become negative!!!!!
 228     htemp->Fit("land",option.Data(),"",range[0],range[1]);
 229     ScaleLangauFun(fitlandau1D,fit.GetParameter(1));
 230     htemp->Fit(fitlandau1D,option.Data(),"",range[0],range[1]); // range for references
 231     return fitlandau1D;
 232 }
 233
 234 void AliTRDNDFast::BuildHistos(){
 235
 236     for(Int_t idim=0;idim<fNDim;idim++){
 237         fHistos[idim]->Reset();
 238         // Fill Histo
 239         Double_t pars[kNpar];
 240         for(Int_t ipar=0;ipar<kNpar;ipar++)pars[ipar]=GetParam(idim,ipar);
 241         // Also Fill overflow bins!!!
 242         for(Int_t ii=1;ii<=fHistos[idim]->GetNbinsX()+1;ii++){
 243             Double_t xval=fHistos[idim]->GetXaxis()->GetBinCenter(ii);
 244             Double_t val=langaufunc(&xval,&pars[0]);
 245             //Double_t val=fFunc[idim]->Eval(xval);
 246             fHistos[idim]->SetBinContent(ii,val);
 247         }
 248         // Normalization
 249         fHistos[idim]->Scale(1./fHistos[idim]->Integral(1,fHistos[idim]->GetNbinsX(),"width"));
 250     }
 251 }
 252
 253 void AliTRDNDFast::Build(Double_t **pars){
 254     // pars[ndim][npar]
 255     for(Int_t idim=0;idim<fNDim;idim++){
 256         for(Int_t ipar=0;ipar<kNpar;ipar++){
 257             fPars[ipar].SetAt(pars[idim][ipar],idim);
 258         }
 259     }
 260     BuildHistos();
 261 }
 262
 263
 264 void AliTRDNDFast::Build(TH1F **hdEdx,TString path){
 265
 266     Double_t range[2];
 267
 268     TCanvas *CANV=new TCanvas("a","a");
 269     if(fNDim==2)CANV->Divide(2,1);
 270     if(fNDim==3)CANV->Divide(2,2);
 271     if(fNDim>6)CANV->Divide(3,3);
 272     // Fit and Extract Parameters
 273     for(Int_t idim=0;idim<fNDim;idim++){
 274         CANV->cd(idim+1);
 275         gPad->SetLogy();
 276         range[0]=hdEdx[idim]->GetXaxis()->GetXmin();
 277         range[1]=hdEdx[idim]->GetXaxis()->GetXmax();
 278         // Norm Histogram
 279         hdEdx[idim]->Scale(1./hdEdx[idim]->Integral(1,hdEdx[idim]->GetNbinsX(),"width"));
 280         // Fit Histogram
 281         fFunc[idim]=FitLandau(Form("fit%d",idim),hdEdx[idim],range,"RMB0");
 282         // Norm Landau
 283         fFunc[idim]->SetParameter(2,fFunc[idim]->GetParameter(2)/fFunc[idim]->Integral(range[0],range[1]));
 284         hdEdx[idim]->DrawCopy();
 285         fFunc[idim]->DrawCopy("same");
 286         // Set Pars
 287         for(Int_t ipar=0;ipar<kNpar;ipar++)fPars[ipar].SetAt(fFunc[idim]->GetParameter(ipar),idim);
 288     }
 289     if(path!="")CANV->Print(Form("%s/%s_Build.pdf",path.Data(),fTitle.Data()));
 290     delete CANV;
 291
 292     BuildHistos();
 293 }
 294
 295 Double_t AliTRDNDFast::Eval(Double_t *point) const{
 296     Double_t val=1;
 297     for(Int_t idim=0;idim<fNDim;idim++){
 298         Int_t bin=fHistos[idim]->GetXaxis()->FindBin(point[idim]);
 299         val*=fHistos[idim]->GetBinContent(bin);
 300     }
 301     return val;
 302 }
 303
 304 void AliTRDNDFast::Random(Double_t *point) const{
 305     for(Int_t idim=0;idim<fNDim;idim++){
 306         point[idim]=fHistos[idim]->GetRandom();
 307     }
 308 }
 309
 310 void AliTRDNDFast::Random(Double_t *point,AliTRDNDFast *nd0,AliTRDNDFast *nd1,Double_t w0,Double_t w1){
 311     for(Int_t idim=0;idim<nd0->fNDim;idim++){
 312         point[idim]=GetRandomInterpolation(nd0->fHistos[idim],nd1->fHistos[idim],w0,w1);
 313     }
 314 }
 315
 316 Int_t AliTRDNDFast::BinarySearchInterpolation(Int_t start,Int_t end,Double_t *a0,Double_t *a1,Double_t w0,Double_t w1,Double_t val){
 317
 318     if((end-start)==1)return start;
 319     Int_t mid=Int_t((end+start)/2);
 320     Double_t valmid=(w0*a0[mid]+w1*a1[mid])/(w0+w1);
 321     if(val>=valmid)return BinarySearchInterpolation(mid,end,a0,a1,w0,w1,val);
 322     return BinarySearchInterpolation(start,mid,a0,a1,w0,w1,val);
 323 }
 324
 325 Double_t AliTRDNDFast::GetRandomInterpolation(TH1F *hist0,TH1F *hist1,Double_t w0,Double_t w1){
 326
 327     // Draw Random Variable
 328     Double_t rand=gRandom->Rndm();
 329
 330     // Get Integral Arrays
 331     Double_t *integral0=hist0->GetIntegral();
 332     Double_t *integral1=hist1->GetIntegral();
 333
 334     Int_t nbinsX=hist0->GetNbinsX();
 335
 336     Int_t ibin=BinarySearchInterpolation(1,nbinsX+1,integral0,integral1,w0,w1,rand);
 337     return hist0->GetXaxis()->GetBinCenter(ibin);
 338 }
 339
 340
 341