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
  37
  38 #include "AliAnalysisTaskFastEmbedding.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDtrack.h"
  41 #include "AliAODEvent.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODJet.h"
  44 #include "AliAODMCParticle.h"
  45
  46 #include "AliLog.h"
  47
  48 ClassImp(AliAnalysisTaskFastEmbedding)
  49
  50 //__________________________________________________________________________
  51 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding()
  52     : AliAnalysisTaskSE()
  53       ,fAODout(0)
  54       ,fAODevent(0)
  55       ,fAODtree(0)
  56       ,fAODfile(0)
  57       ,rndm(0)
  58       ,fAODPathArray(0)
  59       ,fAODPath("AliAOD.root")
  60       ,fTrackBranch("aodExtraTracks")
  61       ,fMCparticlesBranch("aodExtraMCparticles")
  62       ,fJetBranch("")
  63       ,fEntry(0)
  64       ,fEmbedMode(0)
  65       ,fEvtSelecMode(0)
  66       ,fEvtSelMinJetPt(-1)
  67       ,fEvtSelMaxJetPt(-1)
  68       ,fEvtSelMinJetEta(-999.)
  69       ,fEvtSelMaxJetEta( 999.)
  70       ,fEvtSelMinJetPhi(0.)
  71       ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
  72       ,fToyMinNbOfTracks(1)
  73       ,fToyMaxNbOfTracks(1)
  74       ,fToyMinTrackPt(50.)
  75       ,fToyMaxTrackPt(50.)
  76       ,fToyDistributionTrackPt(0.)
  77       ,fToyMinTrackEta(-.5)
  78       ,fToyMaxTrackEta(.5)
  79       ,fToyMinTrackPhi(0.)
  80       ,fToyMaxTrackPhi(2*TMath::Pi())
  81       ,fToyFilterMap(0)
  82       ,fHistList(0)
  83       ,fh1TrackPt(0)
  84       ,fh2TrackEtaPhi(0)
  85       ,fh1TrackN(0)
  86       ,fh1JetPt(0)
  87       ,fh2JetEtaPhi(0)
  88       ,fh1JetN(0)
  89       ,fh1MCTrackPt(0)
  90       ,fh2MCTrackEtaPhi(0)
  91       ,fh1MCTrackN(0)
  92       ,fh1AODfile(0)
  93
  94
  95 {
  96     // default constructor
  97
  98 }
  99
 100
 101 //__________________________________________________________________________
 102 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const char *name)
 103     : AliAnalysisTaskSE(name)
 104       ,fAODout(0)
 105       ,fAODevent(0)
 106       ,fAODtree(0)
 107       ,fAODfile(0)
 108       ,rndm(0)
 109       ,fAODPathArray(0)
 110       ,fAODPath("AliAOD.root")
 111       ,fTrackBranch("aodExtraTracks")
 112       ,fMCparticlesBranch("aodExtraMCparticles")
 113       ,fJetBranch("")
 114       ,fEntry(0)
 115       ,fEmbedMode(0)
 116       ,fEvtSelecMode(0)
 117       ,fEvtSelMinJetPt(-1)
 118       ,fEvtSelMaxJetPt(-1)
 119       ,fEvtSelMinJetEta(-999.)
 120       ,fEvtSelMaxJetEta( 999.)
 121       ,fEvtSelMinJetPhi(0.)
 122       ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
 123       ,fToyMinNbOfTracks(1)
 124       ,fToyMaxNbOfTracks(1)
 125       ,fToyMinTrackPt(50.)
 126       ,fToyMaxTrackPt(50.)
 127       ,fToyDistributionTrackPt(0.)
 128       ,fToyMinTrackEta(-.5)
 129       ,fToyMaxTrackEta(.5)
 130       ,fToyMinTrackPhi(0.)
 131       ,fToyMaxTrackPhi(2*TMath::Pi())
 132       ,fToyFilterMap(0)
 133       ,fHistList(0)
 134       ,fh1TrackPt(0)
 135       ,fh2TrackEtaPhi(0)
 136       ,fh1TrackN(0)
 137       ,fh1JetPt(0)
 138       ,fh2JetEtaPhi(0)
 139       ,fh1JetN(0)
 140       ,fh1MCTrackPt(0)
 141       ,fh2MCTrackEtaPhi(0)
 142       ,fh1MCTrackN(0)
 143       ,fh1AODfile(0)
 144 {
 145     // constructor
 146     DefineOutput(1, TList::Class());
 147 }
 148
 149
 150 //__________________________________________________________________________
 151 AliAnalysisTaskFastEmbedding::AliAnalysisTaskFastEmbedding(const AliAnalysisTaskFastEmbedding &copy)
 152     : AliAnalysisTaskSE()
 153       ,fAODout(copy.fAODout)
 154       ,fAODevent(copy.fAODevent)
 155       ,fAODtree(copy.fAODtree)
 156       ,fAODfile(copy.fAODfile)
 157       ,rndm(copy.rndm)
 158       ,fAODPathArray(copy.fAODPathArray)
 159       ,fAODPath(copy.fAODPath)
 160       ,fTrackBranch(copy.fTrackBranch)
 161       ,fMCparticlesBranch(copy.fMCparticlesBranch)
 162       ,fJetBranch(copy.fJetBranch)
 163       ,fEntry(copy.fEntry)
 164       ,fEmbedMode(copy.fEmbedMode)
 165       ,fEvtSelecMode(copy.fEvtSelecMode)
 166       ,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
 167       ,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
 168       ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
 169       ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
 170       ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
 171       ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
 172       ,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
 173       ,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
 174       ,fToyMinTrackPt(copy.fToyMinTrackPt)
 175       ,fToyMaxTrackPt(copy.fToyMaxTrackPt)
 176       ,fToyDistributionTrackPt(copy.fToyDistributionTrackPt)
 177       ,fToyMinTrackEta(copy.fToyMinTrackEta)
 178       ,fToyMaxTrackEta(copy.fToyMaxTrackEta)
 179       ,fToyMinTrackPhi(copy.fToyMinTrackPhi)
 180       ,fToyMaxTrackPhi(copy.fToyMaxTrackPhi)
 181       ,fToyFilterMap(copy.fToyFilterMap)
 182       ,fHistList(copy.fHistList)
 183       ,fh1TrackPt(copy.fh1TrackPt)
 184       ,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
 185       ,fh1TrackN(copy.fh1TrackN)
 186       ,fh1JetPt(copy.fh1JetPt)
 187       ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
 188       ,fh1JetN(copy.fh1JetN)
 189       ,fh1MCTrackPt(copy.fh1MCTrackPt)
 190       ,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
 191       ,fh1MCTrackN(copy.fh1MCTrackN)
 192       ,fh1AODfile(copy.fh1AODfile)
 193 {
 194     // copy constructor
 195 }
 196
 197
 198 //__________________________________________________________________________
 199 AliAnalysisTaskFastEmbedding& AliAnalysisTaskFastEmbedding::operator=(const AliAnalysisTaskFastEmbedding& o)
 200 {
 201     // assignment
 202
 203     if(this!=&o){
 204         AliAnalysisTaskSE::operator=(o);
 205         fAODout            = o.fAODout;
 206         fAODevent          = o.fAODevent;
 207         fAODtree           = o.fAODtree;
 208         fAODfile           = o.fAODfile;
 209         rndm               = o.rndm;
 210         fAODPathArray       = o.fAODPathArray;
 211         fAODPath           = o.fAODPath;
 212         fTrackBranch       = o.fTrackBranch;
 213         fMCparticlesBranch = o.fMCparticlesBranch;
 214         fJetBranch         = o.fJetBranch;
 215         fEntry             = o.fEntry;
 216         fEmbedMode         = o.fEmbedMode;
 217         fEvtSelecMode      = o.fEvtSelecMode;
 218         fEvtSelMinJetPt    = o.fEvtSelMinJetPt;
 219         fEvtSelMaxJetPt    = o.fEvtSelMaxJetPt;
 220         fEvtSelMinJetEta   = o.fEvtSelMinJetEta;
 221         fEvtSelMaxJetEta   = o.fEvtSelMaxJetEta;
 222         fEvtSelMinJetPhi   = o.fEvtSelMinJetPhi;
 223         fEvtSelMaxJetPhi   = o.fEvtSelMaxJetPhi;
 224         fToyMinNbOfTracks  = o.fToyMinNbOfTracks;
 225         fToyMaxNbOfTracks  = o.fToyMaxNbOfTracks;
 226         fToyMinTrackPt     = o.fToyMinTrackPt;
 227         fToyMaxTrackPt     = o.fToyMaxTrackPt;
 228         fToyDistributionTrackPt = o.fToyDistributionTrackPt;
 229         fToyMinTrackEta    = o.fToyMinTrackEta;
 230         fToyMaxTrackEta    = o.fToyMaxTrackEta;
 231         fToyMinTrackPhi    = o.fToyMinTrackPhi;
 232         fToyMaxTrackPhi    = o.fToyMaxTrackPhi;
 233         fToyFilterMap      = o.fToyFilterMap;
 234         fHistList          = o.fHistList;
 235         fh1TrackPt         = o.fh1TrackPt;
 236         fh2TrackEtaPhi     = o.fh2TrackEtaPhi;
 237         fh1TrackN          = o.fh1TrackN;
 238         fh1JetPt           = o.fh1JetPt;
 239         fh2JetEtaPhi       = o.fh2JetEtaPhi;
 240         fh1JetN            = o.fh1JetN;
 241         fh1MCTrackPt       = o.fh1MCTrackPt;
 242         fh2MCTrackEtaPhi   = o.fh2MCTrackEtaPhi;
 243         fh1MCTrackN        = o.fh1MCTrackN;
 244         fh1AODfile         = o.fh1AODfile;
 245     }
 246
 247     return *this;
 248 }
 249
 250
 251 //__________________________________________________________________________
 252 AliAnalysisTaskFastEmbedding::~AliAnalysisTaskFastEmbedding()
 253 {
 254     // destructor
 255     delete rndm;
 256 }
 257
 258
 259 //__________________________________________________________________________
 260 void AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()
 261 {
 262     // create output objects
 263     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
 264
 265         OpenFile(1);
 266     if(!fHistList) fHistList = new TList();
 267     fHistList->SetOwner(kTRUE);
 268
 269
 270     // set seed
 271     rndm = new TRandom3();
 272     Int_t id = GetJobID();
 273     if(id>-1) rndm->SetSeed(id);
 274     else      rndm->SetSeed();   // a TTUID is generated and used for seed
 275     AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
 276
 277     // embed mode with AOD
 278     Int_t rc = -1;
 279     if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 280
 281        // open input AOD
 282        rc = OpenAODfile();
 283        if(rc<0){
 284            PostData(1, fHistList);
 285            return;
 286            }
 287     } //end: embed mode with AOD
 288
 289
 290     // connect output aod
 291     // create a new branch for extra tracks
 292     fAODout = AODEvent();
 293     if(!fAODout){
 294         AliError("Output AOD not found.");
 295                 PostData(1, fHistList);
 296         return;
 297     }
 298     if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
 299         AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
 300         TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
 301         tracks->SetName(fTrackBranch.Data());
 302         AddAODBranch("TClonesArray", &tracks);
 303     }
 304     // create new branch for extra mcparticle if available as input
 305     if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
 306        AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
 307        TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
 308        mcparticles->SetName(fMCparticlesBranch.Data());
 309        AddAODBranch("TClonesArray", &mcparticles);
 310     }
 311
 312
 313
 314     //qa histograms
 315     Bool_t oldStatus = TH1::AddDirectoryStatus();
 316     TH1::AddDirectory(kFALSE);
 317
 318     fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
 319     fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 320     fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
 321
 322     fHistList->Add(fh1TrackPt);
 323     fHistList->Add(fh2TrackEtaPhi);
 324     fHistList->Add(fh1TrackN);
 325
 326         if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 327
 328             fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 250, 0., 250.);
 329             fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 330             fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
 331
 332             fHistList->Add(fh1JetPt);
 333         fHistList->Add(fh2JetEtaPhi);
 334         fHistList->Add(fh1JetN);
 335     }
 336
 337
 338     if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
 339
 340        fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
 341        fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 342        fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
 343
 344        fHistList->Add(fh1MCTrackPt);
 345        fHistList->Add(fh2MCTrackEtaPhi);
 346        fHistList->Add(fh1MCTrackN);
 347
 348     }
 349
 350     fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 100, 0, 100);
 351     fHistList->Add(fh1AODfile);
 352
 353     // =========== Switch on Sumw2 for all histos ===========
 354     for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 355         TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 356         if (h1){
 357           h1->Sumw2();
 358           continue;
 359         }
 360     }
 361
 362     TH1::AddDirectory(oldStatus);
 363
 364
 365     fh1AODfile->Fill(rc);
 366
 367     PostData(1, fHistList);
 368 }
 369
 370
 371 //__________________________________________________________________________
 372 void AliAnalysisTaskFastEmbedding::Init()
 373 {
 374     // Initialization
 375     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
 376
 377 }
 378
 379
 380 //__________________________________________________________________________
 381 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
 382 {
 383     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
 384
 385     if(!fAODout){
 386        AliError("Need output AOD, but is not connected.");
 387        PostData(1, fHistList);
 388        return;
 389     }
 390
 391     // connect aod out
 392     TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
 393     if(!tracks){
 394         AliError("Extra track branch not found in output.");
 395         PostData(1, fHistList);
 396         return;
 397     }
 398     tracks->Delete();
 399     Int_t nAODtracks=0;
 400
 401
 402     TRef dummy;
 403
 404     // === embed mode with AOD ===
 405     if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 406        if(!fAODevent){
 407           AliError("Need input AOD, but is not connected.");
 408           PostData(1, fHistList);
 409           return;
 410        }
 411
 412        // fetch jets
 413        TClonesArray *aodJets = 0;
 414        if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
 415        else                    aodJets = fAODevent->GetJets();
 416        if(!aodJets){
 417           AliError("Could not find jets in AOD. Check jet branch when indicated.");
 418           PostData(1, fHistList);
 419           return;
 420        }
 421
 422
 423        Int_t nEvents = fAODtree->GetEntries();
 424
 425        Bool_t useEntry = kFALSE;
 426        while(!useEntry){  // protection need, if no event fulfills requierment
 427           if(fEntry>=nEvents){
 428               fEntry=0;
 429               if(!fAODPathArray){
 430                  AliDebug(AliLog::kDebug, "Last event in AOD reached, start from entry 0 again.");
 431               }
 432               else {
 433                  AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
 434
 435                  Int_t rc = OpenAODfile();
 436                  if(rc<0) {
 437                     PostData(1, fHistList);
 438                     return;
 439                  }
 440                  fh1AODfile->Fill(rc);
 441
 442                  // new file => we must use the new jet array
 443                  if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
 444                  else                    aodJets = fAODevent->GetJets();
 445                  if(!aodJets){
 446                    AliError("Could not find jets in AOD. Check jet branch when indicated.");
 447                    PostData(1, fHistList);
 448                    return;
 449                  }
 450               }
 451           }
 452
 453           fAODtree->GetEvent(fEntry);
 454
 455           // jet pt selection
 456           if(fEvtSelecMode==kEventsJetPt){
 457              Int_t nJets = aodJets->GetEntries();
 458              for(Int_t ij=0; ij<nJets; ++ij){
 459                  AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 460                  if(!jet) continue;
 461
 462                  if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 463                     && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 464                     && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 465                     && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
 466                     useEntry = kTRUE;
 467                     break;
 468                  }
 469              }
 470           }
 471
 472           // no selection
 473           if(fEvtSelecMode==kEventsAll){
 474              useEntry = kTRUE;
 475           }
 476
 477           fEntry++;
 478        }
 479        AliInfo(Form("Use entry %d from extra AOD.", fEntry-1));
 480
 481
 482        TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
 483        TClonesArray *mcpartOUT = 0x0;
 484        if(mcpartIN){
 485            mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
 486        if(!mcpartOUT){
 487          AliError("Extra MC particles branch not found in output.");
 488          PostData(1, fHistList);
 489          return;
 490        }
 491            mcpartOUT->Delete();
 492        } else {
 493            AliInfo("No extra MC particles found.");
 494        }
 495
 496
 497        if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
 498            // loop over input tracks
 499            // add to output aod
 500            Int_t nTracks = 0;
 501            Int_t nJets = aodJets->GetEntries();
 502            Int_t nSelectedJets = 0;
 503        AliAODJet *leadJet = 0x0; // used for jet tracks only
 504
 505            if(fEmbedMode==kAODFull){
 506                            nTracks = fAODevent->GetNumberOfTracks();
 507
 508                                for(Int_t ij=0; ij<nJets; ++ij){
 509                        AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 510                        if(!jet) continue;
 511                        if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 512                           && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 513                           && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 514                           && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 515
 516                                                       fh1JetPt->Fill(jet->Pt());
 517                                                           fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
 518                                                           nSelectedJets++;
 519
 520                                            }
 521                    }
 522                                    fh1JetN->Fill(nSelectedJets);
 523                    }
 524
 525            if(fEmbedMode==kAODJetTracks){
 526               // look for leading jet within selection
 527               // get reference tracks for this jet
 528               Float_t leadJetPt = 0.;
 529               for(Int_t ij=0; ij<nJets; ++ij){
 530                   AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 531                   if(!jet) continue;
 532                   if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 533                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 534                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 535                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 536                      if(jet->Pt()>leadJetPt){
 537                          leadJetPt = jet->Pt();
 538                          leadJet = jet;
 539                      }
 540                   }
 541                }
 542                if(leadJet){
 543                            nTracks = leadJet->GetRefTracks()->GetEntriesFast();
 544                                    fh1JetPt->Fill(leadJet->Pt());
 545                    fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
 546                    fh1JetN->Fill(1);
 547                            }
 548            }
 549
 550            fh1TrackN->Fill((Float_t)nTracks);
 551
 552            for(Int_t it=0; it<nTracks; ++it){
 553                AliAODTrack *tmpTr = 0x0;
 554                if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
 555                if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
 556                if(!tmpTr) continue;
 557
 558                tmpTr->SetStatus(AliESDtrack::kEmbedded);
 559
 560                new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 561                dummy = (*tracks)[nAODtracks-1];
 562
 563                fh1TrackPt->Fill(tmpTr->Pt());
 564                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 565            }
 566
 567            if(mcpartIN){
 568                Int_t nMCpart = mcpartIN->GetEntriesFast();
 569
 570                Int_t nPhysicalPrimary=0;
 571                Int_t nAODmcpart=0;
 572                for(Int_t ip=0; ip<nMCpart; ++ip){
 573                    AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
 574
 575                    if(fEmbedMode==kAODJetTracks){
 576                       // jet track? else continue
 577                       // not implemented yet
 578                       continue;
 579                    }
 580
 581                    new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
 582                    dummy = (*mcpartOUT)[nAODmcpart-1];
 583
 584                    if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
 585                       fh1MCTrackPt->Fill(tmpPart->Pt());
 586                       fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
 587                       nPhysicalPrimary++;
 588                    }
 589                }
 590                fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
 591
 592            }
 593        } // end: embed all tracks or jet tracks
 594
 595
 596        if(fEmbedMode==kAODJet4Mom){
 597
 598            // loop over jets
 599            Int_t nJets = aodJets->GetEntries();
 600            fh1TrackN->Fill((Float_t)nJets);
 601            for(Int_t ij=0; ij<nJets; ++ij){
 602                AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 603                if(!jet) continue;
 604                AliAODTrack *tmpTr = (AliAODTrack*)jet;
 605                tmpTr->SetFlags(AliESDtrack::kEmbedded);
 606
 607                new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 608                dummy = (*tracks)[nAODtracks-1];
 609
 610                fh1TrackPt->Fill(tmpTr->Pt());
 611                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 612            }
 613
 614        } // end: embed jets as 4-momenta
 615
 616
 617     } //end: embed mode with AOD
 618
 619
 620     // === embed mode with toy ===
 621     if(fEmbedMode==kToyTracks){
 622         Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
 623
 624         fh1TrackN->Fill((Float_t)nT);
 625
 626         for(Int_t i=0; i<nT; ++i){
 627
 628            Double_t pt = -1.;
 629            if(fToyMinTrackPt!=fToyMaxTrackPt){
 630               if(fToyDistributionTrackPt==0){
 631                  pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
 632               } else {
 633                  while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
 634                     pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
 635                     pt += fToyMinTrackPt;
 636                  }
 637               }
 638            } else {
 639               pt = fToyMinTrackPt;
 640            }
 641            Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
 642            Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
 643            phi = TVector2::Phi_0_2pi(phi);
 644
 645            if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
 646
 647            Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
 648            Float_t mom[3];
 649            mom[0] = pt;
 650            mom[1] = phi;
 651            mom[2] = theta;
 652
 653            Float_t xyz[3];
 654            xyz[0] = -999.;
 655            xyz[1] = -999.;
 656            xyz[2] = -999.;
 657
 658            AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
 659                                               -999,   // label
 660                                               mom,    // momentum[3]
 661                                               kFALSE, // cartesian
 662                                               xyz,    // position
 663                                               kFALSE, // DCA
 664                                               NULL,   // covMatrix,
 665                                               -99,    // charge
 666                                               0,      // itsClusMap
 667                                               NULL,   // pid
 668                                               NULL,   // prodVertex
 669                                               kFALSE, // used for vertex fit
 670                                               kFALSE, // used for prim vtx fit
 671                                               AliAODTrack::kUndef, // type
 672                                               fToyFilterMap,  // select info
 673                                               -999.    // chi2 per NDF
 674                                             );
 675            tmpTr->SetFlags(AliESDtrack::kEmbedded);
 676
 677            new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 678            dummy = (*tracks)[nAODtracks-1];
 679
 680            fh1TrackPt->Fill(pt);
 681            fh2TrackEtaPhi->Fill(eta,phi);
 682
 683            delete tmpTr;
 684         }
 685     } //end: embed mode with toy
 686
 687
 688     PostData(1, fHistList);
 689 }
 690
 691
 692 //__________________________________________________________________________
 693 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
 694 {
 695     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
 696
 697     if(fAODfile && fAODfile->IsOpen())
 698     fAODfile->Close();
 699
 700 }
 701
 702 //__________________________________________________________________________
 703 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
 704 {
 705    Int_t id=-1;
 706
 707    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
 708    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
 709
 710    if(env && strlen(env)){
 711        id= atoi(env);
 712        AliInfo(Form("Job index %d", id));
 713    }
 714    else{
 715        AliInfo("Job index not found. Okay if running locally.");
 716    }
 717
 718    return id;
 719 }
 720
 721 //__________________________________________________________________________
 722
 723 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
 724
 725      Int_t nFiles = fAODPathArray->GetEntries();
 726      Int_t n = rndm->Integer(nFiles);
 727      AliInfo(Form("Select AOD file %d", n));
 728      TObjString *objStr = (TObjString*) fAODPathArray->At(n);
 729      if(!objStr){
 730           AliError("Could not get path of aod file.");
 731           return -1;
 732      }
 733      fAODPath = objStr->GetString();
 734
 735      return n;
 736 }
 737
 738 //__________________________________________________________________________
 739
 740 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
 741
 742     if(trial>3){
 743         AliError("Could not open AOD files, give up ...");
 744         return -1;
 745     }
 746
 747         Int_t rc = 0;
 748         if(fAODPathArray) rc = SelectAODfile();
 749         if(rc<0) return -1;
 750
 751     TDirectory *owd = gDirectory;
 752     if (fAODfile)
 753       fAODfile->Close();
 754     fAODfile = TFile::Open(fAODPath.Data());
 755     owd->cd();
 756     if(!fAODfile){
 757
 758        rc = -1;
 759        if(fAODPathArray){
 760            AliError(Form("Could not open AOD file %s, try another from the list ...", fAODPath.Data()));
 761            rc = OpenAODfile(trial+1);
 762        } else {
 763                AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
 764            }
 765
 766        return rc;
 767     }
 768
 769     fAODtree = (TTree*)fAODfile->Get("aodTree");
 770
 771     if(!fAODtree){
 772        AliError("AOD tree not found.");
 773        return -1;
 774     }
 775
 776     delete fAODevent;
 777     fAODevent = new AliAODEvent();
 778     fAODevent->ReadFromTree(fAODtree);
 779     return rc;  // file position in AOD path array, if array available
 780 }
 781