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
 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    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
 568
 569    if(!fAODout){
 570       AliError("Need output AOD, but is not connected.");
 571       PostData(1, fHistList);
 572       return;
 573    }
 574
 575    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 576    if (!fESD) {
 577       AliError("ESD not available");
 578       PostData(1, fHistList);
 579       return;
 580    }
 581
 582    // -- event selection --
 583    fHistEvtSelection->Fill(1); // number of events before event selection
 584
 585    // physics selection
 586    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 587    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 588    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 589       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 590       fHistEvtSelection->Fill(2);
 591       PostData(1, fHistList);
 592       return;
 593    }
 594
 595    // vertex selection
 596    AliAODVertex* primVtx = fAODout->GetPrimaryVertex();
 597    Int_t nTracksPrim = primVtx->GetNContributors();
 598    if ((nTracksPrim < fMinContribVtx) ||
 599          (primVtx->GetZ() < fVtxZMin) ||
 600          (primVtx->GetZ() > fVtxZMax) ){
 601       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 602       fHistEvtSelection->Fill(3);
 603       PostData(1, fHistList);
 604       return;
 605    }
 606
 607    /** takes wrong value, since jet service tasks runs later
 608    // event class selection (from jet helper task)
 609    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 610    if(fDebug) Printf("Event class %d", eventClass);
 611    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 612    fHistEvtSelection->Fill(4);
 613    PostData(1, fHistList);
 614    return;
 615    }*/
 616
 617    // centrality selection
 618    AliCentrality *cent = 0x0;
 619    Float_t centValue = 0.;
 620    if(fESD) cent = fESD->GetCentrality();
 621    if(cent) centValue = cent->GetCentralityPercentile("V0M");
 622    if(fDebug) printf("centrality: %f\n", centValue);
 623    if (centValue < fCentMin || centValue > fCentMax){
 624       fHistEvtSelection->Fill(4);
 625       PostData(1, fHistList);
 626       return;
 627    }
 628
 629
 630    /*  ** not implemented **
 631    // multiplicity due to input tracks
 632    Int_t nInputTracks = GetNInputTracks();
 633
 634    if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
 635    fHistEvtSelection->Fill(5);
 636    PostData(1, fHistList);
 637    return;
 638    }
 639    */
 640
 641    fHistEvtSelection->Fill(0); // accepted events
 642    // -- end event selection --
 643
 644    // connect aod out
 645    TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
 646    if(!tracks){
 647       AliError("Extra track branch not found in output.");
 648       PostData(1, fHistList);
 649       return;
 650    }
 651    tracks->Delete();
 652    Int_t nAODtracks=0;
 653
 654    TRef dummy;
 655
 656    // === embed mode with AOD ===
 657    if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 658       if(!fAODevent){
 659          AliError("Need input AOD, but is not connected.");
 660          PostData(1, fHistList);
 661          return;
 662       }
 663
 664
 665
 666       Bool_t useEntry = kFALSE;
 667       while(!useEntry){  // protection need, if no event fulfills requierment
 668
 669          fAODEntry++; // go to next event
 670          fCountEvents++;
 671          if(fAODEntry>=fNevents) fAODEntry=0;
 672
 673          // open new aod file
 674          if(fCountEvents>=fNevents){
 675             if(!fAODPathArray){
 676                AliDebug(AliLog::kDebug, "No list of AODs, keep current AOD.");
 677             }
 678             else {
 679                AliDebug(AliLog::kDebug, "Select new AOD file ...");
 680
 681                fFileId = OpenAODfile();
 682                if(fFileId<0) {
 683                   PostData(1, fHistList);
 684                   return;
 685                }
 686                fh1AODfile->Fill(fFileId);
 687                if(fAODEntries!=fNevents){
 688                   AliError("File with list of AODs and file with nb. of entries inconsistent");
 689                   PostData(1, fHistList);
 690                   return;
 691                }
 692             }
 693          }
 694
 695
 696
 697          // get next event
 698          fAODtree->GetEvent(fAODEntry);
 699
 700          // get pt hard
 701          if(mcHeader){
 702             fPtHard = mcHeader->GetPtHard();
 703             GetPtHard(kTRUE,fPtHard); // make value available for other tasks (e.g. jet response task)
 704             AliDebug(AliLog::kDebug,Form("pt hard %.2f", fPtHard));
 705          }
 706          else{
 707             AliDebug(AliLog::kDebug,"No mcHeader");
 708          }
 709          fPtHardBin = GetPtHardBin(fPtHard);
 710
 711          //Printf("pt hard (bin) %.2f (%d), xSection %.2e, trials %.2f",fPtHard, fPtHardBin, fXsection, fAvgTrials);
 712
 713          // fill event variables for each event
 714          fh1Xsec->Fill(fPtHardBin,fXsection);
 715          fh2PtHard->Fill(fPtHardBin,fPtHard);
 716
 717          fh1Trials->Fill(fPtHardBin, fAvgTrials);
 718          fh1TrialsEvtSel->Fill(fPtHardBin);
 719          fh2PtHardTrials->Fill(fPtHardBin,fPtHard,fAvgTrials);
 720
 721          // jet pt selection
 722          if(fEvtSelecMode==kEventsJetPt){
 723             Int_t nJets = fAODJets->GetEntries();
 724             //Printf("nb. jets: %d",nJets);
 725             for(Int_t ij=0; ij<nJets; ++ij){
 726                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
 727                if(!jet) continue;
 728
 729                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 730                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 731                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 732                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
 733                   useEntry = kTRUE;
 734                   break;
 735                }
 736             }
 737          }
 738
 739          // no selection
 740          if(fEvtSelecMode==kEventsAll){
 741             useEntry = kTRUE;
 742          }
 743
 744       }
 745       AliDebug(AliLog::kDebug,Form("Use entry %d from extra AOD.", fAODEntry));
 746
 747       fh2PtHardEvtSel->Fill(fPtHardBin,fPtHard);
 748       fh2AODevent->Fill(fFileId,fAODEntry);
 749
 750       TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
 751       TClonesArray *mcpartOUT = 0x0;
 752       if(mcpartIN){
 753          mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
 754          if(!mcpartOUT){
 755             AliError("Extra MC particles branch not found in output.");
 756             PostData(1, fHistList);
 757             return;
 758          }
 759          mcpartOUT->Delete();
 760       } else {
 761          AliInfo("No extra MC particles found.");
 762       }
 763
 764
 765       if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
 766          // loop over input tracks
 767          // add to output aod
 768          Int_t nTracks = 0;
 769          Int_t nJets = fAODJets->GetEntries();
 770          Int_t nSelectedJets = 0;
 771          AliAODJet *leadJet = 0x0; // used for jet tracks only
 772
 773          if(fEmbedMode==kAODFull){
 774             nTracks = fAODevent->GetNumberOfTracks();
 775
 776             for(Int_t ij=0; ij<nJets; ++ij){
 777                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
 778                if(!jet) continue;
 779                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 780                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 781                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 782                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 783
 784                   fh1JetPt->Fill(jet->Pt());
 785                   fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
 786                   nSelectedJets++;
 787
 788                }
 789             }
 790             fh1JetN->Fill(nSelectedJets);
 791          }
 792
 793          if(fEmbedMode==kAODJetTracks){
 794             // look for leading jet within selection
 795             // get reference tracks for this jet
 796             Float_t leadJetPt = 0.;
 797             for(Int_t ij=0; ij<nJets; ++ij){
 798                AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
 799                if(!jet) continue;
 800                if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 801                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 802                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 803                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 804                   if(jet->Pt()>leadJetPt){
 805                      leadJetPt = jet->Pt();
 806                      leadJet = jet;
 807                   }
 808                }
 809             }
 810             if(leadJet){
 811                nTracks = leadJet->GetRefTracks()->GetEntriesFast();
 812                fh1JetPt->Fill(leadJet->Pt());
 813                fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
 814                fh1JetN->Fill(1);
 815             }
 816          }
 817
 818          fh1TrackN->Fill((Float_t)nTracks);
 819
 820          for(Int_t it=0; it<nTracks; ++it){
 821             AliAODTrack *tmpTr = 0x0;
 822             if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
 823             if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
 824             if(!tmpTr) continue;
 825
 826             tmpTr->SetStatus(AliESDtrack::kEmbedded);
 827
 828             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 829             dummy = (*tracks)[nAODtracks-1];
 830
 831             if(fTrackFilterMap<=0 || tmpTr->TestFilterBit(fTrackFilterMap)){
 832                fh1TrackPt->Fill(tmpTr->Pt());
 833                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 834             }
 835          }
 836
 837          if(mcpartIN){
 838             Int_t nMCpart = mcpartIN->GetEntriesFast();
 839
 840             Int_t nPhysicalPrimary=0;
 841             Int_t nAODmcpart=0;
 842             for(Int_t ip=0; ip<nMCpart; ++ip){
 843                AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
 844
 845                if(fEmbedMode==kAODJetTracks){
 846                   // jet track? else continue
 847                   // not implemented yet
 848                   continue;
 849                }
 850
 851                new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
 852                dummy = (*mcpartOUT)[nAODmcpart-1];
 853
 854                if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
 855                   fh1MCTrackPt->Fill(tmpPart->Pt());
 856                   fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
 857                   nPhysicalPrimary++;
 858                }
 859             }
 860             fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
 861
 862          }
 863       } // end: embed all tracks or jet tracks
 864
 865
 866       if(fEmbedMode==kAODJet4Mom){
 867
 868          // loop over jets
 869          Int_t nJets = fAODJets->GetEntries();
 870          fh1TrackN->Fill((Float_t)nJets);
 871          for(Int_t ij=0; ij<nJets; ++ij){
 872             AliAODJet *jet = dynamic_cast<AliAODJet*>(fAODJets->At(ij));
 873             if(!jet) continue;
 874             AliAODTrack *tmpTr = (AliAODTrack*)jet;
 875             tmpTr->SetFlags(AliESDtrack::kEmbedded);
 876
 877             new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 878             dummy = (*tracks)[nAODtracks-1];
 879
 880             fh1TrackPt->Fill(tmpTr->Pt());
 881             fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 882          }
 883
 884       } // end: embed jets as 4-momenta
 885
 886
 887    } //end: embed mode with AOD
 888
 889
 890    // === embed mode with toy ===
 891    if(fEmbedMode==kToyTracks){
 892       Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
 893
 894       fh1TrackN->Fill((Float_t)nT);
 895
 896       for(Int_t i=0; i<nT; ++i){
 897
 898          Double_t pt = -1.;
 899          if(fToyMinTrackPt!=fToyMaxTrackPt){
 900             if(fToyDistributionTrackPt==0){
 901                pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
 902             } else {
 903                while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
 904                   pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
 905                   pt += fToyMinTrackPt;
 906                }
 907             }
 908          } else {
 909             pt = fToyMinTrackPt;
 910          }
 911          Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
 912          Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
 913          phi = TVector2::Phi_0_2pi(phi);
 914
 915          if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
 916
 917          Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
 918          Float_t mom[3];
 919          mom[0] = pt;
 920          mom[1] = phi;
 921          mom[2] = theta;
 922
 923          Float_t xyz[3];
 924          xyz[0] = -999.;
 925          xyz[1] = -999.;
 926          xyz[2] = -999.;
 927
 928          AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
 929          -999,   // label
 930          mom,    // momentum[3]
 931          kFALSE, // cartesian
 932          xyz,    // position
 933          kFALSE, // DCA
 934          NULL,   // covMatrix,
 935          -99,    // charge
 936          0,      // itsClusMap
 937          NULL,   // pid
 938          NULL,   // prodVertex
 939          kFALSE, // used for vertex fit
 940          kFALSE, // used for prim vtx fit
 941          AliAODTrack::kUndef, // type
 942          fToyFilterMap,  // select info
 943          -999.    // chi2 per NDF
 944          );
 945          tmpTr->SetFlags(AliESDtrack::kEmbedded);
 946
 947          new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 948          dummy = (*tracks)[nAODtracks-1];
 949
 950          fh1TrackPt->Fill(pt);
 951          fh2TrackEtaPhi->Fill(eta,phi);
 952
 953          delete tmpTr;
 954       }
 955    } //end: embed mode with toy
 956
 957
 958    PostData(1, fHistList);
 959 }
 960
 961
 962 //__________________________________________________________________________
 963 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
 964 {
 965    if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
 966
 967    if(fAODfile && fAODfile->IsOpen())
 968    fAODfile->Close();
 969
 970 }
 971
 972 //__________________________________________________________________________
 973 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
 974 {
 975    Int_t id=-1;
 976
 977    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
 978    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
 979
 980    if(env && strlen(env)){
 981       id= atoi(env);
 982       AliInfo(Form("Job index %d", id));
 983    }
 984    else{
 985       AliInfo("Job index not found. Okay if running locally.");
 986    }
 987
 988    return id;
 989 }
 990
 991 //__________________________________________________________________________
 992
 993 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
 994
 995    Int_t nFiles = fAODPathArray->GetEntries();
 996
 997    // choose randomly file and event
 998    Int_t n = -1;
 999    Float_t tmpProp = -1.;
