PWGHF/correlationHF/AliDxHFEParticleSelectionEl.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
   9 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
  10 //*                                                                        *
  11 //* Permission to use, copy, modify and distribute this software and its   *
  12 //* documentation strictly for non-commercial purposes is hereby granted   *
  13 //* without fee, provided that the above copyright notice appears in all   *
  14 //* copies and that both the copyright notice and this permission notice   *
  15 //* appear in the supporting documentation. The authors make no claims     *
  16 //* about the suitability of this software for any purpose. It is          *
  17 //* provided "as is" without express or implied warranty.                  *
  18 //**************************************************************************
  19
  20 /// @file   AliDxHFEParticleSelectionEl.cxx
  21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
  22 /// @date   2012-03-19
  23 /// @brief  D0 selection for D0-HFE correlation
  24 ///
  25 #include "AliDxHFEParticleSelectionEl.h"
  26 #include "AliVParticle.h"
  27 #include "AliVEvent.h"
  28 #include "AliPID.h"
  29 #include "AliPIDResponse.h"
  30 #include "AliHFEcontainer.h"
  31 #include "AliHFEpid.h"
  32 #include "AliHFEpidBase.h"
  33 #include "AliHFEtools.h"
  34 #include "AliHFEcuts.h"
  35 #include "AliAODTrack.h"
  36 #include "AliAnalysisDataSlot.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliCFManager.h"
  40 #include "THnSparse.h"
  41 #include "TH1F.h"
  42 #include "TH2F.h"
  43 #include "TAxis.h"
  44 #include "TList.h"
  45 #include "TObjArray.h"
  46 #include <iostream>
  47 #include <cerrno>
  48 #include <memory>
  49
  50 using namespace std;
  51
  52 /// ROOT macro for the implementation of ROOT specific class methods
  53 ClassImp(AliDxHFEParticleSelectionEl)
  54
  55 AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
  56   : AliDxHFEParticleSelection("electron", opt)
  57   , fPID(NULL)
  58   , fPIDTOF(NULL)
  59   , fElectronProperties(NULL)
  60   , fHistoList(NULL)
  61   , fCutPidList(NULL)
  62   , fPIDResponse(NULL)
  63   , fCuts(NULL)
  64   , fCFM(NULL)
  65 {
  66   // constructor
  67   //
  68   //
  69   //
  70   //
  71 }
  72
  73 AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
  74 {
  75   // destructor
  76   if (fElectronProperties) {
  77     delete fElectronProperties;
  78     fElectronProperties=NULL;
  79   }
  80   if(fHistoList){
  81     delete fHistoList;
  82     fHistoList=NULL;
  83   }
  84   if(fCFM){
  85     delete fCFM;
  86     fCFM=NULL;
  87   }
  88   if(fPIDResponse){
  89     delete fPIDResponse;
  90     fPIDResponse=NULL;
  91   }
  92   // NOTE: external objects fPID and fCuts are not deleted here
  93   fPID=NULL;
  94   fPIDTOF=NULL;
  95   fCuts=NULL;
  96 }
  97
  98 const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
  99   "kRecKineITSTPC",
 100   "kRecPrim",
 101   "kHFEcuts",
 102   "kHFEcutsTOFPID",
 103   "kHFEcutsTPCPID",
 104   "kPIDTOF",
 105   "kPIDTPCTOF",
 106   "Selected e"
 107
 108 };
 109
 110 int AliDxHFEParticleSelectionEl::Init()
 111 {
 112   //
 113   // init function
 114   //
 115   int iResult=0;
 116
 117   // Implicit call to InitControlObjects() before setting up CFM and fCuts
 118   // (if not there)
 119   iResult=AliDxHFEParticleSelection::Init();
 120   if (iResult<0) return iResult;
 121
 122   //--------Initialize correction Framework and Cuts-------------------------
 123   // Consider moving this, either to separate function or
 124   // add a set function for AliCFManager
 125   // Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
 126   AliInfo("Setting up CFM");
 127   fCFM = new AliCFManager;
 128   // the setup of cut objects is done in AliHFEcuts::Initialize
 129   // the ids used by this class must be the same, the code might be broken if
 130   // the sequence in AliHFEcuts::Initialize is changed
 131   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 132   // reset pointers in the CF manager
 133   fCFM->SetNStepParticle(kNcutSteps);
 134   for(Int_t istep = 0; istep < kNcutSteps; istep++) {
 135     fCFM->SetParticleCutsList(istep, NULL);
 136   }
 137   if(!fCuts) {
 138     AliWarning("Cuts not available. Default cuts will be used");
 139     fCuts = new AliHFEcuts;
 140     fCuts->CreateStandardCuts();
 141   }
 142   // TODO: error handling?
 143   fCuts->Initialize(fCFM);
 144
 145   return 0;
 146
 147 }
 148
 149 THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
 150 {
 151   //
 152   // Defines the THnSparse. For now, only calls CreatControlTHnSparse
 153
 154   // here is the only place to change the dimension
 155   const int thnSize = 3;
 156   InitTHnSparseArray(thnSize);
 157   const double Pi=TMath::Pi();
 158   TString name;
 159   name.Form("%s info", GetName());
 160
 161   //                                 0    1       2
 162   //                                 Pt   Phi    Eta
 163   int         thnBins [thnSize] = { 1000,  200, 500};
 164   double      thnMin  [thnSize] = {    0,    0, -1.};
 165   double      thnMax  [thnSize] = {  100, 2*Pi,  1.};
 166   const char* thnNames[thnSize] = { "Pt","Phi","Eta"};
 167
 168   return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
 169 }
 170
 171 int AliDxHFEParticleSelectionEl::InitControlObjects()
 172 {
 173   /// init control and monitoring objects
 174   AliInfo("Electron THnSparse");
 175
 176   fElectronProperties=DefineTHnSparse();
 177   AddControlObject(fElectronProperties);
 178
 179   double dEdxBins[6]={1000,0.,10.,200,0.,200.};
 180   double nSigBins[6]={1000,0.,10.,200,-10.,10.};
 181
 182   fHistoList=new TList;
 183   fHistoList->SetName("HFE Histograms");
 184   fHistoList->SetOwner();
 185
 186   // Histogram storing which cuts have been applied to the tracks
 187   fHistoList->Add(CreateControlHistogram("fWhichCut","effective cut for a rejected particle", kNCutLabels, fgkCutBinNames));
 188
 189   // dEdx plots, TPC signal vs momentum
 190   fHistoList->Add(CreateControl2DHistogram("fdEdx", "dEdx before cuts", dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 191   fHistoList->Add(CreateControl2DHistogram("fdEdxCut", "dEdx after cuts",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 192   fHistoList->Add(CreateControl2DHistogram("fdEdxPidTOF", "dEdx after TOF pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 193   fHistoList->Add(CreateControl2DHistogram("fdEdxPid", "dEdx after pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 194
 195   // nSigmaTPC vs momentum
 196   fHistoList->Add(CreateControl2DHistogram("fnSigTPC", "nSigTPC before cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 197   fHistoList->Add(CreateControl2DHistogram("fnSigTPCCut", "nSigmaTPC after cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 198   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPidTOF", "nSigmaTPC after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 199   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPid", "nSigmaTPC after PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 200
 201   // nSigmaTOF vs momentum
 202   fHistoList->Add(CreateControl2DHistogram("fnSigTOF", "nSigmaTOF before cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 203   fHistoList->Add(CreateControl2DHistogram("fnSigTOFCut", "nSigmaTOF after cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 204   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPidTOF", "nSigmaTOF after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 205   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPid", "nSigmaTOF after PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 206
 207   AddControlObject(fHistoList);
 208
 209   return AliDxHFEParticleSelection::InitControlObjects();
 210 }
 211
 212 int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
 213 {
 214   /// histogram particle properties
 215   if (!p) return -EINVAL;
 216   //if (!fControlObjects) return 0;
 217   if(selectionCode==0) return  0;
 218
 219   if(fElectronProperties && ParticleProperties()) {
 220     FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
 221     fElectronProperties->Fill(ParticleProperties());
 222   }
 223
 224   return 0;
 225 }
 226
 227 int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 228 {
 229   // fill the data array from the particle data
 230   if (!data) return -EINVAL;
 231   AliAODTrack *track=(AliAODTrack*)p;
 232   if (!track) return -ENODATA;
 233   int i=0;
 234   if (dimension!=GetDimTHnSparse()) {
 235     // TODO: think about filling only the available data and throwing a warning
 236     return -ENOSPC;
 237   }
 238   data[i++]=track->Pt();
 239   data[i++]=track->Phi();
 240   data[i++]=track->Eta();
 241   return i;
 242 }
 243
 244
 245 int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
 246 {
 247   /// select El candidates
 248   // TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected.
 249
 250   AliAODTrack *track=(AliAODTrack*)pEl;
 251   fCFM->SetRecEventInfo(pEvent);
 252
 253   ((TH2D*)fHistoList->FindObject("fdEdx"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 254   ((TH2D*)fHistoList->FindObject("fnSigTPC"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 255   ((TH2D*)fHistoList->FindObject("fnSigTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 256
 257   //--------track cut selection-----------------------
 258   //Using AliHFECuts:
 259   // RecKine: ITSTPC cuts
 260   if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
 261     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecKineITSTPC);
 262     return 0;
 263   }
 264
 265   // RecPrim
 266   if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
 267     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecPrim);
 268     return 0;
 269   }
 270
 271   // HFEcuts: ITS layers cuts
 272   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
 273     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsITS);
 274     return 0;
 275   }
 276
 277   // HFE cuts: TOF PID and mismatch flag
 278   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
 279     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTOF);
 280     return 0;
 281   }
 282
 283   // HFE cuts: TPC PID cleanup
 284   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
 285     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTPC);
 286     return 0;
 287   }
 288
 289   // HFEcuts: Nb of tracklets TRD0
 290   //if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 291   ((TH2D*)fHistoList->FindObject("fdEdxCut"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 292   ((TH2D*)fHistoList->FindObject("fnSigTPCCut"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 293   ((TH2D*)fHistoList->FindObject("fnSigTOFCut"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 294
 295
 296   //--------PID selection-----------------------
 297   AliHFEpidObject hfetrack;
 298   hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
 299   hfetrack.SetRecTrack(track);
 300
 301   // TODO: configurable colliding system
 302   //if(IsPbPb()) hfetrack.SetPbPb();
 303   hfetrack.SetPP();
 304
 305   // TODO: Put this into a while-loop instead, looping over the number of pid objects in the cut-list?
 306   // This needs a bit of thinking and finetuning (wrt histogramming)
 307   if(fPIDTOF && fPIDTOF->IsSelected(&hfetrack)) {
 308     ((TH2D*)fHistoList->FindObject("fdEdxPidTOF"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 309     ((TH2D*)fHistoList->FindObject("fnSigTPCPidTOF"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 310     ((TH2D*)fHistoList->FindObject("fnSigTOFPidTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 311   }
 312   else{
 313     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOF);
 314     return 0;
 315   }
 316
 317   //Combined tof & tpc pid
 318   if(fPID && fPID->IsSelected(&hfetrack)) {
 319     AliDebug(3,"Inside FilldPhi, electron is selected");
 320     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
 321     ((TH2D*)fHistoList->FindObject("fdEdxPid"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 322     ((TH2D*)fHistoList->FindObject("fnSigTPCPid"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 323     ((TH2D*)fHistoList->FindObject("fnSigTOFPid"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 324     return 1;
 325   }
 326   else{
 327     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPID);
 328     return 0;
 329   }
 330
 331
 332 }
 333
 334 void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
 335 {
 336   /// set cut objects
 337   if (level==kCutHFE) {
 338     fCuts=dynamic_cast<AliHFEcuts*>(cuts);
 339     if (!fCuts && cuts) {
 340       AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
 341     }
 342     return;
 343   }
 344
 345   if (level==kCutPID) {
 346     fPID=dynamic_cast<AliHFEpid*>(cuts);
 347     if (!fPID && cuts) {
 348       AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 349     }
 350     return;
 351   }
 352
 353   if (level==kCutPIDTOF) {
 354     fPIDTOF=dynamic_cast<AliHFEpid*>(cuts);
 355     if (!fPIDTOF && cuts) {
 356       AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 357     }
 358     return;
 359   }
 360   if(level==kCutList){
 361     fCutPidList=dynamic_cast<TList*>(cuts);
 362     if (!fCutPidList && cuts) {
 363       AliError(Form("cuts object is not of required type TList but %s", cuts->ClassName()));
 364     }
 365     else{
 366       // TODO: Could be done more elegantly, at the moment requires that the cut and pid objects are in
 367       // a specific order..
 368       TObject *obj=NULL;
 369       int iii=0;
 370       TIter next(fCutPidList);
 371       while((obj = next())){
 372         iii++;
 373         if(iii==1) {
 374           fCuts=dynamic_cast<AliHFEcuts*>(obj);
 375           if (!fCuts)
 376             AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
 377         }
 378         if(iii==2){
 379           fPID=dynamic_cast<AliHFEpid*>(obj);
 380           if (!fPID)
 381             AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 382         }
 383         if(iii==3){
 384           fPIDTOF=dynamic_cast<AliHFEpid*>(obj);
 385           if (!fPIDTOF)
 386             AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 387         }
 388       }
 389
 390     }
 391     return;
 392   }
 393 }
 394
 395 //________________________________________________________________________
 396 Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 397 {
 398   // Check single track cuts for a given cut step
 399   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 400   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 401   //if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 402   return kTRUE;
 403 }