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