[u/mrichter/AliRoot.git] / PWG3 / 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.                  *
**************************************************************************/

/* $Id$ */

//
// 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 <TMath.h>
#include <TParticle.h>
#include "TH1.h"
#include "TH2.h"
#include "THnSparse.h"
#include "TAxis.h"

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

#include "AliHFEtools.h"

ClassImp(AliHFEtools)

AliESDpid *AliHFEtools::fgDefaultPID = NULL;
AliAODpidUtil *AliHFEtools::fgDefaultPIDaod = 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;
}

//__________________________________________
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;
}

//_________________________________________
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();
  Double_t *newBins = new Double_t[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);
  delete[] newBins;

  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;
}

//__________________________________________
AliESDpid* AliHFEtools::GetDefaultPID(Bool_t isMC){
  //
  // Get the default PID as singleton instance
  //
  if(!fgDefaultPID){
    fgDefaultPID = new AliESDpid;
    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;
}

//__________________________________________
AliAODpidUtil* AliHFEtools::GetDefaultAODPID(Bool_t isMC){
  //
  // Get the default PID as singleton instance
  //
  if(!fgDefaultPID){
    fgDefaultPIDaod = new AliAODpidUtil;
    Double_t tres = isMC ? 80. : 130.;
    fgDefaultPIDaod->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;
    }
    fgDefaultPIDaod->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
    fgDefaultPIDaod->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 fgDefaultPIDaod;
}

//__________________________________________
void AliHFEtools::DestroyDefaultPID(){
  //
  // Destroy default PID object if existing
  //
  if(fgDefaultPID) delete fgDefaultPID;
  fgDefaultPID = NULL;
  if(fgDefaultPIDaod) delete fgDefaultPIDaod;
  fgDefaultPIDaod = 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;
}
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	**************************************************************************/
27de2dfb	15
	16	/* $Id$ */
	17
70da6c5a	18	//
	19	// Toolkit containing various usefull things
	20	// Usable everywhere in the hfe software package
	21	// For more information see the cxx file
	22	//
	23	// Authors
	24	// All authors of the HFE group
	25	//
	26	#include <TMath.h>
	27	#include <TParticle.h>
	28	#include "TH1.h"
	29	#include "TH2.h"
	30	#include "THnSparse.h"
	31	#include "TAxis.h"
	32
	33	#include "AliAODMCParticle.h"
3a72645a	34	#include "AliAODpidUtil.h"
faee3b18	35	#include "AliESDpid.h"
70da6c5a	36	#include "AliLog.h"
3a72645a	37	#include "AliTPCPIDResponse.h"
faee3b18	38	#include "AliTOFPIDResponse.h"
70da6c5a	39
	40	#include "AliHFEtools.h"
	41
	42	ClassImp(AliHFEtools)
	43
faee3b18	44	AliESDpid *AliHFEtools::fgDefaultPID = NULL;
3a72645a	45	AliAODpidUtil *AliHFEtools::fgDefaultPIDaod = NULL;
69ac0e6f	46	Int_t AliHFEtools::fgLogLevel = 0;
faee3b18	47
70da6c5a	48	//__________________________________________
	49	AliHFEtools::AliHFEtools():
	50	TObject()
	51	{
	52	}
	53
	54	//__________________________________________
	55	Double_t *AliHFEtools::MakeLinearBinning(Int_t nBins, Double_t ymin, Double_t ymax){
	56	//
	57	// Helper function for linearly binned array
	58	//
	59	Double_t *bins = new Double_t[nBins + 1];
	60	Double_t stepsize = (ymax - ymin) / static_cast<Double_t>(nBins);
	61	bins[0] = ymin;
	62	for(Int_t ibin = 1; ibin <= nBins; ibin++)
	63	bins[ibin] = bins[ibin-1] + stepsize;
	64	return bins;
	65	}
	66
	67	//__________________________________________
	68	Double_t *AliHFEtools::MakeLogarithmicBinning(Int_t nBins, Double_t ymin, Double_t ymax){
	69	//
	70	// Helper function for logartimically binned array
	71	//
	72	Double_t *bins = new Double_t[nBins+1];
	73	bins[0] = ymin;
	74	Double_t factor = TMath::Power(ymax/ymin, 1./nBins);
	75	for(Int_t ibin = 1; ibin <= nBins; ibin++){
	76	bins[ibin] = factor * bins[ibin-1];
	77	}
	78	return bins;
	79	}
	80
	81	//_________________________________________
	82	Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
	83
	84	//
	85	// converts the axis (defined by the dimension) of THx or THnSparse
	86	// object to Log scale. Number of bins and bin min and bin max are preserved
	87	//
	88
	89
	90	if(!o){
	91	AliError("Input histogram is null pointer");
	92	return kFALSE;
	93	}
	94
	95	TAxis *axis = 0x0;
	96	if(o->InheritsFrom("TH1")){
bf892a6a	97	TH1 h1 = dynamic_cast<TH1F>(o);
bf892a6a	98	if(h1) axis = h1->GetXaxis();
70da6c5a	99	}
70da6c5a	100	else if(o->InheritsFrom("TH2")){
bf892a6a	101	TH2 h2 = dynamic_cast<TH2F>(o);
	102	if(h2 && 0 == dim){
	103	axis = h2->GetXaxis();
70da6c5a	104	}
bf892a6a	105	else if(h2 && 1 == dim){
bf892a6a	106	axis = h2->GetYaxis();
70da6c5a	107	}
	108	else{
	109	AliError("Only dim = 0 or 1 possible for TH2F");
	110	}
	111	}
	112	else if(o->InheritsFrom("THnSparse")){
bf892a6a	113	THnSparse hs = dynamic_cast<THnSparse>(o);
bf892a6a	114	if(hs) axis = hs->GetAxis(dim);
70da6c5a	115	}
	116	else{
	117	AliError("Type of input object not recognized, please check your code or update this finction");
	118	return kFALSE;
	119	}
	120	if(!axis){
	121	AliError(Form("Axis '%d' could not be identified in the object \n", dim));
	122	return kFALSE;
	123	}
	124
	125	Int_t bins = axis->GetNbins();
	126
	127	Double_t from = axis->GetXmin();
	128	if(from <= 0){
	129	AliError(Form("Log binning not possible for this axis [min = %f]\n", from));
	130	}
	131	Double_t to = axis->GetXmax();
	132	Double_t *newBins = new Double_t[bins+1];
	133	newBins[0] = from;
	134	Double_t factor = TMath::Power(to/from, 1./bins);
	135	for(Int_t i=1; i<=bins; ++i){
	136	newBins[i] = factor * newBins[i-1];
	137	}
	138	axis->Set(bins, newBins);
bf892a6a	139	delete[] newBins;
70da6c5a	140
	141	return kTRUE;
	142	}
	143
	144	//__________________________________________
