3 //**************************************************************************
4 //* This file is property of and copyright by the ALICE Project *
5 //* ALICE Experiment at CERN, All rights reserved. *
7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no> *
8 //* Sedat Altinpinar <Sedat.Altinpinar@cern.ch> *
9 //* Hege Erdal <hege.erdal@gmail.com> *
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 //**************************************************************************
20 /// @file AliDxHFEParticleSelectionEl.cxx
21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
23 /// @brief D0 selection for D0-HFE correlation
25 #include "AliDxHFEParticleSelectionEl.h"
26 #include "AliSelectNonHFE.h"
27 #include "AliReducedParticle.h"
28 #include "AliESDtrackCuts.h"
29 #include "AliVParticle.h"
30 #include "AliVEvent.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"
53 #include "TObjArray.h"
54 #include "AliKFParticle.h"
55 #include "AliKFVertex.h"
62 /// ROOT macro for the implementation of ROOT specific class methods
63 ClassImp(AliDxHFEParticleSelectionEl)
65 AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
66 : AliDxHFEParticleSelection("electron", opt)
70 , fElectronProperties(NULL)
78 , fFinalCutStep(kPIDTOFTPC)
80 , fUseInvMassCut(kNoInvMass)
83 , fImpactParamCutRadial(3)
85 , fSurvivedCutStep(kNotSelected)
86 , fStoreCutStepInfo(kFALSE)
97 AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
100 if (fElectronProperties) {
101 delete fElectronProperties;
102 fElectronProperties=NULL;
116 // NOTE: external objects fPID and fCuts are not deleted here
131 const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
144 int AliDxHFEParticleSelectionEl::Init()
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);
163 fSelNHFE->SetTrackCuts(-3, 3, fTrackCuts);
164 fSelNHFE->SetInvariantMassCut(fInvMassLow);
165 fSelNHFE->SetAlgorithm("KF");
166 fSelNHFE->SetAODanalysis(kTRUE);
169 // Implicit call to InitControlObjects() before setting up CFM and fCuts
171 iResult=AliDxHFEParticleSelection::Init();
172 if (iResult<0) return iResult;
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);
190 AliWarning("Cuts not available. Default cuts will be used");
191 fCuts = new AliHFEcuts;
192 fCuts->CreateStandardCuts();
194 // TODO: error handling?
195 fCuts->Initialize(fCFM);
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);
204 fPIDTPC = new AliHFEpid("hfePidTPC");
205 if(!fPIDTPC->GetNumberOfPIDdetectors()) {
206 fPIDTPC->AddDetector("TPC",1);
208 fPIDTPC->ConfigureTPCdefaultCut(NULL, params, 3.);
209 fPIDTPC->InitializePID();
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));
219 THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse()
222 // Defines the THnSparse. For now, only calls CreatControlTHnSparse
223 // TODO: Add also invariant mass? here and in correlation, to do cut afterwards..
225 // here is the only place to change the dimension
226 const int thnSize = 3;
227 InitTHnSparseArray(thnSize);
228 const double Pi=TMath::Pi();
230 name.Form("%s info", GetName());
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"};
239 return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
242 int AliDxHFEParticleSelectionEl::InitControlObjects()
244 /// init control and monitoring objects
245 AliInfo("Electron THnSparse");
247 fElectronProperties=DefineTHnSparse();
248 AddControlObject(fElectronProperties);
250 fHistoList=new TList;
251 fHistoList->SetName("HFE Histograms");
252 fHistoList->SetOwner();
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.));
262 double dEdxBins[6]={1000,0.,10.,200,0.,200.};
263 double nSigBins[6]={1000,0.,10.,200,-10.,10.};
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.)"));
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.)"));
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.)"));
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));
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.));
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));
300 fSelNHFE->SetHistMass((TH1F*)fHistoList->FindObject("fInvMassULS"));
301 fSelNHFE->SetHistMassBack((TH1F*)fHistoList->FindObject("fInvMassLS"));
303 AddControlObject(fHistoList);
305 return AliDxHFEParticleSelection::InitControlObjects();
308 int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
310 /// histogram particle properties
311 if (!p) return -EINVAL;
312 //if (!fControlObjects) return 0;
313 if(selectionCode==0) return 0;
315 if(fElectronProperties && ParticleProperties()) {
316 FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
317 fElectronProperties->Fill(ParticleProperties());
323 int AliDxHFEParticleSelectionEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
325 // fill the data array from the particle data
326 if (!data) return -EINVAL;
327 AliAODTrack *track=(AliAODTrack*)p;
328 if (!track) return -ENODATA;
330 if (dimension!=GetDimTHnSparse()) {
331 // TODO: think about filling only the available data and throwing a warning
334 data[i++]=track->Pt();
335 data[i++]=track->Phi();
336 data[i++]=track->Eta();
341 TObjArray* AliDxHFEParticleSelectionEl::Select(const AliVEvent* pEvent)
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;
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
358 int nofTracks=pEvent->GetNumberOfTracks();
362 for (int itrack=0; itrack<nofTracks; itrack++) {
364 AliVParticle* track=pEvent->GetTrack(itrack);
365 selectionCode=IsSelected(track,pEvent);
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);
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
383 if(selectedTracks->GetSize()>0)
385 Bool_t* selTrackIndex=new Bool_t[selectedTracks->GetSize()];
387 InvMassFilter(selectedTracks, selTrackIndex);
389 //Only remove electrons if fUseInvMassCut is set
390 if(fUseInvMassCut==kInvMassTwoSelected){
391 for(Int_t k=0; k<selectedTracks->GetSize(); k++)
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)
399 ((TH2D*)fHistoList->FindObject("fdEdxIM"))->Fill(trackIM->GetTPCmomentum(), trackIM->GetTPCsignal());
400 finalTracksIM->Add(trackReduced);
403 HistogramParticleProperties(trackReduced, selectionCode);
406 delete [] selTrackIndex;
411 HistogramEventProperties(AliDxHFEParticleSelection::kHistoNrTracksPrEvent,finalTracks->GetEntries());
412 if(fUseInvMassCut!=kInvMassTwoSelected)
415 return finalTracksIM;
419 int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
421 /// select El candidates
423 AliError("No event information");
426 fSurvivedCutStep=kNotSelected;
428 AliAODTrack *track=(AliAODTrack*)pEl;
429 fCFM->SetRecEventInfo(pEvent);
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());
437 if(fFinalCutStep==kNoCuts){
440 Bool_t acceptTrack=kFALSE;
443 if(TMath::Abs(track->Eta()) < fEtaCut) {
447 const Double_t kBeampiperadius=fImpactParamCutRadial;
448 Double_t dcaD[2]={-999.,-999.}, covD[3]={-999.,-999.,-999.};
450 AliAODEvent *aodevent = (AliAODEvent*)(pEvent);
452 AliDebug(1, "No aod event available\n");
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)) {
465 if(acceptTrack && TMath::Abs(radial) < fImpactParamCutRadial){
470 fSurvivedCutStep=kNoCuts;
472 if(!acceptTrack) return 0;
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);
480 //--------track cut selection-----------------------
482 // RecKine: ITSTPC cuts
483 if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
484 AliDebug(4,"Cut: kStepRecKineITSTPC");
485 ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecKineITSTPC);
488 if(fStoreCutStepInfo) fSurvivedCutStep=kRecKineITSTPC;
489 if(fFinalCutStep==kRecKineITSTPC) {AliDebug(2,"Returns after kStepRecKineITSTPC"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
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)
498 if(fStoreCutStepInfo) fSurvivedCutStep=kRecPrim;
499 if(fFinalCutStep==kRecPrim) {AliDebug(2,"Returns after kRecPrim"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
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)
509 if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsITS;
510 if(fFinalCutStep==kHFEcutsITS) {AliDebug(2,"Returns after kHFEcutsITS "); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
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)
519 if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTOF;
520 if(fFinalCutStep==kHFEcutsTOF) {AliDebug(2,"Returns after kHFEcutsTOF"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
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)
529 if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTPC;
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());
536 if(fFinalCutStep==kHFEcutsTPC) {AliDebug(2,"Returns after track cuts"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
539 //--------PID selection-----------------------
540 AliHFEpidObject hfetrack;
541 hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
542 hfetrack.SetRecTrack(track);
544 // TODO: configurable colliding system
545 // TODO: Check problem with PbPb, for now workaround with setting multiplicity to -1
548 hfetrack.SetCentrality(-1);
550 else hfetrack.SetPP();
551 // hfetrack.SetMulitplicity(ncontribVtx);
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;}
563 ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOF);
566 if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF"); return 0;}
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;}
577 ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTPC);
578 //return 0; //if rejected by TOF, will also be rejected by TPCTOF
580 if(fFinalCutStep==kPIDTPC) {AliDebug(2,"Returns at PIDTTPC"); return 0;}
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());
590 if(fStoreCutStepInfo){
591 fSurvivedCutStep=kPIDTOFTPC;
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)
602 if(fUseInvMassCut==kInvMassSingleSelected)
604 AliDebug(4,"invmass check");
605 fSelNHFE->FindNonHFE(fTrackNum, track, const_cast<AliVEvent*>(pEvent));
606 if(fSelNHFE->IsLS() || fSelNHFE->IsULS())
609 ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kINVMASS);
610 AliDebug(4,"Cut: Invmass");
611 if(!fStoreCutStepInfo) return 0;
614 if(fStoreCutStepInfo) fSurvivedCutStep=kINVMASS;
617 ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
622 void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
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()));
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()));
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()));
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()));
654 // TODO: Could be done more elegantly, at the moment requires that the cut and pid objects are in
655 // a specific order..
658 TIter next(fCutPidList);
659 while((obj = next())){
662 fCuts=dynamic_cast<AliHFEcuts*>(obj);
664 AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
667 fPIDTOFTPC=dynamic_cast<AliHFEpid*>(obj);
669 AliError(Form("(TOFTPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
672 fPIDTOF=dynamic_cast<AliHFEpid*>(obj);
674 AliError(Form("(TOF) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
677 fPIDTPC=dynamic_cast<AliHFEpid*>(obj);
679 AliError(Form("(TPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
688 //________________________________________________________________________
689 Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
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;
698 int AliDxHFEParticleSelectionEl::ParseArguments(const char* arguments)
700 // parse arguments and set internal flags
701 TString strArguments(arguments);
702 auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
703 if (!tokens.get()) return 0;
705 AliInfo(strArguments);
706 TIter next(tokens.get());
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);
726 if(argument.BeginsWith("twoselectedinvmasscut")){
727 fUseInvMassCut=kInvMassTwoSelected;
728 AliInfo("Using Invariant mass cut for two selected particles");
731 if(argument.BeginsWith("useinvmasscut")){
732 fUseInvMassCut=kInvMassSingleSelected;
733 AliInfo("Using Invariant mass cut for single selected particle and looser cuts on partner");
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;
743 if(argument.BeginsWith("impactparamcut=")){
744 argument.ReplaceAll("impactparamcut=", "");
745 fImpactParamCutRadial=argument.Atof();
746 AliInfo(Form("Using impact parameter cut: %f",fImpactParamCutRadial));
749 if(argument.BeginsWith("etacut=")){
750 argument.ReplaceAll("etacut=", "");
751 fEtaCut=argument.Atof();
752 AliInfo(Form("Using Eta cut: %f",fEtaCut));
755 if(argument.BeginsWith("storelastcutstep")){
756 AliInfo("Stores the last cut step");
757 SetStoreLastCutStep(kTRUE);
760 // forwarding of single argument works, unless key-option pairs separated
761 // by blanks are introduced
762 AliDxHFEParticleSelection::ParseArguments(argument);
769 void AliDxHFEParticleSelectionEl::InvMassFilter(TList *elList, Bool_t *selIndx)
772 //Function for getting invariant mass of electron pairs
774 for(int i=0; i<((elList->GetSize())-1); i++)
776 AliAODTrack *trackAsso=(AliAODTrack*)elList->At(i);
777 for(int j=i+1; j<elList->GetSize(); j++)
780 AliAODTrack *track=(AliAODTrack*)elList->At(j);
781 if(trackAsso && track)
784 Double_t mass=-999., width = -999;
785 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
787 Int_t chargeAsso = trackAsso->Charge();
788 Int_t charge = track->Charge();
790 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
791 if(charge>0) fPDGe1 = -11;
792 if(chargeAsso>0) fPDGe2 = -11;
794 if(charge == chargeAsso) fFlagLS = kTRUE;
795 if(charge != chargeAsso) fFlagULS = kTRUE;
797 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
798 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
799 AliKFParticle recg(ge1, ge2);
801 if(recg.GetNDF()<1) continue;
802 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
803 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
805 recg.GetMass(mass,width);
809 // ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
810 ((TH1D*)fHistoList->FindObject("fInvMass2SelLS"))->Fill(mass);
813 ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
820 // ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);
821 ((TH1D*)fHistoList->FindObject("fInvMass2SelULS"))->Fill(mass);
824 ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);