1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Analysis task for Kinematic filtering
18 // Fill AODtrackMC tracks from Kinematic stack
23 #include "TParticle.h"
26 #include "AliAnalysisTaskMCParticleFilter.h"
27 #include "AliAnalysisManager.h"
28 #include "AliAnalysisFilter.h"
29 #include "AliHeader.h"
31 #include "AliMCEvent.h"
32 #include "AliMCEventHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODHandler.h"
37 #include "AliAODVertex.h"
38 #include "AliAODMCParticle.h"
43 ClassImp(AliAnalysisTaskMCParticleFilter)
45 ////////////////////////////////////////////////////////////////////////
47 //____________________________________________________________________
48 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
50 fTrackFilterMother(0x0)
52 // Default constructor
55 //____________________________________________________________________
56 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
57 AliAnalysisTaskSE(name),
58 fTrackFilterMother(0x0)
60 // Default constructor
64 //____________________________________________________________________
65 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
66 AliAnalysisTaskSE(obj),
67 fTrackFilterMother(obj.fTrackFilterMother)
72 //____________________________________________________________________
73 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
75 // if( fTrackFilterMother ) delete fTrackFilterMother;
79 //____________________________________________________________________
80 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
84 AliAnalysisTaskSE::operator=(other);
85 fTrackFilterMother = other.fTrackFilterMother;
90 //____________________________________________________________________
91 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
93 // Create the output container
94 if (OutputTree()&&fTrackFilterMother)
95 OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
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
103 // the branch booking can also go into the AODHandler then
107 TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
108 tca->SetName(AliAODMCParticle::StdBranchName());
109 AddAODBranch("TClonesArray",&tca);
112 AliAODMCHeader *mcHeader = new AliAODMCHeader();
113 Printf("AODMCHeader %p",mcHeader);
114 Printf("AODMCHeader ** %p",&mcHeader);
115 mcHeader->SetName(AliAODMCHeader::StdBranchName());
116 AddAODBranch("AliAODMCHeader",&mcHeader);
120 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
121 AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
123 Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
126 aodH->SetMCEventHandler(mcH);
131 //____________________________________________________________________
132 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
134 // Execute analysis for current event
137 // Fill AOD tracks from Kinematic stack
140 AliAODEvent* aod = AODEvent();
142 AliWarning("No Output Handler connected, doing nothing !") ;
146 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
148 AliWarning("No MC handler Found");
154 // Fill control histos
156 // tmp array for holding the mctracks
157 // need to set mother and duagther __before__
158 // filling in the tree
160 AliMCEvent *mcE = MCEvent();
163 AliWarning("No MC event Found");
167 AliStack* stack = mcE->Stack();
168 Int_t np = mcE->GetNumberOfTracks();
169 Int_t nprim = stack->GetNprimary();
170 // TODO ADD MC VERTEX
172 // We take all real primaries
176 for (Int_t ip = 0; ip < np; ip++){
177 AliMCParticle* mcpart = mcE->GetTrack(ip);
178 TParticle* part = mcpart->Particle();
179 Float_t xv = part->Vx();
180 Float_t yv = part->Vy();
181 Float_t zv = part->Vz();
182 Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
184 Bool_t write = kFALSE;
187 // Select the primary event
189 } else if (part->GetUniqueID() == 4) {
190 // Particles from decay
191 // Check that the decay chain ends at a primary particle
192 TParticle* mother = part;
193 Int_t imo = part->GetFirstMother();
194 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
195 mother = mcE->Stack()->Particle(imo);
196 imo = mother->GetFirstMother();
198 // Select according to pseudorapidity and production point of primary ancestor
199 if (imo < nprim && Select(mcE->Stack()->Particle(imo), rv, zv))write = kTRUE;
200 } else if (part->GetUniqueID() == 5) {
201 // Now look for pair production
202 Int_t imo = part->GetFirstMother();
204 // Select, if the gamma is a primary
207 // Check if the gamma comes from the decay chain of a primary particle
208 TParticle* mother = mcE->Stack()->Particle(imo);
209 imo = mother->GetFirstMother();
210 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
211 mother = mcE->Stack()->Particle(imo);
212 imo = mother->GetFirstMother();
214 // Select according to pseudorapidity and production point
215 if (imo < nprim && Select(mother, rv, zv))
220 else if (part->GetUniqueID() == 13){
222 // Check that we end at a primary particle
223 TParticle* mother = part;
224 Int_t imo = part->GetFirstMother();
225 while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
226 mother = mcE->Stack()->Particle(imo);
227 imo = mother->GetFirstMother();
229 // Select according to pseudorapidity and production point
230 if (imo < nprim && Select(mother, rv, zv))
235 if(mcH)mcH->SelectParticle(ip);
241 // check for charm daughters
242 static int iTaken = 0;
244 static int iCharm = 0;
245 for (Int_t ip = 0; ip < np; ip++){
246 AliMCParticle* mcpart = mcE->GetTrack(ip);
247 TParticle* part = mcpart->Particle();
249 // if((TMath::Abs(part->GetPdgCode())/400)==1){
250 if((TMath::Abs(part->GetPdgCode()))==411){
253 Printf("Decay Mother %s",part->GetPDG()->GetName());
254 Int_t d0 = part->GetFirstDaughter();
255 Int_t d1 = part->GetLastDaughter();
257 for(int id = d0;id <= d1;id++){
258 TParticle* daughter = mcE->Stack()->Particle(id);
259 Printf("Decay Daughter %s",daughter->GetPDG()->GetName());
261 if(mcH->IsParticleSelected(id))iTaken++;
266 Printf("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm);
272 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
274 // Selection accoring to eta of the mother and production point
275 // This has to be refined !!!!!!
277 // Esp if we don't have collisison in the central barrel but e.g. beam gas
280 Float_t eta = part->Eta();
282 // central barrel consider everything in the ITS...
283 // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
284 // larger for V0s in the TPC
285 // if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
287 if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
290 // if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
295 if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
297 // PMD acceptance 2.3 <= eta < = 3.5
298 // if(eta>2.0&&eta<3.7)return kTRUE;