[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEtools.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.                  *
**************************************************************************/
//
// Toolkit containing various usefull things
// Usable everywhere in the hfe software package
// For more information see the cxx file
//
// Authors
//   All authors of the HFE group
//
#include <TArrayD.h>
#include <TMath.h>
#include <TParticle.h>
#include "TH1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "THnSparse.h"
#include "TAxis.h"
#include "TMath.h"
#include "TString.h"

#include "AliAODMCParticle.h"
#include "AliAODpidUtil.h"
#include "AliESDpid.h"
#include "AliLog.h"
#include "AliTPCPIDResponse.h"
#include "AliTOFPIDResponse.h"

#include "AliHFEtools.h"

ClassImp(AliHFEtools)

AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
Int_t AliHFEtools::fgLogLevel = 0;

//__________________________________________
AliHFEtools::AliHFEtools():
  TObject()
{
}

//__________________________________________
Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  //
  // Helper function for linearly binned array
  //
  Double_t *bins = new Double_t[nBins + 1];
  Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  bins[0] = ymin;
  for(Int_t ibin = 1; ibin <= nBins; ibin++)
    bins[ibin] = bins[ibin-1] + stepsize;
  return bins;
}

//__________________________________________
void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
  //
  // Helper function for linearly binned array
  //
  Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
  bins[0] = ymin;
  for(Int_t ibin = 1; ibin <= nBins; ibin++)
    bins[ibin] = bins[ibin-1] + stepsize;
}

//__________________________________________
Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
  //
  // Helper function for logartimically binned array
  //
  Double_t *bins = new Double_t[nBins+1];
  bins[0] = ymin;
  Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
  for(Int_t ibin = 1; ibin <= nBins; ibin++){
    bins[ibin] = factor * bins[ibin-1];
  }
  return bins;
}

//__________________________________________
void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
  //
  // Helper function for logartimically binned array
  //
  bins[0] = ymin;
  Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
  for(Int_t ibin = 1; ibin <= nBins; ibin++)
    bins[ibin] = factor * bins[ibin-1];
}

//_________________________________________
Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){

  // 
  // converts the axis (defined by the dimension) of THx or THnSparse
  // object to Log scale. Number of bins and bin min and bin max are preserved
  //


  if(!o){
    AliError("Input histogram is null pointer");
    return kFALSE;    
  }
  
  TAxis *axis = 0x0;
  if(o->InheritsFrom("TH1")){
    TH1 *h1 = dynamic_cast<TH1F*>(o); 
    if(h1) axis = h1->GetXaxis();
  }
  else if(o->InheritsFrom("TH2")){
    TH2 *h2 = dynamic_cast<TH2F*>(o);
    if(h2 && 0 == dim){
      axis = h2->GetXaxis();
    }
    else if(h2 && 1 == dim){
      axis = h2->GetYaxis();
    }
     else{
       AliError("Only dim = 0 or 1 possible for TH2F");
     }
  }
  else if(o->InheritsFrom("THnSparse")){
    THnSparse *hs = dynamic_cast<THnSparse*>(o);
    if(hs) axis = hs->GetAxis(dim);
  }
  else{
    AliError("Type of input object not recognized, please check your code or update this finction");
    return kFALSE;
  }
  if(!axis){
    AliError(Form("Axis '%d' could not be identified in the object \n", dim));
    return kFALSE;
  }
  
  Int_t bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  if(from <= 0){
    AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
  }
  Double_t to = axis->GetXmax();
  TArrayD newBins(bins+1);
  newBins[0] = from;
  Double_t factor = TMath::Power(to/from, 1./bins);
  for(Int_t i=1; i<=bins; ++i){
    newBins[i] = factor * newBins[i-1];
  }
  axis->Set(bins, newBins.GetArray());

  return kTRUE;
}

//__________________________________________
Float_t AliHFEtools::GetRapidity(const TParticle *part){
  //
  // return rapidity
  //
  Float_t rapidity;
  if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
  else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
  return rapidity;
}

//__________________________________________
Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
  // return rapidity

  Float_t rapidity;        
  if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999; 
  else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz()))); 
  return rapidity;
}