1000    Int_t pendingEvents = fInputEntries-Entry();
1001    //Printf("input entries %d, entry %d, pending events %d", fInputEntries, (Int_t)Entry(), pendingEvents);
1002    while(rndm->Rndm()>tmpProp){
1003       n = rndm->Integer(nFiles);
1004       fAODEntries = fAODEntriesArray->At(n);
1005       Int_t tmpEntries = fAODEntries<pendingEvents ? pendingEvents : fAODEntries;
1006       tmpProp = fAODEntriesMax ? (Float_t)tmpEntries/fAODEntriesMax : 1.;
1007    }
1008    fAODEntry = (Int_t)(rndm->Uniform(fAODEntries));
1009
1010    AliInfo(Form("Select AOD file %d", n));
1011    TObjString *objStr = (TObjString*) fAODPathArray->At(n);
1012    if(!objStr){
1013       AliError("Could not get path of aod file.");
1014       return -1;
1015    }
1016    fAODPath = objStr->GetString();
1017
1018    return n;
1019 }
1020
1021 //__________________________________________________________________________
1022
1023 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
1024
1025    if(fAODPathArray) fFileId = SelectAODfile();
1026    if(fFileId<0) return -1;
1027
1028    TDirectory *owd = gDirectory;
1029    if (fAODfile) fAODfile->Close();
1030    fAODfile = TFile::Open(fAODPath.Data(),"TIMEOUT=180");
1031    owd->cd();
1032    if(!fAODfile){
1033
1034       fFileId = -1;
1035       if(fAODPathArray){
1036          if(trial<=3){
1037             AliError(Form("Could not open AOD file %s (trial %d), try again ...", fAODPath.Data(), trial));
1038             fFileId = OpenAODfile(trial+1);
1039          } else {
1040             AliError(Form("Could not open AOD file %s, give up ...",fAODPath.Data()));
1041             return -1;
1042          }
1043       } else {
1044          AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
1045       }
1046
1047       return fFileId;
1048    }
1049
1050    fAODtree = (TTree*)fAODfile->Get("aodTree");
1051
1052    if(!fAODtree){
1053       AliError("AOD tree not found.");
1054       return -1;
1055    }
1056
1057    /*
1058       fAODtree->SetBranchStatus("*",0);
1059       fAODtree->SetBranchStatus("header",1);
1060       fAODtree->SetBranchStatus("tracks*",1);
1061       fAODtree->SetBranchStatus("jets*",1);
1062       fAODtree->SetBranchStatus("clusters*",1);
1063       fAODtree->SetBranchStatus("mcparticles*",1);
1064       fAODtree->SetBranchStatus("mcHeader*",1);
1065    */
1066
1067    delete fAODevent;
1068    fAODevent = new AliAODEvent();
1069    fAODevent->ReadFromTree(fAODtree);
1070
1071
1072    // fetch header, jets, etc. from new file
1073    fNevents = fAODtree->GetEntries();
1074    mcHeader = (AliAODMCHeader*)fAODevent->FindListObject(AliAODMCHeader::StdBranchName());
1075
1076    if(fJetBranch.Length()) fAODJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
1077    else                    fAODJets = fAODevent->GetJets();
1078    if(!fAODJets){
1079       AliError("Could not find jets in AOD. Check jet branch when indicated.");
1080       return -1;
1081    }
1082
1083    TFile *curfile = fAODtree->GetCurrentFile();
1084    if (!curfile) {
1085       AliError("No current file.");
1086       return -1;
1087    }
1088
1089    Float_t trials = 1.;
1090    Float_t xsec = 0.;
1091    PythiaInfoFromFile(curfile->GetName(),xsec,trials);
1092    fXsection = xsec;
1093
1094    // construct a poor man average trials
1095    Float_t nEntries = (Float_t)fAODtree->GetTree()->GetEntries();
1096    if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
1097
1098    if(fFileId>=0){
1099       AliInfo(Form("Read successfully AOD event from file %d",fFileId));
1100    }
1101
1102    fCountEvents=0; // new file, reset counter
1103
1104    return fFileId;  // file position in AOD path array, if array available
1105 }
1106
1107
1108 //____________________________________________________________________________
1109 Float_t AliAnalysisTaskFastEmbedding::GetPtHard(Bool_t bSet, Float_t newValue){
1110
1111    static Float_t ptHard = -1.;
1112    if(bSet) ptHard = newValue;
1113
1114    return ptHard;
1115 }
1116
1117 //____________________________________________________________________________
1118 Int_t AliAnalysisTaskFastEmbedding::GetPtHardBin(Double_t ptHard){
1119
1120    const Int_t nBins = 10;
1121    Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1122
1123    Int_t bin = -1;
1124    while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1125       bin++;
1126    }
1127
1128    return bin;
1129 }
1130
1131 //____________________________________________________________________________
1132 Bool_t AliAnalysisTaskFastEmbedding::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1133    //
1134    // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1135    // This should provide the path to the AOD/ESD file
1136    // (taken from PWG4/JetTasks/AliAnalysisHelperJetTasks)
1137
1138    TString file(currFile);
1139    fXsec = 0;
1140    fTrials = 1;
1141
1142    if(file.Contains("root_archive.zip#")){
1143       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1144       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1145       file.Replace(pos+1,20,"");
1146    }
1147    else {
1148       // not an archive take the basename....
1149       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1150    }
1151
1152    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
1153    if(!fxsec){
1154       // next trial fetch the histgram file
1155       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1156       if(!fxsec){
1157          // not a severe condition but inciate that we have no information
1158          AliDebug(AliLog::kDebug,Form("Neither pyxsec.root nor pyxsec_hists.root found in %s",file.Data()));
1159          return kFALSE;
1160       }
1161       else{
1162          // find the tlist we want to be independtent of the name so use the Tkey
1163          TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1164          if(!key){
1165             fxsec->Close();
1166             return kFALSE;
1167          }
1168          TList *list = dynamic_cast<TList*>(key->ReadObj());
1169          if(!list){
1170             fxsec->Close();
1171             return kFALSE;
1172          }
1173          fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1174          fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1175          fxsec->Close();
1176       }
1177    } // no tree pyxsec.root
1178    else {
1179       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1180       if(!xtree){
1181          fxsec->Close();
1182          return kFALSE;
1183       }
1184       UInt_t   ntrials  = 0;
1185       Double_t  xsection  = 0;
1186       xtree->SetBranchAddress("xsection",&xsection);
1187       xtree->SetBranchAddress("ntrials",&ntrials);
1188       xtree->GetEntry(0);
1189       fTrials = ntrials;
1190       fXsec = xsection;
1191       fxsec->Close();
1192    }
1193    return kTRUE;
1194 }
1195
1196