[u/mrichter/AliRoot.git] / ITS / AliITSdEdxSamples.cxx

/**************************************************************************
 * Copyright(c) 2009-2012, 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$ */


///////////////////////////////////////////////////////////////////
//                                                               //
// Class to store information for PID with ITS                   //
// and truncated mean computation methods                        //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

#include "AliITSPidParams.h"
#include "AliITSdEdxSamples.h"
#include "AliLog.h"
#include <TMath.h>

ClassImp(AliITSdEdxSamples)

//______________________________________________________________________
AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
  fNSamples(0),
  fClusterMap(0),
  fP(0.),
  fParticleSpecie(0),
  fLayersForPid(0xFFFF)
{
  // Default constructor
  for(Int_t i=0; i<kMaxSamples; i++){
    fdESamples[i]=0.;
    fdxSamples[i]=0.;
    fPAtSample[i]=0.;
  }
}

//______________________________________________________________________
AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) :
  TObject(),
  fNSamples(nSamples),
  fClusterMap(0),
  fP(mom),
  fParticleSpecie(specie),
  fLayersForPid(0xFFFF)
{
  // Standard constructor
  SetdESamples(nSamples,esamples);
  SetdxSamples(nSamples,xsamples);
  SetClusterMapFromdE();
}

//______________________________________________________________________
AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) :
  TObject(),
  fNSamples(source.fNSamples),
  fClusterMap(source.fClusterMap),
  fP(source.fP),
  fParticleSpecie(source.fParticleSpecie),
  fLayersForPid(source.fLayersForPid)
{
  // Copy constructor
  for(Int_t i=0; i<kMaxSamples; i++){
    fdESamples[i]=source.GetdESample(i);
    fdxSamples[i]=source.GetdxSample(i);
    fPAtSample[i]=source.GetMomentumAtSample(i);
  }
}
//______________________________________________________________________
void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){
  // Set the samples

  if(nSamples>kMaxSamples){
    AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples));    
    fNSamples=kMaxSamples;
  }else{
    fNSamples=nSamples;
  }
  for(Int_t i=0; i<fNSamples; i++) fdESamples[i]=samples[i];
  for(Int_t i=fNSamples; i<kMaxSamples; i++) fdESamples[i]=0.;
  return;
}
//______________________________________________________________________
void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){
  // Set the samples

  if(nSamples>kMaxSamples){
    AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples))
    fNSamples=kMaxSamples;
  }else{
    fNSamples=nSamples;
  }
  for(Int_t i=0; i<fNSamples; i++) fdxSamples[i]=samples[i];
  for(Int_t i=fNSamples; i<kMaxSamples; i++) fdxSamples[i]=0.;
  return;
}

