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