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 "AliSelectNonHFE.h"
  27 #include "AliReducedParticle.h"
  28 #include "AliESDtrackCuts.h"
  29 #include "AliVParticle.h"
  30 #include "AliVEvent.h"
  31 #include "AliPID.h"
  32 #include "AliPIDResponse.h"
  33 #include "AliHFEcontainer.h"
  34 #include "AliHFEpid.h"
  35 #include "AliHFEpidBase.h"
  36 #include "AliHFEtools.h"
  37 #include "AliHFEcuts.h"
  38 #include "AliAODTrack.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAnalysisDataSlot.h"
  41 #include "AliAnalysisDataContainer.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliExternalTrackParam.h"
  44 #include "AliAODVertex.h"
  45 #include "AliHFEextraCuts.h"
  46 #include "AliCFManager.h"
  47 #include "THnSparse.h"
  48 #include "AliLog.h"
  49 #include "TH1F.h"
  50 #include "TH2F.h"
  51 #include "TAxis.h"
  52 #include "TList.h"
  53 #include "TObjArray.h"
  54 #include "AliKFParticle.h"
  55 #include "AliKFVertex.h"
  56 #include <iostream>
  57 #include <cerrno>
  58 #include <memory>
  59
  60 using namespace std;
  61
  62 /// ROOT macro for the implementation of ROOT specific class methods
  63 ClassImp(AliDxHFEParticleSelectionEl)
  64
  65 AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
  66   : AliDxHFEParticleSelection("electron", opt)
  67   , fPIDTOFTPC(NULL)
  68   , fPIDTOF(NULL)
  69   , fPIDTPC(NULL)
  70   , fElectronProperties(NULL)
  71   , fHistoList(NULL)
  72   , fCutPidList(NULL)
  73   , fPIDResponse(NULL)
  74   , fCuts(NULL)
  75   , fSelNHFE(NULL)
  76   , fTrackCuts(NULL)
  77   , fCFM(NULL)
  78   , fFinalCutStep(kPIDTOFTPC)
  79   , fInvMassLow(0.1)
  80   , fUseInvMassCut(kNoInvMass)
  81   , fSystem(0)
  82   , fTrackNum(0)
  83   , fImpactParamCutRadial(3)
  84   , fEtaCut(0.8)
  85   , fSurvivedCutStep(kNotSelected)
  86   , fStoreCutStepInfo(kFALSE)
  87
  88 {
  89   // constructor
  90   //
  91   //
  92   //
  93   //
  94   ParseArguments(opt);
  95 }
  96
  97 AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
  98 {
  99   // destructor
 100   if (fElectronProperties) {
 101     delete fElectronProperties;
 102     fElectronProperties=NULL;
 103   }
 104   if(fHistoList){
 105     delete fHistoList;
 106     fHistoList=NULL;
 107   }
 108   if(fCFM){
 109     delete fCFM;
 110     fCFM=NULL;
 111   }
 112   if(fPIDResponse){
 113     delete fPIDResponse;
 114     fPIDResponse=NULL;
 115   }
 116   // NOTE: external objects fPID and fCuts are not deleted here
 117   fPIDTOFTPC=NULL;
 118   fPIDTOF=NULL;
 119   fPIDTPC=NULL;
 120   fCuts=NULL;
 121   if(fSelNHFE){
 122     delete fSelNHFE;
 123     fSelNHFE=NULL;
 124   }
 125   if(fTrackCuts){
 126     delete fTrackCuts;
 127     fTrackCuts=NULL;
 128   }
 129 }
 130
 131 const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
 132   "kRecKineITSTPC",
 133   "kRecPrim",
 134   "kHFEcuts",
 135   "kHFEcutsTOFPID",
 136   "kHFEcutsTPCPID",
 137   "kPIDTOF",
 138   "kPIDTPC",
 139   "kPIDTPCTOF",
 140   "kINVMASS",
 141   "Selected e"
 142 };
 143
 144 int AliDxHFEParticleSelectionEl::Init()
 145 {
 146   //
 147   // init function
 148   //
 149   int iResult=0;
 150   // TODO: think about initializing fSelNHFE and fTrackCuts from the macro
 151   //Initialization of invariant mass cut function (AliSelectNonHFE)
 152   fSelNHFE= new AliSelectNonHFE("IM","IM");
 153   // Cuts for associated track in Inv Mass cut
 154   fTrackCuts = new AliESDtrackCuts();
 155   fTrackCuts->SetAcceptKinkDaughters(kFALSE);
 156   fTrackCuts->SetRequireTPCRefit(kTRUE);
 157   fTrackCuts->SetEtaRange(-0.9,0.9);
 158   fTrackCuts->SetRequireSigmaToVertex(kTRUE);
 159   fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
 160   fTrackCuts->SetMinNClustersTPC(80);
 161   fTrackCuts->SetPtRange(0.3,1e10);
 162
 163   fSelNHFE->SetTrackCuts(-3, 3, fTrackCuts);
 164   fSelNHFE->SetInvariantMassCut(fInvMassLow);
 165   fSelNHFE->SetAlgorithm("KF");
 166   fSelNHFE->SetAODanalysis(kTRUE);
 167
 168
 169   // Implicit call to InitControlObjects() before setting up CFM and fCuts
 170   // (if not there)
 171   iResult=AliDxHFEParticleSelection::Init();
 172   if (iResult<0) return iResult;
 173
 174   //--------Initialize correction Framework and Cuts-------------------------
 175   // Consider moving this, either to separate function or
 176   // add a set function for AliCFManager
 177   // Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
 178   AliInfo("Setting up CFM");
 179   fCFM = new AliCFManager;
 180   // the setup of cut objects is done in AliHFEcuts::Initialize
 181   // the ids used by this class must be the same, the code might be broken if
 182   // the sequence in AliHFEcuts::Initialize is changed
 183   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
 184   // reset pointers in the CF manager
 185   fCFM->SetNStepParticle(kNcutSteps);
 186   for(Int_t istep = 0; istep < kNcutSteps; istep++) {
 187     fCFM->SetParticleCutsList(istep, NULL);
 188   }
 189   if(!fCuts) {
 190     AliWarning("Cuts not available. Default cuts will be used");
 191     fCuts = new AliHFEcuts;
 192     fCuts->CreateStandardCuts();
 193   }
 194   // TODO: error handling?
 195   fCuts->Initialize(fCFM);
 196
 197   //Setting up TPC PID
 198   //Add settings for asymmetric cut on nSigma TPC
 199   //TODO: have this completely set up from addtask
 200   const int paramSize=4;
 201   Double_t params[paramSize];
 202   memset(params, 0, sizeof(Double_t)*paramSize);
 203   params[0]=-1.;
 204   fPIDTPC = new AliHFEpid("hfePidTPC");
 205   if(!fPIDTPC->GetNumberOfPIDdetectors()) {
 206     fPIDTPC->AddDetector("TPC",1);
 207   }
 208   fPIDTPC->ConfigureTPCdefaultCut(NULL, params, 3.);
 209   fPIDTPC->InitializePID();
 210
 211   if(fUseInvMassCut>kNoInvMass)  AliDebug(2,Form("Setting up with invariant mass cut of %f",fInvMassLow));
 212   if(fStoreCutStepInfo) AliDebug(2,"will store cut step info");
 213   AliDebug(2,Form("Stopping selection after step: %d", fFinalCutStep));
 214
 215   return 0;
 216
 217 }
 218
 219 THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
 220 {
 221   //
 222   // Defines the THnSparse. For now, only calls CreatControlTHnSparse
 223   // TODO: Add also invariant mass? here and in correlation, to do cut afterwards..
 224
 225   // here is the only place to change the dimension
 226   const int thnSize = 3;
 227   InitTHnSparseArray(thnSize);
 228   const double Pi=TMath::Pi();
 229   TString name;
 230   name.Form("%s info", GetName());
 231
 232   //                                 0    1       2
 233   //                                 Pt   Phi    Eta
 234   int         thnBins [thnSize] = { 1000,  200, 500};
 235   double      thnMin  [thnSize] = {    0,    0, -1.};
 236   double      thnMax  [thnSize] = {  100, 2*Pi,  1.};
 237   const char* thnNames[thnSize] = { "Pt","Phi","Eta"};
 238
 239   return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
 240 }
 241
 242 int AliDxHFEParticleSelectionEl::InitControlObjects()
 243 {
 244   /// init control and monitoring objects
 245   AliInfo("Electron THnSparse");
 246
 247   fElectronProperties=DefineTHnSparse();
 248   AddControlObject(fElectronProperties);
 249
 250   fHistoList=new TList;
 251   fHistoList->SetName("HFE Histograms");
 252   fHistoList->SetOwner();
 253
 254   // Histogram storing which cuts have been applied to the tracks
 255   // TODO: revise the names of the histograms, starting with 'f'
 256   // confuses the reader to think of class members
 257   fHistoList->Add(CreateControlHistogram("fWhichCut","effective cut for a rejected particle", kNCutLabels, fgkCutBinNames));
 258   fHistoList->Add(CreateControlHistogram("fTPCnClAOD","Nr TPC clusters/track for all tracks", 160, 0., 159.));
 259   fHistoList->Add(CreateControlHistogram("fTPCnClSingleTrackCuts","Nr TPC clusters/track After SingleTrackCuts", 160, 0., 159.));
 260   fHistoList->Add(CreateControlHistogram("fTPCnClTPCTOFPID","Nr TPC clusters/track After TPCTOF PID", 160, 0., 159.));
 261
 262   double dEdxBins[6]={1000,0.,10.,200,0.,200.};
 263   double nSigBins[6]={1000,0.,10.,200,-10.,10.};
 264
 265   // dEdx plots, TPC signal vs momentum
 266   fHistoList->Add(CreateControl2DHistogram("fdEdx", "dEdx before cuts", dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 267   fHistoList->Add(CreateControl2DHistogram("fdEdxCut", "dEdx after cuts",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 268   fHistoList->Add(CreateControl2DHistogram("fdEdxPidTOF", "dEdx after TOF pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 269   fHistoList->Add(CreateControl2DHistogram("fdEdxPidTPC", "dEdx after TPC pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 270   fHistoList->Add(CreateControl2DHistogram("fdEdxPid", "dEdx after pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 271   fHistoList->Add(CreateControl2DHistogram("fdEdxIM", "dEdx after Inv Mass",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 272
 273   // nSigmaTPC vs momentum
 274   fHistoList->Add(CreateControl2DHistogram("fnSigTPC", "nSigTPC before cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 275   fHistoList->Add(CreateControl2DHistogram("fnSigTPCCut", "nSigmaTPC after cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 276   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPidTPC", "nSigmaTPC after TPC PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 277   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPidTOF", "nSigmaTPC after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 278   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPid", "nSigmaTPC after PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 279
 280   // nSigmaTOF vs momentum
 281   fHistoList->Add(CreateControl2DHistogram("fnSigTOF", "nSigmaTOF before cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 282   fHistoList->Add(CreateControl2DHistogram("fnSigTOFCut", "nSigmaTOF after cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 283   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPidTOF", "nSigmaTOF after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 284   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPidTPC", "nSigmaTOF after TPC PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 285   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPid", "nSigmaTOF after PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 286
 287   // Invariant mass LS and ULS without cut
 288   fHistoList->Add(CreateControlHistogram("fInvMassLS", "Invariant mass LS", 1000, 0., 0.5));
 289   fHistoList->Add(CreateControlHistogram("fInvMassULS", "Invariant mass ULS", 1000, 0., 0.5));
 290
 291   // Invariant mass LS and ULS without cut, extended x-axis
 292   fHistoList->Add(CreateControlHistogram("fInvMass2SelLS", "Invariant mass two selected particles LS", 1000, 0., 5.));
 293   fHistoList->Add(CreateControlHistogram("fInvMass2SelULS", "Invariant mass two selected particles ULS", 1000, 0., 5.));
 294
 295   // Invariant mass LS and ULS with cut
 296   // TODO: Remove if remove first try at invariant mass
 297   fHistoList->Add(CreateControlHistogram("fInvMass2SelLScut", "Invariant mass two selected particles LS (cut)", 1000, 0., 0.5));
 298   fHistoList->Add(CreateControlHistogram("fInvMass2SelULScut", "Invariant mass two selected particles ULS (cut)", 1000, 0., 0.5));
 299
 300   fSelNHFE->SetHistMass((TH1F*)fHistoList->FindObject("fInvMassULS"));
 301   fSelNHFE->SetHistMassBack((TH1F*)fHistoList->FindObject("fInvMassLS"));
 302
 303   AddControlObject(fHistoList);
 304
 305   return AliDxHFEParticleSelection::InitControlObjects();
 306 }
 307
 308 int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
 309 {
 310   /// histogram particle properties
 311   if (!p) return -EINVAL;
 312   //if (!fControlObjects) return 0;
 313   if(selectionCode==0) return  0;
 314
 315   if(fElectronProperties && ParticleProperties()) {
 316     FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
 317     fElectronProperties->Fill(ParticleProperties());
 318   }
 319
 320   return 0;
 321 }
 322
 323 int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 324 {
 325   // fill the data array from the particle data
 326   if (!data) return -EINVAL;
 327   AliAODTrack *track=(AliAODTrack*)p;
 328   if (!track) return -ENODATA;
 329   int i=0;
 330   if (dimension!=GetDimTHnSparse()) {
 331     // TODO: think about filling only the available data and throwing a warning
 332     return -ENOSPC;
 333   }
 334   data[i++]=track->Pt();
 335   data[i++]=track->Phi();
 336   data[i++]=track->Eta();
 337   return i;
 338 }
 339
 340
 341 TObjArray* AliDxHFEParticleSelectionEl::Select(const AliVEvent* pEvent)
 342 {
 343   /// create selection from 'Tracks' member of the event,
 344   /// array contains only pointers but does not own the objects
 345   /// object array needs to be deleted by caller
 346   if (!pEvent) return NULL;
 347
 348   TList* selectedTracks=new TList;
 349   fSelNHFE->SetPIDresponse(fPIDResponse);
 350   fPIDTPC->SetPIDResponse(fPIDResponse);
 351   TObjArray* finalTracks=new TObjArray;
 352   if (!finalTracks) return NULL;
 353   finalTracks->SetOwner(kFALSE); // creating new track objects below
 354   TObjArray* finalTracksIM=new TObjArray;
 355   if (!finalTracksIM) return NULL;
 356   finalTracksIM->SetOwner(kFALSE); // creating new track objects below
 357
 358   int nofTracks=pEvent->GetNumberOfTracks();
 359   fTrackNum=0;
 360   int selectionCode=0;
 361   //  int origin
 362   for (int itrack=0; itrack<nofTracks; itrack++) {
 363     selectionCode=0;
 364     AliVParticle* track=pEvent->GetTrack(itrack);
 365     selectionCode=IsSelected(track,pEvent);
 366     fTrackNum++;
 367     // If fUseInvMassCut is set to true, only call HistogramParticleProperties
 368     // (here) for those who didn't pass single track cuts and pid.
 369     // if fUseInvMassCut is not set, then call the function for all tracks
 370     // TODO: Change this when have added invariant mass i thnsparse
 371     //CHANGED FOR OLD INV MASS
 372     if (selectionCode==0) continue;
 373     AliReducedParticle *trackReduced2 =(AliReducedParticle*)CreateParticle(track);
 374     finalTracks->Add(trackReduced2);
 375     if(fUseInvMassCut!=kInvMassTwoSelected || selectionCode==0)
 376       HistogramParticleProperties(trackReduced2, selectionCode);
 377     selectedTracks->Add(track);
 378   }
 379
 380   // Calculating invariant mass for electron pairs with same selection criteria
 381   // TODO: Could be removed if it is verified that looser cuts is better, but keep
 382   // for now
 383   if(selectedTracks->GetSize()>0)
 384     {
 385       Bool_t* selTrackIndex=new Bool_t[selectedTracks->GetSize()];
 386       if (selTrackIndex) {
 387       InvMassFilter(selectedTracks, selTrackIndex);
 388
 389       //Only remove electrons if fUseInvMassCut is set
 390       if(fUseInvMassCut==kInvMassTwoSelected){
 391         for(Int_t k=0; k<selectedTracks->GetSize(); k++)
 392           {
 393             selectionCode=0; //Reset selectionCode here to be based on invariant mass cut
 394             //On the basis that selectedTracks and finalTracks contain same particles
 395             AliAODTrack *trackIM=(AliAODTrack*)selectedTracks->At(k);
 396             AliReducedParticle *trackReduced=(AliReducedParticle*)finalTracks->At(k);
 397             if(selTrackIndex[k]==kTRUE)
 398               {
 399                 ((TH2D*)fHistoList->FindObject("fdEdxIM"))->Fill(trackIM->GetTPCmomentum(), trackIM->GetTPCsignal());
 400                 finalTracksIM->Add(trackReduced);
 401                 selectionCode=1;
 402               }
 403             HistogramParticleProperties(trackReduced, selectionCode);
 404           }
 405       }
 406       delete [] selTrackIndex;
 407       }
 408
 409     }
 410
 411   HistogramEventProperties(AliDxHFEParticleSelection::kHistoNrTracksPrEvent,finalTracks->GetEntries());
 412   if(fUseInvMassCut!=kInvMassTwoSelected)
 413     return finalTracks;
 414   else
 415     return finalTracksIM;
 416
 417 }
 418
 419 int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
 420 {
 421   /// select El candidates
 422   if(!pEvent){
 423     AliError("No event information");
 424     return 0;
 425   }
 426   fSurvivedCutStep=kNotSelected;
 427
 428   AliAODTrack *track=(AliAODTrack*)pEl;
 429   fCFM->SetRecEventInfo(pEvent);
 430
 431   ((TH2D*)fHistoList->FindObject("fdEdx"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 432   ((TH2D*)fHistoList->FindObject("fnSigTPC"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 433   ((TH2D*)fHistoList->FindObject("fnSigTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 434   ((TH1D*)fHistoList->FindObject("fTPCnClAOD"))->Fill(track->GetTPCNcls());
 435
 436
 437   if(fFinalCutStep==kNoCuts){
 438     Float_t radial=999;
 439     //Float_t z=999;
 440     Bool_t acceptTrack=kFALSE;
 441
 442     //Cut on Eta
 443     if(TMath::Abs(track->Eta()) < fEtaCut) {
 444       acceptTrack=kTRUE;
 445     }
 446
 447     const Double_t kBeampiperadius=fImpactParamCutRadial;
 448     Double_t dcaD[2]={-999.,-999.}, covD[3]={-999.,-999.,-999.};
 449
 450     AliAODEvent *aodevent = (AliAODEvent*)(pEvent);
 451     if(!aodevent) {
 452       AliDebug(1, "No aod event available\n");
 453       return 0;
 454     }
 455
 456     AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
 457     if(!vtxAODSkip) return 0;
 458     AliExternalTrackParam etp; etp.CopyFromVTrack(track);
 459     if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
 460       radial = dcaD[0];
 461       //z = dcaD[1];
 462     }
 463
 464     // Cut on z
 465     if(acceptTrack && TMath::Abs(radial) < fImpactParamCutRadial){
 466       acceptTrack=kTRUE;
 467     }
 468
 469     if(acceptTrack)
 470       fSurvivedCutStep=kNoCuts;
 471
 472     if(!acceptTrack) return 0;
 473
 474     // Should only return at this point if you only want to store all tracks
 475     AliDebug(2,"Returns after kNoCuts");
 476     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
 477     return 1;
 478   }
 479
 480   //--------track cut selection-----------------------
 481   //Using AliHFECuts:
 482   // RecKine: ITSTPC cuts
 483   if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
 484     AliDebug(4,"Cut: kStepRecKineITSTPC");
 485     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecKineITSTPC);
 486     return 0;
 487   }
 488   if(fStoreCutStepInfo) fSurvivedCutStep=kRecKineITSTPC;
 489   if(fFinalCutStep==kRecKineITSTPC) {AliDebug(2,"Returns after kStepRecKineITSTPC"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 490
 491   // RecPrim
 492   if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
 493     AliDebug(4,"Cut: kStepRecPrim");
 494     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecPrim);
 495     if(!fStoreCutStepInfo) return 0;
 496     else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
 497   }
 498   if(fStoreCutStepInfo) fSurvivedCutStep=kRecPrim;
 499   if(fFinalCutStep==kRecPrim) {AliDebug(2,"Returns after kRecPrim"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 500
 501   // HFEcuts: ITS layers cuts
 502   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
 503     AliDebug(4,"Cut: kStepHFEcutsITS");
 504     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsITS);
 505     if(!fStoreCutStepInfo) return 0;
 506     else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
 507
 508   }
 509   if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsITS;
 510   if(fFinalCutStep==kHFEcutsITS) {AliDebug(2,"Returns after kHFEcutsITS "); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 511
 512   // HFE cuts: TOF PID and mismatch flag
 513   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
 514     AliDebug(4,"Cut: kStepHFEcutsTOF");
 515     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTOF);
 516     if(!fStoreCutStepInfo) return 0;
 517     else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
 518   }
 519   if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTOF;
 520   if(fFinalCutStep==kHFEcutsTOF) {AliDebug(2,"Returns after kHFEcutsTOF"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 521
 522   // HFE cuts: TPC PID cleanup
 523   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
 524     AliDebug(4,"Cut: kStepHFEcutsTPC");
 525     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTPC);
 526     if(!fStoreCutStepInfo) return 0;
 527     else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
 528   }
 529   if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTPC;
 530
 531   ((TH2D*)fHistoList->FindObject("fdEdxCut"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 532   ((TH2D*)fHistoList->FindObject("fnSigTPCCut"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 533   ((TH2D*)fHistoList->FindObject("fnSigTOFCut"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 534   ((TH1D*)fHistoList->FindObject("fTPCnClSingleTrackCuts"))->Fill(track->GetTPCNcls());
 535
 536   if(fFinalCutStep==kHFEcutsTPC) {AliDebug(2,"Returns after track cuts"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 537
 538
 539   //--------PID selection-----------------------
 540   AliHFEpidObject hfetrack;
 541   hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
 542   hfetrack.SetRecTrack(track);
 543
 544   // TODO: configurable colliding system
 545   // TODO: Check problem with PbPb, for now workaround with setting multiplicity to -1
 546   if(GetSystem()==1) {
 547     hfetrack.SetPbPb();
 548     hfetrack.SetCentrality(-1);
 549   }
 550   else  hfetrack.SetPP();
 551   //  hfetrack.SetMulitplicity(ncontribVtx);
 552
 553   // TODO: Put this into a while-loop instead, looping over the number of pid objects in the cut-list?
 554   // This needs a bit of thinking and finetuning (wrt histogramming)
 555   if(fPIDTOF && fPIDTOF->IsSelected(&hfetrack)) {
 556     ((TH2D*)fHistoList->FindObject("fdEdxPidTOF"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 557     ((TH2D*)fHistoList->FindObject("fnSigTPCPidTOF"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 558     ((TH2D*)fHistoList->FindObject("fnSigTOFPidTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 559     if(fStoreCutStepInfo){fSurvivedCutStep=kPIDTOF; }
 560     if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF");((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 561   }
 562   else{
 563     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOF);
 564   }
 565
 566   if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF"); return 0;}
 567
 568
 569   if(fPIDTPC && fPIDTPC->IsSelected(&hfetrack)) {
 570     ((TH2D*)fHistoList->FindObject("fdEdxPidTPC"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 571     ((TH2D*)fHistoList->FindObject("fnSigTPCPidTPC"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 572     ((TH2D*)fHistoList->FindObject("fnSigTOFPidTPC"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 573     if(fStoreCutStepInfo){fSurvivedCutStep=kPIDTPC; }
 574     if(fFinalCutStep==kPIDTPC) {AliDebug(2,"Returns at PIDTPC"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 575   }
 576   else{
 577     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTPC);
 578     //return 0; //if rejected by TOF, will also be rejected by TPCTOF
 579   }
 580   if(fFinalCutStep==kPIDTPC) {AliDebug(2,"Returns at PIDTTPC"); return 0;}
 581
 582   //Combined tof & tpc pid
 583   if(fPIDTOFTPC && fPIDTOFTPC->IsSelected(&hfetrack)) {
 584     AliDebug(3,"Inside FilldPhi, electron is selected");
 585     ((TH2D*)fHistoList->FindObject("fdEdxPid"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
 586     ((TH2D*)fHistoList->FindObject("fnSigTPCPid"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
 587     ((TH2D*)fHistoList->FindObject("fnSigTOFPid"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
 588     ((TH1D*)fHistoList->FindObject("fTPCnClTPCTOFPID"))->Fill(track->GetTPCNcls());
 589     // Only do this if
 590     if(fStoreCutStepInfo){
 591       fSurvivedCutStep=kPIDTOFTPC;
 592
 593     }
 594   }
 595   else{
 596     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOFTPC);
 597     AliDebug(4,"Cut: kTPCTOFPID");
 598     if(!fStoreCutStepInfo) return 0;
 599     else { return 1; }//return 1 because it passed cuts above, but not this (no need to go further)
 600   }
 601
 602   if(fUseInvMassCut==kInvMassSingleSelected)
 603     {
 604       AliDebug(4,"invmass check");
 605       fSelNHFE->FindNonHFE(fTrackNum, track, const_cast<AliVEvent*>(pEvent));
 606       if(fSelNHFE->IsLS() || fSelNHFE->IsULS())
 607         {
 608           //Not selected
 609           ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kINVMASS);
 610           AliDebug(4,"Cut: Invmass");
 611           if(!fStoreCutStepInfo) return 0;
 612         }
 613       else
 614         if(fStoreCutStepInfo) fSurvivedCutStep=kINVMASS;
 615
 616     }
 617   ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
 618   return 1;
 619
 620 }
 621
 622 void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
 623 {
 624   /// set cut objects
 625   if (level==kCutHFE) {
 626     fCuts=dynamic_cast<AliHFEcuts*>(cuts);
 627     if (!fCuts && cuts) {
 628       AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
 629     }
 630     return;
 631   }
 632
 633   if (level==kCutPIDTOFTPC) {
 634     fPIDTOFTPC=dynamic_cast<AliHFEpid*>(cuts);
 635     if (!fPIDTOFTPC && cuts) {
 636       AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 637     }
 638     return;
 639   }
 640
 641   if (level==kCutPIDTOF) {
 642     fPIDTOF=dynamic_cast<AliHFEpid*>(cuts);
 643     if (!fPIDTOF && cuts) {
 644       AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 645     }
 646     return;
 647   }
 648   if(level==kCutList){
 649     fCutPidList=dynamic_cast<TList*>(cuts);
 650     if (!fCutPidList && cuts) {
 651       AliError(Form("cuts object is not of required type TList but %s", cuts->ClassName()));
 652     }
 653     else{
 654       // TODO: Could be done more elegantly, at the moment requires that the cut and pid objects are in
 655       // a specific order..
 656       TObject *obj=NULL;
 657       int iii=0;
 658       TIter next(fCutPidList);
 659       while((obj = next())){
 660         iii++;
 661         if(iii==1) {
 662           fCuts=dynamic_cast<AliHFEcuts*>(obj);
 663           if (!fCuts)
 664             AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
 665         }
 666         if(iii==2){
 667           fPIDTOFTPC=dynamic_cast<AliHFEpid*>(obj);
 668           if (!fPIDTOFTPC)
 669             AliError(Form("(TOFTPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 670         }
 671         if(iii==3){
 672           fPIDTOF=dynamic_cast<AliHFEpid*>(obj);
 673           if (!fPIDTOF)
 674             AliError(Form("(TOF) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 675         }
 676         /*if(iii=4){
 677           fPIDTPC=dynamic_cast<AliHFEpid*>(obj);
 678           if (!fPIDTPC)
 679             AliError(Form("(TPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
 680             }*/
 681       }
 682
 683     }
 684     return;
 685   }
 686 }
 687
 688 //________________________________________________________________________
 689 Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 690 {
 691   // Check single track cuts for a given cut step
 692   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
 693   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 694   //if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
 695   return kTRUE;
 696 }
 697
 698 int AliDxHFEParticleSelectionEl::ParseArguments(const char* arguments)
 699 {
 700   // parse arguments and set internal flags
 701   TString strArguments(arguments);
 702   auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
 703   if (!tokens.get()) return 0;
 704
 705   AliInfo(strArguments);
 706   TIter next(tokens.get());
 707   TObject* token;
 708   while ((token=next())) {
 709     TString argument=token->GetName();
 710     if(argument.BeginsWith("elreco=")){
 711       argument.ReplaceAll("elreco=", "");
 712       int selectionStep=kPIDTOFTPC; //Default
 713       if(argument.CompareTo("alltracks")==0) selectionStep=kNoCuts;
 714       if(argument.CompareTo("afterreckineitstpc")==0) selectionStep=kRecKineITSTPC;
 715       if(argument.CompareTo("afterrecprim")==0) selectionStep=kRecPrim;
 716       if(argument.CompareTo("afterhfeits")==0) selectionStep=kHFEcutsITS;
 717       if(argument.CompareTo("afterhfetof")==0) selectionStep=kHFEcutsTOF;
 718       if(argument.CompareTo("afterhfetpc")==0) selectionStep=kHFEcutsTPC;
 719       if(argument.CompareTo("aftertrackcuts")==0) selectionStep=kHFEcutsTPC;
 720       if(argument.CompareTo("aftertofpid")==0) selectionStep=kPIDTOF;
 721       if(argument.CompareTo("aftertpcpid")==0) selectionStep=kPIDTPC;
 722       if(argument.CompareTo("afterfullpid")==0) selectionStep=kPIDTOFTPC;
 723       SetFinalCutStep(selectionStep);
 724       continue;
 725     }
 726     if(argument.BeginsWith("twoselectedinvmasscut")){
 727       fUseInvMassCut=kInvMassTwoSelected;
 728       AliInfo("Using Invariant mass cut for two selected particles");
 729       continue;
 730     }
 731     if(argument.BeginsWith("useinvmasscut")){
 732       fUseInvMassCut=kInvMassSingleSelected;
 733       AliInfo("Using Invariant mass cut for single selected particle and looser cuts on partner");
 734       continue;
 735     }
 736     if(argument.BeginsWith("invmasscut=")){
 737       argument.ReplaceAll("invmasscut=", "");
 738       fInvMassLow=argument.Atof();
 739       AliInfo(Form("Using invariant mass cut: %f",fInvMassLow));
 740       //      fUseInvMassCut=kInvMassSingleSelected;
 741       continue;
 742     }
 743     if(argument.BeginsWith("impactparamcut=")){
 744       argument.ReplaceAll("impactparamcut=", "");
 745       fImpactParamCutRadial=argument.Atof();
 746       AliInfo(Form("Using impact parameter cut: %f",fImpactParamCutRadial));
 747       continue;
 748     }
 749     if(argument.BeginsWith("etacut=")){
 750       argument.ReplaceAll("etacut=", "");
 751       fEtaCut=argument.Atof();
 752       AliInfo(Form("Using Eta cut: %f",fEtaCut));
 753       continue;
 754     }
 755     if(argument.BeginsWith("storelastcutstep")){
 756       AliInfo("Stores the last cut step");
 757       SetStoreLastCutStep(kTRUE);
 758       continue;
 759     }
 760     // forwarding of single argument works, unless key-option pairs separated
 761     // by blanks are introduced
 762     AliDxHFEParticleSelection::ParseArguments(argument);
 763   }
 764
 765   return 0;
 766 }
 767
 768
 769 void AliDxHFEParticleSelectionEl::InvMassFilter(TList *elList, Bool_t *selIndx)
 770 {
 771
 772   //Function for getting invariant mass of electron pairs
 773
 774   for(int i=0; i<((elList->GetSize())-1); i++)
 775     {
 776       AliAODTrack *trackAsso=(AliAODTrack*)elList->At(i);
 777       for(int j=i+1; j<elList->GetSize(); j++)
 778         {
 779
 780           AliAODTrack *track=(AliAODTrack*)elList->At(j);
 781           if(trackAsso && track)
 782             {
 783
 784               Double_t mass=-999., width = -999;
 785               Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
 786
 787               Int_t chargeAsso = trackAsso->Charge();
 788               Int_t charge = track->Charge();
 789
 790               Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
 791               if(charge>0) fPDGe1 = -11;
 792               if(chargeAsso>0) fPDGe2 = -11;
 793
 794               if(charge == chargeAsso) fFlagLS = kTRUE;
 795               if(charge != chargeAsso) fFlagULS = kTRUE;
 796
 797               AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
 798               AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
 799               AliKFParticle recg(ge1, ge2);
 800
 801               if(recg.GetNDF()<1) continue;
 802               Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
 803               if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
 804
 805               recg.GetMass(mass,width);
 806
 807               if(fFlagLS) {
 808                 //              if(mass<0.5)
 809                 // ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
 810                 ((TH1D*)fHistoList->FindObject("fInvMass2SelLS"))->Fill(mass);
 811                 if(mass>fInvMassLow)
 812                   {
 813                     ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
 814                     selIndx[i]=kTRUE;
 815                     selIndx[j]=kTRUE;
 816                   }
 817               }
 818               if(fFlagULS) {
 819                 //if(mass<0.5)
 820                 //  ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);
 821                 ((TH1D*)fHistoList->FindObject("fInvMass2SelULS"))->Fill(mass);
 822                 if(mass>fInvMassLow)
 823                   {
 824                     ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);
 825                     selIndx[i]=kTRUE;
 826                     selIndx[j]=kTRUE;
 827                   }
 828
 829               }
 830
 831             }
 832         }
 833     }
 834 }