1 /*************************************************************************
3 * Task for fast embedding *
4 * read extra input from AOD *
6 *************************************************************************/
9 /**************************************************************************
10 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
12 * Author: The ALICE Off-line Project. *
13 * Contributors are mentioned in the code where appropriate. *
15 * Permission to use, copy, modify and distribute this software and its *
16 * documentation strictly for non-commercial purposes is hereby granted *
17 * without fee, provided that the above copyright notice appears in all *
18 * copies and that both the copyright notice and this permission notice *
19 * appear in the supporting documentation. The authors make no claims *
20 * about the suitability of this software for any purpose. It is *
21 * provided "as is" without express or implied warranty. *
22 **************************************************************************/
29 #include <TClonesArray.h>
30 #include <TDirectory.h>
40 #include "AliAnalysisTaskFastEmbedding.h"
41 #include "AliAnalysisManager.h"
42 #include "AliESDEvent.h"
43 #include "AliESDtrack.h"
44 #include "AliAODEvent.h"
45 #include "AliAODTrack.h"
46 #include "AliAODJet.h"
47 #include "AliAODMCParticle.h"
48 #include "AliAODMCHeader.h"
49 #include "AliInputEventHandler.h"
53 ClassImp(AliAnalysisTaskFastEmbedding)
55 //__________________________________________________________________________
56 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
68 ,fAODPath("AliAOD.root")
72 ,fOfflineTrgMask(AliVEvent::kAny)
82 ,fTrackBranch("aodExtraTracks")
83 ,fMCparticlesBranch("aodExtraMCparticles")
92 ,fEvtSelMinJetEta(-999.)
93 ,fEvtSelMaxJetEta( 999.)
95 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
100 ,fToyDistributionTrackPt(0.)
101 ,fToyMinTrackEta(-.5)
104 ,fToyMaxTrackPhi(2*TMath::Pi())
115 ,fHistEvtSelection(0)
134 // default constructor
139 //__________________________________________________________________________
140 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
141 : AliAnalysisTaskSE(name)
152 ,fAODPath("AliAOD.root")
156 ,fOfflineTrgMask(AliVEvent::kAny)
165 ,fNInputTracksMax(-1)
166 ,fTrackBranch("aodExtraTracks")
167 ,fMCparticlesBranch("aodExtraMCparticles")
176 ,fEvtSelMinJetEta(-999.)
177 ,fEvtSelMaxJetEta( 999.)
178 ,fEvtSelMinJetPhi(0.)
179 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
180 ,fToyMinNbOfTracks(1)
181 ,fToyMaxNbOfTracks(1)
184 ,fToyDistributionTrackPt(0.)
185 ,fToyMinTrackEta(-.5)
188 ,fToyMaxTrackPhi(2*TMath::Pi())
199 ,fHistEvtSelection(0)
219 DefineOutput(1, TList::Class());
223 //__________________________________________________________________________
224 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding ©)
225 : AliAnalysisTaskSE()
227 ,fAODout(copy.fAODout)
228 ,fAODevent(copy.fAODevent)
229 ,fAODtree(copy.fAODtree)
230 ,fAODfile(copy.fAODfile)
231 ,mcHeader(copy.mcHeader)
233 ,fInputEntries(copy.fInputEntries)
234 ,fAODPathArray(copy.fAODPathArray)
235 ,fAODEntriesArray(copy.fAODEntriesArray)
236 ,fAODPath(copy.fAODPath)
237 ,fAODEntries(copy.fAODEntries)
238 ,fAODEntriesSum(copy.fAODEntriesSum)
239 ,fAODEntriesMax(copy.fAODEntriesMax)
240 ,fOfflineTrgMask(copy.fOfflineTrgMask)
241 ,fMinContribVtx(copy.fMinContribVtx)
242 ,fVtxZMin(copy.fVtxZMin)
243 ,fVtxZMax(copy.fVtxZMax)
244 ,fEvtClassMin(copy.fEvtClassMin)
245 ,fEvtClassMax(copy.fEvtClassMax)
246 ,fCentMin(copy.fCentMin)
247 ,fCentMax(copy.fCentMax)
248 ,fNInputTracksMin(copy.fNInputTracksMin)
249 ,fNInputTracksMax(copy.fNInputTracksMax)
250 ,fTrackBranch(copy.fTrackBranch)
251 ,fMCparticlesBranch(copy.fMCparticlesBranch)
252 ,fJetBranch(copy.fJetBranch)
253 ,fFileId(copy.fFileId)
254 ,fAODEntry(copy.fAODEntry)
255 ,fCountEvents(copy.fCountEvents)
256 ,fEmbedMode(copy.fEmbedMode)
257 ,fEvtSelecMode(copy.fEvtSelecMode)
258 ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
259 ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
260 ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
261 ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
262 ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
263 ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
264 ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
265 ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
266 ,fToyMinTrackPt(copy.fToyMinTrackPt)
267 ,fToyMaxTrackPt(copy.fToyMaxTrackPt)
268 ,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
269 ,fToyMinTrackEta(copy.fToyMinTrackEta)
270 ,fToyMaxTrackEta(copy.fToyMaxTrackEta)
271 ,fToyMinTrackPhi(copy.fToyMinTrackPhi)
272 ,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
273 ,fToyFilterMap(copy.fToyFilterMap)
274 ,fTrackFilterMap(copy.fTrackFilterMap)
275 ,fNPtHard(copy.fNPtHard)
276 ,fPtHard(copy.fPtHard)
277 ,fPtHardBin(copy.fPtHardBin)
278 ,fAODJets(copy.fAODJets)
279 ,fNevents(copy.fNevents)
280 ,fXsection(copy.fXsection)
281 ,fAvgTrials(copy.fAvgTrials)
282 ,fHistList(copy.fHistList)
283 ,fHistEvtSelection(copy.fHistEvtSelection)
284 ,fh1Xsec(copy.fh1Xsec)
285 ,fh1Trials(copy.fh1Trials)
286 ,fh1TrialsEvtSel(copy.fh1TrialsEvtSel)
287 ,fh2PtHard(copy.fh2PtHard)
288 ,fh2PtHardEvtSel(copy.fh2PtHardEvtSel)
289 ,fh2PtHardTrials(copy.fh2PtHardTrials)
290 ,fh1TrackPt(copy.fh1TrackPt)
291 ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
292 ,fh1TrackN(copy.fh1TrackN)
293 ,fh1JetPt(copy.fh1JetPt)
294 ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
295 ,fh1JetN(copy.fh1JetN)
296 ,fh1MCTrackPt(copy.fh1MCTrackPt)
297 ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
298 ,fh1MCTrackN(copy.fh1MCTrackN)
299 ,fh1AODfile(copy.fh1AODfile)
300 ,fh2AODevent(copy.fh2AODevent)
306 //__________________________________________________________________________
307 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
312 AliAnalysisTaskSE::operator=(o);
315 fAODevent = o.fAODevent;
316 fAODtree = o.fAODtree;
317 fAODfile = o.fAODfile;
318 mcHeader = o.mcHeader;
320 fInputEntries = o.fInputEntries;
321 fAODPathArray = o.fAODPathArray;
322 fAODEntriesArray = o.fAODEntriesArray;
323 fAODPath = o.fAODPath;
324 fAODEntries = o.fAODEntries;
325 fAODEntriesSum = o.fAODEntriesSum;
326 fAODEntriesMax = o.fAODEntriesMax;
327 fOfflineTrgMask = o.fOfflineTrgMask;
328 fMinContribVtx = o.fMinContribVtx;
329 fVtxZMin = o.fVtxZMin;
330 fVtxZMax = o.fVtxZMax;
331 fEvtClassMin = o.fEvtClassMin;
332 fEvtClassMax = o.fEvtClassMax;
333 fCentMin = o.fCentMin;
334 fCentMax = o.fCentMax;
335 fNInputTracksMin = o.fNInputTracksMin;
336 fNInputTracksMax = o.fNInputTracksMax;
337 fTrackBranch = o.fTrackBranch;
338 fMCparticlesBranch = o.fMCparticlesBranch;
339 fJetBranch = o.fJetBranch;
341 fAODEntry = o.fAODEntry;
342 fCountEvents = o.fCountEvents;
343 fEmbedMode = o.fEmbedMode;
344 fEvtSelecMode = o.fEvtSelecMode;
345 fEvtSelMinJetPt = o.fEvtSelMinJetPt;
346 fEvtSelMaxJetPt = o.fEvtSelMaxJetPt;
347 fEvtSelMinJetEta = o.fEvtSelMinJetEta;
348 fEvtSelMaxJetEta = o.fEvtSelMaxJetEta;
349 fEvtSelMinJetPhi = o.fEvtSelMinJetPhi;
350 fEvtSelMaxJetPhi = o.fEvtSelMaxJetPhi;
351 fToyMinNbOfTracks = o.fToyMinNbOfTracks;
352 fToyMaxNbOfTracks = o.fToyMaxNbOfTracks;
353 fToyMinTrackPt = o.fToyMinTrackPt;
354 fToyMaxTrackPt = o.fToyMaxTrackPt;
355 fToyDistributionTrackPt = o.fToyDistributionTrackPt;
356 fToyMinTrackEta = o.fToyMinTrackEta;
357 fToyMaxTrackEta = o.fToyMaxTrackEta;
358 fToyMinTrackPhi = o.fToyMinTrackPhi;
359 fToyMaxTrackPhi = o.fToyMaxTrackPhi;
360 fToyFilterMap = o.fToyFilterMap;
361 fTrackFilterMap = o.fTrackFilterMap;
362 fNPtHard = o.fNPtHard;
364 fPtHardBin = o.fPtHardBin;
365 fAODJets = o.fAODJets;
366 fNevents = o.fNevents;
367 fXsection = o.fXsection;
368 fAvgTrials = o.fAvgTrials;
369 fHistList = o.fHistList;
370 fHistEvtSelection = o.fHistEvtSelection;
372 fh1Trials = o.fh1Trials;
373 fh1TrialsEvtSel = o.fh1TrialsEvtSel;
374 fh2PtHard = o.fh2PtHard;
375 fh2PtHardEvtSel = o.fh2PtHardEvtSel;
376 fh2PtHardTrials = o.fh2PtHardTrials;
377 fh1TrackPt = o.fh1TrackPt;
378 fh2TrackEtaPhi = o.fh2TrackEtaPhi;
379 fh1TrackN = o.fh1TrackN;
380 fh1JetPt = o.fh1JetPt;
381 fh2JetEtaPhi = o.fh2JetEtaPhi;
383 fh1MCTrackPt = o.fh1MCTrackPt;
384 fh2MCTrackEtaPhi = o.fh2MCTrackEtaPhi;
385 fh1MCTrackN = o.fh1MCTrackN;
386 fh1AODfile = o.fh1AODfile;
387 fh2AODevent = o.fh2AODevent;
394 //__________________________________________________________________________
395 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
402 //__________________________________________________________________________
403 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
405 // create output objects
406 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
407 AliLog::SetClassDebugLevel("AliAnalysisTaskFastEmbedding", AliLog::kInfo);
410 if(!fHistList) fHistList = new TList();
411 fHistList->SetOwner(kTRUE);
415 rndm = new TRandom3();
416 Int_t id = GetJobID();
417 if(id>-1) rndm->SetSeed(id);
418 else rndm->SetSeed(); // a TTUID is generated and used for seed
419 AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
423 // embed mode with AOD
424 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
427 fFileId = OpenAODfile();
428 fNevents = 0; // force to open another aod in UserExec()
431 PostData(1, fHistList);
434 } //end: embed mode with AOD
437 // connect output aod
438 // create a new branch for extra tracks
439 fAODout = AODEvent();
441 AliError("Output AOD not found.");
442 PostData(1, fHistList);
445 if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
446 AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
447 TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
448 tracks->SetName(fTrackBranch.Data());
449 AddAODBranch("TClonesArray", &tracks);
451 // create new branch for extra mcparticle if available as input
452 if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
453 AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
454 TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
455 mcparticles->SetName(fMCparticlesBranch.Data());
456 AddAODBranch("TClonesArray", &mcparticles);
462 Bool_t oldStatus = TH1::AddDirectoryStatus();
463 TH1::AddDirectory(kFALSE);
465 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
466 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
467 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
468 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
469 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
470 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
471 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
473 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root;p_{T,hard} bin;<#sigma>",fNPtHard+1,-1.5,fNPtHard-0.5);
474 fh1Trials = new TH1F("fh1Trials","trials (simulation) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
475 fh1TrialsEvtSel = new TH1F("fh1TrialsEvtSel","trials (event selection) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
476 fh2PtHard = new TH2F("fh2PtHard","PYTHIA Pt hard;p_{T,hard} bin;p_{T,hard}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
477 fh2PtHardEvtSel = new TH2F("fh2PtHardEvtSel","PYTHIA Pt hard (event selection);p_{T,hard} bin;p_{T,hard}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
478 fh2PtHardTrials = new TH2F("fh2PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard} bin;#sum{p_{T,hard}}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
480 fHistList->Add(fHistEvtSelection);
481 fHistList->Add(fh1Xsec);
482 fHistList->Add(fh1Trials);
483 fHistList->Add(fh1TrialsEvtSel);
484 fHistList->Add(fh2PtHard);
485 fHistList->Add(fh2PtHardEvtSel);
486 fHistList->Add(fh2PtHardTrials);
488 fh1TrackPt = new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
489 fh2TrackEtaPhi = new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
490 fh1TrackN = new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
492 fHistList->Add(fh1TrackPt);
493 fHistList->Add(fh2TrackEtaPhi);
494 fHistList->Add(fh1TrackN);
496 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
498 fh1JetPt = new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 120, 0., 120.);
499 fh2JetEtaPhi = new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
500 fh1JetN = new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
502 fHistList->Add(fh1JetPt);
503 fHistList->Add(fh2JetEtaPhi);
504 fHistList->Add(fh1JetN);
508 if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
510 fh1MCTrackPt = new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
511 fh2MCTrackEtaPhi = new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
512 fh1MCTrackN = new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
514 fHistList->Add(fh1MCTrackPt);
515 fHistList->Add(fh2MCTrackEtaPhi);
516 fHistList->Add(fh1MCTrackN);
520 fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 2300, -0.5, 2299.5);
521 fh2AODevent = new TH2I("fh1AODevent","selected events;file;event", 2500,-.5,2499.5,5000,-.5,4999.5);
522 fHistList->Add(fh1AODfile);
523 fHistList->Add(fh2AODevent);
525 // =========== Switch on Sumw2 for all histos ===========
526 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
527 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
534 TH1::AddDirectory(oldStatus);
536 PostData(1, fHistList);
540 //__________________________________________________________________________
541 void AliAnalysisTaskFastEmbedding::Init()
544 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
548 //__________________________________________________________________________
549 Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
552 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserNotify()");
554 // get total nb of events in tree (of this subjob)
555 AliInputEventHandler* inputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
556 fInputEntries = inputHandler->GetTree()->GetEntriesFast();
557 AliInfo(Form("Total nb. of events: %d", fInputEntries));
564 //__________________________________________________________________________
565 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
567 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
570 AliError("Need output AOD, but is not connected.");
571 PostData(1, fHistList);
575 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
577 AliError("ESD not available");
578 PostData(1, fHistList);
582 // -- event selection --
583 fHistEvtSelection->Fill(1); // number of events before event selection
586 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
587 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
588 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
589 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
590 fHistEvtSelection->Fill(2);
591 PostData(1, fHistList);
596 AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
597 Int_t nTracksPrim = primVtx->GetNContributors();
598 if ((nTracksPrim < fMinContribVtx) ||
599 (primVtx->GetZ() < fVtxZMin) ||
600 (primVtx->GetZ() > fVtxZMax) ){
601 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
602 fHistEvtSelection->Fill(3);
603 PostData(1, fHistList);
607 /** takes wrong value, since jet service tasks runs later
608 // event class selection (from jet helper task)
609 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
610 if(fDebug) Printf("Event class %d", eventClass);
611 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
612 fHistEvtSelection->Fill(4);
613 PostData(1, fHistList);
617 // centrality selection
618 AliCentrality *cent = 0x0;
619 Float_t centValue = 0.;
620 if(fESD) cent = fESD->GetCentrality();
621 if(cent) centValue = cent->GetCentralityPercentile("V0M");
622 if(fDebug) printf("centrality: %f\n", centValue);
623 if (centValue < fCentMin || centValue > fCentMax){
624 fHistEvtSelection->Fill(4);
625 PostData(1, fHistList);
630 /* ** not implemented **
631 // multiplicity due to input tracks
632 Int_t nInputTracks = GetNInputTracks();
634 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
635 fHistEvtSelection->Fill(5);
636 PostData(1, fHistList);
641 fHistEvtSelection->Fill(0); // accepted events
642 // -- end event selection --
645 TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
647 AliError("Extra track branch not found in output.");
648 PostData(1, fHistList);
656 // === embed mode with AOD ===
657 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
659 AliError("Need input AOD, but is not connected.");
660 PostData(1, fHistList);
666 Bool_t useEntry = kFALSE;
667 while(!useEntry){ // protection need, if no event fulfills requierment
669 fAODEntry++; // go to next event
671 if(fAODEntry>=fNevents) fAODEntry=0;
674 if(fCountEvents>=fNevents){
676 AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
679 AliDebug(AliLog::kDebug, "Select new AOD file ...");
681 fFileId = OpenAODfile();
683 PostData(1, fHistList);
686 fh1AODfile->Fill(fFileId);
687 if(fAODEntries!=fNevents){
688 AliError("File with list of AODs and file with nb. of entries inconsistent");
689 PostData(1, fHistList);
698 fAODtree->GetEvent(fAODEntry);
702 fPtHard = mcHeader->GetPtHard();
703 GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
704 AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
707 AliDebug(AliLog::kDebug,"No mcHeader");
709 fPtHardBin = GetPtHardBin(fPtHard);
711 //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials);
713 // fill event variables for each event
714 fh1Xsec->Fill(fPtHardBin,fXsection);
715 fh2PtHard->Fill(fPtHardBin,fPtHard);
717 fh1Trials->Fill(fPtHardBin, fAvgTrials);
718 fh1TrialsEvtSel->Fill(fPtHardBin);
719 fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
722 if(fEvtSelecMode==kEventsJetPt){
723 Int_t nJets = fAODJets->GetEntries();
724 //Printf("nb. jets: %d",nJets);
725 for(Int_t ij=0; ij<nJets; ++ij){
726 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
729 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
730 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
731 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
732 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
740 if(fEvtSelecMode==kEventsAll){
745 AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
747 fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
748 fh2AODevent->Fill(fFileId,fAODEntry);
750 TClonesArray *mcpartIN = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
751 TClonesArray *mcpartOUT = 0x0;
753 mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
755 AliError("Extra MC particles branch not found in output.");
756 PostData(1, fHistList);
761 AliInfo("No extra MC particles found.");
765 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
766 // loop over input tracks
769 Int_t nJets = fAODJets->GetEntries();
770 Int_t nSelectedJets = 0;
771 AliAODJet *leadJet = 0x0; // used for jet tracks only
773 if(fEmbedMode==kAODFull){
774 nTracks = fAODevent->GetNumberOfTracks();
776 for(Int_t ij=0; ij<nJets; ++ij){
777 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
779 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
780 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
781 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
782 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
784 fh1JetPt->Fill(jet->Pt());
785 fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
790 fh1JetN->Fill(nSelectedJets);
793 if(fEmbedMode==kAODJetTracks){
794 // look for leading jet within selection
795 // get reference tracks for this jet
796 Float_t leadJetPt = 0.;
797 for(Int_t ij=0; ij<nJets; ++ij){
798 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
800 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
801 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
802 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
803 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
804 if(jet->Pt()>leadJetPt){
805 leadJetPt = jet->Pt();
811 nTracks = leadJet->GetRefTracks()->GetEntriesFast();
812 fh1JetPt->Fill(leadJet->Pt());
813 fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
818 fh1TrackN->Fill((Float_t)nTracks);
820 for(Int_t it=0; it<nTracks; ++it){
821 AliAODTrack *tmpTr = 0x0;
822 if(fEmbedMode==kAODFull) tmpTr = fAODevent->GetTrack(it);
823 if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
826 tmpTr->SetStatus(AliESDtrack::kEmbedded);
828 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
829 dummy = (*tracks)[nAODtracks-1];
831 if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
832 fh1TrackPt->Fill(tmpTr->Pt());
833 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
838 Int_t nMCpart = mcpartIN->GetEntriesFast();
840 Int_t nPhysicalPrimary=0;
842 for(Int_t ip=0; ip<nMCpart; ++ip){
843 AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
845 if(fEmbedMode==kAODJetTracks){
846 // jet track? else continue
847 // not implemented yet
851 new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
852 dummy = (*mcpartOUT)[nAODmcpart-1];
854 if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
855 fh1MCTrackPt->Fill(tmpPart->Pt());
856 fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
860 fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
863 } // end: embed all tracks or jet tracks
866 if(fEmbedMode==kAODJet4Mom){
869 Int_t nJets = fAODJets->GetEntries();
870 fh1TrackN->Fill((Float_t)nJets);
871 for(Int_t ij=0; ij<nJets; ++ij){
872 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
874 AliAODTrack *tmpTr = (AliAODTrack*)jet;
875 tmpTr->SetFlags(AliESDtrack::kEmbedded);
877 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
878 dummy = (*tracks)[nAODtracks-1];
880 fh1TrackPt->Fill(tmpTr->Pt());
881 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
884 } // end: embed jets as 4-momenta
887 } //end: embed mode with AOD
890 // === embed mode with toy ===
891 if(fEmbedMode==kToyTracks){
892 Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
894 fh1TrackN->Fill((Float_t)nT);
896 for(Int_t i=0; i<nT; ++i){
899 if(fToyMinTrackPt!=fToyMaxTrackPt){
900 if(fToyDistributionTrackPt==0){
901 pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
903 while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
904 pt = rndm->Exp(fToyDistributionTrackPt); // no power law yet!!
905 pt += fToyMinTrackPt;
911 Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
912 Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
913 phi = TVector2::Phi_0_2pi(phi);
915 if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
917 Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
928 AliAODTrack *tmpTr = new AliAODTrack( -999, // id
939 kFALSE, // used for vertex fit
940 kFALSE, // used for prim vtx fit
941 AliAODTrack::kUndef, // type
942 fToyFilterMap, // select info
943 -999. // chi2 per NDF
945 tmpTr->SetFlags(AliESDtrack::kEmbedded);
947 new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
948 dummy = (*tracks)[nAODtracks-1];
950 fh1TrackPt->Fill(pt);
951 fh2TrackEtaPhi->Fill(eta,phi);
955 } //end: embed mode with toy
958 PostData(1, fHistList);
962 //__________________________________________________________________________
963 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
965 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
967 if(fAODfile && fAODfile->IsOpen())
972 //__________________________________________________________________________
973 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
977 const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
978 //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
980 if(env && strlen(env)){
982 AliInfo(Form("Job index %d", id));
985 AliInfo("Job index not found. Okay if running locally.");
991 //__________________________________________________________________________
993 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
995 Int_t nFiles = fAODPathArray->GetEntries();
997 // choose randomly file and event
999 Float_t tmpProp = -1.;
1000 Int_t pendingEvents = fInputEntries-Entry();
1001 //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1002 while(rndm->Rndm()>tmpProp){
1003 n = rndm->Integer(nFiles);
1004 fAODEntries = fAODEntriesArray->At(n);
1005 Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1006 tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1008 fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1010 AliInfo(Form("Select AOD file %d", n));
1011 TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1013 AliError("Could not get path of aod file.");
1016 fAODPath = objStr->GetString();
1021 //__________________________________________________________________________
1023 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
1025 if(fAODPathArray) fFileId = SelectAODfile();
1026 if(fFileId<0) return -1;
1028 TDirectory *owd = gDirectory;
1029 if (fAODfile) fAODfile->Close();
1030 fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1037 AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1038 fFileId = OpenAODfile(trial+1);
1040 AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1044 AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1050 fAODtree = (TTree*)fAODfile->Get("aodTree");
1053 AliError("AOD tree not found.");
1058 fAODtree->SetBranchStatus("*",0);
1059 fAODtree->SetBranchStatus("header",1);
1060 fAODtree->SetBranchStatus("tracks*",1);
1061 fAODtree->SetBranchStatus("jets*",1);
1062 fAODtree->SetBranchStatus("clusters*",1);
1063 fAODtree->SetBranchStatus("mcparticles*",1);
1064 fAODtree->SetBranchStatus("mcHeader*",1);
1068 fAODevent = new AliAODEvent();
1069 fAODevent->ReadFromTree(fAODtree);
1072 // fetch header, jets, etc. from new file
1073 fNevents = fAODtree->GetEntries();
1074 mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1076 if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1077 else fAODJets = fAODevent->GetJets();
1079 AliError("Could not find jets in AOD. Check jet branch when indicated.");
1083 TFile *curfile = fAODtree->GetCurrentFile();
1085 AliError("No current file.");
1089 Float_t trials = 1.;
1091 PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1094 // construct a poor man average trials
1095 Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1096 if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1099 AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1102 fCountEvents=0; // new file, reset counter
1104 return fFileId; // file position in AOD path array, if array available
1108 //____________________________________________________________________________
1109 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1111 static Float_t ptHard = -1.;
1112 if(bSet) ptHard = newValue;
1117 //____________________________________________________________________________
1118 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1120 const Int_t nBins = 10;
1121 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1124 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1131 //____________________________________________________________________________
1132 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1134 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1135 // This should provide the path to the AOD/ESD file
1136 // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks)
1138 TString file(currFile);
1142 if(file.Contains("root_archive.zip#")){
1143 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1144 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1145 file.Replace(pos+1,20,"");
1148 // not an archive take the basename....
1149 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1152 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1154 // next trial fetch the histgram file
1155 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1157 // not a severe condition but inciate that we have no information
1158 AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1162 // find the tlist we want to be independtent of the name so use the Tkey
1163 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1168 TList *list = dynamic_cast<TList*>(key->ReadObj());
1173 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1174 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1177 } // no tree pyxsec.root
1179 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1185 Double_t xsection = 0;
1186 xtree->SetBranchAddress("xsection",&xsection);
1187 xtree->SetBranchAddress("ntrials",&ntrials);