PWG3/hfe/AliHFEpidITS.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //
  16 // ITS PID class
  17 // checks ITS PID based on ITS dE/dx truncated mean
  18 //
  19 // Authors: Matus Kalisky <matus.kalisky@cern.ch>
  20 //          Markus Fasel <M.Fasel@gsi.de>
  21 //
  22 #include <TClass.h>
  23 #include <TH2F.h>
  24 #include <TList.h>
  25 #include <TMath.h>
  26 #include <TString.h>
  27
  28 //#include "AliAODTrack.h"
  29 //#include "AliAODMCParticle.h"
  30 #include "AliESDtrack.h"
  31 #include "AliLog.h"
  32 #include "AliPID.h"
  33 #include "AliVParticle.h"
  34
  35 #include "AliHFEpidITS.h"
  36
  37 //___________________________________________________________________
  38 AliHFEpidITS::AliHFEpidITS(const Char_t *name):
  39     AliHFEpidBase(name)
  40   , fQAlist(0x0)
  41 {
  42   //
  43   // Default constructor
  44   //
  45 }
  46
  47 //___________________________________________________________________
  48 AliHFEpidITS::AliHFEpidITS(const AliHFEpidITS &ref):
  49     AliHFEpidBase("")
  50   , fQAlist(0x0)
  51 {
  52   //
  53   // Copy constructor
  54   //
  55   ref.Copy(*this);
  56 }
  57
  58 //___________________________________________________________________
  59 AliHFEpidITS &AliHFEpidITS::operator=(const AliHFEpidITS &ref){
  60   //
  61   // Assignment operator
  62   //
  63   if(this != &ref) ref.Copy(*this);
  64   return *this;
  65 }
  66
  67 //___________________________________________________________________
  68 AliHFEpidITS::~AliHFEpidITS(){
  69   //
  70   // Destructor
  71   //
  72   if(fQAlist){
  73     fQAlist->Clear();
  74     delete fQAlist;
  75   }
  76 }
  77
  78 //___________________________________________________________________
  79 void AliHFEpidITS::Copy(TObject &o) const {
  80   //
  81   // Copy function
  82   // Provides a deep copy
  83   //
  84   AliHFEpidITS &target = dynamic_cast<AliHFEpidITS &>(o);
  85
  86   target.fQAlist = dynamic_cast<TList *>(fQAlist->Clone());
  87   AliHFEpidBase::Copy(target);
  88 }
  89
  90 //___________________________________________________________________
  91 Bool_t AliHFEpidITS::InitializePID(){
  92   //
  93   // ITS PID initialization
  94   //
  95   return kTRUE;
  96 }
  97
  98
  99 //___________________________________________________________________
 100 Int_t AliHFEpidITS::IsSelected(AliHFEpidObject* /*track*/){
 101   //
 102   // Does PID decision for ITS
 103   //
 104   return 11;  // @TODO: Implement ITS PID decision
 105 }
 106
 107 //___________________________________________________________________
 108 void AliHFEpidITS::AddQAhistograms(TList *l){
 109   //
 110   // Adding QA histograms for ITS PID
 111   //
 112   // QA histograms are:
 113   // - all Particles ITS signal vs. p (both methods)
 114   // - single particle ITS signal vs. p (both methods)
 115   //
 116   if(!fQAlist) fQAlist = new TList;
 117   fQAlist->SetName("fITSqaHistograms");
 118
 119   // prepare axis
 120   const Int_t kMomentumBins = 41;
 121   const Double_t kPtMin = 0.1;
 122   const Double_t kPtMax = 10.;
 123   const Int_t kSigBins = 300;
 124   Double_t momentumBins[kMomentumBins];
 125   for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
 126     momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
 127
 128   TH2 *histo = NULL;
 129   fQAlist->AddAt((histo = new TH2F("fITSsigV1all", "ITS signal vs. p (all species, Method1):p / GeV/c:ITS signal / a.u.", kMomentumBins - 1, momentumBins, kSigBins, 0., kSigBins)), kITSsigV1);
 130   fQAlist->AddAt((histo = new TH2F("fITSsigV2all", "ITS signal vs. p (all species, Method2):p / GeV/c:ITS signal / a.u.", kMomentumBins - 1, momentumBins, kSigBins, 0., kSigBins)), kITSsigV2);
 131   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 132     fQAlist->AddAt((histo = new TH2F(Form("fITSsigV1%s", AliPID::ParticleName(ispec)), Form("ITS signal vs. p (%s, Method1):p / GeV/c:ITS signal / a.u.", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBins, 0., kSigBins)), 2 * ispec + kHistosSigAll);
 133     fQAlist->AddAt((histo = new TH2F(Form("fITSsigV2%s", AliPID::ParticleName(ispec)), Form("ITS signal vs. p (%s, Method2):p / GeV/c:ITS signal / a.u.", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBins, 0., kSigBins)), 2 * ispec + 1 + kHistosSigAll);
 134   }
 135   l->Add(fQAlist);
 136 }
 137
 138 //___________________________________________________________________
 139 Double_t AliHFEpidITS::GetITSSignalV1(AliVParticle *vtrack, Int_t mcPID){
 140   //
 141   // Calculate the ITS signal according to the mean charge of the clusters
 142   //
 143   if(!TString(vtrack->IsA()->GetName()).CompareTo("AliAODtrack")){
 144     AliError("PID for AODs not implemented yet");
 145     return 0.;
 146   }
 147   AliESDtrack *track = dynamic_cast<AliESDtrack *>(vtrack);
 148   Double_t signal = 0.;
 149 #ifdef TRUNK
 150   Double_t dedx[4];
 151   track->GetITSdEdxSamples(dedx);
 152   signal = TMath::Mean(4, dedx);
 153 #else
 154   signal = track->GetITSsignal();
 155 #endif
 156   Double_t p = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : track->P();
 157   AliDebug(1, Form("Momentum: %f, ITS Signal: %f", p, signal));
 158   if(IsQAon()) FillHistogramsSignalV1(p, signal, mcPID);
 159   return signal;
 160 }
 161
 162 //___________________________________________________________________
 163 Double_t AliHFEpidITS::GetITSSignalV2(AliVParticle *vtrack, Int_t mcPID){
 164   //
 165   // Calculates the ITS signal. Truncated mean is used.
 166   //
 167   if(!TString(vtrack->IsA()->GetName()).CompareTo("AliAODtrack")){
 168     AliError("PID for AODs not implemented yet");
 169     return 0.;
 170   }
 171   AliESDtrack *track = dynamic_cast<AliESDtrack *>(vtrack);
 172   Double_t dedx[4], tmp[4];
 173   Int_t indices[4];
 174   track->GetITSdEdxSamples(tmp);
 175   TMath::Sort(4, tmp, indices);
 176   for(Int_t ien = 0; ien < 4; ien++) dedx[ien] = tmp[indices[ien]];
 177   Double_t signal = TMath::Mean(3, dedx);
 178   Double_t p = track->GetTPCInnerParam() ? track->GetTPCInnerParam()->P() : track->P();
 179   AliDebug(1, Form("Momentum: %f, ITS Signal: %f", p, signal));
 180   if(IsQAon()) FillHistogramsSignalV2(p, signal, mcPID);
 181   return signal;
 182 }
 183
 184 //___________________________________________________________________
 185 void AliHFEpidITS::FillHistogramsSignalV1(Double_t p, Double_t signal, Int_t species){
 186   (dynamic_cast<TH2 *>(fQAlist->At(kITSsigV1)))->Fill(p, signal);
 187   if(species >= 0 && species < AliPID::kSPECIES)
 188     (dynamic_cast<TH2 *>(fQAlist->At(kHistosSigAll + 2 * species)))->Fill(p, signal);
 189 }
 190
 191 //___________________________________________________________________
 192 void AliHFEpidITS::FillHistogramsSignalV2(Double_t p, Double_t signal, Int_t species){
 193   (dynamic_cast<TH2 *>(fQAlist->At(kITSsigV2)))->Fill(p, signal);
 194   if(species >= 0 && species < AliPID::kSPECIES)
 195     (dynamic_cast<TH2 *>(fQAlist->At(kHistosSigAll + 2 * species + 1)))->Fill(p, signal);
 196 }
 197