AliPhysics  d565ceb (d565ceb)
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  // loop over particles
172  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
173 
174  if (partMap) partMap->AddAt(-1, iPart);
175 
176  AliMCParticle* part = static_cast<AliMCParticle*>(mcEvent->GetTrack(iPart));
177 
178  if (!part) continue;
179 
180  Bool_t isPhysPrim = mcEvent->IsPhysicalPrimary(iPart);
181 
182  Int_t flag = 0;
183  if (iPart < nprim) flag |= AliAODMCParticle::kPrimary;
184  if (isPhysPrim) flag |= AliAODMCParticle::kPhysicalPrim;
185  if (mcEvent->IsSecondaryFromWeakDecay(iPart)) flag |= AliAODMCParticle::kSecondaryFromWeakDecay;
186  if (mcEvent->IsSecondaryFromMaterial(iPart)) flag |= AliAODMCParticle::kSecondaryFromMaterial;
187 
188  AliDebugStream(3) << "Particle " << iPart << ": pt = " << part->Pt() << ", PDG = " << part->PdgCode() <<
189  ", mother " << part->GetMother() <<
190  ", kPrimary? " << Bool_t((flag & AliAODMCParticle::kPrimary) != 0) <<
191  ", kPhysicalPrim? " << Bool_t((flag & AliAODMCParticle::kPhysicalPrim) != 0) <<
192  ", kSecondaryFromWeakDecay? " << Bool_t((flag & AliAODMCParticle::kSecondaryFromWeakDecay) != 0) <<
193  ", kSecondaryFromMaterial? " << Bool_t((flag & AliAODMCParticle::kSecondaryFromMaterial) != 0) <<
194  std::endl;
195 
196  AliAODMCParticle *aodPart = new ((*partOut)[nacc]) AliAODMCParticle(part, iPart, flag);
197  aodPart->SetGeneratorIndex(part->GetGeneratorIndex());
198  aodPart->SetStatus(part->Particle()->GetStatusCode());
199  aodPart->SetMCProcessCode(part->Particle()->GetUniqueID());
200 
201  if (!AcceptParticle(aodPart)) continue;
202 
203  if (partMap) partMap->AddAt(nacc, iPart);
204  nacc++;
205  }
206 }
207 
208 //________________________________________________________________________
209 void AliEmcalMCTrackSelector::CopyMCParticles(TClonesArray* partIn, TClonesArray* partOut, AliNamedArrayI* partMap)
210 {
211  // Convert standard MC AOD particles in a new array, and filter if requested.
212 
213  if (!partIn) return;
214 
215  // clear container (normally a null operation as the event should clean it already)
216  partOut->Delete();
217 
218  const Int_t Nparticles = partIn->GetEntriesFast();
219 
220  if (partMap) {
221  // clear particles map
222  partMap->Clear();
223  if (partMap->GetSize() <= Nparticles) partMap->Set(Nparticles*2);
224  }
225 
226  AliDebug(2, Form("Total number of particles = %d", Nparticles));
227 
228  Int_t nacc = 0;
229 
230  // loop over particles
231  for (Int_t iPart = 0; iPart < Nparticles; iPart++) {
232  if (partMap) partMap->AddAt(-1, iPart);
233 
234  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(partIn->At(iPart));
235 
236  if (!AcceptParticle(part)) continue;
237 
238  if (partMap) partMap->AddAt(nacc, iPart);
239 
240  AliAODMCParticle *newPart = new ((*partOut)[nacc]) AliAODMCParticle(*part);
241  newPart->SetGeneratorIndex(part->GetGeneratorIndex());
242 
243  nacc++;
244  }
245 }
246 
247 //________________________________________________________________________
249 {
250  // Determine whether the MC particle is accepted.
251 
252  if (!part) return kFALSE;
253 
254  Int_t partPdgCode = TMath::Abs(part->PdgCode());
255 
256  if (fOnlyHIJING && (part->GetGeneratorIndex() != 0)) return kFALSE;
257 
258  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) return kFALSE;
259 
260  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) return kFALSE;
261 
262  if (fChargedMC && part->Charge() == 0) return kFALSE;
263 
264  if (fOnlyPhysPrim && !part->IsPhysicalPrimary()) return kFALSE;
265 
266  return kTRUE;
267 }
268 
269 //________________________________________________________________________
271 {
272  // Get the pointer to the existing analysis manager via the static access method.
273  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
274  if (!mgr) {
275  ::Error("AddTaskMCTrackSelector", "No analysis manager to connect to.");
276  return NULL;
277  }
278 
279  // Check the analysis type using the event handlers connected to the analysis manager.
280  if (!mgr->GetInputEventHandler()) {
281  ::Error("AddTaskMCTrackSelector", "This task requires an input event handler");
282  return NULL;
283  }
284 
285  // Init the task and do settings
286  TString name("AliEmcalMCTrackSelector_");
287  name += outname;
289  eTask->SetParticlesOutName(outname);
290  eTask->SetRejectNK(nk);
291  eTask->SetChargedMC(ch);
292  eTask->SetEtaMax(etamax);
293  eTask->SetOnlyPhysPrim(physPrim);
294 
295  // Final settings, pass to manager and set the containers
296  mgr->AddTask(eTask);
297 
298  // Create containers for input/output
299  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
300  mgr->ConnectInput (eTask, 0, cinput1 );
301 
302  return eTask;
303 }
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.
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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