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