//__________________________________________
AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
  //
  // Get the default PID as singleton instance
  //
  if(!fgDefaultPID){

    if(isESD) fgDefaultPID = new AliESDpid;
    else fgDefaultPID = new AliAODpidUtil;
    Double_t tres = isMC ? 80. : 130.;
    fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);

    // TPC Bethe Bloch parameters
    Double_t alephParameters[5];
    if(isMC){
      // simulation
      alephParameters[0] = 2.15898e+00/50.;
      alephParameters[1] = 1.75295e+01;
      alephParameters[2] = 3.40030e-09;
      alephParameters[3] = 1.96178e+00;
      alephParameters[4] = 3.91720e+00;
    } else {
      alephParameters[0] = 0.0283086/0.97;
      //alephParameters[0] = 0.0283086;
      alephParameters[1] = 2.63394e+01;
      alephParameters[2] = 5.04114e-11;
      alephParameters[3] = 2.12543e+00;
      alephParameters[4] = 4.88663e+00;
    }
    fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
    fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);

  }
  if(fgLogLevel){
    printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
    printf("Error - Arrrrrrrrr...\n");
    printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
    printf("Error - TOF PID is adapted to Monte Carlo\n");
  }
  return fgDefaultPID;
}


//__________________________________________
void AliHFEtools::DestroyDefaultPID(){
  //
  // Destroy default PID object if existing
  //
  if(fgDefaultPID) delete fgDefaultPID;
  fgDefaultPID = NULL;
}

//__________________________________________
Int_t AliHFEtools::GetPdg(const AliVParticle *track){
  // 
  // Get MC PDG code (only MC particles supported)
  //
  Int_t pdg = 0;
  if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
    const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
    pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
  } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
    const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
    pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
  }
  return pdg;
}

//__________________________________________
Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
  // 
  // Helper function to convert MC PID codes into AliPID codes
  //
  Int_t pid = -1;
  switch(TMath::Abs(pdg)){
    case 11: pid = AliPID::kElectron; break;
    case 13: pid = AliPID::kMuon; break;
    case 211: pid = AliPID::kPion; break;
    case 321: pid = AliPID::kKaon; break;
    case 2212: pid = AliPID::kProton; break;
    default: pid = -1; break;
  };
  return pid;
}

//__________________________________________
TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
{
    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
    // specified. The others bins are calculated automatically. Supports also Graphs
    // with non constant x steps. The axis of the Graph can be exchanged if
    // exchange=kTRUE (modifies the Graph).

    TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
    if( result == 0) return result;

    //------------------------------------------
    // setup the final hist
    Int_t nBinX = g->GetN();
    Double_t* err = g->GetEY();           // error y is still ok even if exchanged
    for(Int_t i = 0; i < nBinX; i ++){
	result->SetBinError(i + 1, err[i]);
    }
    if(exchange){
        AliHFEtools::ExchangeXYGraph(g);        // undo  what has been done in GraphToHist
	AliHFEtools::ExchangeXYGraphErrors(g);  // now exchange including errors
    }

    return result;
}

//__________________________________________
Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
{
    // exchanges x-values and y-values.
    if(g==0) return kFALSE;
    Int_t nbin=g->GetN();
    Double_t x,y;
    for(Int_t i = 0; i < nbin; i ++)
    {
        g->GetPoint(i,x,y);
        g->SetPoint(i,y,x);
    }

    return kTRUE;
}

//__________________________________________
Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
{
    // exchanges x-values and y-values and
    // corresponding errors.
    if(g==0) return kFALSE;
    Int_t nbin=g->GetN();
    Double_t x,y;
    Double_t ex,ey;
    for(Int_t i = 0; i < nbin; i ++)
    {
        g->GetPoint(i,x,y);
        ex = g->GetErrorX(i);
        ey = g->GetErrorY(i);
        g->SetPoint(i,y,x);
        g->SetPointError(i,ey,ex);
    }

    return kTRUE;

}

