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>
41 #include "AliAnalysisTaskFastEmbedding.h"
42 #include "AliAnalysisManager.h"
43 #include "AliESDEvent.h"
44 #include "AliESDtrack.h"
45 #include "AliAODEvent.h"
46 #include "AliAODTrack.h"
47 #include "AliAODJet.h"
48 #include "AliAODMCParticle.h"
49 #include "AliAODMCHeader.h"
50 #include "AliInputEventHandler.h"
54 ClassImp(AliAnalysisTaskFastEmbedding)
56 //__________________________________________________________________________
57 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
69 ,fAODPath("AliAOD.root")
73 ,fOfflineTrgMask(AliVEvent::kAny)
83 ,fTrackBranch("aodExtraTracks")
84 ,fMCparticlesBranch("aodExtraMCparticles")
93 ,fEvtSelMinJetEta(-999.)
94 ,fEvtSelMaxJetEta( 999.)
96 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
102 ,fToyDistributionTrackPt(0.)
103 ,fToyMinTrackEta(-.5)
106 ,fToyMaxTrackPhi(2*TMath::Pi())
117 ,fHistEvtSelection(0)
136 // default constructor
141 //__________________________________________________________________________
142 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
143 : AliAnalysisTaskSE(name)
154 ,fAODPath("AliAOD.root")
158 ,fOfflineTrgMask(AliVEvent::kAny)
167 ,fNInputTracksMax(-1)
168 ,fTrackBranch("aodExtraTracks")
169 ,fMCparticlesBranch("aodExtraMCparticles")
178 ,fEvtSelMinJetEta(-999.)
179 ,fEvtSelMaxJetEta( 999.)
180 ,fEvtSelMinJetPhi(0.)
181 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
183 ,fToyMinNbOfTracks(1)
184 ,fToyMaxNbOfTracks(1)
187 ,fToyDistributionTrackPt(0.)
188 ,fToyMinTrackEta(-.5)
191 ,fToyMaxTrackPhi(2*TMath::Pi())
202 ,fHistEvtSelection(0)
222 DefineOutput(1, TList::Class());
226 //__________________________________________________________________________
227 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding ©)
228 : AliAnalysisTaskSE()
230 ,fAODout(copy.fAODout)
231 ,fAODevent(copy.fAODevent)
232 ,fAODtree(copy.fAODtree)
233 ,fAODfile(copy.fAODfile)
234 ,mcHeader(copy.mcHeader)
236 ,fInputEntries(copy.fInputEntries)
237 ,fAODPathArray(copy.fAODPathArray)
238 ,fAODEntriesArray(copy.fAODEntriesArray)
239 ,fAODPath(copy.fAODPath)
240 ,fAODEntries(copy.fAODEntries)
241 ,fAODEntriesSum(copy.fAODEntriesSum)
242 ,fAODEntriesMax(copy.fAODEntriesMax)
243 ,fOfflineTrgMask(copy.fOfflineTrgMask)
244 ,fMinContribVtx(copy.fMinContribVtx)
245 ,fVtxZMin(copy.fVtxZMin)
246 ,fVtxZMax(copy.fVtxZMax)
247 ,fEvtClassMin(copy.fEvtClassMin)
248 ,fEvtClassMax(copy.fEvtClassMax)
249 ,fCentMin(copy.fCentMin)
250 ,fCentMax(copy.fCentMax)
251 ,fNInputTracksMin(copy.fNInputTracksMin)
252 ,fNInputTracksMax(copy.fNInputTracksMax)
253 ,fTrackBranch(copy.fTrackBranch)
254 ,fMCparticlesBranch(copy.fMCparticlesBranch)
255 ,fJetBranch(copy.fJetBranch)
256 ,fFileId(copy.fFileId)
257 ,fAODEntry(copy.fAODEntry)
258 ,fCountEvents(copy.fCountEvents)
259 ,fEmbedMode(copy.fEmbedMode)
260 ,fEvtSelecMode(copy.fEvtSelecMode)
261 ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
262 ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
263 ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
264 ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
265 ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
266 ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
267 ,fExtraEffPb(copy.fExtraEffPb)
268 ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
269 ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
270 ,fToyMinTrackPt(copy.fToyMinTrackPt)
271 ,fToyMaxTrackPt(copy.fToyMaxTrackPt)
272 ,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
273 ,fToyMinTrackEta(copy.fToyMinTrackEta)
274 ,fToyMaxTrackEta(copy.fToyMaxTrackEta)
275 ,fToyMinTrackPhi(copy.fToyMinTrackPhi)
276 ,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
277 ,fToyFilterMap(copy.fToyFilterMap)
278 ,fTrackFilterMap(copy.fTrackFilterMap)
279 ,fNPtHard(copy.fNPtHard)
280 ,fPtHard(copy.fPtHard)
281 ,fPtHardBin(copy.fPtHardBin)
282 ,fAODJets(copy.fAODJets)
283 ,fNevents(copy.fNevents)
284 ,fXsection(copy.fXsection)
285 ,fAvgTrials(copy.fAvgTrials)
286 ,fHistList(copy.fHistList)
287 ,fHistEvtSelection(copy.fHistEvtSelection)
288 ,fh1Xsec(copy.fh1Xsec)
289 ,fh1Trials(copy.fh1Trials)
290 ,fh1TrialsEvtSel(copy.fh1TrialsEvtSel)
291 ,fh2PtHard(copy.fh2PtHard)
292 ,fh2PtHardEvtSel(copy.fh2PtHardEvtSel)
293 ,fh2PtHardTrials(copy.fh2PtHardTrials)
294 ,fh1TrackPt(copy.fh1TrackPt)
295 ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
296 ,fh1TrackN(copy.fh1TrackN)
297 ,fh1JetPt(copy.fh1JetPt)
298 ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
299 ,fh1JetN(copy.fh1JetN)
300 ,fh1MCTrackPt(copy.fh1MCTrackPt)
301 ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
302 ,fh1MCTrackN(copy.fh1MCTrackN)
303 ,fh1AODfile(copy.fh1AODfile)
304 ,fh2AODevent(copy.fh2AODevent)
310 //__________________________________________________________________________
311 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
316 AliAnalysisTaskSE::operator=(o);
319 fAODevent = o.fAODevent;
320 fAODtree = o.fAODtree;
321 fAODfile = o.fAODfile;
322 mcHeader = o.mcHeader;
324 fInputEntries = o.fInputEntries;
325 fAODPathArray = o.fAODPathArray;
326 fAODEntriesArray = o.fAODEntriesArray;
327 fAODPath = o.fAODPath;
328 fAODEntries = o.fAODEntries;
329 fAODEntriesSum = o.fAODEntriesSum;
330 fAODEntriesMax = o.fAODEntriesMax;
331 fOfflineTrgMask = o.fOfflineTrgMask;
332 fMinContribVtx = o.fMinContribVtx;
333 fVtxZMin = o.fVtxZMin;
334 fVtxZMax = o.fVtxZMax;
335 fEvtClassMin = o.fEvtClassMin;
336 fEvtClassMax = o.fEvtClassMax;
337 fCentMin = o.fCentMin;
338 fCentMax = o.fCentMax;
339 fNInputTracksMin = o.fNInputTracksMin;
340 fNInputTracksMax = o.fNInputTracksMax;
341 fTrackBranch = o.fTrackBranch;
342 fMCparticlesBranch = o.fMCparticlesBranch;
343 fJetBranch = o.fJetBranch;
345 fAODEntry = o.fAODEntry;
346 fCountEvents = o.fCountEvents;
347 fEmbedMode = o.fEmbedMode;
348 fEvtSelecMode = o.fEvtSelecMode;
349 fEvtSelMinJetPt = o.fEvtSelMinJetPt;
350 fEvtSelMaxJetPt = o.fEvtSelMaxJetPt;
351 fEvtSelMinJetEta = o.fEvtSelMinJetEta;
352 fEvtSelMaxJetEta = o.fEvtSelMaxJetEta;
353 fEvtSelMinJetPhi = o.fEvtSelMinJetPhi;
354 fEvtSelMaxJetPhi = o.fEvtSelMaxJetPhi;
355 fExtraEffPb = o.fExtraEffPb;
356 fToyMinNbOfTracks = o.fToyMinNbOfTracks;
357 fToyMaxNbOfTracks = o.fToyMaxNbOfTracks;
358 fToyMinTrackPt = o.fToyMinTrackPt;
359 fToyMaxTrackPt = o.fToyMaxTrackPt;
360 fToyDistributionTrackPt = o.fToyDistributionTrackPt;
361 fToyMinTrackEta = o.fToyMinTrackEta;
362 fToyMaxTrackEta = o.fToyMaxTrackEta;
363 fToyMinTrackPhi = o.fToyMinTrackPhi;
364 fToyMaxTrackPhi = o.fToyMaxTrackPhi;
365 fToyFilterMap = o.fToyFilterMap;
366 fTrackFilterMap = o.fTrackFilterMap;
367 fNPtHard = o.fNPtHard;
369 fPtHardBin = o.fPtHardBin;
370 fAODJets = o.fAODJets;
371 fNevents = o.fNevents;
372 fXsection = o.fXsection;
373 fAvgTrials = o.fAvgTrials;
374 fHistList = o.fHistList;
375 fHistEvtSelection = o.fHistEvtSelection;
377 fh1Trials = o.fh1Trials;
378 fh1TrialsEvtSel = o.fh1TrialsEvtSel;
379 fh2PtHard = o.fh2PtHard;
380 fh2PtHardEvtSel = o.fh2PtHardEvtSel;
381 fh2PtHardTrials = o.fh2PtHardTrials;
382 fh1TrackPt = o.fh1TrackPt;
383 fh2TrackEtaPhi = o.fh2TrackEtaPhi;
384 fh1TrackN = o.fh1TrackN;
385 fh1JetPt = o.fh1JetPt;
386 fh2JetEtaPhi = o.fh2JetEtaPhi;
388 fh1MCTrackPt = o.fh1MCTrackPt;
389 fh2MCTrackEtaPhi = o.fh2MCTrackEtaPhi;
390 fh1MCTrackN = o.fh1MCTrackN;
391 fh1AODfile = o.fh1AODfile;
392 fh2AODevent = o.fh2AODevent;
399 //__________________________________________________________________________
400 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
407 //__________________________________________________________________________
408 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
410 // create output objects
411 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
412 AliLog::SetClassDebugLevel("AliAnalysisTaskFastEmbedding", AliLog::kInfo);
415 if(!fHistList) fHistList = new TList();
416 fHistList->SetOwner(kTRUE);
420 rndm = new TRandom3();
421 Int_t id = GetJobID();
422 if(id>-1) rndm->SetSeed(id);
423 else rndm->SetSeed(); // a TTUID is generated and used for seed
424 AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
428 // embed mode with AOD
429 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
432 fFileId = OpenAODfile();
433 fNevents = 0; // force to open another aod in UserExec()
436 PostData(1, fHistList);
439 } //end: embed mode with AOD
442 // connect output aod
443 // create a new branch for extra tracks
444 fAODout = AODEvent();
446 AliError("Output AOD not found.");
447 PostData(1, fHistList);
450 if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
451 AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
452 TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
453 tracks->SetName(fTrackBranch.Data());
454 AddAODBranch("TClonesArray", &tracks);
456 // create new branch for extra mcparticle if available as input
457 if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
458 AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
459 TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
460 mcparticles->SetName(fMCparticlesBranch.Data());
461 AddAODBranch("TClonesArray", &mcparticles);
467 Bool_t oldStatus = TH1::AddDirectoryStatus();
468 TH1::AddDirectory(kFALSE);
470 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
471 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
472 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
473 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
474 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
475 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
476 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
478 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root;p_{T,hard} bin;<#sigma>",fNPtHard+1,-1.5,fNPtHard-0.5);
479 fh1Trials = new TH1F("fh1Trials","trials (simulation) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
480 fh1TrialsEvtSel = new TH1F("fh1TrialsEvtSel","trials (event selection) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
481 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);
482 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);
483 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);
485 fHistList->Add(fHistEvtSelection);
486 fHistList->Add(fh1Xsec);
487 fHistList->Add(fh1Trials);
488 fHistList->Add(fh1TrialsEvtSel);
489 fHistList->Add(fh2PtHard);
490 fHistList->Add(fh2PtHardEvtSel);
491 fHistList->Add(fh2PtHardTrials);
493 fh1TrackPt = new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
494 fh2TrackEtaPhi = new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
495 fh1TrackN = new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
497 fHistList->Add(fh1TrackPt);
498 fHistList->Add(fh2TrackEtaPhi);
499 fHistList->Add(fh1TrackN);
501 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
503 fh1JetPt = new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 120, 0., 120.);
504 fh2JetEtaPhi = new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
505 fh1JetN = new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
507 fHistList->Add(fh1JetPt);
508 fHistList->Add(fh2JetEtaPhi);
509 fHistList->Add(fh1JetN);
513 if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
515 fh1MCTrackPt = new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
516 fh2MCTrackEtaPhi = new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
517 fh1MCTrackN = new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
519 fHistList->Add(fh1MCTrackPt);
520 fHistList->Add(fh2MCTrackEtaPhi);
521 fHistList->Add(fh1MCTrackN);
525 fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 23, -0.5, 2299.5);
526 fh2AODevent = new TH2I("fh1AODevent","selected events;file;event", 25,-.5,2499.5,50,-.5,4999.5);
527 fHistList->Add(fh1AODfile);
528 fHistList->Add(fh2AODevent);
530 // =========== Switch on Sumw2 for all histos ===========
531 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
532 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
539 TH1::AddDirectory(oldStatus);
541 PostData(1, fHistList);
545 //__________________________________________________________________________
546 void AliAnalysisTaskFastEmbedding::Init()
549 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
553 //__________________________________________________________________________
554 Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
556 // User defined Notify(), called once
557 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserNotify()");
559 // get total nb of events in tree (of this subjob)
560 AliInputEventHandler* inputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
561 fInputEntries = inputHandler->GetTree()->GetEntriesFast();
562 AliInfo(Form("Total nb. of events: %d", fInputEntries));
569 //__________________________________________________________________________
570 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
572 // Analysis, called once per event
573 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
576 AliError("Need output AOD, but is not connected.");
577 PostData(1, fHistList);
581 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
583 AliError("ESD not available");
584 PostData(1, fHistList);
588 // -- event selection --
589 fHistEvtSelection->Fill(1); // number of events before event selection
592 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
593 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
594 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
595 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
596 fHistEvtSelection->Fill(2);
597 PostData(1, fHistList);
602 AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
603 Int_t nTracksPrim = primVtx->GetNContributors();
604 if ((nTracksPrim < fMinContribVtx) ||
605 (primVtx->GetZ() < fVtxZMin) ||
606 (primVtx->GetZ() > fVtxZMax) ){
607 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
608 fHistEvtSelection->Fill(3);
609 PostData(1, fHistList);
613 /** takes wrong value, since jet service tasks runs later
614 // event class selection (from jet helper task)
615 Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
616 if(fDebug) Printf("Event class %d", eventClass);
617 if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
618 fHistEvtSelection->Fill(4);
619 PostData(1, fHistList);
623 // centrality selection
624 AliCentrality *cent = 0x0;
625 Float_t centValue = 0.;
626 if(fESD) cent = fESD->GetCentrality();
627 if(cent) centValue = cent->GetCentralityPercentile("V0M");
628 if(fDebug) printf("centrality: %f\n", centValue);
629 if (centValue < fCentMin || centValue > fCentMax){
630 fHistEvtSelection->Fill(4);
631 PostData(1, fHistList);
636 /* ** not implemented **
637 // multiplicity due to input tracks
638 Int_t nInputTracks = GetNInputTracks();
640 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
641 fHistEvtSelection->Fill(5);
642 PostData(1, fHistList);
647 fHistEvtSelection->Fill(0); // accepted events
648 // -- end event selection --
651 TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
653 AliError("Extra track branch not found in output.");
654 PostData(1, fHistList);
662 // === embed mode with AOD ===
663 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
665 AliError("Need input AOD, but is not connected.");
666 PostData(1, fHistList);
672 Bool_t useEntry = kFALSE;
673 while(!useEntry){ // protection need, if no event fulfills requirement
675 fAODEntry++; // go to next event
677 if(fAODEntry>=fNevents) fAODEntry=0;
680 if(fCountEvents>=fNevents){
682 AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
685 AliDebug(AliLog::kDebug, "Select new AOD file ...");
687 fFileId = OpenAODfile();
689 PostData(1, fHistList);
692 fh1AODfile->Fill(fFileId);
693 if(fAODEntries!=fNevents){
694 AliError("File with list of AODs and file with nb. of entries inconsistent");
695 PostData(1, fHistList);
704 fAODtree->GetEvent(fAODEntry);
708 fPtHard = mcHeader->GetPtHard();
709 GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
710 AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
713 AliDebug(AliLog::kDebug,"No mcHeader");
715 fPtHardBin = GetPtHardBin(fPtHard);
717 //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials);
719 // fill event variables for each event
720 fh1Xsec->Fill(fPtHardBin,fXsection);
721 fh2PtHard->Fill(fPtHardBin,fPtHard);
723 fh1Trials->Fill(fPtHardBin, fAvgTrials);
724 fh1TrialsEvtSel->Fill(fPtHardBin);
725 fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
728 if(fEvtSelecMode==kEventsJetPt){
729 Int_t nJets = fAODJets->GetEntries();
730 //Printf("nb. jets: %d",nJets);
731 for(Int_t ij=0; ij<nJets; ++ij){
732 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
735 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
736 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
737 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
738 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
746 if(fEvtSelecMode==kEventsAll){
751 AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
753 fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
754 fh2AODevent->Fill(fFileId,fAODEntry);
756 TClonesArray *mcpartIN = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
757 TClonesArray *mcpartOUT = 0x0;
759 mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
761 AliError("Extra MC particles branch not found in output.");
762 PostData(1, fHistList);
767 AliInfo("No extra MC particles found.");
771 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
772 // loop over input tracks
775 Int_t nJets = fAODJets->GetEntries();
776 Int_t nSelectedJets = 0;
777 AliAODJet *leadJet = 0x0; // used for jet tracks only
779 if(fEmbedMode==kAODFull){
780 nTracks = fAODevent->GetNumberOfTracks();
782 for(Int_t ij=0; ij<nJets; ++ij){
783 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
785 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
786 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
787 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
788 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
790 fh1JetPt->Fill(jet->Pt());
791 fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
796 fh1JetN->Fill(nSelectedJets);
799 if(fEmbedMode==kAODJetTracks){
800 // look for leading jet within selection
801 // get reference tracks for this jet
802 Float_t leadJetPt = 0.;
803 for(Int_t ij=0; ij<nJets; ++ij){
804 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
806 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
807 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
808 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
809 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
810 if(jet->Pt()>leadJetPt){
811 leadJetPt = jet->Pt();
817 nTracks = leadJet->GetRefTracks()->GetEntriesFast();
818 fh1JetPt->Fill(leadJet->Pt());
819 fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
824 fh1TrackN->Fill((Float_t)nTracks);
826 for(Int_t it=0; it<nTracks; ++it){
827 AliAODTrack *tmpTr = 0x0;
828 if(fEmbedMode==kAODFull) tmpTr = fAODevent->GetTrack(it);
829 if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
831 Double_t rd=rndm->Uniform(0.,1.);
832 if(rd>fExtraEffPb) continue;
833 tmpTr->SetStatus(AliESDtrack::kEmbedded);
835 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
836 dummy = (*tracks)[nAODtracks-1];
838 if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
839 fh1TrackPt->Fill(tmpTr->Pt());
840 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
845 Int_t nMCpart = mcpartIN->GetEntriesFast();
847 Int_t nPhysicalPrimary=0;
849 for(Int_t ip=0; ip<nMCpart; ++ip){
850 AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
851 if(!tmpPart) continue;
852 if(fEmbedMode==kAODJetTracks){
853 // jet track? else continue
854 // not implemented yet
858 if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. && tmpPart->Pt()>0.){
859 new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
860 dummy = (*mcpartOUT)[nAODmcpart-1];
862 if(fDebug>10) printf("added track %d with pT=%.2f to extra branch\n",nAODmcpart,tmpPart->Pt());
864 fh1MCTrackPt->Fill(tmpPart->Pt());
865 fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
870 fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
873 } // end: embed all tracks or jet tracks
876 if(fEmbedMode==kAODJet4Mom){
879 Int_t nJets = fAODJets->GetEntries();
880 fh1TrackN->Fill((Float_t)nJets);
881 for(Int_t ij=0; ij<nJets; ++ij){
882 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
884 AliAODTrack *tmpTr = (AliAODTrack*)jet;
885 tmpTr->SetFlags(AliESDtrack::kEmbedded);
887 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
888 dummy = (*tracks)[nAODtracks-1];
890 fh1TrackPt->Fill(tmpTr->Pt());
891 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
894 } // end: embed jets as 4-momenta
897 } //end: embed mode with AOD
900 // === embed mode with toy ===
901 if(fEmbedMode==kToyTracks){
902 Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
904 fh1TrackN->Fill((Float_t)nT);
906 for(Int_t i=0; i<nT; ++i){
909 if(fToyMinTrackPt!=fToyMaxTrackPt){
910 if(fToyDistributionTrackPt==0){
911 pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
913 while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
914 pt = rndm->Exp(fToyDistributionTrackPt); // no power law yet!!
915 pt += fToyMinTrackPt;
921 Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
922 Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
923 phi = TVector2::Phi_0_2pi(phi);
925 if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
927 Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
938 AliAODTrack *tmpTr = new AliAODTrack( -999, // id
949 kFALSE, // used for vertex fit
950 kFALSE, // used for prim vtx fit
951 AliAODTrack::kUndef, // type
952 fToyFilterMap, // select info
953 -999. // chi2 per NDF
955 tmpTr->SetFlags(AliESDtrack::kEmbedded);
957 new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
958 dummy = (*tracks)[nAODtracks-1];
960 fh1TrackPt->Fill(pt);
961 fh2TrackEtaPhi->Fill(eta,phi);
965 } //end: embed mode with toy
968 PostData(1, fHistList);
972 //__________________________________________________________________________
973 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
976 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
978 if(fAODfile && fAODfile->IsOpen())
983 //__________________________________________________________________________
984 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
986 // gives back the alien subjob id, if available, else -1
990 const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
991 //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
993 if(env && strlen(env)){
995 AliInfo(Form("Job index %d", id));
998 AliInfo("Job index not found. Okay if running locally.");
1004 //__________________________________________________________________________
1005 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
1007 // select an AOD file from fAODPathArray
1009 Int_t nFiles = fAODPathArray->GetEntries();
1011 // choose randomly file and event
1013 Float_t tmpProp = -1.;
1014 Int_t pendingEvents = fInputEntries-Entry();
1015 //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1016 if(fAODEntriesArray){
1017 while(rndm->Rndm()>tmpProp){
1018 n = rndm->Integer(nFiles);
1019 fAODEntries = fAODEntriesArray->At(n);
1020 Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1021 tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1023 fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1026 AliWarning("Number of entries in extra AODs not set!");
1027 n = rndm->Integer(nFiles);
1031 AliInfo(Form("Select AOD file %d", n));
1032 TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1034 AliError("Could not get path of aod file.");
1037 fAODPath = objStr->GetString();
1042 //__________________________________________________________________________
1043 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1045 // select and open an AOD file from fAODPathArray
1047 if(fAODPathArray) fFileId = SelectAODfile();
1048 if(fFileId<0) return -1;
1051 AliInfo("Trying to connect to AliEn ...");
1052 TGrid::Connect("alien://");
1055 TDirectory *owd = gDirectory;
1056 if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
1057 fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1064 AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1065 fFileId = OpenAODfile(trial+1);
1067 AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1071 AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1077 fAODtree = (TTree*)fAODfile->Get("aodTree");
1080 AliError("AOD tree not found.");
1085 fAODtree->SetBranchStatus("*",0);
1086 fAODtree->SetBranchStatus("header",1);
1087 fAODtree->SetBranchStatus("tracks*",1);
1088 fAODtree->SetBranchStatus("jets*",1);
1089 fAODtree->SetBranchStatus("clusters*",1);
1090 fAODtree->SetBranchStatus("mcparticles*",1);
1091 fAODtree->SetBranchStatus("mcHeader*",1);
1095 fAODevent = new AliAODEvent();
1096 fAODevent->ReadFromTree(fAODtree);
1099 // fetch header, jets, etc. from new file
1100 fNevents = fAODtree->GetEntries();
1101 mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1103 if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1104 else fAODJets = fAODevent->GetJets();
1106 AliError("Could not find jets in AOD. Check jet branch when indicated.");
1110 TFile *curfile = fAODtree->GetCurrentFile();
1112 AliError("No current file.");
1116 Float_t trials = 1.;
1118 PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1121 // construct a poor man average trials
1122 Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1123 if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1126 AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1129 fCountEvents=0; // new file, reset counter
1131 return fFileId; // file position in AOD path array, if array available
1135 //____________________________________________________________________________
1136 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1138 // static stored, available for other tasks in train
1140 static Float_t ptHard = -1.;
1141 if(bSet) ptHard = newValue;
1146 //____________________________________________________________________________
1147 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1149 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1151 const Int_t nBins = 10;
1152 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1155 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1162 //____________________________________________________________________________
1163 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1165 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1166 // This should provide the path to the AOD/ESD file
1167 // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks)
1169 TString file(currFile);
1173 if(file.Contains("root_archive.zip#")){
1174 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
1175 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1176 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
1177 file.Replace(pos+1,pos2-pos1,"");
1178 // file.Replace(pos+1,20,"");
1181 // not an archive take the basename....
1182 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1185 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
1187 // next trial fetch the histgram file
1188 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1190 // not a severe condition but inciate that we have no information
1191 AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1195 // find the tlist we want to be independtent of the name so use the Tkey
1196 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1201 TList *list = dynamic_cast<TList*>(key->ReadObj());
1206 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1207 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1210 } // no tree pyxsec.root
1212 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1218 Double_t xsection = 0;
1219 xtree->SetBranchAddress("xsection",&xsection);
1220 xtree->SetBranchAddress("ntrials",&ntrials);