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