Reaction plane angle added (A. Kisiel)
[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"
a3b51fd2 24#include "TParameter.h"
da97a08a 25#include "TString.h"
a3b51fd2 26#include "TList.h"
27#include "TProfile.h"
28#include "TH1F.h"
29
da97a08a 30
31#include "AliAnalysisTaskMCParticleFilter.h"
32#include "AliAnalysisManager.h"
33#include "AliAnalysisFilter.h"
34#include "AliHeader.h"
35#include "AliStack.h"
36#include "AliMCEvent.h"
37#include "AliMCEventHandler.h"
38#include "AliAODEvent.h"
39#include "AliAODHeader.h"
dce1b636 40#include "AliAODMCHeader.h"
da97a08a 41#include "AliAODHandler.h"
42#include "AliAODVertex.h"
43#include "AliAODMCParticle.h"
44
45
46#include "AliLog.h"
47
48ClassImp(AliAnalysisTaskMCParticleFilter)
49
50////////////////////////////////////////////////////////////////////////
51
52//____________________________________________________________________
53AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
54AliAnalysisTaskSE(),
a3b51fd2 55 fTrackFilterMother(0x0),
56 fHistList(0x0)
da97a08a 57{
58 // Default constructor
59}
60
a3b51fd2 61Bool_t AliAnalysisTaskMCParticleFilter::Notify()
62{
63 //
64 // Implemented Notify() to read the cross sections
65 // and number of trials from pyxsec.root
66 //
67 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
68 Double_t xsection = 0;
69 UInt_t ntrials = 0;
70 if(tree){
71 TFile *curfile = tree->GetCurrentFile();
72 if (!curfile) {
73 Error("Notify","No current file");
74 return kFALSE;
75 }
76
77 TString fileName(curfile->GetName());
78 if(fileName.Contains("AliESDs.root")){
79 fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
80 }
81 else if(fileName.Contains("AliAOD.root")){
82 fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
83 }
84 else if(fileName.Contains("galice.root")){
85 // for running with galice and kinematics alone...
86 fileName.ReplaceAll("galice.root", "pyxsec.root");
87 }
88
89
90 TFile *fxsec = TFile::Open(fileName.Data());
91 if(!fxsec){
92 AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
93 // not a severe condition we just do not have the information...
94 return kTRUE;
95 }
96 TTree *xtree = (TTree*)fxsec->Get("Xsection");
97 if(!xtree){
98 Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
99 return kTRUE;
100 }
101 xtree->SetBranchAddress("xsection",&xsection);
102 xtree->SetBranchAddress("ntrials",&ntrials);
103 xtree->GetEntry(0);
104 ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
105 ((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
106
107 }
108 return kTRUE;
109}
110
111
112
da97a08a 113//____________________________________________________________________
114AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
115 AliAnalysisTaskSE(name),
a3b51fd2 116 fTrackFilterMother(0x0),
117 fHistList(0x0)
da97a08a 118{
119 // Default constructor
a3b51fd2 120 DefineOutput(1,TList::Class());
da97a08a 121}
122
a3b51fd2 123/*
da97a08a 124//____________________________________________________________________
125AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
126 AliAnalysisTaskSE(obj),
127 fTrackFilterMother(obj.fTrackFilterMother)
128{
129 // Copy constructor
130}
a3b51fd2 131*/
da97a08a 132//____________________________________________________________________
133AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
134{
135 // if( fTrackFilterMother ) delete fTrackFilterMother;
136}
137
a3b51fd2 138/*
da97a08a 139//____________________________________________________________________
140AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
141{
142// Assignment
143 if(this!=&other) {
144 AliAnalysisTaskSE::operator=(other);
145 fTrackFilterMother = other.fTrackFilterMother;
146 }
147 return *this;
148}
a3b51fd2 149*/
da97a08a 150//____________________________________________________________________
151void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
152{
153 // Create the output container
154 if (OutputTree()&&fTrackFilterMother)
155 OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
d2740fba 156
157 // this part is mainly needed to set the MCEventHandler
158 // to the AODHandler, this will not be needed when
159 // AODHandler goes to ANALYSISalice
160 // setting in the steering macro will not work on proof :(
161 // so we have to do it in a task
162
163 // the branch booking can also go into the AODHandler then
164
dce1b636 165
166 // mcparticles
da97a08a 167 TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
168 tca->SetName(AliAODMCParticle::StdBranchName());
169 AddAODBranch("TClonesArray",&tca);
da97a08a 170
dce1b636 171 // MC header...
172 AliAODMCHeader *mcHeader = new AliAODMCHeader();
dce1b636 173 mcHeader->SetName(AliAODMCHeader::StdBranchName());
174 AddAODBranch("AliAODMCHeader",&mcHeader);
175
176
177
178 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
da97a08a 179 AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
180 if(!aodH){
181 Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
182 return;
183 }
184 aodH->SetMCEventHandler(mcH);
dce1b636 185
a3b51fd2 186
187 // these are histograms, for averaging and summing
188 // do it via histograms to be PROOF-proof
189 // which merges the results from different workers
190 // histos are not saved reside only in memory
191
192
193
194 fHistList = new TList();
195 TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
196 h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
197 fHistList->Add(h1Xsec);
198 TH1F* h1Trials = new TH1F("h1Trials","trials from pyxsec.root",1,0,1);
199 h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
200 fHistList->Add(h1Trials);
201
da97a08a 202
203}
204
205//____________________________________________________________________
206void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
207{
208// Execute analysis for current event
209//
210
211// Fill AOD tracks from Kinematic stack
212
213 // get AliAOD Event
214 AliAODEvent* aod = AODEvent();
215 if (!aod) {
216 AliWarning("No Output Handler connected, doing nothing !") ;
217 return;
218 }
219
220 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
221 if(!mcH){
222 AliWarning("No MC handler Found");
223 return;
224 }
225
226
227 // fetch the output
228 // Fill control histos
229
230 // tmp array for holding the mctracks
231 // need to set mother and duagther __before__
232 // filling in the tree
233
234 AliMCEvent *mcE = MCEvent();
235
236 if(!mcE){
237 AliWarning("No MC event Found");
238 return;
239 }
240
da97a08a 241 Int_t np = mcE->GetNumberOfTracks();
93836e1b 242 Int_t nprim = mcE->GetNumberOfPrimaries();
da97a08a 243 // TODO ADD MC VERTEX
244
245 // We take all real primaries
246
247
da97a08a 248 Int_t j=0;
249 for (Int_t ip = 0; ip < np; ip++){
7aad0c47 250 AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
da97a08a 251 TParticle* part = mcpart->Particle();
da97a08a 252 Float_t xv = part->Vx();
253 Float_t yv = part->Vy();
254 Float_t zv = part->Vz();
255 Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
256
257 Bool_t write = kFALSE;
63892dd0 258
da97a08a 259 if (ip < nprim) {
260 // Select the primary event
261 write = kTRUE;
262 } else if (part->GetUniqueID() == 4) {
263 // Particles from decay
264 // Check that the decay chain ends at a primary particle
93836e1b 265 AliMCParticle* mother = mcpart;
266 Int_t imo = mcpart->GetMother();
da97a08a 267 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
7aad0c47 268 mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b 269 imo = mother->GetMother();
da97a08a 270 }
63892dd0 271 // Select according to pseudorapidity and production point of primary ancestor
7aad0c47 272 if (imo < nprim && Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kTRUE;
da97a08a 273 } else if (part->GetUniqueID() == 5) {
274 // Now look for pair production
93836e1b 275 Int_t imo = mcpart->GetMother();
da97a08a 276 if (imo < nprim) {
277 // Select, if the gamma is a primary
278 write = kTRUE;
279 } else {
280 // Check if the gamma comes from the decay chain of a primary particle
7aad0c47 281 AliMCParticle* mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b 282 imo = mother->GetMother();
da97a08a 283 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
7aad0c47 284 mother = (AliMCParticle*) mcE->GetTrack(imo);
93836e1b 285 imo = mother->GetMother();
da97a08a 286 }
287 // Select according to pseudorapidity and production point
93836e1b 288 if (imo < nprim && Select(mother->Particle(), rv, zv))
da97a08a 289 write = kTRUE;
290 }
291 }
292 /*
293 else if (part->GetUniqueID() == 13){
294 // Evaporation
295 // Check that we end at a primary particle
296 TParticle* mother = part;
297 Int_t imo = part->GetFirstMother();
298 while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
299 mother = mcE->Stack()->Particle(imo);
300 imo = mother->GetFirstMother();
301 }
302 // Select according to pseudorapidity and production point
303 if (imo < nprim && Select(mother, rv, zv))
304 write = kTRUE;
305 }
306 */
307 if (write) {
308 if(mcH)mcH->SelectParticle(ip);
309 j++;
310 }
311 }
312
63892dd0 313
314 // check for charm daughters
315 static int iTaken = 0;
316 static int iAll = 0;
317 static int iCharm = 0;
318 for (Int_t ip = 0; ip < np; ip++){
7aad0c47 319 AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
63892dd0 320 TParticle* part = mcpart->Particle();
321
322 // if((TMath::Abs(part->GetPdgCode())/400)==1){
323 if((TMath::Abs(part->GetPdgCode()))==411){
324 // cases
325 iCharm++;
93836e1b 326 Int_t d0 = mcpart->GetFirstDaughter();
327 Int_t d1 = mcpart->GetLastDaughter();
63892dd0 328 if(d0>0&&d1>0){
329 for(int id = d0;id <= d1;id++){
a3b51fd2 330 // AliMCParticle* daughter = mcE->GetTrack(id);
63892dd0 331 iAll++;
332 if(mcH->IsParticleSelected(id))iTaken++;
333 }
334 }
335 }
336 }
a3c57861 337
338 AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
63892dd0 339
a3b51fd2 340 PostData(1,fHistList);
63892dd0 341
da97a08a 342 return;
343}
344
a3b51fd2 345void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
346 //
347 // Terminate the execution save the PYTHIA x_section to the UserInfo()
348 //
349
350
351 TTree *tree = (TTree*)GetOutputData(0);
352 fHistList = (TList*)GetOutputData(1);
353 if(!fHistList){
354 Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
355 return;
356 }
357
358 if(!tree){
359 Printf("%s:%d Output tree not found",(char*)__FILE__,__LINE__);
360 return;
361 }
362
363 TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
364 TH1F* hTrials = ((TH1F*)(fHistList->FindObject("h1Trials")));
365 if(!hXsec)return;
366 if(!hTrials)return;
367
368 Float_t xsec = hXsec->GetBinContent(1);
369 Float_t trials = hTrials->Integral();
370 AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
371
372 // This only works if the tree is not a special output
373 // in this case it is already written
374 /*
375 TParameter<float> *x = new TParameter<float>("xsection",xsec);
376 tree->GetUserInfo()->Add(x);
377 TParameter<float> *t = new TParameter<float>("trials",trials);
378 tree->GetUserInfo()->Add(t);
379 */
380
381
382}
383
da97a08a 384Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
385{
386 // Selection accoring to eta of the mother and production point
387 // This has to be refined !!!!!!
388
389 // Esp if we don't have collisison in the central barrel but e.g. beam gas
390 // return kTRUE;
391
392 Float_t eta = part->Eta();
393
394 // central barrel consider everything in the ITS...
395 // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
396 // larger for V0s in the TPC
397 // if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
398
63892dd0 399 if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
400
da97a08a 401 // Andreas' Cuts
63892dd0 402 // if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
403
404
da97a08a 405
406 // Muon arm
407 if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
408
409 // PMD acceptance 2.3 <= eta < = 3.5
410 // if(eta>2.0&&eta<3.7)return kTRUE;
411
412 return kFALSE;
413
414}
415