//______________________________________________________________________
void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){
  // Set the samples
  SetdESamples(nSamples,esamples);
  SetdxSamples(nSamples,xsamples);
  for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
  for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.;
  return;
}
//______________________________________________________________________
void AliITSdEdxSamples::SetLayerSample(Int_t iLayer, Bool_t haspoint, Double_t dE, Double_t dx, Double_t p){
  // set info from single layer
  if(haspoint){
    SetPointOnLayer(iLayer);
    fdESamples[iLayer]=dE; 
    fdxSamples[iLayer]=dx; 
    fPAtSample[iLayer]=p;
  }else{
    if(HasPointOnLayer(iLayer)) fClusterMap-=(1<<iLayer);
    fdESamples[iLayer]=0.; 
    fdxSamples[iLayer]=0.; 
    fPAtSample[iLayer]=0.;
       
  }
}
//______________________________________________________________________
Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
  // compute truncated mean 

  Int_t nc=0;
  Double_t dedx[kMaxSamples];
  for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
    Double_t dedxsamp=GetdEdxSample(il);
    if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
      dedx[nc]= dedxsamp;
      nc++;
    }    
  }
  if(nc<1) return 0.;

  Int_t swap; // sort in ascending order
  do {
    swap=0;
    for (Int_t i=0; i<nc-1; i++) {
      if (dedx[i]<=dedx[i+1]) continue;
      Double_t tmp=dedx[i];
      dedx[i]=dedx[i+1]; 
      dedx[i+1]=tmp;
      swap++;
    }
  } while (swap);

  Double_t sumamp=0,sumweight=0;
  Double_t weight[kMaxSamples];
  for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.;
  Int_t lastUsed=(Int_t)(frac*nc+0.00001);
  if(lastUsed==0) lastUsed=1;
  if(lastUsed>kMaxSamples) lastUsed=kMaxSamples;
  for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
  if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5;
  for (Int_t i=0; i<nc; i++) {
    // AliDebug(5,Form("dE/dx %f   weight %f",dedx[i],weight[i]));
    sumamp+= dedx[i]*weight[i];
    sumweight+=weight[i];
  }
  if(sumweight>0.) return sumamp/sumweight;
  else return 0.;

}
//______________________________________________________________________
Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
  // compute generalized mean with k=-2 (used by CMS)
  Int_t nc=0;
  Double_t dedx[kMaxSamples];
  for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
    Double_t dedxsamp=GetdEdxSample(il);
    if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
      dedx[nc]= dedxsamp;
      nc++;      
    }
  }
  if(nc<1) return 0.;

  Double_t weiSum = 0.;
  for (Int_t i=0; i<nc; i++) {
    weiSum+=TMath::Power(dedx[i],-2);
  }
  Double_t wMean=0.;
  if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
  return wMean;

}
//______________________________________________________________________
void  AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
  // compute conditional probablilities
  const Int_t nPart = 3;
  Double_t itsProb[nPart] = {1,1,1}; // p, K, pi

  for(Int_t iS=0; iS<fNSamples; iS++){
    if(!HasPointOnLayer(iS)) continue;
    if(!UseLayerForPid(iS)) continue;
    Int_t iLayer=iS+3; // to match with present ITS
    if(iLayer>6) iLayer=6; // all extra points are treated as SSD
    Float_t dedx = GetdEdxSample(iS);
    if(dedx<mindedx) continue;
    Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
    itsProb[0] *= layProb;

    layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
    if (fP < 0.16) layProb=0.00001;
    itsProb[1] *= layProb;
    
    layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
    itsProb[2] *= layProb;
  }

  // Normalise probabilities
  Double_t sumProb = 0;
  for (Int_t iPart = 0; iPart < nPart; iPart++) {
    sumProb += itsProb[iPart];
  }
  sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions

  for (Int_t iPart = 0; iPart < nPart; iPart++) {
    itsProb[iPart]/=sumProb;
  }
  
  condprob[AliPID::kElectron] = itsProb[2];
  condprob[AliPID::kMuon] = itsProb[2];
  condprob[AliPID::kPion] = itsProb[2];
  condprob[AliPID::kKaon] = itsProb[1];
  condprob[AliPID::kProton] = itsProb[0];
  return;

}
//______________________________________________________________________
void  AliITSdEdxSamples::PrintAll() const{
  // print all the infos
  printf("Particle %d momentum %f GeV/c, number of points %d\n",
	 GetParticleSpecieMC(),
	 fP,
	 GetNumberOfEffectiveSamples());
  for(Int_t iLay=0; iLay<fNSamples; iLay++){
    printf("   Layer %d   Point %d   dE %f keV  dx %f cm  mom %f GeV/c\n",iLay,
	   HasPointOnLayer(iLay),
	   GetdESample(iLay),
	   GetdxSample(iLay),
	   GetMomentumAtSample(iLay));
  }

  printf("Layers used for PID:\n");
  printf("Layer ");
  for(Int_t iLay=0; iLay<fNSamples; iLay++){
    printf("%d ",iLay);
  }
  printf("\n");
  
  printf("Use   ");
  for(Int_t iLay=0; iLay<fNSamples; iLay++){
    printf("%d ",UseLayerForPid(iLay));
  }
  printf("\n");
  printf("Truncated mean = %f\n",GetTruncatedMean());
}
//______________________________________________________________________
void  AliITSdEdxSamples::PrintClusterMap() const{
  // print the cluster map

  printf("Layer ");
  for(Int_t iLay=0; iLay<fNSamples; iLay++){
    printf("%d ",iLay);
  }
  printf("\n");
  
  printf("Point ");
  for(Int_t iLay=0; iLay<fNSamples; iLay++){
    printf("%d ",HasPointOnLayer(iLay));
  }
  printf("\n");
}
Commit	Line	Data
9189f38b	1	/**************************************************************************
	2	* Copyright(c) 2009-2012, 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	/* $Id$ */
	16
	17
	18
	19	///////////////////////////////////////////////////////////////////
	20	// //
	21	// Class to store information for PID with ITS //
	22	// and truncated mean computation methods //
	23	// Origin: F.Prino, Torino, prino@to.infn.it //
	24	// //
	25	///////////////////////////////////////////////////////////////////
	26
	27	#include "AliITSPidParams.h"
	28	#include "AliITSdEdxSamples.h"
	29	#include "AliLog.h"
	30	#include <TMath.h>
	31
	32	ClassImp(AliITSdEdxSamples)
	33
	34	//______________________________________________________________________
	35	AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
	36	fNSamples(0),
