Incrementing class version.
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliDxHFEParticleSelectionMCEl.cxx
CommitLineData
d731501a 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//* Hege Erdal <hege.erdal@gmail.com> *
9//* *
10//* Permission to use, copy, modify and distribute this software and its *
11//* documentation strictly for non-commercial purposes is hereby granted *
12//* without fee, provided that the above copyright notice appears in all *
13//* copies and that both the copyright notice and this permission notice *
14//* appear in the supporting documentation. The authors make no claims *
15//* about the suitability of this software for any purpose. It is *
16//* provided "as is" without express or implied warranty. *
17//**************************************************************************
18
19/// @file AliDxHFEParticleSelectionMCEl.cxx
20/// @author Hege Erdal, Matthias Richter
21/// @date 2012-07-19
22/// @brief MC El selection for D0-HFE correlation
23///
24
25#include "AliDxHFEParticleSelectionMCEl.h"
26#include "AliVParticle.h"
27#include "AliLog.h"
28#include "THnSparse.h"
29#include "AliAODMCParticle.h"
30#include "TH1F.h"
31#include "TAxis.h"
32#include "AliAODTrack.h"
dfe96b90 33#include "AliReducedParticle.h"
d731501a 34#include <iostream>
35#include <cerrno>
36#include <memory>
37
38using namespace std;
39using std::vector;
40
41
42/// ROOT macro for the implementation of ROOT specific class methods
43ClassImp(AliDxHFEParticleSelectionMCEl)
44
45AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
46 : AliDxHFEParticleSelectionEl(opt)
47 , fMCTools()
dfe96b90 48 , fPDGnotMCElectron(NULL)
fad30eb8 49 , fPDGNotHFMother(NULL)
d731501a 50 , fOriginMother(0)
51 , fResultMC(0)
b4779749 52 , fUseKine(kFALSE)
53 , fMotherPDGs(0)
54 , fUseMCReco(kFALSE)
55 , fSelectionStep(AliDxHFEParticleSelectionEl::kNoCuts)
fad30eb8 56 , fStoreCutStepInfo(kFALSE)
d731501a 57{
58 // constructor
59 //
60 //
61 //
62 //
63
64 fMCTools.~AliDxHFEToolsMC();
65
b4779749 66 ParseArguments(opt);
67
68 // This is also checked in AddTasks, but just for safety!
69 if(fUseMCReco && fUseKine) AliFatal("CAN'T SET BOTH usekine AND elmcreco AT THE SAME TIME");
70
d731501a 71 // TODO: argument scan, build tool options accordingly
72 // e.g. set mc mode first/last, skip control histograms
73 TString toolopt("pdg=11 mc-last");
b4779749 74 if(fUseKine) toolopt+=" usekine";
d731501a 75 new (&fMCTools) AliDxHFEToolsMC(toolopt);
76}
77
78//at the moment not used, keep for now (to be used when plotting)
79const char* AliDxHFEParticleSelectionMCEl::fgkPDGMotherBinLabels[]={
80 "d",
81 "u",
82 "s",
83 "c",
84 "b",
85 "gluon",
86 "gamma",
87 "#pi^{0}",
88 "#eta",
89 "proton",
90 "others"
91};
92
dfe96b90 93
94const char* AliDxHFEParticleSelectionMCEl::fgkPDGBinLabels[]={
95 "positron",
96 "electron",
97 "#mu+",
98 "#mu-",
99 "#pi+",
100 "#pi-",
101 "K+",
102 "K-",
103 "proton",
104 "antiproton",
105 "others"
106};
107
d731501a 108AliDxHFEParticleSelectionMCEl::~AliDxHFEParticleSelectionMCEl()
109{
110 // destructor
dfe96b90 111
112 if(fPDGnotMCElectron){
113 delete fPDGnotMCElectron;
114 fPDGnotMCElectron=NULL;
115 }
116
117}
118
119int AliDxHFEParticleSelectionMCEl::Init()
120{
121 //
122 // init function
123 //
b4779749 124
125 // Particles considered HFE background. can be expanded
126 fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGpi0);
127 fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGeta);
128 fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGgamma);
129 fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGJpsi);
130
dfe96b90 131 int iResult=0;
132 iResult=AliDxHFEParticleSelectionEl::Init();
133 if (iResult<0) return iResult;
134
135 // Histo containing PDG of track which was not MC truth electron
b4779749 136 // TODO: Add to list in baseclass
137 fPDGnotMCElectron=CreateControlHistogram("fPDGnotMCElectron","PDG of track not MC truth electron",AliDxHFEToolsMC::kNofPDGLabels,fgkPDGBinLabels);
fad30eb8 138 fPDGNotHFMother=CreateControlHistogram("fPDGNotHFMother","PDG of mother not HF",5000);
dfe96b90 139 AddControlObject(fPDGnotMCElectron);
fad30eb8 140 AddControlObject(fPDGNotHFMother);
dfe96b90 141 return 0;
d731501a 142}
143
dcf83226 144THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
d731501a 145{
146 //
147 // Defines the THnSparse.
148
d731501a 149 const double Pi=TMath::Pi();
150 TString name;
fad30eb8 151 THnSparse* thn=NULL;
d731501a 152 name.Form("%s info", GetName());
153
fad30eb8 154 if(fStoreCutStepInfo){
155 const int thnSizeExt =5;
156
157 InitTHnSparseArray(thnSizeExt);
158
159 // TODO: Redo binning of distributions more?
160 // 0 1 2 3
161 // Pt Phi Eta mother
162 int thnBinsExt[thnSizeExt] = { 100, 100, 100, 15, kNCutLabels };
163 double thnMinExt [thnSizeExt] = { 0, 0, -1., -1.5, kRecKineITSTPC };
164 double thnMaxExt [thnSizeExt] = { 10, 2*Pi, 1., 13.5, kSelected };
165 const char* thnNamesExt[thnSizeExt]={
166 "Pt",
167 "Phi",
168 "Eta",
169 "Mother", //bin==-1: Not MC truth electron
170 "Last survived cut step"
171 };
172 thn=(THnSparse*)CreateControlTHnSparse(name,thnSizeExt,thnBinsExt,thnMinExt,thnMaxExt,thnNamesExt);
173 }
174 else{
175
176 const int thnSize =4;
177 InitTHnSparseArray(thnSize);
178
179 // TODO: Redo binning of distributions more?
180 // 0 1 2 3
181 // Pt Phi Eta mother
182 int thnBins[thnSize] = { 100, 100, 100, 15 };
183 double thnMin [thnSize] = { 0, 0, -1., -1.5 };
184 double thnMax [thnSize] = { 10, 2*Pi, 1., 13.5 };
185 const char* thnNames[thnSize]={
186 "Pt",
187 "Phi",
188 "Eta",
189 "Mother", //bin==-1: Not MC truth electron
190 };
191 thn=(THnSparse*)CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
192
193 }
194 return thn;
d731501a 195}
196
dcf83226 197int AliDxHFEParticleSelectionMCEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
d731501a 198{
199 // fill the data array from the particle data
200 if (!data) return -EINVAL;
fad30eb8 201 if (!p) return -ENODATA;
202 // handle different types of tracks, can be extended
203 AliReducedParticle *trRP=dynamic_cast<AliReducedParticle*>(p);
d731501a 204 int i=0;
dcf83226 205 if (dimension!=GetDimTHnSparse()) {
d731501a 206 // TODO: think about filling only the available data and throwing a warning
207 return -ENOSPC;
208 }
fad30eb8 209 memset(data, 0, dimension*sizeof(data[0]));
210 data[i++]=p->Pt();
211 data[i++]=p->Phi();
212 data[i++]=p->Eta();
213 if (trRP) data[i]=trRP->GetOriginMother();
214 else data[i]=fOriginMother;
215 i++; // take out of conditionals to be save
216 if (i<dimension) {
217 if(fStoreCutStepInfo) data[i]=GetLastSurvivedCutsStep();
218 i++; // take out of conditionals to be save
219 }
d731501a 220 return i;
221}
222
223int AliDxHFEParticleSelectionMCEl::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
224{
225 /// overloaded from AliDxHFEParticleSelection: check particle
226 /// H: Have changed function. Now doing particle selection first, then run MC over
227 /// selected tracks. Could configure it to be configurable, but not sure if it
228 /// is needed.
229 /// Result from normal track selection is returned, result from MC is stored in
230 /// THnSparse.
231
232 int iResult=0;
233 if (!p || !pEvent){
234 return -EINVAL;
235 }
b4779749 236 fOriginMother=-1;
237
238 if(!fUseKine && !fUseMCReco){
d731501a 239 // step 1:
240 // optional MC selection before the particle selection
241 if (fMCTools.MCFirst() && (iResult=CheckMC(p, pEvent))==0) {
242 // histograming?
243 return iResult;
244 }
245
246 // step 2 or 1, depending on sequence:
247 // normal particle selection
248 iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
249 if (fMCTools.MCFirst() || iResult==0) return iResult;
250
b4779749 251
252 }
253
d731501a 254 // step 2, only executed if MC check is last
255 // optional MC selection after the particle selection
dcf83226 256 // result stored to be filled into THnSparse
257 // TODO: strictly speaken the particles should be rejected
258 // if not mc selected, however skip this for the moment, because of
259 // the logic outside
b4779749 260 // This line will run always when running directly over the stack or over MC reconstructed tracks
261 // For MC reconstructed tracks, should consider adding more constrictions on the tracks,
262 // maybe do the track selection first (which means could call AliDxHFEParticleSelectionEl::IsSelected()
263 // and don't do PID
264
fad30eb8 265 if(fUseMCReco) {// && ( fSelectionStep > AliDxHFEParticleSelectionEl::kNoCuts )){
266 // always check the base class method in mode fUseMCReco
b4779749 267 iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
fad30eb8 268 if(iResult == 0) return iResult;
b4779749 269 }
dcf83226 270 fResultMC=CheckMC(p, pEvent);
b4779749 271
272 if(fUseKine || fUseMCReco)
273 return fResultMC;
dcf83226 274
d731501a 275 return iResult;
b4779749 276
d731501a 277}
278
279int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEvent)
280{
281 /// check if MC criteria are fulfilled
282 if (!p || !pEvent){
283 return -EINVAL;
284 }
d731501a 285 int iResult=0;
286
287 if (!fMCTools.IsInitialized() && (iResult=fMCTools.InitMCParticles(pEvent))<0) {
288 // TODO: message? but has to be filtered in order to avoid message flood
289 return 0; // no meaningful filtering on mc possible
290 }
b4779749 291
292 if(fUseKine){
293 // If run on kinematical level, check if particle is electron (only ones who are relevant)
294 // and if pass IsPhysicalPrimary()
295 Int_t test = fMCTools.CheckMCParticle(p);
296 if(test==1) return 0;
297 // Add additional contraints on the electrons?
298 // remove delta e? also remove E<300MeV?
299 // Should also mark dalitz decay and gamma conversion..
300 AliAODMCParticle *mcp=dynamic_cast<AliAODMCParticle*>(p);
1f495e82 301 if(!mcp || !mcp->IsPhysicalPrimary()) return 0;
b4779749 302
303 }
304
dfe96b90 305 int pdgParticle=-1;
b4779749 306 if (!fUseKine && fMCTools.RejectByPDG(p,false, &pdgParticle)) {
d731501a 307 // rejected by pdg
dfe96b90 308 // TODO: Move this to fMCTools???? Can this be part of the statistics in the MC class?
309 fPDGnotMCElectron->Fill(fMCTools.MapPDGLabel(pdgParticle));
d731501a 310 return 0;
311 }
312
313 int pdgMother=0;
b4779749 314 // Find PDG of first mother
d731501a 315 pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetFirstMother);
316
b4779749 317 // Check if first mother is counted as background
318 Bool_t isNotBackground=fMCTools.RejectByPDG(pdgMother,fMotherPDGs);
319
320 if(!isNotBackground){
321 // Set fOriginMother if mother is counted as background
322 // TODO: Could this be done in a more elegant way?
d731501a 323 switch(pdgMother){
324 case(AliDxHFEToolsMC::kPDGpi0): fOriginMother=AliDxHFEToolsMC::kNrOrginMother; break;
325 case(AliDxHFEToolsMC::kPDGeta): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+1; break;
326 case(AliDxHFEToolsMC::kPDGgamma): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+2;break;
327 case(AliDxHFEToolsMC::kPDGJpsi): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+3;break;
328 }
329 }
b4779749 330 else{
331 // If loop over Stack, also checks if First mother is a HF meson
332 Bool_t isHFmeson =fMCTools.TestMotherHFMeson(TMath::Abs(pdgMother));
333
334 if(isHFmeson){
335 // If first mother is a HF meson, loops back to find
336 // original quark + if there was a gluon. Result is
337 // stored in fOriginMother
338 pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetOriginMother);
339 fOriginMother=fMCTools.GetOriginMother();
340 }
341 else{
fad30eb8 342 //NotHFmother
343 fPDGNotHFMother->Fill(pdgMother);
b4779749 344 fOriginMother=AliDxHFEToolsMC::kNrOrginMother+4;
345 }
d731501a 346
b4779749 347 }
d731501a 348
349 return 1;
350}
351
352void AliDxHFEParticleSelectionMCEl::Clear(const char* option)
353{
354 /// clear internal memory
355 fMCTools.Clear(option);
356}
dfe96b90 357
358AliVParticle *AliDxHFEParticleSelectionMCEl::CreateParticle(AliVParticle* track)
359{
360
361 AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),track->Charge(),fOriginMother);
362
363 return part;
364
365}
b4779749 366
367int AliDxHFEParticleSelectionMCEl::ParseArguments(const char* arguments)
368{
369 // parse arguments and set internal flags
370 TString strArguments(arguments);
371 auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
372 if (!tokens.get()) return 0;
373
374 AliInfo(strArguments);
375 TIter next(tokens.get());
376 TObject* token;
377 while ((token=next())) {
378 TString argument=token->GetName();
379 if (argument.BeginsWith("usekine") ){
380 fUseKine=kTRUE;
381 continue;
382 }
383 if (argument.BeginsWith("elmcreco")){
384 fUseMCReco=kTRUE;
385 if(argument.BeginsWith("elmcreco=")){
386 argument.ReplaceAll("elmcreco=", "");
fad30eb8 387 if(argument.CompareTo("alltracks")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kNoCuts;
388 if(argument.CompareTo("afterreckineitstpc")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kRecKineITSTPC;
389 if(argument.CompareTo("afterrecprim")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kRecPrim;
390 if(argument.CompareTo("afterhfeits")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsITS;
391 if(argument.CompareTo("afterhfetof")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTOF;
392 if(argument.CompareTo("afterhfetpc")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTPC;
b4779749 393 if(argument.CompareTo("aftertrackcuts")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTPC;
394 if(argument.CompareTo("aftertofpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOF;
fad30eb8 395 if(argument.CompareTo("aftertpcpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTPC;
b4779749 396 if(argument.CompareTo("afterfullpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOFTPC;
397
398 AliDxHFEParticleSelectionEl::SetFinalCutStep(fSelectionStep);
b4779749 399 }
fad30eb8 400 continue;
401 }
402 if(argument.BeginsWith("storelastcutstep")){
403 AliInfo("Stores the last cut step");
404 fUseMCReco=kTRUE;
405 fStoreCutStepInfo=kTRUE;
406 AliDxHFEParticleSelectionEl::SetStoreLastCutStep(kTRUE);
b4779749 407 continue;
408 }
409 // forwarding of single argument works, unless key-option pairs separated
410 // by blanks are introduced
411 AliDxHFEParticleSelection::ParseArguments(argument);
412 }
413
414 return 0;
415}