Protections for coverity: DIVIDE_BY_ZERO
[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    // User defined Notify(), called once
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    // Analysis, called once per event
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."); 
661          PostData(1, fHistList);
662          return;
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)
781                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
782                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
783                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
784
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));
800                if(!jet) continue;
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          }
818
819          fh1TrackN->Fill((Float_t)nTracks);
820
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)){
833                fh1TrackPt->Fill(tmpTr->Pt());
834                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
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);
960 }
961
962
963 //__________________________________________________________________________
964 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
965 {
966    // terminate
967    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
968
969    if(fAODfile && fAODfile->IsOpen())
970    fAODfile->Close();  
971
972 }
973
974 //__________________________________________________________________________
975 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
976 {
977    // gives back the alien subjob id, if available, else -1
978
979    Int_t id=-1;
980
981    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
982    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
983
984    if(env && strlen(env)){
985       id= atoi(env);
986       AliInfo(Form("Job index %d", id));
987    }
988    else{
989       AliInfo("Job index not found. Okay if running locally.");
990    }
991
992    return id;
993 }
994
995 //__________________________________________________________________________
996 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
997 {
998    // select an AOD file from fAODPathArray
999
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);
1007    if(fAODEntriesArray){
1008    while(rndm->Rndm()>tmpProp){
1009       n = rndm->Integer(nFiles);
1010       fAODEntries = fAODEntriesArray->At(n);
1011       Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1012       tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1013    }
1014    fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1015    } else {
1016       AliWarning("Number of entries in extra AODs not set!");
1017       n = rndm->Integer(nFiles);
1018       fAODEntry = 0;
1019    }
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;
1030 }
1031
1032 //__________________________________________________________________________
1033 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1034 {
1035    // select and open an AOD file from fAODPathArray
1036
1037    if(fAODPathArray) fFileId = SelectAODfile();
1038    if(fFileId<0) return -1;
1039
1040    TDirectory *owd = gDirectory;
1041    if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
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 //____________________________________________________________________________
1121 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1122
1123    // static stored, available for other tasks in train
1124
1125    static Float_t ptHard = -1.;
1126    if(bSet) ptHard = newValue;
1127
1128    return ptHard;
1129 }
1130
1131 //____________________________________________________________________________
1132 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1133
1134    // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1135
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 //____________________________________________________________________________
1148 Bool_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