0d66557a	37	fClusterMap(0),
9189f38b	38	fP(0.),
0d66557a	39	fParticleSpecie(0),
0d66557a	40	fLayersForPid(0xFFFF)
9189f38b	41	{
9189f38b	42	// Default constructor
92329a66	43	for(Int_t i=0; i<kMaxSamples; i++){
	44	fdESamples[i]=0.;
	45	fdxSamples[i]=0.;
cb0a52e0	46	fPAtSample[i]=0.;
cb0a52e0	47	}
9189f38b	48	}
	49
	50	//______________________________________________________________________
92329a66	51	AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) :
9189f38b	52	TObject(),
9189f38b	53	fNSamples(nSamples),
0d66557a	54	fClusterMap(0),
9189f38b	55	fP(mom),
0d66557a	56	fParticleSpecie(specie),
0d66557a	57	fLayersForPid(0xFFFF)
9189f38b	58	{
9189f38b	59	// Standard constructor
92329a66	60	SetdESamples(nSamples,esamples);
92329a66	61	SetdxSamples(nSamples,xsamples);
0d66557a	62	SetClusterMapFromdE();
92329a66	63	}
	64
	65	//______________________________________________________________________
	66	AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) :
	67	TObject(),
	68	fNSamples(source.fNSamples),
0d66557a	69	fClusterMap(source.fClusterMap),
92329a66	70	fP(source.fP),
0d66557a	71	fParticleSpecie(source.fParticleSpecie),
0d66557a	72	fLayersForPid(source.fLayersForPid)
92329a66	73	{
92329a66	74	// Copy constructor
8d91dc44	75	for(Int_t i=0; i<kMaxSamples; i++){
92329a66	76	fdESamples[i]=source.GetdESample(i);
	77	fdxSamples[i]=source.GetdxSample(i);
	78	fPAtSample[i]=source.GetMomentumAtSample(i);
	79	}
9189f38b	80	}
9189f38b	81	//______________________________________________________________________
92329a66	82	void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){
9189f38b	83	// Set the samples
9189f38b	84
92329a66	85	if(nSamples>kMaxSamples){
	86	AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples));
	87	fNSamples=kMaxSamples;
9189f38b	88	}else{
	89	fNSamples=nSamples;
	90	}
92329a66	91	for(Int_t i=0; i<fNSamples; i++) fdESamples[i]=samples[i];
92329a66	92	for(Int_t i=fNSamples; i<kMaxSamples; i++) fdESamples[i]=0.;
9189f38b	93	return;
	94	}
	95	//______________________________________________________________________