//__________________________________________
TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
{
    // Creates a TH1D from TGraph g. The binwidth of the first bin has to
    // specified. The others bins are calculated automatically. Supports also Graphs
    // with non constant x steps. The axis of the Graph can be exchanged if
    // exchange=kTRUE (modifies the Graph).


    TH1D* result = 0;
    if(g == 0)              return result;
    if(firstBinWidth == -1) return result;
    TString myname="";
    myname = g->GetName();
    myname += "_graph";
    if(exchange) AliHFEtools::ExchangeXYGraph(g);

    Int_t nBinX = g->GetN();
    Double_t* x = g->GetX();
    Double_t* y = g->GetY();

    if(nBinX < 1) return result;

    //------------------------------------------
    // create the Matrix for the equation system
    // and init the values

    Int_t nDim = nBinX - 1;
    TMatrixD a(nDim,nDim);
    TMatrixD b(nDim,1);

    Double_t* aA = a.GetMatrixArray();
    Double_t* aB = b.GetMatrixArray();
    memset(aA,0,nDim * nDim * sizeof(Double_t));
    memset(aB,0,nDim * sizeof(Double_t));
    //------------------------------------------

    //------------------------------------------
    // setup equation system
    // width for 1st bin is given therefore
    // we shift bin parameter (column) by one to the left
    // to reduce the matrix size

    Double_t* xAxis = new Double_t [nBinX + 1];
    Double_t* binW  = new Double_t [nBinX ];
    binW[0] = firstBinWidth;

    aB[0] = x[1] - x[0]  - 0.5 * binW[0];
    aA[0] = 0.5;

    for(Int_t col = 1; col < nDim ; col ++)
    {
	Int_t row = col;
	aB[col] = x[col + 1] - x[ col ];
	aA[row * nDim + col - 1 ] = 0.5;
	aA[row * nDim + col     ] = 0.5;
    }
    //------------------------------------------

    //------------------------------------------
    // solve the equations
    a.Invert();
    TMatrixD c = a * b;
    //------------------------------------------

    //------------------------------------------
    // calculate the bin boundaries
    xAxis[0] = x[0] - 0.5 * binW[0];
    memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
    for(Int_t col = 0; col < nBinX ; col ++) {
	xAxis[col + 1] = x[col] + 0.5 * binW[col];
    }
    //------------------------------------------

    //------------------------------------------
    // setup the final hist
    result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
    for(Int_t i = 0; i < nBinX; i ++){
	result->SetBinContent(i + 1, y[i]);
    }
    result->SetMarkerColor(markercolor);
    result->SetMarkerStyle(markerstyle);
    result->SetMarkerSize(markersize);
    //------------------------------------------

    delete [] xAxis;
    delete [] binW;


    return result;
}
Commit	Line	Data
70da6c5a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	//
	16	// Toolkit containing various usefull things
	17	// Usable everywhere in the hfe software package
	18	// For more information see the cxx file
	19	//
	20	// Authors
	21	// All authors of the HFE group
	22	//
0e30407a	23	#include <TArrayD.h>
70da6c5a	24	#include <TMath.h>
	25	#include <TParticle.h>
	26	#include "TH1.h"
a8ef1999	27	#include "TH1D.h"
70da6c5a	28	#include "TH2.h"
a8ef1999	29	#include "TGraph.h"
a8ef1999	30	#include "TGraphErrors.h"
70da6c5a	31	#include "THnSparse.h"
70da6c5a	32	#include "TAxis.h"
a8ef1999	33	#include "TMath.h"
a8ef1999	34	#include "TString.h"
70da6c5a	35
70da6c5a	36	#include "AliAODMCParticle.h"
3a72645a	37	#include "AliAODpidUtil.h"
faee3b18	38	#include "AliESDpid.h"
70da6c5a	39	#include "AliLog.h"
3a72645a	40	#include "AliTPCPIDResponse.h"
faee3b18	41	#include "AliTOFPIDResponse.h"
70da6c5a	42
	43	#include "AliHFEtools.h"
	44
	45	ClassImp(AliHFEtools)
	46
8c1c76e9	47	AliPIDResponse *AliHFEtools::fgDefaultPID = NULL;
69ac0e6f	48	Int_t AliHFEtools::fgLogLevel = 0;
faee3b18	49
70da6c5a	50	//__________________________________________
	51	AliHFEtools::AliHFEtools():
	52	TObject()
	53	{
	54	}
	55
	56	//__________________________________________
	57	Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
	58	//
	59	// Helper function for linearly binned array
	60	//
	61	Double_t *bins = new Double_t[nBins + 1];
	62	Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
	63	bins[0] = ymin;
	64	for(Int_t ibin = 1; ibin <= nBins; ibin++)
	65	bins[ibin] = bins[ibin-1] + stepsize;
	66	return bins;
	67	}
	68
