e131bc4d8d116e22b170520ea970c6c743e6b791
[u/mrichter/AliRoot.git] / 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