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