1 /**************************************************************************
2 * Copyright(c) 2009-2012, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 ///////////////////////////////////////////////////////////////////
21 // Class to store information for PID with ITS //
22 // and truncated mean computation methods //
23 // Origin: F.Prino, Torino, prino@to.infn.it //
25 ///////////////////////////////////////////////////////////////////
27 #include "AliITSPidParams.h"
28 #include "AliITSdEdxSamples.h"
32 ClassImp(AliITSdEdxSamples)
34 //______________________________________________________________________
35 AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
42 // Default constructor
43 for(Int_t i=0; i<kMaxSamples; i++){
50 //______________________________________________________________________
51 AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) :
56 fParticleSpecie(specie),
59 // Standard constructor
60 SetdESamples(nSamples,esamples);
61 SetdxSamples(nSamples,xsamples);
62 SetClusterMapFromdE();
65 //______________________________________________________________________
66 AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) :
68 fNSamples(source.fNSamples),
69 fClusterMap(source.fClusterMap),
71 fParticleSpecie(source.fParticleSpecie),
72 fLayersForPid(source.fLayersForPid)
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);
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;
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);
99 //______________________________________________________________________
100 void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){
103 if(nSamples>kMaxSamples){
104 AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples));
105 fNSamples=kMaxSamples;
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.;
113 //______________________________________________________________________
114 void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){
117 if(nSamples>kMaxSamples){
118 AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples));
119 fNSamples=kMaxSamples;
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.;
128 //______________________________________________________________________
129 void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){
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.;
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
141 SetPointOnLayer(iLayer);
142 fdESamples[iLayer]=dE;
143 fdxSamples[iLayer]=dx;
144 fPAtSample[iLayer]=p;
146 if(HasPointOnLayer(iLayer)) fClusterMap-=(1<<iLayer);
147 fdESamples[iLayer]=0.;
148 fdxSamples[iLayer]=0.;
149 fPAtSample[iLayer]=0.;
153 //______________________________________________________________________
154 Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
155 // compute truncated mean
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){
168 Int_t swap; // sort in ascending order
171 for (Int_t i=0; i<nc-1; i++) {
172 if (dedx[i]<=dedx[i+1]) continue;
173 Double_t tmp=dedx[i];
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];
193 if(sumweight>0.) return sumamp/sumweight;
197 //______________________________________________________________________
198 Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
199 // compute generalized mean with k=-2 (used by CMS)
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){
211 Double_t weiSum = 0.;
212 for (Int_t i=0; i<nc; i++) {
213 weiSum+=TMath::Power(dedx[i],-2);
216 if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
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
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;
236 layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
237 if (fP < 0.16) layProb=0.00001;
238 itsProb[1] *= layProb;
240 layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
241 itsProb[2] *= layProb;
244 // Normalise probabilities
245 Double_t sumProb = 0;
246 for (Int_t iPart = 0; iPart < nPart; iPart++) {
247 sumProb += itsProb[iPart];
249 sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
251 for (Int_t iPart = 0; iPart < nPart; iPart++) {
252 itsProb[iPart]/=sumProb;
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];
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(),
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),
275 GetMomentumAtSample(iLay));
278 printf("Layers used for PID:\n");
280 for(Int_t iLay=0; iLay<fNSamples; iLay++){
286 for(Int_t iLay=0; iLay<fNSamples; iLay++){
287 printf("%d ",UseLayerForPid(iLay));
290 printf("Truncated mean = %f\n",GetTruncatedMean());
292 //______________________________________________________________________
293 void AliITSdEdxSamples::PrintClusterMap() const{
294 // print the cluster map
297 for(Int_t iLay=0; iLay<fNSamples; iLay++){
303 for(Int_t iLay=0; iLay<fNSamples; iLay++){
304 printf("%d ",HasPointOnLayer(iLay));