0e30407a	69	//__________________________________________
	70	void AliHFEtools::FillLinearBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
	71	//
	72	// Helper function for linearly binned array
	73	//
	74	Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
	75	bins[0] = ymin;
	76	for(Int_t ibin = 1; ibin <= nBins; ibin++)
	77	bins[ibin] = bins[ibin-1] + stepsize;
	78	}
	79
70da6c5a	80	//__________________________________________
	81	Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
	82	//
	83	// Helper function for logartimically binned array
	84	//
	85	Double_t *bins = new Double_t[nBins+1];
	86	bins[0] = ymin;
	87	Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
	88	for(Int_t ibin = 1; ibin <= nBins; ibin++){
	89	bins[ibin] = factor * bins[ibin-1];
	90	}
	91	return bins;
	92	}
	93
0e30407a	94	//__________________________________________
	95	void AliHFEtools::FillLogarithmicBinning(TArrayD &bins, Int_t nBins, Double_t ymin, Double_t ymax){
	96	//
	97	// Helper function for logartimically binned array
	98	//
	99	bins[0] = ymin;
	100	Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
	101	for(Int_t ibin = 1; ibin <= nBins; ibin++)
	102	bins[ibin] = factor * bins[ibin-1];
	103	}
	104
70da6c5a	105	//_________________________________________
	106	Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
	107
	108	//
	109	// converts the axis (defined by the dimension) of THx or THnSparse
	110	// object to Log scale. Number of bins and bin min and bin max are preserved
	111	//
	112
	113
	114	if(!o){
	115	AliError("Input histogram is null pointer");
	116	return kFALSE;
	117	}
	118
	119	TAxis *axis = 0x0;
	120	if(o->InheritsFrom("TH1")){
bf892a6a	121	TH1 h1 = dynamic_cast<TH1F>(o);
bf892a6a	122	if(h1) axis = h1->GetXaxis();
70da6c5a	123	}
70da6c5a	124	else if(o->InheritsFrom("TH2")){
bf892a6a	125	TH2 h2 = dynamic_cast<TH2F>(o);
	126	if(h2 && 0 == dim){
	127	axis = h2->GetXaxis();
70da6c5a	128	}
bf892a6a	129	else if(h2 && 1 == dim){
bf892a6a	130	axis = h2->GetYaxis();
70da6c5a	131	}
	132	else{
	133	AliError("Only dim = 0 or 1 possible for TH2F");
	134	}
	135	}
	136	else if(o->InheritsFrom("THnSparse")){
bf892a6a	137	THnSparse hs = dynamic_cast<THnSparse>(o);
bf892a6a	138	if(hs) axis = hs->GetAxis(dim);
70da6c5a	139	}
	140	else{
	141	AliError("Type of input object not recognized, please check your code or update this finction");
	142	return kFALSE;
	143	}
	144	if(!axis){
	145	AliError(Form("Axis '%d' could not be identified in the object \n", dim));
	146	return kFALSE;
	147	}
	148
	149	Int_t bins = axis->GetNbins();
	150
	151	Double_t from = axis->GetXmin();
	152	if(from <= 0){
	153	AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
	154	}
	155	Double_t to = axis->GetXmax();
0e30407a	156	TArrayD newBins(bins+1);
70da6c5a	157	newBins[0] = from;
	158	Double_t factor = TMath::Power(to/from, 1./bins);
	159	for(Int_t i=1; i<=bins; ++i){
	160	newBins[i] = factor * newBins[i-1];
	161	}
0e30407a	162	axis->Set(bins, newBins.GetArray());
70da6c5a	163
	164	return kTRUE;
	165	}
	166
	167	//__________________________________________
3a72645a	168	Float_t AliHFEtools::GetRapidity(const TParticle *part){
70da6c5a	169	//
	170	// return rapidity
	171	//
	172	Float_t rapidity;
	173	if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
	174	else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
	175	return rapidity;
	176	}
	177
	178	//__________________________________________
3a72645a	179	Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
70da6c5a	180	// return rapidity
	181
	182	Float_t rapidity;
	183	if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
	184	else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
	185	return rapidity;
	186	}
	187