92329a66	96	void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){
9189f38b	97	// Set the samples
92329a66	98
	99	if(nSamples>kMaxSamples){
	100	AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples))
	101	fNSamples=kMaxSamples;
	102	}else{
	103	fNSamples=nSamples;
	104	}
	105	for(Int_t i=0; i<fNSamples; i++) fdxSamples[i]=samples[i];
	106	for(Int_t i=fNSamples; i<kMaxSamples; i++) fdxSamples[i]=0.;
	107	return;
	108	}
	109
	110	//______________________________________________________________________
	111	void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){
	112	// Set the samples
	113	SetdESamples(nSamples,esamples);
	114	SetdxSamples(nSamples,xsamples);
9189f38b	115	for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
92329a66	116	for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.;
9189f38b	117	return;
	118	}
	119	//______________________________________________________________________
0d66557a	120	void AliITSdEdxSamples::SetLayerSample(Int_t iLayer, Bool_t haspoint, Double_t dE, Double_t dx, Double_t p){
	121	// set info from single layer
	122	if(haspoint){
	123	SetPointOnLayer(iLayer);
	124	fdESamples[iLayer]=dE;
	125	fdxSamples[iLayer]=dx;
	126	fPAtSample[iLayer]=p;
	127	}else{
	128	if(HasPointOnLayer(iLayer)) fClusterMap-=(1<<iLayer);
	129	fdESamples[iLayer]=0.;
	130	fdxSamples[iLayer]=0.;
	131	fPAtSample[iLayer]=0.;
	132
	133	}
	134	}
	135	//______________________________________________________________________
9189f38b	136	Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
	137	// compute truncated mean
	138
	139	Int_t nc=0;
92329a66	140	Double_t dedx[kMaxSamples];
9189f38b	141	for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
92329a66	142	Double_t dedxsamp=GetdEdxSample(il);
0d66557a	143	if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
92329a66	144	dedx[nc]= dedxsamp;
9189f38b	145	nc++;
92329a66	146	}
9189f38b	147	}
9189f38b	148	if(nc<1) return 0.;
0d66557a	149
9189f38b	150	Int_t swap; // sort in ascending order
	151	do {
	152	swap=0;
	153	for (Int_t i=0; i<nc-1; i++) {
	154	if (dedx[i]<=dedx[i+1]) continue;
	155	Double_t tmp=dedx[i];
	156	dedx[i]=dedx[i+1];
	157	dedx[i+1]=tmp;
	158	swap++;
	159	}
	160	} while (swap);
	161
	162	Double_t sumamp=0,sumweight=0;
92329a66	163	Double_t weight[kMaxSamples];
92329a66	164	for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.;
9189f38b	165	Int_t lastUsed=(Int_t)(frac*nc+0.00001);
9189f38b	166	if(lastUsed==0) lastUsed=1;
92329a66	167	if(lastUsed>kMaxSamples) lastUsed=kMaxSamples;
9189f38b	168	for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
92329a66	169	if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5;
9189f38b	170	for (Int_t i=0; i<nc; i++) {
92329a66	171	// AliDebug(5,Form("dE/dx %f weight %f",dedx[i],weight[i]));
9189f38b	172	sumamp+= dedx[i]*weight[i];
	173	sumweight+=weight[i];
	174	}
	175	if(sumweight>0.) return sumamp/sumweight;
	176	else return 0.;
	177
	178	}
	179	//______________________________________________________________________
	180	Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
	181	// compute generalized mean with k=-2 (used by CMS)
