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