set minimum pT for dicing (M. Verweij)
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskFastEmbedding.cxx
CommitLineData
b725dccf 1/*************************************************************************
2 * *
3 * Task for fast embedding *
4 * read extra input from AOD *
5 * *
6 *************************************************************************/
7
8
9/**************************************************************************
10 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11 * *
12 * Author: The ALICE Off-line Project. *
13 * Contributors are mentioned in the code where appropriate. *
14 * *
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 **************************************************************************/
23
24/* $Id: */
25
26#include <TFile.h>
27#include <TTree.h>
28#include <TChain.h>
29#include <TClonesArray.h>
30#include <TDirectory.h>
31#include <TSystem.h>
32#include <TRef.h>
33#include <TRandom3.h>
34#include <TH1F.h>
35#include <TH2F.h>
31b9d515 36#include <TProfile.h>
37#include <TKey.h>
b725dccf 38
39
40#include "AliAnalysisTaskFastEmbedding.h"
41#include "AliAnalysisManager.h"
31b9d515 42#include "AliESDEvent.h"
46465e39 43#include "AliESDtrack.h"
b725dccf 44#include "AliAODEvent.h"
45#include "AliAODTrack.h"
46#include "AliAODJet.h"
47#include "AliAODMCParticle.h"
31b9d515 48#include "AliAODMCHeader.h"
49#include "AliInputEventHandler.h"
b725dccf 50
51#include "AliLog.h"
52
53ClassImp(AliAnalysisTaskFastEmbedding)
54
55//__________________________________________________________________________
56AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
31b9d515 57: AliAnalysisTaskSE()
58,fESD(0)
59,fAODout(0)
60,fAODevent(0)
61,fAODtree(0)
62,fAODfile(0)
63,mcHeader(0)
64,rndm(0)
65,fInputEntries(0)
66,fAODPathArray(0)
67,fAODEntriesArray(0)
68,fAODPath("AliAOD.root")
69,fAODEntries(-1)
70,fAODEntriesSum(0)
71,fAODEntriesMax(0)
72,fOfflineTrgMask(AliVEvent::kAny)
73,fMinContribVtx(1)
74,fVtxZMin(-8.)
75,fVtxZMax(8.)
76,fEvtClassMin(0)
77,fEvtClassMax(4)
78,fCentMin(0.)
79,fCentMax(100.)
80,fNInputTracksMin(0)
81,fNInputTracksMax(-1)
82,fTrackBranch("aodExtraTracks")
83,fMCparticlesBranch("aodExtraMCparticles")
84,fJetBranch("")
85,fFileId(-1)
86,fAODEntry(0)
87,fCountEvents(-1)
88,fEmbedMode(0)
89,fEvtSelecMode(0)
90,fEvtSelMinJetPt(-1)
91,fEvtSelMaxJetPt(-1)
92,fEvtSelMinJetEta(-999.)
93,fEvtSelMaxJetEta( 999.)
94,fEvtSelMinJetPhi(0.)
95,fEvtSelMaxJetPhi(TMath::Pi()*2.)
96,fToyMinNbOfTracks(1)
97,fToyMaxNbOfTracks(1)
98,fToyMinTrackPt(50.)
99,fToyMaxTrackPt(50.)
100,fToyDistributionTrackPt(0.)
101,fToyMinTrackEta(-.5)
102,fToyMaxTrackEta(.5)
103,fToyMinTrackPhi(0.)
104,fToyMaxTrackPhi(2*TMath::Pi())
105,fToyFilterMap(0)
106,fTrackFilterMap(0)
107,fNPtHard(10)
108,fPtHard(0)
109,fPtHardBin(-1)
110,fAODJets(0x0)
111,fNevents(0)
112,fXsection(0)
113,fAvgTrials(0)
114,fHistList(0)
115,fHistEvtSelection(0)
116,fh1Xsec(0)
117,fh1Trials(0)
118,fh1TrialsEvtSel(0)
119,fh2PtHard(0)
120,fh2PtHardEvtSel(0)
121,fh2PtHardTrials(0)
122,fh1TrackPt(0)
123,fh2TrackEtaPhi(0)
124,fh1TrackN(0)
125,fh1JetPt(0)
126,fh2JetEtaPhi(0)
127,fh1JetN(0)
128,fh1MCTrackPt(0)
129,fh2MCTrackEtaPhi(0)
130,fh1MCTrackN(0)
131,fh1AODfile(0)
132,fh2AODevent(0)
b725dccf 133{
31b9d515 134 // default constructor
b725dccf 135
136}
137
138
139//__________________________________________________________________________
140AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
31b9d515 141: AliAnalysisTaskSE(name)
142,fESD(0)
143,fAODout(0)
144,fAODevent(0)
145,fAODtree(0)
146,fAODfile(0)
147,mcHeader(0)
148,rndm(0)
149,fInputEntries(0)
150,fAODPathArray(0)
151,fAODEntriesArray(0)
152,fAODPath("AliAOD.root")
153,fAODEntries(-1)
154,fAODEntriesSum(0)
155,fAODEntriesMax(0)
156,fOfflineTrgMask(AliVEvent::kAny)
157,fMinContribVtx(1)
158,fVtxZMin(-8.)
159,fVtxZMax(8.)
160,fEvtClassMin(0)
161,fEvtClassMax(4)
162,fCentMin(0.)
163,fCentMax(100.)
164,fNInputTracksMin(0)
165,fNInputTracksMax(-1)
166,fTrackBranch("aodExtraTracks")
167,fMCparticlesBranch("aodExtraMCparticles")
168,fJetBranch("")
169,fFileId(-1)
170,fAODEntry(0)
171,fCountEvents(-1)
172,fEmbedMode(0)
173,fEvtSelecMode(0)
174,fEvtSelMinJetPt(-1)
175,fEvtSelMaxJetPt(-1)
176,fEvtSelMinJetEta(-999.)
177,fEvtSelMaxJetEta( 999.)
178,fEvtSelMinJetPhi(0.)
179,fEvtSelMaxJetPhi(TMath::Pi()*2.)
180,fToyMinNbOfTracks(1)
181,fToyMaxNbOfTracks(1)
182,fToyMinTrackPt(50.)
183,fToyMaxTrackPt(50.)
184,fToyDistributionTrackPt(0.)
185,fToyMinTrackEta(-.5)
186,fToyMaxTrackEta(.5)
187,fToyMinTrackPhi(0.)
188,fToyMaxTrackPhi(2*TMath::Pi())
189,fToyFilterMap(0)
190,fTrackFilterMap(0)
191,fNPtHard(10)
192,fPtHard(0)
193,fPtHardBin(-1)
194,fAODJets(0x0)
195,fNevents(0)
196,fXsection(0)
197,fAvgTrials(0)
198,fHistList(0)
199,fHistEvtSelection(0)
200,fh1Xsec(0)
201,fh1Trials(0)
202,fh1TrialsEvtSel(0)
203,fh2PtHard(0)
204,fh2PtHardEvtSel(0)
205,fh2PtHardTrials(0)
206,fh1TrackPt(0)
207,fh2TrackEtaPhi(0)
208,fh1TrackN(0)
209,fh1JetPt(0)
210,fh2JetEtaPhi(0)
211,fh1JetN(0)
212,fh1MCTrackPt(0)
213,fh2MCTrackEtaPhi(0)
214,fh1MCTrackN(0)
215,fh1AODfile(0)
216,fh2AODevent(0)
b725dccf 217{
31b9d515 218 // constructor
219 DefineOutput(1, TList::Class());
b725dccf 220}
221
222
223//__________________________________________________________________________
224AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy)
31b9d515 225: AliAnalysisTaskSE()
226,fESD(copy.fESD)
227,fAODout(copy.fAODout)
228,fAODevent(copy.fAODevent)
229,fAODtree(copy.fAODtree)
230,fAODfile(copy.fAODfile)
231,mcHeader(copy.mcHeader)
232,rndm(copy.rndm)
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)
b725dccf 301{
31b9d515 302 // copy constructor
b725dccf 303}
304
305
306//__________________________________________________________________________
307AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
308{
31b9d515 309 // assignment
310
311 if(this!=&o){
312 AliAnalysisTaskSE::operator=(o);
313 fESD = o.fESD;
314 fAODout = o.fAODout;
315 fAODevent = o.fAODevent;
316 fAODtree = o.fAODtree;
317 fAODfile = o.fAODfile;
318 mcHeader = o.mcHeader;
319 rndm = o.rndm;
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;
340 fFileId = o.fFileId;
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;
363 fPtHard = o.fPtHard;
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;
371 fh1Xsec = o.fh1Xsec;
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;
382 fh1JetN = o.fh1JetN;
383 fh1MCTrackPt = o.fh1MCTrackPt;
384 fh2MCTrackEtaPhi = o.fh2MCTrackEtaPhi;
385 fh1MCTrackN = o.fh1MCTrackN;
386 fh1AODfile = o.fh1AODfile;
387 fh2AODevent = o.fh2AODevent;
388 }
389
390 return *this;
b725dccf 391}
392
393
394//__________________________________________________________________________
395AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
396{
31b9d515 397 // destructor
398 delete rndm;
b725dccf 399}
400
401
402//__________________________________________________________________________
403void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
404{
31b9d515 405 // create output objects
406 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
407 AliLog::SetClassDebugLevel("AliAnalysisTaskFastEmbedding", AliLog::kInfo);
408
409 OpenFile(1);
410 if(!fHistList) fHistList = new TList();
411 fHistList->SetOwner(kTRUE);
412
413
414 // set seed
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()));
420
421
422
423 // embed mode with AOD
424 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
425
426 // open input AOD
427 fFileId = OpenAODfile();
428 fNevents = 0; // force to open another aod in UserExec()
429 if(fFileId<0){
430 AliError("");
431 PostData(1, fHistList);
432 return;
433 }
434 } //end: embed mode with AOD
435
436
437 // connect output aod
438 // create a new branch for extra tracks
439 fAODout = AODEvent();
440 if(!fAODout){
441 AliError("Output AOD not found.");
442 PostData(1, fHistList);
443 return;
444 }
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);
450 }
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);
457 }
458
459
460
461 //qa histograms
462 Bool_t oldStatus = TH1::AddDirectoryStatus();
463 TH1::AddDirectory(kFALSE);
464
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)");
472
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);
479
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);
487
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.);
491
492 fHistList->Add(fh1TrackPt);
493 fHistList->Add(fh2TrackEtaPhi);
494 fHistList->Add(fh1TrackN);
495
496 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
497
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.);
501
502 fHistList->Add(fh1JetPt);
503 fHistList->Add(fh2JetEtaPhi);
504 fHistList->Add(fh1JetN);
505 }
506
507
508 if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
509
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.);
513
514 fHistList->Add(fh1MCTrackPt);
515 fHistList->Add(fh2MCTrackEtaPhi);
516 fHistList->Add(fh1MCTrackN);
517
518 }
519
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);
524
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));
528 if (h1){
529 h1->Sumw2();
530 continue;
531 }
532 }
533
534 TH1::AddDirectory(oldStatus);
535
536 PostData(1, fHistList);
b725dccf 537}
538
539
540//__________________________________________________________________________
541void AliAnalysisTaskFastEmbedding::Init()
542{
31b9d515 543 // Initialization
544 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
545
546}
547
548//__________________________________________________________________________
549Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
550{
7a88ef30 551 // User defined Notify(), called once
31b9d515 552 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserNotify()");
553
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));
558
559 return kTRUE;
b725dccf 560
b725dccf 561}
562
563
564//__________________________________________________________________________
565void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
566{
7a88ef30 567 // Analysis, called once per event
31b9d515 568 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
569
570 if(!fAODout){
571 AliError("Need output AOD, but is not connected.");
572 PostData(1, fHistList);
573 return;
574 }
575
576 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
577 if (!fESD) {
578 AliError("ESD not available");
579 PostData(1, fHistList);
580 return;
581 }
582
583 // -- event selection --
584 fHistEvtSelection->Fill(1); // number of events before event selection
585
586 // physics 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);
593 return;
594 }
595
596 // vertex selection
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);
605 return;
606 }
607
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);
615 return;
616 }*/
617
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);
627 return;
628 }
629
630
631 /* ** not implemented **
632 // multiplicity due to input tracks
633 Int_t nInputTracks = GetNInputTracks();
634
635 if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
636 fHistEvtSelection->Fill(5);
637 PostData(1, fHistList);
638 return;
639 }
640 */
641
642 fHistEvtSelection->Fill(0); // accepted events
643 // -- end event selection --
644
645 // connect aod out
646 TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
647 if(!tracks){
648 AliError("Extra track branch not found in output.");
649 PostData(1, fHistList);
650 return;
651 }
652 tracks->Delete();
653 Int_t nAODtracks=0;
654
655 TRef dummy;
656
657 // === embed mode with AOD ===
658 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
659 if(!fAODevent){
660 AliError("Need input AOD, but is not connected.");
aa117124 661 PostData(1, fHistList);
662 return;
31b9d515 663 }
664
665
666
667 Bool_t useEntry = kFALSE;
668 while(!useEntry){ // protection need, if no event fulfills requierment
669
670 fAODEntry++; // go to next event
671 fCountEvents++;
672 if(fAODEntry>=fNevents) fAODEntry=0;
673
674 // open new aod file
675 if(fCountEvents>=fNevents){
676 if(!fAODPathArray){
677 AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
678 }
679 else {
680 AliDebug(AliLog::kDebug, "Select new AOD file ...");
681
682 fFileId = OpenAODfile();
683 if(fFileId<0) {
684 PostData(1, fHistList);
685 return;
686 }
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);
691 return;
692 }
693 }
694 }
695
696
697
698 // get next event
699 fAODtree->GetEvent(fAODEntry);
700
701 // get pt hard
702 if(mcHeader){
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));
706 }
707 else{
708 AliDebug(AliLog::kDebug,"No mcHeader");
709 }
710 fPtHardBin = GetPtHardBin(fPtHard);
711
712 //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials);
713
714 // fill event variables for each event
715 fh1Xsec->Fill(fPtHardBin,fXsection);
716 fh2PtHard->Fill(fPtHardBin,fPtHard);
717
718 fh1Trials->Fill(fPtHardBin, fAvgTrials);
719 fh1TrialsEvtSel->Fill(fPtHardBin);
720 fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
721
722 // jet pt selection
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));
728 if(!jet) continue;
729
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)){
734 useEntry = kTRUE;
735 break;
736 }
737 }
738 }
739
740 // no selection
741 if(fEvtSelecMode==kEventsAll){
742 useEntry = kTRUE;
743 }
744
745 }
746 AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
747
748 fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
749 fh2AODevent->Fill(fFileId,fAODEntry);
750
751 TClonesArray *mcpartIN = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
752 TClonesArray *mcpartOUT = 0x0;
753 if(mcpartIN){
754 mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
755 if(!mcpartOUT){
756 AliError("Extra MC particles branch not found in output.");
757 PostData(1, fHistList);
758 return;
759 }
760 mcpartOUT->Delete();
761 } else {
762 AliInfo("No extra MC particles found.");
763 }
764
765
766 if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
767 // loop over input tracks
768 // add to output aod
769 Int_t nTracks = 0;
770 Int_t nJets = fAODJets->GetEntries();
771 Int_t nSelectedJets = 0;
772 AliAODJet *leadJet = 0x0; // used for jet tracks only
773
774 if(fEmbedMode==kAODFull){
775 nTracks = fAODevent->GetNumberOfTracks();
776
777 for(Int_t ij=0; ij<nJets; ++ij){
778 AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
779 if(!jet) continue;
780 if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
46465e39 781 && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
782 && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
783 && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
b725dccf 784
31b9d515 785 fh1JetPt->Fill(jet->Pt());
786 fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
787 nSelectedJets++;
788
789 }
790 }
791 fh1JetN->Fill(nSelectedJets);
792 }
793
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));
1a2bb6d5 800 if(!jet) continue;
31b9d515 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();
807 leadJet = jet;
808 }
809 }
810 }
811 if(leadJet){
812 nTracks = leadJet->GetRefTracks()->GetEntriesFast();
813 fh1JetPt->Fill(leadJet->Pt());
814 fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
815 fh1JetN->Fill(1);
816 }
817 }
b725dccf 818
31b9d515 819 fh1TrackN->Fill((Float_t)nTracks);
b725dccf 820
31b9d515 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));
825 if(!tmpTr) continue;
826
827 tmpTr->SetStatus(AliESDtrack::kEmbedded);
828
829 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
830 dummy = (*tracks)[nAODtracks-1];
831
832 if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
b725dccf 833 fh1TrackPt->Fill(tmpTr->Pt());
834 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
31b9d515 835 }
836 }
837
838 if(mcpartIN){
839 Int_t nMCpart = mcpartIN->GetEntriesFast();
840
841 Int_t nPhysicalPrimary=0;
842 Int_t nAODmcpart=0;
843 for(Int_t ip=0; ip<nMCpart; ++ip){
844 AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
845
846 if(fEmbedMode==kAODJetTracks){
847 // jet track? else continue
848 // not implemented yet
849 continue;
850 }
851
852 new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
853 dummy = (*mcpartOUT)[nAODmcpart-1];
854
855 if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
856 fh1MCTrackPt->Fill(tmpPart->Pt());
857 fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
858 nPhysicalPrimary++;
859 }
860 }
861 fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
862
863 }
864 } // end: embed all tracks or jet tracks
865
866
867 if(fEmbedMode==kAODJet4Mom){
868
869 // loop over jets
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));
874 if(!jet) continue;
875 AliAODTrack *tmpTr = (AliAODTrack*)jet;
876 tmpTr->SetFlags(AliESDtrack::kEmbedded);
877
878 new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
879 dummy = (*tracks)[nAODtracks-1];
880
881 fh1TrackPt->Fill(tmpTr->Pt());
882 fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
883 }
884
885 } // end: embed jets as 4-momenta
886
887
888 } //end: embed mode with AOD
889
890
891 // === embed mode with toy ===
892 if(fEmbedMode==kToyTracks){
893 Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
894
895 fh1TrackN->Fill((Float_t)nT);
896
897 for(Int_t i=0; i<nT; ++i){
898
899 Double_t pt = -1.;
900 if(fToyMinTrackPt!=fToyMaxTrackPt){
901 if(fToyDistributionTrackPt==0){
902 pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
903 } else {
904 while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
905 pt = rndm->Exp(fToyDistributionTrackPt); // no power law yet!!
906 pt += fToyMinTrackPt;
907 }
908 }
909 } else {
910 pt = fToyMinTrackPt;
911 }
912 Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
913 Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
914 phi = TVector2::Phi_0_2pi(phi);
915
916 if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
917
918 Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
919 Float_t mom[3];
920 mom[0] = pt;
921 mom[1] = phi;
922 mom[2] = theta;
923
924 Float_t xyz[3];
925 xyz[0] = -999.;
926 xyz[1] = -999.;
927 xyz[2] = -999.;
928
929 AliAODTrack *tmpTr = new AliAODTrack( -999, // id
930 -999, // label
931 mom, // momentum[3]
932 kFALSE, // cartesian
933 xyz, // position
934 kFALSE, // DCA
935 NULL, // covMatrix,
936 -99, // charge
937 0, // itsClusMap
938 NULL, // pid
939 NULL, // prodVertex
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
945 );
946 tmpTr->SetFlags(AliESDtrack::kEmbedded);
947
948 new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
949 dummy = (*tracks)[nAODtracks-1];
950
951 fh1TrackPt->Fill(pt);
952 fh2TrackEtaPhi->Fill(eta,phi);
953
954 delete tmpTr;
955 }
956 } //end: embed mode with toy
957
958
959 PostData(1, fHistList);
b725dccf 960}
961
962
963//__________________________________________________________________________
964void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
965{
7a88ef30 966 // terminate
31b9d515 967 if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
b725dccf 968
31b9d515 969 if(fAODfile && fAODfile->IsOpen())
970 fAODfile->Close();
b725dccf 971
972}
973
974//__________________________________________________________________________
b725dccf 975Int_t AliAnalysisTaskFastEmbedding::GetJobID()
976{
7a88ef30 977 // gives back the alien subjob id, if available, else -1
978
0ed51e0d 979 Int_t id=-1;
b725dccf 980
0ed51e0d 981 const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
982 //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
b725dccf 983
984 if(env && strlen(env)){
31b9d515 985 id= atoi(env);
986 AliInfo(Form("Job index %d", id));
b725dccf 987 }
988 else{
31b9d515 989 AliInfo("Job index not found. Okay if running locally.");
b725dccf 990 }
991
992 return id;
0ed51e0d 993}
b725dccf 994
995//__________________________________________________________________________
7a88ef30 996Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
997{
998 // select an AOD file from fAODPathArray
b725dccf 999
31b9d515 1000 Int_t nFiles = fAODPathArray->GetEntries();
1001
1002 // choose randomly file and event
1003 Int_t n = -1;
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);
7a88ef30 1007 if(fAODEntriesArray){
31b9d515 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.;
1013 }
1014 fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
7a88ef30 1015 } else {
1016 AliWarning("Number of entries in extra AODs not set!");
1017 n = rndm->Integer(nFiles);
1018 fAODEntry = 0;
1019 }
31b9d515 1020
1021 AliInfo(Form("Select AOD file %d", n));
1022 TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1023 if(!objStr){
1024 AliError("Could not get path of aod file.");
1025 return -1;
1026 }
1027 fAODPath = objStr->GetString();
1028
1029 return n;
1a2bb6d5 1030}
1031
46465e39 1032//__________________________________________________________________________
7a88ef30 1033Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1034{
1035 // select and open an AOD file from fAODPathArray
46465e39 1036
31b9d515 1037 if(fAODPathArray) fFileId = SelectAODfile();
1038 if(fFileId<0) return -1;
1039
1040 TDirectory *owd = gDirectory;
7a88ef30 1041 if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
31b9d515 1042 fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1043 owd->cd();
1044 if(!fAODfile){
1045
1046 fFileId = -1;
1047 if(fAODPathArray){
1048 if(trial<=3){
1049 AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1050 fFileId = OpenAODfile(trial+1);
1051 } else {
1052 AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1053 return -1;
1054 }
1055 } else {
1056 AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1057 }
1058
1059 return fFileId;
1060 }
1061
1062 fAODtree = (TTree*)fAODfile->Get("aodTree");
1063
1064 if(!fAODtree){
1065 AliError("AOD tree not found.");
1066 return -1;
1067 }
1068
1069 /*
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);
1077 */
1078
1079 delete fAODevent;
1080 fAODevent = new AliAODEvent();
1081 fAODevent->ReadFromTree(fAODtree);
1082
1083
1084 // fetch header, jets, etc. from new file
1085 fNevents = fAODtree->GetEntries();
1086 mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1087
1088 if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1089 else fAODJets = fAODevent->GetJets();
1090 if(!fAODJets){
1091 AliError("Could not find jets in AOD. Check jet branch when indicated.");
1092 return -1;
1093 }
1094
1095 TFile *curfile = fAODtree->GetCurrentFile();
1096 if (!curfile) {
1097 AliError("No current file.");
1098 return -1;
1099 }
1100
1101 Float_t trials = 1.;
1102 Float_t xsec = 0.;
1103 PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1104 fXsection = xsec;
1105
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;
1109
1110 if(fFileId>=0){
1111 AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1112 }
1113
1114 fCountEvents=0; // new file, reset counter
1115
1116 return fFileId; // file position in AOD path array, if array available
1117}
1118
1119
1120//____________________________________________________________________________
1121Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1122
7a88ef30 1123 // static stored, available for other tasks in train
1124
31b9d515 1125 static Float_t ptHard = -1.;
1126 if(bSet) ptHard = newValue;
1127
1128 return ptHard;
1a2bb6d5 1129}
b725dccf 1130
31b9d515 1131//____________________________________________________________________________
1132Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1133
7a88ef30 1134 // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1135
31b9d515 1136 const Int_t nBins = 10;
1137 Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1138
1139 Int_t bin = -1;
1140 while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1141 bin++;
1142 }
1143
1144 return bin;
1145}
1146
1147//____________________________________________________________________________
1148Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1149 //
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)
1153
1154 TString file(currFile);
1155 fXsec = 0;
1156 fTrials = 1;
1157
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,"");
1162 }
1163 else {
1164 // not an archive take the basename....
1165 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1166 }
1167
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
1169 if(!fxsec){
1170 // next trial fetch the histgram file
1171 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1172 if(!fxsec){
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()));
1175 return kFALSE;
1176 }
1177 else{
1178 // find the tlist we want to be independtent of the name so use the Tkey
1179 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1180 if(!key){
1181 fxsec->Close();
1182 return kFALSE;
1183 }
1184 TList *list = dynamic_cast<TList*>(key->ReadObj());
1185 if(!list){
1186 fxsec->Close();
1187 return kFALSE;
1188 }
1189 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1190 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1191 fxsec->Close();
1192 }
1193 } // no tree pyxsec.root
1194 else {
1195 TTree *xtree = (TTree*)fxsec->Get("Xsection");
1196 if(!xtree){
1197 fxsec->Close();
1198 return kFALSE;
1199 }
1200 UInt_t ntrials = 0;
1201 Double_t xsection = 0;
1202 xtree->SetBranchAddress("xsection",&xsection);
1203 xtree->SetBranchAddress("ntrials",&ntrials);
1204 xtree->GetEntry(0);
1205 fTrials = ntrials;
1206 fXsec = xsection;
1207 fxsec->Close();
1208 }
1209 return kTRUE;
1210}
1211
1212