Line data Source code
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"
24 : #include "TParameter.h"
25 : #include "TString.h"
26 : #include "TList.h"
27 : #include "TProfile.h"
28 : #include "TH1F.h"
29 :
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 "AliESDInputHandler.h"
39 : #include "AliAODEvent.h"
40 : #include "AliAODHeader.h"
41 : #include "AliAODMCHeader.h"
42 : #include "AliAODHandler.h"
43 : #include "AliAODVertex.h"
44 : #include "AliAODMCParticle.h"
45 : #include "AliCollisionGeometry.h"
46 : #include "AliGenDPMjetEventHeader.h"
47 : #include "AliGenHijingEventHeader.h"
48 : #include "AliGenPythiaEventHeader.h"
49 : #include "AliGenCocktailEventHeader.h"
50 : #include "AliGenEventHeaderTunedPbPb.h"
51 : #include "AliESDtrack.h"
52 : #include "AliAODTrack.h"
53 : #include "AliAODPid.h"
54 : #include "AliESDpid.h"
55 :
56 : #include "AliLog.h"
57 :
58 170 : ClassImp(AliAnalysisTaskMCParticleFilter)
59 :
60 : ////////////////////////////////////////////////////////////////////////
61 :
62 : //____________________________________________________________________
63 : AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
64 0 : AliAnalysisTaskSE(),
65 0 : fTrackFilterMother(0x0),
66 0 : fAODMcHeader(0x0),
67 0 : fAODMcParticles(0x0),
68 0 : fHistList(0x0)
69 0 : {
70 : // Default constructor
71 0 : }
72 :
73 : Bool_t AliAnalysisTaskMCParticleFilter::Notify()
74 : {
75 : //
76 : // Implemented Notify() to read the cross sections
77 : // from pyxsec.root
78 : //
79 2 : AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
80 1 : TTree *tree = mgr->GetTree();
81 1 : Double_t xsection = 0;
82 1 : UInt_t ntrials = 0;
83 1 : if(tree){
84 1 : TFile *curfile = tree->GetCurrentFile();
85 1 : if (!curfile) {
86 0 : Error("Notify","No current file");
87 0 : return kFALSE;
88 : }
89 :
90 1 : TString fileName(curfile->GetName());
91 2 : TString datafile = mgr->GetInputEventHandler()->GetInputFileName();
92 2 : if (fileName.Contains(datafile)) {
93 1 : fileName.ReplaceAll(datafile, "pyxsec.root");
94 : }
95 0 : else if(fileName.Contains("AliESDs.root")){
96 0 : fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
97 : }
98 0 : else if(fileName.Contains("AliAOD.root")){
99 0 : fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
100 : }
101 0 : else if(fileName.Contains("galice.root")){
102 : // for running with galice and kinematics alone...
103 0 : fileName.ReplaceAll("galice.root", "pyxsec.root");
104 : }
105 :
106 :
107 2 : TFile *fxsec = TFile::Open(fileName.Data());
108 1 : if(!fxsec){
109 4 : 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 1 : return kTRUE;
112 : }
113 0 : TTree *xtree = (TTree*)fxsec->Get("Xsection");
114 0 : if(!xtree){
115 0 : AliWarning(Form("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__));
116 0 : return kTRUE;
117 : }
118 0 : xtree->SetBranchAddress("xsection",&xsection);
119 0 : xtree->SetBranchAddress("ntrials",&ntrials);
120 0 : xtree->GetEntry(0);
121 0 : ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
122 1 : }
123 0 : return kTRUE;
124 1 : }
125 :
126 :
127 :
128 : //____________________________________________________________________
129 : AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
130 2 : AliAnalysisTaskSE(name),
131 2 : fTrackFilterMother(0x0),
132 2 : fAODMcHeader(0x0),
133 2 : fAODMcParticles(0x0),
134 2 : fHistList(0x0)
135 10 : {
136 : // Default constructor
137 4 : DefineOutput(1, TList::Class());
138 4 : }
139 :
140 : /*
141 : //____________________________________________________________________
142 : AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
143 : AliAnalysisTaskSE(obj),
144 : fTrackFilterMother(obj.fTrackFilterMother)
145 : {
146 : // Copy constructor
147 : }
148 : */
149 : //____________________________________________________________________
150 : AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
151 0 : {
152 :
153 0 : if(fAODMcHeader){
154 0 : delete fAODMcHeader;
155 : }
156 0 : if(fAODMcParticles){
157 0 : fAODMcParticles->Delete();
158 0 : delete fAODMcParticles;
159 : }
160 0 : }
161 :
162 : /*
163 : //____________________________________________________________________
164 : AliAnalysisTaskMCParticleFilter& 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 : }
173 : */
174 : //____________________________________________________________________
175 : void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
176 : {
177 : // Create the output container
178 :
179 :
180 3 : if (OutputTree()&&fTrackFilterMother)
181 0 : OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
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 :
191 :
192 : // mcparticles
193 2 : fAODMcParticles = new TClonesArray("AliAODMCParticle", 0);
194 1 : fAODMcParticles->SetName(AliAODMCParticle::StdBranchName());
195 1 : AddAODBranch("TClonesArray",&fAODMcParticles);
196 :
197 : // MC header...
198 2 : fAODMcHeader = new AliAODMCHeader();
199 1 : fAODMcHeader->SetName(AliAODMCHeader::StdBranchName());
200 1 : AddAODBranch("AliAODMCHeader",&fAODMcHeader);
201 :
202 1 : AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
203 3 : AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
204 1 : if(!aodH){
205 0 : AliWarning("Could not get AODHandler");
206 0 : return;
207 : }
208 1 : aodH->SetMCEventHandler(mcH);
209 :
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 2 : fHistList = new TList();
219 1 : TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
220 1 : h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
221 1 : fHistList->Add(h1Xsec);
222 1 : TH1F* h1Trials = new TH1F("h1Trials","trials from MC header",1,0,1);
223 1 : h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
224 1 : fHistList->Add(h1Trials);
225 :
226 1 : PostData(1,fHistList);
227 :
228 2 : }
229 :
230 : //____________________________________________________________________
231 : void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
232 : {
233 : // Execute analysis for current event
234 : //
235 :
236 : // Fill AOD tracks from Kinematic stack
237 8 : PostData(1,fHistList);
238 :
239 : // get AliAOD Event
240 4 : AliAODEvent* aod = AODEvent();
241 4 : if (!aod) {
242 0 : AliWarning("No Output Handler connected, doing nothing !") ;
243 0 : return;
244 : }
245 :
246 4 : AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
247 4 : if(!mcH){
248 0 : AliWarning("No MC handler Found");
249 0 : return;
250 : }
251 :
252 12 : AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
253 4 : if(!aodH){
254 0 : AliWarning("Could not get AODHandler");
255 0 : return;
256 : }
257 :
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 4 : AliMCEvent *mcE = MCEvent();
268 :
269 4 : if(!mcE){
270 0 : AliWarning("No MC event Found");
271 0 : return;
272 : }
273 :
274 4 : Int_t np = mcE->GetNumberOfTracks();
275 4 : Int_t nprim = mcE->GetNumberOfPrimaries();
276 : // TODO ADD MC VERTEX
277 :
278 : // Get the proper MC Collision Geometry
279 4 : AliGenEventHeader* mcEH = mcE->GenEventHeader();
280 :
281 12 : AliGenPythiaEventHeader *pyH = dynamic_cast<AliGenPythiaEventHeader*>(mcEH);
282 : AliGenHijingEventHeader *hiH = 0;
283 : AliCollisionGeometry *colG = 0;
284 : AliGenDPMjetEventHeader *dpmH = 0;
285 : AliGenEventHeaderTunedPbPb *tunedH = 0;
286 :
287 : // it can be only one save some casts
288 : // assuming PYTHIA and HIJING are the most likely ones...
289 4 : if(!pyH){
290 12 : hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
291 4 : if(!hiH){
292 12 : dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
293 4 : if(!dpmH){
294 12 : tunedH = dynamic_cast<AliGenEventHeaderTunedPbPb*>(mcEH);
295 4 : }
296 : }
297 : }
298 :
299 4 : if (hiH || dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);
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;
306 4 : if (!colG) {
307 12 : AliGenCocktailEventHeader *ccEH = dynamic_cast<AliGenCocktailEventHeader *>(mcEH);
308 4 : if (ccEH) {
309 4 : TList *genHeaders = ccEH->GetHeaders();
310 40 : for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
311 64 : if(!pyH)pyH = dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
312 64 : if(!hiH)hiH = dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
313 64 : if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
314 64 : if(!dpmH)dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
315 : }
316 4 : }
317 4 : }
318 :
319 : // take the trials from the p+p event
320 4 : if(hiH)ntrials = hiH->Trials();
321 4 : if(dpmH)ntrials = dpmH->Trials();
322 4 : if(pyH)ntrials = pyH->Trials();
323 4 : if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials);
324 :
325 :
326 :
327 :
328 4 : if (colG) {
329 0 : fAODMcHeader->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
330 0 : AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
331 0 : }
332 4 : else if(tunedH) {
333 0 : fAODMcHeader->SetReactionPlaneAngle(tunedH->GetPsi2());
334 0 : fAODMcHeader->SetCrossSection(tunedH->GetCentrality());
335 0 : }
336 :
337 :
338 : Int_t j=0;
339 4306 : for (Int_t ip = 0; ip < np; ip++){
340 2149 : AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
341 2149 : TParticle* part = mcpart->Particle();
342 2149 : Float_t xv = part->Vx();
343 2149 : Float_t yv = part->Vy();
344 2149 : Float_t zv = part->Vz();
345 2149 : Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
346 :
347 : Bool_t write = kFALSE;
348 :
349 2149 : if (ip < nprim) {
350 : // Select the primary event
351 : write = kTRUE;
352 2149 : } else if (part->GetUniqueID() == kPDecay) {
353 : // Particles from decay
354 : // Check that the decay chain ends at a primary particle
355 : AliMCParticle* mother = mcpart;
356 66 : Int_t imo = mcpart->GetMother();
357 335 : while((imo >= nprim) && (mother->Particle()->GetUniqueID() == kPDecay)) {
358 57 : mother = (AliMCParticle*) mcE->GetTrack(imo);
359 57 : imo = mother->GetMother();
360 : }
361 : // Select according to pseudorapidity and production point of primary ancestor
362 100 : 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
364 2037 : } else if (part->GetUniqueID() == kPPair) {
365 : // Now look for pair production
366 419 : Int_t imo = mcpart->GetMother();
367 419 : if (imo < nprim) {
368 : // Select, if the gamma is a primary
369 : write = kTRUE;
370 46 : } else {
371 : // Check if the gamma comes from the decay chain of a primary particle
372 373 : AliMCParticle* mother = (AliMCParticle*) mcE->GetTrack(imo);
373 373 : imo = mother->GetMother();
374 1155 : while((imo >= nprim) && (mother->Particle()->GetUniqueID() == kPDecay)) {
375 24 : mother = (AliMCParticle*) mcE->GetTrack(imo);
376 24 : imo = mother->GetMother();
377 : }
378 : // Select according to pseudorapidity and production point
379 409 : if (imo < nprim && Select(mother->Particle(), rv, zv))
380 8 : write = kTRUE;
381 : }
382 419 : }
383 :
384 2149 : if (write) {
385 400 : if(mcH)mcH->SelectParticle(ip);
386 200 : j++;
387 200 : }
388 : }
389 :
390 : /*
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 :
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);
416 : }
417 : }
418 : }
419 : }// if charm
420 : AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
421 : }
422 : */
423 :
424 :
425 4 : aodH->StoreMCParticles();
426 :
427 : return;
428 4 : }
429 :
430 : void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
431 : //
432 : // Terminate the execution save the PYTHIA x_section to the UserInfo()
433 : //
434 :
435 :
436 2 : fHistList = (TList*)GetOutputData(1);
437 1 : if(!fHistList){
438 1 : Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
439 1 : return;
440 : }
441 :
442 0 : TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
443 0 : TH1F* hTrials = ((TH1F*)(fHistList->FindObject("h1Trials")));
444 0 : if(!hXsec)return;
445 0 : if(!hTrials)return;
446 :
447 0 : Float_t xsec = hXsec->GetBinContent(1);
448 0 : Float_t trials = hTrials->Integral();
449 0 : AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
450 :
451 1 : }
452 :
453 : Bool_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 72 : 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 :
468 44 : if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;
469 :
470 : // Andreas' Cuts
471 : // if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;
472 :
473 :
474 :
475 : // Muon arm
476 56 : 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 28 : return kFALSE;
482 :
483 36 : }
484 :
485 : void AliAnalysisTaskMCParticleFilter::PrintMCParticle(const AliMCParticle *mcp,Int_t np){
486 :
487 0 : const TParticle *p = mcp->Particle();
488 0 : Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughte\
489 : r2 %d ",
490 0 : np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter \
491 0 : (0),p->GetDaughter(1));
492 0 : Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
493 0 : p->Print();
494 0 : Printf("---------------------------------------");
495 0 : }
|
