[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFECorrectSpectrumBase.cxx


/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
// Class for spectrum correction
// Subtraction of hadronic background, Unfolding of the data and
// Renormalization done here
// The following containers have to be set:
//  - Correction framework container for real data
//  - Correction framework container for MC (Efficiency Map)
//  - Correction framework container for background coming from data
//  - Correction framework container for background coming from MC
//
//  Author: 
//            Raphaelle Bailhache <R.Bailhache@gsi.de>
//            Markus Fasel <M.Fasel@gsi.de>
//

#include <TArrayD.h>
#include <TH1.h>
#include <TList.h>
#include <TObjArray.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TMath.h>
#include <TAxis.h>
#include <TGraphErrors.h>
#include <TFile.h>
#include <TPad.h>
#include <TH2D.h>
#include <TF1.h>

#include "AliPID.h"
#include "AliCFContainer.h"
#include "AliCFDataGrid.h"
#include "AliCFEffGrid.h"
#include "AliCFGridSparse.h"
#include "AliCFUnfolding.h"
#include "AliLog.h"

#include "AliHFECorrectSpectrumBase.h"
#include "AliHFEcuts.h"
#include "AliHFEcontainer.h"
#include "AliHFEtools.h"

ClassImp(AliHFECorrectSpectrumBase)

//____________________________________________________________
AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const char *name):
  TNamed(name, ""),
  fCFContainers(new TObjArray(kNbCFContainers)),
  fCorrelation(NULL),
  fEfficiencyFunction(NULL),
  fEtaSelected(kFALSE),
  fSetSmoothing(kFALSE),
  fNbDimensions(1),
  fNEvents(0),
  fStepMC(-1),
  fStepTrue(-1),
  fStepData(-1),
  fStepBeforeCutsV0(-1),
  fStepAfterCutsV0(-1),
  fStepGuessedUnfolding(-1),
  fNumberOfIterations(10),
  fChargeChoosen(kAllCharge),
  fTestCentralityLow(-1),
  fTestCentralityHigh(-1)
{
  //
  // Default constructor
  //

  memset(fEtaRange, 0, sizeof(Double_t) * 2);
  memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
 
}
//____________________________________________________________
AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref):
  TNamed(ref),
  fCFContainers(NULL),
  fCorrelation(NULL),
  fEfficiencyFunction(NULL),
  fEtaSelected(ref.fEtaSelected),
  fSetSmoothing(ref.fSetSmoothing),
  fNbDimensions(ref.fNbDimensions),
  fNEvents(ref.fNEvents),
  fStepMC(ref.fStepMC),
  fStepTrue(ref.fStepTrue),
  fStepData(ref.fStepData),
  fStepBeforeCutsV0(ref.fStepBeforeCutsV0),
  fStepAfterCutsV0(ref.fStepAfterCutsV0),
  fStepGuessedUnfolding(ref.fStepGuessedUnfolding),
  fNumberOfIterations(ref.fNumberOfIterations),
  fChargeChoosen(ref.fChargeChoosen),
  fTestCentralityLow(ref.fTestCentralityLow),
  fTestCentralityHigh(ref.fTestCentralityHigh)
{
  //
  // Copy constructor
  //
  ref.Copy(*this);

}
//____________________________________________________________
AliHFECorrectSpectrumBase &AliHFECorrectSpectrumBase::operator=(const AliHFECorrectSpectrumBase &ref){
  //
  // Assignment operator
  //
  if(this == &ref) 
    ref.Copy(*this);
  return *this;
}
//____________________________________________________________
void AliHFECorrectSpectrumBase::Copy(TObject &o) const {
  // 
  // Copy into object o
  //
  AliHFECorrectSpectrumBase &target = dynamic_cast<AliHFECorrectSpectrumBase &>(o);
  target.fCFContainers = fCFContainers;
  target.fCorrelation = fCorrelation;
  target.fEfficiencyFunction = fEfficiencyFunction;
  target.fEtaSelected = fEtaSelected;
  target.fSetSmoothing = fSetSmoothing;
  target.fNbDimensions = fNbDimensions;
  target.fNEvents = fNEvents;
  target.fStepMC = fStepMC;
  target.fStepTrue = fStepTrue;
  target.fStepData = fStepData;
  target.fStepBeforeCutsV0 = fStepBeforeCutsV0;
  target.fStepAfterCutsV0 = fStepAfterCutsV0;
  target.fStepGuessedUnfolding = fStepGuessedUnfolding;
  target.fNumberOfIterations = fNumberOfIterations;
  target.fChargeChoosen = fChargeChoosen;
  target.fTestCentralityLow = fTestCentralityLow;
  target.fTestCentralityHigh = fTestCentralityHigh;
  target.fEtaRange[0] = fEtaRange[0];
  target.fEtaRange[1] = fEtaRange[1];
  target.fEtaRangeNorm[0] = fEtaRangeNorm[0];
  target.fEtaRangeNorm[1] = fEtaRangeNorm[1];

}