3a72645a	145	Float_t AliHFEtools::GetRapidity(const TParticle *part){
70da6c5a	146	//
	147	// return rapidity
	148	//
	149	Float_t rapidity;
	150	if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999;
	151	else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz())));
	152	return rapidity;
	153	}
	154
	155	//__________________________________________
3a72645a	156	Float_t AliHFEtools::GetRapidity(const AliAODMCParticle *part){
70da6c5a	157	// return rapidity
	158
	159	Float_t rapidity;
	160	if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999;
	161	else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz())));
	162	return rapidity;
	163	}
	164
faee3b18	165	//__________________________________________
	166	AliESDpid* AliHFEtools::GetDefaultPID(Bool_t isMC){
	167	//
	168	// Get the default PID as singleton instance
	169	//
	170	if(!fgDefaultPID){
	171	fgDefaultPID = new AliESDpid;
	172	Double_t tres = isMC ? 80. : 130.;
	173	fgDefaultPID->GetTOFResponse().SetTimeResolution(tres);
	174
	175	// TPC Bethe Bloch parameters
	176	Double_t alephParameters[5];
	177	if(isMC){
	178	// simulation
	179	alephParameters[0] = 2.15898e+00/50.;
	180	alephParameters[1] = 1.75295e+01;
	181	alephParameters[2] = 3.40030e-09;
	182	alephParameters[3] = 1.96178e+00;
	183	alephParameters[4] = 3.91720e+00;
	184	} else {
	185	alephParameters[0] = 0.0283086/0.97;
	186	//alephParameters[0] = 0.0283086;
	187	alephParameters[1] = 2.63394e+01;
	188	alephParameters[2] = 5.04114e-11;
	189	alephParameters[3] = 2.12543e+00;
	190	alephParameters[4] = 4.88663e+00;
	191	}
	192	fgDefaultPID->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
67fe7bd0	193	fgDefaultPID->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
faee3b18	194
	195	}
	196	if(fgLogLevel){
67fe7bd0	197	printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
	198	printf("Error - Arrrrrrrrr...\n");
	199	printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
	200	printf("Error - TOF PID is adapted to Monte Carlo\n");
faee3b18	201	}
	202	return fgDefaultPID;
	203	}
	204
