JETAN/AliAnalysisTaskFastEmbedding.cxx

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