//____________________________________________________________
AliHFECorrectSpectrumBase::~AliHFECorrectSpectrumBase(){
  //
  // Destructor
  //
  if(fCFContainers) delete fCFContainers;
 
}
//__________________________________________________________________________________
TGraphErrors *AliHFECorrectSpectrumBase::Normalize(THnSparse * const spectrum) const {
  //
  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
  // Give the final pt spectrum to be compared
  //

  if(fNEvents > 0) {

    TH1D* projection = spectrum->Projection(0);
    CorrectFromTheWidth(projection);
    TGraphErrors *graphError = NormalizeTH1(projection);
    return graphError;
  
  }
    
  return 0x0;
  

}
//__________________________________________________________________________________
TGraphErrors *AliHFECorrectSpectrumBase::Normalize(AliCFDataGrid * const spectrum) const {
  //
  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
  // Give the final pt spectrum to be compared
  //
  if(fNEvents > 0) {

    TH1D* projection = (TH1D *) spectrum->Project(0);
    CorrectFromTheWidth(projection);
    TGraphErrors *graphError = NormalizeTH1(projection);

    return graphError;
    
  }

  return 0x0;
  

}
//__________________________________________________________________________________
TGraphErrors *AliHFECorrectSpectrumBase::NormalizeTH1(TH1 *input) const {
  //
  // Normalize the spectrum to 1/(2*Pi*p_{T})*dN/dp_{T} (GeV/c)^{-2}
  // Give the final pt spectrum to be compared
  //
  Double_t chargecoefficient = 0.5;
  if(fChargeChoosen != 0) chargecoefficient = 1.0;

  Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
  printf("Normalizing Eta Range %f\n", etarange);
  printf("Number of events in Normalisation: %d\n", fNEvents);
  AliDebug(3, Form("charge coefficient: %f\n", chargecoefficient));
  if(fNEvents > 0) {

    TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
    Double_t p = 0, dp = 0; Int_t point = 1;
    Double_t n = 0, dN = 0;
    Double_t nCorr = 0, dNcorr = 0;
    //Double_t errdN = 0, errdp = 0;
    Double_t errdN = 0;
    for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
      point = ibin - input->GetXaxis()->GetFirst();
      p = input->GetXaxis()->GetBinCenter(ibin);
      //dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
      n = input->GetBinContent(ibin);
      AliDebug(6, Form("p: %f, n: %e\n", p, n));
      dN = input->GetBinError(ibin);
      // New point
      nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * 1./(2. * TMath::Pi() * p) * n;
      errdN = 1./(2. * TMath::Pi() * p);
      //errdp = 1./(2. * TMath::Pi() * p*p) * n;
      //dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
      dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN);
      
      spectrumNormalized->SetPoint(point, p, nCorr);
      spectrumNormalized->SetPointError(point, dp, dNcorr);
    }
    spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
    spectrumNormalized->SetMarkerStyle(22);
    spectrumNormalized->SetMarkerColor(kBlue);
    spectrumNormalized->SetLineColor(kBlue);

    return spectrumNormalized;
    
  }

  return 0x0;
  

}
//____________________________________________________________
void AliHFECorrectSpectrumBase::SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type){
  //
  // Set the container for a given step to the 
  //
  if(!fCFContainers) fCFContainers = new TObjArray(kNbCFContainers);
  fCFContainers->AddAt(cont, type);
}