3a72645a	205	//__________________________________________
	206	AliAODpidUtil* AliHFEtools::GetDefaultAODPID(Bool_t isMC){
	207	//
	208	// Get the default PID as singleton instance
	209	//
	210	if(!fgDefaultPID){
	211	fgDefaultPIDaod = new AliAODpidUtil;
	212	Double_t tres = isMC ? 80. : 130.;
	213	fgDefaultPIDaod->GetTOFResponse().SetTimeResolution(tres);
	214
	215	// TPC Bethe Bloch parameters
	216	Double_t alephParameters[5];
	217	if(isMC){
	218	// simulation
	219	alephParameters[0] = 2.15898e+00/50.;
	220	alephParameters[1] = 1.75295e+01;
	221	alephParameters[2] = 3.40030e-09;
	222	alephParameters[3] = 1.96178e+00;
	223	alephParameters[4] = 3.91720e+00;
	224	} else {
	225	alephParameters[0] = 0.0283086/0.97;
	226	//alephParameters[0] = 0.0283086;
	227	alephParameters[1] = 2.63394e+01;
	228	alephParameters[2] = 5.04114e-11;
	229	alephParameters[3] = 2.12543e+00;
	230	alephParameters[4] = 4.88663e+00;
	231	}
	232	fgDefaultPIDaod->GetTPCResponse().SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]);
	233	fgDefaultPIDaod->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
	234	}
	235	if(fgLogLevel){
	236	printf("Error - You are using the default PID: You should use the PID coming from the tender\n");
	237	printf("Error - Arrrrrrrrr...\n");
	238	printf("Error - Please rethink your program logic. Using default PID is really dangerous\n");
	239	printf("Error - TOF PID is adapted to Monte Carlo\n");
	240	}
	241	return fgDefaultPIDaod;
	242	}
	243
faee3b18	244	//__________________________________________
	245	void AliHFEtools::DestroyDefaultPID(){
	246	//
	247	// Destroy default PID object if existing
	248	//
	249	if(fgDefaultPID) delete fgDefaultPID;
	250	fgDefaultPID = NULL;
3a72645a	251	if(fgDefaultPIDaod) delete fgDefaultPIDaod;
	252	fgDefaultPIDaod = NULL;
	253	}
	254
	255	//__________________________________________
	256	Int_t AliHFEtools::GetPdg(const AliVParticle *track){
	257	//
	258	// Get MC PDG code (only MC particles supported)
	259	//
	260	Int_t pdg = 0;
	261	if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
	262	const AliMCParticle mctrack = dynamic_cast<const AliMCParticle >(track);
e3ae862b	263	pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
3a72645a	264	} else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
3a72645a	265	const AliAODMCParticle aodmctrack = dynamic_cast<const AliAODMCParticle >(track);
e3ae862b	266	pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
3a72645a	267	}
	268	return pdg;
	269	}
	270
	271	//__________________________________________
	272	Int_t AliHFEtools::PDG2AliPID(Int_t pdg){
	273	//
	274	// Helper function to convert MC PID codes into AliPID codes
	275	//
	276	Int_t pid = -1;
	277	switch(TMath::Abs(pdg)){
	278	case 11: pid = AliPID::kElectron; break;
	279	case 13: pid = AliPID::kMuon; break;
	280	case 211: pid = AliPID::kPion; break;
	281	case 321: pid = AliPID::kKaon; break;
	282	case 2212: pid = AliPID::kProton; break;
	283	default: pid = -1; break;
	284	};
	285	return pid;
faee3b18	286	}