JETAN/AliAnalysisTaskFastEmbedding.cxx

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