AliPhysics  vAN-20150630 (513c479)
 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  fSpecialPDG(0),
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  fSpecialPDG(0),
54  fParticlesMapName(""),
55  fInit(kFALSE),
56  fParticlesIn(0),
57  fParticlesOut(0),
58  fParticlesMap(0),
59  fEvent(0),
60  fMC(0),
61  fIsESD(kFALSE),
62  fDisabled(kFALSE)
63 {
64  // Constructor.
65 }
66 
67 //________________________________________________________________________
69 {
70  // Destructor.
71 }
72 
73 //________________________________________________________________________
75 {
76  // Create my user objects.
77 }
78 
79 //________________________________________________________________________
81 {
82  // Main loop, called for each event.
83 
84  if (fDisabled) return;
85 
86  if (!fInit) {
87  fEvent = InputEvent();
88  if (!fEvent) {
89  AliError("Could not retrieve event! Returning");
90  return;
91  }
92 
93  if (fEvent->InheritsFrom("AliESDEvent")) fIsESD = kTRUE;
94  else fIsESD = kFALSE;
95 
96  TObject *obj = fEvent->FindListObject(fParticlesOutName);
97  if (obj) { // the output array is already present in the array!
98  AliError(Form("The output array %s is already present in the event! Task will be disabled.", fParticlesOutName.Data()));
99  fDisabled = kTRUE;
100  return;
101  }
102  else { // copy the array from the standard ESD/AOD collections, and filter if requested
103 
104  fParticlesOut = new TClonesArray("AliAODMCParticle"); // the output will always be of AliAODMCParticle, regardless of the input
106  fEvent->AddObject(fParticlesOut);
107 
109  fParticlesMapName += "_Map";
110 
111  if (fEvent->FindListObject(fParticlesMapName)) {
112  AliError(Form("The output array map %s is already present in the event! Task will be disabled.", fParticlesMapName.Data()));
113  fDisabled = kTRUE;
114  return;
115  }
116  else {
117  fParticlesMap = new AliNamedArrayI(fParticlesMapName, 99999);
118  fEvent->AddObject(fParticlesMap);
119  }
120 
121  if (!fIsESD) {
122  fParticlesIn = static_cast<TClonesArray*>(InputEvent()->FindListObject(AliAODMCParticle::StdBranchName()));
123  if (!fParticlesIn) {
124  AliError("Could not retrieve AOD MC particles! Task will be disabled.");
125  fDisabled = kTRUE;
126  return;
127  }
128  TClass *cl = fParticlesIn->GetClass();
129  if (!cl->GetBaseClass("AliAODMCParticle")) {
130  AliError(Form("%s: Collection %s does not contain AliAODMCParticle! Task will be disabled.", GetName(), AliAODMCParticle::StdBranchName()));
131  fDisabled = kTRUE;
132  fParticlesIn = 0;
133  return;
134  }
135  }
136  }
137 
138  fMC = MCEvent();
139  if (!fMC) {
140  AliError("Could not retrieve MC event! Returning");
141  fDisabled = kTRUE;
142  return;
143  }
144 
145  fInit = kTRUE;
146  }
147 
148  if (fIsESD) ConvertMCParticles();
149  else CopyMCParticles();
150 }
151 
152 //________________________________________________________________________
154 {
155  // Convert MC particles in MC AOD particles.
156 
157  // clear container (normally a null operation as the event should clean it already)
158  fParticlesOut->Delete();
159 
160  // clear particles map
161  fParticlesMap->Clear();
162 
163  const Int_t Nparticles = fMC->GetNumberOfTracks();
164  const Int_t nprim = fMC->GetNumberOfPrimaries();
165 
166  if (fParticlesMap->GetSize() <= Nparticles) fParticlesMap->Set(Nparticles*2);
167 
168  // loop over particles
169  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
170 
171  fParticlesMap->AddAt(-1, iPart);
172 
173  AliMCParticle* part = static_cast<AliMCParticle*>(fMC->GetTrack(iPart));
174 
175  if (!part) continue;
176 
177  Int_t partPdgCode = TMath::Abs(part->PdgCode());
178  Int_t genIndex = part->GetGeneratorIndex();
179 
180  if (fOnlyHIJING && genIndex != 0) continue;
181 
182  AliDebug(10, Form("Particle %d: generator index %d", iPart, part->GetGeneratorIndex()));
183 
184  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) continue;
185 
186  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) continue;
187 
188  Bool_t isSpecialPdg = (fSpecialPDG != 0 && partPdgCode == fSpecialPDG && iPart < nprim);
189 
190  if (isSpecialPdg) {
191  AliDebug(2, Form("Including particle %d (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f)",
192  iPart, partPdgCode, part->Pt(), part->Eta(), part->Phi()));
193  }
194 
195  if (fChargedMC && part->Charge() == 0 && !isSpecialPdg) continue;
196 
197  Bool_t isPhysPrim = fMC->IsPhysicalPrimary(iPart);
198  if (fOnlyPhysPrim && !isPhysPrim && !isSpecialPdg) continue;
199 
200  if (!CheckSpecialPDGDaughter(iPart)) continue;
201 
202  fParticlesMap->AddAt(nacc, iPart);
203 
204  Int_t flag = 0;
205  if (iPart < nprim) flag |= AliAODMCParticle::kPrimary;
206  if (isPhysPrim) flag |= AliAODMCParticle::kPhysicalPrim;
207  if (fMC->IsSecondaryFromWeakDecay(iPart)) flag |= AliAODMCParticle::kSecondaryFromWeakDecay;
208  if (fMC->IsSecondaryFromMaterial(iPart)) flag |= AliAODMCParticle::kSecondaryFromMaterial;
209 
210  AliAODMCParticle *aodPart = new ((*fParticlesOut)[nacc]) AliAODMCParticle(part, iPart, flag);
211  aodPart->SetGeneratorIndex(genIndex);
212  aodPart->SetStatus(part->Particle()->GetStatusCode());
213  aodPart->SetMCProcessCode(part->Particle()->GetUniqueID());
214 
215  nacc++;
216  }
217 }
218 
219 //________________________________________________________________________
221 {
222  // Convert standard MC AOD particles in a new array, and filter if requested.
223 
224  if (!fParticlesIn) return;
225 
226  // clear container (normally a null operation as the event should clean it already)
227  fParticlesOut->Delete();
228 
229  // clear particles map
230  fParticlesMap->Clear();
231 
232  const Int_t Nparticles = fParticlesIn->GetEntriesFast();
233  const Int_t nprim = fMC->GetNumberOfPrimaries();
234 
235  if (fParticlesMap->GetSize() <= Nparticles) fParticlesMap->Set(Nparticles*2);
236 
237  AliDebug(2, Form("Total number of particles = %d", Nparticles));
238 
239  // loop over particles
240  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
241 
242  fParticlesMap->AddAt(-1, iPart);
243 
244  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(fParticlesIn->At(iPart));
245 
246  if (!part) continue;
247 
248  Int_t partPdgCode = TMath::Abs(part->PdgCode());
249 
250  AliDebug(10, Form("Particle %d: generator index %d", iPart, part->GetGeneratorIndex()));
251 
252  if (fOnlyHIJING && (part->GetGeneratorIndex() != 0)) continue;
253 
254  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) continue;
255 
256  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) continue;
257 
258  Bool_t isSpecialPdg = (fSpecialPDG != 0 && partPdgCode == fSpecialPDG && iPart < nprim);
259 
260  if (isSpecialPdg && AliLog::GetDebugLevel("AliEmcalMCTrackSelector","AliEmcalMCTrackSelector") >= 2) {
261  Int_t ndaugh = part->GetNDaughters();
262  AliDebug(2, Form("Including particle %d (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f, n daughters = %d)",
263  iPart, part->PdgCode(), part->Pt(), part->Eta(), part->Phi(), ndaugh));
264 
265  for (Int_t idaugh = 0; idaugh < ndaugh; idaugh++) {
266  Int_t posDaugh = part->GetDaughter(idaugh);
267  AliAODMCParticle *daugh = static_cast<AliAODMCParticle*>(fParticlesIn->At(posDaugh));
268  if (daugh) {
269  AliDebug(2, Form("Daughter %d: i = %d, PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f",
270  idaugh, posDaugh, daugh->PdgCode(), daugh->Pt(), daugh->Eta(), daugh->Phi()));
271  }
272  }
273  }
274 
275  if (fChargedMC && part->Charge() == 0 && !isSpecialPdg) continue;
276 
277  if (fOnlyPhysPrim && !part->IsPhysicalPrimary() && !isSpecialPdg) continue;
278 
279  if (!CheckSpecialPDGDaughter(part, nprim)) continue;
280 
281  fParticlesMap->AddAt(nacc, iPart);
282 
283  AliAODMCParticle *newPart = new ((*fParticlesOut)[nacc]) AliAODMCParticle(*part);
284  newPart->SetGeneratorIndex(part->GetGeneratorIndex());
285 
286  nacc++;
287  }
288 }
289 
290 //________________________________________________________________________
291 Bool_t AliEmcalMCTrackSelector::CheckSpecialPDGDaughter(AliAODMCParticle* part, Int_t nprim)
292 {
293  // skip particle if it's a daughter of a "special" PDG particle
294  if (fSpecialPDG == 0) return kTRUE;
295 
296  AliAODMCParticle* pm = part;
297  Int_t imo = -1;
298  while (pm != 0) {
299  imo = pm->GetMother();
300  if (imo < 0) break;
301  pm = static_cast<AliAODMCParticle*>(fParticlesIn->At(imo));
302  if (TMath::Abs(pm->GetPdgCode()) == fSpecialPDG && imo < nprim) {
303  AliDebug(2, Form("Rejecting particle (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f) daughter of %d (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f)",
304  part->PdgCode(), part->Pt(), part->Eta(), part->Phi(), imo, pm->PdgCode(), pm->Pt(), pm->Eta(), pm->Phi()));
305  return kFALSE;
306  }
307  }
308  return kTRUE;
309 }
310 
311 //________________________________________________________________________
313 {
314  // skip particle if it's a daughter of a "special" PDG particle
315  if (fSpecialPDG == 0) return kTRUE;
316 
317  AliStack* stack = fMC->Stack();
318  TParticle* part = stack->Particle(iPart);
319  TParticle* pm = part;
320  Int_t imo = -1;
321  while (pm != 0) {
322  imo = pm->GetFirstMother();
323  if (imo < 0) break;
324  pm = stack->Particle(imo);
325  if (TMath::Abs(pm->GetPdgCode()) == fSpecialPDG && imo < stack->GetNprimary()) {
326  AliDebug(2, Form("Rejecting particle (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f) daughter of %d (PDG = %d, pT = %.3f, eta = %.3f, phi = %.3f)",
327  part->GetPdgCode(), part->Pt(), part->Eta(), part->Phi(), imo, pm->GetPdgCode(), pm->Pt(), pm->Eta(), pm->Phi()));
328  return kFALSE;
329  }
330  }
331  return kTRUE;
332 }
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
Bool_t CheckSpecialPDGDaughter(AliAODMCParticle *part, Int_t nprim)
AliVEvent * fEvent
particle index/label
Bool_t fInit
name of the particle map
TClonesArray * fParticlesIn
true = task initialized
ClassImp(AliEmcalMCTrackSelector) AliEmcalMCTrackSelector