1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Analysis task for Kinematic filtering
18 // Fill AODtrackMC tracks from Kinematic stack
23 #include "TParticle.h"
24 #include "TParameter.h"
31 #include "AliAnalysisTaskMCParticleFilter.h"
32 #include "AliAnalysisManager.h"
33 #include "AliAnalysisFilter.h"
34 #include "AliHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliMCEventHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODHeader.h"
40 #include "AliAODMCHeader.h"
41 #include "AliAODHandler.h"
42 #include "AliAODVertex.h"
43 #include "AliAODMCParticle.h"
44 #include "AliCollisionGeometry.h"
45 #include "AliGenDPMjetEventHeader.h"
46 #include "AliGenHijingEventHeader.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliGenCocktailEventHeader.h"
52 ClassImp(AliAnalysisTaskMCParticleFilter)
54 ////////////////////////////////////////////////////////////////////////
56 //____________________________________________________________________
57 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
59 fTrackFilterMother(0x0),
62 // Default constructor
65 Bool_t AliAnalysisTaskMCParticleFilter::Notify()
68 // Implemented Notify() to read the cross sections
71 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
72 Double_t xsection = 0;
75 TFile *curfile = tree->GetCurrentFile();
77 Error("Notify","No current file");
81 TString fileName(curfile->GetName());
82 if(fileName.Contains("AliESDs.root")){
83 fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
85 else if(fileName.Contains("AliAOD.root")){
86 fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
88 else if(fileName.Contains("galice.root")){
89 // for running with galice and kinematics alone...
90 fileName.ReplaceAll("galice.root", "pyxsec.root");
94 TFile *fxsec = TFile::Open(fileName.Data());
96 AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
97 // not a severe condition we just do not have the information...
100 TTree *xtree = (TTree*)fxsec->Get("Xsection");
102 Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
105 xtree->SetBranchAddress("xsection",&xsection);
106 xtree->SetBranchAddress("ntrials",&ntrials);
108 ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
115 //____________________________________________________________________
116 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
117 AliAnalysisTaskSE(name),
118 fTrackFilterMother(0x0),
121 // Default constructor
122 DefineOutput(1,TList::Class());
126 //____________________________________________________________________
127 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
128 AliAnalysisTaskSE(obj),
129 fTrackFilterMother(obj.fTrackFilterMother)
134 //____________________________________________________________________
135 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
137 // if( fTrackFilterMother ) delete fTrackFilterMother;
141 //____________________________________________________________________
142 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
146 AliAnalysisTaskSE::operator=(other);
147 fTrackFilterMother = other.fTrackFilterMother;
152 //____________________________________________________________________
153 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
155 // Create the output container
156 if (OutputTree()&&fTrackFilterMother)
157 OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
159 // this part is mainly needed to set the MCEventHandler
160 // to the AODHandler, this will not be needed when
161 // AODHandler goes to ANALYSISalice
162 // setting in the steering macro will not work on proof :(
163 // so we have to do it in a task
165 // the branch booking can also go into the AODHandler then
169 TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
170 tca->SetName(AliAODMCParticle::StdBranchName());
171 AddAODBranch("TClonesArray",&tca);
174 AliAODMCHeader *mcHeader = new AliAODMCHeader();
175 mcHeader->SetName(AliAODMCHeader::StdBranchName());
176 AddAODBranch("AliAODMCHeader",&mcHeader);
180 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
181 AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
183 Printf("%s:&d Could not get AODHandler",(char*)__FILE__,__LINE__);
186 aodH->SetMCEventHandler(mcH);
189 // these are histograms, for averaging and summing
190 // do it via histograms to be PROOF-proof
191 // which merges the results from different workers
192 // histos are not saved reside only in memory
196 fHistList = new TList();
197 TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
198 h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
199 fHistList->Add(h1Xsec);
200 TH1F* h1Trials = new TH1F("h1Trials","trials from MC header",1,0,1);
201 h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
202 fHistList->Add(h1Trials);
207 //____________________________________________________________________
208 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
210 // Execute analysis for current event
213 // Fill AOD tracks from Kinematic stack
216 AliAODEvent* aod = AODEvent();
218 AliWarning("No Output Handler connected, doing nothing !") ;
222 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
224 AliWarning("No MC handler Found");
230 // Fill control histos
232 // tmp array for holding the mctracks
233 // need to set mother and duagther __before__
234 // filling in the tree
236 AliMCEvent *mcE = MCEvent();
239 AliWarning("No MC event Found");
243 Int_t np = mcE->GetNumberOfTracks();
244 Int_t nprim = mcE->GetNumberOfPrimaries();
245 // TODO ADD MC VERTEX
247 AliAODMCHeader *aodMCHo = (AliAODMCHeader *) aod->FindListObject("mcHeader");
248 // Get the proper MC Collision Geometry
249 AliGenEventHeader* mcEH = mcE->GenEventHeader();
251 AliGenPythiaEventHeader *pyH = dynamic_cast<AliGenPythiaEventHeader*>(mcEH);
252 AliGenHijingEventHeader *hiH = 0;
253 AliCollisionGeometry *colG = 0;
254 AliGenDPMjetEventHeader *dpmH = 0;
255 // it can be only one save some casts
256 // assuming PYTHIA and HIJING are the most likely ones...
258 hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
260 colG = dynamic_cast<AliCollisionGeometry *>(mcEH);
261 if(!colG)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
266 // fetch the trials on a event by event basis, not from pyxsec.root otherwise
267 // we will get a problem when running on proof since Notify may be called
268 // more than once per file
269 // consider storing this information in the AOD output via AliAODHandler
272 AliGenCocktailEventHeader *ccEH = dynamic_cast<AliGenCocktailEventHeader *>(mcEH);
274 TList *genHeaders = ccEH->GetHeaders();
275 for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
276 if(!pyH)dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
277 if(!hiH)dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
278 if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
279 if(!dpmH)dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
284 // take the trials from the p+p event
285 if(hiH)ntrials = hiH->Trials();
286 if(dpmH)ntrials = dpmH->Trials();
287 if(pyH)ntrials = pyH->Trials();
288 if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
294 aodMCHo->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
295 printf("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle());
300 // check varibales for charm need all daughters
301 static int iTaken = 0;
303 static int iCharm = 0;
307 for (Int_t ip = 0; ip < np; ip++){
308 AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
309 TParticle* part = mcpart->Particle();
310 Float_t xv = part->Vx();
311 Float_t yv = part->Vy();
312 Float_t zv = part->Vz();
313 Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
315 Bool_t write = kFALSE;
318 // Select the primary event
320 } else if (part->GetUniqueID() == 4) {
321 // Particles from decay
322 // Check that the decay chain ends at a primary particle
323 AliMCParticle* mother = mcpart;
324 Int_t imo = mcpart->GetMother();
325 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
326 mother = (AliMCParticle*) mcE->GetTrack(imo);
327 imo = mother->GetMother();
329 // Select according to pseudorapidity and production point of primary ancestor
330 if (imo < nprim && Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kTRUE;
331 } else if (part->GetUniqueID() == 5) {
332 // Now look for pair production
333 Int_t imo = mcpart->GetMother();
335 // Select, if the gamma is a primary
338 // Check if the gamma comes from the decay chain of a primary particle
339 AliMCParticle* mother = (AliMCParticle*) mcE->GetTrack(imo);
340 imo = mother->GetMother();
341 while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
342 mother = (AliMCParticle*) mcE->GetTrack(imo);
343 imo = mother->GetMother();
345 // Select according to pseudorapidity and production point
346 if (imo < nprim && Select(mother->Particle(), rv, zv))
351 else if (part->GetUniqueID() == 13){
353 // Check that we end at a primary particle
354 TParticle* mother = part;
355 Int_t imo = part->GetFirstMother();
356 while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
357 mother = mcE->Stack()->Particle(imo);
358 imo = mother->GetFirstMother();
360 // Select according to pseudorapidity and production point
361 if (imo < nprim && Select(mother, rv, zv))
366 if(mcH)mcH->SelectParticle(ip);
369 // debug info to check fro charm daugthers
370 if((TMath::Abs(part->GetPdgCode()))==411){
372 Int_t d0 = mcpart->GetFirstDaughter();
373 Int_t d1 = mcpart->GetLastDaughter();
375 for(int id = d0;id <= d1;id++){
377 if(mcH->IsParticleSelected(id))iTaken++;
384 AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
386 PostData(1,fHistList);
391 void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
393 // Terminate the execution save the PYTHIA x_section to the UserInfo()
397 TTree *tree = (TTree*)GetOutputData(0);
398 fHistList = (TList*)GetOutputData(1);
400 Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
405 Printf("%s:%d Output tree not found",(char*)__FILE__,__LINE__);
409 TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
410 TH1F* hTrials = ((TH1F*)(fHistList->FindObject("h1Trials")));
414 Float_t xsec = hXsec->GetBinContent(1);
415 Float_t trials = hTrials->Integral();
416 AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
418 // This only works if the tree is not a special output
419 // in this case it is already written
421 TParameter<float> *x = new TParameter<float>("xsection",xsec);
422 tree->GetUserInfo()->Add(x);
423 TParameter<float> *t = new TParameter<float>("trials",trials);
424 tree->GetUserInfo()->Add(t);
430 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
432 // Selection accoring to eta of the mother and production point
433 // This has to be refined !!!!!!
435 // Esp if we don't have collisison in the central barrel but e.g. beam gas
438 Float_t eta = part->Eta();
440 // central barrel consider everything in the ITS...
441 // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
442 // larger for V0s in the TPC
443 // if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
445 if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
448 // if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
453 if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
455 // PMD acceptance 2.3 <= eta < = 3.5
456 // if(eta>2.0&&eta<3.7)return kTRUE;