faee3b18	188	//__________________________________________
8c1c76e9	189	AliPIDResponse* AliHFEtools::GetDefaultPID(Bool_t isMC, Bool_t isESD){
faee3b18	190	//
	191	// Get the default PID as singleton instance
	192	//
	193	if(!fgDefaultPID){
8c1c76e9	194
	195	if(isESD) fgDefaultPID = new AliESDpid;
	196	else fgDefaultPID = new AliAODpidUtil;
faee3b18	197	Double_t tres = isMC ? 80. : 130.;
	198	fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
	199
	200	// TPC Bethe Bloch parameters
	201	Double_t alephParameters[5];
	202	if(isMC){
	203	// simulation
	204	alephParameters[0] = 2.15898e+00/50.;
	205	alephParameters[1] = 1.75295e+01;
	206	alephParameters[2] = 3.40030e-09;
	207	alephParameters[3] = 1.96178e+00;
	208	alephParameters[4] = 3.91720e+00;
	209	} else {
	210	alephParameters[0] = 0.0283086/0.97;
	211	//alephParameters[0] = 0.0283086;
	212	alephParameters[1] = 2.63394e+01;
	213	alephParameters[2] = 5.04114e-11;
	214	alephParameters[3] = 2.12543e+00;
	215	alephParameters[4] = 4.88663e+00;
	216	}
	217	fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
67fe7bd0	218	fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
faee3b18	219
	220	}
	221	if(fgLogLevel){
67fe7bd0	222	printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
	223	printf("Error - Arrrrrrrrr...\n");
	224	printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
	225	printf("Error - TOF PID is adapted to Monte Carlo\n");
faee3b18	226	}
	227	return fgDefaultPID;
	228	}
	229
3a72645a	230
faee3b18	231	//__________________________________________
	232	void AliHFEtools::DestroyDefaultPID(){
	233	//
	234	// Destroy default PID object if existing
	235	//
	236	if(fgDefaultPID) delete fgDefaultPID;
	237	fgDefaultPID = NULL;
3a72645a	238	}
	239
	240	//__________________________________________
	241	Int_t AliHFEtools::GetPdg(const AliVParticle *track){
	242	//
	243	// Get MC PDG code (only MC particles supported)
	244	//
	245	Int_t pdg = 0;
	246	if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
	247	const AliMCParticle mctrack = dynamic_cast<const AliMCParticle >(track);
e3ae862b	248	pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
3a72645a	249	} else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
3a72645a	250	const AliAODMCParticle aodmctrack = dynamic_cast<const AliAODMCParticle >(track);
e3ae862b	251	pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
3a72645a	252	}
	253	return pdg;
	254	}
	255
	256	//__________________________________________
	257	Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
	258	//
	259	// Helper function to convert MC PID codes into AliPID codes
	260	//
	261	Int_t pid = -1;
	262	switch(TMath::Abs(pdg)){
	263	case 11: pid = AliPID::kElectron; break;
	264	case 13: pid = AliPID::kMuon; break;
	265	case 211: pid = AliPID::kPion; break;
	266	case 321: pid = AliPID::kKaon; break;
	267	case 2212: pid = AliPID::kProton; break;
	268	default: pid = -1; break;
	269	};
	270	return pid;
