589d0b233efd5308a526327c86cf8c92a59f0951
[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 #include <TGrid.h>
39
40
41 #include "AliAnalysisTaskFastEmbedding.h"
42 #include "AliAnalysisManager.h"
43 #include "AliESDEvent.h"
44 #include "AliESDtrack.h"
45 #include "AliAODEvent.h"
46 #include "AliAODTrack.h"
47 #include "AliAODJet.h"
48 #include "AliAODMCParticle.h"
49 #include "AliAODMCHeader.h"
50 #include "AliInputEventHandler.h"
51
52 #include "AliLog.h"
53
54 ClassImp(AliAnalysisTaskFastEmbedding)
55
56 //__________________________________________________________________________
57 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
58 : AliAnalysisTaskSE()
59 ,fESD(0)
60 ,fAODout(0)
61 ,fAODevent(0)
62 ,fAODtree(0)
63 ,fAODfile(0)
64 ,mcHeader(0)
65 ,rndm(0)
66 ,fInputEntries(0)
67 ,fAODPathArray(0)
68 ,fAODEntriesArray(0)
69 ,fAODPath("AliAOD.root")
70 ,fAODEntries(-1)
71 ,fAODEntriesSum(0)
72 ,fAODEntriesMax(0)
73 ,fOfflineTrgMask(AliVEvent::kAny)
74 ,fMinContribVtx(1)
75 ,fVtxZMin(-8.)
76 ,fVtxZMax(8.)
77 ,fEvtClassMin(0)
78 ,fEvtClassMax(4)
79 ,fCentMin(0.)
80 ,fCentMax(100.)
81 ,fNInputTracksMin(0)
82 ,fNInputTracksMax(-1)
83 ,fTrackBranch("aodExtraTracks")
84 ,fMCparticlesBranch("aodExtraMCparticles")
85 ,fJetBranch("")
86 ,fFileId(-1)
87 ,fAODEntry(0)
88 ,fCountEvents(-1)
89 ,fEmbedMode(0)
90 ,fEvtSelecMode(0)
91 ,fEvtSelMinJetPt(-1)
92 ,fEvtSelMaxJetPt(-1)
93 ,fEvtSelMinJetEta(-999.)
94 ,fEvtSelMaxJetEta( 999.)
95 ,fEvtSelMinJetPhi(0.)
96 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
97 ,fToyMinNbOfTracks(1)
98 ,fToyMaxNbOfTracks(1)
99 ,fToyMinTrackPt(50.)
100 ,fToyMaxTrackPt(50.)
101 ,fToyDistributionTrackPt(0.)
102 ,fToyMinTrackEta(-.5)
103 ,fToyMaxTrackEta(.5)
104 ,fToyMinTrackPhi(0.)
105 ,fToyMaxTrackPhi(2*TMath::Pi())
106 ,fToyFilterMap(0)
107 ,fTrackFilterMap(0)
108 ,fNPtHard(10)
109 ,fPtHard(0)
110 ,fPtHardBin(-1)
111 ,fAODJets(0x0)
112 ,fNevents(0)
113 ,fXsection(0)
114 ,fAvgTrials(0)
115 ,fHistList(0)
116 ,fHistEvtSelection(0)
117 ,fh1Xsec(0)
118 ,fh1Trials(0)
119 ,fh1TrialsEvtSel(0)
120 ,fh2PtHard(0)
121 ,fh2PtHardEvtSel(0)
122 ,fh2PtHardTrials(0)
123 ,fh1TrackPt(0)
124 ,fh2TrackEtaPhi(0)
125 ,fh1TrackN(0)
126 ,fh1JetPt(0)
127 ,fh2JetEtaPhi(0)
128 ,fh1JetN(0)
129 ,fh1MCTrackPt(0)
130 ,fh2MCTrackEtaPhi(0)
131 ,fh1MCTrackN(0)
132 ,fh1AODfile(0)
133 ,fh2AODevent(0)
134 {
135    // default constructor
136
137 }
138
139
140 //__________________________________________________________________________
141 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
142 : AliAnalysisTaskSE(name)
143 ,fESD(0)
144 ,fAODout(0)
145 ,fAODevent(0)
146 ,fAODtree(0)
147 ,fAODfile(0)
148 ,mcHeader(0)
149 ,rndm(0)
150 ,fInputEntries(0)
151 ,fAODPathArray(0)
152 ,fAODEntriesArray(0)
153 ,fAODPath("AliAOD.root")
154 ,fAODEntries(-1)
155 ,fAODEntriesSum(0)
156 ,fAODEntriesMax(0)
157 ,fOfflineTrgMask(AliVEvent::kAny)
158 ,fMinContribVtx(1)
159 ,fVtxZMin(-8.)
160 ,fVtxZMax(8.)
161 ,fEvtClassMin(0)
162 ,fEvtClassMax(4)
163 ,fCentMin(0.)
164 ,fCentMax(100.)
165 ,fNInputTracksMin(0)
166 ,fNInputTracksMax(-1)
167 ,fTrackBranch("aodExtraTracks")
168 ,fMCparticlesBranch("aodExtraMCparticles")
169 ,fJetBranch("")
170 ,fFileId(-1)
171 ,fAODEntry(0)
172 ,fCountEvents(-1)
173 ,fEmbedMode(0)
174 ,fEvtSelecMode(0)
175 ,fEvtSelMinJetPt(-1)
176 ,fEvtSelMaxJetPt(-1)
177 ,fEvtSelMinJetEta(-999.)
178 ,fEvtSelMaxJetEta( 999.)
179 ,fEvtSelMinJetPhi(0.)
180 ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
181 ,fToyMinNbOfTracks(1)
182 ,fToyMaxNbOfTracks(1)
183 ,fToyMinTrackPt(50.)
184 ,fToyMaxTrackPt(50.)
185 ,fToyDistributionTrackPt(0.)
186 ,fToyMinTrackEta(-.5)
187 ,fToyMaxTrackEta(.5)
188 ,fToyMinTrackPhi(0.)
189 ,fToyMaxTrackPhi(2*TMath::Pi())
190 ,fToyFilterMap(0)
191 ,fTrackFilterMap(0)
192 ,fNPtHard(10)
193 ,fPtHard(0)
194 ,fPtHardBin(-1)
195 ,fAODJets(0x0)
196 ,fNevents(0)
197 ,fXsection(0)
198 ,fAvgTrials(0)
199 ,fHistList(0)
200 ,fHistEvtSelection(0)
201 ,fh1Xsec(0)
202 ,fh1Trials(0)
203 ,fh1TrialsEvtSel(0)
204 ,fh2PtHard(0)
205 ,fh2PtHardEvtSel(0)
206 ,fh2PtHardTrials(0)
207 ,fh1TrackPt(0)
208 ,fh2TrackEtaPhi(0)
209 ,fh1TrackN(0)
210 ,fh1JetPt(0)
211 ,fh2JetEtaPhi(0)
212 ,fh1JetN(0)
213 ,fh1MCTrackPt(0)
214 ,fh2MCTrackEtaPhi(0)
215 ,fh1MCTrackN(0)
216 ,fh1AODfile(0)
217 ,fh2AODevent(0)
218 {
219    // constructor
220    DefineOutput(1, TList::Class());
221 }
222
223
224 //__________________________________________________________________________
225 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy)
226 : AliAnalysisTaskSE()
227 ,fESD(copy.fESD)
228 ,fAODout(copy.fAODout)
229 ,fAODevent(copy.fAODevent)
230 ,fAODtree(copy.fAODtree)
231 ,fAODfile(copy.fAODfile)
232 ,mcHeader(copy.mcHeader)
233 ,rndm(copy.rndm)
234 ,fInputEntries(copy.fInputEntries)
235 ,fAODPathArray(copy.fAODPathArray)
236 ,fAODEntriesArray(copy.fAODEntriesArray)
237 ,fAODPath(copy.fAODPath)
238 ,fAODEntries(copy.fAODEntries)
239 ,fAODEntriesSum(copy.fAODEntriesSum)
240 ,fAODEntriesMax(copy.fAODEntriesMax)
241 ,fOfflineTrgMask(copy.fOfflineTrgMask)
242 ,fMinContribVtx(copy.fMinContribVtx)
243 ,fVtxZMin(copy.fVtxZMin)
244 ,fVtxZMax(copy.fVtxZMax)
245 ,fEvtClassMin(copy.fEvtClassMin)
246 ,fEvtClassMax(copy.fEvtClassMax)
247 ,fCentMin(copy.fCentMin)
248 ,fCentMax(copy.fCentMax)
249 ,fNInputTracksMin(copy.fNInputTracksMin)
250 ,fNInputTracksMax(copy.fNInputTracksMax)
251 ,fTrackBranch(copy.fTrackBranch)
252 ,fMCparticlesBranch(copy.fMCparticlesBranch)
253 ,fJetBranch(copy.fJetBranch)
254 ,fFileId(copy.fFileId)
255 ,fAODEntry(copy.fAODEntry)
256 ,fCountEvents(copy.fCountEvents)
257 ,fEmbedMode(copy.fEmbedMode)
258 ,fEvtSelecMode(copy.fEvtSelecMode)
259 ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
260 ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
261 ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
262 ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
263 ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
264 ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
265 ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
266 ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
267 ,fToyMinTrackPt(copy.fToyMinTrackPt)
268 ,fToyMaxTrackPt(copy.fToyMaxTrackPt)
269 ,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
270 ,fToyMinTrackEta(copy.fToyMinTrackEta)
271 ,fToyMaxTrackEta(copy.fToyMaxTrackEta)
272 ,fToyMinTrackPhi(copy.fToyMinTrackPhi)
273 ,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
274 ,fToyFilterMap(copy.fToyFilterMap)
275 ,fTrackFilterMap(copy.fTrackFilterMap)
276 ,fNPtHard(copy.fNPtHard)
277 ,fPtHard(copy.fPtHard)
278 ,fPtHardBin(copy.fPtHardBin)
279 ,fAODJets(copy.fAODJets)
280 ,fNevents(copy.fNevents)
281 ,fXsection(copy.fXsection)
282 ,fAvgTrials(copy.fAvgTrials)
283 ,fHistList(copy.fHistList)
284 ,fHistEvtSelection(copy.fHistEvtSelection)
285 ,fh1Xsec(copy.fh1Xsec)
286 ,fh1Trials(copy.fh1Trials)
287 ,fh1TrialsEvtSel(copy.fh1TrialsEvtSel)
288 ,fh2PtHard(copy.fh2PtHard)
289 ,fh2PtHardEvtSel(copy.fh2PtHardEvtSel)
290 ,fh2PtHardTrials(copy.fh2PtHardTrials)
291 ,fh1TrackPt(copy.fh1TrackPt)
292 ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
293 ,fh1TrackN(copy.fh1TrackN)
294 ,fh1JetPt(copy.fh1JetPt)
295 ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
296 ,fh1JetN(copy.fh1JetN)
297 ,fh1MCTrackPt(copy.fh1MCTrackPt)
298 ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
299 ,fh1MCTrackN(copy.fh1MCTrackN)
300 ,fh1AODfile(copy.fh1AODfile)
301 ,fh2AODevent(copy.fh2AODevent)
302 {
303    // copy constructor
304 }
305
306
307 //__________________________________________________________________________
308 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
309 {
310    // assignment
311
312    if(this!=&o){
313       AliAnalysisTaskSE::operator=(o);
314       fESD               = o.fESD;
315       fAODout            = o.fAODout;
316       fAODevent          = o.fAODevent;
317       fAODtree           = o.fAODtree;
318       fAODfile           = o.fAODfile;
319       mcHeader           = o.mcHeader;
320       rndm               = o.rndm;
321       fInputEntries      = o.fInputEntries;
322       fAODPathArray      = o.fAODPathArray;
323       fAODEntriesArray   = o.fAODEntriesArray;
324       fAODPath           = o.fAODPath;
325       fAODEntries        = o.fAODEntries;
326       fAODEntriesSum     = o.fAODEntriesSum;
327       fAODEntriesMax     = o.fAODEntriesMax;
328       fOfflineTrgMask    = o.fOfflineTrgMask;
329       fMinContribVtx     = o.fMinContribVtx;
330       fVtxZMin           = o.fVtxZMin;
331       fVtxZMax           = o.fVtxZMax;
332       fEvtClassMin       = o.fEvtClassMin;
333       fEvtClassMax       = o.fEvtClassMax;
334       fCentMin           = o.fCentMin;
335       fCentMax           = o.fCentMax;
336       fNInputTracksMin   = o.fNInputTracksMin;
337       fNInputTracksMax   = o.fNInputTracksMax;
338       fTrackBranch       = o.fTrackBranch;
339       fMCparticlesBranch = o.fMCparticlesBranch;
340       fJetBranch         = o.fJetBranch;
341       fFileId            = o.fFileId;
342       fAODEntry          = o.fAODEntry;
343       fCountEvents       = o.fCountEvents;
344       fEmbedMode         = o.fEmbedMode;
345       fEvtSelecMode      = o.fEvtSelecMode;
346       fEvtSelMinJetPt    = o.fEvtSelMinJetPt;
347       fEvtSelMaxJetPt    = o.fEvtSelMaxJetPt;
348       fEvtSelMinJetEta   = o.fEvtSelMinJetEta;
349       fEvtSelMaxJetEta   = o.fEvtSelMaxJetEta;
350       fEvtSelMinJetPhi   = o.fEvtSelMinJetPhi;
351       fEvtSelMaxJetPhi   = o.fEvtSelMaxJetPhi;
352       fToyMinNbOfTracks  = o.fToyMinNbOfTracks;
353       fToyMaxNbOfTracks  = o.fToyMaxNbOfTracks;
354       fToyMinTrackPt     = o.fToyMinTrackPt;
355       fToyMaxTrackPt     = o.fToyMaxTrackPt;
356       fToyDistributionTrackPt = o.fToyDistributionTrackPt;
357       fToyMinTrackEta    = o.fToyMinTrackEta;
358       fToyMaxTrackEta    = o.fToyMaxTrackEta;
359       fToyMinTrackPhi    = o.fToyMinTrackPhi;
360       fToyMaxTrackPhi    = o.fToyMaxTrackPhi;
361       fToyFilterMap      = o.fToyFilterMap;
362       fTrackFilterMap    = o.fTrackFilterMap;
363       fNPtHard           = o.fNPtHard;
364       fPtHard            = o.fPtHard;
365       fPtHardBin         = o.fPtHardBin;
366       fAODJets           = o.fAODJets;
367       fNevents           = o.fNevents;
368       fXsection          = o.fXsection;
369       fAvgTrials         = o.fAvgTrials;
370       fHistList          = o.fHistList;
371       fHistEvtSelection  = o.fHistEvtSelection;
372       fh1Xsec            = o.fh1Xsec;
373       fh1Trials          = o.fh1Trials;
374       fh1TrialsEvtSel    = o.fh1TrialsEvtSel;
375       fh2PtHard          = o.fh2PtHard;
376       fh2PtHardEvtSel    = o.fh2PtHardEvtSel;
377       fh2PtHardTrials    = o.fh2PtHardTrials;
378       fh1TrackPt         = o.fh1TrackPt;
379       fh2TrackEtaPhi     = o.fh2TrackEtaPhi;
380       fh1TrackN          = o.fh1TrackN;
381       fh1JetPt           = o.fh1JetPt;
382       fh2JetEtaPhi       = o.fh2JetEtaPhi;
383       fh1JetN            = o.fh1JetN;
384       fh1MCTrackPt       = o.fh1MCTrackPt;
385       fh2MCTrackEtaPhi   = o.fh2MCTrackEtaPhi;
386       fh1MCTrackN        = o.fh1MCTrackN;
387       fh1AODfile         = o.fh1AODfile;
388       fh2AODevent        = o.fh2AODevent;
389    }
390
391    return *this;
392 }
393
394
395 //__________________________________________________________________________
396 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
397 {
398    // destructor
399    delete rndm;
400 }
401
402
403 //__________________________________________________________________________
404 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
405 {
406    // create output objects
407    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
408    AliLog::SetClassDebugLevel("AliAnalysisTaskFastEmbedding", AliLog::kInfo);
409
410    OpenFile(1);
411    if(!fHistList) fHistList = new TList();
412    fHistList->SetOwner(kTRUE);
413    
414
415    // set seed
416    rndm = new TRandom3();
417    Int_t id = GetJobID();
418    if(id>-1) rndm->SetSeed(id);
419    else      rndm->SetSeed();   // a TTUID is generated and used for seed
420    AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
421
422
423
424    // embed mode with AOD
425    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
426
427       // open input AOD
428       fFileId = OpenAODfile();
429       fNevents = 0; // force to open another aod in UserExec()
430       if(fFileId<0){
431          AliError("");
432          PostData(1, fHistList);
433          return;
434       }
435    } //end: embed mode with AOD
436
437
438    // connect output aod
439    // create a new branch for extra tracks
440    fAODout = AODEvent();
441    if(!fAODout){
442       AliError("Output AOD not found.");
443       PostData(1, fHistList);
444       return;
445    }
446    if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
447       AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
448       TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
449       tracks->SetName(fTrackBranch.Data());
450       AddAODBranch("TClonesArray", &tracks);
451    }
452    // create new branch for extra mcparticle if available as input
453    if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
454       AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
455       TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
456       mcparticles->SetName(fMCparticlesBranch.Data());
457       AddAODBranch("TClonesArray", &mcparticles);
458    }
459
460
461
462    //qa histograms
463    Bool_t oldStatus = TH1::AddDirectoryStatus();
464    TH1::AddDirectory(kFALSE);
465    
466    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
467    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
468    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
469    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
470    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
471    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
472    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
473    
474    fh1Xsec         = new TProfile("fh1Xsec","xsec from pyxsec.root;p_{T,hard} bin;<#sigma>",fNPtHard+1,-1.5,fNPtHard-0.5);
475    fh1Trials       = new TH1F("fh1Trials","trials (simulation) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
476    fh1TrialsEvtSel = new TH1F("fh1TrialsEvtSel","trials (event selection) from pyxsec.root;p_{T,hard} bin;#sum{ntrials}",fNPtHard+1,-1.5,fNPtHard-0.5);
477    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);
478    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);
479    fh2PtHardTrials = new TH2F("fh2PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard} bin;#sum{p_{T,hard}}",fNPtHard+1,-1.5,fNPtHard-0.5,350,-.5,349.5);
480    
481    fHistList->Add(fHistEvtSelection);
482    fHistList->Add(fh1Xsec);
483    fHistList->Add(fh1Trials);
484    fHistList->Add(fh1TrialsEvtSel);
485    fHistList->Add(fh2PtHard);
486    fHistList->Add(fh2PtHardEvtSel);
487    fHistList->Add(fh2PtHardTrials);
488
489    fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
490    fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
491    fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
492
493    fHistList->Add(fh1TrackPt);
494    fHistList->Add(fh2TrackEtaPhi);
495    fHistList->Add(fh1TrackN);
496    
497    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
498       
499       fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 120, 0., 120.);
500       fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
501       fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
502       
503       fHistList->Add(fh1JetPt);
504       fHistList->Add(fh2JetEtaPhi);
505       fHistList->Add(fh1JetN);
506    }
507
508
509    if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){ 
510
511       fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
512       fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
513       fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
514       
515       fHistList->Add(fh1MCTrackPt);
516       fHistList->Add(fh2MCTrackEtaPhi);
517       fHistList->Add(fh1MCTrackN);
518       
519    }
520    
521    fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 2300, -0.5, 2299.5);
522    fh2AODevent = new TH2I("fh1AODevent","selected events;file;event", 2500,-.5,2499.5,5000,-.5,4999.5);
523    fHistList->Add(fh1AODfile);
524    fHistList->Add(fh2AODevent);
525
526    // =========== Switch on Sumw2 for all histos ===========
527    for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
528       TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
529       if (h1){
530          h1->Sumw2();
531          continue;
532       }
533    }
534
535    TH1::AddDirectory(oldStatus);
536
537    PostData(1, fHistList);
538 }
539
540
541 //__________________________________________________________________________
542 void AliAnalysisTaskFastEmbedding::Init()
543 {
544    // Initialization
545    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
546
547 }
548
549 //__________________________________________________________________________
550 Bool_t AliAnalysisTaskFastEmbedding::UserNotify()
551 {
552    // User defined Notify(), called once
553    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserNotify()");
554
555    // get total nb of events in tree (of this subjob)
556    AliInputEventHandler* inputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
557    fInputEntries = inputHandler->GetTree()->GetEntriesFast();
558    AliInfo(Form("Total nb. of events: %d", fInputEntries));
559
560    return kTRUE;
561
562 }
563
564
565 //__________________________________________________________________________
566 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
567 {
568    // Analysis, called once per event
569    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
570
571    if(!fAODout){
572       AliError("Need output AOD, but is not connected."); 
573       PostData(1, fHistList);
574       return;
575    }
576
577    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
578    if (!fESD) {
579       AliError("ESD not available");
580       PostData(1, fHistList);
581       return;
582    }
583
584    // -- event selection --
585    fHistEvtSelection->Fill(1); // number of events before event selection
586
587    // physics selection
588    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
589    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
590    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
591       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
592       fHistEvtSelection->Fill(2);
593       PostData(1, fHistList);
594       return;
595    } 
596
597    // vertex selection
598    AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
599    Int_t nTracksPrim = primVtx->GetNContributors();
600    if ((nTracksPrim < fMinContribVtx) ||
601          (primVtx->GetZ() < fVtxZMin) ||
602          (primVtx->GetZ() > fVtxZMax) ){
603       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
604       fHistEvtSelection->Fill(3);
605       PostData(1, fHistList);
606       return;
607    }
608
609    /** takes wrong value, since jet service tasks runs later
610    // event class selection (from jet helper task)
611    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
612    if(fDebug) Printf("Event class %d", eventClass);
613    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
614    fHistEvtSelection->Fill(4);
615    PostData(1, fHistList);
616    return;
617    }*/
618
619    // centrality selection
620    AliCentrality *cent = 0x0;
621    Float_t centValue = 0.; 
622    if(fESD) cent = fESD->GetCentrality();
623    if(cent) centValue = cent->GetCentralityPercentile("V0M");
624    if(fDebug) printf("centrality: %f\n", centValue);
625    if (centValue < fCentMin || centValue > fCentMax){
626       fHistEvtSelection->Fill(4);
627       PostData(1, fHistList);
628       return;
629    }
630
631
632    /*  ** not implemented **
633    // multiplicity due to input tracks
634    Int_t nInputTracks = GetNInputTracks();
635
636    if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
637    fHistEvtSelection->Fill(5);
638    PostData(1, fHistList);
639    return;
640    }
641    */
642
643    fHistEvtSelection->Fill(0); // accepted events  
644    // -- end event selection --
645
646    // connect aod out
647    TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
648    if(!tracks){
649       AliError("Extra track branch not found in output.");
650       PostData(1, fHistList);
651       return;
652    }
653    tracks->Delete();
654    Int_t nAODtracks=0;
655
656    TRef dummy;
657
658    // === embed mode with AOD ===
659    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
660       if(!fAODevent){
661          AliError("Need input AOD, but is not connected."); 
662          PostData(1, fHistList);
663          return;
664       }
665
666
667
668       Bool_t useEntry = kFALSE;
669       while(!useEntry){  // protection need, if no event fulfills requirement
670
671          fAODEntry++; // go to next event 
672          fCountEvents++;
673          if(fAODEntry>=fNevents) fAODEntry=0; 
674
675          // open new aod file
676          if(fCountEvents>=fNevents){ 
677             if(!fAODPathArray){
678                AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
679             } 
680             else {
681                AliDebug(AliLog::kDebug, "Select new AOD file ...");
682
683                fFileId = OpenAODfile();
684                if(fFileId<0) {
685                   PostData(1, fHistList);
686                   return;
687                }
688                fh1AODfile->Fill(fFileId);
689                if(fAODEntries!=fNevents){
690                   AliError("File with list of AODs and file with nb. of entries inconsistent");
691                   PostData(1, fHistList);
692                   return;
693                }
694             }
695          }
696
697
698
699          // get next event
700          fAODtree->GetEvent(fAODEntry);
701
702          // get pt hard
703          if(mcHeader){
704             fPtHard = mcHeader->GetPtHard();
705             GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
706             AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
707          }
708          else{
709             AliDebug(AliLog::kDebug,"No mcHeader");
710          }
711          fPtHardBin = GetPtHardBin(fPtHard);
712
713          //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials); 
714
715          // fill event variables for each event
716          fh1Xsec->Fill(fPtHardBin,fXsection);
717          fh2PtHard->Fill(fPtHardBin,fPtHard);
718
719          fh1Trials->Fill(fPtHardBin, fAvgTrials);
720          fh1TrialsEvtSel->Fill(fPtHardBin);
721          fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
722
723          // jet pt selection
724          if(fEvtSelecMode==kEventsJetPt){
725             Int_t nJets = fAODJets->GetEntries();
726             //Printf("nb. jets: %d",nJets);
727             for(Int_t ij=0; ij<nJets; ++ij){
728                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
729                if(!jet) continue;
730
731                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
732                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
733                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
734                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
735                   useEntry = kTRUE;
736                   break;
737                } 
738             }
739          }
740
741          // no selection
742          if(fEvtSelecMode==kEventsAll){
743             useEntry = kTRUE;
744          }
745
746       }
747       AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
748
749       fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
750       fh2AODevent->Fill(fFileId,fAODEntry);
751
752       TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
753       TClonesArray *mcpartOUT = 0x0;
754       if(mcpartIN){
755          mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
756          if(!mcpartOUT){
757             AliError("Extra MC particles branch not found in output.");
758             PostData(1, fHistList);
759             return;
760          }
761          mcpartOUT->Delete();
762       } else {
763          AliInfo("No extra MC particles found.");
764       }
765
766
767       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
768          // loop over input tracks
769          // add to output aod
770          Int_t nTracks = 0;
771          Int_t nJets = fAODJets->GetEntries();
772          Int_t nSelectedJets = 0;
773          AliAODJet *leadJet = 0x0; // used for jet tracks only
774
775          if(fEmbedMode==kAODFull){
776             nTracks = fAODevent->GetNumberOfTracks();
777
778             for(Int_t ij=0; ij<nJets; ++ij){
779                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
780                if(!jet) continue;
781                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
782                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
783                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
784                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
785
786                   fh1JetPt->Fill(jet->Pt());
787                   fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
788                   nSelectedJets++;
789
790                }
791             }                              
792             fh1JetN->Fill(nSelectedJets);
793          }
794
795          if(fEmbedMode==kAODJetTracks){
796             // look for leading jet within selection
797             // get reference tracks for this jet
798             Float_t leadJetPt = 0.;
799             for(Int_t ij=0; ij<nJets; ++ij){
800                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
801                if(!jet) continue;
802                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
803                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
804                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
805                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
806                   if(jet->Pt()>leadJetPt){
807                      leadJetPt = jet->Pt();
808                      leadJet = jet;
809                   } 
810                }
811             }
812             if(leadJet){
813                nTracks = leadJet->GetRefTracks()->GetEntriesFast();
814                fh1JetPt->Fill(leadJet->Pt());
815                fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
816                fh1JetN->Fill(1);                                   
817             }
818          }
819
820          fh1TrackN->Fill((Float_t)nTracks);
821
822          for(Int_t it=0; it<nTracks; ++it){
823             AliAODTrack *tmpTr = 0x0;
824             if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
825             if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
826             if(!tmpTr) continue;
827
828             tmpTr->SetStatus(AliESDtrack::kEmbedded);
829
830             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
831             dummy = (*tracks)[nAODtracks-1];
832
833             if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
834                fh1TrackPt->Fill(tmpTr->Pt());
835                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
836             }
837          }
838
839          if(mcpartIN){
840             Int_t nMCpart = mcpartIN->GetEntriesFast();
841
842             Int_t nPhysicalPrimary=0;
843             Int_t nAODmcpart=0;
844             for(Int_t ip=0; ip<nMCpart; ++ip){
845                AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
846                if(!tmpPart) continue;
847                if(fEmbedMode==kAODJetTracks){
848                   // jet track? else continue
849                   // not implemented yet
850                   continue;
851                } 
852
853                if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99.  && tmpPart->Pt()>0.){
854                  new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
855                  dummy = (*mcpartOUT)[nAODmcpart-1];
856                  
857                  fh1MCTrackPt->Fill(tmpPart->Pt());
858                  fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
859                  nPhysicalPrimary++;
860                }
861             }
862             fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
863
864          }
865       } // end: embed all tracks or jet tracks
866
867
868       if(fEmbedMode==kAODJet4Mom){
869
870          // loop over jets
871          Int_t nJets = fAODJets->GetEntries();
872          fh1TrackN->Fill((Float_t)nJets);
873          for(Int_t ij=0; ij<nJets; ++ij){
874             AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
875             if(!jet) continue;
876             AliAODTrack *tmpTr = (AliAODTrack*)jet;
877             tmpTr->SetFlags(AliESDtrack::kEmbedded);
878
879             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
880             dummy = (*tracks)[nAODtracks-1]; 
881
882             fh1TrackPt->Fill(tmpTr->Pt());
883             fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
884          }
885
886       } // end: embed jets as 4-momenta
887
888
889    } //end: embed mode with AOD
890
891
892    // === embed mode with toy ===
893    if(fEmbedMode==kToyTracks){
894       Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
895
896       fh1TrackN->Fill((Float_t)nT);
897
898       for(Int_t i=0; i<nT; ++i){
899
900          Double_t pt = -1.;
901          if(fToyMinTrackPt!=fToyMaxTrackPt){
902             if(fToyDistributionTrackPt==0){
903                pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
904             } else {
905                while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
906                   pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
907                   pt += fToyMinTrackPt;
908                }
909             }
910          } else {
911             pt = fToyMinTrackPt;
912          }
913          Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
914          Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
915          phi = TVector2::Phi_0_2pi(phi);
916
917          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
918
919          Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
920          Float_t mom[3];
921          mom[0] = pt;
922          mom[1] = phi;
923          mom[2] = theta;
924
925          Float_t xyz[3];
926          xyz[0] = -999.;
927          xyz[1] = -999.;
928          xyz[2] = -999.;
929
930          AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
931          -999,   // label
932          mom,    // momentum[3]
933          kFALSE, // cartesian
934          xyz,    // position
935          kFALSE, // DCA
936          NULL,   // covMatrix,
937          -99,    // charge
938          0,      // itsClusMap
939          NULL,   // pid 
940          NULL,   // prodVertex
941          kFALSE, // used for vertex fit
942          kFALSE, // used for prim vtx fit
943          AliAODTrack::kUndef, // type
944          fToyFilterMap,  // select info
945          -999.    // chi2 per NDF
946          );
947          tmpTr->SetFlags(AliESDtrack::kEmbedded);
948
949          new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
950          dummy = (*tracks)[nAODtracks-1];
951
952          fh1TrackPt->Fill(pt);
953          fh2TrackEtaPhi->Fill(eta,phi);
954
955          delete tmpTr;
956       }
957    } //end: embed mode with toy
958
959
960    PostData(1, fHistList);
961 }
962
963
964 //__________________________________________________________________________
965 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
966 {
967    // terminate
968    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
969
970    if(fAODfile && fAODfile->IsOpen())
971    fAODfile->Close();  
972
973 }
974
975 //__________________________________________________________________________
976 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
977 {
978    // gives back the alien subjob id, if available, else -1
979
980    Int_t id=-1;
981
982    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
983    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
984
985    if(env && strlen(env)){
986       id= atoi(env);
987       AliInfo(Form("Job index %d", id));
988    }
989    else{
990       AliInfo("Job index not found. Okay if running locally.");
991    }
992
993    return id;
994 }
995
996 //__________________________________________________________________________
997 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
998 {
999    // select an AOD file from fAODPathArray
1000
1001    Int_t nFiles = fAODPathArray->GetEntries();
1002
1003    // choose randomly file and event
1004    Int_t n = -1;
1005    Float_t tmpProp = -1.;
1006    Int_t pendingEvents = fInputEntries-Entry();
1007    //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1008    if(fAODEntriesArray){
1009      while(rndm->Rndm()>tmpProp){
1010        n = rndm->Integer(nFiles);
1011        fAODEntries = fAODEntriesArray->At(n);
1012        Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1013        tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1014      }
1015      fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1016    }
1017    else {
1018       AliWarning("Number of entries in extra AODs not set!");
1019       n = rndm->Integer(nFiles);
1020       fAODEntry = 0;
1021    }
1022
1023    AliInfo(Form("Select AOD file %d", n));
1024    TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1025    if(!objStr){
1026       AliError("Could not get path of aod file.");
1027       return -1;
1028    }
1029    fAODPath = objStr->GetString();
1030
1031    return n;
1032 }
1033
1034 //__________________________________________________________________________
1035 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1036 {
1037    // select and open an AOD file from fAODPathArray
1038
1039    if(fAODPathArray) fFileId = SelectAODfile();
1040    if(fFileId<0) return -1;
1041
1042    if (!gGrid) {
1043      AliInfo("Trying to connect to AliEn ...");
1044      TGrid::Connect("alien://");
1045    }
1046
1047    TDirectory *owd = gDirectory;
1048    if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
1049    fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1050    owd->cd();
1051    if(!fAODfile){
1052
1053       fFileId = -1;
1054       if(fAODPathArray){
1055          if(trial<=3){ 
1056             AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1057             fFileId = OpenAODfile(trial+1);
1058          } else {
1059             AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1060             return -1;
1061          }
1062       } else {
1063          AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1064       }
1065
1066       return fFileId;
1067    }
1068
1069    fAODtree = (TTree*)fAODfile->Get("aodTree");
1070
1071    if(!fAODtree){
1072       AliError("AOD tree not found.");
1073       return -1;
1074    }
1075
1076    /*
1077       fAODtree->SetBranchStatus("*",0);
1078       fAODtree->SetBranchStatus("header",1);
1079       fAODtree->SetBranchStatus("tracks*",1);
1080       fAODtree->SetBranchStatus("jets*",1);
1081       fAODtree->SetBranchStatus("clusters*",1);
1082       fAODtree->SetBranchStatus("mcparticles*",1);
1083       fAODtree->SetBranchStatus("mcHeader*",1);
1084    */
1085
1086    delete fAODevent;
1087    fAODevent = new AliAODEvent();
1088    fAODevent->ReadFromTree(fAODtree);
1089
1090
1091    // fetch header, jets, etc. from new file
1092    fNevents = fAODtree->GetEntries();
1093    mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1094
1095    if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1096    else                    fAODJets = fAODevent->GetJets();
1097    if(!fAODJets){
1098       AliError("Could not find jets in AOD. Check jet branch when indicated.");
1099       return -1;
1100    }
1101
1102    TFile *curfile = fAODtree->GetCurrentFile();
1103    if (!curfile) {
1104       AliError("No current file.");
1105       return -1;
1106    }
1107
1108    Float_t trials = 1.;
1109    Float_t xsec = 0.;
1110    PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1111    fXsection = xsec;
1112
1113    // construct a poor man average trials 
1114    Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1115    if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1116
1117    if(fFileId>=0){
1118       AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1119    }
1120
1121    fCountEvents=0; // new file, reset counter
1122
1123    return fFileId;  // file position in AOD path array, if array available
1124 }
1125
1126
1127 //____________________________________________________________________________
1128 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1129
1130    // static stored, available for other tasks in train
1131
1132    static Float_t ptHard = -1.;
1133    if(bSet) ptHard = newValue;
1134
1135    return ptHard;
1136 }
1137
1138 //____________________________________________________________________________
1139 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1140
1141    // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1142
1143    const Int_t nBins = 10;
1144    Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1145
1146    Int_t bin = -1;
1147    while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1148       bin++;
1149    }
1150
1151    return bin;
1152 }
1153
1154 //____________________________________________________________________________
1155 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1156    //
1157    // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1158    // This should provide the path to the AOD/ESD file
1159    // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks) 
1160
1161    TString file(currFile);  
1162    fXsec = 0;
1163    fTrials = 1;
1164
1165    if(file.Contains("root_archive.zip#")){
1166      Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
1167      Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1168      Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
1169       file.Replace(pos+1,pos2-pos1,"");
1170       //      file.Replace(pos+1,20,"");
1171    }
1172    else {
1173       // not an archive take the basename....
1174       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1175    }
1176
1177    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
1178    if(!fxsec){
1179       // next trial fetch the histgram file
1180       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1181       if(!fxsec){
1182          // not a severe condition but inciate that we have no information
1183          AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1184          return kFALSE;
1185       }
1186       else{
1187          // find the tlist we want to be independtent of the name so use the Tkey
1188          TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1189          if(!key){
1190             fxsec->Close();
1191             return kFALSE;
1192          }
1193          TList *list = dynamic_cast<TList*>(key->ReadObj());
1194          if(!list){
1195             fxsec->Close();
1196             return kFALSE;
1197          }
1198          fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1199          fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1200          fxsec->Close();
1201       }
1202    } // no tree pyxsec.root
1203    else {
1204       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1205       if(!xtree){
1206          fxsec->Close();
1207          return kFALSE;
1208       }
1209       UInt_t   ntrials  = 0;
1210       Double_t  xsection  = 0;
1211       xtree->SetBranchAddress("xsection",&xsection);
1212       xtree->SetBranchAddress("ntrials",&ntrials);
1213       xtree->GetEntry(0);
1214       fTrials = ntrials;
1215       fXsec = xsection;
1216       fxsec->Close();
1217    }
1218    return kTRUE;
1219 }
1220
1221