Adding includes for EMCAL_Utils and OADB PATH (A. Shabetai)
[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     // 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