Switching from CMAKE_SOURCE_DIR to AliRoot_SOURCE_DIR
[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)  {
829               tmpTr = dynamic_cast<AliAODTrack*>(fAODevent->GetTrack(it));
830               if(!tmpTr) AliFatal("Not a standard AOD");
831             }
832             if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
833             if(!tmpTr) continue;
834             Double_t rd=rndm->Uniform(0.,1.);
835             if(rd>fExtraEffPb) continue; 
836             tmpTr->SetStatus(AliESDtrack::kEmbedded);
837
838             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr); 
839             dummy = (*tracks)[nAODtracks-1];
840
841             if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
842                fh1TrackPt->Fill(tmpTr->Pt());
843                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
844             }
845          }
846
847          if(mcpartIN){
848             Int_t nMCpart = mcpartIN->GetEntriesFast();
849
850             Int_t nPhysicalPrimary=0;
851             Int_t nAODmcpart=0;
852             for(Int_t ip=0; ip<nMCpart; ++ip){
853                AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
854                if(!tmpPart) continue;
855                if(fEmbedMode==kAODJetTracks){
856                   // jet track? else continue
857                   // not implemented yet
858                   continue;
859                } 
860
861                if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99.  && tmpPart->Pt()>0.){
862                  new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
863                  dummy = (*mcpartOUT)[nAODmcpart-1];
864
865                  if(fDebug>10) printf("added track %d with pT=%.2f to extra branch\n",nAODmcpart,tmpPart->Pt());
866                  
867                  fh1MCTrackPt->Fill(tmpPart->Pt());
868                  fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
869                  nPhysicalPrimary++;
870
871                }
872             }
873             fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
874
875          }
876       } // end: embed all tracks or jet tracks
877
878
879       if(fEmbedMode==kAODJet4Mom){
880
881          // loop over jets
882          Int_t nJets = fAODJets->GetEntries();
883          fh1TrackN->Fill((Float_t)nJets);
884          for(Int_t ij=0; ij<nJets; ++ij){
885             AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
886             if(!jet) continue;
887             AliAODTrack *tmpTr = (AliAODTrack*)jet;
888             tmpTr->SetFlags(AliESDtrack::kEmbedded);
889
890             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
891             dummy = (*tracks)[nAODtracks-1]; 
892
893             fh1TrackPt->Fill(tmpTr->Pt());
894             fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
895          }
896
897       } // end: embed jets as 4-momenta
898
899
900    } //end: embed mode with AOD
901
902
903    // === embed mode with toy ===
904    if(fEmbedMode==kToyTracks){
905       Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
906
907       fh1TrackN->Fill((Float_t)nT);
908
909       for(Int_t i=0; i<nT; ++i){
910
911          Double_t pt = -1.;
912          if(fToyMinTrackPt!=fToyMaxTrackPt){
913             if(fToyDistributionTrackPt==0){
914                pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
915             } else {
916                while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
917                   pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
918                   pt += fToyMinTrackPt;
919                }
920             }
921          } else {
922             pt = fToyMinTrackPt;
923          }
924          Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
925          Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
926          phi = TVector2::Phi_0_2pi(phi);
927
928          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
929
930          Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
931          Float_t mom[3];
932          mom[0] = pt;
933          mom[1] = phi;
934          mom[2] = theta;
935
936          Float_t xyz[3];
937          xyz[0] = -999.;
938          xyz[1] = -999.;
939          xyz[2] = -999.;
940
941          AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
942          -999,   // label
943          mom,    // momentum[3]
944          kFALSE, // cartesian
945          xyz,    // position
946          kFALSE, // DCA
947          NULL,   // covMatrix,
948          -99,    // charge
949          0,      // itsClusMap
950          //         NULL,   // pid 
951          NULL,   // prodVertex
952          kFALSE, // used for vertex fit
953          kFALSE, // used for prim vtx fit
954          AliAODTrack::kUndef, // type
955          fToyFilterMap,  // select info
956          -999.    // chi2 per NDF
957          );
958          tmpTr->SetFlags(AliESDtrack::kEmbedded);
959
960          new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
961          dummy = (*tracks)[nAODtracks-1];
962
963          fh1TrackPt->Fill(pt);
964          fh2TrackEtaPhi->Fill(eta,phi);
965
966          delete tmpTr;
967       }
968    } //end: embed mode with toy
969
970
971    PostData(1, fHistList);
972 }
973
974
975 //__________________________________________________________________________
976 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
977 {
978    // terminate
979    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
980
981    if(fAODfile && fAODfile->IsOpen())
982    fAODfile->Close();  
983
984 }
985
986 //__________________________________________________________________________
987 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
988 {
989    // gives back the alien subjob id, if available, else -1
990
991    Int_t id=-1;
992
993    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
994    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
995
996    if(env && strlen(env)){
997       id= atoi(env);
998       AliInfo(Form("Job index %d", id));
999    }
1000    else{
1001       AliInfo("Job index not found. Okay if running locally.");
1002    }
1003
1004    return id;
1005 }
1006
1007 //__________________________________________________________________________
1008 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile()
1009 {
1010    // select an AOD file from fAODPathArray
1011
1012    Int_t nFiles = fAODPathArray->GetEntries();
1013
1014    // choose randomly file and event
1015    Int_t n = -1;
1016    Float_t tmpProp = -1.;
1017    Int_t pendingEvents = fInputEntries-Entry();
1018    //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1019    if(fAODEntriesArray){
1020      while(rndm->Rndm()>tmpProp){
1021        n = rndm->Integer(nFiles);
1022        fAODEntries = fAODEntriesArray->At(n);
1023        Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1024        tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1025      }
1026      fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1027    }
1028    else {
1029       AliWarning("Number of entries in extra AODs not set!");
1030       n = rndm->Integer(nFiles);
1031       fAODEntry = 0;
1032    }
1033
1034    AliInfo(Form("Select AOD file %d", n));
1035    TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1036    if(!objStr){
1037       AliError("Could not get path of aod file.");
1038       return -1;
1039    }
1040    fAODPath = objStr->GetString();
1041
1042    return n;
1043 }
1044
1045 //__________________________________________________________________________
1046 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial)
1047 {
1048    // select and open an AOD file from fAODPathArray
1049
1050    if(fAODPathArray) fFileId = SelectAODfile();
1051    if(fFileId<0) return -1;
1052
1053    if (!gGrid) {
1054      AliInfo("Trying to connect to AliEn ...");
1055      TGrid::Connect("alien://");
1056    }
1057
1058    TDirectory *owd = gDirectory;
1059    if (fAODfile && fAODfile->IsOpen()) fAODfile->Close();
1060    fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1061    owd->cd();
1062    if(!fAODfile){
1063
1064       fFileId = -1;
1065       if(fAODPathArray){
1066          if(trial<=3){ 
1067             AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1068             fFileId = OpenAODfile(trial+1);
1069          } else {
1070             AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1071             return -1;
1072          }
1073       } else {
1074          AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1075       }
1076
1077       return fFileId;
1078    }
1079
1080    fAODtree = (TTree*)fAODfile->Get("aodTree");
1081
1082    if(!fAODtree){
1083       AliError("AOD tree not found.");
1084       return -1;
1085    }
1086
1087    /*
1088       fAODtree->SetBranchStatus("*",0);
1089       fAODtree->SetBranchStatus("header",1);
1090       fAODtree->SetBranchStatus("tracks*",1);
1091       fAODtree->SetBranchStatus("jets*",1);
1092       fAODtree->SetBranchStatus("clusters*",1);
1093       fAODtree->SetBranchStatus("mcparticles*",1);
1094       fAODtree->SetBranchStatus("mcHeader*",1);
1095    */
1096
1097    delete fAODevent;
1098    fAODevent = new AliAODEvent();
1099    fAODevent->ReadFromTree(fAODtree);
1100
1101
1102    // fetch header, jets, etc. from new file
1103    fNevents = fAODtree->GetEntries();
1104    mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1105
1106    if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1107    else                    fAODJets = fAODevent->GetJets();
1108    if(!fAODJets){
1109       AliError("Could not find jets in AOD. Check jet branch when indicated.");
1110       return -1;
1111    }
1112
1113    TFile *curfile = fAODtree->GetCurrentFile();
1114    if (!curfile) {
1115       AliError("No current file.");
1116       return -1;
1117    }
1118
1119    Float_t trials = 1.;
1120    Float_t xsec = 0.;
1121    PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1122    fXsection = xsec;
1123
1124    // construct a poor man average trials 
1125    Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1126    if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1127
1128    if(fFileId>=0){
1129       AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1130    }
1131
1132    fCountEvents=0; // new file, reset counter
1133
1134    return fFileId;  // file position in AOD path array, if array available
1135 }
1136
1137
1138 //____________________________________________________________________________
1139 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1140
1141    // static stored, available for other tasks in train
1142
1143    static Float_t ptHard = -1.;
1144    if(bSet) ptHard = newValue;
1145
1146    return ptHard;
1147 }
1148
1149 //____________________________________________________________________________
1150 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1151
1152    // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1153
1154    const Int_t nBins = 10;
1155    Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1156
1157    Int_t bin = -1;
1158    while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1159       bin++;
1160    }
1161
1162    return bin;
1163 }
1164
1165 //____________________________________________________________________________
1166 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1167    //
1168    // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1169    // This should provide the path to the AOD/ESD file
1170    // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks) 
1171
1172    TString file(currFile);  
1173    fXsec = 0;
1174    fTrials = 1;
1175
1176    if(file.Contains("root_archive.zip#")){
1177      Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
1178      Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1179      Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
1180       file.Replace(pos+1,pos2-pos1,"");
1181       //      file.Replace(pos+1,20,"");
1182    }
1183    else {
1184       // not an archive take the basename....
1185       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1186    }
1187
1188    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
1189    if(!fxsec){
1190       // next trial fetch the histgram file
1191       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1192       if(!fxsec){
1193          // not a severe condition but inciate that we have no information
1194          AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1195          return kFALSE;
1196       }
1197       else{
1198          // find the tlist we want to be independtent of the name so use the Tkey
1199          TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1200          if(!key){
1201             fxsec->Close();
1202             return kFALSE;
1203          }
1204          TList *list = dynamic_cast<TList*>(key->ReadObj());
1205          if(!list){
1206             fxsec->Close();
1207             return kFALSE;
1208          }
1209          fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1210          fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1211          fxsec->Close();
1212       }
1213    } // no tree pyxsec.root
1214    else {
1215       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1216       if(!xtree){
1217          fxsec->Close();
1218          return kFALSE;
1219       }
1220       UInt_t   ntrials  = 0;
1221       Double_t  xsection  = 0;
1222       xtree->SetBranchAddress("xsection",&xsection);
1223       xtree->SetBranchAddress("ntrials",&ntrials);
1224       xtree->GetEntry(0);
1225       fTrials = ntrials;
1226       fXsec = xsection;
1227       fxsec->Close();
1228    }
1229    return kTRUE;
1230 }
1231
1232