Example macro for filtering MC info to the AOD on proof, added comment
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.cxx
CommitLineData
da97a08a 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//
17// Analysis task for Kinematic filtering
18// Fill AODtrackMC tracks from Kinematic stack
19//
20
21#include <TChain.h>
22#include <TFile.h>
23#include "TParticle.h"
24#include "TString.h"
25
26#include "AliAnalysisTaskMCParticleFilter.h"
27#include "AliAnalysisManager.h"
28#include "AliAnalysisFilter.h"
29#include "AliHeader.h"
30#include "AliStack.h"
31#include "AliMCEvent.h"
32#include "AliMCEventHandler.h"
33#include "AliAODEvent.h"
34#include "AliAODHeader.h"
35#include "AliAODHandler.h"
36#include "AliAODVertex.h"
37#include "AliAODMCParticle.h"
38
39
40#include "AliLog.h"
41
42ClassImp(AliAnalysisTaskMCParticleFilter)
43
44////////////////////////////////////////////////////////////////////////
45
46//____________________________________________________________________
47AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
48AliAnalysisTaskSE(),
49 fTrackFilterMother(0x0)
50{
51 // Default constructor
52}
53
54//____________________________________________________________________
55AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
56 AliAnalysisTaskSE(name),
57 fTrackFilterMother(0x0)
58{
59 // Default constructor
60}
61
62
63//____________________________________________________________________
64AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
65 AliAnalysisTaskSE(obj),
66 fTrackFilterMother(obj.fTrackFilterMother)
67{
68 // Copy constructor
69}
70
71//____________________________________________________________________
72AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
73{
74 // if( fTrackFilterMother ) delete fTrackFilterMother;
75}
76
77
78//____________________________________________________________________
79AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
80{
81// Assignment
82 if(this!=&other) {
83 AliAnalysisTaskSE::operator=(other);
84 fTrackFilterMother = other.fTrackFilterMother;
85 }
86 return *this;
87}
88
89//____________________________________________________________________
90void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
91{
92 // Create the output container
93 if (OutputTree()&&fTrackFilterMother)
94 OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba 95
96 // this part is mainly needed to set the MCEventHandler
97 // to the AODHandler, this will not be needed when
98 // AODHandler goes to ANALYSISalice
99 // setting in the steering macro will not work on proof :(
100 // so we have to do it in a task
101
102 // the branch booking can also go into the AODHandler then
103
da97a08a 104 TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
105 tca->SetName(AliAODMCParticle::StdBranchName());
106 AddAODBranch("TClonesArray",&tca);
107 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
108
109 AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
110 if(!aodH){
111 Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
112 return;
113 }
114 aodH->SetMCEventHandler(mcH);
115 // TODO ADD MC VERTEX
116
117}
118
119//____________________________________________________________________
120void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
121{
122// Execute analysis for current event
123//
124
125// Fill AOD tracks from Kinematic stack
126
127 // get AliAOD Event
128 AliAODEvent* aod = AODEvent();
129 if (!aod) {
130 AliWarning("No Output Handler connected, doing nothing !") ;
131 return;
132 }
133
134 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
135 if(!mcH){
136 AliWarning("No MC handler Found");
137 return;
138 }
139
140
141 // fetch the output
142 // Fill control histos
143
144 // tmp array for holding the mctracks
145 // need to set mother and duagther __before__
146 // filling in the tree
147
148 AliMCEvent *mcE = MCEvent();
149
150 if(!mcE){
151 AliWarning("No MC event Found");
152 return;
153 }
154
155 AliStack* stack = mcE->Stack();
156 Int_t np = mcE->GetNumberOfTracks();
157 Int_t nprim = stack->GetNprimary();
158 // TODO ADD MC VERTEX
159
160 // We take all real primaries
161
162
163
164 Int_t j=0;
165 for (Int_t ip = 0; ip < np; ip++){
166 AliMCParticle* mcpart = mcE->GetTrack(ip);
167 TParticle* part = mcpart->Particle();
168
169 Float_t xv = part->Vx();
170 Float_t yv = part->Vy();
171 Float_t zv = part->Vz();
172 Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
173
174 Bool_t write = kFALSE;
175 Int_t flag = 0;
176
177
178 if (ip < nprim) {
179 // Select the primary event
180 write = kTRUE;
181 } else if (part->GetUniqueID() == 4) {
182 // Particles from decay
183 // Check that the decay chain ends at a primary particle
184 TParticle* mother = part;
185 Int_t imo = part->GetFirstMother();
186 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
187 mother = mcE->Stack()->Particle(imo);
188 imo = mother->GetFirstMother();
189 }
190 // Select according to pseudorapidity and production point
191 if (imo < nprim && Select(mother, rv, zv))
192 write = kTRUE;
193 } else if (part->GetUniqueID() == 5) {
194 // Now look for pair production
195 continue;
196 Int_t imo = part->GetFirstMother();
197 if (imo < nprim) {
198 // Select, if the gamma is a primary
199 write = kTRUE;
200 } else {
201 // Check if the gamma comes from the decay chain of a primary particle
202 TParticle* mother = mcE->Stack()->Particle(imo);
203 imo = mother->GetFirstMother();
204 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
205 mother = mcE->Stack()->Particle(imo);
206 imo = mother->GetFirstMother();
207 }
208 // Select according to pseudorapidity and production point
209 if (imo < nprim && Select(mother, rv, zv))
210 write = kTRUE;
211 }
212 }
213 /*
214 else if (part->GetUniqueID() == 13){
215 // Evaporation
216 // Check that we end at a primary particle
217 TParticle* mother = part;
218 Int_t imo = part->GetFirstMother();
219 while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
220 mother = mcE->Stack()->Particle(imo);
221 imo = mother->GetFirstMother();
222 }
223 // Select according to pseudorapidity and production point
224 if (imo < nprim && Select(mother, rv, zv))
225 write = kTRUE;
226 }
227 */
228 if (write) {
229 if(mcH)mcH->SelectParticle(ip);
230 j++;
231 }
232 }
233
234 return;
235}
236
237Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
238{
239 // Selection accoring to eta of the mother and production point
240 // This has to be refined !!!!!!
241
242 // Esp if we don't have collisison in the central barrel but e.g. beam gas
243 // return kTRUE;
244
245 Float_t eta = part->Eta();
246
247 // central barrel consider everything in the ITS...
248 // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
249 // larger for V0s in the TPC
250 // if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
251
252 // Andreas' Cuts
253 if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
254
255 // Muon arm
256 if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
257
258 // PMD acceptance 2.3 <= eta < = 3.5
259 // if(eta>2.0&&eta<3.7)return kTRUE;
260
261 return kFALSE;
262
263}
264