AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 "AliVEvent.h"
11 #include "AliMCEvent.h"
12 #include "AliAODMCParticle.h"
13 #include "AliMCParticle.h"
14 #include "AliStack.h"
15 #include "AliNamedArrayI.h"
16 
17 #include "AliLog.h"
18 
20 
21 //________________________________________________________________________
23  AliAnalysisTaskSE("AliEmcalMCTrackSelector"),
24  fParticlesOutName("MCParticlesSelected"),
25  fOnlyPhysPrim(kTRUE),
26  fRejectNK(kFALSE),
27  fChargedMC(kFALSE),
28  fOnlyHIJING(kFALSE),
29  fEtaMax(1),
30  fParticlesMapName(""),
31  fInit(kFALSE),
32  fParticlesIn(0),
33  fParticlesOut(0),
34  fParticlesMap(0),
35  fEvent(0),
36  fMC(0),
37  fIsESD(kFALSE),
38  fDisabled(kFALSE)
39 {
40  // Constructor.
41 }
42 
43 //________________________________________________________________________
45  AliAnalysisTaskSE(name),
46  fParticlesOutName("MCParticlesSelected"),
47  fOnlyPhysPrim(kTRUE),
48  fRejectNK(kFALSE),
49  fChargedMC(kFALSE),
50  fOnlyHIJING(kFALSE),
51  fEtaMax(1),
52  fParticlesMapName(""),
53  fInit(kFALSE),
54  fParticlesIn(0),
55  fParticlesOut(0),
56  fParticlesMap(0),
57  fEvent(0),
58  fMC(0),
59  fIsESD(kFALSE),
60  fDisabled(kFALSE)
61 {
62  // Constructor.
63 }
64 
65 //________________________________________________________________________
67 {
68  // Destructor.
69 }
70 
71 //________________________________________________________________________
73 {
74  // Create my user objects.
75 }
76 
77 //________________________________________________________________________
79 {
80  // Main loop, called for each event.
81 
82  if (fDisabled) return;
83 
84  if (!fInit) {
85  fEvent = InputEvent();
86  if (!fEvent) {
87  AliError("Could not retrieve event! Returning");
88  return;
89  }
90 
91  if (fEvent->InheritsFrom("AliESDEvent")) fIsESD = kTRUE;
92  else fIsESD = kFALSE;
93 
94  TObject *obj = fEvent->FindListObject(fParticlesOutName);
95  if (obj) { // the output array is already present in the array!
96  AliError(Form("The output array %s is already present in the event! Task will be disabled.", fParticlesOutName.Data()));
97  fDisabled = kTRUE;
98  return;
99  }
100  else { // copy the array from the standard ESD/AOD collections, and filter if requested
101 
102  fParticlesOut = new TClonesArray("AliAODMCParticle"); // the output will always be of AliAODMCParticle, regardless of the input
104  fEvent->AddObject(fParticlesOut);
105 
107  fParticlesMapName += "_Map";
108 
109  if (fEvent->FindListObject(fParticlesMapName)) {
110  AliError(Form("The output array map %s is already present in the event! Task will be disabled.", fParticlesMapName.Data()));
111  fDisabled = kTRUE;
112  return;
113  }
114  else {
115  fParticlesMap = new AliNamedArrayI(fParticlesMapName, 99999);
116  fEvent->AddObject(fParticlesMap);
117  }
118 
119  if (!fIsESD) {
120  fParticlesIn = static_cast<TClonesArray*>(InputEvent()->FindListObject(AliAODMCParticle::StdBranchName()));
121  if (!fParticlesIn) {
122  AliError("Could not retrieve AOD MC particles! Task will be disabled.");
123  fDisabled = kTRUE;
124  return;
125  }
126  TClass *cl = fParticlesIn->GetClass();
127  if (!cl->GetBaseClass("AliAODMCParticle")) {
128  AliError(Form("%s: Collection %s does not contain AliAODMCParticle! Task will be disabled.", GetName(), AliAODMCParticle::StdBranchName()));
129  fDisabled = kTRUE;
130  fParticlesIn = 0;
131  return;
132  }
133  }
134  }
135 
136  fMC = MCEvent();
137  if (!fMC) {
138  AliError("Could not retrieve MC event! Returning");
139  fDisabled = kTRUE;
140  return;
141  }
142 
143  fInit = kTRUE;
144  }
145 
146  if (fIsESD) ConvertMCParticles();
147  else CopyMCParticles();
148 }
149 
150 //________________________________________________________________________
152 {
153  // Convert MC particles in MC AOD articles.
154 
155  // clear container (normally a null operation as the event should clean it already)
156  fParticlesOut->Delete();
157 
158  // clear particles map
159  fParticlesMap->Clear();
160 
161  const Int_t Nparticles = fMC->GetNumberOfTracks();
162  const Int_t nprim = fMC->GetNumberOfPrimaries();
163 
164  if (fParticlesMap->GetSize() <= Nparticles) fParticlesMap->Set(Nparticles*2);
165 
166  // loop over particles
167  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
168 
169  fParticlesMap->AddAt(-1, iPart);
170 
171  AliMCParticle* part = static_cast<AliMCParticle*>(fMC->GetTrack(iPart));
172 
173  if (!part) continue;
174 
175  Bool_t isPhysPrim = fMC->IsPhysicalPrimary(iPart);
176 
177  Int_t flag = 0;
178  if (iPart < nprim) flag |= AliAODMCParticle::kPrimary;
179  if (isPhysPrim) flag |= AliAODMCParticle::kPhysicalPrim;
180  if (fMC->IsSecondaryFromWeakDecay(iPart)) flag |= AliAODMCParticle::kSecondaryFromWeakDecay;
181  if (fMC->IsSecondaryFromMaterial(iPart)) flag |= AliAODMCParticle::kSecondaryFromMaterial;
182 
183  AliAODMCParticle *aodPart = new ((*fParticlesOut)[nacc]) AliAODMCParticle(part, iPart, flag);
184  aodPart->SetGeneratorIndex(part->GetGeneratorIndex());
185  aodPart->SetStatus(part->Particle()->GetStatusCode());
186  aodPart->SetMCProcessCode(part->Particle()->GetUniqueID());
187 
188  if (!AcceptParticle(aodPart)) continue;
189 
190  fParticlesMap->AddAt(nacc, iPart);
191  nacc++;
192  }
193 }
194 
195 //________________________________________________________________________
197 {
198  // Convert standard MC AOD particles in a new array, and filter if requested.
199 
200  if (!fParticlesIn) return;
201 
202  // clear container (normally a null operation as the event should clean it already)
203  fParticlesOut->Delete();
204 
205  // clear particles map
206  fParticlesMap->Clear();
207 
208  const Int_t Nparticles = fParticlesIn->GetEntriesFast();
209 
210  if (fParticlesMap->GetSize() <= Nparticles) fParticlesMap->Set(Nparticles*2);
211 
212  AliDebug(2, Form("Total number of particles = %d", Nparticles));
213 
214  Int_t nacc = 0;
215 
216  // loop over particles
217  for (Int_t iPart = 0; iPart < Nparticles; iPart++) {
218  fParticlesMap->AddAt(-1, iPart);
219 
220  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(fParticlesIn->At(iPart));
221 
222  if (!AcceptParticle(part)) continue;
223 
224  fParticlesMap->AddAt(nacc, iPart);
225 
226  AliAODMCParticle *newPart = new ((*fParticlesOut)[nacc]) AliAODMCParticle(*part);
227  newPart->SetGeneratorIndex(part->GetGeneratorIndex());
228 
229  nacc++;
230  }
231 }
232 
233 //________________________________________________________________________
234 Bool_t AliEmcalMCTrackSelector::AcceptParticle(AliAODMCParticle* part) const
235 {
236  // Determine whether the MC particle is accepted.
237 
238  if (!part) return kFALSE;
239 
240  Int_t partPdgCode = TMath::Abs(part->PdgCode());
241 
242  if (fOnlyHIJING && (part->GetGeneratorIndex() != 0)) return kFALSE;
243 
244  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) return kFALSE;
245 
246  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) return kFALSE;
247 
248  if (fChargedMC && part->Charge() == 0) return kFALSE;
249 
250  if (fOnlyPhysPrim && !part->IsPhysicalPrimary()) return kFALSE;
251 
252  return kTRUE;
253 }
TClonesArray * fParticlesOut
particle array in (AOD)
Bool_t fIsESD
MC event (ESD)
void UserExec(Option_t *option)
Bool_t fDisabled
ESD or AOD analysis.
AliNamedArrayI * fParticlesMap
particle array out
AliVEvent * fEvent
particle index/label
Bool_t fInit
name of the particle map
TClonesArray * fParticlesIn
true = task initialized
ClassImp(AliEmcalMCTrackSelector) AliEmcalMCTrackSelector
virtual Bool_t AcceptParticle(AliAODMCParticle *part) const