]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/FASTJETAN/AliAnalysisTaskJetClusterKine.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
FASTJETAN module
[u/mrichter/AliRoot.git] / JETAN / FASTJETAN / AliAnalysisTaskJetClusterKine.cxx
CommitLineData
44b92289 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
72using std::vector;
73
74ClassImp(AliAnalysisTaskJetClusterKine)
75
76AliAnalysisTaskJetClusterKine::~AliAnalysisTaskJetClusterKine(){
77 //
78 // Destructor
79 //
80
81 delete fRef;
82
83 if(fTCAJetsOut)fTCAJetsOut->Delete();
84 delete fTCAJetsOut;
85
86}
87
88//_____________________________________________________________________
89
90AliAnalysisTaskJetClusterKine::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
148AliAnalysisTaskJetClusterKine::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
209Bool_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
219void 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
404void 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//_________________________________________________________________________
411void AliAnalysisTaskJetClusterKine::LocalInit(){
412 // MC handler
413 //
414 // Initialization
415 //
416
417 if(fDebug > 1) printf("AnalysisTaskJetClusterKine::LocalInit() \n");
418 Init();
419}
420
421//_________________________________________________________________________
422void 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
690void AliAnalysisTaskJetClusterKine::Terminate(Option_t* /*option*/){
691 //
692 // Terminate analysis
693 //
694 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
695
696
697}
698
699//_____________________________________________________________________________________
700
701Int_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