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     // connect output aod
 266     // create a new branch for extra tracks
 267     fAODout = AODEvent();
 268     if(!fAODout){
 269         AliError("Output AOD not found.");
 270         return;
 271     }
 272     if(!fAODout->FindListObject(fTrackBranch.Data()) && strlen(fTrackBranch.Data())){
 273         AliInfo(Form("Add AOD branch %s", fTrackBranch.Data()));
 274         TClonesArray *tracks = new TClonesArray("AliAODTrack",0);
 275         tracks->SetName(fTrackBranch.Data());
 276         AddAODBranch("TClonesArray", &tracks);
 277     }
 278     // create new branch for extra mcparticle if available as input
 279     if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
 280        AliInfo(Form("Add AOD branch %s", fMCparticlesBranch.Data()));
 281        TClonesArray *mcparticles = new TClonesArray("AliAODMCParticle",0);
 282        mcparticles->SetName(fMCparticlesBranch.Data());
 283        AddAODBranch("TClonesArray", &mcparticles);
 284     }
 285
 286
 287
 288
 289     //qa histograms
 290
 291     OpenFile(1);
 292     if(!fHistList) fHistList = new TList();
 293     fHistList->SetOwner(kTRUE);
 294
 295     Bool_t oldStatus = TH1::AddDirectoryStatus();
 296     TH1::AddDirectory(kFALSE);
 297
 298     fh1TrackPt      =  new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
 299     fh2TrackEtaPhi  =  new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 300     fh1TrackN       =  new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
 301
 302     fHistList->Add(fh1TrackPt);
 303     fHistList->Add(fh2TrackEtaPhi);
 304     fHistList->Add(fh1TrackN);
 305
 306         if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 307
 308             fh1JetPt        =  new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 250, 0., 250.);
 309             fh2JetEtaPhi    =  new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 310             fh1JetN         =  new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
 311
 312             fHistList->Add(fh1JetPt);
 313         fHistList->Add(fh2JetEtaPhi);
 314         fHistList->Add(fh1JetN);
 315     }
 316
 317
 318     if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
 319
 320        fh1MCTrackPt      =  new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
 321        fh2MCTrackEtaPhi  =  new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
 322        fh1MCTrackN       =  new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
 323
 324        fHistList->Add(fh1MCTrackPt);
 325        fHistList->Add(fh2MCTrackEtaPhi);
 326        fHistList->Add(fh1MCTrackN);
 327
 328     }
 329
 330     fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 100, 0, 100);
 331     fHistList->Add(fh1AODfile);
 332
 333     // =========== Switch on Sumw2 for all histos ===========
 334     for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 335         TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 336         if (h1){
 337           h1->Sumw2();
 338           continue;
 339         }
 340     }
 341
 342     TH1::AddDirectory(oldStatus);
 343
 344
 345     // set seed
 346     rndm = new TRandom3();
 347     Int_t id = GetJobID();
 348     if(id>-1) rndm->SetSeed(id);
 349     else      rndm->SetSeed();   // a TTUID is generated and used for seed
 350     AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
 351
 352     // embed mode with AOD
 353     if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 354
 355        // open input AOD
 356        Int_t rc = OpenAODfile();
 357        if(rc<0) return;
 358        fh1AODfile->Fill(rc);
 359
 360     } //end: embed mode with AOD
 361
 362
 363
 364    PostData(1, fHistList);
 365 }
 366
 367
 368 //__________________________________________________________________________
 369 void AliAnalysisTaskFastEmbedding::Init()
 370 {
 371     // Initialization
 372     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
 373
 374 }
 375
 376
 377 //__________________________________________________________________________
 378 void AliAnalysisTaskFastEmbedding::UserExec(Option_t *)
 379 {
 380     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserExec()");
 381
 382     if(!fAODout){
 383        AliError("Need output AOD, but is not connected.");
 384        PostData(1, fHistList);
 385        return;
 386     }
 387
 388     // connect aod out
 389     TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
 390     if(!tracks){
 391         AliError("Extra track branch not found in output.");
 392         PostData(1, fHistList);
 393         return;
 394     }
 395     tracks->Delete();
 396     Int_t nAODtracks=0;
 397
 398
 399     TRef dummy;
 400
 401     // === embed mode with AOD ===
 402     if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
 403        if(!fAODevent){
 404           AliError("Need input AOD, but is not connected.");
 405           PostData(1, fHistList);
 406           return;
 407        }
 408
 409        // fetch jets
 410        TClonesArray *aodJets = 0;
 411        if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
 412        else                    aodJets = fAODevent->GetJets();
 413        if(!aodJets){
 414           AliError("Could not find jets in AOD. Check jet branch when indicated.");
 415           PostData(1, fHistList);
 416           return;
 417        }
 418
 419
 420        Int_t nEvents = fAODtree->GetEntries();
 421
 422        Bool_t useEntry = kFALSE;
 423        while(!useEntry){  // protection need, if no event fulfills requierment
 424           if(fEntry>=nEvents){
 425               fEntry=0;
 426               if(!fAODPathArray){
 427                  AliDebug(AliLog::kDebug, "Last event in AOD reached, start from entry 0 again.");
 428               }
 429               else {
 430                  AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
 431
 432                  Int_t rc = OpenAODfile();
 433                  if(rc<0) {
 434                     PostData(1, fHistList);
 435                     return;
 436                  }
 437                  fh1AODfile->Fill(rc);
 438
 439                  // new file => we must use the new jet array
 440                  if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
 441                  else                    aodJets = fAODevent->GetJets();
 442                  if(!aodJets){
 443                    AliError("Could not find jets in AOD. Check jet branch when indicated.");
 444                    PostData(1, fHistList);
 445                    return;
 446                  }
 447               }
 448           }
 449
 450           fAODtree->GetEvent(fEntry);
 451
 452           // jet pt selection
 453           if(fEvtSelecMode==kEventsJetPt){
 454              Int_t nJets = aodJets->GetEntries();
 455              for(Int_t ij=0; ij<nJets; ++ij){
 456                  AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 457                  if(!jet) continue;
 458
 459                  if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 460                     && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 461                     && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 462                     && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
 463                     useEntry = kTRUE;
 464                     break;
 465                  }
 466              }
 467           }
 468
 469           // no selection
 470           if(fEvtSelecMode==kEventsAll){
 471              useEntry = kTRUE;
 472           }
 473
 474           fEntry++;
 475        }
 476        AliInfo(Form("Use entry %d from extra AOD.", fEntry-1));
 477
 478
 479        TClonesArray *mcpartIN  = (TClonesArray*)(fAODevent->FindListObject("mcparticles"));
 480        TClonesArray *mcpartOUT = 0x0;
 481        if(mcpartIN){
 482            mcpartOUT = (TClonesArray*)(fAODout->FindListObject(fMCparticlesBranch.Data()));
 483            mcpartOUT->Delete();
 484        } else {
 485            AliInfo("No extra MC particles found.");
 486        }
 487
 488
 489        if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
 490            // loop over input tracks
 491            // add to output aod
 492            Int_t nTracks = 0;
 493            Int_t nJets = aodJets->GetEntries();
 494            Int_t nSelectedJets = 0;
 495        AliAODJet *leadJet = 0x0; // used for jet tracks only
 496
 497            if(fEmbedMode==kAODFull){
 498                            nTracks = fAODevent->GetNumberOfTracks();
 499
 500                                for(Int_t ij=0; ij<nJets; ++ij){
 501                        AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 502                        if(!jet) continue;
 503                        if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 504                           && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 505                           && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 506                           && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 507
 508                                                       fh1JetPt->Fill(jet->Pt());
 509                                                           fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
 510                                                           nSelectedJets++;
 511
 512                                            }
 513                    }
 514                                    fh1JetN->Fill(nSelectedJets);
 515                    }
 516
 517            if(fEmbedMode==kAODJetTracks){
 518               // look for leading jet within selection
 519               // get reference tracks for this jet
 520               Float_t leadJetPt = 0.;
 521               for(Int_t ij=0; ij<nJets; ++ij){
 522                   AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 523                   if(!jet) continue;
 524                   if(   (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
 525                      && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
 526                      && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
 527                      && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
 528                      if(jet->Pt()>leadJetPt){
 529                          leadJetPt = jet->Pt();
 530                          leadJet = jet;
 531                      }
 532                   }
 533                }
 534                if(leadJet){
 535                            nTracks = leadJet->GetRefTracks()->GetEntriesFast();
 536                                    fh1JetPt->Fill(leadJet->Pt());
 537                    fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
 538                    fh1JetN->Fill(1);
 539                            }
 540            }
 541
 542            fh1TrackN->Fill((Float_t)nTracks);
 543
 544            for(Int_t it=0; it<nTracks; ++it){
 545                AliAODTrack *tmpTr = 0x0;
 546                if(fEmbedMode==kAODFull)      tmpTr = fAODevent->GetTrack(it);
 547                if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
 548                if(!tmpTr) continue;
 549
 550                tmpTr->SetStatus(AliESDtrack::kEmbedded);
 551
 552                new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 553                dummy = (*tracks)[nAODtracks-1];
 554
 555                fh1TrackPt->Fill(tmpTr->Pt());
 556                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 557            }
 558
 559            if(mcpartIN){
 560                Int_t nMCpart = mcpartIN->GetEntriesFast();
 561
 562                Int_t nPhysicalPrimary=0;
 563                Int_t nAODmcpart=0;
 564                for(Int_t ip=0; ip<nMCpart; ++ip){
 565                    AliAODMCParticle *tmpPart = (AliAODMCParticle*) mcpartIN->At(ip);
 566
 567                    if(fEmbedMode==kAODJetTracks){
 568                       // jet track? else continue
 569                       // not implemented yet
 570                       continue;
 571                    }
 572
 573                    new((*mcpartOUT)[nAODmcpart++]) AliAODMCParticle(*tmpPart);
 574                    dummy = (*mcpartOUT)[nAODmcpart-1];
 575
 576                    if(tmpPart->IsPhysicalPrimary() && tmpPart->Charge()!=0. && tmpPart->Charge()!=-99. ){
 577                       fh1MCTrackPt->Fill(tmpPart->Pt());
 578                       fh2MCTrackEtaPhi->Fill(tmpPart->Eta(), tmpPart->Phi());
 579                       nPhysicalPrimary++;
 580                    }
 581                }
 582                fh1MCTrackN->Fill((Float_t)nPhysicalPrimary);
 583
 584            }
 585        } // end: embed all tracks or jet tracks
 586
 587
 588        if(fEmbedMode==kAODJet4Mom){
 589
 590            // loop over jets
 591            Int_t nJets = aodJets->GetEntries();
 592            fh1TrackN->Fill((Float_t)nJets);
 593            for(Int_t ij=0; ij<nJets; ++ij){
 594                AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
 595                if(!jet) continue;
 596                AliAODTrack *tmpTr = (AliAODTrack*)jet;
 597                tmpTr->SetFlags(AliESDtrack::kEmbedded);
 598
 599                new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 600                dummy = (*tracks)[nAODtracks-1];
 601
 602                fh1TrackPt->Fill(tmpTr->Pt());
 603                fh2TrackEtaPhi->Fill(tmpTr->Eta(), tmpTr->Phi());
 604            }
 605
 606        } // end: embed jets as 4-momenta
 607
 608
 609     } //end: embed mode with AOD
 610
 611
 612     // === embed mode with toy ===
 613     if(fEmbedMode==kToyTracks){
 614         Int_t nT = (Int_t)(rndm->Uniform(fToyMinNbOfTracks, fToyMaxNbOfTracks)+0.5);
 615
 616         fh1TrackN->Fill((Float_t)nT);
 617
 618         for(Int_t i=0; i<nT; ++i){
 619
 620            Double_t pt = -1.;
 621            if(fToyMinTrackPt!=fToyMaxTrackPt){
 622               if(fToyDistributionTrackPt==0){
 623                  pt = rndm->Uniform(fToyMinTrackPt, fToyMaxTrackPt);
 624               } else {
 625                  while(pt<fToyMinTrackPt||pt>fToyMaxTrackPt){
 626                     pt = rndm->Exp(fToyDistributionTrackPt);   // no power law yet!!
 627                     pt += fToyMinTrackPt;
 628                  }
 629               }
 630            } else {
 631               pt = fToyMinTrackPt;
 632            }
 633            Double_t eta = rndm->Uniform(fToyMinTrackEta,fToyMaxTrackEta);
 634            Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
 635            phi = TVector2::Phi_0_2pi(phi);
 636
 637            if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
 638
 639            Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
 640            Float_t mom[3];
 641            mom[0] = pt;
 642            mom[1] = phi;
 643            mom[2] = theta;
 644
 645            Float_t xyz[3];
 646            xyz[0] = -999.;
 647            xyz[1] = -999.;
 648            xyz[2] = -999.;
 649
 650            AliAODTrack *tmpTr = new AliAODTrack( -999,   // id
 651                                               -999,   // label
 652                                               mom,    // momentum[3]
 653                                               kFALSE, // cartesian
 654                                               xyz,    // position
 655                                               kFALSE, // DCA
 656                                               NULL,   // covMatrix,
 657                                               -99,    // charge
 658                                               0,      // itsClusMap
 659                                               NULL,   // pid
 660                                               NULL,   // prodVertex
 661                                               kFALSE, // used for vertex fit
 662                                               kFALSE, // used for prim vtx fit
 663                                               AliAODTrack::kUndef, // type
 664                                               fToyFilterMap,  // select info
 665                                               -999.    // chi2 per NDF
 666                                             );
 667            tmpTr->SetFlags(AliESDtrack::kEmbedded);
 668
 669            new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
 670            dummy = (*tracks)[nAODtracks-1];
 671
 672            fh1TrackPt->Fill(pt);
 673            fh2TrackEtaPhi->Fill(eta,phi);
 674
 675            delete tmpTr;
 676         }
 677     } //end: embed mode with toy
 678
 679
 680     PostData(1, fHistList);
 681 }
 682
 683
 684 //__________________________________________________________________________
 685 void AliAnalysisTaskFastEmbedding::Terminate(Option_t *)
 686 {
 687     if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Terminate()");
 688
 689     if(fAODfile && fAODfile->IsOpen())
 690     fAODfile->Close();
 691
 692 }
 693
 694 //__________________________________________________________________________
 695 Int_t AliAnalysisTaskFastEmbedding::GetJobID()
 696 {
 697    Int_t id=-1;
 698
 699    const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
 700    //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
 701
 702    if(env && strlen(env)){
 703        id= atoi(env);
 704        AliInfo(Form("Job index %d", id));
 705    }
 706    else{
 707        AliInfo("Job index not found. Okay if running locally.");
 708    }
 709
 710    return id;
 711 }
 712
 713 //__________________________________________________________________________
 714
 715 Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
 716
 717      Int_t nFiles = fAODPathArray->GetEntries();
 718      Int_t n = rndm->Integer(nFiles);
 719      AliInfo(Form("Select AOD file %d", n));
 720      TObjString *objStr = (TObjString*) fAODPathArray->At(n);
 721      if(!objStr){
 722           AliError("Could not get path of aod file.");
 723           return -1;
 724      }
 725      fAODPath = objStr->GetString();
 726
 727      return n;
 728 }
 729
 730 //__________________________________________________________________________
 731
 732 Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
 733
 734     if(trial>3){
 735         AliError("Could not open AOD files, give up ...");
 736         return -1;
 737     }
 738
 739         Int_t rc = 0;
 740         if(fAODPathArray) rc = SelectAODfile();
 741         if(rc<0) return -1;
 742
 743     TDirectory *owd = gDirectory;
 744     if (fAODfile)
 745       fAODfile->Close();
 746     fAODfile = TFile::Open(fAODPath.Data());
 747     owd->cd();
 748     if(!fAODfile){
 749
 750        rc = -1;
 751        if(fAODPathArray){
 752            AliError(Form("Could not open AOD file %s, try another from the list ...", fAODPath.Data()));
 753            rc = OpenAODfile(trial+1);
 754        } else {
 755                AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
 756            }
 757
 758        return rc;
 759     }
 760
 761     fAODtree = (TTree*)fAODfile->Get("aodTree");
 762
 763     if(!fAODtree){
 764        AliError("AOD tree not found.");
 765        return -1;
 766     }
 767
 768     delete fAODevent;
 769     fAODevent = new AliAODEvent();
 770     fAODevent->ReadFromTree(fAODtree);
 771     return rc;  // file position in AOD path array, if array available
 772 }
 773