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