Coverity fix
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskFastEmbedding.cxx
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>
36 #include <TProfile.h>
37 #include <TKey.h>
38
39
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"
50
51 #include "AliLog.h"
52
53 ClassImp(AliAnalysisTaskFastEmbedding)
54
55 //__________________________________________________________________________
56 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
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)
133 {
134    // default constructor
135
136 }
137
138
139 //__________________________________________________________________________
140 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
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)
217 {
218    // constructor
219    DefineOutput(1, TList::Class());
220 }
221
222
223 //__________________________________________________________________________
224 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy)
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)
301 {
302    // copy constructor
303 }
304
305
306 //__________________________________________________________________________
307 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
308 {
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;
391 }
392
393
394 //__________________________________________________________________________
395 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
396 {
397    // destructor
398    delete rndm;
399 }
400
401
402 //__________________________________________________________________________
403 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
404 {
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);
537 }
538
539
540 //__________________________________________________________________________
541 void AliAnalysisTaskFastEmbedding::Init()
542 {
543    // Initialization
544    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
545
546 }
547
548 //__________________________________________________________________________
549 Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
550 {
551
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;
560
561 }
562
563
564 //__________________________________________________________________________
565 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
566 {
567    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
568
569    if(!fAODout){
570       AliError("Need output AOD, but is not connected."); 
571       PostData(1, fHistList);
572       return;
573    }
574
575    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
576    if (!fESD) {
577       AliError("ESD not available");
578       PostData(1, fHistList);
579       return;
580    }
581
582    // -- event selection --
583    fHistEvtSelection->Fill(1); // number of events before event selection
584
585    // physics selection
586    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
587    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
588    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
589       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
590       fHistEvtSelection->Fill(2);
591       PostData(1, fHistList);
592       return;
593    } 
594
595    // vertex selection
596    AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
597    Int_t nTracksPrim = primVtx->GetNContributors();
598    if ((nTracksPrim < fMinContribVtx) ||
599          (primVtx->GetZ() < fVtxZMin) ||
600          (primVtx->GetZ() > fVtxZMax) ){
601       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
602       fHistEvtSelection->Fill(3);
603       PostData(1, fHistList);
604       return;
605    }
606
607    /** takes wrong value, since jet service tasks runs later
608    // event class selection (from jet helper task)
609    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
610    if(fDebug) Printf("Event class %d", eventClass);
611    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
612    fHistEvtSelection->Fill(4);
613    PostData(1, fHistList);
614    return;
615    }*/
616
617    // centrality selection
618    AliCentrality *cent = 0x0;
619    Float_t centValue = 0.; 
620    if(fESD) cent = fESD->GetCentrality();
621    if(cent) centValue = cent->GetCentralityPercentile("V0M");
622    if(fDebug) printf("centrality: %f\n", centValue);
623    if (centValue < fCentMin || centValue > fCentMax){
624       fHistEvtSelection->Fill(4);
625       PostData(1, fHistList);
626       return;
627    }
628
629
630    /*  ** not implemented **
631    // multiplicity due to input tracks
632    Int_t nInputTracks = GetNInputTracks();
633
634    if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
635    fHistEvtSelection->Fill(5);
636    PostData(1, fHistList);
637    return;
638    }
639    */
640
641    fHistEvtSelection->Fill(0); // accepted events  
642    // -- end event selection --
643
644    // connect aod out
645    TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
646    if(!tracks){
647       AliError("Extra track branch not found in output.");
648       PostData(1, fHistList);
649       return;
650    }
651    tracks->Delete();
652    Int_t nAODtracks=0;
653
654    TRef dummy;
655
656    // === embed mode with AOD ===
657    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
658       if(!fAODevent){
659          AliError("Need input AOD, but is not connected."); 
660          PostData(1, fHistList);
661          return;
662       }
663
664
665
666       Bool_t useEntry = kFALSE;
667       while(!useEntry){  // protection need, if no event fulfills requierment
668
669          fAODEntry++; // go to next event 
670          fCountEvents++;
671          if(fAODEntry>=fNevents) fAODEntry=0; 
672
673          // open new aod file
674          if(fCountEvents>=fNevents){ 
675             if(!fAODPathArray){
676                AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
677             } 
678             else {
679                AliDebug(AliLog::kDebug, "Select new AOD file ...");
680
681                fFileId = OpenAODfile();
682                if(fFileId<0) {
683                   PostData(1, fHistList);
684                   return;
685                }
686                fh1AODfile->Fill(fFileId);
687                if(fAODEntries!=fNevents){
688                   AliError("File with list of AODs and file with nb. of entries inconsistent");
689                   PostData(1, fHistList);
690                   return;
691                }
692             }
693          }
694
695
696
697          // get next event
698          fAODtree->GetEvent(fAODEntry);
699
700          // get pt hard
701          if(mcHeader){
702             fPtHard = mcHeader->GetPtHard();
703             GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
704             AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
705          }
706          else{
707             AliDebug(AliLog::kDebug,"No mcHeader");
708          }
709          fPtHardBin = GetPtHardBin(fPtHard);
710
711          //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials); 
712
713          // fill event variables for each event
714          fh1Xsec->Fill(fPtHardBin,fXsection);
715          fh2PtHard->Fill(fPtHardBin,fPtHard);
716
717          fh1Trials->Fill(fPtHardBin, fAvgTrials);
718          fh1TrialsEvtSel->Fill(fPtHardBin);
719          fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
720
721          // jet pt selection
722          if(fEvtSelecMode==kEventsJetPt){
723             Int_t nJets = fAODJets->GetEntries();
724             //Printf("nb. jets: %d",nJets);
725             for(Int_t ij=0; ij<nJets; ++ij){
726                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
727                if(!jet) continue;
728
729                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
730                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
731                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
732                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
733                   useEntry = kTRUE;
734                   break;
735                } 
736             }
737          }
738
739          // no selection
740          if(fEvtSelecMode==kEventsAll){
741             useEntry = kTRUE;
742          }
743
744       }
745       AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
746
747       fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
748       fh2AODevent->Fill(fFileId,fAODEntry);
749
750       TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
751       TClonesArray *mcpartOUT = 0x0;
752       if(mcpartIN){
753          mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
754          if(!mcpartOUT){
755             AliError("Extra MC particles branch not found in output.");
756             PostData(1, fHistList);
757             return;
758          }
759          mcpartOUT->Delete();
760       } else {
761          AliInfo("No extra MC particles found.");
762       }
763
764
765       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
766          // loop over input tracks
767          // add to output aod
768          Int_t nTracks = 0;
769          Int_t nJets = fAODJets->GetEntries();
770          Int_t nSelectedJets = 0;
771          AliAODJet *leadJet = 0x0; // used for jet tracks only
772
773          if(fEmbedMode==kAODFull){
774             nTracks = fAODevent->GetNumberOfTracks();
775
776             for(Int_t ij=0; ij<nJets; ++ij){
777                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
778                if(!jet) continue;
779                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
780                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
781                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
782                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
783
784                   fh1JetPt->Fill(jet->Pt());
785                   fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
786                   nSelectedJets++;
787
788                }
789             }                              
790             fh1JetN->Fill(nSelectedJets);
791          }
792
793          if(fEmbedMode==kAODJetTracks){
794             // look for leading jet within selection
795             // get reference tracks for this jet
796             Float_t leadJetPt = 0.;
797             for(Int_t ij=0; ij<nJets; ++ij){
798                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
799                if(!jet) continue;
800                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
801                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
802                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
803                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
804                   if(jet->Pt()>leadJetPt){
805                      leadJetPt = jet->Pt();
806                      leadJet = jet;
807                   } 
808                }
809             }
810             if(leadJet){
811                nTracks = leadJet->GetRefTracks()->GetEntriesFast();
812                fh1JetPt->Fill(leadJet->Pt());
813                fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
814                fh1JetN->Fill(1);                                   
815             }
816          }
817
818          fh1TrackN->Fill((Float_t)nTracks);
819
820          for(Int_t it=0; it<nTracks; ++it){
821             AliAODTrack *tmpTr = 0x0;
822             if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
823             if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
824             if(!tmpTr) continue;
825
826             tmpTr->SetStatus(AliESDtrack::kEmbedded);
827
828             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
829             dummy = (*tracks)[nAODtracks-1];
830
831             if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
832                fh1TrackPt->Fill(tmpTr->Pt());
833                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
834             }
835          }
836
837          if(mcpartIN){
838             Int_t nMCpart = mcpartIN->GetEntriesFast();
839
840             Int_t nPhysicalPrimary=0;
841             Int_t nAODmcpart=0;
842             for(Int_t ip=0; ip<nMCpart; ++ip){
843                AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
844
845                if(fEmbedMode==kAODJetTracks){
846                   // jet track? else continue
847                   // not implemented yet
848                   continue;
849                } 
850
851                new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
852                dummy = (*mcpartOUT)[nAODmcpart-1];
853
854                if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
855                   fh1MCTrackPt->Fill(tmpPart->Pt());
856                   fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
857                   nPhysicalPrimary++;
858                }
859             }
860             fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
861
862          }
863       } // end: embed all tracks or jet tracks
864
865
866       if(fEmbedMode==kAODJet4Mom){
867
868          // loop over jets
869          Int_t nJets = fAODJets->GetEntries();
870          fh1TrackN->Fill((Float_t)nJets);
871          for(Int_t ij=0; ij<nJets; ++ij){
872             AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
873             if(!jet) continue;
874             AliAODTrack *tmpTr = (AliAODTrack*)jet;
875             tmpTr->SetFlags(AliESDtrack::kEmbedded);
876
877             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
878             dummy = (*tracks)[nAODtracks-1]; 
879
880             fh1TrackPt->Fill(tmpTr->Pt());
881             fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
882          }
883
884       } // end: embed jets as 4-momenta
885
886
887    } //end: embed mode with AOD
888
889
890    // === embed mode with toy ===
891    if(fEmbedMode==kToyTracks){
892       Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
893
894       fh1TrackN->Fill((Float_t)nT);
895
896       for(Int_t i=0; i<nT; ++i){
897
898          Double_t pt = -1.;
899          if(fToyMinTrackPt!=fToyMaxTrackPt){
900             if(fToyDistributionTrackPt==0){
901                pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
902             } else {
903                while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
904                   pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
905                   pt += fToyMinTrackPt;
906                }
907             }
908          } else {
909             pt = fToyMinTrackPt;
910          }
911          Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
912          Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
913          phi = TVector2::Phi_0_2pi(phi);
914
915          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
916
917          Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
918          Float_t mom[3];
919          mom[0] = pt;
920          mom[1] = phi;
921          mom[2] = theta;
922
923          Float_t xyz[3];
924          xyz[0] = -999.;
925          xyz[1] = -999.;
926          xyz[2] = -999.;
927
928          AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
929          -999,   // label
930          mom,    // momentum[3]
931          kFALSE, // cartesian
932          xyz,    // position
933          kFALSE, // DCA
934          NULL,   // covMatrix,
935          -99,    // charge
936          0,      // itsClusMap
937          NULL,   // pid 
938          NULL,   // prodVertex
939          kFALSE, // used for vertex fit
940          kFALSE, // used for prim vtx fit
941          AliAODTrack::kUndef, // type
942          fToyFilterMap,  // select info
943          -999.    // chi2 per NDF
944          );
945          tmpTr->SetFlags(AliESDtrack::kEmbedded);
946
947          new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
948          dummy = (*tracks)[nAODtracks-1];
949
950          fh1TrackPt->Fill(pt);
951          fh2TrackEtaPhi->Fill(eta,phi);
952
953          delete tmpTr;
954       }
955    } //end: embed mode with toy
956
957
958    PostData(1, fHistList);
959 }
960
961
962 //__________________________________________________________________________
963 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
964 {
965    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
966
967    if(fAODfile && fAODfile->IsOpen())
968    fAODfile->Close();  
969
970 }
971
972 //__________________________________________________________________________
973 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
974 {
975    Int_t id=-1;
976
977    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
978    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
979
980    if(env && strlen(env)){
981       id= atoi(env);
982       AliInfo(Form("Job index %d", id));
983    }
984    else{
985       AliInfo("Job index not found. Okay if running locally.");
986    }
987
988    return id;
989 }
990
991 //__________________________________________________________________________
992
993 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
994
995    Int_t nFiles = fAODPathArray->GetEntries();
996
997    // choose randomly file and event
998    Int_t n = -1;
999    Float_t tmpProp = -1.;
1000    Int_t pendingEvents = fInputEntries-Entry();
1001    //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1002    while(rndm->Rndm()>tmpProp){
1003       n = rndm->Integer(nFiles);
1004       fAODEntries = fAODEntriesArray->At(n);
1005       Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1006       tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1007    }
1008    fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1009
1010    AliInfo(Form("Select AOD file %d", n));
1011    TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1012    if(!objStr){
1013       AliError("Could not get path of aod file.");
1014       return -1;
1015    }
1016    fAODPath = objStr->GetString();
1017
1018    return n;
1019 }
1020
1021 //__________________________________________________________________________
1022
1023 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
1024
1025    if(fAODPathArray) fFileId = SelectAODfile();
1026    if(fFileId<0) return -1;
1027
1028    TDirectory *owd = gDirectory;
1029    if (fAODfile) fAODfile->Close();
1030    fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1031    owd->cd();
1032    if(!fAODfile){
1033
1034       fFileId = -1;
1035       if(fAODPathArray){
1036          if(trial<=3){ 
1037             AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1038             fFileId = OpenAODfile(trial+1);
1039          } else {
1040             AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1041             return -1;
1042          }
1043       } else {
1044          AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1045       }
1046
1047       return fFileId;
1048    }
1049
1050    fAODtree = (TTree*)fAODfile->Get("aodTree");
1051
1052    if(!fAODtree){
1053       AliError("AOD tree not found.");
1054       return -1;
1055    }
1056
1057    /*
1058       fAODtree->SetBranchStatus("*",0);
1059       fAODtree->SetBranchStatus("header",1);
1060       fAODtree->SetBranchStatus("tracks*",1);
1061       fAODtree->SetBranchStatus("jets*",1);
1062       fAODtree->SetBranchStatus("clusters*",1);
1063       fAODtree->SetBranchStatus("mcparticles*",1);
1064       fAODtree->SetBranchStatus("mcHeader*",1);
1065    */
1066
1067    delete fAODevent;
1068    fAODevent = new AliAODEvent();
1069    fAODevent->ReadFromTree(fAODtree);
1070
1071
1072    // fetch header, jets, etc. from new file
1073    fNevents = fAODtree->GetEntries();
1074    mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1075
1076    if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1077    else                    fAODJets = fAODevent->GetJets();
1078    if(!fAODJets){
1079       AliError("Could not find jets in AOD. Check jet branch when indicated.");
1080       return -1;
1081    }
1082
1083    TFile *curfile = fAODtree->GetCurrentFile();
1084    if (!curfile) {
1085       AliError("No current file.");
1086       return -1;
1087    }
1088
1089    Float_t trials = 1.;
1090    Float_t xsec = 0.;
1091    PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1092    fXsection = xsec;
1093
1094    // construct a poor man average trials 
1095    Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1096    if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1097
1098    if(fFileId>=0){
1099       AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1100    }
1101
1102    fCountEvents=0; // new file, reset counter
1103
1104    return fFileId;  // file position in AOD path array, if array available
1105 }
1106
1107
1108 //____________________________________________________________________________
1109 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1110
1111    static Float_t ptHard = -1.;
1112    if(bSet) ptHard = newValue;
1113
1114    return ptHard;
1115 }
1116
1117 //____________________________________________________________________________
1118 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1119
1120    const Int_t nBins = 10;
1121    Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1122
1123    Int_t bin = -1;
1124    while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1125       bin++;
1126    }
1127
1128    return bin;
1129 }
1130
1131 //____________________________________________________________________________
1132 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1133    //
1134    // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1135    // This should provide the path to the AOD/ESD file
1136    // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks) 
1137
1138    TString file(currFile);  
1139    fXsec = 0;
1140    fTrials = 1;
1141
1142    if(file.Contains("root_archive.zip#")){
1143       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1144       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1145       file.Replace(pos+1,20,"");
1146    }
1147    else {
1148       // not an archive take the basename....
1149       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1150    }
1151
1152    TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1153    if(!fxsec){
1154       // next trial fetch the histgram file
1155       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1156       if(!fxsec){
1157          // not a severe condition but inciate that we have no information
1158          AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1159          return kFALSE;
1160       }
1161       else{
1162          // find the tlist we want to be independtent of the name so use the Tkey
1163          TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1164          if(!key){
1165             fxsec->Close();
1166             return kFALSE;
1167          }
1168          TList *list = dynamic_cast<TList*>(key->ReadObj());
1169          if(!list){
1170             fxsec->Close();
1171             return kFALSE;
1172          }
1173          fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1174          fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1175          fxsec->Close();
1176       }
1177    } // no tree pyxsec.root
1178    else {
1179       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1180       if(!xtree){
1181          fxsec->Close();
1182          return kFALSE;
1183       }
1184       UInt_t   ntrials  = 0;
1185       Double_t  xsection  = 0;
1186       xtree->SetBranchAddress("xsection",&xsection);
1187       xtree->SetBranchAddress("ntrials",&ntrials);
1188       xtree->GetEntry(0);
1189       fTrials = ntrials;
1190       fXsec = xsection;
1191       fxsec->Close();
1192    }
1193    return kTRUE;
1194 }
1195
1196