AliPhysics  608b256 (608b256)
AliEmcalMCTrackSelector.cxx
Go to the documentation of this file.
1 /**************************************************************************************
2  * Copyright (C) 2012, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  **************************************************************************************/
28 
29 #include <iostream>
30 
31 #include <TClonesArray.h>
32 
33 #include "AliAnalysisManager.h"
34 #include "AliVEvent.h"
35 #include "AliMCEvent.h"
36 #include "AliAODMCParticle.h"
37 #include "AliMCParticle.h"
38 #include "AliStack.h"
39 #include "AliNamedArrayI.h"
40 
41 #include "AliLog.h"
42 
44 
47  fParticlesOutName("MCParticlesSelected"),
48  fOnlyPhysPrim(kTRUE),
49  fRejectNK(kFALSE),
50  fChargedMC(kFALSE),
51  fOnlyHIJING(kFALSE),
52  fRejectPhotonMothers(false),
53  fEtaMax(1),
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 }
65 
67  AliAnalysisTaskSE(name),
68  fParticlesOutName("MCParticlesSelected"),
69  fOnlyPhysPrim(kTRUE),
70  fRejectNK(kFALSE),
71  fChargedMC(kFALSE),
72  fOnlyHIJING(kFALSE),
73  fRejectPhotonMothers(false),
74  fEtaMax(1),
75  fParticlesMapName(""),
76  fInit(kFALSE),
77  fParticlesIn(0),
78  fParticlesOut(0),
79  fParticlesMap(0),
80  fEvent(0),
81  fMC(0),
82  fIsESD(kFALSE),
83  fDisabled(kFALSE)
84 {
85 }
86 
88 {
89  if (fDisabled) return;
90 
91  if (!fInit) {
92  fEvent = InputEvent();
93  if (!fEvent) {
94  AliErrorStream() << "Could not retrieve event! Returning" << std::endl;
95  return;
96  }
97 
98  if (fEvent->InheritsFrom("AliESDEvent")) fIsESD = kTRUE;
99  else fIsESD = kFALSE;
100 
101  TObject *obj = fEvent->FindListObject(fParticlesOutName);
102  if (obj) { // the output array is already present in the array!
103  AliErrorStream() << "The output array " << fParticlesOutName << " is already present in the event! Task will be disabled." << std::endl;
104  fDisabled = kTRUE;
105  return;
106  }
107  else { // copy the array from the standard ESD/AOD collections, and filter if requested
108 
109  fParticlesOut = new TClonesArray("AliAODMCParticle"); // the output will always be of AliAODMCParticle, regardless of the input
111  fEvent->AddObject(fParticlesOut);
112 
114  fParticlesMapName += "_Map";
115 
116  if (fEvent->FindListObject(fParticlesMapName)) {
117  AliErrorStream() << "The output array map " << fParticlesMapName << " is already present in the event! Task will be disabled." << std::endl;
118  fDisabled = kTRUE;
119  return;
120  }
121  else {
123  fEvent->AddObject(fParticlesMap);
124  }
125 
126  if (!fIsESD) {
127  fParticlesIn = static_cast<TClonesArray*>(InputEvent()->FindListObject(AliAODMCParticle::StdBranchName()));
128  if (!fParticlesIn) {
129  AliErrorStream() << "Could not retrieve AOD MC particles! Task will be disabled." << std::endl;
130  fDisabled = kTRUE;
131  return;
132  }
133  TClass *cl = fParticlesIn->GetClass();
134  if (!cl->GetBaseClass("AliAODMCParticle")) {
135  AliErrorStream() << GetName() << ": Collection " << AliAODMCParticle::StdBranchName() << " %s does not contain AliAODMCParticle! Task will be disabled." << std::endl;
136  fDisabled = kTRUE;
137  fParticlesIn = 0;
138  return;
139  }
140  }
141  }
142 
143  fMC = MCEvent();
144  if (!fMC) {
145  AliErrorStream() << "Could not retrieve MC event! Returning" << std::endl;
146  fDisabled = kTRUE;
147  return;
148  }
149 
150  fInit = kTRUE;
151  }
152 
155 }
156 
157 void AliEmcalMCTrackSelector::ConvertMCParticles(AliMCEvent* mcEvent, TClonesArray* partOut, AliNamedArrayI* partMap)
158 {
159 
160  // clear container (normally a null operation as the event should clean it already)
161  partOut->Delete();
162 
163  const Int_t Nparticles = mcEvent->GetNumberOfTracks();
164 
165  if (partMap) {
166  // clear particles map
167  partMap->Clear();
168  if (partMap->GetSize() <= Nparticles) partMap->Set(Nparticles*2);
169  }
170 
171  const Int_t nprim = mcEvent->GetNumberOfPrimaries();
172 
173  AliDebugStream(3) << "Number of particles: " << Nparticles << std::endl;
174 
175  // loop over particles
176  for (Int_t iPart = 0, nacc = 0; iPart < Nparticles; iPart++) {
177 
178  if (partMap) partMap->AddAt(-1, iPart);
179 
180  AliMCParticle* part = static_cast<AliMCParticle*>(mcEvent->GetTrack(iPart));
181 
182  if (!part) continue;
183 
184  Bool_t isPhysPrim = mcEvent->IsPhysicalPrimary(iPart);
185 
186  Int_t flag = 0;
187  if (iPart < nprim) flag |= AliAODMCParticle::kPrimary;
188  if (isPhysPrim) flag |= AliAODMCParticle::kPhysicalPrim;
189  if (mcEvent->IsSecondaryFromWeakDecay(iPart)) flag |= AliAODMCParticle::kSecondaryFromWeakDecay;
190  if (mcEvent->IsSecondaryFromMaterial(iPart)) flag |= AliAODMCParticle::kSecondaryFromMaterial;
191 
192  // Reject particle if is a photon and mother of another photon
193  if (fRejectPhotonMothers) {
194  if(TMath::Abs(part->PdgCode()) == 22){
195  bool hasPhotonDaughter = false;
196  // check if particle has photon daughter
197  if(part->GetDaughterFirst() > -1 && part->GetDaughterLast() > -1) {
198  for(int idaughter = part->GetDaughterFirst(); idaughter <= part->GetDaughterLast(); idaughter++) {
199  AliMCParticle *daughter = static_cast<AliMCParticle *>(mcEvent->GetTrack(idaughter));
200  if(TMath::Abs(daughter->PdgCode()) == 22) {
201  hasPhotonDaughter = true;
202  break;
203  }
204  }
205  }
206  if(hasPhotonDaughter){
207  AliDebugStream(2) << "Found photon which is the mother of another photon, not selecting ..." << std::endl;
208  continue;
209  }
210  }
211  }
212 
213  AliDebugStream(3) << "Particle " << iPart << ": pt = " << part->Pt() << ", PDG = " << part->PdgCode() <<
214  ", mother " << part->GetMother() <<
215  ", kPrimary? " << Bool_t((flag & AliAODMCParticle::kPrimary) != 0) <<
216  ", kPhysicalPrim? " << Bool_t((flag & AliAODMCParticle::kPhysicalPrim) != 0) <<
217  ", kSecondaryFromWeakDecay? " << Bool_t((flag & AliAODMCParticle::kSecondaryFromWeakDecay) != 0) <<
218  ", kSecondaryFromMaterial? " << Bool_t((flag & AliAODMCParticle::kSecondaryFromMaterial) != 0) <<
219  ", nacc = " << nacc <<
220  ", iPart = " << iPart <<
221  std::endl;
222 
223  // Do not put rejected particles on the output container
224  AliAODMCParticle parttmp(part, iPart, flag);
225  parttmp.SetGeneratorIndex(part->GetGeneratorIndex());
226  parttmp.SetStatus(part->Particle()->GetStatusCode());
227  parttmp.SetMCProcessCode(part->Particle()->GetUniqueID());
228 
229  if (!AcceptParticle(&parttmp)) continue;
230 
231  new ((*partOut)[nacc]) AliAODMCParticle(parttmp);
232  if (partMap) partMap->AddAt(nacc, iPart);
233  nacc++;
234  }
235 }
236 
237 void AliEmcalMCTrackSelector::CopyMCParticles(TClonesArray* partIn, TClonesArray* partOut, AliNamedArrayI* partMap)
238 {
239 
240  if (!partIn) return;
241  AliDebugStream(5) << "Particles in: " << partIn->GetName() << ", out: " << partOut->GetName() << std::endl;
242 
243  // clear container (normally a null operation as the event should clean it already)
244  partOut->Delete();
245 
246  const Int_t Nparticles = partIn->GetEntriesFast();
247 
248  if (partMap) {
249  // clear particles map
250  partMap->Clear();
251  if (partMap->GetSize() <= Nparticles) partMap->Set(Nparticles*2);
252  }
253 
254  AliDebugStream(2) << "Total number of particles = " << Nparticles << std::endl;
255 
256  Int_t nacc = 0;
257 
258  // loop over particles
259  for (Int_t iPart = 0; iPart < Nparticles; iPart++) {
260  if (partMap) partMap->AddAt(-1, iPart);
261 
262  AliAODMCParticle* part = static_cast<AliAODMCParticle*>(partIn->At(iPart));
263 
264  if (!AcceptParticle(part)) continue;
265 
266  // Reject particle if is a photon and mother of another photon
267  if (fRejectPhotonMothers) {
268  if(TMath::Abs(part->PdgCode()) == 22){
269  bool hasPhotonDaughter = false;
270  // check if particle has photon daughter
271  if(part->GetDaughterFirst() > -1 && part->GetDaughterLast() > -1) {
272  for(int idaughter = part->GetDaughterFirst(); idaughter <= part->GetDaughterLast(); idaughter++) {
273  AliAODMCParticle *daughter = static_cast<AliAODMCParticle *>(partIn->At(idaughter));
274  if(TMath::Abs(daughter->PdgCode()) == 22) {
275  hasPhotonDaughter = true;
276  break;
277  }
278  }
279  }
280  if(hasPhotonDaughter){
281  AliDebugStream(2) << "Found photon which is the mother of another photon, not selecting ..." << std::endl;
282  continue;
283  }
284  }
285  }
286 
287  if (partMap) partMap->AddAt(nacc, iPart);
288 
289  AliAODMCParticle *newPart = new ((*partOut)[nacc]) AliAODMCParticle(*part);
290  newPart->SetGeneratorIndex(part->GetGeneratorIndex());
291 
292  nacc++;
293  }
294  AliDebugStream(5) << "Particles in: " << partIn->GetEntries() << ", out: " << partOut->GetEntries() << std::endl;
295 }
296 
297 Bool_t AliEmcalMCTrackSelector::AcceptParticle(AliAODMCParticle* part) const
298 {
299  if (!part) return kFALSE;
300 
301  Int_t partPdgCode = TMath::Abs(part->PdgCode());
302 
303  if (fOnlyHIJING && (part->GetGeneratorIndex() != 0)) return kFALSE;
304 
305  if (fEtaMax > 0. && TMath::Abs(part->Eta()) > fEtaMax) return kFALSE;
306 
307  if (fRejectNK && (partPdgCode == 130 || partPdgCode == 2112)) return kFALSE;
308 
309  if (fChargedMC && part->Charge() == 0) return kFALSE;
310 
311  if (fOnlyPhysPrim && !part->IsPhysicalPrimary()) return kFALSE;
312 
313  return kTRUE;
314 }
315 
317 {
318  // Get the pointer to the existing analysis manager via the static access method.
319  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
320  if (!mgr) {
321  ::Error("AddTaskMCTrackSelector", "No analysis manager to connect to.");
322  return NULL;
323  }
324 
325  // Check the analysis type using the event handlers connected to the analysis manager.
326  if (!mgr->GetInputEventHandler()) {
327  ::Error("AddTaskMCTrackSelector", "This task requires an input event handler");
328  return NULL;
329  }
330 
331  // Init the task and do settings
332  TString name("AliEmcalMCTrackSelector_");
333  name += outname;
335  eTask->SetParticlesOutName(outname);
336  eTask->SetRejectNK(nk);
337  eTask->SetChargedMC(ch);
338  eTask->SetEtaMax(etamax);
339  eTask->SetOnlyPhysPrim(physPrim);
340 
341  // Final settings, pass to manager and set the containers
342  mgr->AddTask(eTask);
343 
344  // Create containers for input/output
345  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
346  mgr->ConnectInput (eTask, 0, cinput1 );
347 
348  return eTask;
349 }
void ConvertMCParticles(AliMCEvent *mcEvent, TClonesArray *partOut, AliNamedArrayI *partMap=0)
Convert MC particles in MC AOD articles (for ESD analysis).
Bool_t fOnlyHIJING
true = only HIJING particles
double Double_t
Definition: External.C:58
Bool_t fChargedMC
true = only charged particles
Double_t fEtaMax
maximum eta to accept particles
TClonesArray * fParticlesOut
! particle array out
Bool_t fIsESD
! ESD or AOD analysis
Bool_t fOnlyPhysPrim
true = only physical primary particles
void UserExec(Option_t *option)
Main event loop.
Bool_t fDisabled
! Disable task if a problem occurs at initialization
Bool_t fRejectPhotonMothers
Reject photons that are mothers of other photons.
TString fParticlesMapName
! name of the particle map
Bool_t fRejectNK
true = reject K_0^L and neutrons
void SetParticlesOutName(const char *name)
Set the name of the output container.
AliMCEvent * fMC
! MC event (ESD)
int Int_t
Definition: External.C:63
AliNamedArrayI * fParticlesMap
! particle index/label
AliEmcalMCTrackSelector()
Dummy constructor.
void CopyMCParticles(TClonesArray *partIn, TClonesArray *partOut, AliNamedArrayI *partMap=0)
Convert standard MC AOD particles in a new array, and filter if requested (for AOD analysis)...
Class to select particles in MC events. Salvatore Aiola, Yale Univeristy.
void SetOnlyPhysPrim(Bool_t s)
Select only physical primary particles.
void Clear(Option_t *option="")
const Double_t etamax
void SetEtaMax(Double_t e)
Set the eta acceptance.
void SetRejectNK(Bool_t r=kTRUE)
Reject neutrons and K0long particles.
TString fParticlesOutName
name of output particle array
const char Option_t
Definition: External.C:48
void SetChargedMC(Bool_t c=kTRUE)
Select only charged particles.
static AliEmcalMCTrackSelector * AddTaskMCTrackSelector(TString outname="mcparticles", Bool_t nk=kFALSE, Bool_t ch=kFALSE, Double_t etamax=1, Bool_t physPrim=kTRUE)
Create new AliEmcalMCTrackSelector task and add it to the analysis manager.
Bool_t fInit
! true = task initialized
bool Bool_t
Definition: External.C:53
TClonesArray * fParticlesIn
! particle array in (AOD)
virtual Bool_t AcceptParticle(AliAODMCParticle *part) const
Check whether paricle is selected.