]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG1/AliTreeDraw.cxx
Preliminary files for CMake
[u/mrichter/AliRoot.git] / PWG1 / AliTreeDraw.cxx
index 9592e767fadc2316d755e1ecf49871188bf5d81a..1ab865cd201a5733a55a25f8f5f30e07e4f854a8 100644 (file)
@@ -49,6 +49,8 @@ marian.ivanov@cern.ch
 #include "TGeometry.h"
 #include "TPolyLine3D.h"
 #include "TPolyMarker3D.h"
+#include "TObjString.h"
+
 
 //ALIROOT includes
 #include "AliTrackPointArray.h"
@@ -349,6 +351,109 @@ TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFit
   return hRes;
 }
 
+TH1F* AliTreeDraw::CreateResHistoI(TH2F* hRes2, TH1F **phMean, Int_t integ,  Bool_t drawBinFits)
+{
+  TVirtualPad* currentPad = gPad;
+  TAxis* axis = hRes2->GetXaxis();
+  Int_t nBins = axis->GetNbins();
+  Bool_t overflowBinFits = kFALSE;
+  TH1F* hRes, *hMean;
+  if (axis->GetXbins()->GetSize()){
+    hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
+    hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
+  }
+  else{
+    hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
+    hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
+
+  }
+  hRes->SetStats(false);
+  hRes->SetOption("E");
+  hRes->SetMinimum(0.);
+  //
+  hMean->SetStats(false);
+  hMean->SetOption("E");
+  // create the fit function
+  TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
+  
+  fitFunc->SetLineWidth(2);
+  fitFunc->SetFillStyle(0);
+  // create canvas for fits
+  TCanvas* canBinFits = NULL;
+  Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
+  Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
+  Int_t ny = (nPads-1) / nx + 1;
+  if (drawBinFits) {
+    canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
+    if (canBinFits) delete canBinFits;
+    canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
+    canBinFits->Divide(nx, ny);
+  }
+
+  // loop over x bins and fit projection
+  Int_t dBin = ((overflowBinFits) ? 1 : 0);
+  for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
+    if (drawBinFits) canBinFits->cd(bin + dBin);
+    Int_t bin0=TMath::Max(bin-integ,0);
+    Int_t bin1=TMath::Min(bin+integ,nBins);
+    TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
+    //    
+    if (hBin->GetEntries() > 5) {
+      fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
+      hBin->Fit(fitFunc,"s");
+      Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
+
+      if (sigma > 0.){
+       hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
+       hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
+      }
+      else{
+       hRes->SetBinContent(bin, 0.);
+       hMean->SetBinContent(bin,0);
+      }
+      hRes->SetBinError(bin, fitFunc->GetParError(2));
+      hMean->SetBinError(bin, fitFunc->GetParError(1));
+      
+      //
+      //
+
+    } else {
+      hRes->SetBinContent(bin, 0.);
+      hRes->SetBinError(bin, 0.);
+      hMean->SetBinContent(bin, 0.);
+      hMean->SetBinError(bin, 0.);
+    }
+    
+
+    if (drawBinFits) {
+      char name[256];
+      if (bin == 0) {
+       sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
+      } else if (bin == nBins+1) {
+       sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
+      } else {
+       sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
+               axis->GetTitle(), axis->GetBinUpEdge(bin));
+      }
+      canBinFits->cd(bin + dBin);
+      hBin->SetTitle(name);
+      hBin->SetStats(kTRUE);
+      hBin->DrawCopy("E");
+      canBinFits->Update();
+      canBinFits->Modified();
+      canBinFits->Update();
+    }
+    
+    delete hBin;
+  }
+
+  delete fitFunc;
+  currentPad->cd();
+  *phMean = hMean;
+  return hRes;
+}
+
 
 
 
@@ -385,3 +490,68 @@ void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
    }
 }
 
+
+
+
+TString* AliTreeDraw::FitPlane(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix, Int_t start, Int_t stop){
+   //
+   // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
+   // returns chi2, fitParam and covMatrix
+   // returns TString with fitted formula
+   //
+    
+   TString formulaStr(formula); 
+   TString drawStr(drawCommand);
+   TString cutStr(cuts);
+      
+   formulaStr.ReplaceAll("++", "~");
+   TObjArray* formulaTokens = formulaStr.Tokenize("~"); 
+   Int_t dim = formulaTokens->GetEntriesFast();
+   
+   fitParam.ResizeTo(dim);
+   covMatrix.ResizeTo(dim,dim);
+   
+   TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
+   fitter->StoreData(kTRUE);   
+   fitter->ClearPoints();
+   
+   Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
+   if (entries == -1) return new TString("An ERROR has occured during fitting!");
+   Double_t **values = new Double_t*[dim+1] ; 
+   
+   for (Int_t i = 0; i < dim + 1; i++){
+      Int_t centries = 0;
+      if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
+      else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
+      
+      if (entries != centries) return new TString("An ERROR has occured during fitting!");
+      values[i] = new Double_t[entries];
+      memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t)); 
+   }
+   
+   // add points to the fitter
+   for (Int_t i = 0; i < entries; i++){
+      Double_t x[1000];
+      for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
+      fitter->AddPoint(x, values[dim][i], 1);
+   }
+
+   fitter->Eval();
+   fitter->GetParameters(fitParam);
+   fitter->GetCovarianceMatrix(covMatrix);
+   chi2 = fitter->GetChisquare();
+   chi2 = chi2;
+   
+   TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
+   
+   for (Int_t iparam = 0; iparam < dim; iparam++) {
+     returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
+     if (iparam < dim-1) returnFormula.Append("+");
+   }
+   returnFormula.Append(" )");
+   delete formulaTokens;
+   delete fitter;
+   delete[] values;
+   return preturnFormula;
+}
+