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()
551 // User defined Notify(), called once
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 // Analysis, called once per event
568 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
571 AliError("Need output AOD, but is not connected.");
572 PostData(1, fHistList);
576 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
578 AliError("ESD not available");
579 PostData(1, fHistList);
583 // -- event selection --
584 fHistEvtSelection->Fill(1); // number of events before event selection
587 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
588 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
589 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
590 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
591 fHistEvtSelection->Fill(2);
592 PostData(1, fHistList);
597 AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
598 Int_t nTracksPrim = primVtx->GetNContributors();
599 if ((nTracksPrim < fMinContribVtx) ||
600 (primVtx->GetZ() < fVtxZMin) ||
601 (primVtx->GetZ() > fVtxZMax) ){
602 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
603 fHistEvtSelection->Fill(3);
604 PostData(1, fHistList);
608 /** takes wrong value, since jet service tasks runs later
609 // event class selection (from jet helper task)
610 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
611 if(fDebug) Printf("Event class %d", eventClass);
612 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
613 fHistEvtSelection->Fill(4);
614 PostData(1, fHistList);
618 // centrality selection
619 AliCentrality *cent = 0x0;
620 Float_t centValue = 0.;
621 if(fESD) cent = fESD->GetCentrality();
622 if(cent) centValue = cent->GetCentralityPercentile("V0M");
623 if(fDebug) printf("centrality: %f\n", centValue);
624 if (centValue < fCentMin || centValue > fCentMax){
625 fHistEvtSelection->Fill(4);
626 PostData(1, fHistList);
631 /* ** not implemented **
632 // multiplicity due to input tracks
633 Int_t nInputTracks = GetNInputTracks();
635 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
636 fHistEvtSelection->Fill(5);
637 PostData(1, fHistList);
642 fHistEvtSelection->Fill(0); // accepted events
643 // -- end event selection --
646 TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
648 AliError("Extra track branch not found in output.");
649 PostData(1, fHistList);
657 // === embed mode with AOD ===
658 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
660 AliError("Need input AOD, but is not connected.");
661 PostData(1, fHistList);
667 Bool_t useEntry = kFALSE;
668 while(!useEntry){ // protection need, if no event fulfills requierment
670 fAODEntry++; // go to next event
672 if(fAODEntry>=fNevents) fAODEntry=0;
675 if(fCountEvents>=fNevents){
677 AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
680 AliDebug(AliLog::kDebug, "Select new AOD file ...");
682 fFileId = OpenAODfile();
684 PostData(1, fHistList);
687 fh1AODfile->Fill(fFileId);
688 if(fAODEntries!=fNevents){
689 AliError("File with list of AODs and file with nb. of entries inconsistent");
690 PostData(1, fHistList);
699 fAODtree->GetEvent(fAODEntry);
703 fPtHard = mcHeader->GetPtHard();
704 GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
705 AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
708 AliDebug(AliLog::kDebug,"No mcHeader");
710 fPtHardBin = GetPtHardBin(fPtHard);
712 //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials);
714 // fill event variables for each event
715 fh1Xsec->Fill(fPtHardBin,fXsection);
716 fh2PtHard->Fill(fPtHardBin,fPtHard);
718 fh1Trials->Fill(fPtHardBin, fAvgTrials);
719 fh1TrialsEvtSel->Fill(fPtHardBin);
720 fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
723 if(fEvtSelecMode==kEventsJetPt){
724 Int_t nJets = fAODJets->GetEntries();
725 //Printf("nb. jets: %d",nJets);
726 for(Int_t ij=0; ij<nJets; ++ij){
727 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
730 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
731 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
732 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
733 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
741 if(fEvtSelecMode==kEventsAll){
746 AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
748 fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
749 fh2AODevent->Fill(fFileId,fAODEntry);
751 TClonesArray *mcpartIN = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
752 TClonesArray *mcpartOUT = 0x0;
754 mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
756 AliError("Extra MC particles branch not found in output.");
757 PostData(1, fHistList);
762 AliInfo("No extra MC particles found.");
766 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
767 // loop over input tracks
770 Int_t nJets = fAODJets->GetEntries();
771 Int_t nSelectedJets = 0;
772 AliAODJet *leadJet = 0x0; // used for jet tracks only
774 if(fEmbedMode==kAODFull){
775 nTracks = fAODevent->GetNumberOfTracks();
777 for(Int_t ij=0; ij<nJets; ++ij){
778 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
780 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
781 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
782 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
783 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
785 fh1JetPt->Fill(jet->Pt());
786 fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
791 fh1JetN->Fill(nSelectedJets);
794 if(fEmbedMode==kAODJetTracks){
795 // look for leading jet within selection
796 // get reference tracks for this jet
797 Float_t leadJetPt = 0.;
798 for(Int_t ij=0; ij<nJets; ++ij){
799 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
801 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
802 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
803 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
804 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
805 if(jet->Pt()>leadJetPt){
806 leadJetPt = jet->Pt();
812 nTracks = leadJet->GetRefTracks()->GetEntriesFast();
813 fh1JetPt->Fill(leadJet->Pt());
814 fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
819 fh1TrackN->Fill((Float_t)nTracks);
821 for(Int_t it=0; it<nTracks; ++it){
822 AliAODTrack *tmpTr = 0x0;
823 if(fEmbedMode==kAODFull) tmpTr = fAODevent->GetTrack(it);
824 if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
827 tmpTr->SetStatus(AliESDtrack::kEmbedded);
829 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
830 dummy = (*tracks)[nAODtracks-1];
832 if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
833 fh1TrackPt->Fill(tmpTr->Pt());
834 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
839 Int_t nMCpart = mcpartIN->GetEntriesFast();
841 Int_t nPhysicalPrimary=0;
843 for(Int_t ip=0; ip<nMCpart; ++ip){
844 AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
846 if(fEmbedMode==kAODJetTracks){
847 // jet track? else continue
848 // not implemented yet
852 new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
853 dummy = (*mcpartOUT)[nAODmcpart-1];
855 if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
856 fh1MCTrackPt->Fill(tmpPart->Pt());
857 fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
861 fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
864 } // end: embed all tracks or jet tracks
867 if(fEmbedMode==kAODJet4Mom){
870 Int_t nJets = fAODJets->GetEntries();
871 fh1TrackN->Fill((Float_t)nJets);
872 for(Int_t ij=0; ij<nJets; ++ij){
873 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
875 AliAODTrack *tmpTr = (AliAODTrack*)jet;
876 tmpTr->SetFlags(AliESDtrack::kEmbedded);
878 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
879 dummy = (*tracks)[nAODtracks-1];
881 fh1TrackPt->Fill(tmpTr->Pt());
882 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
885 } // end: embed jets as 4-momenta
888 } //end: embed mode with AOD
891 // === embed mode with toy ===
892 if(fEmbedMode==kToyTracks){
893 Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
895 fh1TrackN->Fill((Float_t)nT);
897 for(Int_t i=0; i<nT; ++i){
900 if(fToyMinTrackPt!=fToyMaxTrackPt){
901 if(fToyDistributionTrackPt==0){
902 pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
904 while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
905 pt = rndm->Exp(fToyDistributionTrackPt); // no power law yet!!
906 pt += fToyMinTrackPt;
912 Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
913 Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
914 phi = TVector2::Phi_0_2pi(phi);
916 if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
918 Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
929 AliAODTrack *tmpTr = new AliAODTrack( -999, // id
940 kFALSE, // used for vertex fit
941 kFALSE, // used for prim vtx fit
942 AliAODTrack::kUndef, // type
943 fToyFilterMap, // select info
944 -999. // chi2 per NDF
946 tmpTr->SetFlags(AliESDtrack::kEmbedded);
948 new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
949 dummy = (*tracks)[nAODtracks-1];
951 fh1TrackPt->Fill(pt);
952 fh2TrackEtaPhi->Fill(eta,phi);
956 } //end: embed mode with toy
959 PostData(1, fHistList);
963 //__________________________________________________________________________
964 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
967 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
969 if(fAODfile && fAODfile->IsOpen())
974 //__________________________________________________________________________
975 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
977 // gives back the alien subjob id, if available, else -1
981 const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
982 //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
984 if(env && strlen(env)){
986 AliInfo(Form("Job index %d", id));
989 AliInfo("Job index not found. Okay if running locally.");
995 //__________________________________________________________________________
996 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
998 // select an AOD file from fAODPathArray
1000 Int_t nFiles = fAODPathArray->GetEntries();
1002 // choose randomly file and event
1004 Float_t tmpProp = -1.;
1005 Int_t pendingEvents = fInputEntries-Entry();
1006 //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1007 if(fAODEntriesArray){
1008 while(rndm->Rndm()>tmpProp){
1009 n = rndm->Integer(nFiles);
1010 fAODEntries = fAODEntriesArray->At(n);
1011 Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1012 tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1014 fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1016 AliWarning("Number of entries in extra AODs not set!");
1017 n = rndm->Integer(nFiles);
1021 AliInfo(Form("Select AOD file %d", n));
1022 TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1024 AliError("Could not get path of aod file.");
1027 fAODPath = objStr->GetString();
1032 //__________________________________________________________________________
1033 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1035 // select and open an AOD file from fAODPathArray
1037 if(fAODPathArray) fFileId = SelectAODfile();
1038 if(fFileId<0) return -1;
1040 TDirectory *owd = gDirectory;
1041 if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
1042 fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1049 AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1050 fFileId = OpenAODfile(trial+1);
1052 AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1056 AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1062 fAODtree = (TTree*)fAODfile->Get("aodTree");
1065 AliError("AOD tree not found.");
1070 fAODtree->SetBranchStatus("*",0);
1071 fAODtree->SetBranchStatus("header",1);
1072 fAODtree->SetBranchStatus("tracks*",1);
1073 fAODtree->SetBranchStatus("jets*",1);
1074 fAODtree->SetBranchStatus("clusters*",1);
1075 fAODtree->SetBranchStatus("mcparticles*",1);
1076 fAODtree->SetBranchStatus("mcHeader*",1);
1080 fAODevent = new AliAODEvent();
1081 fAODevent->ReadFromTree(fAODtree);
1084 // fetch header, jets, etc. from new file
1085 fNevents = fAODtree->GetEntries();
1086 mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1088 if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1089 else fAODJets = fAODevent->GetJets();
1091 AliError("Could not find jets in AOD. Check jet branch when indicated.");
1095 TFile *curfile = fAODtree->GetCurrentFile();
1097 AliError("No current file.");
1101 Float_t trials = 1.;
1103 PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1106 // construct a poor man average trials
1107 Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1108 if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1111 AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1114 fCountEvents=0; // new file, reset counter
1116 return fFileId; // file position in AOD path array, if array available
1120 //____________________________________________________________________________
1121 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1123 // static stored, available for other tasks in train
1125 static Float_t ptHard = -1.;
1126 if(bSet) ptHard = newValue;
1131 //____________________________________________________________________________
1132 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1134 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1136 const Int_t nBins = 10;
1137 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1140 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1147 //____________________________________________________________________________
1148 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1150 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1151 // This should provide the path to the AOD/ESD file
1152 // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks)
1154 TString file(currFile);
1158 if(file.Contains("root_archive.zip#")){
1159 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1160 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1161 file.Replace(pos+1,20,"");
1164 // not an archive take the basename....
1165 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1168 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
1170 // next trial fetch the histgram file
1171 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1173 // not a severe condition but inciate that we have no information
1174 AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1178 // find the tlist we want to be independtent of the name so use the Tkey
1179 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1184 TList *list = dynamic_cast<TList*>(key->ReadObj());
1189 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1190 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1193 } // no tree pyxsec.root
1195 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1201 Double_t xsection = 0;
1202 xtree->SetBranchAddress("xsection",&xsection);
1203 xtree->SetBranchAddress("ntrials",&ntrials);