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