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