JETAN/AliAnalysisTaskFastEmbedding.cxx

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