9189f38b	182	Int_t nc=0;
92329a66	183	Double_t dedx[kMaxSamples];
9189f38b	184	for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
92329a66	185	Double_t dedxsamp=GetdEdxSample(il);
0d66557a	186	if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
92329a66	187	dedx[nc]= dedxsamp;
92329a66	188	nc++;
9189f38b	189	}
	190	}
	191	if(nc<1) return 0.;
	192
	193	Double_t weiSum = 0.;
	194	for (Int_t i=0; i<nc; i++) {
	195	weiSum+=TMath::Power(dedx[i],-2);
	196	}
	197	Double_t wMean=0.;
	198	if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
	199	return wMean;
	200
	201	}
	202	//______________________________________________________________________
	203	void AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
	204	// compute conditional probablilities
	205	const Int_t nPart = 3;
	206	Double_t itsProb[nPart] = {1,1,1}; // p, K, pi
	207
	208	for(Int_t iS=0; iS<fNSamples; iS++){
0d66557a	209	if(!HasPointOnLayer(iS)) continue;
0d66557a	210	if(!UseLayerForPid(iS)) continue;
9189f38b	211	Int_t iLayer=iS+3; // to match with present ITS
9189f38b	212	if(iLayer>6) iLayer=6; // all extra points are treated as SSD
92329a66	213	Float_t dedx = GetdEdxSample(iS);
92329a66	214	if(dedx<mindedx) continue;
9189f38b	215	Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
	216	itsProb[0] *= layProb;
	217
	218	layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
	219	if (fP < 0.16) layProb=0.00001;
	220	itsProb[1] *= layProb;
	221
	222	layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
	223	itsProb[2] *= layProb;
	224	}
	225
	226	// Normalise probabilities
	227	Double_t sumProb = 0;
	228	for (Int_t iPart = 0; iPart < nPart; iPart++) {
	229	sumProb += itsProb[iPart];
	230	}
	231	sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
	232
	233	for (Int_t iPart = 0; iPart < nPart; iPart++) {
	234	itsProb[iPart]/=sumProb;
	235	}
	236
	237	condprob[AliPID::kElectron] = itsProb[2];
	238	condprob[AliPID::kMuon] = itsProb[2];
	239	condprob[AliPID::kPion] = itsProb[2];
	240	condprob[AliPID::kKaon] = itsProb[1];
	241	condprob[AliPID::kProton] = itsProb[0];
	242	return;
	243
	244	}
0d66557a	245	//______________________________________________________________________
	246	void AliITSdEdxSamples::PrintAll() const{
	247	// print all the infos
	248	printf("Particle %d momentum %f GeV/c, number of points %d\n",
	249	GetParticleSpecieMC(),
	250	fP,
	251	GetNumberOfEffectiveSamples());
	252	for(Int_t iLay=0; iLay<fNSamples; iLay++){
	253	printf(" Layer %d Point %d dE %f keV dx %f cm mom %f GeV/c\n",iLay,
	254	HasPointOnLayer(iLay),
	255	GetdESample(iLay),
	256	GetdxSample(iLay),
	257	GetMomentumAtSample(iLay));
	258	}
	259
	260	printf("Layers used for PID:\n");
	261	printf("Layer ");
	262	for(Int_t iLay=0; iLay<fNSamples; iLay++){
	263	printf("%d ",iLay);
	264	}
	265	printf("\n");
	266
	267	printf("Use ");
	268	for(Int_t iLay=0; iLay<fNSamples; iLay++){
	269	printf("%d ",UseLayerForPid(iLay));
	270	}
	271	printf("\n");
	272	printf("Truncated mean = %f\n",GetTruncatedMean());
	273	}
	274	//______________________________________________________________________
	275	void AliITSdEdxSamples::PrintClusterMap() const{
	276	// print the cluster map
	277
	278	printf("Layer ");
	279	for(Int_t iLay=0; iLay<fNSamples; iLay++){
	280	printf("%d ",iLay);
	281	}
	282	printf("\n");
	283
	284	printf("Point ");
	285	for(Int_t iLay=0; iLay<fNSamples; iLay++){
	286	printf("%d ",HasPointOnLayer(iLay));
	287	}
	288	printf("\n");
	289	}