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