//____________________________________________________________
AliCFContainer *AliHFECorrectSpectrumBase::GetContainer(AliHFECorrectSpectrumBase::CFContainer_t type){
  //
  // Get Correction Framework Container for given type
  //
  if(!fCFContainers) return NULL;
  return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
}
//____________________________________________________________
AliCFContainer *AliHFECorrectSpectrumBase::GetSlicedContainer(AliCFContainer *container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) {
  //
  // Slice bin for a given source of electron
  // nDim is the number of dimension the corrections are done
  // dimensions are the definition of the dimensions
  // source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
  // positivenegative if we want to keep positive (1) or negative (0) or both (-1)
  // centrality (-1 means we do not cut on centrality)
  //
  const double kVerySmall = 1e-5;
  
  Double_t *varMin = new Double_t[container->GetNVar()],
           *varMax = new Double_t[container->GetNVar()];

  Double_t *binLimits;
  for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
    
    binLimits = new Double_t[container->GetNBins(ivar)+1];
    container->GetBinLimits(ivar,binLimits);
    varMin[ivar] = binLimits[0];
    varMax[ivar] = binLimits[container->GetNBins(ivar)];
    // source
    if(ivar == 4){
      if((source>= 0) && (source<container->GetNBins(ivar))) {
	      varMin[ivar] = binLimits[source];
	      varMax[ivar] = binLimits[source];
      }     
    }
    // charge
    if(ivar == 3) {
      if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
    }
    // eta
    if(ivar == 1){
      for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++) 
        AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic))); 
      if(fEtaSelected){
        varMin[ivar] = fEtaRange[0];
        varMax[ivar] = fEtaRange[1];
      }
    }
    if(fEtaSelected){
      fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
      fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
      AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[1]));
    }
    // centrality
    if(ivar == 5){
	    if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar)) && doCentralityProjection) {
	      varMin[ivar] = binLimits[centralitylow]+ kVerySmall;
	      varMax[ivar] = binLimits[centralityhigh]+ kVerySmall;

	      TAxis *axistest = container->GetAxis(5,0);
	      AliDebug(1, Form("Number of bin in centrality direction %d\n",axistest->GetNbins()));
	      AliDebug(1, Form("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]));
	      Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
	      Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
	      AliDebug(1, Form("Low centrality %f and high centrality %f\n",lowcentrality,highcentrality));
	    }
    }
    
    // Protect varmax against overflow 
    if(TMath::Abs(varMax[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
      AliInfo("Protection against overflow bin");
      varMax[ivar] -= kVerySmall;
    }
    if(varMax[ivar] > binLimits[container->GetNBins(ivar)]){
      AliError("Upper limit exceeds allowed range");
      varMax[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
    }

    // Protect varmin against overflow 
    if(TMath::Abs(varMin[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
      AliInfo("Protection against overflow bin");
      varMin[ivar] -= kVerySmall;
    }
    if(varMin[ivar] > binLimits[container->GetNBins(ivar)]){
      AliError("Upper limit exceeds allowed range");
      varMin[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
    }
    AliDebug(1, Form("variable %d: Settting limits to %f and %f\n", ivar, varMin[ivar], varMax[ivar]));
    delete[] binLimits;
   
  }

  AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
  delete[] varMin; delete[] varMax;

  return k;

}

//_________________________________________________________________________
THnSparseF *AliHFECorrectSpectrumBase::GetSlicedCorrelation(THnSparseF *correlationmatrix, Int_t nDim, Int_t *dimensions,Chargetype_t charge,Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) const {
  //
  // Slice correlation
  //

  Int_t ndimensions = correlationmatrix->GetNdimensions();
  //printf("Number of dimension %d correlation map\n",ndimensions);
  if(ndimensions < (2*nDim)) {
    AliError("Problem in the dimensions");
    return NULL;
  }
  
  // Cut in centrality is centrality > -1
  if((5+((Int_t)(ndimensions/2.))) < ndimensions) {
    if((centralitylow >=0) && (centralityhigh >=0)) {
      
      TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
      TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
      
      Int_t bins0 = axiscentrality0->GetNbins();
      Int_t bins1 = axiscentrality1->GetNbins();
      
      AliDebug(1, Form("Number of centrality bins: %d and %d\n",bins0,bins1));
      if(bins0 != bins1) {
	      AliError("Problem in the dimensions");
	      return NULL;
      }
      
      if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0) && doCentralityProjection) {
	      axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
	      axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
	
	      Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
	      Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
	      Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
	      Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
	      AliDebug(1,Form("0 Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0));
	      AliDebug(1,Form("1 Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1));
	
      }    
    }
  }

  // Cut in eta > -1
  if(fEtaSelected){
    if((1+((Int_t)(ndimensions/2.))) < ndimensions) {
      
      TAxis *axiseta0 = correlationmatrix->GetAxis(1);
      TAxis *axiseta1 = correlationmatrix->GetAxis(1+((Int_t)(ndimensions/2.)));
      
      Int_t bins0 = axiseta0->GetNbins();
      Int_t bins1 = axiseta1->GetNbins();
      
      AliDebug(1, Form("Number of eta bins: %d and %d\n",bins0,bins1));
      if(bins0 != bins1) {
	AliError("Problem in the dimensions");
	return NULL;
      }
      
      axiseta0->SetRangeUser(fEtaRange[0],fEtaRange[1]);
      axiseta1->SetRangeUser(fEtaRange[0],fEtaRange[1]);
	
	Double_t loweta0 = axiseta0->GetBinLowEdge(axiseta0->FindBin(fEtaRange[0]));
	Double_t higheta0 = axiseta0->GetBinUpEdge(axiseta0->FindBin(fEtaRange[1]));
	Double_t loweta1 = axiseta1->GetBinLowEdge(axiseta1->FindBin(fEtaRange[0]));
	Double_t higheta1 = axiseta1->GetBinUpEdge(axiseta1->FindBin(fEtaRange[1]));
	AliInfo(Form("0 Low eta %f and high eta %f\n",loweta0,higheta0));
	AliInfo(Form("1 Low eta %f and high eta %f\n",loweta1,higheta1));
	
    }    
  }

  // Cut in charge
  if(charge != kAllCharge) {
    if((3+((Int_t)(ndimensions/2.))) < ndimensions) {
      
      TAxis *axischarge0 = correlationmatrix->GetAxis(3);
      TAxis *axischarge1 = correlationmatrix->GetAxis(3+((Int_t)(ndimensions/2.)));
      
      Int_t bins0 = axischarge0->GetNbins();
      Int_t bins1 = axischarge1->GetNbins();
      
      AliDebug(1, Form("Number of charge bins: %d and %d\n",bins0,bins1));
      if(bins0 != bins1) {
	AliError("Problem in the dimensions");
	return NULL;
      }
      
      axischarge0->SetRangeUser(charge,charge);
      axischarge1->SetRangeUser(charge,charge);
      
      Double_t lowcharge0 = axischarge0->GetBinLowEdge(axischarge0->FindBin(charge));
      Double_t highcharge0 = axischarge0->GetBinUpEdge(axischarge0->FindBin(charge));
      Double_t lowcharge1 = axischarge1->GetBinLowEdge(axischarge1->FindBin(charge));
      Double_t highcharge1 = axischarge1->GetBinUpEdge(axischarge1->FindBin(charge));
      AliInfo(Form("0 Low charge %f and high charge %f\n",lowcharge0,highcharge0));
      AliInfo(Form("1 Low charge %f and high charge %f\n",lowcharge1,highcharge1));
      
    }
  }
  

  Int_t ndimensionsContainer = (Int_t) ndimensions/2;
  
  Int_t *dim = new Int_t[nDim*2];
  for(Int_t iter=0; iter < nDim; iter++){
    dim[iter] = dimensions[iter];
    dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
  }
    
  THnSparseF *k = (THnSparseF *) correlationmatrix->Projection(nDim*2,dim);

  delete[] dim; 
  return k;
  
}
//___________________________________________________________________________
void AliHFECorrectSpectrumBase::CorrectFromTheWidth(TH1D *h1) const {
  //
  // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
  //

  TAxis *axis = h1->GetXaxis();
  Int_t nbinX = h1->GetNbinsX();

  for(Int_t i = 1; i <= nbinX; i++) {

    Double_t width = axis->GetBinWidth(i);
    Double_t content = h1->GetBinContent(i);
    Double_t error = h1->GetBinError(i); 
    h1->SetBinContent(i,content/width); 
    h1->SetBinError(i,error/width);
  }

}

//___________________________________________________________________________
void AliHFECorrectSpectrumBase::CorrectStatErr(AliCFDataGrid *backgroundGrid) const { 
  //
  // Correct statistical error
  //

  TH1D *h1 = (TH1D*)backgroundGrid->Project(0);
  Int_t nbinX = h1->GetNbinsX();
  Int_t bins[1];
  for(Long_t i = 1; i <= nbinX; i++) {
    bins[0] = i;
    Float_t content = h1->GetBinContent(i);
    if(content>0){
      Float_t error = TMath::Sqrt(content);
      backgroundGrid->SetElementError(bins, error);
    }
  }
}
//_________________________________________________________________________
TObject* AliHFECorrectSpectrumBase::GetSpectrum(const AliCFContainer * const c, Int_t step) {
  AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
  return data;
}
//_________________________________________________________________________
TObject* AliHFECorrectSpectrumBase::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
  // 
  // Create efficiency grid and calculate efficiency
  // of step to step0
  //
  TString name("eff");
  name += step;
  name+= step0;
  AliCFEffGrid* eff = new AliCFEffGrid((const char*)name,"",*c);
  eff->CalculateEfficiency(step,step0);
  return eff;
}
Commit	Line	Data
959ea9d8	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16	//
	17	// Class for spectrum correction
	18	// Subtraction of hadronic background, Unfolding of the data and
	19	// Renormalization done here
	20	// The following containers have to be set:
	21	// - Correction framework container for real data
	22	// - Correction framework container for MC (Efficiency Map)
	23	// - Correction framework container for background coming from data
	24	// - Correction framework container for background coming from MC
	25	//
	26	// Author:
	27	// Raphaelle Bailhache <R.Bailhache@gsi.de>
	28	// Markus Fasel <M.Fasel@gsi.de>
	29	//
	30
	31	#include <TArrayD.h>
	32	#include <TH1.h>
	33	#include <TList.h>
	34	#include <TObjArray.h>
	35	#include <TROOT.h>
	36	#include <TCanvas.h>
	37	#include <TLegend.h>
	38	#include <TStyle.h>
	39	#include <TMath.h>
	40	#include <TAxis.h>
	41	#include <TGraphErrors.h>
	42	#include <TFile.h>
	43	#include <TPad.h>
	44	#include <TH2D.h>
	45	#include <TF1.h>
	46
	47	#include "AliPID.h"
	48	#include "AliCFContainer.h"
	49	#include "AliCFDataGrid.h"
	50	#include "AliCFEffGrid.h"
	51	#include "AliCFGridSparse.h"
	52	#include "AliCFUnfolding.h"
	53	#include "AliLog.h"
	54
	55	#include "AliHFECorrectSpectrumBase.h"
	56	#include "AliHFEcuts.h"
	57	#include "AliHFEcontainer.h"
	58	#include "AliHFEtools.h"
	59
	60	ClassImp(AliHFECorrectSpectrumBase)
	61
	62	//____________________________________________________________
	63	AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const char *name):
	64	TNamed(name, ""),
63bdf450	65	fCFContainers(new TObjArray(kNbCFContainers)),
959ea9d8	66	fCorrelation(NULL),
	67	fEfficiencyFunction(NULL),
	68	fEtaSelected(kFALSE),
	69	fSetSmoothing(kFALSE),
	70	fNbDimensions(1),
	71	fNEvents(0),
	72	fStepMC(-1),
	73	fStepTrue(-1),
	74	fStepData(-1),
	75	fStepBeforeCutsV0(-1),
	76	fStepAfterCutsV0(-1),
	77	fStepGuessedUnfolding(-1),
	78	fNumberOfIterations(10),
	79	fChargeChoosen(kAllCharge),
	80	fTestCentralityLow(-1),
	81	fTestCentralityHigh(-1)
	82	{
	83	//
	84	// Default constructor
	85	//
	86
	87	memset(fEtaRange, 0, sizeof(Double_t) * 2);
	88	memset(fEtaRangeNorm, 0, sizeof(Double_t) * 2);
	89
	90	}
	91	//____________________________________________________________
	92	AliHFECorrectSpectrumBase::AliHFECorrectSpectrumBase(const AliHFECorrectSpectrumBase &ref):
	93	TNamed(ref),
	94	fCFContainers(NULL),
	95	fCorrelation(NULL),
	96	fEfficiencyFunction(NULL),
	97	fEtaSelected(ref.fEtaSelected),
	98	fSetSmoothing(ref.fSetSmoothing),
	99	fNbDimensions(ref.fNbDimensions),
	100	fNEvents(ref.fNEvents),
	101	fStepMC(ref.fStepMC),
	102	fStepTrue(ref.fStepTrue),
	103	fStepData(ref.fStepData),
	104	fStepBeforeCutsV0(ref.fStepBeforeCutsV0),
	105	fStepAfterCutsV0(ref.fStepAfterCutsV0),
	106	fStepGuessedUnfolding(ref.fStepGuessedUnfolding),
	107	fNumberOfIterations(ref.fNumberOfIterations),
	108	fChargeChoosen(ref.fChargeChoosen),
	109	fTestCentralityLow(ref.fTestCentralityLow),
	110	fTestCentralityHigh(ref.fTestCentralityHigh)
	111	{
	112	//
	113	// Copy constructor
	114	//
	115	ref.Copy(*this);
	116
	117	}
	118	//____________________________________________________________
	119	AliHFECorrectSpectrumBase &AliHFECorrectSpectrumBase::operator=(const AliHFECorrectSpectrumBase &ref){
	120	//
	121	// Assignment operator
	122	//
	123	if(this == &ref)
	124	ref.Copy(*this);
	125	return *this;
	126	}
	127	//____________________________________________________________
	128	void AliHFECorrectSpectrumBase::Copy(TObject &o) const {
	129	//
130	// Copy into object o
131	//
132	AliHFECorrectSpectrumBase &target = dynamic_cast<AliHFECorrectSpectrumBase &>(o);
133	target.fCFContainers = fCFContainers;
134	target.fCorrelation = fCorrelation;
135	target.fEfficiencyFunction = fEfficiencyFunction;
136	target.fEtaSelected = fEtaSelected;
137	target.fSetSmoothing = fSetSmoothing;
138	target.fNbDimensions = fNbDimensions;
139	target.fNEvents = fNEvents;
140	target.fStepMC = fStepMC;
141	target.fStepTrue = fStepTrue;
142	target.fStepData = fStepData;
143	target.fStepBeforeCutsV0 = fStepBeforeCutsV0;
144	target.fStepAfterCutsV0 = fStepAfterCutsV0;
145	target.fStepGuessedUnfolding = fStepGuessedUnfolding;
146	target.fNumberOfIterations = fNumberOfIterations;
147	target.fChargeChoosen = fChargeChoosen;
148	target.fTestCentralityLow = fTestCentralityLow;
149	target.fTestCentralityHigh = fTestCentralityHigh;
150	target.fEtaRange[0] = fEtaRange[0];
151	target.fEtaRange[1] = fEtaRange[1];
152	target.fEtaRangeNorm[0] = fEtaRangeNorm[0];
153	target.fEtaRangeNorm[1] = fEtaRangeNorm[1];
154
155	}
156
157	//____________________________________________________________
158	AliHFECorrectSpectrumBase::~AliHFECorrectSpectrumBase(){
159	//
160	// Destructor
161	//
162	if(fCFContainers) delete fCFContainers;
163
164	}
165	//__________________________________________________________________________________
166	TGraphErrors AliHFECorrectSpectrumBase::Normalize(THnSparse const spectrum) const {
167	//
168	// Normalize the spectrum to 1/(2Pip_{T})*dN/dp_{T} (GeV/c)^{-2}
169	// Give the final pt spectrum to be compared
170	//
171
172	if(fNEvents > 0) {
173
174	TH1D* projection = spectrum->Projection(0);
175	CorrectFromTheWidth(projection);
176	TGraphErrors *graphError = NormalizeTH1(projection);
177	return graphError;
178
179	}
180
181	return 0x0;
182
183
184	}
185	//__________________________________________________________________________________
186	TGraphErrors AliHFECorrectSpectrumBase::Normalize(AliCFDataGrid const spectrum) const {
187	//
188	// Normalize the spectrum to 1/(2Pip_{T})*dN/dp_{T} (GeV/c)^{-2}
189	// Give the final pt spectrum to be compared
190	//
191	if(fNEvents > 0) {
192
193	TH1D* projection = (TH1D *) spectrum->Project(0);
194	CorrectFromTheWidth(projection);
195	TGraphErrors *graphError = NormalizeTH1(projection);
196
197	return graphError;
198
199	}
200
201	return 0x0;
202
203
204	}
205	//__________________________________________________________________________________
206	TGraphErrors AliHFECorrectSpectrumBase::NormalizeTH1(TH1 input) const {
207	//
208	// Normalize the spectrum to 1/(2Pip_{T})*dN/dp_{T} (GeV/c)^{-2}
209	// Give the final pt spectrum to be compared
210	//
211	Double_t chargecoefficient = 0.5;
212	if(fChargeChoosen != 0) chargecoefficient = 1.0;
213
214	Double_t etarange = fEtaSelected ? fEtaRangeNorm[1] - fEtaRangeNorm[0] : 1.6;
215	printf("Normalizing Eta Range %f\n", etarange);
7483e78e	216	printf("Number of events in Normalisation: %d\n", fNEvents);
7483e78e	217	AliDebug(3, Form("charge coefficient: %f\n", chargecoefficient));
959ea9d8	218	if(fNEvents > 0) {
	219
	220	TGraphErrors *spectrumNormalized = new TGraphErrors(input->GetNbinsX());
	221	Double_t p = 0, dp = 0; Int_t point = 1;
	222	Double_t n = 0, dN = 0;
	223	Double_t nCorr = 0, dNcorr = 0;
7483e78e	224	//Double_t errdN = 0, errdp = 0;
7483e78e	225	Double_t errdN = 0;
959ea9d8	226	for(Int_t ibin = input->GetXaxis()->GetFirst(); ibin <= input->GetXaxis()->GetLast(); ibin++){
	227	point = ibin - input->GetXaxis()->GetFirst();
	228	p = input->GetXaxis()->GetBinCenter(ibin);
7483e78e	229	//dp = input->GetXaxis()->GetBinWidth(ibin)/2.;
959ea9d8	230	n = input->GetBinContent(ibin);
7483e78e	231	AliDebug(6, Form("p: %f, n: %e\n", p, n));
959ea9d8	232	dN = input->GetBinError(ibin);
	233	// New point
	234	nCorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * 1./(2. * TMath::Pi() * p) * n;
	235	errdN = 1./(2. * TMath::Pi() * p);
7483e78e	236	//errdp = 1./(2. * TMath::Pi() * pp) n;
	237	//dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN dN + errdp errdp * dp *dp);
	238	dNcorr = chargecoefficient * 1./etarange * 1./(Double_t)(fNEvents) * TMath::Sqrt(errdN * errdN * dN *dN);
959ea9d8	239
	240	spectrumNormalized->SetPoint(point, p, nCorr);
	241	spectrumNormalized->SetPointError(point, dp, dNcorr);
	242	}
	243	spectrumNormalized->GetXaxis()->SetTitle("p_{T} [GeV/c]");
	244	spectrumNormalized->GetYaxis()->SetTitle("#frac{1}{2 #pi p_{T}} #frac{dN}{dp_{T}} / [GeV/c]^{-2}");
	245	spectrumNormalized->SetMarkerStyle(22);
	246	spectrumNormalized->SetMarkerColor(kBlue);
	247	spectrumNormalized->SetLineColor(kBlue);
	248
	249	return spectrumNormalized;
	250
	251	}
	252
	253	return 0x0;
	254
	255
	256	}
	257	//____________________________________________________________
	258	void AliHFECorrectSpectrumBase::SetContainer(AliCFContainer *cont, AliHFECorrectSpectrumBase::CFContainer_t type){
	259	//
	260	// Set the container for a given step to the
	261	//
63bdf450	262	if(!fCFContainers) fCFContainers = new TObjArray(kNbCFContainers);
959ea9d8	263	fCFContainers->AddAt(cont, type);
	264	}
	265
	266	//____________________________________________________________
	267	AliCFContainer *AliHFECorrectSpectrumBase::GetContainer(AliHFECorrectSpectrumBase::CFContainer_t type){
	268	//
	269	// Get Correction Framework Container for given type
	270	//
	271	if(!fCFContainers) return NULL;
	272	return dynamic_cast<AliCFContainer *>(fCFContainers->At(type));
	273	}
	274	//____________________________________________________________
7483e78e	275	AliCFContainer AliHFECorrectSpectrumBase::GetSlicedContainer(AliCFContainer container, Int_t nDim, Int_t *dimensions,Int_t source,Chargetype_t charge, Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) {
959ea9d8	276	//
	277	// Slice bin for a given source of electron
	278	// nDim is the number of dimension the corrections are done
	279	// dimensions are the definition of the dimensions
	280	// source is if we want to keep only one MC source (-1 means we don't cut on the MC source)
	281	// positivenegative if we want to keep positive (1) or negative (0) or both (-1)
	282	// centrality (-1 means we do not cut on centrality)
	283	//
7483e78e	284	const double kVerySmall = 1e-5;
959ea9d8	285
	286	Double_t *varMin = new Double_t[container->GetNVar()],
	287	*varMax = new Double_t[container->GetNVar()];
	288
	289	Double_t *binLimits;
	290	for(Int_t ivar = 0; ivar < container->GetNVar(); ivar++){
	291
	292	binLimits = new Double_t[container->GetNBins(ivar)+1];
	293	container->GetBinLimits(ivar,binLimits);
	294	varMin[ivar] = binLimits[0];
	295	varMax[ivar] = binLimits[container->GetNBins(ivar)];
	296	// source
	297	if(ivar == 4){
	298	if((source>= 0) && (source<container->GetNBins(ivar))) {
	299	varMin[ivar] = binLimits[source];
	300	varMax[ivar] = binLimits[source];
	301	}
	302	}
	303	// charge
	304	if(ivar == 3) {
	305	if(charge != kAllCharge) varMin[ivar] = varMax[ivar] = charge;
	306	}
	307	// eta
	308	if(ivar == 1){
	309	for(Int_t ic = 1; ic <= container->GetAxis(1,0)->GetLast(); ic++)
	310	AliDebug(1, Form("eta bin %d, min %f, max %f\n", ic, container->GetAxis(1,0)->GetBinLowEdge(ic), container->GetAxis(1,0)->GetBinUpEdge(ic)));
	311	if(fEtaSelected){
	312	varMin[ivar] = fEtaRange[0];
	313	varMax[ivar] = fEtaRange[1];
	314	}
	315	}
	316	if(fEtaSelected){
	317	fEtaRangeNorm[0] = container->GetAxis(1,0)->GetBinLowEdge(container->GetAxis(1,0)->FindBin(fEtaRange[0]));
	318	fEtaRangeNorm[1] = container->GetAxis(1,0)->GetBinUpEdge(container->GetAxis(1,0)->FindBin(fEtaRange[1]));
7483e78e	319	AliInfo(Form("Normalization done in eta range [%f,%f]\n", fEtaRangeNorm[0], fEtaRangeNorm[1]));
959ea9d8	320	}
	321	// centrality
	322	if(ivar == 5){
7483e78e	323	if((centralitylow>= 0) && (centralitylow<container->GetNBins(ivar)) && (centralityhigh>= 0) && (centralityhigh<container->GetNBins(ivar)) && doCentralityProjection) {
	324	varMin[ivar] = binLimits[centralitylow]+ kVerySmall;
	325	varMax[ivar] = binLimits[centralityhigh]+ kVerySmall;
	326
	327	TAxis *axistest = container->GetAxis(5,0);
	328	AliDebug(1, Form("Number of bin in centrality direction %d\n",axistest->GetNbins()));
	329	AliDebug(1, Form("Project from %f to %f\n",binLimits[centralitylow],binLimits[centralityhigh]));
	330	Double_t lowcentrality = axistest->GetBinLowEdge(axistest->FindBin(binLimits[centralitylow]));
	331	Double_t highcentrality = axistest->GetBinUpEdge(axistest->FindBin(binLimits[centralityhigh]));
	332	AliDebug(1, Form("Low centrality %f and high centrality %f\n",lowcentrality,highcentrality));
	333	}
959ea9d8	334	}
959ea9d8	335
7483e78e	336	// Protect varmax against overflow
	337	if(TMath::Abs(varMax[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
	338	AliInfo("Protection against overflow bin");
	339	varMax[ivar] -= kVerySmall;
	340	}
	341	if(varMax[ivar] > binLimits[container->GetNBins(ivar)]){
	342	AliError("Upper limit exceeds allowed range");
	343	varMax[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
	344	}
	345
	346	// Protect varmin against overflow
	347	if(TMath::Abs(varMin[ivar] - binLimits[container->GetNBins(ivar)]) < kVerySmall){
	348	AliInfo("Protection against overflow bin");
	349	varMin[ivar] -= kVerySmall;
	350	}
	351	if(varMin[ivar] > binLimits[container->GetNBins(ivar)]){
	352	AliError("Upper limit exceeds allowed range");
	353	varMin[ivar] = binLimits[container->GetNBins(ivar)] - kVerySmall;
	354	}
	355	AliDebug(1, Form("variable %d: Settting limits to %f and %f\n", ivar, varMin[ivar], varMax[ivar]));
959ea9d8	356	delete[] binLimits;
7483e78e	357
959ea9d8	358	}
7483e78e	359
959ea9d8	360	AliCFContainer *k = container->MakeSlice(nDim, dimensions, varMin, varMax);
	361	delete[] varMin; delete[] varMax;
	362
	363	return k;
	364
	365	}
	366
	367	//_________________________________________________________________________
7483e78e	368	THnSparseF AliHFECorrectSpectrumBase::GetSlicedCorrelation(THnSparseF correlationmatrix, Int_t nDim, Int_t *dimensions,Chargetype_t charge,Int_t centralitylow, Int_t centralityhigh, Bool_t doCentralityProjection) const {
959ea9d8	369	//
	370	// Slice correlation
	371	//
	372
	373	Int_t ndimensions = correlationmatrix->GetNdimensions();
	374	//printf("Number of dimension %d correlation map\n",ndimensions);
	375	if(ndimensions < (2*nDim)) {
	376	AliError("Problem in the dimensions");
	377	return NULL;
	378	}
	379
	380	// Cut in centrality is centrality > -1
	381	if((5+((Int_t)(ndimensions/2.))) < ndimensions) {
	382	if((centralitylow >=0) && (centralityhigh >=0)) {
	383
	384	TAxis *axiscentrality0 = correlationmatrix->GetAxis(5);
	385	TAxis *axiscentrality1 = correlationmatrix->GetAxis(5+((Int_t)(ndimensions/2.)));
	386
	387	Int_t bins0 = axiscentrality0->GetNbins();
	388	Int_t bins1 = axiscentrality1->GetNbins();
	389
	390	AliDebug(1, Form("Number of centrality bins: %d and %d\n",bins0,bins1));
	391	if(bins0 != bins1) {
7483e78e	392	AliError("Problem in the dimensions");
7483e78e	393	return NULL;
959ea9d8	394	}
959ea9d8	395
7483e78e	396	if((centralitylow>= 0) && (centralitylow<bins0) && (centralityhigh>= 0) && (centralityhigh<bins0) && doCentralityProjection) {
	397	axiscentrality0->SetRangeUser(centralitylow,centralityhigh);
	398	axiscentrality1->SetRangeUser(centralitylow,centralityhigh);
959ea9d8	399
7483e78e	400	Double_t lowcentrality0 = axiscentrality0->GetBinLowEdge(axiscentrality0->FindBin(centralitylow));
	401	Double_t highcentrality0 = axiscentrality0->GetBinUpEdge(axiscentrality0->FindBin(centralityhigh));
	402	Double_t lowcentrality1 = axiscentrality1->GetBinLowEdge(axiscentrality1->FindBin(centralitylow));
	403	Double_t highcentrality1 = axiscentrality1->GetBinUpEdge(axiscentrality1->FindBin(centralityhigh));
	404	AliDebug(1,Form("0 Low centrality %f and high centrality %f\n",lowcentrality0,highcentrality0));
	405	AliDebug(1,Form("1 Low centrality %f and high centrality %f\n",lowcentrality1,highcentrality1));
959ea9d8	406
	407	}
	408	}
	409	}
	410
	411	// Cut in eta > -1
	412	if(fEtaSelected){
	413	if((1+((Int_t)(ndimensions/2.))) < ndimensions) {
	414
	415	TAxis *axiseta0 = correlationmatrix->GetAxis(1);
	416	TAxis *axiseta1 = correlationmatrix->GetAxis(1+((Int_t)(ndimensions/2.)));
	417
	418	Int_t bins0 = axiseta0->GetNbins();
	419	Int_t bins1 = axiseta1->GetNbins();
	420
	421	AliDebug(1, Form("Number of eta bins: %d and %d\n",bins0,bins1));
	422	if(bins0 != bins1) {
	423	AliError("Problem in the dimensions");
	424	return NULL;
	425	}
	426
	427	axiseta0->SetRangeUser(fEtaRange[0],fEtaRange[1]);
	428	axiseta1->SetRangeUser(fEtaRange[0],fEtaRange[1]);
	429
	430	Double_t loweta0 = axiseta0->GetBinLowEdge(axiseta0->FindBin(fEtaRange[0]));
	431	Double_t higheta0 = axiseta0->GetBinUpEdge(axiseta0->FindBin(fEtaRange[1]));
	432	Double_t loweta1 = axiseta1->GetBinLowEdge(axiseta1->FindBin(fEtaRange[0]));
	433	Double_t higheta1 = axiseta1->GetBinUpEdge(axiseta1->FindBin(fEtaRange[1]));
	434	AliInfo(Form("0 Low eta %f and high eta %f\n",loweta0,higheta0));
	435	AliInfo(Form("1 Low eta %f and high eta %f\n",loweta1,higheta1));
	436
	437	}
	438	}
	439
	440	// Cut in charge
	441	if(charge != kAllCharge) {
	442	if((3+((Int_t)(ndimensions/2.))) < ndimensions) {
	443
	444	TAxis *axischarge0 = correlationmatrix->GetAxis(3);
	445	TAxis *axischarge1 = correlationmatrix->GetAxis(3+((Int_t)(ndimensions/2.)));
	446
	447	Int_t bins0 = axischarge0->GetNbins();
	448	Int_t bins1 = axischarge1->GetNbins();
	449
	450	AliDebug(1, Form("Number of charge bins: %d and %d\n",bins0,bins1));
	451	if(bins0 != bins1) {
	452	AliError("Problem in the dimensions");
	453	return NULL;
	454	}
	455
	456	axischarge0->SetRangeUser(charge,charge);
	457	axischarge1->SetRangeUser(charge,charge);
	458
	459	Double_t lowcharge0 = axischarge0->GetBinLowEdge(axischarge0->FindBin(charge));
	460	Double_t highcharge0 = axischarge0->GetBinUpEdge(axischarge0->FindBin(charge));
	461	Double_t lowcharge1 = axischarge1->GetBinLowEdge(axischarge1->FindBin(charge));
	462	Double_t highcharge1 = axischarge1->GetBinUpEdge(axischarge1->FindBin(charge));
	463	AliInfo(Form("0 Low charge %f and high charge %f\n",lowcharge0,highcharge0));
	464	AliInfo(Form("1 Low charge %f and high charge %f\n",lowcharge1,highcharge1));
	465
	466	}
	467	}
	468
	469
470	Int_t ndimensionsContainer = (Int_t) ndimensions/2;
471
472	Int_t dim = new Int_t[nDim2];
473	for(Int_t iter=0; iter < nDim; iter++){
474	dim[iter] = dimensions[iter];
475	dim[iter+nDim] = ndimensionsContainer + dimensions[iter];
476	}
477
478	THnSparseF k = (THnSparseF ) correlationmatrix->Projection(nDim*2,dim);
479
480	delete[] dim;
481	return k;
482
483	}
484	//___________________________________________________________________________
485	void AliHFECorrectSpectrumBase::CorrectFromTheWidth(TH1D *h1) const {
486	//
487	// Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
488	//
489
490	TAxis *axis = h1->GetXaxis();
491	Int_t nbinX = h1->GetNbinsX();
492
493	for(Int_t i = 1; i <= nbinX; i++) {
494
495	Double_t width = axis->GetBinWidth(i);
496	Double_t content = h1->GetBinContent(i);
497	Double_t error = h1->GetBinError(i);
498	h1->SetBinContent(i,content/width);
499	h1->SetBinError(i,error/width);
500	}
501
502	}
503
504	//___________________________________________________________________________
505	void AliHFECorrectSpectrumBase::CorrectStatErr(AliCFDataGrid *backgroundGrid) const {
506	//
507	// Correct statistical error
508	//
509
510	TH1D h1 = (TH1D)backgroundGrid->Project(0);
511	Int_t nbinX = h1->GetNbinsX();
512	Int_t bins[1];
513	for(Long_t i = 1; i <= nbinX; i++) {
514	bins[0] = i;
515	Float_t content = h1->GetBinContent(i);
516	if(content>0){
517	Float_t error = TMath::Sqrt(content);
518	backgroundGrid->SetElementError(bins, error);
519	}
520	}
521	}
522	//_________________________________________________________________________
523	TObject* AliHFECorrectSpectrumBase::GetSpectrum(const AliCFContainer * const c, Int_t step) {
524	AliCFDataGrid* data = new AliCFDataGrid("data","",*c, step);
525	return data;
526	}
527	//_________________________________________________________________________
528	TObject* AliHFECorrectSpectrumBase::GetEfficiency(const AliCFContainer * const c, Int_t step, Int_t step0){
529	//
530	// Create efficiency grid and calculate efficiency
531	// of step to step0
532	//
533	TString name("eff");
534	name += step;
535	name+= step0;
536	AliCFEffGrid* eff = new AliCFEffGrid((const char)name,"",c);
537	eff->CalculateEfficiency(step,step0);
538	return eff;
539	}