]>
Commit | Line | Data |
---|---|---|
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 | ||
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 |