[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);
  }
}
//_____________________________________________________________________________
AliITSdEdxSamples& AliITSdEdxSamples::operator=(const AliITSdEdxSamples &source){
  // Assignment operator
 if(this==&source) return *this;
  ((TObject *)this)->operator=(source);
  fNSamples = source.fNSamples;
  fClusterMap = source.fClusterMap;
  fP = source.fP;
  fParticleSpecie = source.fParticleSpecie;
  fLayersForPid = source.fLayersForPid;
  for(Int_t i=0; i<kMaxSamples; i++){
    fdESamples[i]=source.GetdESample(i);
    fdxSamples[i]=source.GetdxSample(i);
    fPAtSample[i]=source.GetMomentumAtSample(i);
  }
  return *this;
}

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