Fixing a little overlap. Adding flags to perform material budget studies (Mario)
[u/mrichter/AliRoot.git] / ITS / AliITSdEdxSamples.cxx
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 }