AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalMCTrackSelector.cxx
Go to the documentation of this file.
1 //
2 // Class to select particles in MC events.
3 //
4 // Author: S. Aiola
5 
7 
8 #include <TClonesArray.h>
9 
10 #include "AliAnalysisManager.h"
11 #include "AliVEvent.h"
12 #include "AliMCEvent.h"
13 #include "AliAODMCParticle.h"
14 #include "AliMCParticle.h"
15 #include "AliStack.h"
16 #include "AliNamedArrayI.h"
17 
18 #include "AliLog.h"
19 
21 
22 //________________________________________________________________________
25  fParticlesOutName("MCParticlesSelected"),
26  fOnlyPhysPrim(kTRUE),
27  fRejectNK(kFALSE),
28  fChargedMC(kFALSE),
29  fOnlyHIJING(kFALSE),
30  fEtaMax(1),
31  fParticlesMapName(""),
32  fInit(kFALSE),
33  fParticlesIn(0),
34  fParticlesOut(0),
35  fParticlesMap(0),
36  fEvent(0),
37  fMC(0),
38  fIsESD(kFALSE),
39  fDisabled(kFALSE)
40 {
41  // Constructor.
42 }
43 
44 //________________________________________________________________________
46  AliAnalysisTaskSE(name),
47  fParticlesOutName("MCParticlesSelected"),
48  fOnlyPhysPrim(kTRUE),
49  fRejectNK(kFALSE),
50  fChargedMC(kFALSE),
51  fOnlyHIJING(kFALSE),
52  fEtaMax(1),
53  fParticlesMapName(""),
54  fInit(kFALSE),
55  fParticlesIn(0),
56  fParticlesOut(0),
57  fParticlesMap(0),
58  fEvent(0),
59  fMC(0),
60  fIsESD(kFALSE),
61  fDisabled(kFALSE)
62 {
63  // Constructor.
64 }
65 
66 //________________________________________________________________________
68 {
69  // Destructor.
70 }
71 
72 //________________________________________________________________________
74 {
75  // Create my user objects.
76 }
77 
78 //________________________________________________________________________
80 {
81  // Main loop, called for each event.
82 
83  if (fDisabled) return;
84 
85  if (!fInit) {
86  fEvent = InputEvent();
87  if (!fEvent) {
88  AliError("Could not retrieve event! Returning");
89  return;
90  }
91 
92  if (fEvent->InheritsFrom("AliESDEvent")) fIsESD = kTRUE;
93  else fIsESD = kFALSE;
94 
95  TObject *obj = fEvent->FindListObject(fParticlesOutName);
96  if (obj) { // the output array is already present in the array!
97  AliError(Form("The output array %s is already present in the event! Task will be disabled.", fParticlesOutName.Data()));
98  fDisabled = kTRUE;
99  return;
100  }
101  else { // copy the array from the standard ESD/AOD collections, and filter if requested
102 
103  fParticlesOut = new TClonesArray("AliAODMCParticle"); // the output will always be of AliAODMCParticle, regardless of the input
105  fEvent->AddObject(fParticlesOut);
106 
108  fParticlesMapName += "_Map";
109 
110  if (fEvent->FindListObject(fParticlesMapName)) {
111  AliError(Form("The output array map %s is already present in the event! Task will be disabled.", fParticlesMapName.Data()));
112  fDisabled = kTRUE;
113  return;
114  }
115  else {
117  fEvent->AddObject(fParticlesMap);
118  }
119 
120  if (!fIsESD) {
121  fParticlesIn = static_cast<TClonesArray*>(InputEvent()->FindListObject(AliAODMCParticle::StdBranchName()));
122  if (!fParticlesIn) {
123  AliError("Could not retrieve AOD MC particles! Task will be disabled.");
124  fDisabled = kTRUE;
125  return;
126  }
127  TClass *cl = fParticlesIn->GetClass();
128  if (!cl->GetBaseClass("AliAODMCParticle")) {
129  AliError(Form("%s: Collection %s does not contain AliAODMCParticle! Task will be disabled.", GetName(), AliAODMCParticle::StdBranchName()));
130  fDisabled = kTRUE;
131  fParticlesIn = 0;
132  return;
133  }
134  }
135  }
136 
137  fMC = MCEvent();
138  if (!fMC) {
139  AliError("Could not retrieve MC event! Returning");
140  fDisabled = kTRUE;
141  return;
142  }
143 
144  fInit = kTRUE;
145  }
146 
149 }
150 
151 //________________________________________________________________________
152 void AliEmcalMCTrackSelector::ConvertMCParticles(AliMCEvent* mcEvent, TClonesArray* partOut, AliNamedArrayI* partMap)
153 {
154  // Convert MC particles in MC AOD articles.
155 
156  // clear container (normally a null operation as the event should clean it already)
157  partOut->Delete();
158 
159  const Int_t Nparticles = mcEvent->GetNumberOfTracks();
160 
161  if (partMap) {
162  // clear particles map
163  partMap->Clear();
164  if (partMap->GetSize() <= Nparticles) partMap->Set(Nparticles*2);
165  }
166 
167  const Int_t nprim = mcEvent->GetNumberOfPrimaries();
168 
169  // loop over particles
170  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
171 
172  if (partMap) partMap->AddAt(-1, iPart);
173 
174  AliMCParticle* part = static_cast<AliMCParticle*>(mcEvent->GetTrack(iPart));
175 
176  if (!part) continue;
177 
178  Bool_t isPhysPrim = mcEvent->IsPhysicalPrimary(iPart);
179 
180  Int_t flag = 0;
181  if (iPart < nprim) flag |= AliAODMCParticle::kPrimary;
182  if (isPhysPrim) flag |= AliAODMCParticle::kPhysicalPrim;
183  if (mcEvent->IsSecondaryFromWeakDecay(iPart)) flag |= AliAODMCParticle::kSecondaryFromWeakDecay;
184  if (mcEvent->IsSecondaryFromMaterial(iPart)) flag |= AliAODMCParticle::kSecondaryFromMaterial;
185 
186  AliAODMCParticle *aodPart = new ((*partOut)[nacc]) AliAODMCParticle(part, iPart, flag);
187  aodPart->SetGeneratorIndex(part->GetGeneratorIndex());
188  aodPart->SetStatus(part->Particle()->GetStatusCode());
189  aodPart->SetMCProcessCode(part->Particle()->GetUniqueID());
190 
191  if (!AcceptParticle(aodPart)) continue;
192 
193  if (partMap) partMap->AddAt(nacc, iPart);
194  nacc++;
195  }
196 }
197 
198 //________________________________________________________________________
199 void AliEmcalMCTrackSelector::CopyMCParticles(TClonesArray* partIn, TClonesArray* partOut, AliNamedArrayI* partMap)
200 {
201  // Convert standard MC AOD particles in a new array, and filter if requested.
202 
203  if (!partIn) return;
204 
205  // clear container (normally a null operation as the event should clean it already)
206  partOut->Delete();
207 
208  const Int_t Nparticles = partIn->GetEntriesFast();
209 
210  if (partMap) {
211  // clear particles map
212  partMap->Clear();
213  if (partMap->GetSize() <= Nparticles) partMap->Set(Nparticles*2);
214  }
215 
216  AliDebug(2, Form("Total number of particles = %d", Nparticles));
217 
218  Int_t nacc = 0;
219 
220  // loop over particles
221  for (Int_t iPart = 0; iPart < Nparticles; iPart++) {
222  if (partMap) partMap->AddAt(-1, iPart);
223 
224  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(partIn->At(iPart));
225 
226  if (!AcceptParticle(part)) continue;
227 
228  if (partMap) partMap->AddAt(nacc, iPart);
229 
230  AliAODMCParticle *newPart = new ((*partOut)[nacc]) AliAODMCParticle(*part);
231  newPart->SetGeneratorIndex(part->GetGeneratorIndex());
232 
233  nacc++;
234  }
235 }
236 
237 //________________________________________________________________________
238 Bool_t AliEmcalMCTrackSelector::AcceptParticle(AliAODMCParticle* part) const
239 {
240  // Determine whether the MC particle is accepted.
241 
242  if (!part) return kFALSE;
243 
244  Int_t partPdgCode = TMath::Abs(part->PdgCode());
245 
246  if (fOnlyHIJING && (part->GetGeneratorIndex() != 0)) return kFALSE;
247 
248  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) return kFALSE;
249 
250  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) return kFALSE;
251 
252  if (fChargedMC && part->Charge() == 0) return kFALSE;
253 
254  if (fOnlyPhysPrim && !part->IsPhysicalPrimary()) return kFALSE;
255 
256  return kTRUE;
257 }
258 
259 //________________________________________________________________________
261 {
262  // Get the pointer to the existing analysis manager via the static access method.
263  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
264  if (!mgr) {
265  ::Error("AddTaskMCTrackSelector", "No analysis manager to connect to.");
266  return NULL;
267  }
268 
269  // Check the analysis type using the event handlers connected to the analysis manager.
270  if (!mgr->GetInputEventHandler()) {
271  ::Error("AddTaskMCTrackSelector", "This task requires an input event handler");
272  return NULL;
273  }
274 
275  // Init the task and do settings
276  TString name("AliEmcalMCTrackSelector_");
277  name += outname;
279  eTask->SetParticlesOutName(outname);
280  eTask->SetRejectNK(nk);
281  eTask->SetChargedMC(ch);
282  eTask->SetEtaMax(etamax);
283  eTask->SetOnlyPhysPrim(physPrim);
284 
285  // Final settings, pass to manager and set the containers
286  mgr->AddTask(eTask);
287 
288  // Create containers for input/output
289  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
290  mgr->ConnectInput (eTask, 0, cinput1 );
291 
292  return eTask;
293 }
void ConvertMCParticles(AliMCEvent *mcEvent, TClonesArray *partOut, AliNamedArrayI *partMap=0)
double Double_t
Definition: External.C:58
TClonesArray * fParticlesOut
particle array in (AOD)
Bool_t fIsESD
MC event (ESD)
void UserExec(Option_t *option)
Bool_t fDisabled
ESD or AOD analysis.
void SetParticlesOutName(const char *name)
int Int_t
Definition: External.C:63
AliNamedArrayI * fParticlesMap
particle array out
void CopyMCParticles(TClonesArray *partIn, TClonesArray *partOut, AliNamedArrayI *partMap=0)
void Clear(Option_t *option="")
AliVEvent * fEvent
particle index/label
const Double_t etamax
void SetRejectNK(Bool_t r=kTRUE)
const char Option_t
Definition: External.C:48
void SetChargedMC(Bool_t c=kTRUE)
static AliEmcalMCTrackSelector * AddTaskMCTrackSelector(TString outname="mcparticles", Bool_t nk=kFALSE, Bool_t ch=kFALSE, Double_t etamax=1, Bool_t physPrim=kTRUE)
Bool_t fInit
name of the particle map
bool Bool_t
Definition: External.C:53
TClonesArray * fParticlesIn
true = task initialized
ClassImp(AliEmcalMCTrackSelector) AliEmcalMCTrackSelector
virtual Bool_t AcceptParticle(AliAODMCParticle *part) const