CMake target for TENDER.
[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"
dce1b636 35#include "AliAODMCHeader.h"
da97a08a 36#include "AliAODHandler.h"
37#include "AliAODVertex.h"
38#include "AliAODMCParticle.h"
39
40
41#include "AliLog.h"
42
43ClassImp(AliAnalysisTaskMCParticleFilter)
44
45////////////////////////////////////////////////////////////////////////
46
47//____________________________________________________________________
48AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
49AliAnalysisTaskSE(),
50 fTrackFilterMother(0x0)
51{
52 // Default constructor
53}
54
55//____________________________________________________________________
56AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
57 AliAnalysisTaskSE(name),
58 fTrackFilterMother(0x0)
59{
60 // Default constructor
61}
62
63
64//____________________________________________________________________
65AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
66 AliAnalysisTaskSE(obj),
67 fTrackFilterMother(obj.fTrackFilterMother)
68{
69 // Copy constructor
70}
71
72//____________________________________________________________________
73AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
74{
75 // if( fTrackFilterMother ) delete fTrackFilterMother;
76}
77
78
79//____________________________________________________________________
80AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
81{
82// Assignment
83 if(this!=&other) {
84 AliAnalysisTaskSE::operator=(other);
85 fTrackFilterMother = other.fTrackFilterMother;
86 }
87 return *this;
88}
89
90//____________________________________________________________________
91void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
92{
93 // Create the output container
94 if (OutputTree()&&fTrackFilterMother)
95 OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba 96
97 // this part is mainly needed to set the MCEventHandler
98 // to the AODHandler, this will not be needed when
99 // AODHandler goes to ANALYSISalice
100 // setting in the steering macro will not work on proof :(
101 // so we have to do it in a task
102
103 // the branch booking can also go into the AODHandler then
104
dce1b636 105
106 // mcparticles
da97a08a 107 TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
108 tca->SetName(AliAODMCParticle::StdBranchName());
109 AddAODBranch("TClonesArray",&tca);
da97a08a 110
dce1b636 111 // MC header...
112 AliAODMCHeader *mcHeader = new AliAODMCHeader();
dce1b636 113 mcHeader->SetName(AliAODMCHeader::StdBranchName());
114 AddAODBranch("AliAODMCHeader",&mcHeader);
115
116
117
118 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
da97a08a 119 AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
120 if(!aodH){
121 Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
122 return;
123 }
124 aodH->SetMCEventHandler(mcH);
dce1b636 125
da97a08a 126
127}
128
129//____________________________________________________________________
130void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
131{
132// Execute analysis for current event
133//
134
135// Fill AOD tracks from Kinematic stack
136
137 // get AliAOD Event
138 AliAODEvent* aod = AODEvent();
139 if (!aod) {
140 AliWarning("No Output Handler connected, doing nothing !") ;
141 return;
142 }
143
144 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
145 if(!mcH){
146 AliWarning("No MC handler Found");
147 return;
148 }
149
150
151 // fetch the output
152 // Fill control histos
153
154 // tmp array for holding the mctracks
155 // need to set mother and duagther __before__
156 // filling in the tree
157
158 AliMCEvent *mcE = MCEvent();
159
160 if(!mcE){
161 AliWarning("No MC event Found");
162 return;
163 }
164
da97a08a 165 Int_t np = mcE->GetNumberOfTracks();
93836e1b 166 Int_t nprim = mcE->GetNumberOfPrimaries();
da97a08a 167 // TODO ADD MC VERTEX
168
169 // We take all real primaries
170
171
da97a08a 172 Int_t j=0;
173 for (Int_t ip = 0; ip < np; ip++){
174 AliMCParticle* mcpart = mcE->GetTrack(ip);
175 TParticle* part = mcpart->Particle();
da97a08a 176 Float_t xv = part->Vx();
177 Float_t yv = part->Vy();
178 Float_t zv = part->Vz();
179 Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
180
181 Bool_t write = kFALSE;
63892dd0 182
da97a08a 183 if (ip < nprim) {
184 // Select the primary event
185 write = kTRUE;
186 } else if (part->GetUniqueID() == 4) {
187 // Particles from decay
188 // Check that the decay chain ends at a primary particle
93836e1b 189 AliMCParticle* mother = mcpart;
190 Int_t imo = mcpart->GetMother();
da97a08a 191 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
93836e1b 192 mother = mcE->GetTrack(imo);
193 imo = mother->GetMother();
da97a08a 194 }
63892dd0 195 // Select according to pseudorapidity and production point of primary ancestor
93836e1b 196 if (imo < nprim && Select(mcE->GetTrack(imo)->Particle(), rv, zv))write = kTRUE;
da97a08a 197 } else if (part->GetUniqueID() == 5) {
198 // Now look for pair production
93836e1b 199 Int_t imo = mcpart->GetMother();
da97a08a 200 if (imo < nprim) {
201 // Select, if the gamma is a primary
202 write = kTRUE;
203 } else {
204 // Check if the gamma comes from the decay chain of a primary particle
93836e1b 205 AliMCParticle* mother = mcE->GetTrack(imo);
206 imo = mother->GetMother();
da97a08a 207 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
93836e1b 208 mother = mcE->GetTrack(imo);
209 imo = mother->GetMother();
da97a08a 210 }
211 // Select according to pseudorapidity and production point
93836e1b 212 if (imo < nprim && Select(mother->Particle(), rv, zv))
da97a08a 213 write = kTRUE;
214 }
215 }
216 /*
217 else if (part->GetUniqueID() == 13){
218 // Evaporation
219 // Check that we end at a primary particle
220 TParticle* mother = part;
221 Int_t imo = part->GetFirstMother();
222 while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
223 mother = mcE->Stack()->Particle(imo);
224 imo = mother->GetFirstMother();
225 }
226 // Select according to pseudorapidity and production point
227 if (imo < nprim && Select(mother, rv, zv))
228 write = kTRUE;
229 }
230 */
231 if (write) {
232 if(mcH)mcH->SelectParticle(ip);
233 j++;
234 }
235 }
236
63892dd0 237
238 // check for charm daughters
239 static int iTaken = 0;
240 static int iAll = 0;
241 static int iCharm = 0;
242 for (Int_t ip = 0; ip < np; ip++){
243 AliMCParticle* mcpart = mcE->GetTrack(ip);
244 TParticle* part = mcpart->Particle();
245
246 // if((TMath::Abs(part->GetPdgCode())/400)==1){
247 if((TMath::Abs(part->GetPdgCode()))==411){
248 // cases
249 iCharm++;
93836e1b 250 Int_t d0 = mcpart->GetFirstDaughter();
251 Int_t d1 = mcpart->GetLastDaughter();
63892dd0 252 if(d0>0&&d1>0){
253 for(int id = d0;id <= d1;id++){
93836e1b 254 AliMCParticle* daughter = mcE->GetTrack(id);
63892dd0 255 iAll++;
256 if(mcH->IsParticleSelected(id))iTaken++;
257 }
258 }
259 }
260 }
a3c57861 261
262 AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
63892dd0 263
264
da97a08a 265 return;
266}
267
268Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
269{
270 // Selection accoring to eta of the mother and production point
271 // This has to be refined !!!!!!
272
273 // Esp if we don't have collisison in the central barrel but e.g. beam gas
274 // return kTRUE;
275
276 Float_t eta = part->Eta();
277
278 // central barrel consider everything in the ITS...
279 // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
280 // larger for V0s in the TPC
281 // if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
282
63892dd0 283 if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
284
da97a08a 285 // Andreas' Cuts
63892dd0 286 // if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
287
288
da97a08a 289
290 // Muon arm
291 if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
292
293 // PMD acceptance 2.3 <= eta < = 3.5
294 // if(eta>2.0&&eta<3.7)return kTRUE;
295
296 return kFALSE;
297
298}
299