Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetClusterKine.cxx
1 // **************************************
2 // Task used for jet finding in the Kine Train (generation and analysis on the fly, no detector effects)
3 // Output is stored in an exchance container
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom3.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include <TChain.h>
28 #include <TRefArray.h>
29 #include <TFile.h>
30 #include <TKey.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3F.h>
34 #include <TProfile.h>
35 #include <TF1.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include  "TDatabasePDG.h"
40 #include <TGrid.h>
41
42 #include "AliAnalysisTaskJetClusterKine.h"
43 #include "AliAnalysisManager.h"
44 #include "AliJetFinder.h"
45 #include "AliJetHeader.h"
46 #include "AliJetReader.h"
47 #include "AliESDEvent.h"
48 #include "AliAODEvent.h"
49 #include "AliAODHandler.h"
50 #include "AliAODExtension.h"
51 #include "AliAODTrack.h"
52 #include "AliAODJet.h"
53 #include "AliVParticle.h" //FK//
54 #include "AliAODMCParticle.h"
55 #include "AliMCEventHandler.h"
56 #include "AliMCEvent.h"
57 #include "AliStack.h"
58 #include "AliGenEventHeader.h" //FK//
59 #include "AliGenPythiaEventHeader.h"
60 #include "AliJetKineReaderHeader.h"
61 #include "AliGenCocktailEventHeader.h"
62 #include "AliInputEventHandler.h"
63 #include "AliAODJetEventBackground.h"
64
65 #include "fastjet/PseudoJet.hh"
66 #include "fastjet/ClusterSequenceArea.hh"
67 #include "fastjet/AreaDefinition.hh"
68 #include "fastjet/JetDefinition.hh"
69 // get info on how fastjet was configured
70 #include "fastjet/config.h"
71
72 using std::vector;
73
74 ClassImp(AliAnalysisTaskJetClusterKine)
75
76 AliAnalysisTaskJetClusterKine::~AliAnalysisTaskJetClusterKine(){
77   //
78   // Destructor
79   //
80
81   delete fRef;
82
83   if(fTCAJetsOut)fTCAJetsOut->Delete();
84   delete fTCAJetsOut;
85   
86 }
87
88 //_____________________________________________________________________
89
90 AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine(): 
91   AliAnalysisTaskSE(),
92   fMcEvent(0x0),  
93   fMcHandler(0x0),
94   fRef(new TRefArray),
95   fTrackTypeGen(kTrackKineCharged), //Kine charged? 
96   fAvgTrials(1),
97   fTrackEtaWindow(0.9),    
98   fTrackPtCut(0.),                                                      
99   fJetOutputMinPt(0.150),
100   fMaxTrackPtInJet(100.),
101   fVtxZCut(10.0),
102   fNonStdBranch(""),
103   fOutContainer(kNoOutput), //FF//
104   fNonStdFile(""),
105   fRparam(1.0), 
106   fAlgorithm(fastjet::kt_algorithm),
107   fStrategy(fastjet::Best),
108   fRecombScheme(fastjet::BIpt_scheme),
109   fAreaType(fastjet::active_area), 
110   fGhostArea(0.01),
111   fActiveAreaRepeats(1),
112   fGhostEtamax(1.5),
113   fTCAJetsOut(0x0),
114   fh1Xsec(0x0),
115   fh1Trials(0x0),
116   fh1PtHard(0x0),
117   fh1PtHardNoW(0x0),  
118   fh1PtHardTrials(0x0),
119   fh1NJetsGen(0x0),
120   fh1NConstGen(0x0),
121   fh1NConstLeadingGen(0x0),
122   fh1PtJetsGenIn(0x0),
123   fh1PtJetsLeadingGenIn(0x0),
124   fh1PtJetConstGen(0x0),
125   fh1PtJetConstLeadingGen(0x0),
126   fh1PtTracksGenIn(0x0),
127   fh1Nch(0x0),
128   fh1Z(0x0), 
129   fh2NConstPt(0x0),
130   fh2NConstLeadingPt(0x0),
131   fh2JetPhiEta(0x0),
132   fh2LeadingJetPhiEta(0x0),
133   fh2JetEtaPt(0x0),
134   fh2LeadingJetEtaPt(0x0),
135   fh2TrackEtaPt(0x0),
136   fh2JetsLeadingPhiEta(0x0),
137   fh2JetsLeadingPhiPt(0x0),
138   fh2JetsLeadingPhiPtW(0x0),
139   fHistList(0x0)  
140 {
141   //
142   // Constructor
143   //
144 }
145
146 //_____________________________________________________________________
147
148 AliAnalysisTaskJetClusterKine::AliAnalysisTaskJetClusterKine(const char* name):
149   AliAnalysisTaskSE(name),
150   fMcEvent(0x0),  
151   fMcHandler(0x0),
152   fRef(new TRefArray),
153   fTrackTypeGen(kTrackKineCharged),  //kine charged?
154   fAvgTrials(1),
155   fTrackEtaWindow(0.9),    
156   fTrackPtCut(0.),                                                      
157   fJetOutputMinPt(0.150),
158   fMaxTrackPtInJet(100.),
159   fVtxZCut(10.0),
160   fNonStdBranch(""),
161   fOutContainer(kNoOutput),//FF//
162   fNonStdFile(""),
163   fRparam(1.0), 
164   fAlgorithm(fastjet::kt_algorithm),
165   fStrategy(fastjet::Best),
166   fRecombScheme(fastjet::BIpt_scheme),
167   fAreaType(fastjet::active_area), 
168   fGhostArea(0.01),
169   fActiveAreaRepeats(1),
170   fGhostEtamax(1.5),
171   fTCAJetsOut(0x0),
172   fh1Xsec(0x0),
173   fh1Trials(0x0),
174   fh1PtHard(0x0),
175   fh1PtHardNoW(0x0),  
176   fh1PtHardTrials(0x0),
177   fh1NJetsGen(0x0),
178   fh1NConstGen(0x0),
179   fh1NConstLeadingGen(0x0),
180   fh1PtJetsGenIn(0x0),
181   fh1PtJetsLeadingGenIn(0x0),
182   fh1PtJetConstGen(0x0),
183   fh1PtJetConstLeadingGen(0x0),
184   fh1PtTracksGenIn(0x0),
185   fh1Nch(0x0),
186   fh1Z(0x0), 
187   fh2NConstPt(0x0),
188   fh2NConstLeadingPt(0x0),
189   fh2JetPhiEta(0x0),
190   fh2LeadingJetPhiEta(0x0),
191   fh2JetEtaPt(0x0),
192   fh2LeadingJetEtaPt(0x0),
193   fh2TrackEtaPt(0x0),
194   fh2JetsLeadingPhiEta(0x0),
195   fh2JetsLeadingPhiPt(0x0),
196   fh2JetsLeadingPhiPtW(0x0),
197   fHistList(0x0)
198 {
199   //
200   // named ctor
201   //
202   DefineOutput(1, TList::Class());  
203   DefineOutput(2, TClonesArray::Class());  
204 }
205
206
207 //_____________________________________________________________________
208
209 Bool_t AliAnalysisTaskJetClusterKine::Notify(){
210   //
211   // Implemented Notify() to read the cross sections
212   // and number of trials from pyxsec.root
213   // 
214   return kTRUE;
215 }
216
217 //_____________________________________________________________________
218
219 void AliAnalysisTaskJetClusterKine::UserCreateOutputObjects(){
220
221    //
222    // Create the output container
223    //
224
225    // Connect the AOD
226
227
228    if(fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
229
230
231    if(fNonStdBranch.Length()!=0){
232       // only create the output branch if we have a name
233       // Create a new branch for jets...
234       //  -> cleared in the UserExec....
235       // here we can also have the case that the brnaches are written to a separate file
236       fTCAJetsOut = new TClonesArray("AliAODJet", 0);
237       fTCAJetsOut->SetName(fNonStdBranch.Data());
238       if(fOutContainer==kAODBranch){   //FF//
239          AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
240       }
241
242     
243       //if(fNonStdFile.Length()!=0){
244         // 
245         // case that we have an AOD extension we need to fetch the jets from the extended output
246         // we identify the extension aod event by looking for the branchname
247         //AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
248         // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
249         //fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
250       //}
251    }
252
253
254    if(!fHistList) fHistList = new TList();
255    fHistList->SetOwner(kTRUE);
256    PostData(1, fHistList); // post data in any case once
257
258
259    if(fOutContainer==kExchCont){
260       fTCAJetsOut->SetOwner();
261       PostData(2, fTCAJetsOut); //FF// post data in any case once
262    }
263
264    Bool_t oldStatus = TH1::AddDirectoryStatus();
265    TH1::AddDirectory(kFALSE);
266
267
268    //
269    //  Histogram
270     
271    const Int_t nBinPt = 100;
272    Double_t binLimitsPt[nBinPt+1];
273    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
274       if(iPt == 0){
275          binLimitsPt[iPt] = 0.0;
276       }else {// 1.0
277          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
278       }
279    }
280   
281    const Int_t nBinPhi = 90;
282    Double_t binLimitsPhi[nBinPhi+1];
283    for(Int_t iPhi = 0; iPhi<=nBinPhi; iPhi++){
284        if(iPhi==0){
285           binLimitsPhi[iPhi] = -1.*TMath::Pi();
286        }else{
287           binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
288        }
289     }
290    
291     const Int_t nBinEta = 40;
292     Double_t binLimitsEta[nBinEta+1];
293     for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
294        if(iEta==0){
295           binLimitsEta[iEta] = -2.0;
296        }else{
297           binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
298      }
299    }
300    
301    const int nChMax = 5000;
302    
303    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
304    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
305   
306    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
307    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
308    
309    
310    fh1NJetsGen = new TH1F("fh1NJetsGen","N reconstructed jets",120,-0.5,119.5);
311    
312    fh1NConstGen = new TH1F("fh1NConstGen","# jet constituents",120,-0.5,119.5);
313    fh1NConstLeadingGen = new TH1F("fh1NConstLeadingGen","jet constituents",120,-0.5,119.5);
314    
315    
316    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
317    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
318    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
319    
320    fh1PtJetsGenIn  = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
321    fh1PtJetsLeadingGenIn = new TH1F("fh1PtJetsLeadingGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
322    fh1PtJetConstGen = new TH1F("fh1PtJetsConstGen","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
323    fh1PtJetConstLeadingGen = new TH1F("fh1PtJetsConstLeadingGen","Gen jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
324    fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("Gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
325    fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
326    
327    
328    fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
329    
330    
331  
332    fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
333    fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
334  
335  
336    fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
337                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
338    fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
339                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
340  
341    fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
342                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
343    fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
344                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
345  
346    fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
347                             nBinEta,binLimitsEta,nBinPt,binLimitsPt);
348  
349  
350  
351    fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
352                                   nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
353    fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
354                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
355
356    fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
357                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
358
359
360   
361
362    const Int_t saveLevel = 3; // large save level more histos
363    if(saveLevel>0){
364       fHistList->Add(fh1Xsec);
365       fHistList->Add(fh1Trials);
366       
367       fHistList->Add(fh1NJetsGen);
368       fHistList->Add(fh1NConstGen);
369       fHistList->Add(fh1NConstLeadingGen);
370       fHistList->Add(fh1PtJetsGenIn);
371       fHistList->Add(fh1PtJetsLeadingGenIn);
372       fHistList->Add(fh1PtTracksGenIn);
373       fHistList->Add(fh1PtJetConstGen);
374       fHistList->Add(fh1PtJetConstLeadingGen);
375       fHistList->Add(fh1Nch);
376       fHistList->Add(fh1Z);
377       
378       fHistList->Add(fh2NConstPt);
379       fHistList->Add(fh2NConstLeadingPt);
380       fHistList->Add(fh2JetPhiEta);
381       fHistList->Add(fh2LeadingJetPhiEta);
382       fHistList->Add(fh2JetEtaPt);
383       fHistList->Add(fh2LeadingJetEtaPt);
384       fHistList->Add(fh2TrackEtaPt);
385       fHistList->Add(fh2JetsLeadingPhiEta);
386       fHistList->Add(fh2JetsLeadingPhiPt);
387       fHistList->Add(fh2JetsLeadingPhiPtW);
388    }
389
390    // =========== Switch on Sumw2 for all histos ===========
391    for(Int_t i=0; i<fHistList->GetEntries(); ++i){
392       TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
393       if(h1){
394          h1->Sumw2();
395          continue;
396       }
397       THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
398       if(hn) hn->Sumw2();
399    }
400    TH1::AddDirectory(oldStatus);
401 }
402 //_________________________________________________________________________
403
404 void AliAnalysisTaskJetClusterKine::Init(){
405    // MC handler
406    if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Init() \n");
407    fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); //FK//
408 }
409
410 //_________________________________________________________________________
411 void AliAnalysisTaskJetClusterKine::LocalInit(){
412    // MC handler
413    //
414    // Initialization
415    //
416
417    if(fDebug > 1) printf("AnalysisTaskJetClusterKine::LocalInit() \n");
418    Init();
419 }
420
421 //_________________________________________________________________________
422 void AliAnalysisTaskJetClusterKine::UserExec(Option_t* /*option*/){
423    // handle and reset the output jet branch 
424   
425    if(fDebug > 1) printf("AliAnalysisTaskJetClusterKine::UserExec \n");
426    if(fTCAJetsOut) fTCAJetsOut->Delete(); //clean TClonesArray
427
428    //
429    // Execute analysis for current event
430    //
431    Init();
432    if(fMcHandler){
433       fMcEvent = fMcHandler->MCEvent(); 
434    }else{
435       if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Handler() fMcHandler=NULL\n");
436       PostData(1, fHistList);
437
438       if(fOutContainer==kExchCont){
439          PostData(2, fTCAJetsOut); //FF//
440       }
441
442       return;
443    }
444    if(!fMcEvent){
445       if(fDebug > 1) printf("AnalysisTaskJetClusterKine::Exec()   fMcEvent=NULL \n");
446       PostData(1, fHistList);
447
448       if(fOutContainer==kExchCont){
449          PostData(2, fTCAJetsOut); //FF// 
450       }
451
452       return;
453    }
454
455    const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
456    Float_t zVtx = vtxMC->GetZ();
457    if(TMath::Abs(zVtx) > fVtxZCut){ //vertex cut
458       PostData(1, fHistList);
459
460       if(fOutContainer==kExchCont){
461          PostData(2, fTCAJetsOut); //FF// 
462       }
463
464       return;
465    }
466
467    fh1Z->Fill(zVtx);
468  
469    fh1Trials->Fill("#sum{ntrials}",1);
470    if(fDebug > 10) Printf("%s:%d",(char*)__FILE__,__LINE__);
471
472    // we simply fetch the mc particles as a list of AliVParticles
473
474    TList genParticles; //list of particles  with above min pT and  good eta range
475
476    Int_t nParticles = GetListOfTracks(&genParticles, fTrackTypeGen); //fill the list
477    fh1Nch->Fill(nParticles);
478
479    if(fDebug>2) Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__, nParticles, genParticles.GetEntries());
480
481    // find the jets....
482
483    vector<fastjet::PseudoJet> inputParticlesKine; 
484
485    for(int i = 0; i < genParticles.GetEntries(); i++){  //loop over generated particles
486       AliVParticle *vp = (AliVParticle*) genParticles.At(i);
487       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
488       jInp.set_user_index(i);
489       inputParticlesKine.push_back(jInp);
490
491       if(fTCAJetsOut){
492          if(i == 0){
493             fRef->Delete(); // make sure to delete before placement new...
494             new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
495          }
496          fRef->Add(vp);  //TRefArray does not work with toy model ...
497       }
498    }//generated particles
499
500    if(inputParticlesKine.size()==0){ //FK//
501       if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
502       PostData(1, fHistList);
503
504       if(fOutContainer==kExchCont){
505          PostData(2, fTCAJetsOut); //FF// 
506       }
507
508       return;
509    } //FK//
510  
511    // run fast jet
512    // employ setters for these...
513    // now create the object that holds info about ghosts                        
514
515    fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
516    fastjet::AreaType areaType =   fastjet::active_area;
517    fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
518    fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
519    fastjet::ClusterSequenceArea clustSeq(inputParticlesKine, jetDef, areaDef); //FK//
520   
521    //range where to compute background
522    Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
523    phiMin = 0;
524    phiMax = 2*TMath::Pi();
525    rapMax = fGhostEtamax - fRparam;
526    rapMin = - fGhostEtamax + fRparam;
527    fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
528  
529
530    const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
531    const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
532
533  
534    fh1NJetsGen->Fill(sortedJets.size()); //the number of found jets
535
536    // loop over all jets an fill information, first on is the leading jet
537
538    Int_t nGen     = inclusiveJets.size();
539    if(inclusiveJets.size()>0){
540
541       //leading Jet
542       AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
543       Double_t area = clustSeq.area(sortedJets[0]);
544       leadingJet.SetEffArea(area,0);
545       Float_t pt = leadingJet.Pt();
546       Float_t phi = leadingJet.Phi();
547       if(phi<0) phi += TMath::Pi()*2.;    
548       Float_t eta = leadingJet.Eta();
549
550       //inclusive jet
551       Int_t nAodOutJets = 0;
552       AliAODJet *aodOutJet = NULL;
553
554
555       // correlation of leading jet with tracks
556   //FK//    TIterator *particleIter = genParticles.MakeIterator();//FK//
557   //FK//    particleIter->Reset();
558   //FK//    AliVParticle *tmpParticle = NULL;
559    //FK//   while((tmpParticle = (AliVParticle*)(particleIter->Next()))){
560    //FK//      Float_t tmpPt = tmpParticle->Pt();
561     //FK//     // correlation
562     //FK//     Float_t tmpPhi =  tmpParticle->Phi();     
563    //FK//      if(tmpPhi<0) tmpPhi+= TMath::Pi()*2.;    
564    //FK//      Float_t dPhi = phi - tmpPhi;
565      //FK//    if(dPhi > TMath::Pi())    dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
566     //FK//     if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();  //-pi,pi     
567      //FK//    fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
568      //FK//    fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
569     //FK//  }  
570    
571       TLorentzVector vecareab;
572       for(int j = 0; j < nGen;j++){ //loop over inclusive jets
573          AliAODJet tmpGenJet (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
574          aodOutJet = NULL;
575          Float_t tmpPt = tmpGenJet.Pt(); //incl jet pt 
576       
577          if((tmpPt > fJetOutputMinPt) && fTCAJetsOut){// cut on the non-background subtracted...
578             aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpGenJet);
579             aodOutJet->GetRefTracks()->Clear();
580             Double_t area1 = clustSeq.area(sortedJets[j]);
581             aodOutJet->SetEffArea(area1,0);
582             fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);  
583             vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());     
584             aodOutJet->SetVectorAreaCharged(&vecareab);
585          }
586
587          fh1PtJetsGenIn->Fill(tmpPt); //incl jet Pt
588          // Fill Spectra with constituentsemacs
589          const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
590
591          fh1NConstGen->Fill(constituents.size());  //number of constituents
592          fh2NConstPt->Fill(tmpPt,constituents.size()); //number of constituents vs jet pt 
593
594          // loop over constiutents and fill spectrum
595          AliVParticle *partLead = 0;
596          Float_t ptLead = -1;
597
598          for(unsigned int ic = 0; ic < constituents.size(); ic++){
599             AliVParticle *part = (AliVParticle*) genParticles.At(constituents[ic].user_index());
600             if(!part) continue;
601             fh1PtJetConstGen->Fill(part->Pt()); //pt of constituents
602
603             if(aodOutJet){
604                aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));//FK//
605
606                if(part->Pt()>fMaxTrackPtInJet){
607                   aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
608                }
609             }
610
611             if(part->Pt()>ptLead){
612                ptLead = part->Pt();
613                partLead = part;
614             }
615
616             if(j==0) fh1PtJetConstLeadingGen->Fill(part->Pt()); //pt of leading jet constituents
617          }
618
619          if(partLead){
620             if(aodOutJet){
621                //set pT of leading constituent of jet
622                aodOutJet->SetPtLeading(partLead->Pt());
623             }
624          }
625     
626          // correlation
627          Float_t tmpPhi =  tmpGenJet.Phi(); //incl jet phi
628          Float_t tmpEta =  tmpGenJet.Eta(); //incl jet eta
629
630          if(tmpPhi<0) tmpPhi += TMath::Pi()*2.;        
631          fh2JetPhiEta->Fill(tmpGenJet.Phi(),tmpEta);
632          fh2JetEtaPt->Fill(tmpEta,tmpPt);
633
634          if(j==0){   //leading jet
635             fh1PtJetsLeadingGenIn->Fill(tmpPt); //leading jet pt
636             fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
637             fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
638             fh1NConstLeadingGen->Fill(constituents.size());  //number of constituents in leading jet
639             fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
640             continue;
641          }
642
643          Float_t dPhi = phi - tmpPhi;
644          if(dPhi > TMath::Pi())     dPhi = dPhi - 2.*TMath::Pi(); //-pi,pi
645          if(dPhi<(-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi(); //-pi.pi     
646          Float_t dEta = eta - tmpGenJet.Eta();
647          fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
648          fh2JetsLeadingPhiPt->Fill(dPhi,pt);
649          fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
650       }// loop over reconstructed jets
651     //delete particleIter;
652    }//number of jets>0
653
654  
655    // fill track information
656    Int_t nTrackOver = genParticles.GetSize();
657    // do the same for tracks and jets
658
659    if(nTrackOver>0){
660       TIterator *particleIter = genParticles.MakeIterator();
661       AliVParticle *tmpGen = 0;
662       particleIter->Reset();
663
664       while((tmpGen = (AliVParticle*)(particleIter->Next()))){
665          Float_t tmpPt = tmpGen->Pt();
666          Float_t tmpEta = tmpGen->Eta();
667          fh1PtTracksGenIn->Fill(tmpPt);
668          fh2TrackEtaPt->Fill(tmpEta,tmpPt);
669       }  
670       delete particleIter;
671    }
672
673
674
675
676
677    if(fDebug > 2){
678       if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
679    }
680    PostData(1, fHistList);
681  
682   if(fOutContainer==kExchCont){
683       PostData(2, fTCAJetsOut); //FF// 
684    }
685
686
687 }
688 //____________________________________________________________________________
689
690 void AliAnalysisTaskJetClusterKine::Terminate(Option_t* /*option*/){
691   //
692   // Terminate analysis
693   //
694   if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
695
696     
697 }
698
699 //_____________________________________________________________________________________
700
701 Int_t  AliAnalysisTaskJetClusterKine::GetListOfTracks(TList *list, Int_t type){
702
703   //
704   // get list of tracks/particles for different types
705   //
706
707    if(fDebug>2) Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
708
709    Int_t iCount = 0; //number of particles 
710
711    if(type ==  kTrackKineAll || type == kTrackKineCharged){
712       if(! fMcEvent) return iCount;
713
714       //we want to have alivpartilces so use get track
715       for(int it = 0; it < fMcEvent->GetNumberOfTracks(); ++it){
716          if(!fMcEvent->IsPhysicalPrimary(it)) continue;
717
718          AliMCParticle* part = (AliMCParticle*)fMcEvent->GetTrack(it);
719          if(TMath::Abs(part->Eta()) > fTrackEtaWindow) continue;   //apply eta cut
720          if(part->Pt() < fTrackPtCut) continue;  //apply pT cut
721
722          if(type == kTrackKineAll){  // full jets
723             list->Add(part);
724             iCount++;
725
726          }else if(type == kTrackKineCharged){  //charged jets
727             if(part->Particle()->GetPDG()->Charge()==0) continue;
728             list->Add(part);
729             iCount++;
730          }
731       }
732    }
733   
734    list->Sort();  //?//
735
736    return iCount; //number of particles in the list
737 }
738
739