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 }