870452f9ddd1d86e4eba89b3167a9870d4c54a89
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.cxx
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
42 ClassImp(AliAnalysisTaskMCParticleFilter)
43
44 ////////////////////////////////////////////////////////////////////////
45
46 //____________________________________________________________________
47 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
48 AliAnalysisTaskSE(),
49   fTrackFilterMother(0x0)
50 {
51   // Default constructor
52 }
53
54 //____________________________________________________________________
55 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
56     AliAnalysisTaskSE(name),
57     fTrackFilterMother(0x0)
58 {
59   // Default constructor
60 }
61
62
63 //____________________________________________________________________
64 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
65     AliAnalysisTaskSE(obj),
66     fTrackFilterMother(obj.fTrackFilterMother)
67 {
68   // Copy constructor
69 }
70
71 //____________________________________________________________________
72 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
73 {
74   //  if( fTrackFilterMother ) delete fTrackFilterMother;
75 }
76
77
78 //____________________________________________________________________
79 AliAnalysisTaskMCParticleFilter& 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 //____________________________________________________________________
90 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
91 {
92   // Create the output container
93     if (OutputTree()&&fTrackFilterMother) 
94         OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
95     // how is this is reset cleared in the UserExec....
96     // Can this be handled by the framework?
97     TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
98     tca->SetName(AliAODMCParticle::StdBranchName());
99     AddAODBranch("TClonesArray",&tca);
100     AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
101
102     AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
103     if(!aodH){
104       Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
105       return;
106     }
107     aodH->SetMCEventHandler(mcH);
108     // TODO ADD MC VERTEX
109     
110 }
111
112 //____________________________________________________________________
113 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
114 {
115 // Execute analysis for current event
116 //
117
118 // Fill AOD tracks from Kinematic stack
119     
120   // get AliAOD Event 
121   AliAODEvent* aod = AODEvent();
122   if (!aod) {
123       AliWarning("No Output Handler connected, doing nothing !") ;
124       return;
125   }
126
127   AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
128   if(!mcH){ 
129     AliWarning("No MC handler Found");
130     return;
131   }
132
133
134   // fetch the output 
135   // Fill control histos
136
137   // tmp array for holding the mctracks
138   // need to set mother and duagther __before__
139   // filling in the tree
140
141   AliMCEvent *mcE = MCEvent();
142
143   if(!mcE){
144     AliWarning("No MC event Found");
145     return;
146   }
147
148   AliStack* stack = mcE->Stack();
149   Int_t np    = mcE->GetNumberOfTracks();
150   Int_t nprim = stack->GetNprimary();
151   // TODO ADD MC VERTEX
152   
153   // We take all real primaries
154
155
156             
157   Int_t j=0;
158   for (Int_t ip = 0; ip < np; ip++){
159     AliMCParticle* mcpart =  mcE->GetTrack(ip);
160     TParticle* part = mcpart->Particle();
161     
162     Float_t xv = part->Vx();
163     Float_t yv = part->Vy();
164     Float_t zv = part->Vz();
165     Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
166       
167     Bool_t write = kFALSE;
168     Int_t flag = 0;
169
170
171     if (ip < nprim) {
172       // Select the primary event
173       write = kTRUE;
174     } else if (part->GetUniqueID() == 4) {
175       // Particles from decay
176       // Check that the decay chain ends at a primary particle
177       TParticle* mother = part;
178       Int_t imo = part->GetFirstMother();
179       while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
180         mother =  mcE->Stack()->Particle(imo);
181         imo =  mother->GetFirstMother();
182       }
183       // Select according to pseudorapidity and production point 
184         if (imo < nprim && Select(mother, rv, zv)) 
185           write = kTRUE;
186     } else if (part->GetUniqueID() == 5) {
187       // Now look for pair production
188       continue;
189       Int_t imo = part->GetFirstMother();
190       if (imo < nprim) {
191         // Select, if the gamma is a primary
192         write = kTRUE;
193       } else {
194         // Check if the gamma comes from the decay chain of a primary particle
195         TParticle* mother = mcE->Stack()->Particle(imo);
196         imo = mother->GetFirstMother();
197         while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
198           mother =  mcE->Stack()->Particle(imo);
199           imo =  mother->GetFirstMother();
200         }
201         // Select according to pseudorapidity and production point 
202         if (imo < nprim && Select(mother, rv, zv)) 
203           write = kTRUE;
204       }
205     }
206     /*
207     else if (part->GetUniqueID() == 13){
208       // Evaporation
209       // Check that we end at a primary particle
210       TParticle* mother = part;
211       Int_t imo = part->GetFirstMother();
212       while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
213         mother =  mcE->Stack()->Particle(imo);
214         imo =  mother->GetFirstMother();
215       }
216       // Select according to pseudorapidity and production point 
217         if (imo < nprim && Select(mother, rv, zv)) 
218           write = kTRUE;
219     }
220     */    
221     if (write) {
222       if(mcH)mcH->SelectParticle(ip);
223       j++;
224     }
225   }
226
227   return;
228 }
229
230 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
231 {
232   // Selection accoring to eta of the mother and production point
233   // This has to be refined !!!!!!
234
235   // Esp if we don't have collisison in the central barrel but e.g. beam gas
236   //  return kTRUE;
237
238   Float_t eta = part->Eta();
239
240   // central barrel consider everything in the ITS...
241   // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
242   // larger for V0s in the TPC
243   //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
244
245   // Andreas' Cuts
246   if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   
247
248   // Muon arm
249   if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
250
251   // PMD acceptance 2.3 <= eta < = 3.5
252   //  if(eta>2.0&&eta<3.7)return kTRUE; 
253
254   return kFALSE;
255  
256 }
257