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