CMake target for TENDER.
[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 "AliAODMCHeader.h"
36 #include "AliAODHandler.h"
37 #include "AliAODVertex.h"
38 #include "AliAODMCParticle.h"
39
40
41 #include "AliLog.h"
42
43 ClassImp(AliAnalysisTaskMCParticleFilter)
44
45 ////////////////////////////////////////////////////////////////////////
46
47 //____________________________________________________________________
48 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
49 AliAnalysisTaskSE(),
50   fTrackFilterMother(0x0)
51 {
52   // Default constructor
53 }
54
55 //____________________________________________________________________
56 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
57     AliAnalysisTaskSE(name),
58     fTrackFilterMother(0x0)
59 {
60   // Default constructor
61 }
62
63
64 //____________________________________________________________________
65 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
66     AliAnalysisTaskSE(obj),
67     fTrackFilterMother(obj.fTrackFilterMother)
68 {
69   // Copy constructor
70 }
71
72 //____________________________________________________________________
73 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
74 {
75   //  if( fTrackFilterMother ) delete fTrackFilterMother;
76 }
77
78
79 //____________________________________________________________________
80 AliAnalysisTaskMCParticleFilter& 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 //____________________________________________________________________
91 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
92 {
93   // Create the output container
94     if (OutputTree()&&fTrackFilterMother) 
95         OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
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
105
106     // mcparticles
107     TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
108     tca->SetName(AliAODMCParticle::StdBranchName());
109     AddAODBranch("TClonesArray",&tca);
110
111     // MC header...
112     AliAODMCHeader *mcHeader = new AliAODMCHeader();
113     mcHeader->SetName(AliAODMCHeader::StdBranchName());
114     AddAODBranch("AliAODMCHeader",&mcHeader);    
115
116     
117
118     AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
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);
125
126     
127 }
128
129 //____________________________________________________________________
130 void 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
165   Int_t np    = mcE->GetNumberOfTracks();
166   Int_t nprim = mcE->GetNumberOfPrimaries();
167   // TODO ADD MC VERTEX
168   
169   // We take all real primaries
170
171
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();
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;
182     
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
189       AliMCParticle* mother = mcpart;
190       Int_t imo = mcpart->GetMother();
191       while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
192         mother = mcE->GetTrack(imo);
193         imo =  mother->GetMother();
194       }
195       // Select according to pseudorapidity and production point of primary ancestor
196       if (imo < nprim && Select(mcE->GetTrack(imo)->Particle(), rv, zv))write = kTRUE;         
197     } else if (part->GetUniqueID() == 5) {
198       // Now look for pair production
199       Int_t imo = mcpart->GetMother();
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
205         AliMCParticle* mother = mcE->GetTrack(imo);
206         imo = mother->GetMother();
207         while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
208           mother =  mcE->GetTrack(imo);
209           imo =  mother->GetMother();
210         }
211         // Select according to pseudorapidity and production point 
212         if (imo < nprim && Select(mother->Particle(), rv, zv)) 
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
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++;
250       Int_t d0 =  mcpart->GetFirstDaughter();
251       Int_t d1 =  mcpart->GetLastDaughter();
252       if(d0>0&&d1>0){
253         for(int id = d0;id <= d1;id++){
254           AliMCParticle* daughter =  mcE->GetTrack(id);
255           iAll++;
256           if(mcH->IsParticleSelected(id))iTaken++;
257         }
258       }
259     }
260   }
261
262   AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
263
264
265   return;
266 }
267
268 Bool_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
283   if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;   
284
285   // Andreas' Cuts
286   //  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   
287
288
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