faee3b18	271	}
a8ef1999	272
	273	//__________________________________________
	274	TH1D* AliHFEtools::GraphErrorsToHist(TGraphErrors* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
	275	{
	276	// Creates a TH1D from TGraph g. The binwidth of the first bin has to
	277	// specified. The others bins are calculated automatically. Supports also Graphs
	278	// with non constant x steps. The axis of the Graph can be exchanged if
	279	// exchange=kTRUE (modifies the Graph).
	280
	281	TH1D* result = GraphToHist(g,firstBinWidth,exchange, markerstyle,markercolor,markersize);
	282	if( result == 0) return result;
	283
	284	//------------------------------------------
	285	// setup the final hist
	286	Int_t nBinX = g->GetN();
	287	Double_t* err = g->GetEY(); // error y is still ok even if exchanged
	288	for(Int_t i = 0; i < nBinX; i ++){
	289	result->SetBinError(i + 1, err[i]);
	290	}
	291	if(exchange){
	292	AliHFEtools::ExchangeXYGraph(g); // undo what has been done in GraphToHist
	293	AliHFEtools::ExchangeXYGraphErrors(g); // now exchange including errors
	294	}
	295
	296	return result;
	297	}
	298
	299	//__________________________________________
	300	Bool_t AliHFEtools::ExchangeXYGraph(TGraph* g)
	301	{
	302	// exchanges x-values and y-values.
	303	if(g==0) return kFALSE;
	304	Int_t nbin=g->GetN();
	305	Double_t x,y;
	306	for(Int_t i = 0; i < nbin; i ++)
	307	{
	308	g->GetPoint(i,x,y);
	309	g->SetPoint(i,y,x);
	310	}
	311
	312	return kTRUE;
	313	}
	314
	315	//__________________________________________
	316	Bool_t AliHFEtools::ExchangeXYGraphErrors(TGraphErrors* g)
	317	{
	318	// exchanges x-values and y-values and
	319	// corresponding errors.
	320	if(g==0) return kFALSE;
	321	Int_t nbin=g->GetN();
	322	Double_t x,y;
	323	Double_t ex,ey;
	324	for(Int_t i = 0; i < nbin; i ++)
	325	{
	326	g->GetPoint(i,x,y);
	327	ex = g->GetErrorX(i);
	328	ey = g->GetErrorY(i);
	329	g->SetPoint(i,y,x);
	330	g->SetPointError(i,ey,ex);
	331	}
	332
	333	return kTRUE;
	334
	335	}
336
337	//__________________________________________
338	TH1D* AliHFEtools::GraphToHist(TGraph* g, Double_t firstBinWidth, Bool_t exchange, Int_t markerstyle,Int_t markercolor,Float_t markersize)
339	{
340	// Creates a TH1D from TGraph g. The binwidth of the first bin has to
341	// specified. The others bins are calculated automatically. Supports also Graphs
342	// with non constant x steps. The axis of the Graph can be exchanged if
343	// exchange=kTRUE (modifies the Graph).
344
345
346	TH1D* result = 0;
347	if(g == 0) return result;
348	if(firstBinWidth == -1) return result;
349	TString myname="";
350	myname = g->GetName();
351	myname += "_graph";
352	if(exchange) AliHFEtools::ExchangeXYGraph(g);
353
354	Int_t nBinX = g->GetN();
355	Double_t* x = g->GetX();
356	Double_t* y = g->GetY();
357
358	if(nBinX < 1) return result;
359
360	//------------------------------------------
361	// create the Matrix for the equation system
362	// and init the values
363
364	Int_t nDim = nBinX - 1;
365	TMatrixD a(nDim,nDim);
366	TMatrixD b(nDim,1);
367
368	Double_t* aA = a.GetMatrixArray();
369	Double_t* aB = b.GetMatrixArray();
370	memset(aA,0,nDim * nDim * sizeof(Double_t));
371	memset(aB,0,nDim * sizeof(Double_t));
372	//------------------------------------------
373
374	//------------------------------------------
375	// setup equation system
376	// width for 1st bin is given therefore
377	// we shift bin parameter (column) by one to the left
378	// to reduce the matrix size
379
380	Double_t* xAxis = new Double_t [nBinX + 1];
381	Double_t* binW = new Double_t [nBinX ];
382	binW[0] = firstBinWidth;
383
384	aB[0] = x[1] - x[0] - 0.5 * binW[0];
385	aA[0] = 0.5;
386
387	for(Int_t col = 1; col < nDim ; col ++)
388	{
389	Int_t row = col;
390	aB[col] = x[col + 1] - x[ col ];
391	aA[row * nDim + col - 1 ] = 0.5;
392	aA[row * nDim + col ] = 0.5;
393	}
394	//------------------------------------------
395
396	//------------------------------------------
397	// solve the equations
398	a.Invert();
399	TMatrixD c = a * b;
400	//------------------------------------------
401
402	//------------------------------------------
403	// calculate the bin boundaries
404	xAxis[0] = x[0] - 0.5 * binW[0];
405	memcpy(&binW[1],c.GetMatrixArray(),nDim * sizeof(Double_t));
406	for(Int_t col = 0; col < nBinX ; col ++) {
407	xAxis[col + 1] = x[col] + 0.5 * binW[col];
408	}
409	//------------------------------------------
410
411	//------------------------------------------
412	// setup the final hist
413	result = new TH1D(myname.Data(),myname.Data(),nBinX, xAxis);
414	for(Int_t i = 0; i < nBinX; i ++){
415	result->SetBinContent(i + 1, y[i]);
416	}
417	result->SetMarkerColor(markercolor);
418	result->SetMarkerStyle(markerstyle);
419	result->SetMarkerSize(markersize);
420	//------------------------------------------
421
422	delete [] xAxis;
423	delete [] binW;
424
425
426	return result;
427	}